text:
##===----------------------------------------------------------------------===##
##                                     CLEAN.JL                               ##
##                                                                            ##
##===----------------------------------------------------------------------===##
##                                                                            ##
##   This file cleans the Vectorize directory to provide a fresh slate        ##
##   before building or as a debugging step.                                  ##
##                                                                            ##
##===----------------------------------------------------------------------===##

# erase old function definitions to prevent issues
currdir = @__FILE__
pkgdir = currdir[1:end-13]

# Clean up directory status
@static if Sys.iswindows()
    run(`cmd /C del $(pkgdir)src\\Functions.jl`)
    run(`cmd /C copy NUL Functions.jl`)
    run(`cmd /C move Functions.jl $(pkgdir)src`)
else
    run(`rm -rf $(pkgdir)src/Functions.jl`)
    run(`touch $(pkgdir)src/Functions.jl`)
end


text:
""""
    ∇st(st, at,maxDepth=-1)
Take a gradient w.r.t. every coordinate, and store in a scattered array. The shape is (space variable, pathIndices..., outputLocation). For example, if the input signal has size 128
"""
function ∇st(st, at, maxDepth=-1)
    y, back = pullback(st, at)
    ∇ = scattered(map(x -> zeros(size(at)[1:ndims(y)]...,
            x[ndims(y):end-1]...), size(y)))
    if maxDepth < 0
        maxDepth = length(∇)
    end
    for m = 0:maxDepth
        axes(y[2])[1:end-1]
        for is in Iterators.product(axes(y[m])[1:end-1]...)
            grad = scattered(map(x -> zeros(size(x)...), y.output))
            p = pathLocs(m, is)
            grad[p] = 1
            ∇[m][:, is...] = back(grad)[1][:, 1, 1]
        end
    end
    return ∇
end

"""
    plotFirstLayer1D(j,origLoc,cline=:darkrainbow)
"""
function plotFirstLayer1D(j, origLoc, original, cline=:darkrainbow)

    space = plot(origLoc[1][:, :, j], line_z=(1:size(origLoc[1], 2))',
        legend=false, colorbar=true, color=cline,
        title="first layer gradient wavelet $j varying location")
    org = plot([original original], line_z=([-20; 1:size(origLoc[1], 2)...])', legend=false, colorbar=true, color=cline) # todo: including the colorbar here is a real hack to get them to line up
    ∇h = heatmap(origLoc[1][:, 1:end, j]', xlabel="space",
        ylabel="wavelet location", title="First layer gradient j=$j")
    ∇̂h = heatmap(log.(abs.(rfft(origLoc[1][:, 1:end, j], 1)) .^ 2)', xlabel="frequency",
        ylabel="wavelet location", title="Log Power Frequency domain j=$j")
    l = Plots.@layout [a; b{0.1h}; [b c]]
    plot(space, org, ∇h, ∇̂h, layout=l)
end

"""
    gifFirstLayer(origLoc, saveTo="tmp.gif", fps = 2)
"""
function gifFirstLayer(origLoc, firstSig, saveTo="gradientFigures/tmp.gif", fps=2)
    anim = Animation()
    for j = 1:size(origLoc[1])[end]
        plotFirstLayer1D(j, origLoc, firstSig)
        frame(anim)
    end
    gif(anim, saveTo, fps=2)
end

meanWave(wave) = sum(real.(range(0, stop=1, length=size(wave, 1)) .* wave), dims=1) ./ sum(real.(wave), dims=1)

function plotSecondLayer1D(loc, origLoc, wave1, wave2, original=false, subsamSz=(128, 85,), c=:thermal, lastDiagFreq=true)
    waveUsed = real.(ifftshift(irfft(wave1[:, loc[2]], subsamSz[1] * 2)))
    l1wave = plot(waveUsed, legend=false, titlefontsize=8, title="layer 1 ($(loc[2]))")
    annotate!(size(waveUsed, 1) * 5 / 6, maximum(waveUsed), Plots.text("freq = $(meanWave(wave1)[loc[2]])"[1:13], 5))
    waveUsed = real.(ifftshift(irfft(wave2[:, loc[1]], subsamSz[2] * 2)))
    l2wave = plot(waveUsed, legend=false, titlefontsize=8, title="layer 2 ($(loc[1]))")
    annotate!(size(waveUsed, 1) * 5 / 6, maximum(waveUsed),
        Plots.text("freq = $(meanWave(wave2)[loc[1]])"[1:13], 5))
    if original != false
        org = plot([original original], line_z=(1:size(origLoc[1], 2))', legend=false, colorbar=true, color=:darkrainbow) # todo: including the colorbar here is a real hack to get them to line up
    end
    ∇heat = heatmap(origLoc[2][:, :, loc...]', ylabel="wavelet location",
        xlabel="space", color=c, title="second layer gradient ")
    if lastDiagFreq
        ∇plt = heatmap(log.(abs.(rfft(origLoc[2][:, :, loc...], 1))), xlabel="wavelet location", ylabel="frequency", color=c, title="varying location, path $(loc[2:-1:1])")
    else
        ∇plt = plot(origLoc[2][:, :, loc...], line_z=(28:-1:1)', legend=false, color=:cividis, title="varying location, path $(loc[2:-1:1])")
    end

    if original != false
        l = Plots.@layout [[a; b{0.1h}; [c d]] e]
        return plot(∇heat, org, l1wave, l2wave, ∇plt, layout=l)
    else
        l = Plots.@layout [[a; [b c]] d]
        return plot(∇heat, l1wave, l2wave, ∇plt, layout=l)
    end
end

"""
    plotSecondLayer(stw; title="Second Layer results", xVals=-1, yVals=-1, logPower=true, toHeat=nothing, c=cgrad(:viridis, [0,.9]), threshold=0, linePalette=:greys, minLog=NaN, kwargs...)
TODO fix the similarity of these names.
xVals and yVals give the spacing of the grid, as it doesn't seem to be done
correctly by default. xVals gives the distance from the left and the right
as a tuple, while yVals gives the distance from the top and the bottom,
also as a tuple. Default values are `xVals = (.037, .852), yVals = (.056, .939)`, or if you have no title, use `xVals = (.0105, .882), yVals = (.056, .939)`
If you have no colorbar, set `xVals = (.0015, .997), yVals = (.002, .992)`
In the case that arbitrary space has been introduced, if you have a title, use `xVals = (.037, .852), yVals = (.056, .939)`, or if you have no title, use `xVals = (.0105, .882), yVals = (.056, .939)`
"""
function plotSecondLayer(stw::ScatteredOut, st; kwargs...)
    secondLayerRes = stw[2]
    if ndims(secondLayerRes) > 3
        return plotSecondLayer(secondLayerRes[:, :, :, 1], st; kwargs...)
    else
        return plotSecondLayer(secondLayerRes, st; kwargs...)
    end
end

function plotSecondLayer(stw, st; title="Second Layer results", xVals=-1, yVals=-1, logPower=true, toHeat=nothing, c=cgrad(:viridis, [0, 0.9]), threshold=0, freqsigdigits=3, linePalette=:greys, minLog=NaN, subClims=(Inf, -Inf), δt=1000, firstFreqSpacing=nothing, secondFreqSpacing=nothing, transp=true, labelRot=30, xlabel=nothing, ylabel=nothing, frameTypes=:box, miniFillAlpha=0.5, kwargs...)
    n, m = size(stw)[2:3]
    freqs = getMeanFreq(st, δt)
    freqs = map(x -> round.(x, sigdigits=freqsigdigits), freqs)[1:2]
    gr(size=2.5 .* (280, 180))
    if !(typeof(c) <: PlotUtils.ContinuousColorGradient)
        c = cgrad(c)
    end
    if toHeat == nothing
        toHeat = [norm(stw[:, i, j, 1]) for i = 1:n, j = 1:m]
    end
    if firstFreqSpacing == nothing
        firstFreqSpacing = 1:m
    end
    if secondFreqSpacing == nothing
        secondFreqSpacing = 1:n
    end
    # transp means transpose so that the second layer frequency is along the x-axis
    if transp
        xTicksFreq = (secondFreqSpacing, freqs[2][secondFreqSpacing])
        yTicksFreq = (firstFreqSpacing, freqs[1][firstFreqSpacing])
        freqs = reverse(freqs)
        nTmp = n
        n = m
        m = nTmp
        toHeat = toHeat'
        xInd = 2
        yInd = 1
        stw = permutedims(stw, (1, 3, 2, (4:ndims(stw))...))
    else
        xTicksFreq = (firstFreqSpacing, freqs[1][firstFreqSpacing])
        yTicksFreq = (secondFreqSpacing, freqs[2][secondFreqSpacing])
        (1:m, freqs[1])
        xInd = 1
        yInd = 2
    end
    if xVals == -1 && title == ""
        xVals = (0.0105, 0.882)
    elseif xVals == -1
        xVals = (0.002, 0.880)
    end
    Δx = xVals[2] - xVals[1]
    xrange = range(xVals[1] + Δx / m - Δx / (m + 3),
        stop=xVals[2] - 2 * Δx / (m + 3) + Δx / m, length=m)
    if yVals == -1
        yVals = (0.0, 0.995)
    end
    Δy = yVals[2] - yVals[1]
    yrange = range(yVals[1] + Δy / n - Δy / (n + 3), stop=yVals[2] - 2 * Δy / (n + 3) + Δy / n, length=n)
    if logPower
        toHeat = log10.(toHeat)
        if !isnan(minLog)
            toHeat = max.(minLog, toHeat)
        end
    end
    bottom = min(minimum(toHeat), subClims[1])
    top = max(subClims[2], maximum(toHeat))
    totalRange = top - bottom
    # substitute given x and y labels if needed
    if isnothing(xlabel)
        xlabel = "Layer $(xInd) frequency (Hz)"
    end
    if isnothing(ylabel)
        ylabel = "Layer $(yInd) frequency (Hz)"
    end

    if title == ""
        plt = heatmap(toHeat; yticks=yTicksFreq, xticks=xTicksFreq, tick_direction=:out, rotation=labelRot,
            xlabel=xlabel, ylabel=ylabel, c=c, clims=(bottom, top), kwargs...)
    else
        plt = heatmap(toHeat; yticks=yTicksFreq, xticks=xTicksFreq, tick_direction=:out, rotation=labelRot,
            title=title, xlabel=xlabel, ylabel=ylabel, c=c, clims=(bottom, top), kwargs...)
    end
    nPlot = 2
    for i in 1:n, j in 1:m
        if maximum(abs.(stw[:, i, j, :])) > threshold
            plt = plot!(stw[:, i, j, :], legend=false, subplot=nPlot,
                bg_inside=c[(toHeat[i, j]-bottom)/totalRange],
                ticks=nothing, palette=linePalette, frame=frameTypes,
                inset=(1, bbox(xrange[j], yrange[i], Δx / (m + 10),
                    Δy / (n + 10), :bottom, :left)), fill=0, fillalpha=miniFillAlpha)
            nPlot += 1
        end
    end
    plt
end

function jointPlot(thingToPlot, thingName, cSymbol, St; sharedColorScaling=:exp, targetExample=1, δt=1000, freqigdigits=3, sharedColorbar=true, extraPlot=nothing, allPositive=false, logPower=false, kwargs...)
    if sharedColorbar
        clims = (min(minimum.(thingToPlot)...), max(maximum.(thingToPlot)...))
        climszero = clims
        climsfirst = clims
        climssecond = clims
        toHeat = [norm(thingToPlot[2][:, i, j, 1], Inf) for i = 1:size(thingToPlot[2], 2), j = 1:size(thingToPlot[2], 3)]
    else
        climszero = (min(minimum.(thingToPlot[0])...), max(maximum.(thingToPlot[0])...))
        climsfirst = (min(minimum.(thingToPlot[1])...), max(maximum.(thingToPlot[1])...))
        climssecond = (min(minimum.(thingToPlot[2])...), max(maximum.(thingToPlot[2])...))
        toHeat = [norm(thingToPlot[2][:, i, j, 1], Inf) for i = 1:size(thingToPlot[2], 2), j = 1:size(thingToPlot[2], 3)]
    end
    firstLay = thingToPlot[1][:, :, targetExample]'
    zeroLay = thingToPlot[0][:, :, targetExample]'
    toHeat[toHeat.==0] .= -Inf    # we would like zeroes to not actually render
    firstLay[firstLay.==0] .= -Inf # for either layer

    # adjust the other parts to be log if logPower is true
    if logPower && allPositive
        absThing = map(x -> abs.(x), (thingToPlot[0], thingToPlot[1], thingToPlot[2]))
        clims = (min(minimum.(absThing)...), max(maximum.(absThing)...))
        firstLay = log10.(firstLay)
        zeroLay = log10.(abs.(zeroLay))
        climszero = log10.(clims)
        climsfirst = log10.(clims)
        climssecond = log10.(clims)
    elseif logPower
        error("not currently plotting log power and negative values")
    end

    if allPositive
        c = cgrad(cSymbol, scale=sharedColorScaling)
        cSecond = cgrad(cSymbol)
    else
        zeroAt = -climssecond[1] / (climssecond[2] - climssecond[1]) # set the mid color switch to zero
        c = cgrad(cSymbol, [0, zeroAt], scale=sharedColorScaling)
        cSecond = cgrad(cSymbol, [0, zeroAt])
    end

    # define the spatial locations as they correspond to the input
    spaceLocs = range(1, size(St)[1], length=length(zeroLay))

    p2 = plotSecondLayer(thingToPlot[2][:, :, :, targetExample], St; title="Second Layer", toHeat=toHeat, logPower=logPower, c=c, clims=climssecond, subClims=climssecond, cbar=false, xVals=(0.000, 0.993), yVals=(0.0, 0.994), transp=true, xlabel="", kwargs...)
    freqs = getMeanFreq(St, δt)
    freqs = map(x -> round.(x, sigdigits=freqigdigits), freqs)
    p1 = heatmap(firstLay, c=c, title="First Layer", clims=climsfirst, cbar=false, yticks=((1:size(firstLay, 1)), ""), xticks=((1:size(firstLay, 2)), ""), bottom_margin=-10Plots.px)
    p0 = heatmap(spaceLocs, 1:1, zeroLay, c=c, xlabel="time (ms)\nZeroth Layer", clims=climszero, cbar=false, yticks=nothing, top_margin=-10Plots.px, bottom_margin=10Plots.px)
    colorbarOnly = scatter([0, 0], [0, 1], zcolor=[0, 3], clims=climssecond, xlims=(1, 1.1), xshowaxis=false, yshowaxis=false, label="", c=c, grid=false, framestyle=:none)
    if extraPlot == nothing
        # extraPlot = scatter([0,0], [0,1], legend=false, grid=false, x="Layer 2 frequency", foreground_color_subplot=:white, top_margin=-10Plots.px, showaxis=false, yticks=nothing)
        extraPlot = plot(xlabel="Layer 2 frequency (Hz)", grid=false, xticks=(1:5, ""), showaxis=false, yticks=nothing, bottom_margin=-10Plots.px, top_margin=-20Plots.px)
        # extraPlot = heatmap(zeroLay, c=c, xlabel="location\nZeroth Layer", clims=climszero, cbar=false, yticks=nothing, top_margin=-10Plots.px, bottom_margin=10Plots.px)
    end
    titlePlot = plot(title=thingName, grid=false, showaxis=false, xticks=nothing, yticks=nothing, bottom_margin=-10Plots.px)
    lay = Plots.@layout [o{0.00001h}; [[a b; c{0.1h} d{0.1h}] b{0.04w}]]
    plot(titlePlot, p2, p1, extraPlot, p0, colorbarOnly, layout=lay)
end

"""
    addNextPath(addTo,addFrom)
    addNextPath(addFrom)

adds the next nonzero path in addFrom not present in addTo, ordered by layer, and
then the standard order for julia (dimension 1, then 2, 3,etc). If only handed `addFrom`
it makes a new pathLoc only containing the first non-empty path.

presently it only works for pathLocs whose indices are Boolean Arrays
"""
function addNextPath(addTo::pathLocs{m}, addFrom) where {m}
    shouldWeDo = foldl(checkShouldAdd, zip(addTo.indices, addFrom.indices), init=Tuple{}())
    inds = map(makeNext, addTo.indices, addFrom.indices, shouldWeDo)
    return pathLocs{m}(inds)
end
function checkShouldAdd(didPrevs, current)
    if !any(didPrevs) && current[1] != nothing
        addLoc = findfirst(current[1] .!= current[2])
        if addLoc != nothing
            return (didPrevs..., true)
        else
            return (didPrevs..., false)
        end
    else
        return (didPrevs..., false)
    end
end
function makeNext(current, fillWith, shouldDo)
    if shouldDo
        addLoc = findfirst(fillWith .!= current)
        newVersion = copy(current)
        newVersion[addLoc] = true
        return newVersion
    else
        return current
    end
end

addNextPath(addFrom::pathLocs{m}) where {m} = pathLocs{m}(foldl(makeSingle, addFrom.indices, init=Tuple{}()))

makeSingle(prev, x::Nothing) = (prev..., nothing)
makeSingle(x::Nothing) = (nothing,)
function makeSingle(x::BitArray)
    almostNull = falses(size(x))
    almostNull[findfirst(x)] = true
    return (prev..., almostNull)
end
function makeSingle(prev, x::BitArray)
    almostNull = falses(size(x))
    # only set this layer to true if all the previous are nothing
    if typeof(prev) <: Tuple{Vararg{Nothing}} # || !any.(any.(prev[prev .!=nothing]))
        almostNull[findfirst(x)] = true
    end
    return (prev..., almostNull)
end

ext: jl


text:
using AtomsBase
using InteratomicPotentials
using StaticArrays
using Unitful
using UnitfulAtomic

position = @SVector [1.0, 0.0, 0.0]
element = :Ar
atom1 = Atom(element, position * u"Å")

ϵ = 1.0
σ = 0.25
rcutoff = 2.0
lj = LennardJones(ϵ, σ, rcutoff, [:Ar, :H])

@test isa(atom1, Atom)
@test isa(lj, EmpiricalPotential)

atom2 = Atom(element, 0.25 * position * u"Å")
atom3 = Atom(element, 0.5 * position * u"Å")
atom4 = Atom(element, 0.75 * position * u"Å")
atom5 = Atom(:H, 1.1 * position * u"Å")
atom6 = Atom(element, 2.0 * position * u"Å")
# atoms = [atom1, atom2, atom3, atom4, atom5, atom6]
atoms = [atom3, atom4, atom5]
box = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
bcs = [Periodic(), Periodic(), Periodic()]
system = FlexibleSystem(atoms, box * u"Å", bcs)
n = length(system)
nnlist = neighborlist(system, rcutoff)
@test isa(nnlist, NeighborList)

true_r = [
    [SVector{3}([0.25, 0.0, 0.0]), SVector{3}([-0.40, 0.0, 0.0])],
    [SVector{3}([0.35, 0.0, 0.0])],
    []
]
nnlist_r = nnlist.r

for (ri, ri_true) in zip(nnlist_r, true_r)
    for (rij, rij_true) in zip(ri, ri_true)
        @test isapprox(sum(rij + rij_true), 0.0, atol = eps(2.0))
    end
end

@test isa(potential_energy(system, lj), AbstractFloat)
@test isa(force(system, lj), AbstractVector{<:SVector{3,<:AbstractFloat}})

@test isa(virial(system, lj), AbstractFloat)
@test isa(virial_stress(system, lj), SVector{6,<:AbstractFloat})

lj2 = LennardJones(2*ϵ, 2*σ, rcutoff, [:Ar, :H])

sum_lj = lj+lj2
mixed_potentials = [lj+lj2, lj-lj2, 2.0*lj, -2.0*lj, 2.0*lj - 1.0*lj2, lj2 / 2.0]
true_energies = [potential_energy(system, lj) + potential_energy(system, lj2),
                    potential_energy(system, lj) - potential_energy(system, lj2),
                    2.0*potential_energy(system, lj),
                    -2.0*potential_energy(system, lj),
                    2.0*potential_energy(system, lj) - 1.0potential_energy(system, lj2),
                    potential_energy(system, lj2)/2.0]

true_forces = [force(system, lj) + force(system, lj2),
                    force(system, lj) - force(system, lj2),
                    2.0*force(system, lj),
                    -2.0*force(system, lj),
                    2.0*force(system, lj) - 1.0force(system, lj2),
                    force(system, lj2)/2.0]

true_virial = [virial(system, lj) + virial(system, lj2),
                    virial(system, lj) - virial(system, lj2),
                    2.0*virial(system, lj),
                    -2.0*virial(system, lj),
                    2.0*virial(system, lj) - 1.0virial(system, lj2),
                    virial(system, lj2)/2.0]

for (result_lj, te, tf, tv) in zip(mixed_potentials, true_energies, true_forces, true_virial)
    @test isa(result_lj, MixedPotential)
    @test isa(potential_energy(system, result_lj), AbstractFloat)
    @test isapprox(potential_energy(system, result_lj), te, rtol = 1e-6)
    @test isa(force(system, result_lj), AbstractVector{<:SVector{3,<:AbstractFloat}})
    @test isapprox(force(system, result_lj), tf, rtol = 1e-6)
    @test isa(virial(system, result_lj), AbstractFloat)
    @test isapprox(virial(system, result_lj), tv, rtol = 1e-6)

    e_f = energy_and_force(system, result_lj)
    @test isapprox(e_f.e, te, rtol = 1e-6)
    @test isapprox(e_f.f, tf, rtol = 1e-6)
end


text:
#=

    Reading input

=#
function read_input()

    pass=[]
    policy=[]
    letter=[]
    open("./input", "r") do f
        while ! eof(f)
            s::String = readline(f)
            entry = split(s," ")
            push!(pass,entry[3])
            push!(letter,entry[2][1])

            policy_split = split(entry[1],"-")
            policy_tuple = (parse(Int32,policy_split[1]),
                            parse(Int32,policy_split[2]))
            push!(policy,policy_tuple)
        end
    end

    return pass, policy, letter
end

input = read_input()
pass = input[1]
policy = input[2]
letter = input[3]

#=

--- Day 2: Password Philosophy ---

Your flight departs in a few days from the coastal airport; the easiest way down to
the coast from here is via toboggan.

The shopkeeper at the North Pole Toboggan Rental Shop is having a bad day. "Something's
wrong with our computers; we can't log in!" You ask if you can take a look.

Their password database seems to be a little corrupted: some of the passwords wouldn't
have been allowed by the Official Toboggan Corporate Policy that was in effect when they
were chosen.

To try to debug the problem, they have created a list (your puzzle input) of passwords
(according to the corrupted database) and the corporate policy when that password was set.

For example, suppose you have the following list:

1-3 a: abcde
1-3 b: cdefg
2-9 c: ccccccccc

Each line gives the password policy and then the password. The password policy indicates
the lowest and highest number of times a given letter must appear for the password to be
valid. For example, 1-3 a means that the password must contain a at least 1 time and at
most 3 times.

In the above example, 2 passwords are valid. The middle password, cdefg, is not; it contains
no instances of b, but needs at least 1. The first and third passwords are valid: they contain
one a or nine c, both within the limits of their respective policies.

How many passwords are valid according to their policies?

=#

function is_valid(policy, letter, password)
    letter_count = count(==(letter), password)
    if letter_count<policy[1]
        return false
    elseif letter_count>policy[2]
        return false
    else
        return true
    end
end

function count_valid_passwords(policy, letter, pass)
    count_valid = 0
    for i in 1:size(pass,1)
        if is_valid(policy[i], letter[i], pass[i])
            count_valid+=1
        end
    end
    return count_valid
end

println(size(policy,1),", ",size(letter,1), ", ", size(pass,1))
println("Valid passwords found: ", count_valid_passwords(policy, letter, pass))

#=
--- Part Two ---

While it appears you validated the passwords correctly, they don't seem to be
what the Official Toboggan Corporate Authentication System is expecting.

The shopkeeper suddenly realizes that he just accidentally explained the password
policy rules from his old job at the sled rental place down the street! The Official
Toboggan Corporate Policy actually works a little differently.

Each policy actually describes two positions in the password, where 1 means the
first character, 2 means the second character, and so on. (Be careful; Toboggan
Corporate Policies have no concept of "index zero"!) Exactly one of these positions
must contain the given letter. Other occurrences of the letter are irrelevant for
the purposes of policy enforcement.

Given the same example list from above:

    1-3 a: abcde is valid: position 1 contains a and position 3 does not.
    1-3 b: cdefg is invalid: neither position 1 nor position 3 contains b.
    2-9 c: ccccccccc is invalid: both position 2 and position 9 contain c.

How many passwords are valid according to the new interpretation of the policies?
=#

function is_valid(policy, letter, password)
    in_fst_pos = (password[policy[1]] == letter)
    in_snd_pos = (password[policy[2]] == letter)

    if in_fst_pos
        if in_snd_pos
            return false
        else
            return true
        end
    elseif in_snd_pos
        return true
    else
        return false
    end
end

function count_valid_passwords(policy, letter, pass)
    count_valid = 0
    for i in 1:size(pass,1)
        if is_valid(policy[i], letter[i], pass[i])
            count_valid+=1
        end
    end
    return count_valid
end

println("Valid passwords found with second policy: ", count_valid_passwords(policy, letter, pass))


text:
function adaptative_time_integration_2_points(r::Rheology,p::Params,S_i,σ₃,Dⁱ,Dᵒ,ϵ̇11,ϵ̇ⁱξη,Δt,θ,tspan ; damage_growth_out=true)
    #unpack
    time_maxiter = p.solver.time_maxiter
    flags = p.flags

    # fill first values
    σᵒᵢⱼ_i, σⁱᵢⱼ_i, ϵᵒᵢⱼ_i, ϵⁱᵢⱼ_i = initialize_state_var_D(r,p,S_i,σ₃,Dⁱ,Dᵒ,θ ; coords=:band)
    Ttensors = eltype(σᵒᵢⱼ_i)

    # initialize containers
    t_vec = Float64[tspan[1]]
    S_vec = Float64[S_i]
    σⁱᵢⱼ_vec = SymmetricTensor{2,3,Ttensors}[σⁱᵢⱼ_i]
    σᵒᵢⱼ_vec = SymmetricTensor{2,3,Ttensors}[σᵒᵢⱼ_i]
    ϵⁱᵢⱼ_vec = SymmetricTensor{2,3,Ttensors}[ϵⁱᵢⱼ_i]
    Dⁱ_vec = Float64[Dⁱ]
    Dᵒ_vec = Float64[Dᵒ]

    # initialize values
    tsim::Float64 = tspan[1]
    Snext::Float64 = S_i
    σⁱᵢⱼnext::typeof(σⁱᵢⱼ_i) = σⁱᵢⱼ_i
    σᵒᵢⱼnext::typeof(σᵒᵢⱼ_i) = σᵒᵢⱼ_i
    ϵⁱᵢⱼnext::typeof(ϵⁱᵢⱼ_i) = ϵⁱᵢⱼ_i
    Dⁱnext::Float64 = Dⁱ
    Dᵒnext::Float64 = Dᵒ
    Δt_next::Float64 = Δt
    Δt_used::Float64 = Δt

    flags.bifurcation = (Dⁱ==Dᵒ) ? false : true
    i = 1 # iter counter
    while tsim < tspan[2]
      set_print_flag!(p,i,tsim)
      flags.print && print_time_iteration(i,tsim)

      # solving procedure depends on the sign of S derivative
      if flags.bifurcation
        Snext, σⁱᵢⱼnext, σᵒᵢⱼnext, ϵⁱᵢⱼnext, Dⁱnext, Dᵒnext, Δt_used, Δt_next = adaptative_solve_2_points(r,p,Snext,σ₃,σⁱᵢⱼnext,σᵒᵢⱼnext,ϵⁱᵢⱼnext,Dᵒnext,Dⁱnext,ϵ̇ⁱξη,θ,-1,Δt_next ; damage_growth_out)
      else
        σⁱᵢⱼnext, ϵⁱᵢⱼnext, Dⁱnext, Δt_used, Δt_next = adaptative_solve_1_point(r,p,σⁱᵢⱼnext,ϵⁱᵢⱼnext,Dⁱnext,ϵ̇11,θ,Δt)
        σᵒᵢⱼnext = σⁱᵢⱼnext
        Dᵒnext = Dⁱnext
        σᵒᵢⱼnext_principal = principal_coords(σᵒᵢⱼnext,θ)
        Snext = σᵒᵢⱼnext_principal[1,1]/σᵒᵢⱼnext_principal[2,2]
      end

      if !flags.bifurcation # activate bifurcation procedure if S starts to decrease or if the derivative is zero.
        if Snext-S_vec[end] <= 0
          flags.bifurcation = true
          Snext, σⁱᵢⱼnext, σᵒᵢⱼnext, ϵⁱᵢⱼnext, Dⁱnext, Dᵒnext, Δt_used, Δt_next = adaptative_solve_2_points(r,p,Snext,σ₃,σⁱᵢⱼnext,σᵒᵢⱼnext,ϵⁱᵢⱼnext,Dᵒnext,Dⁱnext,ϵ̇ⁱξη,θ,-1,Δt_next ; damage_growth_out)
        end
      end

      # test :
      p.flags.nan && println("Nan flag true in main function, should return")
      # update tsim if no nans
      !flags.nan ? (tsim += Δt_used) : nothing

      set_save_flag!(p,i,tsim)
      if flags.save
        push!(S_vec,Snext)
        push!(σⁱᵢⱼ_vec,σⁱᵢⱼnext)
        push!(σᵒᵢⱼ_vec,σᵒᵢⱼnext)
        push!(ϵⁱᵢⱼ_vec,ϵⁱᵢⱼnext)
        push!(Dⁱ_vec,Dⁱnext)
        push!(Dᵒ_vec,Dᵒnext)
        push!(t_vec,tsim)
      end

      i += 1

      # break loop under conditions
      flags.nan && break
      if !isnothing(time_maxiter)
        (length(t_vec)==time_maxiter+1) && break
      end
      D_max = 0.999
      (Dⁱnext > D_max) && (@info("Dⁱ reached $D_max, ending simulation") ; break)
      (Dᵒnext > D_max) && (@info("Dᵒ reached $D_max, ending simulation") ; break)
    end
    return t_vec, S_vec, σⁱᵢⱼ_vec, σᵒᵢⱼ_vec, ϵⁱᵢⱼ_vec, Dⁱ_vec, Dᵒ_vec
  end

  function initialize_state_var_D(r,p,S_i,σ₃,Dⁱ,Dᵒ,θ ; coords=:band)
    ps = p.solver
    # get outside state to be used as reference
    σᵒᵢⱼ_i_principal = build_principal_stress_tensor(r,S_i,σ₃,Dᵒ ; abstol=1e-15) # takes care of the plane strain constraint by solving non linear out of plane strain wrt σ₃₃ using Newton algorithm
    ϵᵒᵢⱼ_i_principal = compute_ϵij(r,Dᵒ,σᵒᵢⱼ_i_principal)
    σᵒᵢⱼ_i = band_coords(σᵒᵢⱼ_i_principal,θ)
    ϵᵒᵢⱼ_i = band_coords(ϵᵒᵢⱼ_i_principal,θ)

    σⁱᵢⱼ_i_principal = build_principal_stress_tensor(r,S_i,σ₃,Dⁱ ; abstol=1e-15) # takes care of the plane strain constraint by solving non linear out of plane strain wrt σ₃₃ using Newton algorithm
    ϵⁱᵢⱼ_i_principal = compute_ϵij(r,Dⁱ,σⁱᵢⱼ_i_principal)
    σⁱᵢⱼ_i = band_coords(σⁱᵢⱼ_i_principal,θ)
    ϵⁱᵢⱼ_i = band_coords(ϵⁱᵢⱼ_i_principal,θ)

    # insert stress and strain componant constrained by stress continuity and strain compatibility
    σⁱᵢⱼ_guess = insert_into(σⁱᵢⱼ_i, (σᵒᵢⱼ_i[2,2],σᵒᵢⱼ_i[1,2]), ((2,2),(1,2)))
    ϵⁱᵢⱼ_guess = insert_into(ϵⁱᵢⱼ_i, ϵᵒᵢⱼ_i[1,1], (1,1))

    # initialize first guess for u = [ σξξ, σoop, ϵηη, ϵξη ]
    u = Vec(σⁱᵢⱼ_i[1,1], σⁱᵢⱼ_i[3,3], ϵⁱᵢⱼ_i[2,2], ϵⁱᵢⱼ_i[1,2])
    result = DiffResults.JacobianResult(u)

    for i in 1:ps.newton_maxiter
      # get residual and its gradient with respect to u
      ForwardDiff.jacobian!(result, u -> residual_initialize(r,Dⁱ,ϵⁱᵢⱼ_guess,σⁱᵢⱼ_guess,u), u)
      ∇res = DiffResults.jacobian(result)
      res  = DiffResults.value(result)

      # update u with Newton algo
      δu = - ∇res\res
      u = u + δu
      #@debug "δu = $δu"
      #@debug "typeof(u) = $(typeof(u))"

      (norm(res) <= ps.newton_abstol) && break

      (i == ps.newton_maxiter) && @debug("ending norm res = $(norm(res))")#("max newton iteration reached ($i), residual still higher than abstol with $(norm(res))")
    end
    σξξ, σoop, ϵηη, ϵξη = u
    σⁱᵢⱼ = insert_into(σⁱᵢⱼ_guess, (σξξ, σoop), ((1,1),(3,3)))
    ϵⁱᵢⱼ = insert_into(ϵⁱᵢⱼ_guess, (ϵηη,ϵξη), ((2,2),(1,2)))
    if coords == :principal
      σᵒᵢⱼ = σᵒᵢⱼ_i_principal
      ϵᵒᵢⱼ = ϵᵒᵢⱼ_i_principal
      σⁱᵢⱼ = principal_coords(σⁱᵢⱼ,θ)
      ϵⁱᵢⱼ = principal_coords(ϵⁱᵢⱼ,θ)
    elseif coords == :band
      σᵒᵢⱼ = σᵒᵢⱼ_i
      ϵᵒᵢⱼ = ϵᵒᵢⱼ_i
    else
      @error "`coords` keyword argument $coords is not recognized. Use `:principal` or `:band` instead, `:band` being the default"
    end
    return σᵒᵢⱼ, σⁱᵢⱼ, ϵᵒᵢⱼ, ϵⁱᵢⱼ
  end

  function residual_initialize(r,Dⁱ,ϵⁱᵢⱼ_guess,σⁱᵢⱼ_guess,u)
    σξξ, σoop, ϵηη, ϵξη = u
    σⁱᵢⱼ = insert_into(σⁱᵢⱼ_guess, (σξξ, σoop), ((1,1),(3,3)))
    ϵⁱᵢⱼ = insert_into(ϵⁱᵢⱼ_guess, (ϵηη,ϵξη), ((2,2),(1,2)))
    ϵⁱᵢⱼ_analytical = compute_ϵij(r,Dⁱ,σⁱᵢⱼ)
    Δϵⁱᵢⱼ = ϵⁱᵢⱼ_analytical - ϵⁱᵢⱼ
    res = Vec(Δϵⁱᵢⱼ[1,1], Δϵⁱᵢⱼ[2,2], Δϵⁱᵢⱼ[3,3], Δϵⁱᵢⱼ[1,2])
    return res
  end

  function time_integration_2_points(r,p,S_i,σ₃,D_i,ϵ̇11,ϵ̇ⁱξη,Δt,θ,tspan ; damage_growth_out=true)
    flags = p.flags
    #initialize vectors
    t_vec = tspan[1]:Δt:tspan[2]
    σⁱᵢⱼ_vec = Vector{SymmetricTensor{2,3}}(undef,length(t_vec))
    σᵒᵢⱼ_vec = Vector{SymmetricTensor{2,3}}(undef,length(t_vec))
    ϵⁱᵢⱼ_vec = similar(σⁱᵢⱼ_vec)
    Dⁱ_vec = Vector{Float64}(undef,length(t_vec))
    Dᵒ_vec = Vector{Float64}(undef,length(t_vec))
    S_vec = Vector{Float64}(undef,length(t_vec))

    # fill first values
    σᵢⱼ_i = build_principal_stress_tensor(r,S_i,σ₃,D_i ; abstol=1e-15) # takes care of the plane strain constraint by solving non linear out of plane strain wrt σ₃₃ using Newton algorithm
    ϵᵢⱼ_i = compute_ϵij(r,D_i,σᵢⱼ_i)

    σⁱᵢⱼ_vec[begin] = band_coords(σᵢⱼ_i,θ)
    σᵒᵢⱼ_vec[begin] = band_coords(σᵢⱼ_i,θ)
    ϵⁱᵢⱼ_vec[begin] = band_coords(ϵᵢⱼ_i,θ) # test that stress-strain relation is still valid after rotation !!
    Dⁱ_vec[begin] = D_i
    Dᵒ_vec[begin] = D_i
    S_vec[begin] = S_i

    last_iter = 0
    for i in 2:length(t_vec)
      set_print_flag!(p,i,t_vec[i-1])
      flags.print && print_time_iteration(i,t_vec[i-1])

      # solving procedure depends on the sign of S derivative
      if flags.bifurcation
        S_vec[i], σⁱᵢⱼ_vec[i], σᵒᵢⱼ_vec[i], ϵⁱᵢⱼ_vec[i], Dⁱ_vec[i], Dᵒ_vec[i], _ = solve_2_points(r,p,S_vec[i-1],σ₃,σⁱᵢⱼ_vec[i-1],σᵒᵢⱼ_vec[i-1],ϵⁱᵢⱼ_vec[i-1],Dᵒ_vec[i-1],Dⁱ_vec[i-1],ϵ̇ⁱξη,θ,Δt ; damage_growth_out)
      else
        σⁱᵢⱼ_vec[i], ϵⁱᵢⱼ_vec[i], Dⁱ_vec[i], _ = solve_1_point(r,p,σⁱᵢⱼ_vec[i-1],ϵⁱᵢⱼ_vec[i-1],Dⁱ_vec[i-1],ϵ̇11,θ,Δt)
        σᵒᵢⱼ_vec[i] = σⁱᵢⱼ_vec[i]
        Dᵒ_vec[i] = Dⁱ_vec[i]
        σᵒᵢⱼ_principal = principal_coords(σᵒᵢⱼ_vec[i],θ)
        S_vec[i] = σᵒᵢⱼ_principal[1,1]/σᵒᵢⱼ_principal[2,2]
      end

      if !flags.bifurcation # activate bifurcation procedure if S starts to decrease or if the derivative is zero.
        (S_vec[i]-S_vec[i-1] <= 0) && (flags.bifurcation = true)
        if flags.bifurcation
          S_vec[i], σⁱᵢⱼ_vec[i], σᵒᵢⱼ_vec[i], ϵⁱᵢⱼ_vec[i], Dⁱ_vec[i], Dᵒ_vec[i], _ = solve_2_points(r,p,S_vec[i-1],σ₃,σⁱᵢⱼ_vec[i-1],σᵒᵢⱼ_vec[i-1],ϵⁱᵢⱼ_vec[i-1],Dᵒ_vec[i-1],Dⁱ_vec[i-1],ϵ̇ⁱξη,θ,Δt ; damage_growth_out)
        end
      end
      if Dⁱ_vec[i] > 0.999
        last_iter = i
        break
      end
    end
    li = last_iter
    return t_vec[1:li], S_vec[1:li], σⁱᵢⱼ_vec[1:li], σᵒᵢⱼ_vec[1:li], ϵⁱᵢⱼ_vec[1:li], Dⁱ_vec[1:li], Dᵒ_vec[1:li]
  end


text:
using RecursiveArrayTools;
using DataStructures;
using JSON;
using TimerOutputs;
const to = TimerOutput()

include("jsonhandler.jl");
include("astar.jl");
include("bfs.jl");
include("ids.jl");
include("tree.jl");
include("utils.jl");

max_stars_repo_name: brice-allen/npuz

max_issues_repo_name: brice-allen/npuz

max_forks_repo_name: brice-allen/npuz


text:
using DashBootstrapComponents, DashHtmlComponents

simple_jumbotron = dbc_jumbotron([
    html_h1("Jumbotron", className="display-3"),
    html_p(
        "Use a jumbotron to call attention to " * "featured content or information.",
        className="lead",
    ),
    html_hr(className="my-2"),
    html_p(
        "Jumbotrons use utility classes for typography and " *
        "spacing to suit the larger container.",
    ),
    html_p(dbc_button("Learn more", color="primary"), className="lead"),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
module Enzyme

export autodiff, autodiff_deferred
export Const, Active, Duplicated, DuplicatedNoNeed

"""
    abstract type Annotation{T}

Abstract type for [`autodiff`](@ref) function argument wrappers like
[`Const`](@ref), [`Active`](@ref) and [`Duplicated`](@ref).
"""
abstract type Annotation{T} end

"""
    struct Const{T} <: Annotation{T}

Constructor: `Const(x)`

Mark a function argument `x` of [`autodiff`](@ref) as constant,
Enzyme will not auto-differentiate in respect `Const` arguments.
"""
struct Const{T} <: Annotation{T}
    val::T
end

# To deal with Const(Int) and prevent it to go to `Const{DataType}(T)`
Const(::Type{T}) where T = Const{Type{T}}(T)

"""
    struct Active{T} <: Annotation{T}

Constructor: `Active(x)`

Mark a function argument `x` of [`autodiff`](@ref) as active,
Enzyme will auto-differentiate in respect `Active` arguments.

!!! note

    Enzyme gradients with respect to integer values are zero.
    [`Active`](@ref) will automatically convert plain integers to floating
    point values, but cannot do so for integer values in tuples and structs.
"""
struct Active{T} <: Annotation{T}
    val::T
end

Active(i::Integer) = Active(float(i))

"""
    struct Duplicated{T} <: Annotation{T}

Constructor: `Duplicated(x, ∂f_∂x)`

Mark a function argument `x` of [`autodiff`](@ref) as duplicated, Enzyme will
auto-differentiate in respect to such arguments, with `dx` acting as an
accumulator for gradients (so ``\\partial f / \\partial x`` will be *added to*)
`∂f_∂x`.
"""
struct Duplicated{T} <: Annotation{T}
    val::T
    dval::T
end

struct DuplicatedNoNeed{T} <: Annotation{T}
    val::T
    dval::T
end

Base.eltype(::Type{<:Annotation{T}}) where T = T

function guess_activity(T)
    if T <: AbstractFloat || T <: Complex{<:AbstractFloat}
        return Active{T}
    elseif T <: AbstractArray
        return Duplicated{T}
    else
        return Const{T}
    end
end

import LLVM

include("api.jl")
include("logic.jl")
include("typeanalysis.jl")
include("typetree.jl")
include("utils.jl")
include("compiler.jl")

# @inline annotate() = ()
# @inline annotate(arg::A, args::Vararg{Any, N}) where {A<:Annotation, N} = (arg, annotate(args...)...)
# @inline annotate(arg, args::Vararg{Any, N}) where N = (Const(arg), annotate(args...)...)

@inline function annotate(args::Vararg{Any, N}) where N
    ntuple(Val(N)) do i
        Base.@_inline_meta
        arg = @inbounds args[i]
        if arg isa Annotation
            return arg
        else
            return Const(arg)
        end
    end
end

prepare_cc() = ()
prepare_cc(arg::Duplicated, args...) = (arg.val, arg.dval, prepare_cc(args...)...)
prepare_cc(arg::DuplicatedNoNeed, args...) = (arg.val, arg.dval, prepare_cc(args...)...)
prepare_cc(arg::Annotation, args...) = (arg.val, prepare_cc(args...)...)

"""
    autodiff(f, Activity, args...)

Auto-differentiate function `f` at arguments `args`.

Limitations:

* `f` may only return a `Real` (of a built-in/primitive type) or `nothing`,
  not an array, struct, `BigFloat`, etc. To handle vector-valued return
  types, use a mutating `f!` that returns `nothing` and stores it's return
  value in one of the arguments, which must be wrapped in a
  [`Duplicated`](@ref).

* `f` may not allocate memory, this restriction is likely to be removed in
  future versions. Technically it can currently allocate memory directly, but not in a function called by `f`.

`args` may be numbers, arrays, structs of numbers, structs of arrays and so
on. Enzyme will only differentiate in respect to arguments that are wrapped
in an [`Active`](@ref) (for immutable arguments like primitive types and
structs thereof) or [`Duplicated`](@ref) (for mutable arguments like arrays,
`Ref`s and structs thereof). Non-annotated arguments will automatically be
treated as [`Const`](@ref).

`Activity` is the Activity of the return value, it may be `Const` or `Active`.

Example:

```jldoctest
using Enzyme

a = 4.2
b = [2.2, 3.3]; ∂f_∂b = zero(b)
c = 55; d = 9

f(a, b, c, d) = a * √(b[1]^2 + b[2]^2) + c^2 * d^2
∂f_∂a, ∂f_∂d = autodiff(f, Active, Active(a), Duplicated(b, ∂f_∂b), c, Active(d))

# output

(3.966106403010388, 54450.0)
```

here, `autodiff` returns a tuple
``(\\partial f/\\partial a, \\partial f/\\partial d)``,
while ``\\partial f/\\partial b`` will be *added to* `∂f_∂b` (but not returned).
`c` will be treated as `Const(c)`.

!!! note

    Enzyme gradients with respect to integer values are zero.
    [`Active`](@ref) will automatically convert plain integers to floating
    point values, but cannot do so for integer values in tuples and structs.
"""
@inline function autodiff(f::F, ::Type{A}, args...) where {F, A<:Annotation}
    args′  = annotate(args...)
    tt′    = Tuple{map(Core.Typeof, args′)...}
    if A <: Active
        tt    = Tuple{map(T->eltype(Core.Typeof(T)), args′)...}
        rt = Core.Compiler.return_type(f, tt)
        if !allocatedinline(rt)
            forward, adjoint = Enzyme.Compiler.thunk(f, #=df=#nothing, Duplicated{rt}, tt′, #=Split=# Val(true))
            res = forward(args′...)
            tape = res[1]
            if res[3] isa Base.RefValue
                res[3][] += one(eltype(typeof(res[3])))
            else
                res[3] += one(eltype(typeof(res[3])))
            end
            return adjoint(args′..., tape)
        end
    elseif A <: Duplicated
        throw(ErrorException("Duplicated Returns not yet handled"))
    end
    thunk = Enzyme.Compiler.thunk(f, #=df=#nothing, A, tt′, #=Split=# Val(false))
    rt = eltype(Compiler.return_type(thunk))
    if A <: Active
        args′ = (args′..., one(rt))
    end
    thunk(args′...)
end

@inline function autodiff(dupf::Duplicated{F}, ::Type{A}, args...) where {F, A<:Annotation}
    args′  = annotate(args...)
    tt′    = Tuple{map(Core.Typeof, args′)...}
    thunk = Enzyme.Compiler.thunk(#=f=#dupf.val, #=df=#dupf.dval, A, tt′, #=Split=# Val(false))
    if A <: Active
        rt = eltype(Compiler.return_type(thunk))
        args′ = (args′..., one(rt))
    end
    thunk(args′...)
end

"""
    autodiff(f, args...)

Like [`autodiff`](@ref) but will try to guess the activity of the return value.
"""
@inline function autodiff(f::F, args...) where {F}
    args′ = annotate(args...)
    tt′   = Tuple{map(Core.Typeof, args′)...}
    tt    = Tuple{map(T->eltype(Core.Typeof(T)), args′)...}
    rt    = Core.Compiler.return_type(f, tt)
    A     = guess_activity(rt)
    autodiff(f, A, args′...)
end

@inline function autodiff(dupf::Duplicated{F}, args...) where {F}
    args′ = annotate(args...)
    tt′   = Tuple{map(Core.Typeof, args′)...}
    tt    = Tuple{map(T->eltype(Core.Typeof(T)), args′)...}
    rt    = Core.Compiler.return_type(dupf.val, tt)
    A     = guess_activity(rt)
    autodiff(dupf, A, args′...)
end

"""
    autodiff_deferred(f, Activity, args...)

Same as [`autodiff`](@ref) but uses deferred compilation to support usage in GPU
code, as well as high-order differentiation.
"""
@inline function autodiff_deferred(f::F, ::Type{A}, args...) where {F, A<:Annotation}
    args′ = annotate(args...)
    tt′   = Tuple{map(Core.Typeof, args′)...}
    if A isa UnionAll
        tt = Tuple{map(T->eltype(Core.Typeof(T)), args′)...}
        rt = Core.Compiler.return_type(f, tt)
        rt = A{rt}
    else
        @assert A isa DataType
        rt = A
    end

    if eltype(rt) == Union{}
        error("Return type inferred to be Union{}. Giving up.")
    end

    ptr   = Compiler.deferred_codegen(Val(f), Val(tt′), Val(rt))
    thunk = Compiler.CombinedAdjointThunk{F, rt, tt′, Nothing}(f, ptr, #=df=#nothing)
    if rt <: Active
        args′ = (args′..., one(eltype(rt)))
    elseif A <: Duplicated
        throw(ErrorException("Duplicated Returns not yet handled"))
    end
    thunk(args′...)
end

"""
    autodiff_deferred(f, args...)

Like [`autodiff_deferred`](@ref) but will try to guess the activity of the return value.
"""
@inline function autodiff_deferred(f::F, args...) where {F}
    args′ = annotate(args...)
    tt′   = Tuple{map(Core.Typeof, args′)...}
    tt    = Tuple{map(T->eltype(Core.Typeof(T)), args′)...}
    rt    = Core.Compiler.return_type(f, tt)
    rt    = guess_activity(rt)
    autodiff_deferred(f, rt, args′...)
end

using Adapt
Adapt.adapt_structure(to, x::Duplicated) = Duplicated(adapt(to, x.val), adapt(to, x.dval))
Adapt.adapt_structure(to, x::DuplicatedNoNeed) = DuplicatedNoNeed(adapt(to, x.val), adapt(to, x.dval))
Adapt.adapt_structure(to, x::Const) = Const(adapt(to, x.val))
Adapt.adapt_structure(to, x::Active) = Active(adapt(to, x.val))

end # module


text:
# ==============================================================================
# Random linear switching system
#
# system type: hybrid system with linear dynamics
# state dimension: 5
# input dimension: 1
# modes: 5
# transitions: 5
#
# The model is taken from [1].
#
# [1] https://ths.rwth-aachen.de/research/projects/hypro/5-dimensional-switching-linear-system/
# ==============================================================================
using HybridSystems, MathematicalSystems, LazySets, SymEngine
using SX: readsxmodel, _get_coeffs

function linearSwitching_model()
    n = 5  # state dimension
    U = Interval(-1.0, 1.0)  # common input domain
    m = dim(U)  # input dimension
    ε = 1e-6  # auxiliary bloating of guards for ensuring intersection

    # hybrid automaton
    HA = LightAutomaton(5)

    # SpaceEx model
    file = @relpath "SpaceEx/model.xml"
    model = readsxmodel(file, raw_dict=true)
    variables = convert.(Basic, [f.args[1].args[1] for f in model["flows"][1]])
    inputs = [convert(Basic, :u)]
    load_dynamics(loc) = _get_coeffs(model["flows"][loc], n, m, variables, inputs)

    # mode 1
    A, B, _ = load_dynamics(1)
    X = HalfSpace([-1.0, 0.0, 0.0, 0.0, 0.0], -3.0 + ε)  # x1 ≥ 3
    q1 = @system(x' = Ax + Bu, x ∈ X, u ∈ U)

    # mode 2
    A, B, c = load_dynamics(2)
    X = HalfSpace([-1.0, 0.0, 0.0, 0.0, 0.0], -2.0 + ε)  # x1 ≥ 2
    q2 = @system(x' = Ax + Bu, x ∈ X, u ∈ U)

    # mode 3
    A, B, _ = load_dynamics(3)
    X = HalfSpace([-1.0, 0.0, 0.0, 0.0, 0.0], -1.0 + ε)  # x1 ≥ 1
    q3 = @system(x' = Ax + Bu, x ∈ X, u ∈ U)

    # mode 4
    A, B, _ = load_dynamics(4)
    X = HalfSpace([-1.0, 0.0, 0.0, 0.0, 0.0], 0.0 + ε)  # x1 ≥ 0
    q4 = @system(x' = Ax + Bu, x ∈ X, u ∈ U)

    # mode 5
    A, B, _ = load_dynamics(5)
    X = HalfSpace([1.0, 0.0, 0.0, 0.0, 0.0], 1.0 + ε)  # x1 ≤ 1
    q5 = @system(x' = Ax + Bu, x ∈ X, u ∈ U)

    # transition mode 1 → mode 2
    add_transition!(HA, 1, 2, 1)
    G = HPolyhedron([HalfSpace([1.0, 0.0, 0.0, 0.0, 0.0], 3.0 + ε),
                     HalfSpace([-1.0, 0.0, 0.0, 0.0, 0.0], -3.0 + ε)])  # x1 = 3
    R12 = ConstrainedIdentityMap(5, G)

    # transition mode 2 → mode 3
    add_transition!(HA, 2, 3, 2)
    G = HPolyhedron([HalfSpace([1.0, 0.0, 0.0, 0.0, 0.0], 2.0 + ε),
                     HalfSpace([-1.0, 0.0, 0.0, 0.0, 0.0], -2.0 + ε)])  # x1 = 2
    R23 = ConstrainedIdentityMap(5, G)

    # transition mode 3 → mode 4
    add_transition!(HA, 3, 4, 3)
    G = HPolyhedron([HalfSpace([1.0, 0.0, 0.0, 0.0, 0.0], 1.0 + ε),
                     HalfSpace([-1.0, 0.0, 0.0, 0.0, 0.0], -1.0 + ε)])  # x1 = 1
    R34 = ConstrainedIdentityMap(5, G)

    # transition mode 4 → mode 5
    add_transition!(HA, 4, 5, 4)
    G = HPolyhedron([HalfSpace([1.0, 0.0, 0.0, 0.0, 0.0], 0.0 + ε),
                     HalfSpace([-1.0, 0.0, 0.0, 0.0, 0.0], -0.0 + ε)])  # x1 = 0
    R45 = ConstrainedIdentityMap(5, G)

    # transition mode 5 → mode 1
    add_transition!(HA, 5, 1, 5)
    G = HPolyhedron([HalfSpace([1.0, 0.0, 0.0, 0.0, 0.0], 1.0 + ε),
                     HalfSpace([-1.0, 0.0, 0.0, 0.0, 0.0], -1.0 + ε)])  # x1 = 1
    R51 = ConstrainedIdentityMap(5, G)

    # hybrid system
    S = HybridSystem(HA, [q1, q2, q3, q4, q5], [R12, R23, R34, R45, R51],
                     fill(AutonomousSwitching(), 5))

    return S
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
function mapping_factor(obj::IonosphericMappingFunction, arg0::RealFieldElement)
    return jcall(obj, "mappingFactor", RealFieldElement, (RealFieldElement,), arg0)
end

function mapping_factor(obj::IonosphericMappingFunction, arg0::jdouble)
    return jcall(obj, "mappingFactor", jdouble, (jdouble,), arg0)
end


text:
export
    LSFI

"""
    function LSFI(query_vars;
        increment = 10, start_size = increment, max_iterations = 100, start_depth = 1)

A lazy inference algorithm that uses variable elimination at every step.

# Arguments
- `query_vars`: Variables that can be queried after each `refine` step
- `increment`: The increment to range size on every iteration
- `start_size`: The starting range size
- `max_iterations`: The maximum number of refinement steps
- `start_depth`: The depth of recursive expansion in the first iteration
"""
function LSFI(query_vars;
        increment = 10, start_size = increment, max_iterations = 100, start_depth = 1)
    # query_vars = [get_node(net, :out)]
    maker(size, depth) = VE(query_vars; depth = depth, bounds = true, range_size =size)
    return LazyInference(maker; increment = increment, start_size = start_size,
        max_iterations = max_iterations, start_depth = start_depth)
end


text:
# Simple "print"-like rendering, use "{ ... }" brackets for compactness
function Base.show(io::IO, inds::AbstractIndices)
    limit = get(io, :limit, false) ? Int64(10) : typemax(Int64)
    comma = false
    print(io, "{")
    for i in inds
        if comma
            print(io, ", ")
        end
        if limit == 0
            print(io, "…")
            break
        end
        show(io, i)
        comma = true
        limit -= 1
    end
    print(io, "}")
end

function Base.show(io::IO, d::AbstractDictionary)
    limit = get(io, :limit, false) ? Int64(10) : typemax(Int64)
    comma = false
    print(io, "{")
    for (i, v) in pairs(d)
        if comma
            print(io, ", ")
        end
        if limit == 0
            print(io, "…")
            break
        end
        show(io, i)
        print(io, " │ ")
        show(io, v)
        comma = true
        limit -= 1
    end
    print(io, "}")
end

# The "display"-like rendering
function Base.show(io::IO, ::MIME"text/plain", d::AbstractIndices)
    if isempty(d)
        print(io, "0-element $(typeof(d))")
        return
    end

    # Designed to be efficient for very large sets of unknown lengths

    n_lines = max(Int64(3), get(io, :limit, false) ? Int64(displaysize(io)[1] - 4) : typemax(Int64))
    n_cols = max(Int64(2), get(io, :limit, false) ? Int64(displaysize(io)[2] - 1) : typemax(Int64))
    n_lines_top = n_lines ÷ Int64(2)
    n_lines_bottom = n_lines - n_lines_top

    # First we collect strings of all the relevant elements
    top_ind_strs = Vector{String}()
    bottom_ind_strs = Vector{String}()

    top_lines = Int64(1)
    top_full = false
    top_last_index = Base.RefValue{keytype(d)}()
    for i in keys(d)
        push!(top_ind_strs, sprint(show, i, context = io, sizehint = 0))
        top_lines += 1
        if top_lines > n_lines_top
            top_full = true
            top_last_index[] = i
            break
        end
    end

    bottom_lines = Int64(1)
    bottom_full = false
    if top_full
        for i in Iterators.reverse(keys(d))
            if bottom_full
                if isequal(i, top_last_index[])
                    bottom_full = false # overide this, we don't need the ⋮
                else
                    bottom_ind_strs[end] = "⋮"
                end
                break
            end

            if isequal(i, top_last_index[])
                # Already rendered, we are finished
                break
            end

            push!(bottom_ind_strs, sprint(show, i, context = io, sizehint = 0))
            bottom_lines += 1
            if bottom_lines > n_lines_bottom
                bottom_full = true
                # We check the next element to see if this one should be a ⋮
            end
        end
        ind_strs = vcat(top_ind_strs, reverse(bottom_ind_strs))
    else
        ind_strs = top_ind_strs
    end

    if Base.IteratorSize(d) === Base.SizeUnknown()
        if bottom_full
            print(io, "Greater than $(length(ind_strs))-element $(typeof(d))")
        else
            print(io, "$(length(ind_strs))-element $(typeof(d))")
        end
    else
        print(io, "$(length(d))-element $(typeof(d))")
    end

    # Now find padding sizes
    max_ind_width = maximum(textwidth, ind_strs)
    if max_ind_width + 1 > n_cols
        shrink_to!(ind_strs, n_cols)
    end

    for ind_str in ind_strs
        print(io, "\n ")
        print(io, ind_str)
    end
end

# The "display"-like rendering
function Base.show(io::IO, ::MIME"text/plain", d::AbstractDictionary)
    if isempty(d)
        print(io, "0-element $(typeof(d))")
        return
    end

    # Designed to be efficient for very large sets of unknown lengths

    n_lines = max(Int64(3), get(io, :limit, false) ? Int64(displaysize(io)[1] - 4) : typemax(Int64))
    n_cols = max(Int64(8), get(io, :limit, false) ? Int64(displaysize(io)[2] - 4) : typemax(Int64))
    n_lines_top = n_lines ÷ Int64(2)
    n_lines_bottom = n_lines - n_lines_top

    # First we collect strings of all the relevant elements
    top_ind_strs = Vector{String}()
    top_val_strs = Vector{String}()
    bottom_val_strs = Vector{String}()
    bottom_ind_strs = Vector{String}()

    top_lines = Int64(1)
    top_full = false
    top_last_index = Base.RefValue{keytype(d)}()
    for i in keys(d)
        push!(top_ind_strs, sprint(show, i, context = io, sizehint = 0))
        if isassigned(d, i)
            push!(top_val_strs, sprint(show, d[i], context = io, sizehint = 0))
        else
            push!(top_val_strs, "#undef")
        end
        top_lines += 1
        if top_lines > n_lines_top
            top_full = true
            top_last_index[] = i
            break
        end
    end

    bottom_lines = Int64(1)
    bottom_full = false
    if top_full
        for i in Iterators.reverse(keys(d))
            if bottom_full
                if isequal(i, top_last_index[])
                    bottom_full = false # overide this, we don't need the ⋮
                else
                    bottom_ind_strs[end] = "⋮"
                    bottom_val_strs[end] = "⋮"
                end
                break
            end

            if isequal(i, top_last_index[])
                # Already rendered, we are finished
                break
            end

            push!(bottom_ind_strs, sprint(show, i, context = io, sizehint = 0))
            if isassigned(d, i)
                push!(bottom_val_strs, sprint(show, d[i], context = io, sizehint = 0))
            else
                push!(bottom_val_strs, "#undef")
            end

            bottom_lines += 1
            if bottom_lines > n_lines_bottom
                bottom_full = true
                # We check the next element to see if this one should be a ⋮
            end
        end
        ind_strs = vcat(top_ind_strs, reverse(bottom_ind_strs))
        val_strs = vcat(top_val_strs, reverse(bottom_val_strs))
    else
        ind_strs = top_ind_strs
        val_strs = top_val_strs
    end

    if Base.IteratorSize(d) === Base.SizeUnknown()
        if bottom_full
            print(io, "Greater than $(length(ind_strs))-element $(typeof(d))")
        else
            print(io, "$(length(ind_strs))-element $(typeof(d))")
        end
    else
        print(io, "$(length(d))-element $(typeof(d))")
    end

    # Now find padding sizes
    max_ind_width = maximum(textwidth, ind_strs)
    max_val_width = maximum(textwidth, val_strs)
    if max_ind_width + max_val_width > n_cols
        val_width = max_val_width
        ind_width = max_ind_width
        while ind_width + val_width > n_cols
            if ind_width > val_width
                ind_width -= 1
            else
                val_width -= 1
            end
        end
        if ind_width != max_ind_width
            shrink_to!(ind_strs, ind_width)
        end
        if val_width != max_val_width
            shrink_to!(val_strs, val_width)
        end
    else
        ind_width = max_ind_width
    end

    for (ind_str, val_str) in zip(ind_strs, val_strs)
        print(io, "\n ")
        print(io, " " ^ max(0, ind_width - textwidth(ind_str)))
        print(io, ind_str)
        print(io, " │ ")
        print(io, val_str)
    end
end

function shrink_to!(strs, width)
    for i in keys(strs)
        str = strs[i]
        if textwidth(str) > width
            new_str = ""
            w = 0
            for c in str
                new_w = textwidth(c)
                if new_w + w < width
                    new_str = new_str * c
                    w += new_w
                else
                    new_str = new_str * "…"
                    break
                end
            end
            strs[i] = new_str
        end
    end
end

# TODO fix `repr`


text:
# represent a sum of tensor networks
struct NetworkSum
  nodes::Array
end

# TODO: add caching intermediates here
function run(net_sum::NetworkSum, feed_dict::Dict; kwargs...)
  return ITensors.sum(compute_graph(net_sum.nodes, feed_dict; kwargs...); kwargs...)
end

inner(net_sum::NetworkSum, n2) = NetworkSum([inner(n1, n2) for n1 in net_sum.nodes])

Base.push!(net_sum::NetworkSum, node) = Base.push!(net_sum.nodes, node)


text:
"""
A linear combination of Laguerre functions.
"""
struct LaguerreFunction{T} <: AbstractLaguerre{T}
    # representation of a linear Combination of Laguerre function
    # it's just Laguerre polynomial multiplied by exp(-x/2)
    orders::Vector{Int} # order of Laguerre polynomials.
    coeffs::Vector{T} # coefficients of Laguerre polynomials.
    len::Int # length of orders
    function LaguerreFunction(orders::Vector{Int}, coeffs::Vector{T}) where T
        length(orders) == length(coeffs) || throw(DimensionMismatch("orders and coeffs are of different length."))
        all(orders .>= 0) || throw(DomainError("negative orders not allowed."))
        # arrange in a ascending order.
        perm = sortperm(orders)
        orders = orders[perm]
        coeffs = coeffs[perm]
        return new{T}(orders, coeffs, length(orders))
    end
end

function LaguerreFunction(order::Int, coef::Number)
    # handle a single Laguerre Function.
    return LaguerreFunction([order], [coef])
end

function LaguerreFunction(order::Int)
    # A single Laguerre Function with coefficient 1.
    return LaguerreFunction(order, 1)
end

function Base.:*(la::LaguerrePolynomial{T}, lf::LaguerreFunction{S}) where {T,S}
    # A product of LaguerrePolynomial and LaguerreFunction can be represented by
    # LaguerreFunctions.
    # construct a corresponding LaguerrePolynomial.
    la2 = LaguerrePolynomial(lf.orders, lf.coeffs)
    prod = la * la2
    return LaguerreFunction(prod.orders, prod.coeffs)
end

Base.:*(lf::LaguerreFunction{T}, la::LaguerrePolynomial{S}) where {T,S} = la * lf


text:
function load_data(obj::DataLoader, arg0::InputStream, arg1::JString)
    return jcall(obj, "loadData", void, (InputStream, JString), arg0, arg1)
end

function still_accepts_data(obj::DataLoader)
    return jcall(obj, "stillAcceptsData", jboolean, ())
end


text:
function addPackedVariable(navAbilityClient::NavAbilityClient, client::Client, variable)::String
    response = navAbilityClient.mutate(MutationOptions(
        "addVariable",
        MUTATION_ADDVARIABLE,
        Dict(
            "variable" => Dict(
                "client" => client,
                "packedData" => json(variable)
            )
        )
    )) |> fetch
    rootData = JSON.parse(response.Data)
    if haskey(rootData, "errors")
        throw("Error: $(rootData["errors"])")
    end
    data = get(rootData,"data",nothing)
    if data === nothing return "Error" end
    addVariable = get(data,"addVariable","Error")
    return addVariable
end

function addVariable(navAbilityClient::NavAbilityClient, client::Client, variable::Variable)
    return @async addPackedVariable(navAbilityClient, client, variable)
end

function getVariableEvent(navAbilityClient::NavAbilityClient, client::Client, label::String)::Dict{String,Any}
    response = navAbilityClient.query(QueryOptions(
        "sdk_get_variable",
        """
            $GQL_FRAGMENT_VARIABLES
            $GQL_GETVARIABLE
        """,
        Dict(
            "label" => label,
            "userId" => client.userId,
            "robotId" => client.robotId,
            "sessionId" => client.sessionId
        )
    )) |> fetch
    rootData = JSON.parse(response.Data)
    if haskey(rootData, "errors")
        throw("Error: $(rootData["errors"])")
    end
    data = get(rootData,"data",nothing)
    if data === nothing return Dict() end
    user = get(data,"users",[])
    if size(user)[1] < 1 return Dict() end
    robots = get(user[1],"robots",[])
    if size(robots)[1] < 1 return Dict() end
    sessions = get(robots[1],"sessions",[])
    if size(sessions)[1] < 1 return Dict() end
    variables = get(sessions[1],"variables",[])
    if size(variables)[1] < 1 return Dict() end
    return variables[1]
end

getVariable(navAbilityClient::NavAbilityClient, client::Client, label::String) = @async getVariableEvent(navAbilityClient, client, label)

function getVariablesEvent(navAbilityClient::NavAbilityClient, client::Client; detail::QueryDetail = SKELETON)::Vector{Dict{String,Any}}
    response = navAbilityClient.query(QueryOptions(
        "sdk_get_variables",
        """
            $GQL_FRAGMENT_VARIABLES
            $GQL_GETVARIABLES
        """,
        Dict(
            "userId" => client.userId,
            "robotId" => client.robotId,
            "sessionId" => client.sessionId,
            "fields_summary" => detail === SUMMARY || detail === FULL,
            "fields_full" => detail === FULL,
        )
    )) |> fetch
    rootData = JSON.parse(response.Data)
    if haskey(rootData, "errors")
        throw("Error: $(rootData["errors"])")
    end
    data = get(rootData,"data",nothing)
    if data === nothing return [] end
    user = get(data,"users",[])
    if size(user)[1] < 1 return [] end
    robots = get(user[1],"robots",[])
    if size(robots)[1] < 1 return [] end
    sessions = get(robots[1],"sessions",[])
    if size(sessions)[1] < 1 return [] end
    return get(sessions[1],"variables",[])
end

getVariables(navAbilityClient::NavAbilityClient, client::Client; detail::QueryDetail = SKELETON) = @async getVariablesEvent(navAbilityClient, client; detail)

function listVariables(navAbilityClient::NavAbilityClient, client::Client)
    @async begin
        variables = getVariables(navAbilityClient,client) |> fetch
        map(v -> v["label"], variables)
    end
end

function ls(navAbilityClient::NavAbilityClient, client::Client)
    return listVariables(navAbilityClient,client)
end


text:
using Pkg;
Pkg.activate(".");
Pkg.instantiate();
using SolidModeling;

Pkg.test();


text:
###############################################################################
#Calculation of the number of the ions given a concentration and the MM of the#
###############################################################################

#
# Ainda preciso checar como fazer para o caso em que carga é diferente de zero
#

# checar diferente métodos para calcular a concentração

export prot_elec

# One IL + Protein + Water
function prot_elec(data; cube=false, charge=false)

  protein = data.protein
  MMP = data.MMP

  cation1 = data.cation
  anion1  = data.anion

  MM1  = data.MM
  c    = data.c

  vprot = (MMP / 6.02e23) * 1e-3

  ### decidir aqui a melhore forma de calcular as caixas

#  if cube==false
#    lx,ly,lz = measure_prot(protein,30)
#  elseif cube==true
    lx = 94
    ly = 94
    lz = 94
#  end

  sides = [lx, ly, lz]

  vtotal = lx*ly*lz*1.e-27                                      # volume in L
  vsol = - vprot + vtotal                                         # solution volume = volume of the box  - volume of the protein / L
  nions = round(Int64, vsol * c * 6.02e23)                      # number of ions
  vions = nions * MM1 * 1.e-3 / (6.02e23)                       # volume occupied by the ions
  nwater = round(Int64,(vsol - vions) * 6.02e23 / (18*1.e-3))   # number of water molecules to fill the box

  if charge==false
    return nions, nwater, sides
  else
    return nions, nwater, charge, sides
  end

end

export sol_elec
# One IL + Water
function sol_elec(data)

  cation1 = data.cation
  anion1  = data.anion

  MM1  = data.MM
  c    = data.c

  lx = 50
  ly = 50
  lz = 50

  sides = [lx, ly, lz]

  vsol   = lx*ly*lz*1.e-27                                        # Volume (total=solution) / L
  nions  = round(Int64, vsol * c * 6.02e23)                      # number of ions
  vions  = nions * MM1 * 1.e-3 / (6.02e23)                       # volume occupied by the ions
  nwater = round(Int64,(vsol - vions) * 6.02e23 / (18*1.e-3))   # number of water molecules to fill the box

  return nions, nwater, sides

end

export prot_melec
# two IL + Protein + Water
function prot_melec(data; cube=false, charge=false)

  protein = data.protein
  MMP = data.MMP

  cation1 = data.cation
  anion1  = data.anion1
  anion2  = data.anion2

  MM1  = data.MM1
  MM2  = data.MM2
  c    = data.c

  vprot = (MMP / 6.02e23) * 1e-3
#  if cube==false
#    lx,ly,lz = measure_prot(protein,30)
#  elseif cube==true
    lx = 94
    ly = 94
    lz = 94
#  end

  sides = [lx, ly, lz]

  vtotal = lx*ly*lz*1.e-27                                      # volume in L
  vsol = - vprot + vtotal                                         # solution volume = volume of the box  - volume of the protein / L
  nions = round(Int64, vsol * c * 6.02e23)                      # number of ions
  vions = nions * MM1 * 1.e-3 / (6.02e23) + nions * MM2 * 1.e-3 / (6.02e23)                       # volume occupied by the ions
  nwater = round(Int64,(vsol - vions) * 6.02e23 / (18*1.e-3))   # number of water molecules to fill the box

  nan = round(Int,(nions/2)) # same number of molecules for each anion = ncation/2
  # adjustment of the nions - If ncat % 2 != 0 -> ncat = ncat + 1 - ncat must be 2 times greater than nan
  ncat = 2 * nan

  if charge==false
    return ncat, nan, nwater, sides
  else
    return ncat, nan, charge, nwater, sides
  end

end

export sol_melec

# two IL +  Water
function sol_melec(data; cube=false, charge=false)
  cation1 = data.cation
  anion1  = data.anion1
  anion2  = data.anion2

  MM1  = data.MM1
  MM2  = data.MM2
  c    = data.c

  lx = 50
  ly = 50
  lz = 50

  sides = [lx, ly, lz]

  vsol = lx*ly*lz*1.e-27                                      # volume in L
  nions = round(Int64, vsol * c * 6.02e23)                      # number of ions
  vions = nions * MM1 * 1.e-3 / (6.02e23) + nions * MM2 * 1.e-3 / (6.02e23)                       # volume occupied by the ions
  nwater = round(Int64,(vsol - vions) * 6.02e23 / (18*1.e-3))   # number of water molecules to fill the box

  ncat = nions
  nan = round(Int,(nions/2)) # same number of molecules for each anion = ncation/2

  return ncat, nan, nwater, sides

end


text:
#ubox method.

"""
    splitintwo(x::Valid) takes a valid and bifurcates it into a upper and lower
    half, based on their posit representation.
"""
function splitintwo(x::Valid{N,ES}) where {N, ES}
    if (x.upper == next(x.lower))
        tile(x.lower), tile(x.upper)
    elseif (isallreals(x))
        Vnum{N,ES}(Inf) → maxneg(Vnum{N,ES}), zero(Vnum{N,ES}) → maxpos(Vnum{N,ES})
    elseif (roundsinf(x))
        if !isfinite(x.lower)
            Valid{N,ES}(Inf), minneg(Vnum{N,ES}) → x.upper
        elseif !isfinite(x.upper)
            Valid{N,ES}(Inf), x.lower → maxpos(Vnum{N,ES})
        else
            x.lower → Vnum{N,ES}(Inf), minneg(Vnum{N,ES}) → x.upper
        end
    else
        middle = reinterpret(Vnum{N,ES}, (@u x.lower) >> 1 + (@u x.upper) >> 1) |> __round
        x.lower → middle, next(middle) → x.upper
    end
end

function splitinthree(x::Valid{N,ES}) where {N, ES}
    if (x.upper == next(next(x.lower)))
        tile(x.lower), x.lower |> next |> tile, tile(x.upper)
    elseif (isallreals(x))
        Vnum{N,ES}(Inf) → maxneg(Vnum{N,ES}), Vnum{N,ES} |> zero |> tile, minpos(Vnum{N,ES}) → maxpos(Vnum{N,ES})
    elseif (roundsinf(x))
        if !isfinite(x.lower)
            (Valid{N,ES}(Inf), splitintwo(minneg(Vnum{N,ES}) → x.upper)...)
        elseif !isfinite(x.upper)
            (splitintwo(x.lower → maxpos(Vnum{N,ES}))..., Valid{N,ES}(Inf), )
        else
            x.lower → maxneg(Vnum{N,ES}), Inf |> Vnum{N,ES} |> tile, minneg(Vnum{N,ES}) → x.upper
        end
    else
        middle1 = reinterpret(Vnum{N,ES}, (3 * ((@s x.lower) >> 2) + ((@s x.upper) >> 2))) |> __round
        middle2 = reinterpret(Vnum{N,ES}, (((@s x.lower) >> 2) + 3 * ((@s x.upper) >> 2))) |> __round
        x.lower → middle1, next(middle1) → prev(middle2), middle2 → x.upper
    end
end

"""
    vwidth(x::Valid{N,ES})

    describe the width of the valid number.  This is the number of tiles.
"""
vwidth(x::Valid{N,ES}) where {N,ES} = ((@u x.upper) - (@u x.lower)) >> (64 - N) + 1

function merge_contiguous!(v1::Array{Valid{N,ES},J}, v2::Array{Valid{N,ES},J}) where J where {N,ES}
    dim = 0
    for idx = 1:length(v1)
        dim = idx
        (v1[idx] != v2[idx]) && break
    end
    (dim == 0) && return false
    (next(v1[dim].upper) == v2[dim].lower) || (next(v2[dim].upper) == v1[dim].lower) || return false
    for idx = (dim + 1):length(v1)
        v1[idx] != v2[idx] && return false
    end

    if next(v1[dim].upper) == v2[dim].lower
        v1[dim] = v1[dim].lower → v2[dim].upper
    else
        v1[dim] = v2[dim].lower → v1[dim].upper
    end
    true
end

function coalesce!(a::Vector{Array{Valid{N,ES},J}}, b::Vector{Array{Valid{N,ES},J}}) where J where {N,ES}
    for b_item in b
        coalesce!(a, b_item)
    end
    #coalesce!(a)
    return a
end
function coalesce!(a::Vector{Array{Valid{N,ES},J}}, b::Array{Valid{N,ES},J}) where J where {N,ES}
    merged = false
    for idx = 1:length(a)
        if merge_contiguous!(a[idx],b)
            merged = true
            break
        end
    end
    merged || push!(a, b)
end

function coalesce_2!(a::Vector{Array{Valid{N,ES},J}}) where J where {N,ES}
    #pop the top two values and attempt to coalesce them.
    if length(a) >= 2
        if merge_contiguous!(a[end-1], a[end])
            pop!(a)
        end
    end
end

function coalesce_top!(a::Vector{Array{Valid{N,ES},J}}) where J where {N,ES}
    while length(a) >= 2
        if merge_contiguous!(a[end-1], a[end])
            pop!(a)
        else
            break
        end
    end
end

"""
    ufilter(predicate, dataset::Vector{Any})
"""
function ufilter(pred::Function, x::Valid{N,ES}, ::Type{Val{verbose}} = Val{false}) where {N, ES, verbose}
    if pred(x)
        verbose && println("interval $x checks out.")
        if istile(x)
            return Valid{N,ES}[x]
        else
            (upper, lower) = splitintwo(x)
            vcat(ufilter(pred, lower), ufilter(pred, upper))
        end
    else
        return Valid{N,ES}[]
    end
end
function ufilter(pred::Function, x::Array{Valid{N,ES},J}) where J where {N, ES}
    if pred(x)
        #check to see if we're all tiles.
        mapreduce(istile, &, true, x) && return Array{Valid{N,ES},J}[x]
        #find the thing with the biggest index.
        (_, idx) = x .|> vwidth |> findmax
        upper_array = copy(x)
        lower_array = copy(x)

        lower_value, upper_value = x[idx] |> splitintwo

        lower_array[idx] = lower_value
        upper_array[idx] = upper_value

        coalesce!(ufilter(pred, lower_array), ufilter(pred, upper_array))
    else
        return Array{Valid{N,ES},J}[]
    end
end

function ufilter_bfs(pred::Function, x::Array{Valid{N,ES},J}) where J where {N,ES}
    #consider implementing this as a linked list type for better performance.
    queue = Array{Valid{N,ES},J}[]
    results = Array{Valid{N,ES},J}[]
    unshift!(queue, x)
    iter = 0
    while length(queue) > 0
        #examine the first element in the queue.
        this_box = shift!(queue)
        if pred(this_box)
            if mapreduce(istile, &, true, this_box)
                coalesce!(results, this_box)
            else
                #find the widest index in the array.
                (_, idx) = this_box .|> vwidth |> findmax
                upper_array = copy(this_box)
                lower_array = copy(this_box)

                #split this index in the array into two.
                lower_value, upper_value = this_box[idx] |> splitintwo

                lower_array[idx] = lower_value
                upper_array[idx] = upper_value

                #push both arrays onto the check queue.
                unshift!(queue, lower_array)
                unshift!(queue, upper_array)
            end
        end
    end
    results
end

#a temporary version of this, needs to take into account roundsinf.
function Base.union(a::Valid{N,ES},b::Valid{N,ES}) where {N,ES}
    return min(a.lower, b.lower) → max(a.upper, b.upper)
end

function merge_into!(a::Array{Valid{N,ES},J}, b::Array{Valid{N,ES},J}) where J where {N,ES}
    #just assign when a in a.
    isempty(a[1]) && (a[:] = b[:]; return)
    for idx = 1:length(a)
        a[idx] = union(a[idx], b[idx])
    end
end

function Base.in(a::Valid{N,ES}, b::Valid{N,ES}) where {N,ES}
    return a.lower >= b.lower && a.upper <= b.upper
end

function ufilter_dfs(pred::Function, x::Array{Valid{N,ES},J}) where J where {N,ES}
    #consider implementing this as a linked list type for better performance.
    queue = Array{Valid{N,ES},J}[]
    bbox = [Valid{N,ES}(∅) for idx in 1:length(x)]
    push!(queue, x)
    iter = 0
    while length(queue) > 0
        #examine the first element in the queue.
        this_box = pop!(queue)

        this_box in bbox && continue

        if pred(this_box)
            if mapreduce(istile, &, true, this_box)
                merge_into!(bbox, this_box)
            else
                #find the widest index in the array.
                (_, idx) = this_box .|> vwidth |> findmax

                this_width = this_box[idx] |> vwidth

                if this_width >= 3
                    upper_array = copy(this_box)
                    middle_array = copy(this_box)
                    lower_array = copy(this_box)

                    #split this index in the array into two.
                    lower_value, middle_value, upper_value = this_box[idx] |> splitinthree

                    lower_array[idx] = lower_value
                    middle_array[idx] = middle_value
                    upper_array[idx] = upper_value

                    #push all arrays onto the check queue, but put the middle last.
                    push!(queue, middle_array)
                    push!(queue, lower_array)
                    push!(queue, upper_array)
                else
                    upper_array = copy(this_box)
                    lower_array = copy(this_box)

                    #split this index in the array into two.
                    lower_value, upper_value = this_box[idx] |> splitintwo

                    lower_array[idx] = lower_value
                    upper_array[idx] = upper_value

                    #push both arrays onto the check queue.
                    push!(queue, lower_array)
                    push!(queue, upper_array)
                end
            end
        end
    end
    bbox
end

export ufilter, ufilter_bfs, ufilter_dfs

max_issues_repo_name: MohHizzani/SigmoidNumbers

max_forks_repo_name: MohHizzani/SigmoidNumbers


text:
function legendrepolynomial(x,p)

n = length(x);
v0 = zeros(Float64,n,1);
v1 = ones(Float64,n,1);

x = 2.0*reshape(x,n,1) .- 1.0;

if (p==0)
    f = v1;
    fx = v0;
elseif (p==1)
    f = [v1 x];
    fx = [v0 v1];
elseif (p==2)
    f = [v1 x 0.5*(3*(x.^2) .- 1)];
    fx = [v0 v1 3*x];
elseif (p==3)
    f = [v1 x 0.5*(3*(x.^2) .- 1) 0.5*(5*(x.^3) .- 3*x)];
    fx = [v0 v1 3*x 7.5*(x.^2) .- 1.5];
elseif (p==4)
    f = [v1 x 0.5*(3*(x.^2) .- 1) 0.5*(5*(x.^3) .- 3*x) (1.0/8.0)*(35*(x.^4) .- 30*(x.^2) .+ 3)];
    fx = [v0 v1 3*x (7.5*(x.^2) .- 1.5) (1.0/8.0)*(140*(x.^3) .- 60*x)];
else
    error("polynomial degree must be less than 5");
end

fx = 2.0*fx;

return f, fx;

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
module Monitoring

export get_gpu_state, add_logging_to_pyfile, run_task

function get_gpu_state(gpu_id::Integer)
    i = 0
    output = readstring(open(pipeline(`nvidia-smi`, `cat`))[1])
    for line in split(output, "\n")
        if contains(line, "%")
            if i == gpu_id
                if parse(Int, split(line)[13][1:end-1]) > 5
                    return 0 # Busy
                else
                    return 1 # Free
                end
            else
                i += 1
            end
        end
    end
    error("The selected gpu id wasn't found.")
end

function update_logs(log_pos, file)
    text = join(readlines(seek(open(file), log_pos)), "")
    return (text, log_pos + length(text))
end

function add_logging_to_pyfile(file)
    # This function will append changed stdin to python file.
    # So the file will write not to terminal but to file.
    prepend = "import sys\nimport functools\nsys.stdout = open(\"log.txt\", \"w\")\nprint = functools.partial(print, flush=True)"
    data = join(readlines(open(file, "r")), "")
    if startswith(data, prepend)
        return
    end
    rm(file)
    fw = open(file, "w")
    write(fw, string(prepend,data))
    close(fw)
end

function run_task(command::Cmd)::Base.Process
    return spawn(command)
end

end # module

max_stars_repo_name: Gollor/JuliaAutomation

max_issues_repo_name: Gollor/JuliaAutomation

max_forks_repo_name: Gollor/JuliaAutomation


text:
mutable struct DCModel
    Y::Vector{Int32} # N-by-1, choices of individuals
    X1::Array{Float64,3} # N-J-K1: indiv-prod-char
    X2::Matrix{Float64} # J-K2: product char
    D::Matrix{Float64} # N-L: demographics
    C::Vector{Vector{Int64}} # N: choice set of each individual,
                             # can be actual constraints or random draws to facilitate est

    # Parameters
    hetero_preference::Bool # unobservable hetero preference

    α::Vector{Float64} # K1: coef for X1
    Π::Matrix{Float64} # K2-L: coef for X2 interact with demo
    ξ::Vector{Float64} # J-by-1: product FE, normalize the first commodity to 0
    σ::Vector{Float64} # K2-by-1: s.d. for unobs. hetero. preference for Π

    optResults # optimization results
end

function DCModel(Y, X1::Array{T,3}, X2::Matrix{T}, D::Matrix{T}; hetero_preference) where T
    Y = Int.(Y)
    N = length(Y)
    N2, J, K1 = size(X1)
    J2, K2 = size(X2)
    N3, L = size(D)

    @assert J==J2
    @assert N==N2==N3

    C = Vector{Vector{eltype(Y)}}(undef, N)

    # Initialize coef
    α = zeros(T, K1)
    Π = zeros(T, K2, L)
    ξ = zeros(T, J-1)
    if hetero_preference
        σ = zeros(T, K2)
    else
        σ = zeros(T, 0)
    end

    optResults = nothing

    DCModel(Y,X1,X2,D,C, hetero_preference, α, Π, ξ, σ, optResults)
end

function updateCoef!(model::DCModel)
    @assert !isnothing(model.optResults)
    N,J,K1,K2,L = get_dim(model.X1, model.X2, model.D)

    res = Optim.minimizer(model.optResults)
    model.α = res[1:K1]
    model.Π = reshape(res[K1+1:K1+K2*L], K2, L)
    model.ξ = res[K1+K2*L+1:K1+K2*L+J-1] # ξ is only J-1 dim
    if model.hetero_preference
        model.σ = res[K1+K2*L+J+1:end]
    end
    return nothing
end


text:
module Database

using MySQL

const HOST = "127.0.0.1"
const USER = "root"
const PASS = "root"
const DB = "six_degrees"
const PORT = 8889

const CONN = MySQL.connect(HOST, USER, PASS, db = DB, port = PORT)

export CONN

disconnect() = MySQL.disconnect(CONN)

atexit(disconnect)

end

max_stars_repo_name: PacktPublishing/Julia-Projects

max_issues_repo_name: PacktPublishing/Julia-Projects

max_forks_repo_name: PacktPublishing/Julia-Projects


text:
#######################
# The Lagrange basis
#######################

using Polynomials: Poly, polyint
using BasisFunctions: ChebyshevInterval, PolynomialBasis
using BasisFunctions: hasderivative, hasantiderivative, ordering, support

const LagrangeInterval = ChebyshevInterval
const LagrangeIndex = NativeIndex{:lagrange}

"""
A basis of the Lagrange polynomials `l_i(x) = ∏_(j,i≠j) (x - ξ^j) / (ξ^i - ξ^j)`
on the interval [-1,+1].
"""
struct Lagrange{S,T} <: PolynomialBasis{T}
    n :: Int
    nodes :: ScatteredGrid{T}

    denom::Vector{T}    # denominator of Lagrange polynomials
    diffs::Matrix{T}    # inverse of differences between nodes
    vdminv::Matrix{T}   # inverse Vandermonde matrix

    function Lagrange{S,T}(nodes::ScatteredGrid{T}) where {S,T}
        local p::T

        local ξ = nodes.points
        local n = length(ξ)

        @assert minimum(ξ) ≥ leftendpoint(support(nodes))
        @assert maximum(ξ) ≤ rightendpoint(support(nodes))

        denom = ones(n)
        diffs = zeros(n,n)

        for i in eachindex(ξ)
            p = 1
            for j in eachindex(ξ)
                diffs[i,j] = 1 / (ξ[i] - ξ[j])
                if i ≠ j
                    denom[i] *= diffs[i,j]
                end
            end
        end

        new(n, nodes, denom, diffs, vandermonde_matrix_inverse(ξ))
    end
end

Lagrange{S}(nodes::ScatteredGrid{T}) where {S,T} = Lagrange{S,T}(nodes)
Lagrange(nodes::ScatteredGrid) = Lagrange{:custom}(nodes)
Lagrange{S}(ξ::Vector) where {S} = Lagrange{S}(ScatteredGrid(ξ, LagrangeInterval()))
Lagrange(ξ::Vector) = Lagrange(ScatteredGrid(ξ, LagrangeInterval()))
Lagrange{S}(n::Int) where {S} = Lagrange{S}(get_nodes(Val(S), n))
Lagrange{S,T}(n::Int) where {S,T} = Lagrange{S}(get_nodes(Val(S), n, T))

const LagrangeGaussLegendre = Lagrange{:gauss_legendre}
const LagrangeLobattoLegendre = Lagrange{:lobatto_legendre}
const LagrangeGaussChebyshev = Lagrange{:gauss_chebyshev}
const LagrangeLobattoChebyshev = Lagrange{:lobatto_chebyshev}

# Convenience constructor: map the Lagrange basis to the interval [a,b]
Lagrange{S}(n::Int, a::Number, b::Number) where {S} = rescale(Lagrange{S}(n), a, b)

nodes(b::Lagrange)  = b.nodes.points
nnodes(b::Lagrange) = b.n
degree(b::Lagrange) = b.n-1

BasisFunctions.native_index(b::Lagrange, idx) = LagrangeIndex(idx)
BasisFunctions.linear_index(b::Lagrange, idx) = LagrangeIndex(idx)
BasisFunctions.ordering(b::Lagrange) = Base.OneTo(nnodes(b))
BasisFunctions.support(b::Lagrange) = support(b.nodes)

BasisFunctions.interpolation_grid(b::Lagrange{S,T}) where {S,T} = ScatteredGrid(get_nodes(Val(S), length(b), T), LagrangeInterval())

BasisFunctions.similar(::Lagrange{S}, ::Type{T}, n::Int) where {S,T} = Lagrange{S,T}(n)

Base.size(b::Lagrange) = (b.n,)

function BasisFunctions.unsafe_eval_element(b::Lagrange{S,T}, idx::LagrangeIndex, x::T) where {S,T}
    local y::T = 1
    local ξ = nodes(b)
    for i in 1:length(ξ)
        i ≠ idx ? y *= x - ξ[i] : nothing
    end
    y * b.denom[value(idx)]
end

function BasisFunctions.unsafe_eval_element_derivative(b::Lagrange{S,T}, idx::LagrangeIndex, x::T) where {S,T}
    local y::T = 0
    local z::T
    local ξ = nodes(b)
    local d = b.diffs

    for l in 1:nnodes(b)
        if l ≠ idx
            z = d[idx,l]
            for i in 1:nnodes(b)
                i ≠ idx && i ≠ l ? z *= (x - ξ[i]) * d[idx,i] : nothing
            end
            y += z
        end
    end
    y
end

function unsafe_eval_element_antiderivative(b::Lagrange{S,T}, idx::LagrangeIndex, x::T) where {S,T}
    local y = zero(nodes(b))
    y[value(idx)] = 1
    lint = polyint(Poly(b.vdminv*y))
    return lint(x) - lint(leftendpoint(support(b)))
end


text:
function DiffEqBase.__solve(prob::DiffEqBase.AbstractSteadyStateProblem,
                            alg::SteadyStateDiffEqAlgorithm,args...;
                            abstol=1e-8,kwargs...)

  if typeof(prob.u0) <: Number
    u0 = [prob.u0]
  else
    u0 = vec(deepcopy(prob.u0))
  end

  sizeu = size(prob.u0)
  p = prob.p

  if typeof(prob) <: SteadyStateProblem
    if !isinplace(prob) && (typeof(prob.u0)<:AbstractVector || typeof(prob.u0)<:Number)
      f! = (du,u) -> (du[:] = prob.f(u,p,Inf); nothing)
    elseif !isinplace(prob) && typeof(prob.u0)<:AbstractArray
      f! = (du,u) -> (du[:] = vec(prob.f(reshape(u, sizeu),p,Inf)); nothing)
    elseif typeof(prob.u0)<:AbstractVector
      f! = (du,u) -> (prob.f(du,u,p,Inf); nothing)
    else # Then it's an in-place function on an abstract array
      f! = (du,u) -> (prob.f(reshape(du, sizeu),
                      reshape(u, sizeu),p,Inf);
                      du=vec(du); nothing)
    end
  elseif typeof(prob) <: NonlinearProblem
    if !isinplace(prob) && (typeof(prob.u0)<:AbstractVector || typeof(prob.u0)<:Number)
      f! = (du,u) -> (du[:] = prob.f(u,p); nothing)
    elseif !isinplace(prob) && typeof(prob.u0)<:AbstractArray
      f! = (du,u) -> (du[:] = vec(prob.f(reshape(u, sizeu),p)); nothing)
    elseif typeof(prob.u0)<:AbstractVector
      f! = (du,u) -> (prob.f(du,u,p); nothing)
    else # Then it's an in-place function on an abstract array
      f! = (du,u) -> (prob.f(reshape(du, sizeu),
                      reshape(u, sizeu),p);
                      du=vec(du); nothing)
    end
  end

  # du = similar(u)
  # f = (u) -> (f!(du,u); du) # out-of-place version

  if typeof(alg) <: SSRootfind
    u = alg.nlsolve(f!,u0,abstol)
    resid = similar(u)
    f!(resid,u)
    DiffEqBase.build_solution(prob,alg,u,resid;retcode = :Success)
  else
    error("Algorithm not recognized")
  end
end

function DiffEqBase.__solve(prob::DiffEqBase.AbstractSteadyStateProblem,
                            alg::DynamicSS,args...;save_everystep=false,
                            save_start=false,kwargs...)

  tspan = alg.tspan isa Tuple ? alg.tspan : convert.(eltype(prob.u0),(zero(alg.tspan), alg.tspan))
  if typeof(prob) <: SteadyStateProblem
    f = prob.f
  elseif typeof(prob) <: NonlinearProblem
    if isinplace(prob)
      f = (du,u,p,t) -> prob.f(du,u,p)
    else
      f = (u,p,t) -> prob.f(u,p)
    end
  end

  _prob = ODEProblem(f,prob.u0,tspan,prob.p)
  sol = solve(_prob,alg.alg,args...;kwargs...,
              callback=TerminateSteadyState(alg.abstol,alg.reltol),
              save_everystep=save_everystep,save_start=save_start)
  if isinplace(prob)
    du = similar(sol.u[end])
    f(du, sol.u[end], prob.p, sol.t[end])
  else
    du = f(sol.u[end], prob.p, sol.t[end])
  end
  array_condition() = all(abs(d) <= abstol || abs(d) <= reltol*abs(u)
                          for (d,abstol, reltol, u)
                          in zip(du, Iterators.cycle(alg.abstol), Iterators.cycle(alg.reltol), sol.u[end]))
  broadcast_condition() = all((abs.(du) .<= alg.abstol) .| (abs.(du) .<= alg.reltol .* abs.(sol.u[end])))
  if sol.retcode == :Terminated && (typeof(sol.u[end]) <: Array ? array_condition() : broadcast_condition())
    _sol = DiffEqBase.build_solution(prob,alg,sol.u[end],du;retcode = :Success)
  else
    _sol = DiffEqBase.build_solution(prob,alg,sol.u[end],du;retcode = :Failure)
  end
  _sol
end


text:
struct RollingWindow
    WindowSize::Int
    observations::Int
    maxiter::Int
    skip::Int
end

"""
    RollingWindow(samples::Int,windowsize::Int)

Creates a RollingWindow iterator from a number of `samples` and a static `windowsize`. The iterator can be used
in for loops to iteratively return indices of a dynamic rolling window.
"""
RollingWindow(samples::Int,windowsize::Int) = RollingWindow(windowsize, samples,samples - windowsize + 1, 1)

"""
    RollingWindow(samples::Int,windowsize::Int,skip::Int)

Creates a RollingWindow iterator from a number of `samples` and a static `windowsize` where every iteration `skip` steps are skipped.
The iterator can be used in for loops to iteratively return indices of a dynamic rolling window.
"""
RollingWindow(samples::Int,windowsize::Int,skip::Int) = RollingWindow(windowsize, samples,samples - windowsize + 1, skip)

function Base.iterate( iter::RollingWindow, state = 1 )
    window = 0
    if state <= iter.maxiter
        window = collect( state : ( state + iter.WindowSize - 1) )
    else
        return nothing
    end
    return ( window, state + iter.skip  )
end

mutable struct ewma
    lambda::Float64
    center::Float64
    lastval::Float64
    rv::RunningVar
end

"""
    EWMA(Initial::Float64, Lambda::Float64) = ewma(Lambda, Initial, Initial, RunningVar(Initial))

Constructs an exponentially weighted moving average object from an initial scalar property value `Initial` and
the decay parameter `Lambda`. This defaults the center value to be the initial value.
"""
EWMA(Initial::Float64, Lambda::Float64) = ewma(Lambda, Initial, Initial, RunningVar(Initial))

"""
    EWMA(Initial::Float64, Lambda::Float64) = ewma(Lambda, Initial, Initial, RunningVar(Initial))

Constructs an exponentially weighted moving average object from an vector of scalar property values `Initial` and
the decay parameter `Lambda`. This computes the running statistcs neccesary for creating the EWMA model using the
interval provided and updates the center value to the mean of the provided values.
"""
function EWMA(Initial::Array, Lambda::Float64)
    burnin = EWMA(Initial[1], Lambda)#Call constructor above
    len = length(Initial)
    for sample in 1:len
        burnin(Initial[sample])
    end
    burnin.center = burnin.rv.m.mu
    return burnin
end

"""
    EWMA(P::ewma)(New; train = true)

Provides an EWMA score for a `New` scalar value. If ```train == true``` the model is updated to include this new value.
"""
function (P::ewma)(New; train = true)
    P.lastval = ( P.lambda * New ) + ( 1 - P.lambda ) * P.lastval
    if train == true
        Update!(P.rv, New)
    end
    return P.lastval
end

"""
    ChangeCenter(P::ewma, new::Float64)

This is a convenience function to update the center of a `P` EWMA model, to a `new` scalar value.
"""
function ChangeCenter(P::ewma, new::Float64)
    P.center .= new
end

"""
    Variance(P::ewma)

This function returns the EWMA control variance.
"""
Variance(P::ewma) = (P.lambda / (2.0 - P.lambda) ) * Variance(P.rv)

"""
    Limits(P::ewma; k = 3.0)

This function returns the upper and lower control limits with a `k` span of variance for an EWMA object `P`.
"""
Limits(P::ewma; k = 3.0) = (P.center + (k * sqrt( Variance( P ) ) ), P.center - (k * sqrt( Variance( P ) ) )  )

mutable struct NaiveForecaster
    lastval::Float64
end
"""
    NaiveForecaster( univariate::Array )

Makes a forecaster model that simply predicts the last value it has seen for all future values.
"""
NaiveForecast(univariate::Array) = NaiveForecaster( last( univariate ) )

"""
    update!( model::NaiveForecaster, newdata::Float64 )

Makes a forecast model that simply predicts the last value it has seen for all future values.
"""
update!( model::NaiveForecaster, newdata::Float64 ) = model.lastval .= newdata

"""
    update( model::NaiveForecaster, newdata::Float64 )

returns an updated forecast model that simply predicts the last value it has seen for all future values.
"""
update(model::NaiveForecaster, newdata::Float64) = return NaiveForecaster( newdata )

"""
    (nf::NaiveForecaster)()

Predicts with a naive forecaster model.
"""
(nf::NaiveForecaster)() = nf.lastval

mutable struct SimpleAverage
    runmean::RunningMean
end

"""
    SimpleAverage( univariate )
"""
SimpleAverage(initial::Float64) = SimpleAverage( RunningMean( initial) )

"""
    SimpleAverage( univariate::Array )
"""
function SimpleAverage(univariate::Array)
    rm = RunningMean( univariate[1])
    for i in 2:length(univariate)
        Update!( rm, univariate[ i ] )
    end
    return SimpleAverage( rm )
end

"""
    Update( model::naiveforecaster, newdata::Float64 )
"""
function Update(model::SimpleAverage, newdata::Float64)
    return Update!( model.runmean, newdata )
end

"""
    (sa::SimpleAverage)()

Predicts with a naive forecaster model.
"""
(sa::SimpleAverage)() = sa.runmean.mu

struct EchoStateNetwork
    Winput::Array
    W::Array
    Woutput::Array
    States::Array
    LastState::Array
    L2::Float64
    alpha::Float64
    bias
end

"""
    EchoStateNetwork(X, Y, Reservoir = 1000;
                        L2 = 1e-8, alpha = 0.25, SpectralRadius = 1.00, Sparsity = 0.99, Noise = -1.00,
                        bias = true, burnin = 0)

    Creates an echo state network model/object.
"""
function EchoStateNetwork(X, Y, Reservoir = 1000;
                        L2 = 1e-8, alpha = 0.25, SpectralRadius = 1.00, Sparsity = 0.99, Noise = -1.00,
                        bias = true, burnin = 0)
    X = forceMatrix( X ) ; Y = forceMatrix( Y )#Ensure inputs are Array2's
    (Obs, Vars) = size( X )
    #Allocate memory for the collected states matrix
    States = zeros(bias + Vars + Reservoir, Obs - burnin - 1)
    Predictions = zeros( size( Y ) )
    #Initialize input weight matrix and recurrent weight matrix
    Winput = rand( Reservoir, bias + Vars) .* -2 .+ 1
    #Generate Resevoir matrix with random uniform values from -1 to +1
    W = rand( Reservoir, Reservoir ) .* -2 .+ 1
    #Create a 'mask' from the Bernouli distribution (binomial distribution size = 1)
    Dropout = rbinomial(1.0 - Sparsity, Reservoir, Reservoir )
    #Now we need to introduce sparsity in our reservoir connections based on the mask
    W = W .* Dropout
    #Normalize and set spectral radius of reservoir weight matrix
    W = W .* ( SpectralRadius / abs( LinearAlgebra.eigen( W ).values[ 1 ] ) )
    #Run the reservoir with the data and collect X
    state = zeros(Reservoir, 1)
    #Go through each time step and compute the model
    for t in 1 : (Obs - 1)
        u = bias ? vcat( 1, X[t,:] ) : X[t,:]
        state = ((1.0 - alpha) * state) .+ (alpha * tanh.( Winput * u .+ W * state ))
        #If we are past the burnin number, start saving the model
        if t > burnin
            States[ :, t - burnin ] = vcat(u, state)
        end
    end
    if Noise > 0.0
        VS = size(S)[2]
        for obs in 1:Obs
            mini = reduce( min, States[obs,:] )
            maxi = reduce( max, States[obs,:] )
            States[obs,:] .+= States[obs,:] .* ( rand( 1, VS ) .* Noise .- (Noise/2.0) )
        end
    end
    #Train the output layer from the random reservoir via least squares
    Woutput = Y[ (burnin + 2) : end, :]' * States' * Base.inv( States * States' .+ L2 * LinearAlgebra.Diagonal(ones(bias + Vars + Reservoir) ) )
    return EchoStateNetwork( Winput, W, Woutput, States, state, L2, alpha, bias )
end

"""
    TuneRidge(Y, model::EchoStateNetwork, L2, burnin = 0 )

    This function modifies an ESN model's output weights without reconstructing the dynamic resevoir.
    Basically this is a time saver!
"""
function TuneRidge(Y, model::EchoStateNetwork, L2, burnin = 0 )
    model.L2 .= L2
    model.Woutput .= Y[ (burnin + 2) : end]' * model.States' * Base.inv( model.States * model.States' .+ L2 * diag(ones(bias + Vars + Reservoir) ) )
end

"""
    PredictFn(model::EchoStateNetwork, X; UseLastState = true, burnin = 0)

    This function uses X data to predict Y values from a stored ESN model.
    Burnin can be used here, but it simply rejects the first `N` samples from storage/allocation.
    UseLastState will employ the last known state from the training set as the first state estimate (reccomended)
"""
function PredictFn(model::EchoStateNetwork, X; UseLastState = true, burnin = 0)
    (Obs,Vars) = size(X)
    #Run the trained ESN using last state as starting point
    state = UseLastState ? model.LastState : state = zeros( size( model.Reservoir ) )
    Prediction = zeros( size(model.Woutput)[1], Obs - burnin)

    u = X[1,:]
    for t in 1 : Obs
        u = model.bias ? vcat( 1, X[t,:]) : X[t,:]
        state = ((1.0 - model.alpha ) * state) .+ (model.alpha * tanh.( model.Winput * u .+ model.W * state ) )
        if t > burnin
            Prediction[ : , t - burnin ] = model.Woutput * vcat(u, state)
        end
        if t < Obs
            u = X[t + 1,:] #Load next sample for prediction
        end
    end#End for

    return Prediction
end


text:
module SBML

using SBML_jll, Libdl, Pkg
using SparseArrays

include("structs.jl")
include("version.jl")
include("readsbml.jl")
include("math.jl")
include("utils.jl")

sbml = (sym::Symbol) -> dlsym(SBML_jll.libsbml_handle, sym)

export SBMLVersion,
    readSBML, Model, UnitPart, Species, Reaction, getS, getLBs, getUBs, getOCs

end # module


text:
module QuantumReservoirComputing

using DifferentialEquations, LightGraphs, LinearAlgebra

include("input_handler.jl")
export InputSignal

include("ginzburg_landau.jl")
export GLLayer

end # module


text:
module RunExample

using Test

@time @testset "SubTriangulation" begin
  include("SubTriangulation.jl")

  filename = joinpath(@__DIR__,"..","test/data/47076.stl")
  SubTriangulation.main(filename,n=50)

  filename = SubTriangulation.download(37384)
  SubTriangulation.main(filename,n=50)
  rm(filename)

end

@time @testset "Poisson" begin
  include("Poisson.jl")

  filename = joinpath(@__DIR__,"..","test/data/293137.stl")
  Poisson.main(filename,n=20)
end

end # module


text:
abstract type TimeError <: Exception end

"""
    AmbiguousTimeError(local_datetime, timezone)

The provided local datetime is ambiguous within the specified timezone. Typically occurs on
daylight saving time transitions which "fall back" causing duplicate hour long period.
"""
struct AmbiguousTimeError <: TimeError
    local_dt::DateTime
    timezone::TimeZone
end

function Base.showerror(io::IO, e::AmbiguousTimeError)
    print(io, "AmbiguousTimeError: ")
    print(io, "Local DateTime $(e.local_dt) is ambiguous within $(string(e.timezone))")
end

"""
    NonExistentTimeError(local_datetime, timezone)

The provided local datetime is does not exist within the specified timezone. Typically
occurs on daylight saving time transitions which "spring forward" causing an hour long
period to be skipped.
"""
struct NonExistentTimeError <: TimeError
    local_dt::DateTime
    timezone::TimeZone
end

function Base.showerror(io::IO, e::NonExistentTimeError)
    print(io, "NonExistentTimeError: ")
    print(io, "Local DateTime $(e.local_dt) does not exist within $(string(e.timezone))")
end

"""
    UnhandledTimeError(timezone)

The timezone calculation occurs beyond the last calculated transition.
"""
struct UnhandledTimeError <: TimeError
    tz::VariableTimeZone
end

function Base.showerror(io::IO, e::UnhandledTimeError)
    print(io, "UnhandledTimeError: ")
    print(io, "TimeZone $(string(e.tz)) does not handle dates on or after $(e.tz.cutoff) UTC")
end


text:
# =============================================================================
# Max volume subject to superficial Area
#  ___
# |___|
#
# Luis Lopez Diaz - llopeza@uninorte.edu.co
# =============================================================================
using JuMP, Ipopt, Gurobi;
# Input data
As = 73;                                        # Max superficial area
X = 8;                                          # Max width
Y = 20;                                         # Max lenght
ml = Model(with_optimizer(Ipopt.Optimizer));    # Optimizer
@variable(ml, 0 ≤ x);                           # x dimension
@variable(ml, 0 ≤ y);                           # y dimension
@variable(ml, 0 ≤ z);                           # z dimension
@NLconstraint(ml, 2*(x*y + x*z + y*z) == As);   # constraint
@NLobjective(ml, Max, x*y*z);                   # Objetive function
optimize!(ml);                                  # Optimize
println("Length: ", value(x))
println("Width : ", value(y))
println("Height: ", value(z))


text:
includet("simulation_model.jl")
using StatsBase, GLM, DataFrames, LinearAlgebra, Plots

function estimateOptimalPolicy()
    players = SCTypeData(10000)
    d = 10000
    βs = rand(Uniform(0.5, 1.1), d)

    μs = zeros(d, 2)
    σs = zeros(d, 4)

    for i in 1:d
        data = ObservedData(players, [0, βs[i]])
        σs[i, :] = reshape(cov([data.y data.x]), (1, 4))
        μs[i, 1] = mean(data.y)
        μs[i, 2] = mean(data.x)
     end

    display(plot(βs, σs[:, 4], seriestype=:scatter))
    μ_coef = [coef(lm(@formula(μ ~ β + β^2), DataFrame(μ=μs[:, i], β=βs))) for i in 1:2]
    σ_coef = [coef(lm(@formula(σ ~ β + β^2), DataFrame(σ=σs[:, i], β=βs))) for i in 1:4]
    function loss(β)
        μb = [ dot(μ_coef[i], [1, β, β^2]) for i in 1:2]
        σb = reshape([ dot(σ_coef[i],[1, β, β^2]) for i in 1:4], (2,2))
        samples = rand(MvNormal(μb, σb), 1000000)
        β0 = mean(samples[1, :] .- samples[2, :].*β)
        return mean((samples[1, :] .- β0 .- samples[2, :]*β).^2)
    end

    βstar = Optim.optimize(loss, 0.0, 1.0).minimizer
    return βstar
end

#then we estimate a linear function for mean and variance of the x variables

# then we optimize the loss function when we integrate over a normal distribution

max_stars_repo_name: evanmunro/personalized-policy

max_issues_repo_name: evanmunro/personalized-policy

max_forks_repo_name: evanmunro/personalized-policy


text:
using Huggett890, Test

function euler_test()
    @testset "Euler" begin
        m = init_model();
        t = 2;  # test other ages +++++
        k = 5.0;
        eIdx = 3;
        kPrimeTomorrow = init_test_worker_solution(eIdx + 1);

        kPrime = Huggett890.solve_one_point(m, t, k, eIdx, kPrimeTomorrow);
        dev = Huggett890.euler_dev_one_point(m, t, k, eIdx, kPrime, kPrimeTomorrow);
        @test kPrime >= Huggett890.kprime_min(m.kGrid, t);
        @test kPrime <= Huggett890.kprime_max(m.kGrid, t);
        @test abs(dev) < 1e-4
        # dev2 = Huggett890.euler_dev_one_point2(m, t, k, kPrime, kPrimeTomorrow);
        # @test abs(dev2) < 1e-4

        # Euler needs more direct test +++++
    end
end

function solve_worker_test()
    @testset "Solve worker" begin
        m = init_model();

        # Just a syntax test. Need to solve t and t+1 for substantive testing.
        t = 3; # test other ages +++++
        k = 5.0;
        eIdx = 2;
        nEff = 4;
        kPrimeTomorrow = init_test_worker_solution(nEff);
        sol = Huggett890.solve_one_period(m, kPrimeTomorrow, t);

        Huggett890.solve_worker(m);
        # +++++  @test Huggett890.check_retired_solution(m, sol)
    end
end

function solve_period_test()
    @testset "One period" begin
        m = init_model();
        t = 2;
        k = 5.0;
        nEff = 4;
        eIdx = 2;
        kPrimeTomorrow = init_test_worker_solution(nEff);
        kPrimeFct, cFct = hg.solve_one_efficiency(m, t, eIdx, kPrimeTomorrow);
        kPrime = kPrimeFct(k);
        eMu = hg.e_uprime_c(m, t, eIdx, kPrime, kPrimeTomorrow);
        c = cFct(k);
        muC = hg.marg_utility(m.util, c);
        dev2 = muC / eMu - hg.betar(m);
        # Surprisingly imprecise +++
        @test abs(dev2) < 1e-2

        eff = hg.efficiency(m.endow, eIdx);
        kPrimeMin, kPrimeMax = hg.kprime_range(m, t, k, eff);
        @test kPrimeMax > kPrimeMin
        kPrime2 = hg.solve_one_point(m, t, k, eIdx, kPrimeTomorrow);
        @test isapprox(kPrime2, kPrime)
    end
end

function solve_test()
    @testset "Solve" begin
        m = init_model();
        sol = solve_worker(m);
        @test sol isa Solution
    end
end

@testset "Worker" begin
    euler_test();
    solve_period_test();
    solve_worker_test();
    solve_test();
end

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