text:
"""
```
    GammaCorrection <: AbstractHistogramAdjustmentAlgorithm
    GammaCorrection(; gamma = 1)

    adjust_histogram([T,] img, f::GammaCorrection)
    adjust_histogram!([out,] img, f::GammaCorrection)
```

Returns a gamma corrected image.

# Details

Gamma correction is a non-linear  transformation given by the relation
```math
f(x) = x^\\gamma \\quad \\text{for} \\; x \\in \\mathbb{R}, \\gamma > 0.
```
It is called a *power law* transformation because one quantity varies as a power
of another quantity.

Gamma correction has historically been used to preprocess
an image to compensate for the fact that the intensity of light generated by a
physical device is not usually a linear function of the applied signal but
instead follows a power law [1]. For example, for many Cathode Ray Tubes (CRTs) the
emitted light intensity on the display is approximately equal to the voltage
raised to the power of γ, where γ ∈ [1.8, 2.8]. Hence preprocessing a raw image with
an exponent of 1/γ  would have ensured a linear response to brightness.

Research in psychophysics has also established an [empirical  power law
](https://en.wikipedia.org/wiki/Stevens%27s_power_law)  between light intensity and perceptual
brightness. Hence, gamma correction often serves as a useful image enhancement
tool.

# Options

Various options for the parameters of the `adjust_histogram` function and the
`Gamma` type are described in more detail below.

## Choices for `img`

The function can handle a variety of input types. The returned
image depends on the input type.

For colored images, the input is converted to YIQ type and the Y channel is
gamma corrected. This is the combined with the I and Q channels and the
resulting image converted to the same type as the input.

## Choice for `gamma`

The `gamma` value must be a non-zero positive number. A `gamma` value less than
one will yield a brighter image whereas a value greater than one will produce a
darker image. If left unspecified a default value of one is assumed.

# Example

```julia
using ImageContrastAdjustment, ImageView

# Create an example image consisting of a linear ramp of intensities.
n = 32
intensities = 0.0:(1.0/n):1.0
img = repeat(intensities, inner=(20,20))'

# Brighten the dark tones.
imgadj = adjust_histogram( img, GammaCorrection(gamma = 1/2))

# Display the original and adjusted image.
imshow(img)
imshow(imgadj)
```

# References
1. W. Burger and M. J. Burge. *Digital Image Processing*. Texts in Computer Science, 2016. [doi:10.1007/978-1-4471-6684-9](https://doi.org/10.1007/978-1-4471-6684-9)
"""
Base.@kwdef struct GammaCorrection{T <: Real} <: AbstractHistogramAdjustmentAlgorithm
    gamma::T = 1.0
end

function (f::GammaCorrection)(out::GenericGrayImage, img::GenericGrayImage)
    out .= img
    correct_gamma!(out, f.gamma)
    return out
end

# TODO Expand this to more generic gray color types.
function correct_gamma!(img::AbstractArray{Gray{T}},  gamma::Real) where T <: FixedPointNumbers.Normed
    γ = Float64(gamma)
    # Create a lookup-table for the gamma transformation of the grayvalues.
    raw_type = FixedPointNumbers.rawtype(T)
    table = zeros(T, typemax(raw_type) + 1)
    for i in zero(raw_type):typemax(raw_type)
        table[i + 1] = T((i / typemax(raw_type))^γ)
    end
    # Map the pixels to their new grayvalues.
    map!(img,img) do p
        table[p.val.i + 1]
    end
end

function correct_gamma!(img::GenericGrayImage,  gamma::Real)
    γ = Float64(gamma)
    T = eltype(img)
    map!(img,img) do val
        if isnan(val)
            return val
        else
            return  T <: Integer ? round(Int, val^γ) : val^γ
        end
    end
end

function (f::GammaCorrection)(out::AbstractArray{<:Color3}, img::AbstractArray{<:Color3})
    T = eltype(out)
    yiq = convert.(YIQ, img)
    yiq_view = channelview(yiq)
    adjust_histogram!(view(yiq_view,1,:,:), f)
    out .= convert.(T, yiq)
end

(f::GammaCorrection)(out::GenericGrayImage, img::AbstractArray{<:Color3}) =
    f(out, of_eltype(Gray, img))


text:
module Extend

using Reexport: @reexport

# Only exports things one typically dispatches on when extending NaiveNASlib with new functions, e.g. abstract types and traits

# Vertex types
@reexport using ..NaiveNASlib: AbstractVertex, InputVertex, InputSizeVertex, CompVertex, MutationVertex, OutputsVertex

# Mutation traits
@reexport using ..NaiveNASlib: trait, MutationTrait, DecoratingTrait, Immutable, MutationSizeTrait, SizeAbsorb, SizeTransparent, SizeStack, SizeInvariant, FixedSizeTrait

# Size strategies
@reexport using ..NaiveNASlib: AbstractΔSizeStrategy, AbstractJuMPΔSizeStrategy, DecoratingJuMPΔSizeStrategy, AbstractAfterΔSizeStrategy

@reexport using ..NaiveNASlib: AbstractAlignSizeStrategy
@reexport using ..NaiveNASlib: AbstractConnectStrategy

@reexport using ..NaiveNASlib: base, parselect, vertex

end

ext: jl


text:
module VQE4NPP

using Yao, YaoExtensions
using Flux
using Random
using Statistics

include("TrainFunc.jl")
include("OtherFunc.jl")

end


text:
using MEngProject
using Documenter

makedocs(;
    modules=[MEngProject],
    authors="nialltc <39237476+nialltc@users.noreply.github.com> and contributors",
    repo="https://github.com/nialltc/MEngProject.jl/blob/{commit}{path}#L{line}",
    sitename="MEngProject.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://nialltc.github.io/MEngProject.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/nialltc/MEngProject.jl",
)


text:
@testset "construct" begin
   @test Double64(one(Double64)) === one(Double64)
   @test Double64(one(Double32)) === one(Double64)
   @test Double64(one(Double16)) === one(Double64)
   @test Double32(one(Double64)) === one(Double32)
   @test Double32(one(Double32)) === one(Double32)
   @test Double32(one(Double16)) === one(Double32)
   @test Double16(one(Double64)) === one(Double16)
   @test Double16(one(Double32)) === one(Double16)
   @test Double16(one(Double16)) === one(Double16)

   @test Double64(one(Float64)) === one(Double64)
   @test Double64(one(Float32)) === one(Double64)
   @test Double64(one(Float16)) === one(Double64)
   @test Double32(one(Float64)) === one(Double32)
   @test Double32(one(Float32)) === one(Double32)
   @test Double32(one(Float16)) === one(Double32)
   @test Double16(one(Float64)) === one(Double16)
   @test Double16(one(Float32)) === one(Double16)
   @test Double16(one(Float16)) === one(Double16)

   @test DoubleFloat(one(Double64)) === one(Double64)
   @test DoubleFloat(one(Double32)) === one(Double32)
   @test DoubleFloat(one(Double16)) === one(Double16)

   @test DoubleFloat(1.0) === one(Double64)
   @test DoubleFloat(1.0f0) === one(Double32)
   @test DoubleFloat(Float16(1.0)) === one(Double16)

   @test ComplexDF64(1, 1) === ComplexDF64(1.0, 1.0)
   @test ComplexDF32(1, 1) === ComplexDF32(1.0f0, 1.0f0)
   @test ComplexDF16(1, 1) === ComplexDF16(Float16(1.0), Float16(1.0))
   @test ComplexDF64(1//1, -1//1) === ComplexDF64(1.0, -1.0)
   @test ComplexDF32(1//1, -1//1) === ComplexDF32(1.0f0, -1.0f0)
   @test ComplexDF16(1//1, -1//1) === ComplexDF16(Float16(1.0), -Float16(1.0))
end

@testset "intraconstruct" begin
   @test Double64(Double32(1)) === Double64(1)
   @test Double64(Double16(1)) === Double64(1)
   @test Double32(Double16(1)) === Double32(1)

   @test Double64(1.0) === Double64(1)
   @test Double64(1.0f0) === Double64(1)
   @test Double64(Float16(1.0)) === Double64(1)

   @test Double32(1.0) === Double32(1)
   @test Double32(1.0f0) === Double32(1)
   @test Double32(Float16(1.0)) === Double32(1)

   @test Double16(1.0) == Double16(1)
   @test Double16(1.0f0) == Double16(1)
   @test Double16(Float16(1.0)) == Double16(1)

   @test Double64(Int64(1)) === Double64(1.0)
   @test Double32(Int64(1)) === Double32(1.0)
   @test Double16(Int64(1)) === Double16(1.0)

   @test Double64(Int32(1)) === Double64(1.0)
   @test Double32(Int32(1)) === Double32(1.0)
   @test Double16(Int32(1)) === Double16(1.0)

   @test Double64(Int16(1)) === Double64(1.0)
   @test Double32(Int16(1)) === Double32(1.0)
   @test Double16(Int16(1)) === Double16(1.0)
end

@testset "tuple construct" begin
   tup = (1.0e3, 1.0e-2)
   @test Double64(tup) == DoubleFloat(Float64(HI(tup)), Float64(LO(tup)))
   @test Double32(tup) == DoubleFloat(Float32(HI(tup)), Float32(LO(tup)))
   @test Double16(tup) == DoubleFloat(Float16(HI(tup)), Float16(LO(tup)))
end

@testset "tuple deconstruct" begin
   tup = (1.0e3, 1.0e-2)
   @test HI(tup) == tup[1]
   @test LO(tup) == tup[2]
   @test HILO(tup) == tup
end


text:
@testset "Utils" begin
    # rowrange
    n = 23
    v = rand(n)
    for m in 1:9
        step = floor(Int, n/m)
        for i in 1:m
            @test QuantizedArrays.rowrange(n, m, i) ==
                step*(i-1)+1 : step*i
        end
    end

    # quantized_eltype
    for T in subtypes(Unsigned)
        k = rand(typemin(T):typemax(T))
        @test QuantizedArrays.quantized_eltype(k+1) == T
    end

    # Quantization type
    for qt in [:sample, :pq, :opq, :rvq]
        qt == :rvq && @test QuantizedArrays.quantization_type(qt) isa QuantizedArrays.AdditiveQuantization
        qt != :rvq && @test QuantizedArrays.quantization_type(qt) isa QuantizedArrays.OrthogonalQuantization
    end
end


text:
struct UncertainDepthValueDataset <: UncertainIndexValueDataset
    indices::Vector{AbstractUncertainValue}
    values::Vector{AbstractUncertainValue}
end

export UncertainDepthValueDataset


text:
module EMNIST_Tests
    using Test
    using DataDeps
    using MLDatasets

    if parse(Bool, get(ENV, "CI", "false"))
        @info "CI detected: skipping dataset download"
    else
        data_dir = withenv("DATADEPS_ALWAY_ACCEPT"=>"true") do
            datadep"EMNIST"
        end

        @testset "EMNIST" begin
            X_train = EMNIST.Balanced.traindata()
            X_test = EMNIST.Balanced.testdata()
            Y_train = EMNIST.Balanced.trainlabels()
            Y_test = EMNIST.Balanced.testlabels()
            @test X_train isa Array{UInt8,3}
            @test X_test isa Array{UInt8,3}
            @test Y_train isa Array{Float64,2}
            @test Y_test isa Array{Float64,2}
            @test size(X_train) == (112800, 28, 28)
            @test size(X_test) == (18800, 28, 28)
            @test size(Y_train) == (112800, 1)
            @test size(Y_test) == (18800, 1)

            X_train = EMNIST.ByClass.traindata()
            X_test = EMNIST.ByClass.testdata()
            Y_train = EMNIST.ByClass.trainlabels()
            Y_test = EMNIST.ByClass.testlabels()
            @test X_train isa Array{UInt8,3}
            @test X_test isa Array{UInt8,3}
            @test Y_train isa Array{Float64,2}
            @test Y_test isa Array{Float64,2}
            @test size(X_train) == (697932, 28, 28)
            @test size(X_test) == (116323, 28, 28)
            @test size(Y_train) == (697932, 1)
            @test size(Y_test) == (116323, 1)

            X_train = EMNIST.ByMerge.traindata()
            X_test = EMNIST.ByMerge.testdata()
            Y_train = EMNIST.ByMerge.trainlabels()
            Y_test = EMNIST.ByMerge.testlabels()
            @test X_train isa Array{UInt8,3}
            @test X_test isa Array{UInt8,3}
            @test Y_train isa Array{Float64,2}
            @test Y_test isa Array{Float64,2}
            @test size(X_train) == (697932, 28, 28)
            @test size(X_test) == (116323, 28, 28)
            @test size(Y_train) == (697932, 1)
            @test size(Y_test) == (116323, 1)

            X_train = EMNIST.Digits.traindata()
            X_test = EMNIST.Digits.testdata()
            Y_train = EMNIST.Digits.trainlabels()
            Y_test = EMNIST.Digits.testlabels()
            @test X_train isa Array{UInt8,3}
            @test X_test isa Array{UInt8,3}
            @test Y_train isa Array{Float64,2}
            @test Y_test isa Array{Float64,2}
            @test size(X_train) == (240000, 28, 28)
            @test size(X_test) == (40000, 28, 28)
            @test size(Y_train) == (240000, 1)
            @test size(Y_test) == (40000, 1)

            X_train = EMNIST.Letters.traindata()
            X_test = EMNIST.Letters.testdata()
            Y_train = EMNIST.Letters.trainlabels()
            Y_test = EMNIST.Letters.testlabels()
            @test X_train isa Array{UInt8,3}
            @test X_test isa Array{UInt8,3}
            @test Y_train isa Array{Float64,2}
            @test Y_test isa Array{Float64,2}
            @test size(X_train) == (124800, 28, 28)
            @test size(X_test) == (20800, 28, 28)
            @test size(Y_train) == (124800, 1)
            @test size(Y_test) == (20800, 1)

            X_train = EMNIST.MNIST.traindata()
            X_test = EMNIST.MNIST.testdata()
            Y_train = EMNIST.MNIST.trainlabels()
            Y_test = EMNIST.MNIST.testlabels()
            @test X_train isa Array{UInt8,3}
            @test X_test isa Array{UInt8,3}
            @test Y_train isa Array{Float64,2}
            @test Y_test isa Array{Float64,2}
            @test size(X_train) == (60000, 28, 28)
            @test size(X_test) == (10000, 28, 28)
            @test size(Y_train) == (60000, 1)
            @test size(Y_test) == (10000, 1)
        end
    end

end


text:
"""
Generate a history of requests from given MDP.

This method is used to genereate sequences of requests for the experiments.
"""
function simulate_trace(m::PMDP,rng::AbstractRNG)
    hr = HistoryRecorder(max_steps=selling_period_end(m), rng = rng)
    reject = FunctionPolicy(s->1000.)
    history = simulate(hr, m, reject)
end

"""
Apply BSON specific fixes to loaded trace data
"""
function load_traces(fname)
    data = load(fname)

    # fix type of traces
    traces = data[:traces]
    @assert all(y->typeof(y)==typeof(traces[1]), traces) "traces contain elements of different types!"
    data[:traces] = [traces...] # this repackages the loaded array of typr Array{Any,1} into Array{SimHistory, 1}

    # fix type of pp_params
    pp_params = data[:pp_params]
    data[:pp_params] = Dict(pp_params...)

    return Dict(data)
end


text:
using MicrobiomeAnalysis
using Statistics # , SummarizedExperiments
using Test

@testset "MicrobiomeAnalysis.jl" begin

    se = SummarizedExperiments.exampleobject(40, 20)
    foo_assay = assay(se, "foo")
    bar_assay = assay(se, "bar")

    log10_assay = log10(foo_assay, 1)
    pa_assay = pa(bar_assay, 1)
    relabund_assay = relabund(foo_assay)
    z_assay1 = ztransform(bar_assay)
    clr_assay1 = clr(foo_assay, 1)

    transform!(se, "foo", log10, 1, output_name = "log10_assay")
    transform!(se, "bar", pa, 1, output_name = "pa_assay")
    transform!(se, "foo", relabund, output_name = "relabund_assay")
    z_assay2 = transform(se, "bar", ztransform)
    clr_assay2 = transform(se, "foo", clr, 1)

    # are the rounded means equal to each other up to a tolerance of 1 digit?
    @test round(mean(log10_assay), digits = 1) == round(mean(assay(se, "log10_assay")), digits = 1)
    # is each element in A equal to each element in B?
    @test pa_assay == assay(se, "pa_assay")
    # are the rounded means equal to each other up to a tolerance of 15 digits?
    @test round(mean(relabund_assay), digits = 15) == round(mean(assay(se, "relabund_assay")), digits = 15)
    # is each element in A equal to each element in B?
    @test z_assay1 == z_assay2
    # is the absolute value of the difference between the means less than the tolerance?
    @test abs(mean(clr_assay1) - mean(clr_assay2)) < 10e-16

end


text:
using Distributions, Plots; pyplot()

alpha, beta = 8, 2
prior(lam) = pdf(Gamma(alpha, 1/beta), lam)
data = [2,1,0,0,1,0,2,2,5,2,4,0,3,2,5,0]

like(lam) = *([pdf(Poisson(lam),x) for x in data]...)
posteriorUpToK(lam) = like(lam)*prior(lam)

delta = 10^-4.
lamRange = 0:delta:10
K = sum([posteriorUpToK(lam)*delta for lam in lamRange])
posterior(lam) = posteriorUpToK(lam)/K

bayesEstimate = sum([lam*posterior(lam)*delta for lam in lamRange])

newAlpha, newBeta = alpha + sum(data), beta + length(data)
closedFormBayesEstimate = mean(Gamma(newAlpha, 1/newBeta))

println("Computational Bayes Estimate: ", bayesEstimate)
println("Closed form Bayes Estimate: ", closedFormBayesEstimate)

plot(lamRange, prior.(lamRange),
	c=:blue, label="Prior distribution")
plot!(lamRange, posterior.(lamRange),
	c=:red, label="Posterior distribution",
	xlims=(0, 10), ylims=(0, 1.2),
    	xlabel=L"\lambda",ylabel="Density")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# tools to help manipulate graphs

using JEMSS
import LightGraphs
using SparseArrays

travelModeString(i::Int) = "mode_$i"

function checkNodeIndices(nodes::Vector{Node})
	result = all(i -> i == nodes[i].index, 1:length(nodes)) # true if ok
	result || @warn("Node indices are not 1:n")
	return result
end

function checkArcIndices(arcs::Vector{Arc})
	result = all(i -> i == arcs[i].index, 1:length(arcs)) # true if ok
	result || @warn("Arc indices are not 1:n")
	return result
end

# check arc.fromNodeIndex and arc.toNodeIndex values are in expected range
function checkArcNodeIndices(numNodes::Int, arcs::Vector{Arc})
	result = all(i -> (1 <= arcs[i].fromNodeIndex <= numNodes) && (1 <= arcs[i].toNodeIndex <= numNodes), 1:length(arcs)) # true if ok
	result || @warn("Arc from/to node indices are not all within 1:numNodes")
	return result
end
checkArcNodeIndices(nodes::Vector{Node}, arcs::Vector{Arc}) = checkArcNodeIndices(length(nodes), arcs)

# get number of travel modes for an arc
function getNumTravelModes(arc::Arc)
	i = 1
	while haskey(arc.fields, travelModeString(i))
		i += 1
	end
	numModes = i-1
	@assert(numModes >= 1, "Arcs should contain at least one travel mode field, e.g. '$(travelModeString(1))'")
	return numModes
end
getNumTravelModes(arcs::Vector{Arc}) = getNumTravelModes(arcs[1]) # assume all arcs have same fields

# get travel times for each travel mode of an arc
function getArcTravelTimes(arc::Arc)
	numModes = getNumTravelModes(arcs)
	travelTimes = Vector{Float}(undef, numModes)
	for i = 1:numModes
		travelTimes[i] = arc.fields[travelModeString(i)]
	end
	return travelTimes
end

# return vector of travel times for arcs,
# first dimension is for travel modes, second is for arcs
function getArcsTravelTimes(arcs::Vector{Arc})
	numArcs = length(arcs)
	numModes = getNumTravelModes(arcs)
	travelTimes = Array{Float,2}(undef, numModes, numArcs)
	for i = 1:numModes
		modeString = travelModeString(i)
		for j = 1:numArcs
			travelTimes[i,j] = arcs[j].fields[modeString]
		end
	end
	return travelTimes
end

function nodesAddField!(nodes::Vector{Node}, field::String; default::Any = nothing, overwrite::Bool = false)
	if !overwrite && haskey(nodes[1].fields, field) return end
	for node in nodes
		node.fields[field] = default
	end
end

function arcsAddField!(arcs::Vector{Arc}, field::String; default::Any = nothing, overwrite::Bool = false)
	if !overwrite && haskey(arcs[1].fields, field) return end
	for arc in arcs
		arc.fields[field] = default
	end
end

nodesDeleteField!(nodes::Vector{Node}, field::String) = for node in nodes delete!(node.fields, field) end
arcsDeleteField!(arcs::Vector{Arc}, field::String) = for arc in arcs delete!(arc.fields, field) end

function JEMSS.writeArcsFile(arcsFilename::String, arcs::Vector{Arc}, arcForm::String)
	travelTimes = getArcsTravelTimes(arcs)
	writeArcsFile(arcsOutputFilename, arcs, travelTimes, arcForm)
end

function createLightGraph(numNodes::Int, arcs::Vector{Arc})
	lightGraph = LightGraphs.DiGraph(numNodes)
	for arc in arcs
		LightGraphs.add_edge!(lightGraph, arc.fromNodeIndex, arc.toNodeIndex)
	end
	return lightGraph
end
createLightGraph(nodes::Vector{Node}, arcs::Vector{Arc}) = createLightGraph(length(nodes), arcs)

# tag the elements (nodes and arcs) that contribute to the largest strongly connected component
# tagging is done by adding a field with name given by kwarg `lsccHeader`
function graphTagLargestComponentElts!(nodes::Vector{Node}, arcs::Vector{Arc};
	lsccHeader::String = "in_largest_component")

	# shorthand
	numNodes = length(nodes)
	numArcs = length(arcs)

	lightGraph = createLightGraph(nodes, arcs)

	# find strongly connected components; note that Graphs.strongly_connected_components() did not work correctly (the output was not what I expected), so LightGraphs.strongly_connected_components() is needed instead
	components = LightGraphs.strongly_connected_components(lightGraph)
	largestComponent = components[argmax([length(c) for c in components])] # node indices of largest strongly connected component

	# tag nodes and arcs in largest strongly connected component
	tagNode = fill(false, numNodes) # tagNode[i] = true if nodes[i] in largest component
	tagArc = fill(false, numArcs) # tagArc[i] = true if arcs[i] in largest component
	for i in largestComponent
		tagNode[i] = true
	end
	for (i,arc) in enumerate(arcs)
		if tagNode[arc.fromNodeIndex] && tagNode[arc.toNodeIndex]
			tagArc[i] = true
		end
	end

	# tag by adding fields to nodes and arcs
	for i = 1:numNodes
		nodes[i].fields[lsccHeader] = tagNode[i]
	end
	for i = 1:numArcs
		arcs[i].fields[lsccHeader] = tagArc[i]
	end

	return
end

# keep nodes and arcs where filters return true
# indices will be renumbered
# arcs with either (to/from) node removed are not automatically removed
function graphRemoveElts!(nodes::Vector{Node}, arcs::Vector{Arc};
	nodeFilter::Function = (x->true), arcFilter::Function = (x->true))

	# before filtering, create mapping from old indices to new
	keepNodeIndices = findall(nodeFilter, nodes)
	nodeNewIndex = fill(nullIndex, length(nodes))
	nodeNewIndex[keepNodeIndices] = 1:length(keepNodeIndices) # nodeNewIndex[i] gives new node index for old nodes[i]

	# filter nodes, renumber nodes from 1 to n
	filter!(nodeFilter, nodes)
	for (i,node) in enumerate(nodes)
		node.index = i
	end

	# filter arcs, renumber nodes from 1 to n
	filter!(arcFilter, arcs)
	for (i,arc) in enumerate(arcs)
		arc.index = i
	end

	# connect arcs to nodes
	# arc.fromNodeIndex or arc.toNodeIndex will become nullIndex if node does not exist
	for arc in arcs
		arc.fromNodeIndex = nodeNewIndex[arc.fromNodeIndex]
		arc.toNodeIndex = nodeNewIndex[arc.toNodeIndex]
	end

	# check that all needed nodes exist for filtered arcs
	checkArcNodeIndices(nodes, arcs)

	return
end

# remove elements outside of largest strongly connected component
function graphKeepLargestComponent!(nodes::Vector{Node}, arcs::Vector{Arc};
	lsccHeader::String = "in_largest_component")
	graphTagLargestComponentElts!(nodes, arcs; lsccHeader = lsccHeader)
	nodeFilter(node::Node) = node.fields[lsccHeader]
	arcFilter(arc::Arc) = arc.fields[lsccHeader]
	graphRemoveElts!(nodes, arcs, nodeFilter = nodeFilter, arcFilter = arcFilter)
	nodesDeleteField!(nodes, lsccHeader)
	arcsDeleteField!(arcs, lsccHeader)
	return
end

# remove arcs from graph that are connected to/from nodes outside of range
function graphRemoveDisconnectedArcs!(nodes::Vector{Node}, arcs::Vector{Arc})
	numNodes = length(nodes)
	arcFilter(arc::Arc) = (1 <= arc.fromNodeIndex <= numNodes) && (1 <= arc.toNodeIndex <= numNodes)
	graphRemoveElts!(nodes, arcs, arcFilter = arcFilter)
end

# return dict with tuple of node indices (from, to) as keys, and vector of arc indices as values
function graphNodePairArcIndices(arcs::Vector{Arc})
	nodePairArcIndices = Dict{Tuple{Int,Int}, Vector{Int}}() # nodePairArcIndices[i,j] gives indices of arcs going from node i to node j
	for arc in arcs
		i = arc.fromNodeIndex
		j = arc.toNodeIndex
		if !haskey(nodePairArcIndices, (i,j))
			nodePairArcIndices[(i,j)] = []
		end
		push!(nodePairArcIndices[(i,j)], arc.index)
	end
	return nodePairArcIndices
end

# find duplicate arcs (same 'to' and 'from' nodes),
# return vector of vector of indices of duplicate arcs
function graphFindDuplicateArcs(arcs::Vector{Arc})
	nodePairArcIndices = graphNodePairArcIndices(arcs)
	duplicateArcsSets = []
	for (nodePair, arcIndices) in nodePairArcIndices
		if length(arcIndices) > 1
			push!(duplicateArcsSets, arcIndices)
		end
	end
	return duplicateArcsSets
end

graphContainsDuplicateArcs(arcs::Vector{Arc}) = graphFindDuplicateArcs(arcs) != []
function checkGraphContainsNoDuplicateArcs(arcs::Vector{Arc})
	result = !graphContainsDuplicateArcs(arcs) # true if ok
	result || @warn("There are duplicate arcs.")
	return result
end

# merge arc2 into arc1, keeping minimum travel time (from both arcs) for each travel mode in arc1
# arc2 will need to be removed in separate step
function mergeDuplicateArcs!(arc1::Arc, arc2::Arc)
	@assert(arc1.fromNodeIndex == arc2.fromNodeIndex)
	@assert(arc1.toNodeIndex == arc2.toNodeIndex)
	numModes = getNumTravelModes(arc1)
	for i = 1:numModes
		modeString = travelModeString(i)
		arc1.fields[modeString] = min(arc1.fields[modeString], arc2.fields[modeString])
	end
	# assume that arc.distance is same for both arcs
end
mergeDuplicateArcs!(arcs::Vector{Arc}, i::Int, j::Int) = mergeDuplicateArcs!(arcs[i], arcs[j])

# merge any duplicate arcs into one
# will take the shortest travel time from any merged arcs, for each travel mode
function graphMergeDuplicateArcs!(nodes::Vector{Node}, arcs::Vector{Arc};
	mergeArcsFunction::Function = mergeDuplicateArcs!, mergeHeader::String = "merge_result")

	arcsAddField!(arcs, mergeHeader, default = "unmerged")

	# find which arcs are duplicates
	duplicateArcsSets = graphFindDuplicateArcs(arcs) # duplicateArcsSets[i] has a vector of indices of duplicate arcs
	for duplicateArcs in duplicateArcsSets
		i = duplicateArcs[1]
		arcs[i].fields[mergeHeader] = "merged"
		for j in duplicateArcs[2:end]
			mergeArcsFunction(arcs, i, j)
			arcs[j].fields[mergeHeader] = "merged_in"
		end
	end

	# remove duplicate arcs
	arcFilter(arc::Arc) = arc.fields[mergeHeader] != "merged_in"
	graphRemoveElts!(nodes, arcs, arcFilter = arcFilter)
end

# Divide an arc into even sections a given number of times, 0 divides makes no change.
# Only some nodes and arcs fields are updated (node offRoadAccess, arc travel times and osm weight), the rest are unchanged but may need changing.
function graphDivideArc!(nodes::Vector{Node}, arcs::Vector{Arc}; arcIndex::Int = 0, numDivides::Int = 0)
	if numDivides <= 0 return end
	@assert(1 <= arcIndex <= length(arcs))

	# shorthand
	numNodes = length(nodes)
	numArcs = length(arcs)
	arc = arcs[arcIndex]
	fromNode = nodes[arc.fromNodeIndex]
	toNode = nodes[arc.toNodeIndex]

	# keep track of which nodes have been added, and which arcs have been divided
	addedHeader = "added"
	divArcHeader = "divided_arc_index"
	nodesAddField!(nodes, addedHeader; default = false)
	arcsAddField!(arcs, divArcHeader; default = nullIndex)

	# change some fields of arc
	for i = 1:getNumTravelModes(arc)
		arc.fields[travelModeString(i)] /= numDivides + 1
	end
	if haskey(arc.fields, "osm_weight")
		arc.fields["osm_weight"] /= numDivides + 1
	end
	arc.distance /= numDivides + 1
	arc.fields[divArcHeader] = arc.index

	tempNode = deepcopy(fromNode) # could copy any node
	tempNode.fields[addedHeader] = true
	tempArc = deepcopy(arc)
	# some fields of these new nodes and arcs will have incorrect/irrelevant information and should be ignored

	for i = 1:numDivides
		# add node
		newNode = deepcopy(tempNode)
		newNode.index = numNodes + i
		newNode.location = linearInterpLocation(fromNode.location, toNode.location, 0.0, 1.0, i/(numDivides+1))
		newNode.offRoadAccess = fromNode.offRoadAccess && toNode.offRoadAccess
		push!(nodes, newNode)

		# add arc
		newArc = deepcopy(tempArc)
		newArc.index = numArcs + i
		newArc.fromNodeIndex = newNode.index
		newArc.toNodeIndex = newNode.index + 1
		push!(arcs, newArc)
	end
	arc.toNodeIndex = numNodes + 1
	arcs[end].toNodeIndex = toNode.index
end

# Divide arcs into multiple arcs (with added nodes) so that the the maximum arc travel time <= maxArcTravelTime.
function graphDivideArcs!(nodes::Vector{Node}, arcs::Vector{Arc}; maxArcTravelTime::Float = Inf)
	if maxArcTravelTime == Inf return end
	@assert(maxArcTravelTime > 0)

	nodePairArcIndex = Dict{Tuple{Int,Int},Int}()
	for arc in arcs
		@assert(!haskey(nodePairArcIndex, (arc.fromNodeIndex, arc.toNodeIndex)), arc) # duplicate arcs not allowed
		nodePairArcIndex[arc.fromNodeIndex, arc.toNodeIndex] = arc.index
	end

	travelTimes = getArcsTravelTimes(arcs) # values before dividing arcs, actual values will change (but not in this var)
	arcsDivided = fill(false, length(arcs))
	for (node1, node2) in keys(nodePairArcIndex)
		if arcsDivided[nodePairArcIndex[node1, node2]] continue end # arc has already been divided
		numNodes = length(nodes) # keep track of number of nodes before current divide
		if !haskey(nodePairArcIndex, (node2, node1)) # no arc in opposite direction
			i = nodePairArcIndex[node1, node2]
			numDivides = ceil(Int, maximum(travelTimes[:,i]) / maxArcTravelTime) - 1
			if numDivides < 1 continue end
			graphDivideArc!(nodes, arcs; arcIndex = i, numDivides = numDivides)
			arcsDivided[i] = true
		else
			# there are two arcs in opposite directions, need to split both simultaneously
			arc1 = nodePairArcIndex[node1, node2]
			arc2 = nodePairArcIndex[node2, node1]
			numDivides = ceil(Int, maximum(travelTimes[:,[arc1,arc2]]) / maxArcTravelTime) - 1 # look at maximum travel time of both arcs
			if numDivides < 1 continue end
			graphDivideArc!(nodes, arcs; arcIndex = arc1, numDivides = numDivides)
			graphDivideArc!(nodes, arcs; arcIndex = arc2, numDivides = numDivides)
			arcsDivided[arc1] = true
			arcsDivided[arc2] = true

			# have created extra nodes when dividing arc2; need to remove these and instead use nodes that were created when dividing arc1
			for i = 1:numDivides
				pop!(nodes)
			end
			for i = 1:numDivides
				arc = arcs[end-i+1]
				arc.fromNodeIndex = numNodes + i
				i == 1 || (arc.toNodeIndex = numNodes + i - 1)
			end
			arcs[arc2].toNodeIndex = numNodes + numDivides
		end
	end

	# check max arc travel time
	travelTimes = getArcsTravelTimes(arcs)
	@assert(maximum(travelTimes) <= maxArcTravelTime)

	checkNodeIndices(nodes)
	checkArcIndices(arcs)
	checkArcNodeIndices(nodes, arcs)
	checkGraphContainsNoDuplicateArcs(arcs)
end

# Tag all nodes and arcs that are used in any shortest path:
# from an origin node to a chosen node,
# from a chosen node to a destination node, or
# from an origin node to a destination node
# This assumes that any duplicate arcs have been removed.
# This version requires LightGraphs.dijkstra_shortest_paths function to allow distmx::Dict, for sake of reducing run time
function graphTagSpElts!(nodes::Vector{Node}, arcs::Vector{Arc};
	spHeader::String = "in_a_sp", chosenNodes::Vector{Int} = Int[], originNodes::Vector{Int} = Int[], destNodes::Vector{Int} = Int[], dijkstraSpDistmxAllowDict::Bool = false)

	@assert(!graphContainsDuplicateArcs(arcs), "Cannot handle duplicate arcs")
	if dijkstraSpDistmxAllowDict @info("Setting `dijkstraSpDistmxAllowDict` = true requires LightGraphs.dijkstra_shortest_paths function to be changed, to allow arg `distmx` to be type Dict") end

	# shorthand
	numArcs = length(arcs)
	numNodes = length(nodes)

	# keep track of all tagged nodes and arcs
	nodeTagged = fill(false, numNodes)
	arcTagged = fill(false, numArcs)

	# to be used, and reset, for each SP tree:
	tagNode = Vector{Bool}(undef, numNodes)
	tagArc = Vector{Bool}(undef, numArcs)

	nodePairArcIndex = Dict{Tuple{Int,Int},Int}()
	revNodePairArcIndex = Dict{Tuple{Int,Int},Int}()
	for arc in arcs
		nodePairArcIndex[arc.fromNodeIndex, arc.toNodeIndex] = arc.index
		revNodePairArcIndex[arc.toNodeIndex, arc.fromNodeIndex] = arc.index
	end

	# # using sparse matrices takes much longer to populate than using dict
	# nodePairArcIndex = spzeros(Int, numNodes, numNodes)
	# for arc in arcs
		# nodePairArcIndex[arc.fromNodeIndex, arc.toNodeIndex] = arc.index
	# end
	# revNodePairArcIndex = transpose(nodePairArcIndex)

	lightGraph = createLightGraph(nodes, arcs)
	revLightGraph = LightGraphs.reverse(lightGraph)

	travelTimes = getArcsTravelTimes(arcs)
	numModes = size(travelTimes,1)

	# create distance matrix, to be reused
	println("Setting arc time matrix / dict")
	if dijkstraSpDistmxAllowDict
		# using this with LightGraphs requires dijkstra_shortest_paths function distmx arg type to be changed
		arcTimes = Dict{Tuple{Int,Int},Float}()
		revArcTimes = Dict{Tuple{Int,Int},Float}()
		for arc in arcs
			if mod(arc.index, 1000) == 0 print("\r\tarc: ", arc.index) end
			arcTimes[arc.fromNodeIndex, arc.toNodeIndex] = nullTime # to be set later
			revArcTimes[arc.toNodeIndex, arc.fromNodeIndex] = nullTime # to be set later
		end
		println()
	else
		# use sparse matrix
		# this is much slower to populate than a dict, but works with unedited version of LightGraphs.dijkstra_shortest_paths
		arcTimes = spzeros(Float, numNodes, numNodes)
		revArcTimes = spzeros(Float, numNodes, numNodes)
		for arc in arcs
			if mod(arc.index, 1000) == 0 print("\r\tarc: ", arc.index) end
			arcTimes[arc.fromNodeIndex, arc.toNodeIndex] = nullTime # to be set later
			revArcTimes[arc.toNodeIndex, arc.fromNodeIndex] = nullTime # to be set later
		end
		println()
	end

	# Tag all nodes and arcs (by setting vars tagNode and tagArc) that are in a shortest path from the originNode to any node in destNodeSets
	# To work with destination node (instead of origin), provide a reversed lightGraph and transposed arcTimes and nodePairArcIndex
	function graphTagSpTreeElts!(lightGraph::LightGraphs.DiGraph{Int}, originNode::Int, destNodeSets::Vector{Vector{Int}}, arcTimes::Union{SparseMatrixCSC{Float,Int},Dict{Tuple{Int,Int},Float}}, nodePairArcIndex::Dict{Tuple{Int,Int},Int})
		# global tagNode, tagArc
		i = originNode # shorthand
		spData = LightGraphs.dijkstra_shortest_paths(lightGraph, i, arcTimes)
		spPreds = spData.parents # spPreds[j] gives index of node before j on shortest path from i to j
		tagNode[:] .= false
		tagNode[i] = true
		tagArc[:] .= false
		for nodeIndices in destNodeSets
			for j in nodeIndices
				while !tagNode[j]
					tagNode[j] = true
					k = j
					j = spPreds[j]
					tagArc[nodePairArcIndex[j,k]] = true
				end
			end
		end
		nodeTagged .|= tagNode
		arcTagged .|= tagArc
		return
	end

	println("Calculating shortest paths")
	for modeIndex in 1:numModes

		println("\tmode: $modeIndex")

		# set distance matrix values
		for arc in arcs
			arcTimes[arc.fromNodeIndex, arc.toNodeIndex] = travelTimes[modeIndex, arc.index]
			revArcTimes[arc.toNodeIndex, arc.fromNodeIndex] = travelTimes[modeIndex, arc.index]
		end

		for (i,originNode) in enumerate(originNodes)
			print("\r\torigin node: $i (of $(length(originNodes)))")
			# tag nodes in a shortest path from node originNode to any chosen/dest nodes
			graphTagSpTreeElts!(lightGraph, originNode, [chosenNodes, destNodes], arcTimes, nodePairArcIndex)
		end
		println()

		for (i,destNode) in enumerate(destNodes)
			print("\r\tdest node: $i (of $(length(destNodes)))")
			# tag nodes in a shortest path from any chosen/origin node to node destNode
			graphTagSpTreeElts!(revLightGraph, destNode, [chosenNodes, originNodes], revArcTimes, revNodePairArcIndex)
		end
		println()

	end

	# use vars `nodeTagged` and `arcTagged` to set value of new field for nodes and arcs
	for i = 1:numNodes
		nodes[i].fields[spHeader] = nodeTagged[i]
	end
	for i = 1:numArcs
		arcs[i].fields[spHeader] = arcTagged[i]
	end

	return
end


text:
module SupplyChainModels

include("EOQ.jl")
include("entities.jl")

a = EOQ.eoq(ct=500, c=50, demand=2000, ce=50*0.25)
a.Q

# Costs
a.HoldingCost
a.OrderingCost
a.TRC()
a.TC()

# Policy
a.policy()
a.policy_for("Q", 500)
a.policy_for("T", 3.0)

for i in (1:10^4)
    c = rand(45:55)
    b = EOQ.eoq(ct=rand(450:550), c=c, demand=rand(1950:2050), ce=c*0.25)
    b.TC()
end

end # Module EOQ

# No puedo correr los códigos acá (fuera del módulo) ---> Investigar y arreglar


text:
abstract type JavascriptSerializer <: IO end

"""
    io_object(io::JavascriptSerializer)
gets the IO object of a JavascriptSerializer
"""
io_object(io::JavascriptSerializer) = io.io

struct JSONSerializer{T <: IO} <: JavascriptSerializer
    io::T
end

# Overload base functions for any object inheriting from JavascriptSerializer
for func in (:write, :print, :println)
    @eval Base.$(func)(io::JavascriptSerializer, args...) = $(func)(io_object(io), args...)
end
Base.print(io::JavascriptSerializer, arg::Union{SubString{String}, String}) = print(io_object(io), arg)
Base.print(io::JavascriptSerializer, arg::Char) = print(io_object(io), arg)

Base.write(io::JSServe.JavascriptSerializer, x::Array) = write(JSServe.io_object(io), x)
Base.write(io::JSServe.JavascriptSerializer, arg::UInt8) = Base.write(io_object(io), arg)

function serialize_websocket(io::IO, message)
    # right now, we simply use
    str = sprint() do io2
        serialize_string(JSONSerializer(io2), message)
    end
    write(io, str)
end

"""
Function used to serialize objects interpolated into a JS string, e.g.
via:
```julia
var = "hi"
js"console.log(\$var)"
```
When sent to the frontend, the final javascript string will be constructed as:
```julia
"console.log(" * serialize_string(var) * ")"
```
"""
serialize_string(@nospecialize(object)) = sprint(io-> serialize_string(io, object))
serialize_string(io::IO, @nospecialize(object)) = JSON3.write(io, object)
serialize_string(io::IO, x::Observable) = print(io, "'", x.id, "'")

# Handle interpolation into javascript
function serialize_string(jsio::JSONSerializer, jss::JSCode)
    # serialize as string
    serialize_string(jsio, serialize_string(jss))
end

serialize_string(io::IO, x::JSString) = print(io, x.source)
function serialize_string(io::IO, jss::JSCode)
    for elem in jss.source
        serialize_string(io, elem)
    end
end
function serialize_string(io::IO, jsss::AbstractVector{JSCode})
    for jss in jsss
        serialize_string(io, jss)
        println(io)
    end
end

struct AsJSON{T}
    x::T
end

serialize_string(io::JSONSerializer, js::AsJSON) = serialize_string(io, js.x)

function serialize_string(io::IO, js::AsJSON)
    serialize_string(JSONSerializer(io), js.x)
end

# Since there is no easy way to mess with JSON3 printer, we print
# Vectors & Dictionaries ourselves, so that we cann apply serialize_string recursively
function serialize_string(io::IO, vector::Union{AbstractVector, Tuple})
    print(io, '[')
    for (i, element) in enumerate(vector)
        serialize_string(io, element)
        (i == length(vector)) || print(io, ',') # dont print for last item
    end
    print(io, ']')
end

function serialize_string(io::IO, dict::Union{NamedTuple, AbstractDict})
    print(io, '{')
    for (i, (key, value)) in enumerate(pairs(dict))
        serialize_string(io, key); print(io, ':')
        serialize_string(io, value)
        (i == length(dict)) || print(io, ',') # dont print for last item
    end
    print(io, '}')
end

# function deserialize_js(::JSObject)
#     # deserialize with the function get_heap_object
#     return js"get_heap_object"
# end


text:
using Dates
using GitHub
using Pkg
using Printf
using RegistryCI
using Test
using TimeZones

const AutoMerge = RegistryCI.AutoMerge

const path = joinpath(DEPOT_PATH[1], "registries", "General")
RegistryCI.test(path)


text:
struct SimpleTauLeaping <: DiffEqBase.DEAlgorithm end
struct RegularSSA <: DiffEqBase.DEAlgorithm end

function DiffEqBase.solve(jump_prob::JumpProblem,alg::SimpleTauLeaping;
                          seed = nothing,
                          dt = error("dt is required for SimpleTauLeaping"))

  @assert isempty(jump_prob.jump_callback.continuous_callbacks)
  @assert isempty(jump_prob.jump_callback.discrete_callbacks)
  prob = jump_prob.prob
  seed === nothing ? rng = Xorshifts.Xoroshiro128Plus() :
                    rng = Xorshifts.Xoroshiro128Plus(seed)

  rj = jump_prob.regular_jump
  rate = rj.rate
  c = rj.c
  dc = zero(rj.dc)
  fill!(dc,0)
  rate_cache = zeros(size(dc,2))

  u0 = copy(prob.u0)
  tspan = prob.tspan
  p = prob.p
  mark = nothing

  n = Int((tspan[2] - tspan[1])/dt) + 1
  u = Vector{typeof(prob.u0)}(undef,n)
  u[1] = u0
  t = [tspan[1] + i*dt for i in 0:n-1]

  counts = zeros(Int,size(dc,2))
  update = similar(u0)

  rj.constant_c && c(dc,u0,p,tspan[1],mark)

  for i in 2:n
      uprev = u[i-1]
      tprev = t[i-1]
      rate(rate_cache,uprev,p,tprev)
      rate_cache .*= dt
      counts .= pois_rand.((rng,), rate_cache)
      !rj.constant_c && c(dc,uprev,p,tprev,mark)
      mul!(update,dc,counts)
      u[i] = uprev .+ update
  end

  sol = DiffEqBase.build_solution(prob,alg,t,u,
                       calculate_error = false,
                       interp = DiffEqBase.ConstantInterpolation(t,u))
end

function DiffEqBase.solve(jump_prob::JumpProblem,alg::RegularSSA)

  @assert isempty(jump_prob.jump_callback.continuous_callbacks)
  @assert isempty(jump_prob.jump_callback.discrete_callbacks)
  prob = jump_prob.prob

  rj = jump_prob.regular_jump
  rate = rj.rate
  c = rj.c
  dc = zero(rj.dc)
  fill!(dc,0)
  rate_cache = zeros(size(dc,2))
  rate_sum = similar(rate_cache)

  u0 = copy(prob.u0)
  tspan = prob.tspan
  p = prob.p
  mark = nothing

  u =[u0]
  t = [tspan[1]]

  rj.constant_c && c(dc,u0,p,tspan[1],mark)

  curt = tspan[1]
  while curt < tspan[2]
      rate(rate_cache,u[end],p,curt)
      cumsum!(rate_sum,rate_cache)
      ttnj = randexp()/rate_sum[end]
      r = rand()
      rate_sum ./= rate_sum[end]
      i = searchsortedfirst(rate_sum,r)
      !rj.constant_c && c(dc,u[end],p,curt,mark)
      # Can instead be dc*[0,0,...,1,...,0,0] if sparse
      # https://github.com/JuliaLang/julia/issues/13438
      unext = u[end] .+ @view dc[:,i]
      curt += ttnj
      push!(t,curt)
      push!(u,unext)
  end

  t[end] = tspan[2]

  sol = DiffEqBase.build_solution(prob,alg,t,u,
                       calculate_error = false,
                       interp = DiffEqBase.ConstantInterpolation(t,u))
end

export SimpleTauLeaping, RegularSSA


text:
abstract type StaticTimeSeries <: TimeSeriesData end

Base.length(ts::StaticTimeSeries) = length(get_data(ts))
get_initial_timestamp(ts::StaticTimeSeries) = TimeSeries.timestamp(get_data(ts))[1]
get_count(ts::StaticTimeSeries) = 1


text:
using Test
using DFTK
using SpecialFunctions

include("testcases.jl")

@testset "Smearing functions" begin
    for m in DFTK.Smearing.smearing_methods
        @test DFTK.Smearing.occupation(m(), -Inf) == 1
        @test DFTK.Smearing.occupation(m(), Inf) == 0
        x = .04
        ε = 1e-8
        @test abs((DFTK.Smearing.occupation(m(), x+ε) - DFTK.Smearing.occupation(m(), x))/ε -
                  DFTK.Smearing.occupation_derivative(m(), x)) < 1e-4

        # entropy functions should satisfy s' = x f'
        sprime = (DFTK.Smearing.entropy(m(), x+ε) - DFTK.Smearing.entropy(m(), x))/ε
        fprime = (DFTK.Smearing.occupation(m(), x+ε) - DFTK.Smearing.occupation(m(), x))/ε
        @test abs(sprime - x*fprime) < 1e-4
    end
end

@testset "Smearing for insulators" begin
    Ecut = 5
    n_bands = 10
    fft_size = [15, 15, 15]

    # Emulate an insulator ... prepare energy levels
    energies = [zeros(n_bands) for k in silicon.kcoords]
    n_occ = div(silicon.n_electrons, 2)
    n_k = length(silicon.kcoords)
    for ik in 1:n_k
        energies[ik] = sort(rand(n_bands))
        energies[ik][n_occ+1:end] .+= 2
    end
    εHOMO = maximum(energies[ik][n_occ] for ik in 1:n_k)
    εLUMO = minimum(energies[ik][n_occ + 1] for ik in 1:n_k)

    # Occupation for zero temperature
    model = Model(silicon.lattice; n_electrons=silicon.n_electrons, temperature=0.0, smearing=nothing)
    basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.ksymops; fft_size=fft_size)
    occupation0, εF0 = find_occupation_bandgap(basis, energies)
    @test εHOMO < εF0 < εLUMO
    @test sum(basis.kweights .* sum.(occupation0)) ≈ model.n_electrons

    # See that the electron count still works if we add temperature
    Ts = (0, 1e-6, .1, 1.0)
    for T in Ts, meth in DFTK.Smearing.smearing_methods
        model = Model(silicon.lattice; n_electrons=silicon.n_electrons, temperature=T, smearing=meth())
        basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.ksymops; fft_size=fft_size)
        occs, _ = find_occupation(basis, energies)
        @test sum(basis.kweights .* sum.(occs)) ≈ model.n_electrons
    end

    # See that the occupation is largely uneffected with only a bit of temperature
    Ts = (0, 1e-6, 1e-4)
    for T in Ts, meth in DFTK.Smearing.smearing_methods
        model = Model(silicon.lattice; n_electrons=silicon.n_electrons, temperature=T, smearing=meth())
        basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.ksymops; fft_size=fft_size)
        occupation, _= find_occupation(basis, energies)

        for ik in 1:n_k
            @test all(isapprox.(occupation[ik], occupation0[ik], atol=1e-2))
        end
    end
end

@testset "Smearing for metals" begin
    testcase = magnesium
    Ecut = 5
    fft_size = [15, 15, 15]
    kgrid_size = [2, 3, 4]

    # Emulate a metal ...
    energies = [[-0.08063210585291,  0.11227915155236, 0.13057816014162, 0.57672256037074],
                [ 0.09509047528102,  0.09538152469111, 0.27197836572013, 0.28750689088845],
                [-0.00144586520885,  0.18640677556553, 0.19603060374450, 0.24422060327989],
                [ 0.05693643182609,  0.16919740718547, 0.24190245274401, 0.25674283154835],
                [-0.06756541677784,  0.03381889875058, 0.23162853469956, 0.50981867707851],
                [ 0.10685980948954,  0.10728887405642, 0.20784971952147, 0.20786603845828],
                [ 0.01122399002894,  0.11011069317735, 0.24016826005369, 0.30770620467001],
                [ 0.06925846412968,  0.16087157153058, 0.19146746736359, 0.27463770659603],
                [-0.02937886574534, -0.02937886574483, 0.36206906745747, 0.36206906745749],
                [ 0.13314087354890,  0.13314087354890, 0.15834732772541, 0.15834732772541],
                [ 0.04869672986772,  0.04869672986772, 0.27749728805752, 0.27749728805768],
                [ 0.10585630776222,  0.10585630776223, 0.22191839818805, 0.22191839818822]]

    spec = ElementPsp(testcase.atnum, psp=load_psp(testcase.psp))
    symops = DFTK.symmetry_operations(testcase.lattice,
                                      [spec => testcase.positions])
    kcoords, ksymops = bzmesh_ir_wedge(kgrid_size, symops)

    n_bands = length(energies[1])
    n_k = length(kcoords)
    @assert n_k == length(energies)

    parameters = (
        (DFTK.Smearing.FermiDirac,         0.01, 0.17251898225370),
        (DFTK.Smearing.FermiDirac,         0.02, 0.17020763046058),
        (DFTK.Smearing.FermiDirac,         0.03, 0.16865552281082),
        (DFTK.Smearing.MethfesselPaxton1, 0.01, 0.16917895217084),
        (DFTK.Smearing.MethfesselPaxton1, 0.02, 0.17350869020891),
        (DFTK.Smearing.MethfesselPaxton1, 0.03, 0.17395190342809),
    )

    for (meth, Tsmear, εF_ref) in parameters
        model = Model(silicon.lattice, n_electrons=testcase.n_electrons;
                      temperature=Tsmear, smearing=meth())
        basis = PlaneWaveBasis(model, Ecut, kcoords, ksymops; fft_size=fft_size)
        occupation, εF = find_occupation(basis, energies)

        @test sum(basis.kweights .* sum.(occupation)) ≈ model.n_electrons
        @test εF ≈ εF_ref
    end
end


text:
#=
    Main test suite, which runs regression & comparison tests.
=#

using Test

@testset verbose=true "All c_sw tests"  begin
    @testset verbose=true "Regression" begin
        @test include("regression_tests.jl")
    end
    @testset verbose=true "Comparison " begin
        @test include("compare_tests.jl")
    end
end


text:
#using Pkg
#pkg"activate ."

#push!(LOAD_PATH, pwd(),"C")

include("Gopts.jl")
using .Gopts
include("Juleva.jl")
using .Juleva

# Avoid ambiguity
fpi = convert(Float64,π)

# Parameters
x = Array{Float64,1}(undef,3)
y = Array{Float64,1}(undef,3)

xs = [0.; 0.; 0.]
xl = [0.; 0.; 0.]
xu = [fpi; fpi; fpi]

function integrand(x::Ptr{Float64},dim::Ptr{Float64},params::Ptr{Float64})
  A = 1.0 / (pi * pi * pi)
  return 3*A / (1.0 - cos(unsafe_load(x,1))*cos(unsafe_load(x,2))*cos(unsafe_load(x,3)))::Cdouble
end

# Init options
g_opts = GOpts()
g_opts.algo = 2
g_opts.dim = 3
g_opts.iter = 10
g_opts.retRes = pointer(x)
g_opts.retErr = pointer(y)
g_opts.x_s = pointer(xs)
g_opts.x_l = pointer(xl)
g_opts.x_u = pointer(xu)
g_opts.func = @cfunction(integrand,Cdouble,(Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}))

println(g_opts.dim)

res = optimize(g_opts, [0.1,0.2,0.3])

println(res)

max_stars_repo_name: denisbertini/juleva

max_issues_repo_name: denisbertini/juleva

max_forks_repo_name: denisbertini/juleva


text:
function part1()
    tot = 0
    @inbounds for line in eachline("../data/input8.txt")
        _, final = split(line, " | ")
        segs = split(final)
        tot += count(x -> length(x) in (2, 4, 7, 3), segs)
    end
    return tot
end

function part2()
    tot = 0
    sum2digit = Dict(42 => 0, 17 => 1, 34 => 2, 39 => 3, 30 => 4, 37 => 5, 41 => 6, 25 => 7, 49 => 8, 45 => 9)

    @inbounds for line in eachline("../data/input8.txt")
        initial, final = split(line, " | ")

        freqs = Dict()
        for c in "abcdefg"
            freqs[c] = count(==(c), initial)
        end

        segs_o = split(final)
        for (i, seg) in enumerate(segs_o)
            digit = sum2digit[sum(freqs[c] for c in seg)]
            tot += digit * 10^(4-i)
        end
    end
    return tot
end

max_stars_repo_name: lucaferranti/AoC

max_issues_repo_name: lucaferranti/AoC

max_forks_repo_name: lucaferranti/AoC


text:
module DOS

using Interpolations
using DelimitedFiles
using Statistics
using Zygote

export DOSData, get_dos

"""
    DOSData(dos_file; Ef=0, cut_energy=false)

Struct for storing density of states (DOS) information. Can be constructed from a CSV file.

# Fields
* interp_func: function that returns interpolated values of DOS for any energy between E_min and E_max
* average_value: Average value of DOS over this interval
* E_min, E_max: bounds over which this DOS is defined
"""
struct DOSData
    interp_func
    average_value::Float64
    E_min::Float64
    E_max::Float64
end

DOSData(dos_file; Ef=0, cut_energy=false) = DOSData(get_dos(dos_file; Ef, cut_energy)...)

# pretty printing
Base.show(io::IO, dd::DOSData) = print(io, "DOSData: avg value $(round(dd.average_value, sigdigits=3)) from energy $(round(dd.E_min, sigdigits=2)) to $(round(dd.E_max, sigdigits=2))")

(dd::DOSData)(E::Real) = dd.interp_func(E)

"""
    get_dos(dos_file; Ef=0, cut_energy=false)

Retrieve density of states data from a file at path `dos_file`, assuming that the Fermi level is at energy `Ef`.

Returns a callable interpolated DOS, the average value of DOS, and the lower and upper bound energies.

# Notes
* `dos_file` is assumed to contain two colums, energy and DOS value
* If `cut_energy==true`, return data in a symmetric range of energies around `Ef` (this is sometimes useful for nice-looking plots). Note that the code presumes the cutoff for this comes from the upper bound. This would not be hard to fix, I just have not done so as yet.
"""
function get_dos(dos_file::String; kwargs...)
    dos_data = Zygote.ignore() do
        x = readdlm(dos_file, Float64, skipstart=1)
        x
    end
    get_dos(dos_data; kwargs...)
end
function get_dos(dd::Matrix; Ef=0, cut_energy=false)
    dos_data = Zygote.ignore() do
        x = dd
        # recenter so Ef=0
        x[:,1] = x[:,1] .- Ef
        x
    end
    E_max = dos_data[end,1]
    if cut_energy #TODO: fix this so that it works if magnitude of upper bound is smaller
        # for now, we'll pick a symmetric energy range about 0 (i.e. the new Ef)
        len_keep = sum(dos_data[:,1] .> -E_max)
        # cut off
        dos_data = dos_data[end-len_keep:end,:]
    end
    E_min = dos_data[1,1]
    average_dos = mean(dos_data[:,2]) # for whole structure
    # interpolate
    E_step = mean(dos_data[2:end,1].-dos_data[1:end-1,1])
    #dos_interp = scale(interpolate(dos_data[:,2], BSpline(Linear())), range(E_min, E_max+0.0001, step=E_step))
    dos_interp = scale(interpolate(dos_data[:,2], BSpline(Linear())), range(E_min, E_max, step=E_step))
    return dos_interp, average_dos, E_min, E_max
end

end


text:
mutable struct SixelAllocator
    ptr::Ptr{C.sixel_allocator_t}
    function SixelAllocator()
        # libsixel checks if the pointer is NULL during initialization so we use a non-zero
        # "invalid" address value to bypass the NULL check
        ptr_ref = Ref{Ptr{C.sixel_allocator_t}}(1)
        status = C.sixel_allocator_new(ptr_ref, C_NULL, C_NULL, C_NULL, C_NULL)
        check_status(status)
        new(ptr_ref[])
    end
end

mutable struct SixelOutput
    ptr::Ptr{C.sixel_output_t}
    function SixelOutput(fn_write, priv; allocator::SixelAllocator=SixelAllocator())
        # libsixel checks if the pointer is NULL during initialization so we use a non-zero
        # "invalid" address value to bypass the NULL check
        ptr_ref = Ref{Ptr{C.sixel_output_t}}(1)
        status = C.sixel_output_new(ptr_ref, fn_write, priv, allocator.ptr)
        check_status(status)
        new(ptr_ref[])
    end
end

mutable struct SixelDither
    ptr::Ptr{C.sixel_dither_t}
    function SixelDither(
            data, width, height, pixelformat, quality_mode,
            method_for_largest=C.LARGE_NORM,
            method_for_rep=C.REP_CENTER_BOX;
            allocator::SixelAllocator=SixelAllocator())
        # libsixel checks if the pointer is NULL during initialization so we use a non-zero
        # "invalid" address value to bypass the NULL check
        ptr_ref = Ref{Ptr{C.sixel_dither_t}}(1)
        ncolors = -1 # libsixel internally sets -1 to length(palette)
        status = C.sixel_dither_new(ptr_ref, ncolors, allocator.ptr)
        check_status(status)
        ptr = ptr_ref[]
        C.sixel_dither_initialize(ptr,
            data,
            Int32(width),
            Int32(height),
            Int32(pixelformat),
            Int32(method_for_largest),
            Int32(method_for_rep),
            Int32(quality_mode)
        )
        new(ptr)
    end
end

check_status(status) = status == 0 || error("libsixel error code: $status")


text:
## ----------  Parameter Transformations

"""
	ParamTransformation

Abstract type for transforming parameters into bounded values (guesses) and reverse.

Define a concrete type with its own parameters and method
```julia
    transform_param(tr :: ParamTransformation, p :: Param)
```
"""
abstract type ParamTransformation{F1 <: Real} end

lb(::ParamTransformation{F1}) where F1 = F1(TransformationLb);
ub(::ParamTransformation{F1}) where F1 = F1(TransformationUb);

"""
	LinearTransformation

Default linear transformation into default interval [1, 2].
`F1` encodes the element type. The bounds are always scalar.
"""
struct LinearTransformation{F1} <: ParamTransformation{F1}
    # lb :: F1 = one(F1)
    # ub :: F1 = F1(2.0)
end

# ----------------


text:
using Test
using TimeMachine

@testset "dtw 1" begin
    # dist = Distances.Euclidean()
    local cost = dtw([1, 2, 4, 1, 2], [1, 2, 4, 1, 2], (a, b) -> abs(a - b))
    @test 0.0 == cost

    cost = dtw([2, 6, 5, 3, 2], [1, 2, 4, 1, 2], (a, b) -> abs(a - b))

    @test 1 == cost
end

@testset "dtw 2" begin
    # dist = Distances.SqEuclidean()
    local a = [1, 1, 1, 2, 4, 6, 5, 5, 5, 4, 4, 3, 1, 1, 1]
    local b = [1, 1, 2, 4, 6, 6, 6, 5, 4, 4, 4, 3, 3, 3, 1]
    local cost = dtw(a, b, (a, b) -> (a - b)^2)

    @test 0.0 == cost

    a[end] += 2
    cost = dtw(a, b, (a, b) -> (a - b)^2)
    @test (2.0 / 11.0) == cost

    a = collect(1:10)
    b = a .+ 1
    cost = dtw(a, b, (a, b) -> (a - b)^2)
    @test (2.0 / 11.0) == cost
end


text:
#!/usr/bin/env julia

using SphericalFourierBesselDecompositions
SFB = SphericalFourierBesselDecompositions

using Test
using LinearAlgebra

@testset "Mixing matrices" begin
    @testset "No window" begin
        # Note: inaccuracies in this test are dominated by inaccuracies in healpix.
        rmin = 500.0
        rmax = 1000.0
        amodes = SFB.AnlmModes(3, 5, rmin, rmax)
        wmodes = SFB.ConfigurationSpaceModes(rmin, rmax, 1000, amodes.nside)
        win = SFB.make_window(wmodes, :separable)

        # wmix
        wmix = SFB.calc_wmix(win, wmodes, amodes)
        diff = wmix - I
        for i′=1:size(diff,2), i=1:size(diff,1)
            @test isapprox(diff[i,i′], 0, atol=1e-2)
        end
        @test isapprox(wmix, I, atol=1e-2)

        # wmix_negm
        wmix_negm = SFB.calc_wmix(win, wmodes, amodes, neg_m=true)
        for i′=1:size(diff,2), i=1:size(diff,1)
            n, l, m = SFB.getnlm(amodes, i)
            n′, l′, m′ = SFB.getnlm(amodes, i′)
            correct = (n==n′ && l==l′ && m==-m′) ? 1 : 0
            @test isapprox(wmix_negm[i,i′], correct, atol=1e-2)
        end

        # M
        cmodes = SFB.ClnnModes(amodes, Δnmax=1)
        mmix1 = SFB.power_win_mix(wmix, wmix_negm, cmodes)
        @test isapprox(mmix1, I, atol=1e-2)
        M = SFB.power_win_mix(win, wmodes, cmodes)
        @test M ≈ mmix1
        @test isapprox(M, I, atol=1e-2)

        # CWlnn
        pk(k) = 1e4 * (k/1e-2)^(-3.1)
        Clnn = SFB.gen_Clnn_theory(pk, cmodes)
        CnlmNLM = SFB.Clnn2CnlmNLM(Clnn, cmodes)
        CWlnn1 = SFB.sum_m_lmeqLM(wmix * CnlmNLM * wmix', cmodes)
        CWlnn2 = M * Clnn
        @test CWlnn1 ≈ CWlnn2

        # Nshot
        nbar = 3e-4
        Nshotobs = SFB.win_lnn(win, wmodes, cmodes) ./ nbar
        Nshotobs2 = SFB.sum_m_lmeqLM(wmix, cmodes) ./ nbar
        @test isapprox(Nshotobs, Nshotobs2, rtol=1e-2)
    end

    # Inhomogeneous masks
    @testset "Inhomogeneous Window sep & insep" begin
        rmin = 500.0
        rmax = 1000.0
        amodes = SFB.AnlmModes(2, 5, rmin, rmax)
        wmodes = SFB.ConfigurationSpaceModes(rmin, rmax, 1000, amodes.nside)

        win1 = SFB.make_window(wmodes, :separable)
        @. win1.mask = rand()
        sidx = (1:length(win1.mask))[win1.mask .> 0]
        sidx = 1:(length(sidx) ÷ 2)
        @. win1.mask[sidx] = 0.5 * win1.mask[sidx]

        win2 = [win1[i,j] for i=1:size(win1,1), j=1:size(win1,2)]
        @show sum(win1.mask) ./ length(win1.mask)
        @show extrema(win1.mask)
        @show typeof(win1) typeof(win2)
        @show size(win1) size(win2)

        # wmix
        wmix1 = SFB.calc_wmix(win1, wmodes, amodes)
        wmix1_negm = SFB.calc_wmix(win1, wmodes, amodes, neg_m=true)
        wmix2 = SFB.calc_wmix(win2, wmodes, amodes)
        wmix2_negm = SFB.calc_wmix(win2, wmodes, amodes, neg_m=true)
        @show typeof(wmix1) typeof(wmix2)
        @test wmix1 == wmix2

        # M
        cmodes = SFB.ClnnModes(amodes, Δnmax=1)
        M = SFB.power_win_mix(wmix1, wmix1_negm, cmodes)
        M1 = SFB.power_win_mix(win1, wmodes, cmodes)
        M2 = SFB.power_win_mix(win2, wmodes, cmodes)
        @show extrema(M .- M1)
        @show extrema(M .- M2)
        @test M ≈ M1
        @test M ≈ M2

        # CWlnn
        pk(k) = 1e4 * (k/1e-2)^(-3.1)
        Clnn = SFB.gen_Clnn_theory(pk, cmodes)
        CnlmNLM = SFB.Clnn2CnlmNLM(Clnn, cmodes)
        CWlnn1 = SFB.sum_m_lmeqLM(wmix1 * CnlmNLM * wmix1', cmodes)
        CWlnn2 = M * Clnn
        @show extrema((CWlnn1 .- CWlnn2) ./ CWlnn1)
        @test CWlnn1 ≈ CWlnn2 rtol=1e-2

        # Nshot
        nbar = 3e-4
        Nshotobs1 = SFB.win_lnn(win1, wmodes, cmodes) ./ nbar
        Nshotobs2 = SFB.win_lnn(win2, wmodes, cmodes) ./ nbar
        @test Nshotobs1 ≈ Nshotobs2

        # Bandpower binning
        fsky = sum(win1[1,:]) / size(win1,2)
        Δℓ = round(Int, 1 / fsky)
        w̃, v = SFB.bandpower_binning_weights(cmodes; Δℓ=Δℓ, Δn=1)
        bcmodes = SFB.ClnnBinnedModes(w̃, v, cmodes)
        Nmix1 = SFB.power_win_mix(win1, w̃, v, wmodes, bcmodes)
        Nmix2 = SFB.power_win_mix(win2, w̃, v, wmodes, bcmodes)
        @show extrema(Nmix1 .- Nmix2)
        @test Nmix1 ≈ Nmix2
        w̃M1 = SFB.power_win_mix(win1, w̃, I, wmodes, bcmodes)
        Mv1 = SFB.power_win_mix(win1, I, v, wmodes, bcmodes)
        w̃M2 = SFB.power_win_mix(win2, w̃, I, wmodes, bcmodes)
        Mv2 = SFB.power_win_mix(win2, I, v, wmodes, bcmodes)
        @show extrema(w̃M1 .- w̃M2)
        @test w̃M1 ≈ w̃M2
        @test inv(Nmix1) * w̃M1 ≈ inv(Nmix2) * w̃M2
        @test Mv1 * inv(Nmix1) ≈ Mv2 * inv(Nmix2)
    end

    # complex windows
    win_descriptions = [
                        (:ang_75,),
                        (:ang_75, :radial),
                        (:ang_half, :radial),
                        (:ang_quarter, :radial),
                        (:separable, :ang_75,),
                        (:separable, :ang_75, :radial),
                        (:separable, :ang_half, :radial),
                        (:separable, :ang_quarter, :radial),
                      ]
    @testset "Window $(win_features...)" for win_features in win_descriptions
        rmin = 500.0
        rmax = 1000.0
        amodes = SFB.AnlmModes(2, 5, rmin, rmax)
        wmodes = SFB.ConfigurationSpaceModes(rmin, rmax, 1000, amodes.nside)
        win = SFB.make_window(wmodes, win_features...)

        # wmix
        wmix = SFB.calc_wmix(win, wmodes, amodes)
        wmix_negm = SFB.calc_wmix(win, wmodes, amodes, neg_m=true)

        # M
        cmodes = SFB.ClnnModes(amodes, Δnmax=1)
        mmix1 = SFB.power_win_mix(wmix, wmix_negm, cmodes)
        M = SFB.power_win_mix(win, wmodes, cmodes)
        @test M ≈ mmix1

        # CWlnn
        pk(k) = 1e4 * (k/1e-2)^(-3.1)
        Clnn = SFB.gen_Clnn_theory(pk, cmodes)
        CnlmNLM = SFB.Clnn2CnlmNLM(Clnn, cmodes)
        CWlnn1 = SFB.sum_m_lmeqLM(wmix * CnlmNLM * wmix', cmodes)
        CWlnn2 = M * Clnn
        @test CWlnn1 ≈ CWlnn2

        # Nshot
        nbar = 3e-4
        Nshotobs = SFB.win_lnn(win, wmodes, cmodes) ./ nbar
        Nshotobs2 = SFB.sum_m_lmeqLM(wmix, cmodes) ./ nbar
        @test isapprox(Nshotobs, Nshotobs2, rtol=1e-2)

        # Bandpower binning basics
        fsky = sum(win[1,:]) / size(win,2)
        Δℓ = round(Int, 1 / fsky)
        w̃, v = SFB.bandpower_binning_weights(cmodes; Δℓ=Δℓ, Δn=1)
        N = w̃ * M * v
        w = inv(N) * w̃ * M
        ṽ = M * v * inv(N)
        @test w * v ≈ I  # these really just test linear algebra
        @test w̃ * ṽ ≈ I  # these really just test linear algebra

        # Bandpower binning optimizations
        bcmodes = SFB.ClnnBinnedModes(w̃, v, cmodes)
        Nmix = SFB.power_win_mix(win, w̃, v, wmodes, bcmodes)
        @test Nmix ≈ N
        w̃M = SFB.power_win_mix(win, w̃, I, wmodes, bcmodes)
        Mv = SFB.power_win_mix(win, I, v, wmodes, bcmodes)
        @test w̃M ≈ w̃ * M
        @test inv(Nmix) * w̃M ≈ w
        @test Mv * inv(Nmix) ≈ ṽ
    end
end

# vim: set sw=4 et sts=4 :


text:
@testset "rank/file attack" begin
    occ = [0,1,0,0,0,0,0,0,
           0,0,0,0,0,1,0,0,
           0,1,0,0,0,1,0,0,
           0,0,0,1,0,0,0,0,
           0,0,0,0,1,0,0,0,
           0,0,0,1,0,0,0,0,
           0,0,0,0,0,1,0,0,
           0,1,0,0,0,0,1,0]

    att = [0,0,1,0,0,0,0,0,
           0,0,1,0,0,0,0,0,
           0,0,1,0,0,0,0,0,
           0,0,1,0,0,0,0,0,
           1,1,0,1,0,0,0,0,
           0,0,1,0,0,0,0,0,
           0,0,1,0,0,0,0,0,
           0,0,1,0,0,0,0,0]

    edg = [0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,1,0,0,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0,
           0,0,0,0,0,0,0,0]
    occ_uint = Bobby.cvt_to_uint(BitArray(occ))
    ui = Bobby.INT2UINT[35]
    m, e = Bobby.orthogonal_attack(occ_uint, ui)
    @test all(Int.(Bobby.cvt_to_bitarray(m)) .== att)
    @test all(Int.(Bobby.cvt_to_bitarray(e)) .== edg)
end

@testset "rook fen" begin
    @test test_fen("k2r1r2/8/8/8/8/p3p2p/P3P2P/R3K3 w KQkq - 0 1", 1, [3])
    @test test_fen("k2r1r2/8/8/8/8/p3p2p/P3P2P/4K2R w KQkq - 0 1", 1, [2])
    @test test_fen("k2r1r2/8/8/8/8/p3p2p/P3P2P/R3K2R w KQkq - 0 1", 1, [5])
end


text:
@doc raw"""
`SDE`: Stratonovich Stochastic Differential Equation

Defines a stochastic differential initial value problem
```math
\begin{align*}
\dq (t) &= v(t, q(t)) \, dt + B(t, q(t)) \circ dW , & q(t_{0}) &= q_{0} ,
\end{align*}
```
with drift vector field ``v``, diffusion matrix ``B``,
initial conditions ``q_{0}``, the dynamical variable ``q``
taking values in ``\mathbb{R}^{d}``, and the m-dimensional Wiener process W

### Fields

* `d`:  dimension of dynamical variable ``q`` and the vector field ``v``
* `m`:  dimension of the Wiener process
* `n`:  number of initial conditions
* `ns`: number of sample paths
* `v`:  function computing the deterministic vector field
* `B`:  function computing the d x m diffusion matrix
* `t₀`: initial time
* `q₀`: initial condition for dynamical variable ``q`` (may be a random variable itself)

The functions `v` and `B`, providing the drift vector field and diffusion matrix,
`v(t, q, v)` and `B(t, q, B; col=0)`, where `t` is the current time, `q` is the
current solution vector, and `v` and `B` are the variables which hold the result
of evaluating the vector field ``v`` and the matrix ``B`` on `t` and `q` (if col==0),
or the column col of the matrix B (if col>0).

### Example

```julia
    function v(λ, t, q, v)
        v[1] = λ*q[1]
        v[2] = λ*q[2]
    end

    function B(μ, t, q, B; col=0)
        if col==0 #whole matrix
            B[1,1] = μ*q[1]
            B[2,1] = μ*q[2]
        elseif col==1
            #just first column
        end
    end

    t₀ = 0.
    q₀ = [1., 1.]
    λ  = 2.
    μ  = 1.

    v_sde = (t, q, v) -> v(λ, t, q, v)
    B_sde = (t, q, B) -> B(μ, t, q, B)

    sde = SDE(v_sde, B_sde, t₀, q₀)
```
"""
struct SDE{dType <: Number, tType <: Number, vType <: Function, BType <: Function, N} <: Equation{dType, tType}
    d::Int
    m::Int
    n::Int
    ns::Int
    v::vType
    B::BType
    t₀::tType
    q₀::Array{dType, N}           #Initial condition: N=1 - single deterministic, N=2 - single random or multiple deterministic, N=3 - multiple random
    periodicity::Vector{dType}

    function SDE{dType,tType,vType,BType,N}(d, m, n, ns, v, B, t₀, q₀; periodicity=[]) where {dType <: Number, tType <: Number, vType <: Function, BType <: Function, N}

        @assert dType == eltype(q₀)
        @assert tType == typeof(t₀)
        @assert ndims(q₀) == N
        @assert d == size(q₀,1)

        if ns==1 && n==1
            # single sample path and single initial condition, therefore N=1
            @assert N == 1
        elseif ns>1 && n==1
            # single initial condition, but may be random (N=2) or deterministic (N=1)
            @assert N ∈ (1,2)
            if N==2
                @assert ns == size(q₀,2)
            end
        elseif ns==1 && n>1
            # multiple  deterministic initial conditions, so N=2
            @assert N == 2
            @assert n == size(q₀,2)
        elseif ns>1 && n>1
            # either multiple random initial conditions (N=3),
            # or multiple deterministic initial conditions (N=2)
            @assert N ∈ (2,3)

            if N==2
                @assert n == size(q₀,2)
            else
                @assert ns == size(q₀,2)
                @assert n == size(q₀,3)
            end
        end

        if !(length(periodicity) == d)
            periodicity = zeros(dType, d)
        end

        new(d, m, n, ns, v, B, t₀, q₀, periodicity)
    end
end

function SDE(m::Int, ns::Int, v::VT, B::BT, t₀::TT, q₀::DenseArray{DT,1}; periodicity=[]) where {DT,TT,VT,BT}
    # A 1D array q₀ contains a single deterministic initial condition, so n=1, but we still need to specify
    # the number of sample paths ns
    SDE{DT, TT, VT, BT, 1}(size(q₀, 1), m, 1, ns, v, B, t₀, q₀, periodicity=periodicity)
end

# A 2-dimensional matrix q0 can represent a single random initial condition with ns>1 and n=1,
# or a set of deterministic initial conditions with n>1 (for which we can have both ns=1 and ns>1)
# The function below assumes q0 to represent a single random initial condition (n=1, ns=size(q₀, 2))
function SDE(m::Int, v::VT, B::BT, t₀::TT, q₀::DenseArray{DT,2}; periodicity=[]) where {DT,TT,VT,BT}
    SDE{DT, TT, VT, BT, 2}(size(q₀, 1), m, 1, size(q₀, 2), v, B, t₀, q₀, periodicity=periodicity)
end

# On the other hand, the function below assumes q₀ represents multiple deterministic initial conditions
# (n=size(q₀, 2)), but these initial conditions may be run an arbitrary number ns of sample paths, so ns has to be explicitly specified
function SDE(m::Int, ns::Int, v::VT, B::BT, t₀::TT, q₀::DenseArray{DT,2}; periodicity=[]) where {DT,TT,VT,BT}
    SDE{DT, TT, VT, BT, 2}(size(q₀, 1), m, size(q₀, 2), ns, v, B, t₀, q₀, periodicity=periodicity)
end

# # OLD FUNCTION FOR A 2D matrix q₀
# # The argument IC specifies whether there are multiple deterministic initial conditions
# # (IC=true, so n>1) or a single random one (default IC=false, so n=1)
# function SDE(m::Int, v::VT, B::BT, t₀::TT, q₀::DenseArray{DT,2}; periodicity=[], IC=false) where {DT,TT,VT,BT}
#     if IC==true
#         SDE{DT, TT, VT, BT, 2}(size(q₀, 1), m, size(q₀, 2), 1, v, B, t₀, q₀, periodicity=periodicity)
#     else
#         SDE{DT, TT, VT, BT, 2}(size(q₀, 1), m, 1, size(q₀, 2), v, B, t₀, q₀, periodicity=periodicity)
#     end
# end

function SDE(m::Int, v::VT, B::BT, t₀::TT, q₀::DenseArray{DT,3}; periodicity=[]) where {DT,TT,VT,BT}
    # A 3D array q₀ contains multiple random initial condition, so n=size(q₀,3) and ns=size(q₀,2)
    SDE{DT, TT, VT, BT, 3}(size(q₀, 1), m, size(q₀,3), size(q₀,2), v, B, t₀, q₀, periodicity=periodicity)
end

function SDE(m::Int, ns::Int, v::VT, B::BT, q₀::DenseArray{DT,1}; periodicity=[]) where {DT,VT,BT}
    SDE(m, ns, v, B, zero(DT), q₀, periodicity=periodicity)
end

# Assumes q0 represents a single random initial condition (n=1, ns=size(q₀, 2))
function SDE(m::Int, v::VT, B::BT, q₀::DenseArray{DT,2}; periodicity=[]) where {DT,VT,BT}
    SDE(m, v, B, zero(DT), q₀, periodicity=periodicity)
end

# Assumes q₀ represents multiple deterministic initial conditions (n=size(q₀, 2))
function SDE(m::Int, ns::Int, v::VT, B::BT, q₀::DenseArray{DT,2}; periodicity=[]) where {DT,VT,BT}
    SDE(m, ns, v, B, zero(DT), q₀, periodicity=periodicity)
end

function SDE(m::Int, v::VT, B::BT, q₀::DenseArray{DT,3}; periodicity=[]) where {DT,VT,BT}
    SDE(m, v, B, zero(DT), q₀, periodicity=periodicity)
end

Base.hash(sde::SDE, h::UInt) = hash(sde.d, hash(sde.m, hash(sde.n, hash(sde.ns, hash(sde.v, hash(sde.B, hash(sde.t₀, hash(sde.q₀, hash(sde.periodicity, h)))))))))

Base.:(==)(sde1::SDE, sde2::SDE) = (
                                sde1.d == sde2.d
                             && sde1.m == sde2.m
                             && sde1.n == sde2.n
                             && sde1.ns == sde2.ns
                             && sde1.v == sde2.v
                             && sde1.B == sde2.B
                             && sde1.t₀ == sde2.t₀
                             && sde1.q₀ == sde2.q₀
                             && sde1.periodicity == sde2.periodicity)

function Base.similar(sde::SDE{DT,TT,VT,UT}, q₀::DenseArray{DT,1}) where {DT, TT, VT, UT}
    similar(sde, sde.t₀, q₀)
end

# Assumes q0 represents a single random initial condition (n=1, ns=size(q₀, 2))
function Base.similar(sde::SDE{DT,TT,VT,UT}, q₀::DenseArray{DT,2}) where {DT, TT, VT, UT}
    similar(sde, sde.t₀, q₀)
end

function Base.similar(sde::SDE{DT,TT,VT,UT}, q₀::DenseArray{DT,3}) where {DT, TT, VT, UT}
    similar(sde, sde.t₀, q₀)
end

function Base.similar(sde::SDE{DT,TT,VT,BT}, t₀::TT, q₀::DenseArray{DT,1}) where {DT, TT, VT, BT}
    @assert sde.d == size(q₀,1)
    SDE(sde.m, sde.ns, sde.v, sde.B, t₀, q₀, periodicity=sde.periodicity)
end

# Assumes q0 represents a single random initial condition (n=1, ns=size(q₀, 2))
function Base.similar(sde::SDE{DT,TT,VT,BT}, t₀::TT, q₀::DenseArray{DT,2}) where {DT, TT, VT, BT}
    @assert sde.d == size(q₀,1)
    SDE(sde.m, sde.ns, sde.v, sde.B, t₀, q₀, periodicity=sde.periodicity)
end

function Base.similar(sde::SDE{DT,TT,VT,BT}, t₀::TT, q₀::DenseArray{DT,3}) where {DT, TT, VT, BT}
    @assert sde.d == size(q₀,1)
    SDE(sde.m, sde.ns, sde.v, sde.B, t₀, q₀, periodicity=sde.periodicity)
end

function Base.similar(sde::SDE{DT,TT,VT,UT}, q₀::DenseArray{DT,1}, ns::Int) where {DT, TT, VT, UT}
    similar(sde, sde.t₀, q₀, ns)
end

# Assumes q₀ represents multiple deterministic initial conditions (n=size(q₀, 2))
function Base.similar(sde::SDE{DT,TT,VT,UT}, q₀::DenseArray{DT,2}, ns::Int) where {DT, TT, VT, UT}
    similar(sde, sde.t₀, q₀, ns)
end

function Base.similar(sde::SDE{DT,TT,VT,BT}, t₀::TT, q₀::DenseArray{DT,1}, ns::Int) where {DT, TT, VT, BT}
    @assert sde.d == size(q₀,1)
    SDE(sde.m, ns, sde.v, sde.B, t₀, q₀, periodicity=sde.periodicity)
end

# Assumes q₀ represents multiple deterministic initial conditions (n=size(q₀, 2))
function Base.similar(sde::SDE{DT,TT,VT,BT}, t₀::TT, q₀::DenseArray{DT,2}, ns::Int) where {DT, TT, VT, BT}
    @assert sde.d == size(q₀,1)
    SDE(sde.m, ns, sde.v, sde.B, t₀, q₀, periodicity=sde.periodicity)
end

Base.ndims(sde::SDE) = sde.d


text:
## Extend convert to handle JuMP containers
# Array -> SparseAxisArray
function Base.convert(::Type{JuMPC.SparseAxisArray}, arr::Array)
    data = Dict(Tuple(k) => arr[k] for k in CartesianIndices(arr))
    return JuMPC.SparseAxisArray(data)
end

# DenseAxisArray -> SparseAxisArray
function Base.convert(::Type{JuMPC.SparseAxisArray},
                      arr::JuMPC.DenseAxisArray)
    data = Dict(k.I => arr[k] for k in keys(arr))
    return JuMPC.SparseAxisArray(data)
end

# Convert numbers to jump objects
function Base.convert(::Type{JuMP.AbstractJuMPScalar}, c::Number)
    return zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef}) + c
end

# Make workaround for keys of containers
_keys(a::JuMPC.SparseAxisArray) = keys(a.data)
_keys(a::AbstractArray) = keys(a)

## Define efficient function to check if all elements in array are equal
# method found at https://stackoverflow.com/questions/47564825/check-if-all-the-elements-of-a-julia-array-are-equal/47578613
@inline function _allequal(x::AbstractArray)::Bool
    length(x) < 2 && return true
    e1 = first(x)
    for value in x
        isequal(value, e1) || return false
    end
    return true
end

# Extend comparison for JuMP.VariableInfo
function Base.:(==)(info1::JuMP.VariableInfo, info2::JuMP.VariableInfo)::Bool
    return info1.has_lb == info2.has_lb &&
           (!info1.has_lb || info1.lower_bound == info2.lower_bound) &&
           info1.has_ub == info2.has_ub &&
           (!info1.has_ub || info1.upper_bound == info2.upper_bound) &&
           info1.has_fix == info2.has_fix &&
           (!info1.has_fix || info1.fixed_value == info2.fixed_value) &&
           info1.has_start == info2.has_start &&
           (!info1.has_start || info1.start == info2.start) &&
           info1.binary == info2.binary && info1.integer == info2.integer
end

## Convert JuMP variable info to only use Float64
# Change needed
function _make_float_info(info::JuMP.VariableInfo
    )::JuMP.VariableInfo{Float64, Float64, Float64, Float64}
    return JuMP.VariableInfo{Float64, Float64, Float64, Float64}(
        info.has_lb, info.lower_bound, info.has_ub, info.upper_bound,
        info.has_fix, info.fixed_value, info.has_start, info.start,
        info.binary, info.integer)
end

# No change needed
function _make_float_info(
    info::JuMP.VariableInfo{Float64, Float64, Float64, Float64}
    )::JuMP.VariableInfo{Float64, Float64, Float64, Float64}
    return info
end

# Make a string of NamedTuple arguments
function _kwargs2string(nt::NamedTuple)::String
    if length(nt) == 1
        return string(nt)[2:end-2]
    else
        return string(nt)[2:end-1]
    end
end


text:
module TestSetExtensions

using Distributed, Test, DeepDiffs
using Test: AbstractTestSet, DefaultTestSet, FallbackTestSet
using Test: Result, Fail, Error, Pass
export ExtendedTestSet, @includetests

"""
Includes the given test files, given as a list without their ".jl" extensions.
If none are given it will scan the directory of the calling file and include all
the julia files.
"""
macro includetests(testarg...)
    if length(testarg) == 0
        tests = []
    elseif length(testarg) == 1
        tests = testarg[1]
    else
        error("@includetests takes zero or one argument")
    end

    rootfile = "$(__source__.file)"
    mod = __module__

    quote
        tests = $tests
        rootfile = $rootfile

        if length(tests) == 0
            tests = readdir(dirname(rootfile))
            tests = filter(f->endswith(f, ".jl") && f!= basename(rootfile), tests)
        else
            tests = map(f->string(f, ".jl"), tests)
        end

        println();

        for test in tests
            print(splitext(test)[1], ": ")
            Base.include($mod, test)
            println()
        end
    end
end

struct ExtendedTestSet{T<:AbstractTestSet} <: AbstractTestSet
    wrapped::T

    ExtendedTestSet{T}(desc) where {T} = new(T(desc))
    ExtendedTestSet{FallbackTestSet}(desc) = new(FallbackTestSet())
end

struct FailDiff <: Result
    result::Fail
end

struct ExtendedTestSetException <: Exception
    msg::AbstractString
end

function ExtendedTestSet(desc; wrap=DefaultTestSet)
    ExtendedTestSet{wrap}(desc)
end

function Test.record(ts::ExtendedTestSet{T}, res::Fail) where {T}
    println("\n=====================================================")
    Test.record(ts.wrapped, res)
    print("=====================================================\n")
end

function Test.record(ts::ExtendedTestSet{DefaultTestSet}, res::Fail)
    if Distributed.myid() == 1
        println("\n=====================================================")
        printstyled(ts.wrapped.description, ": "; color = :white)

        if res.test_type === :test
            try
                test_expr = if isa(res.data, Expr)
                    res.data
                elseif isa(res.data, String)
                    Meta.parse(res.data)
                end

                if test_expr.head === :call && test_expr.args[1] === Symbol("==")
                    dd = if isa(test_expr.args[2], String) && isa(test_expr.args[3], String)
                        deepdiff(test_expr.args[2], test_expr.args[3])
                    elseif test_expr.args[2].head === :vect && test_expr.args[3].head === :vect
                        deepdiff(test_expr.args[2].args, test_expr.args[3].args)
                    elseif test_expr.args[2].head === :call && test_expr.args[3].head === :call &&
                            test_expr.args[2].args[1].head === :curly && test_expr.args[3].args[1].head === :curly
                        deepdiff(Base.eval(test_expr.args[2].args), Base.eval(test_expr.args[3].args))
                    end

                    if ! isa(dd, DeepDiffs.SimpleDiff)
                        # The test was an comparison between things we can diff,
                        # so display the diff
                        printstyled("Test Failed\n"; bold = true, color = Base.error_color())
                        println("  Expression: ", res.orig_expr)
                        printstyled("\nDiff:\n"; color = Base.info_color())
                        display(dd)
                        println()
                    else
                        # fallback to the default printing if we don't have a pretty diff
                        print(res)
                    end
                end
            catch ex
                print(res)
            end
        else
            # fallback to the default printing for non-comparisons
            print(res)
        end

        Base.show_backtrace(stdout, Test.scrub_backtrace(backtrace()))
        # show_backtrace doesn't print a trailing newline
        println("\n=====================================================")
    end
    push!(ts.wrapped.results, res)
    res, backtrace()
end

function Test.record(ts::ExtendedTestSet{T}, res::Error) where {T}
    # Ignore errors generated from failed FallbackTestSet
    if occursin(r"^Test.FallbackTestSetException", res.value) ||
           (occursin(r"^TestSetExtensions.ExtendedTestSetException", res.value) &&
            occursin("FallbackTestSetException occurred", res.value))
        throw(ExtendedTestSetException("FallbackTestSetException occurred"))
    end

    println("\n=====================================================")
    Test.record(ts.wrapped, res)
    print("=====================================================\n")
end

function Test.record(ts::ExtendedTestSet{T}, res::Pass) where {T}
    printstyled("."; color = :green)
    Test.record(ts.wrapped, res)
    res
end

Test.record(ts::ExtendedTestSet{T}, res) where {T} = Test.record(ts.wrapped, res)

# When recording DefaultTestSet results to an ExtendedTestSet{FallbackTestSet},
# throw an exception if there are any failures or errors in the DefaultTestSet.
#
# Note: this method is only needed for backward compatibility with Julia<=1.3
function Test.record(ts::ExtendedTestSet{FallbackTestSet}, res::DefaultTestSet)
    # Check for failures and errors
    passes, fails, errors, broken, c_passes, c_fails, c_errors, c_broken =
        Test.get_test_counts(res)
    if (fails > 0) || (errors > 0)
        throw(ExtendedTestSetException("Failure or error occurred in DefaultTestSet nested within FallbackTestSet."))
    end

    return res
end

function Test.finish(ts::ExtendedTestSet{T}) where {T}
    Test.get_testset_depth() == 0 && print("\n\n")
    Test.finish(ts.wrapped)
end

end # module


text:
using Latexify
using Test

arr = [1 2; 3 4]

@test latexify(arr) == replace(
raw"\begin{equation}
\left[
\begin{array}{cc}
1 & 2 \\
3 & 4 \\
\end{array}
\right]
\end{equation}
", "\r\n"=>"\n")

# Latexify.@generate_test latexify(arr; env=:inline)
@test latexify(arr; env = :inline) == replace(
raw"$\left[
\begin{array}{cc}
1 & 2 \\
3 & 4 \\
\end{array}
\right]$", "\r\n"=>"\n")

# Latexify.@generate_test latexify(arr; env=:equation)
@test latexify(arr; env = :equation) == replace(
raw"\begin{equation}
\left[
\begin{array}{cc}
1 & 2 \\
3 & 4 \\
\end{array}
\right]
\end{equation}
", "\r\n"=>"\n")

# Latexify.@generate_test latexify(arr; env=:bracket)
@test latexify(arr; env = :bracket) == replace(
raw"\[
\left[
\begin{array}{cc}
1 & 2 \\
3 & 4 \\
\end{array}
\right]\]
", "\r\n"=>"\n")

# Latexify.@generate_test latexify(arr; env=:raw)
@test latexify(arr; env = :raw) == replace(
raw"\left[
\begin{array}{cc}
1 & 2 \\
3 & 4 \\
\end{array}
\right]", "\r\n"=>"\n")

arr = [1,2,:(x/y),4]

@test latexify(arr) == replace(
raw"\begin{equation}
\left[
\begin{array}{c}
1 \\
2 \\
\frac{x}{y} \\
4 \\
\end{array}
\right]
\end{equation}
", "\r\n"=>"\n")

@test latexify(arr; transpose=true) == replace(
raw"\begin{equation}
\left[
\begin{array}{cccc}
1 & 2 & \frac{x}{y} & 4 \\
\end{array}
\right]
\end{equation}
", "\r\n"=>"\n")

@test latexify((1.0, 2), (3, 4)) == replace(
raw"\begin{equation}
\left[
\begin{array}{cc}
1.0 & 3 \\
2 & 4 \\
\end{array}
\right]
\end{equation}
", "\r\n"=>"\n")

@test latexify(((1.0, 2), (3, 4))) == replace(
raw"\begin{equation}
\left[
\begin{array}{cc}
1.0 & 3 \\
2 & 4 \\
\end{array}
\right]
\end{equation}
", "\r\n"=>"\n")

@test latexify(:(x = [1 2] * [1, 2] * [1 2; 3 4])) == replace(
raw"$x = \left[
\begin{array}{cc}
1 & 2 \\
\end{array}
\right] \cdot \left[
\begin{array}{c}
1 \\
2 \\
\end{array}
\right] \cdot \left[
\begin{array}{cc}
1 & 2 \\
3 & 4 \\
\end{array}
\right]$", "\r\n"=>"\n")

@test latexify(:(x = $arr)) == replace(
raw"$x = \left[
\begin{array}{c}
1 \\
2 \\
\frac{x}{y} \\
4 \\
\end{array}
\right]$", "\r\n"=>"\n")


text:
s = open("input") do file
    readlines(file);
end

nums = map(x->parse(Int, x), s);

sort!(nums);

function GetDist(nums)
    counts = [0, 0, 0];
    last = 0;
    for i = 1:length(nums)
        dif = nums[i] - last;
        counts[dif] += 1;
        last = nums[i];
    end

    counts[3] += 1;

    return counts;
end

dist = GetDist(nums);

dist[1] * dist[3]