text:
function morsefreq(c::CWT{W,T,Morse,N}) where {W,T,N}

    # measures of frequency for generalized Morse wavelet. [with F. Rekibi]
    # the output returns the modal or peak

    # For be=0, the "wavelet" becomes an analytic lowpass filter

    # Lilly and Olhede (2009).  Higher-order properties of analytic wavelets.
    # IEEE Trans. Sig. Proc., 57 (1), 146--160.

    ga = c.waveType.ga;
    be = c.waveType.be;

    fm = exp.((log.(be) - log.(ga)) ./ ga);

    if sum(be.==0) != 0 && size(fm) == ()
        fm = (log(2))^(1 / ga);
    elseif sum(be.==0) != 0 && size(fm) != ()
        fm[be.==0] = (log(2))^(1 / ga[be.==0]);
    end

    fm = fm / (2 * pi);

    return fm

end

"""
    nOctaves, totalWavelets, sRanges, sWidth = getNWavelets(n1,c)

utility for understanding the spacing of the wavelets. `sRanges` is a list of
the s values used in each octave. sWidth is a list of the corresponding
variance adjustments
"""
function getNWavelets(n1, c::CWT)
    nOctaves = getNOctaves(n1, c)
    n,nSpace = setn(n1,c)
    isAve = !(typeof(c.averagingType) <: NoAve)
    if nOctaves ≤ c.averagingLength + getMinScaling(c) # there isn't enough space for anything but averaging.
        totalWavelets = isAve
        sRanges = [1.0]
        sWidth = [1.0]
        return nOctaves, totalWavelets, sRanges, sWidth
    end
    sRange = 2 .^ (polySpacing(nOctaves, c))
    totalWavelets = round(Int, length(sRange) + isAve)
    sWidth = varianceAdjust(c,totalWavelets,nOctaves)
    return nOctaves, totalWavelets, sRange, sWidth
end

"""
different wavelet familes need to end at a different number of octaves because they have different tail behavior
"""
getNOctaves(n1,c::CWT{W,T, M, N}) where {W, T, N, M} = log2(n1>>1+1) + c.extraOctaves
# choose the number of octaves so the last mean, which is at s*σ[1]
# is 3 standard devations away from the end
getNOctaves(n1,c::CWT{W,T, Morlet, N}) where {W, T, N} = log2((n1>>1+1)/(c.σ[1]+3)) + c.extraOctaves
getNOctaves(n1,c::CWT{W,T, <:Paul, N}) where {W, T, N} = log2((n1>>1+1)/(2c.α+5)) + c.extraOctaves
# choose the number of octaves so the last mean is 5 standard deviations from the end
function getNOctaves(n1,c::CWT{W,T, <:Dog, N}) where {W, T, N}
    μ = getMean(c)
    σ = getStd(c)
    log2(n1>>1/(μ+5σ)) + c.extraOctaves
end
# choose the number of octaves so the smallest support is twice the qmf
getNOctaves(n1,c::CWT{W,T, <:ContOrtho, N}) where {W, T, N} = log2(n1) - 2 - log2(length(qmf(c.waveType))) + c.extraOctaves
getNOctaves(n1,c::CWT{W,T, Morse, N}) where {W, T, N} = log2((n1>>1+1)/(morsefreq(c) + 1)) + c.extraOctaves
# getNOctaves(n1,c::CWT{W,T, Morse, N}) where {W, T, N} = 4 + c.extraOctaves

"""
As with the last octave, different wavelet families have different space decay rates, and in the case of symmetric or zero padding we don't want wavelets that bleed across the boundary from the center.
"""
getMinScaling(c::CWT{W,T,M,N}) where {W,T,N,M} = 0 # by default all scales are allowed (all of the orthogonal transforms)
getMinScaling(c::CWT{W,T,<:Morlet,N}) where {W,T,N} = 1/(c.β)^.8 # morlet is slightly too large at the boundaries by default
getMinScaling(c::CWT{W,T,<:Paul,N}) where {W,T,N} = 2/(2c.α+1) # Paul presents some difficulties, as the decay changes quickly (like 1/t^(α+1))
getMinScaling(c::CWT{W,T,<:Dog,N}) where {W,T,N} = 2 # like morlet, the decay for Dog is exponential and consistent across derivatives
# getMinScaling(c::CWT{W,T,<:Morse,N}) where {W,T,N,M} = log2(2 * morsefreq(c))
getMinScaling(c::CWT{W,T,<:Morse,N}) where {W,T,N,M} = 5 * morsefreq(c)

function varianceAdjust(this::CWT{W,T, M, N}, totalWavelets, nOct) where {W,T,N, M}
    # increases the width of the wavelets by σ[i] = (1+a(total-i)ᴾ)σₛ
    # which is a polynomial of order p going from 1 at the highest frequency to
    # sqrt(p) at the lowest
    β = this.β
    x = exp(nOct)/(100+exp(nOct)) # the fewer octaves there are, the smaller the width adjustment we need
    a = (β^x-1)/(totalWavelets-1)^β
    sWidth = 1 .+a .*(totalWavelets .- (1:totalWavelets)).^β
    if any(isnan.(sWidth))
        return [1]
    end
    return sWidth
end

function polySpacing(nOct, c)
    a = max(getMinScaling(c) + c.averagingLength, 0)
    O = nOct
    β = c.β
    Q = c.Q
    if O ≤ a
        # averagingLength and the min scaling are too large for anything to be done
        return [1.0]
    elseif 0 < O - a ≤ 1
        # there's only one octave, just return the linear rate of Q
        return range(a, O, length=1+round(Int, Q))
    end
    # x is the index, y is the scale
    # y= b*x^(1/β), solve for a and b with
    # (x₀,y₀)=(1, aveLength + minScale)
    # dy/dx = 1/s, so the Quality factor gives the slope at the last frequency
    b = (β/Q)^(1/β) * (O)^((β-1)/β)
    # the point x so that the second condition holds
    lastWavelet = Q * (O)/β
    # the point so that the first wavelet is at a
    firstWavelet = (a/b)^β
    # step size so that there are actually Q wavelets in the last octave
    startOfLastOctave = ((nOct-1)/b)^β
    stepSize = (lastWavelet - startOfLastOctave)/Q
    samplePoints = range(firstWavelet, stop=lastWavelet,
                         step = stepSize)#, length = round(Int, O * Q^(1/p)))
    return b .* (samplePoints).^(1/β)
end

# a utility to just get the start, stop, and step size used in polySpacing. Only used for explanatory purposes
function genSamplePoints(nOct, c)
    a = getMinScaling(c) + c.averagingLength
    O = nOct
    β = c.β
    Q = c.Q
    if O ≤ a
        # averagingLength and the min scaling are too large for anything to be done
        return [1.0]
    elseif 0 < O - a ≤ 1
        # there's only one octave, just return the linear rate of Q
        return range(a, O, length=1+round(Int, Q))
    end
    # x is the index, y is the scale
    # y= b*x^(1/β), solve for a and b with
    # (x₀,y₀)=(1, aveLength + minScale)
    # dy/dx = 1/s, so the Quality factor gives the slope at the last frequency
    b = (β/Q)^(1/β) * (O)^((β-1)/β)
    # the point x so that the second condition holds
    lastWavelet = Q * (O)/β
    # the point so that the first wavelet is at a
    firstWavelet = (a/b)^β
    # step size so that there are actually Q wavelets in the last octave
    startOfLastOctave = ((nOct-1)/b)^β
    stepSize = (lastWavelet - startOfLastOctave)/Q
    firstWavelet, lastWavelet, stepSize
end

"""

adjust the length of the storage based on the boundary conditions
"""
function setn(n1, c)
    if boundaryType(c) <: ZPBoundary
        base2 = ceil(Int,log2(n1));   # power of 2 nearest to n1
        nSpace = 2^(base2)
        n = nSpace>>1 + 1
    elseif boundaryType(c) <: SymBoundary
        # n1+1 rather than just n1 because this is going to be used in an rfft
        # for real data
        n = n1 + 1
        nSpace = 2*n1
    elseif boundaryType(c) <: PerBoundary
        n= n1>>1 + 1
        nSpace = n1
    else
        error("No such boundary as $(boundaryType(c))")
    end
    return n, nSpace
end

# computing averaging function utils

# this makes sure the averaging function is real and positive
adjust(c::CWT) = 1
adjust(c::CWT{W, T, <:Dog}) where {W,T} = 1/(im)^(c.α)

function locationShift(c::CWT{W, T, <:Morlet, N}, s, ω, sWidth) where {W,T,N}
        s0 = 3*c.σ[1] *s*sWidth
        ω_shift = ω .+ c.σ[1] * s0
    return (s0, ω_shift)
end

function locationShift(c::CWT{W, T, <:Dog, N}, s, ω, sWidth) where {W,T,N}
    μNoS = getMean(c)
    μLast = μNoS * (2s)
    s0 = μLast / 2 / getStd(c)
    μ = μNoS * s0
    ω_shift = ω .+ μ
    return (s0, ω_shift)
end

function locationShift(c::CWT{W, T, <:Morse, N}, s, ω, sWidth) where {W,T,N}
        # s0 = c.waveType.cf * s * sWidth
        s0 = morsefreq(c) * s * sWidth
        # ω_shift = ω .+ c.waveType.cf * s0
        ω_shift = ω .+ s0
    return (s0, ω_shift)
end

"""
get the mean of the mother wavelet where s=1
"""
function getMean(c::CWT{W, T, <:Dog},s=1) where {W,T}
    Gm1 = gamma((c.α+1)/2)
    Gm2 = gamma((c.α+2)/2)
    Gm2 / Gm1 * sqrt(2) * s^2
end
getMean(c::CWT{W,T,<:Paul}, s=1) where {W,T} = (c.α+1)*s
function getMean(c::CWT{W, T, <:Morlet},s=1) where {W,T}
    return s*c.σ[1]
end
function getMean(c::CWT{W, T, <:Morse},s=1) where {W,T}
    #return s*c.waveType.cf
    return s*morsefreq(c)
end
"""
    getStd(c::CWT{W, T, <:Dog}, s=1) where {W,T}
get the standard deviation of the mother wavelet
"""
getStd(c::CWT{W, T, <:Dog}, s=1) where {W,T} = sqrt(c.α + 1 - getMean(c)^2)*s^(3/2)
getStd(c::CWT{W, T, <:Paul},s=1) where {W, T} = sqrt((c.α+2)*(c.α+1)) * s
getStd(c::CWT{W, T, <:Morlet},s=1) where {W,T} = (s * c.β^.8)/sqrt(2)
getStd(c::CWT{W, T, <:Morse},s=1) where {W,T} = (s * c.β^.8)/sqrt(2)

function locationShift(c::CWT{W, T, <:Paul, N}, s, ω, sWidth) where {W,T,N}
    s0 = s*sqrt(c.α+1)
    ω_shift = ω .+ (c.α .+ 1) * s0
    return (s0, ω_shift)
end

function getUpperBound(c::CWT{W, T, <:Morlet, N}, s) where {W,T,N}
    return c.σ[1] * s
end

function getUpperBound(c::CWT{W, T, <:Dog, N}, s) where {W,T,N}
    return sqrt(2)* s * gamma((c.α+2)/2) / gamma((c.α+1)/2)
end

function getUpperBound(c::CWT{W, T, <:Paul, N}, s) where {W,T,N}
    return (c.α + 1) * s
end

"""
    arrayOfFreqs = getMeanFreq(Ŵ, δt=1000)
Calculate each of the mean frequencies of a collection of analytic or real wavelets Ŵ. assumes a sampling rate of 2kHz, so the maximum frequency is 1kHz. Change δt to adjust to your problem.
"""
function getMeanFreq(Ŵ, δt=2000)
    eachNorm = [norm(w,1) for w in eachslice(Ŵ,dims=ndims(Ŵ))]'
    freqs = range(0,δt/2, length= size(Ŵ,1))
    dropdims(sum(freqs .* abs.(Ŵ), dims=1) ./ eachNorm,dims=1)
end

# orthogonal wavelet utils
function getContWaveFromOrtho(c,N)
    depth, ψLen = calcDepth(c,N)
    φ, ψ, t = makewavelet(c.waveType.o, depth)
    return φ, ψ, ψLen
end

"""
    nIters, sigLength = calcDepth(w,N)
given a `CWT` with an orthogonal wavelet type, calculate the depth `nIters`
necessary to get a resulting wavelet longer than `N`. `sigLength` is the
resulting length.
"""
function calcDepth(w, N)
    origLen = length(qmf(w.waveType))
    netLen = origLen
    nIters = 0
    if N < netLen
        return nIters
    end
    while N > netLen - 2^nIters + 1
        netLen = nextLength(netLen, origLen)
        nIters += 1
    end
    return (nIters, netLen - 2^nIters + 1)
end
# the actual length calculator, as this is recursive
nextLength(curr,orig) = curr*2 + orig - 1 # upsample then convolve
# pad a vector `v` to have length `N`. Or chop off from both sides if its too long
function padTo(v, N)
    nDiff = N - length(v)
    if nDiff > 0 # there should be more entries than there are
        cat(v, zeros(eltype(v), max(N-length(v),0)), dims = 1)
    elseif nDiff ==0
        return v
    else # there are too many, so chop off half the difference from each side
       inds = ceil(Int, 1-nDiff/2):length(v)+ceil(Int,nDiff/2)
       return v[inds]
    end
end

# create interpolater for the orthogonal cases
genInterp(ψ) = interpolate(ψ, BSpline(Quadratic(Reflect(OnGrid()))))

"""
    β = computeDualWeights(Ŵ)
compute the weight given to each wavelet so that in the Fourier domain, the sum across wavelets is as close to 1 at every frequency below the peak of the last wavelet (that is `β' .*abs.(Ŵ[1:lastFreq,:]) ≈ ones(lastFreq)` in an ℓ^2 sense)
"""
function computeDualWeights(Ŵ, wav)
    @views lastReasonableFreq = computeLastFreq(Ŵ[:,end], wav)
    Wdag = pinv(Ŵ[1:lastReasonableFreq,:])
    β = conj.(Wdag * ones(size(Wdag,2)))'
    return β
end
# function computeDualWeights(Ŵ, wav::CWT{W,T,<:ContOrtho,N}, n,naive=true) where {W,T,N}
#     β = computeNaiveDualWeights(Ŵ, wav, n)
#     return β
# end
computeLastFreq(ŵ, wav) = length(ŵ) - 3 # we have constructed the wavelets so the last frequency is <1%, so trying to get it to 1 will fail
computeLastFreq(ŵ, wav::CWT{W,T,<:Morlet}) where {W,T} = findlast(abs.(ŵ) .> 1/4 * maximum(abs.(ŵ))) # the gaussian decay makes this way too large by the end

"""
    computeNaiveDualWeights(Ŵ, wav, n1)
Compute the dual weights using the scaling amounts and the normalization power p. Primarily used when the least squares version is poorly constructed. When the least squares version doesn't perform well, this is also likely to have poor reconstruction, but it won't give extremely negative or oscillatory weights like the least squares version.
"""
function computeNaiveDualWeights(Ŵ, wav, n1)
    _, _, sRange, _ = ContinuousWavelets.getNWavelets(n1, wav)
    isAve = !(wav.averagingType isa NoAve)
    p = wav.p
    # every wavelet is multiplied by 1/s^(1/p), so we need to multiply by that to normalize in frequency. The extra +1 is to convert from uniform in scale to uniform in frequency
    sToTheP = sRange' .^ (1/p+1)
    if isAve && size(Ŵ,2)>1
        # averaging needs to have maximum value of 1
        aveMultiplier = 1/norm(Ŵ[:,1],Inf)
        dual = sum(Ŵ[:,2:end] .* sToTheP, dims=2)
        renormalize = norm(dual,Inf)
        dual ./=renormalize
        #dual += Ŵ[:,1] .* sToTheP[1]
        β = cat(aveMultiplier, 2 .* sToTheP./renormalize..., dims=2)
    elseif size(Ŵ,2)>1
        dual = sum(Ŵ .* sToTheP, dims=2)
        β = sToTheP ./ norm(dual, Inf)
        β[2:end] .* 2
    else
        # there's only the averaging
        β = [1/norm(Ŵ[:,1],Inf)]
    end
    return β
end

"""
    dualCover, dualNorm = getDualCoverage(n,cWav, invType)
get the sum of the weights and its deviation from 1
"""
function getDualCoverage(n,cWav, invType)
    Ŵ = computeWavelets(n,cWav)[1]
    if invType isa DualFrames
        dualCover = sum(conj.(Ŵ) .* explicitConstruction(Ŵ),dims=2)
        return dualCover, norm(dualCover .- 1)
    elseif invType isa PenroseDelta
        β = computeDualWeights(Ŵ, cWav)
    elseif invType isa NaiveDelta
        β = computeNaiveDualWeights(Ŵ, cWav, n)
    end
    dualCover = sum(conj.(β) .* Ŵ, dims=2)
    dualCover, norm(dualCover .- 1)
end

function dualDeviance(n, cWav, naive)
    getDualCoverage(n,cWav,naive)
end

"""
Explicitly construct the canonical dual frame elements for a translation invariant frame. This is likely to end poorly due to the low representation at high frequencies.
"""
function explicitConstruction(Ŵ)
    Ŵdual = conj.(Ŵ ./ [norm(ŵ)^2 for ŵ in eachslice(Ŵ,dims=1)])
end


text:
using Plots, LaTeXStrings,JLD
function neighbor_checking(network_,dim,i,j)
    neighbors_list=[]
    if i+1!=dim+1
        push!(neighbors_list,network_[i+1,j])
    end
    if i-1!=0
        push!(neighbors_list,network_[i-1,j])
    end
    if j+1!=dim+1
        push!(neighbors_list,network_[i,j+1])
    end
    if j-1!=0
        push!(neighbors_list, network_[i,j-1])
    end
    return neighbors_list
end
function on_or_blocked(network_,dim,i,j,p)
    if p>rand()
        network_[i,j]=1
    else
        network_[i,j]=-1
    end
    return network_
end
p=0.55
dim=50
network_=zeros(Int,dim,dim)
network_[rand(1:dim),rand(1:dim)]=1 #the first entry
for step in 1:(dim^2)
    operation_list=[]
    for i in 1:dim
        for j in 1:dim
            if network_[i,j]==0 && (1 in neighbor_checking(network_,dim,i,j))
                push!(operation_list,[i,j])
            end
        end
    end
    if length(operation_list)>0
        for i in operation_list
            network_=on_or_blocked(network_,dim,i[1],i[2],p)
        end
    elseif length(operation_list)==0
        break
    end
end
heatmap(network_)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter4\\Fig\\4.7_p=0.5.png")

ext: jl


text:
using Pluto
import Pluto: run_reactive!, WorkspaceManager

Pluto.set_ENV_defaults()

@testset "Workspace manager" begin
# basic functionality is already tested by the reactivity tests

    @testset "Multiple notebooks" begin
        fakeclientA = Client(:fakeA, nothing)
        fakeclientB = Client(:fakeB, nothing)
        Pluto.connectedclients[fakeclientA.id] = fakeclientA
        Pluto.connectedclients[fakeclientB.id] = fakeclientB

        notebookA = Notebook(joinpath(tempdir(), "test.jl"), [
            Cell("x = 3")
        ])
        fakeclientA.connected_notebook = notebookA

        notebookB = Notebook(joinpath(tempdir(), "test.jl"), [
            Cell("x")
        ])
        fakeclientB.connected_notebook = notebookB

        @test_nowarn run_reactive!(notebookA, notebookA.cells[1])
        @test_nowarn run_reactive!(notebookB, notebookB.cells[1])

        @test notebookB.cells[1].errored == true
    end
end


text: step_up() = while !step() step_up() end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function interpLinear(xpt::AbstractArray,ypt::AbstractArray,x)
    if x >= minimum(xpt) && x <= maximum(xpt)
        idx = findfirst(t->t>= x, xpt) -1
        idx == 0 ? idx += 1 : nothing
        θ = (x - xpt[idx])/ (xpt[idx+1] - xpt[idx])
        (1-θ)*ypt[idx] + θ*ypt[idx+1]
    else
        0
    end
end

function odeTK(du, u, p,t)
    du[1] = p[3] * (interpLinear(p[1],p[2],t)-u[1])
end

function odeSD(du, u, p,t)
    du[1] = p[3] * (interpLinear(p[1],p[2],t)-u[1])
    du[2] = p[4] + maximum([u[1] - p[5], 0]) * p[6]
end

function logLogisticLaw(x, α, β)
    1 /(1 +(x/α)^(-β))
end

"""
    Solve Toxicokinetics Models with a single parameter

    \$\\frac{dx(t)}{dt} = k_d (conc(t) - x(t))\$

**Fields**
    - `tps` -- time vector
    - `conc` -- exposure vector
    - `kd` -- parameter, a scalar

Return a Tuple with `Array{Float64,1}` of th same length as `tps`

# Example
```julia
julia> myTK = runTK([0,1,2,3], [0,1,20,2], 0.5);
julia> myTK.TK
4-element Array{Float64,1}:
 0.0
 0.21800446293645104
 4.570824232907027
 6.801199145490319
```
"""
function runTK(tps, conc, kd)
    prob = DiffEqBase.ODEProblem(odeTK,[0.0],(0.0, maximum(tps)),[tps, conc, kd])
    _saveat = tps === nothing ? Float64[] : tps
    sol = DiffEqBase.solve(prob ; saveat = _saveat)
    TK = [sol.u[j][1] for j in 1:length(sol)]
    return (TK=TK,)
end

"""
    Solve Toxicokinetics Models with MCMC parameters

    \$\\frac{dx(t)}{dt} = k_d (conc(t) - x(t))\$

**Fields**

    - `tps` -- time vector
    - `conc` -- exposure vector
    - `kd` -- parameter, scalar

Return an `Array{Array{Float64,1},1}`.

# Example
```julia
julia> myTK = runTK_MCMC([0,1,2,3], [0,1,20,2], [0.5, 2.4]);

julia> myTK.TK
2-element Array{Array{Float64,1},1}:
 [0.0, 0.21800446293645104, 4.570824232907027, 6.801199145490319]
 [0.0, 0.6019164649918195, 12.77507909562312, 8.171280475528322]
```
"""
function runTK_MCMC(tps, conc, kd)
    TK = Array{Array{Float64,1},1}(undef, length(kd))
    @inbounds for i in eachindex(kd)
        prob = DiffEqBase.ODEProblem(odeTK,[0.0],(0.0, maximum(tps)),[tps, conc, kd[i]])
        _saveat = tps === nothing ? Float64[] : tps
        sol = DiffEqBase.solve(prob ; saveat = _saveat)
        TK[i] = [sol.u[j][1] for j in 1:length(sol)]
    end
    return (TK=TK,)
end

"""
Solve ODE for Stochastic Death TKTD model

    \$\\frac{dx(t)}{dt} = k_d (conc(t) - x(t))\$

    \$\\frac{dy(t)}{dt} = k_d (x(t) - z) + h_b\$

    \$ pSurv(t) = \\exp{(-y(t))}\$

**Fields**

- `tps` -- time vector
- `conc` -- exposure vector
- `kd` -- parameters
- `hb` -- background mortality
- `z` -- threshold
- `kk` -- killing rate

Return a Tuple of `Array{Float64,1}`.

# Example
```julia
julia> mySD = runSD([0,1,2,3], [0,1,20,2], 0.5, 2.4, 16.3 ,10.0);

julia> mySD.TK
4-element Array{Float64,1}:
 0.0
 0.21089023114020433
 4.575797874578747
 6.801834882676336

julia> mySD.TD
4-element Array{Float64,1}:
 1.0
 0.09071795328941247
 0.008229747049020037
 0.0007465858083766806
```
"""
function runSD(tps, conc, kd, hb, z, kk)
    prob = DiffEqBase.ODEProblem(odeSD,[0.0, 0.0],(0.0, maximum(tps)),[tps, conc, kd, hb, z, kk])
    _saveat = tps === nothing ? Float64[] : tps
    sol = DiffEqBase.solve(prob ; saveat = _saveat)
    TK = [sol.u[j][1] for j in 1:length(sol)]
    TD = [exp(-sol.u[j][2]) for j in 1:length(sol)]
    return (TK=TK, TD=TD,)
end

"""
Solve ODE for Stochastic Death TKTD model

    \$\\frac{dx(t)}{dt} = k_d (conc(t) - x(t))\$

    \$\\frac{dy(t)}{dt} = k_d (x(t) - z) + h_b\$

    \$ pSurv(t) = \\exp{(-y(t))}\$

**Fields**

- `tps` -- time vector
- `conc` -- exposure vector
- `kd` -- parameters
- `hb` -- background mortality
- `z` -- threshold
- `kk` -- killing rate

Return a Tuple of `Array{Array{Float64,1},1}`.

# Example
```julia
julia> mySD = runSD_MCMC([0,1,2,3], [0,1,20,2], [0.5, 2.4],[2.5, 2.1],[10.5, 16.0],[5.0, 2.4]);

julia> mySD.TK
2-element Array{Array{Float64,1},1}:
 [0.0, 0.21089023114020433, 4.575797874578747, 6.801834882676336]
 [0.0, 0.6309497946830392, 12.689036317855038, 8.153727962011798]

julia> mySD.TD
2-element Array{Array{Float64,1},1}:
 [1.0, 0.08208499862389876, 0.006737946999085473, 0.0005530843701478346]
 [1.0, 0.12245642825298195, 0.014995576820477731, 0.0018363047770289104]
```

"""

function runSD_MCMC(tps, conc, kd, hb, z, kk)
    TK = Array{Array{Float64,1},1}(undef, length(kd))
    TD = Array{Array{Float64,1},1}(undef, length(kd))
    @inbounds for i in eachindex(kd)
        prob = DiffEqBase.ODEProblem(odeSD,[0.0, 0.0],(0.0, maximum(tps)),[tps, conc, kd[i], hb[i], z[i], kk[i]])
        _saveat = tps === nothing ? Float64[] : tps
        sol = DiffEqBase.solve(prob ; saveat = _saveat)
        TK[i] = [sol.u[j][1] for j in 1:length(sol)]
        TD[i] = [exp(- sol.u[j][2]) for j in 1:length(sol)]
    end
    return (TK=TK, TD=TD)
end

"""
Solve ODE for Individual Tolerance TKTD model

    \$\\frac{dx(t)}{dt} = k_d (conc(t) - x(t))\$

    \$ pSurv(t) = \\exp{(-h_b t)} (1-logLogistic(max(x(t), alpha,beta))  \$

**Fields**

- `tps` -- time vector
- `conc` -- exposure vector
- `kd` -- toxicokinetics parameter
- `hb` -- background mortality
- `alpha` -- median of log-logistic distribution
- `beta` -- shape of the log-logistic distribution

Return an `Array{Float64,1}`.

# Example
```julia
julia> myIT = runIT([0,1,2,3], [0,1,20,2], 0.5, 2.1, 10.5 , 2.4);

julia> myIT.TK
4-element Array{Float64,1}:
 0.0
 0.21800446293645104
 4.570824232907027
 6.801199145490319

julia> myIT.TD
4-element Array{Float64,1}:
 1.0
 0.12244522344410577
 0.013201809910381767
 0.001357555235565695
```

"""
function runIT(tps, conc, kd, hb, alpha, beta)
    prob = DiffEqBase.ODEProblem(odeTK,[0.0],(0.0, maximum(tps)),[tps, conc, kd])
    _saveat = tps === nothing ? Float64[] : tps
    sol = DiffEqBase.solve(prob ; saveat = _saveat)

    pSurv = Array{Float64}(undef,length(tps))
    Dmax_tmp = accumulate(max, sol.u)

    TK = [sol.u[j][1] for j in 1:length(sol)]
    Dmax_tmp = accumulate(max, sol.u)
    TD = [ exp(-hb*tps[j])*(1-logLogisticLaw(Dmax_tmp[j][1], alpha,beta)) for j in 1:length(sol)]

    return (TK=TK, TD=TD)
end

"""
Solve ODE for Individual Tolerance TKTD model

    \$\\frac{dx(t)}{dt} = k_d (conc(t) - x(t))\$

    \$ pSurv(t) = \\exp{(-h_b t)} (1-logLogistic(max(x(t), alpha,beta))  \$

**Fields**

- `tps` -- time vector
- `conc` -- exposure vector
- `kd` -- toxicokinetics parameter
- `hb` -- background mortality
- `alpha` -- median of log-logistic distribution
- `beta` -- shape of the log-logistic distribution

Return an `Array{Array{Float64,1},1}`.

# Example
```julia
julia> myIT = runIT_MCMC([0,1,2,3], [0,1,20,2], [0.5, 2.4],[2.5, 2.1],[10.5, 16.0],[5.0, 2.4]);

julia> myIT.TK
2-element Array{Array{Float64,1},1}:
 [0.0, 0.21800446293645104, 4.570824232907027, 6.801199145490319]
 [0.0, 0.6019164649918195, 12.77507909562312, 8.171280475528322]

julia> myIT.TD
2-element Array{Array{Float64,1},1}:
 [1.0, 0.0820849983072016, 0.006634237831805787, 0.0004964761672208183]
 [1.0, 0.12240978217836299, 0.009475165156808458, 0.0011602948822098684]
```

"""

function runIT_MCMC(tps, conc, kd, hb, alpha, beta)
    TK = Array{Array{Float64,1},1}(undef, length(kd))
    TD = Array{Array{Float64,1},1}(undef, length(kd))
    @inbounds for i in eachindex(kd)
        prob = DiffEqBase.ODEProblem(odeTK,[0.0],(0.0, maximum(tps)),[tps, conc, kd[i]])
        _saveat = tps === nothing ? Float64[] : tps
        sol = DiffEqBase.solve(prob ; saveat = _saveat)
        TK[i] = [sol.u[j][1] for j in 1:length(sol)]
        Dmax_tmp = accumulate(max, sol.u)
        TD[i] = [ exp(-hb[i]*tps[j])*(1-logLogisticLaw(Dmax_tmp[j][1], alpha[i],beta[i])) for j in 1:length(sol)]
    end
    return (TK=TK, TD=TD)
end

# struct ExposureProfile{T<:Number,X<:Number}
#     time::Array{T,1}
#     exposure::Array{X,1}
#     function ExposureProfile(time::Array{T,1},
#                              exposure::Array{X,1}) where {T,X}
#         if length(time) != length(exposure)
#             throw(ArgumentError("length of time and exposure differ"))
#         end
#         if length(time) != length(unique(time))
#             throw(ArgumentError("values in time vector must be unique"))
#         end
#         if sort(time) != time
#             throw(ArgumentError("time vector must be sorted"))
#         end
#         new{T,X}(time, exposure)
#     end
# end
# ExposureProfile([1.0],[2.0])

# ExposureProfile([1.0,2,3],[2.0,3.0,2.0])
# ExposureProfile([1,2,3],[2.0,3.0,2.0])
# ExposureProfile([1,2,3],[2.0,3.0,2.0]).time
# ExposureProfile([1,2,3],[2.0,3.0,2.0]).exposure
# # should return false
# ExposureProfile([1,3,3],[2.0,3.0,2.0])
# ExposureProfile([1,3,2],[2.0,3.0,2.0])
# ExposureProfile([1.0],[2.0,3.0])


text:
# Recreates the analysis from the original data files with different parameters

datapath = "dataset"
experimentdirnames = ["Mdm2-p53", "MEG3-p53"]

samplenames = ["A", "B", "C", "D"]

nreplicates = 3
nsamples = 4
ncells = 10

outputdir = "output/optimize"
mkpath(outputdir)
outputdatapath = joinpath(outputdir, "results_optimize.jld2")

using Distributed
currentworkers = addprocs(exeflags = "--project")
@everywhere using SMLMAssociationAnalysis_NCB
using Printf
using LocalizationMicroscopy
using Statistics
using FileIO

experimentresults = Vector{Vector{Vector{ResultOptimizing}}}[]
for experimentdirname ∈ experimentdirnames
    println("Starting experiment $experimentdirname.")
    experimentpath = joinpath(datapath, experimentdirname)
    replicateresults = Vector{Vector{ResultOptimizing}}[]
    for i ∈ 1:nreplicates
        sampleresults = Vector{ResultOptimizing}[]
        println("    Starting replicate $i.")
        replicatepath = joinpath(experimentpath, "Replicate $i")
        for samplename ∈ samplenames
            results = ResultOptimizing[]
            println("        Starting sample $samplename.")
            for j ∈ 1:ncells
                println("            Starting cell $j.")
                cellpath = joinpath(replicatepath, "$samplename $(Printf.@sprintf("%03i", j)).bin.txt")
                localizations = loadlocalizations(cellpath, LocalizationMicroscopy.nikonelementstext)
                # account for variances in data collection
                if experimentdirname == experimentdirnames[2] && samplename ∈ samplenames[3:4]
                    ch1_name = "561"
                else
                    ch1_name = "647"
                end
                ch2_name = "488"
                if experimentdirname == experimentdirnames[2] &&
                   ((i == 3 && samplename ∈ samplenames[1:3]) || (i == 2 && samplename == samplenames[2]))
                    ch1_startframe = 1
                    ch2_startframe = 15001
                else
                    ch1_startframe = 1
                    ch2_startframe = 11001
                end
                ch1 = getmolecules.(
                    Ref(localizations),
                    Ref(ch1_name),
                    Ref(ch1_startframe),
                    11000,
                    100,
                    10,
                    34.2,
                    [0, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000],
                    0
                )
                ch2 = getmolecules.(
                    Ref(localizations),
                    Ref(ch2_name),
                    Ref(ch2_startframe),
                    11000,
                    100,
                    10,
                    34.2,
                    [0, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000],
                    0
                )

                ch1_molecules, ch1_localizations = getmolecules(
                    localizations,
                    ch1_name,
                    ch1_startframe,
                    11000,
                    100,
                    10,
                    34.2,
                    500,
                    200,
                )
                ch2_molecules, ch2_localizations = getmolecules(
                    localizations,
                    ch2_name,
                    ch2_startframe,
                    11000,
                    100,
                    10,
                    34.2,
                    500,
                    200,
                )

                ch1_neighbors, ch2_neighbors, distances = exclusivenearestneighbors(ch1_molecules, ch2_molecules)

                percentileranks_by_localdensity = montecarloaffinity.(
                    Ref(ch1_molecules),
                    Ref(ch2_molecules),
                    Ref(ch1_neighbors),
                    Ref(ch2_neighbors),
                    Ref(distances),
                    200:200:2000,
                    10000
                )

                if length(distances) == 0
                    mediandistance = NaN
                else
                    mediandistance = median(distances)
                end
                println("                Done:$(length(distances)) neighbors from $(length(ch1_molecules)) and $(length(ch2_molecules)) molecules, $(length(ch1_localizations)) and $(length(ch2_localizations)) localizations; median distance $mediandistance")
                ch1_data = ChannelDataOptimizing(ch1_name, first.(ch1), ch1_neighbors)
                ch2_data = ChannelDataOptimizing(ch2_name, first.(ch2), ch2_neighbors)
                result = ResultOptimizing(
                    "",
                    experimentdirname,
                    i,
                    samplename,
                    j,
                    [ch1_data, ch2_data],
                    distances,
                    mediandistance,
                    percentileranks_by_localdensity
                )
                push!(results, result)
            end
            push!(sampleresults, results)
        end
        push!(replicateresults, sampleresults)
    end
    push!(experimentresults, replicateresults)
end

save(outputdatapath, "experimentresults", experimentresults)

rmprocs(currentworkers)


text:
c = one(PermCycs)
p = one(PermList)
s = one(PermSparse)
m = one(PermMat)

@test null(c) == null(p) == null(s)

for el in ( :c, :p, :s, :m)
    @eval begin
        @test isid($el)
        @test issym($el)
        @test isposdef($el)
        @test iseven($el)
    end
end

##

p = PermList([3,1,2])
c = cycles(p)
s = psparse(p)
m = matrix(p)

for el in ( :c, :p, :s, :m)
    @eval begin
        @test ! isid($el)
        @test ! issym($el)
        @test ! isposdef($el)
        @test inv($el) == transpose($el)
        @test rank($el) == plength($el)
        @test det($el) == sign($el)
        @test iseven($el)
    end
end

p = PermList([2,7,8,3,1,9,6,4,5,10])
@test pivtopermlist(topiv(p)) == p

q = PermList([5,6,9,1,4,7,8,3,2,10])
q1 = pivtopermlist([5,6,9,5,5,7,8,9,9,10])
@test q1 == q


text:
# ---
# title: 307. Range Sum Query - Mutable
# id: problem307
# author: zhwang
# date: 2022-03-12
# difficulty: Medium
# categories: Binary Indexed Tree, Segment Tree
# link: <https://leetcode.com/problems/range-sum-query-mutable/description/>
# hidden: true
# ---
#
# Given an integer array _nums_ , find the sum of the elements between indices
# _i_ and _j_ ( _i_ ≤ _j_ ), inclusive.
#
# The _update(i, val)_ function modifies _nums_ by updating the element at index
# _i_ to _val_.
#
# **Example:**
#
#
#
#     Given nums = [1, 3, 5]
#
#     sumRange(0, 2) -> 9
#     update(1, 2)
#     sumRange(0, 2) -> 8
#
#
#
#
# **Constraints:**
#
#   * The array is only modifiable by the _update_ function.
#   * You may assume the number of calls to _update_ and _sumRange_ function is distributed evenly.
#   * `0 <= i <= j <= nums.length - 1`
#
#
## @lc code=start
using LeetCode

mutable struct SegmentTree{T<:Real}
    n::Int
    tree::Vector{T}

    function SegmentTree(nums::Vector{T}) where T
        n = length(nums)
        tree = append!(zeros(Int, n - 1), nums)
        for i in (n - 1):-1:1
            tree[i] = tree[i << 1] + tree[i << 1 | 1]
        end
        new{T}(n, tree)
    end
end

function update!(ST::SegmentTree, ind::Int, val::Int)::Nothing
    tree = ST.tree
    ind += ST.n - 1
    delta = val - tree[ind]
    while ind > 0
        tree[ind] += delta
        ind >>= 1
    end
end

function sum_range(ST::SegmentTree, left::Int, right::Int)::Int
    left += ST.n - 1
    right += ST.n - 1
    res, tree = 0, ST.tree
    while left <= right
        if isodd(left) ## right child
            res += tree[left] ## record outside value
            left += 1
        end
        left >>= 1
        if iseven(right) ## left child
            res += tree[right] ## record outside value
            right -= 1
        end
        right >>= 1
    end
    res
end
## @lc code=end

## @lc test=start
@testset "307.range-sum-query-mutable.jl" begin
    @testset "307 case 1" begin
        ST = SegmentTree([1, 3, 5])
        @test sum_range(ST, 1, 3) == 9
        update!(ST, 2, 2)
        @test sum_range(ST, 1, 3) == 8
    end
end
## @lc test=end

## @lc info=start
# link: [solution to 307]()
## @lc info=end


text:
# Author: Lukasz Bratos

# Wektory do testowania, odpowiednio dla Float64 i Float32 oraz prawidlowa odpowiedz
x64 = [2.718281828,-3.141592654,1.414213562,0.5772156649,0.3010299957]
y64 = [1486.2497,878366.9879,-22.37492,4773714.647,0.000185049]
x32 = map(a -> convert(Float32, a), x64)
y32 = map(a -> convert(Float32, a), y64)
answer = -1.00657107000000e-11

# Algorytm 'w przod'. Przyjmuje wektor x oraz y i rozmiar wektorow n.
function forward(x, y, n)
    s = 0.0
    for i in 1:n
        s = s + x[i] * y[i]
    end
    return s
end

# Algorytm 'w tyl'. Przyjmuje wektor x oraz y i rozmiar wektorow n.
function backward(x, y, n)
    s = 0.0
    for i in reverse(1:n)
        s = s + x[i] * y[i]
    end
    return s
end

# Algorytm 'od najwiekszego do najmniejszego'. Przyjmuje wektor x oraz y, typ zmiennopozycyjny w ktorym wykonywane sa obliczenia oraz rozmiar wektorow n.
function fromTheBiggestToTheSmallest(x, y, type, n)
    partials = Array{type}(undef, 5)
    for i in 1:n
        partials[i] = x[i] * y[i]
    end

    pos = sort(filter(x -> x > 0, partials), rev=true)
    neg = sort(filter(x -> x < 0, partials))

    return foldl(+, pos) + foldl(+, neg)
end

# Algorytm 'od najmniejszego do najwiekszego'. Przyjmuje wektor x oraz y, typ zmiennopozycyjny w ktorym wykonywane sa obliczenia oraz rozmiar wektorow n.
function fromTheSmallestToTheBiggest(x, y, type, n)
    partials = Array{type}(undef, 5)
    for i in 1:n
        partials[i] = x[i] * y[i]
    end

    pos = sort(filter(x -> x > 0, partials))
    neg = sort(filter(x -> x < 0, partials), rev=true)

    return foldl(+, pos) + foldl(+, neg)
end

println("Wynik dla 'w przod' i Float64: ", forward(x64, y64, 5))
println("Wynik dla 'w przod' i Float32: ", forward(x32, y32, 5))

println("Wynik dla 'w tyl' i Float64: ", backward(x64, y64, 5))
println("Wynik dla 'w tyl' i Float32: ", backward(x32, y32, 5))

println("Wynik dla najwiekszy-najmniejszy i Float64: ", fromTheBiggestToTheSmallest(x64, y64, Float64, 5))
println("Wynik dla najwiekszy-najmniejszy i Float32: ", fromTheBiggestToTheSmallest(x32, y32, Float32, 5))

println("Wynik dla najmniejszy-najwiekszy i Float64: ", fromTheSmallestToTheBiggest(x64, y64, Float64, 5))
println("Wynik dla najmniejszy-najwiekszy i Float32: ", fromTheSmallestToTheBiggest(x32, y32, Float32, 5))


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Order

import ..@__MODULE__, ..parentmodule
const Base = parentmodule(@__MODULE__)
import .Base:
    AbstractVector, @propagate_inbounds, isless, identity, getindex,
    +, -, !, &, <, |

## notions of element ordering ##

export # not exported by Base
    Ordering, Forward, Reverse, Backward,
    By, Less, Perm,
    ReverseOrdering, ForwardOrdering,
    DirectOrdering

abstract type Ordering end

struct Less{T<:Function} <: Ordering
    is_less::T
end

struct By{T,O<:Ordering} <: Ordering
    by::T
    is_less::O
end

struct Reverse{O<:Ordering} <: Ordering
    is_less::O
end

struct Perm{O<:Ordering,V<:AbstractVector} <: Ordering
    order::O
    data::V
end

# Simplify some things
Reverse(ord::Reverse) = ord.is_less
By(by::typeof(identity), ord::Ordering) = ord

(o::Less)(a, b)    = o.is_less(a, b)
(o::By)(a, b)      = o.is_less(o.by(a), o.by(b))
(o::Reverse)(a, b) = o.is_less(b, a)

const Forward = Less(isless)
const ForwardOrdering = typeof(Forward)
const Backward = Reverse(Forward)
const BackwardOrdering = typeof(Backward)
const DirectOrdering = Union{ForwardOrdering,BackwardOrdering}

_ord(lt, by, order::ForwardOrdering) = By(by,Less(lt))
_ord(lt, by, order::BackwardOrdering) = Reverse(By(by,Less(lt)))

ord(lt, by, rev::Nothing, order::Ordering=Forward) = _ord(lt, by, order)

function ord(lt, by, rev::Bool, order::Ordering=Forward)
    o = _ord(lt, by, order)
    return rev ? Reverse(o) : o
end

end

max_stars_repo_name: NymanLauri/julia

max_issues_repo_name: NymanLauri/julia

max_forks_repo_name: NymanLauri/julia


text:
export computeViscosity_TauII,
    computeViscosity_EpsII, computeViscosity_TauII!, computeViscosity_EpsII!

"""
    compute viscosity given strain rate 2nd invariant

    τ = 2ηε -> η = τ/2/ε
"""
function computeViscosity_EpsII(εII, v, args)
    τII = computeCreepLaw_TauII(εII, v, args) # gives
    η = 0.5 * τII / εII
    return η
end

"""
    compute viscosity given strain rate 2nd invariant

    τ = 2ηε -> η = τ/2/ε
"""
@inline function computeViscosity_TauII(τII, v, args)
    εII = computeCreepLaw_EpsII(τII, v, args)
    η = 0.5 * τII / εII
    return η
end

@inline function computeViscosity_TauII(τII::T, v::Tuple, args) where {T}
    return computeViscosity(computeViscosity_TauII, τII, v, args, Val(length(v)))
end
@inline function computeViscosity_EpsII(εII::T, v::Tuple, args) where {T}
    return computeViscosity(computeViscosity_EpsII, εII, v, args, Val(length(v)))
end

@inline @generated function computeViscosity(
    fn::F, CII::T, v::Tuple, args::NamedTuple, ::Val{N}
) where {F,T,N}
    quote
        η = 0.0
        Base.Cartesian.@nexprs $N i -> η += 1 / fn(CII, v[i], args)
        return 1 / η
    end
end

@inline function computeViscosity_TauII!(η::AbstractArray{T,nDim}, τII::AbstractArray{T,nDim}, v::Tuple, args; cutoff=(1e16, 1e25)) where {T, nDim}
    Threads.@threads for I in eachindex(τII)
        η[I] = max(
            cutoff[1], min(cutoff[2], computeViscosity(
                    computeViscosity_TauII,
                    τII[I],
                    v,
                    (; zip(keys(args), getindex.(values(args), I))...),
                    Val(length(v)),
                )
            )
        )
    end
end

@inline function computeViscosity_EpsII!(η::AbstractArray{T,nDim}, τII::AbstractArray{T,nDim}, v::Tuple, args; cutoff=(1e16, 1e25)) where {T, nDim}
    Threads.@threads for I in eachindex(τII)
        η[I] =  max(
            cutoff[1], min(cutoff[2], computeViscosity(
                computeViscosity_EpsII,
                    τII[I],
                    v,
                    (; zip(keys(args), getindex.(values(args), I))...),
                    Val(length(v)),
                )
            )
        )
    end
end

# OPTION 1
strainCircuit(TauII, v, args) = strainCircuit(TauII, v, args, Val(length(v)))

@inline @generated function strainCircuit(TauII, v, args, ::Val{N}; n=1) where {N}
    quote
        c = 0.0
        Base.Cartesian.@nexprs $N i ->
            c += if v[i] isa Tuple
                1 / strainCircuit(TauII, v[i], args, Val(length(v[i])); n=-1)
            else
                computeCreepLaw_EpsII(TauII, v[i], args)^n
            end
        return c
    end
end

function viscosityCircuit_TauII(τII::T, v, args) where {T}
    εII = strainCircuit(τII, v, args)
    return η = 2 * εII / τII
end

# v = (DiffusionCreep(), DislocationCreep())
# v = (DiffusionCreep(), (DiffusionCreep(), DislocationCreep()))

# EpsII = 1e-10
# TauII = 20e6
# P, T, f, d = 20e6, 300.0, 1.0, 1.0
# args = (P=P, T=T, f=f, d=d)

# @benchmark viscosityCircuit_TauII($TauII, $v, $args)
# @benchmark strainCircuit($TauII, $v, $args)

# strainCircuit(TauII, v, args)

# computeCreepLaw_TauII(EpsII, DiffusionCreep(), args)
# @btime computeCreepLaw_TauII($EpsII, $(DiffusionCreep()), $args)

# computeCreepLaw_TauII(EpsII, DislocationCreep(), args)
# computeCreepLaw_EpsII(TauII, DiffusionCreep(), args)
# computeCreepLaw_EpsII(TauII, DislocationCreep(), args)

# v = (DiffusionCreep(), DislocationCreep())

# computeViscosity_EpsII2(EpsII, v, args)
# computeViscosity_EpsII(EpsII, v, args)

# @benchmark computeViscosity_EpsII2($EpsII, $v, $args)
# @benchmark computeViscosity_EpsII($EpsII, $v, $args)


text:
module CourseMatch

	VERSION = VersionNumber(0,0,1)

using JuMP, Gurobi, CSV, DataFrames

# includes
include("demand.jl")
include("student.jl")

# exports

end # module

max_stars_repo_name: ScPo-CompEcon/CourseMatch

max_issues_repo_name: ScPo-CompEcon/CourseMatch

max_forks_repo_name: ScPo-CompEcon/CourseMatch


text:
# This is the module file
module ThreeBodyProblem

using OrdinaryDiffEq
using LinearAlgebra
using PolynomialRoots
using RecipesBase
using SPICE

include("constants.jl")
export JD, JY, AU, G

include("parameters.jl")
export  System,
        Body,
        set_system,
        sun_mercury,
        sun_venus,
        sun_earth,
        earth_moon,
        sun_mars,
        sun_jupiter,
        sun_saturn,
        sun_uranus,
        sun_neptune,
        jupiter_europa,
        saturn_enceladus,
        earth_moon_sun,
        SUN,
        MERCURY,
        VENUS,
        EARTH,
        MARS,
        JUPITER,
        SATURN,
        URANUS,
        NEPTUNE,
        PLUTO,
        MOON,
        PHOBOS,
        DEIMOS,
        IO,
        EUROPA,
        GANYMEDE,
        CALLISTO,
        AMALTHEA,
        MIMAS,
        ENCELADUS,
        TETHYS,
        DIONE,
        RHEA,
        TITAN,
        HYPERION,
        IAPETUS,
        PHOEBE,
        ARIEL,
        UMBRIEL,
        TITANIA,
        OBERON,
        MIRANDA,
        TRITON,
        NEREID,
        PROTEUS,
        CHARON,
        NIX,
        HYDRA,
        KERBEROS,
        STYX

include("frames.jl")
export  rot2inert!,
        rot2inert,
        inert2rot!,
        inert2rot,
        r2latlon,
        enu2ecef,
        ecef2enu,
        ecef2eci,
        eci2ecef,
        eci2sci,
        sci2eci,
        dimensionalize!,
        dimensionalize,
        nondimensionalize!,
        nondimensionalize

include("dynamics.jl")
export  R2BPdynamics!,
        R2BPdynamics,
        CR3BPdynamics!,
        CR3BPdynamics,
        CR3BPstm!,
        CR3BPstm,
        CR3BPinert!,
        CR3BPinert,
        CWdynamics!,
        CWdynamics,
        BCPdynamics!,
        BCPdynamics,
        BCPstm!,
        BCPstm

include("orbitalelements.jl")
export  AP2a,
        cart2oe,
        oe2cart,
        nu2E,
        E2nu,
        E2M,
        M2E,
        azel2cart,
        cart2azel

include("util.jl")
export  computeR1R2,
        computer1r2,
        computeL1,
        computeL2,
        computeL3,
        computeL4,
        computeL5,
        computeLpts,
        computeUeff,
        computeC,
        rotx,
        rotxd,
        roty,
        rotyd,
        rotz,
        rotzd,
        rotlatlon,
        date2mjd,
        mjd2gmst,
        wrapto360,
        wrapto180,
        wrapto2pi,
        wraptopi,
        date2str,
        deserno_sphere,
        deserno_hemisphere,
        spherical_ring,
        gammaL

include("plot.jl")
export  circle,
        sphere
        # seczoom,
        # plot_prim!,
        # plot_prim,
        # plot_sec!,
        # plot_sec,
        # plot_Lpts!,
        # plot_Lpts

include("special.jl")
export  invariant_manifolds,
        rich3,
        differential_corrector

end # module


text:
# Plot recipes

# Recipe for Connectivity Histogram (Cardinality Histogram)
@recipe function f(W::SpatialWeights)
    xguide --> "Number of Neighbors"
    legend --> false
    seriestype --> :histogram
    cardinalities(W)
end

# Recipe for Moran plot
@recipe function f(x::Vector{Float64}, W::SpatialWeights)

    Wx = slag(W, x)

    zstandardize  = get(plotattributes, :standardize, true)
    if zstandardize
        z = standardize(ZScoreTransform, x)
        Wz = standardize(ZScoreTransform, Wx)
    else
        z = x
        Wz = Wx
    end

    title --> "Moran Scatterplot"
    xguide --> "Attribute"
    yguide --> "Spatial Lag of " * plotattributes[:xguide]
    legend --> false
    grid --> false

    # Variable and spatial lag of variable
    @series begin
        seriestype := :scatter
        smooth := true
        z, Wz
    end

    # Vertical line at variable mean
    zmean = mean(z)
    @series begin
        seriestype := :vline
        linestyle := :dash
        seriescolor --> :red
        [zmean]
    end

    # Horizontal line at W * variable mean
    Wzmean = mean(Wz)
    @series begin
        seriestype := :hline
        linestyle := :dash
        seriescolor --> :red
        [Wzmean]
    end

end

# Recipe for Global Spatial Autocorrelation reference distribution
@recipe function plot(x::AbstractGlobalSpatialAutocorrelation)

    S = score(x)
    Sperms = scoreperms(x)

    legend --> false

    # Histogram of permutated values
    @series begin
        seriestype := :histogram
        Sperms
    end

    # Vertical line at Moran's i
    @series begin
        seriestype := :vline
        seriescolor --> :red
        [S]
    end

end


text:
module ParametersTests

using Base.Test
using BenchmarkTools
using BenchmarkTools: Parameters

@test BenchmarkTools.DEFAULT_PARAMETERS == Parameters()

p = Parameters(seconds = 1, gctrial = false)
oldseconds = BenchmarkTools.DEFAULT_PARAMETERS.seconds
oldgctrial = BenchmarkTools.DEFAULT_PARAMETERS.gctrial
BenchmarkTools.DEFAULT_PARAMETERS.seconds = p.seconds
BenchmarkTools.DEFAULT_PARAMETERS.gctrial = p.gctrial
@test p == Parameters()
@test Parameters(p; evals = 3, time_tolerance = .32) == Parameters(evals = 3, time_tolerance = .32)
BenchmarkTools.DEFAULT_PARAMETERS.seconds = oldseconds
BenchmarkTools.DEFAULT_PARAMETERS.gctrial = oldgctrial

p = Parameters(seconds = 1, gctrial = false, samples = 2, evals = 2, overhead = 42,
               gcsample = false, time_tolerance = 0.043, memory_tolerance = 0.15)
oldseconds = BenchmarkTools.DEFAULT_PARAMETERS.seconds
oldgctrial = BenchmarkTools.DEFAULT_PARAMETERS.gctrial
old_time_tolerance = BenchmarkTools.DEFAULT_PARAMETERS.time_tolerance
old_memory_tolerance = BenchmarkTools.DEFAULT_PARAMETERS.memory_tolerance
oldsamples = BenchmarkTools.DEFAULT_PARAMETERS.samples
oldevals = BenchmarkTools.DEFAULT_PARAMETERS.evals
oldoverhead = BenchmarkTools.DEFAULT_PARAMETERS.overhead
oldgcsample = BenchmarkTools.DEFAULT_PARAMETERS.gcsample
BenchmarkTools.DEFAULT_PARAMETERS.seconds = p.seconds
BenchmarkTools.DEFAULT_PARAMETERS.gctrial = p.gctrial
BenchmarkTools.DEFAULT_PARAMETERS.time_tolerance = p.time_tolerance
BenchmarkTools.DEFAULT_PARAMETERS.memory_tolerance = p.memory_tolerance
BenchmarkTools.DEFAULT_PARAMETERS.samples = p.samples
BenchmarkTools.DEFAULT_PARAMETERS.evals = p.evals
BenchmarkTools.DEFAULT_PARAMETERS.overhead = p.overhead
BenchmarkTools.DEFAULT_PARAMETERS.gcsample = p.gcsample
@test p == Parameters()
@test p == Parameters(p)
BenchmarkTools.DEFAULT_PARAMETERS.seconds = oldseconds
BenchmarkTools.DEFAULT_PARAMETERS.gctrial = oldgctrial
BenchmarkTools.DEFAULT_PARAMETERS.time_tolerance = old_time_tolerance
BenchmarkTools.DEFAULT_PARAMETERS.memory_tolerance = old_memory_tolerance
BenchmarkTools.DEFAULT_PARAMETERS.samples = oldsamples
BenchmarkTools.DEFAULT_PARAMETERS.evals = oldevals
BenchmarkTools.DEFAULT_PARAMETERS.overhead = oldoverhead
BenchmarkTools.DEFAULT_PARAMETERS.gcsample = oldgcsample

end # module


text:
mutable struct PQPolygonSDKUserModel

    # data about the user -
    email::String
    apikey::String

    PQPolygonSDKUserModel() = new()
end

# what is the base type for all endpoint calls?
abstract type AbstractPolygonEndpointModel end

mutable struct PolygonGroupedDailyEndpointModel <: AbstractPolygonEndpointModel

    # data -
    date::Date
    adjusted::Bool
    apikey::String

    # constructor -
    PolygonGroupedDailyEndpointModel() = new()
end

mutable struct PolygonDailyOpenCloseEndpointModel <: AbstractPolygonEndpointModel

    # data -
    ticker::String
    adjusted::Bool
    date::Date
    apikey::String

    # constructor -
    PolygonDailyOpenCloseEndpointModel() = new()
end

mutable struct PolygonPreviousCloseEndpointModel <: AbstractPolygonEndpointModel

    # data -
    ticker::String
    adjusted::Bool
    apikey::String

    # constructor -
    PolygonPreviousCloseEndpointModel() = new()
end

mutable struct PolygonAggregatesEndpointModel <: AbstractPolygonEndpointModel

    # data -
    adjusted::Bool
    limit::Int64
    sortdirection::String
    ticker::String
    to::Date
    from::Date
    multiplier::Int64
    timespan::String
    apikey::String

    # constructor -
    PolygonAggregatesEndpointModel() = new()
end

mutable struct PolygonTickerNewsEndpointModel <: AbstractPolygonEndpointModel

    # data -
    ticker::String
    published_utc::Union{Date, Nothing}
    order::String
    limit::Int64
    sort::Union{String, Nothing}
    apikey::String

    # constructor -
    PolygonTickerNewsEndpointModel() = new()
end

mutable struct PolygonTickerDetailsEndpointModel <: AbstractPolygonEndpointModel #ycpan

    # data -
    apikey::String
    ticker::String
    date::String

    # constructor -
    PolygonTickerDetailsEndpointModel() = new()
end

mutable struct PolygonOptionsContractReferenceEndpoint <: AbstractPolygonEndpointModel

    # data -
    ticker::Union{Nothing,String}
    underlying_ticker::Union{Nothing, String}
    contract_type::String
    expiration_date::Date
    limit::Int64
    order::String
    sort::String
    apikey::String

    # constructor -
    PolygonOptionsContractReferenceEndpoint() = new()
end

mutable struct PolygonMarketHolidaysEndpointModel <: AbstractPolygonEndpointModel #ycpan

    # data -
    apikey::String

    # constructor -
    PolygonMarketHolidaysEndpointModel() = new()
end

mutable struct PolygonExchangesEndpointModel <: AbstractPolygonEndpointModel #ycpan

    # data -
    apikey::String
    asset_class::String
    locale::String

    # constructor -
    PolygonExchangesEndpointModel() = new()
end

mutable struct PolygonStockSplitsEndpointModel <: AbstractPolygonEndpointModel #ycpan

    # data -
    apikey::String
    ticker::String
    execution_date::String
    reverse_split::String
    order::String
    limit::String
    sort::String

    # constructor -
    PolygonStockSplitsEndpointModel() = new()
end

mutable struct PolygonMarketStatusEndpointModel <: AbstractPolygonEndpointModel #ycpan

    # data -
    apikey::String

    # constructor -
    PolygonMarketStatusEndpointModel() = new()
end

mutable struct PolygonDividendsEndpointModel <: AbstractPolygonEndpointModel #ycpan

    # data -
    apikey::String
    ticker::String
    ex_dividend_date::String
    record_date::String
    declaration_date::String
    pay_date::String
    frequency::String
    cash_amount::String
    dividend_type::String
    order::String
    limit::String
    sort::String

    # constructor -
    PolygonDividendsEndpointModel() = new()
end

mutable struct PolygonTickersEndpointModel <: AbstractPolygonEndpointModel #ycpan

    # data -
    apikey::String
    ticker::String
    type::String
    market::String
    exchange::String
    cusip::String
    cik::String
    date::String
    search::String
    active::Bool
    sort::String
    order::String
    limit::String

    # constructor -
    PolygonTickersEndpointModel() = new()
end

mutable struct PolygonConditionsEndpointModel <: AbstractPolygonEndpointModel #ycpan

    # data -
    apikey::String
    asset_class::String
    data_type::String
    id::String
    sip::String
    order::String
    limit::String
    sort::String

    # constructor -
    PolygonConditionsEndpointModel() = new()
end

mutable struct PolygonStockFinancialsEndpointModel <: AbstractPolygonEndpointModel #ycpan

    # data -
    apikey::String
    ticker::String
    cik::String
    company_name::String
    sic::String
    filing_date::String
    period_of_report_date::String
    timeframe::String
    include_sources::String
    order::String
    limit::String
    sort::String

    # constructor -
    PolygonStockFinancialsEndpointModel() = new()
end

mutable struct PolygonTickerTypesEndpointModel <: AbstractPolygonEndpointModel #ycpan

    # data -
    apikey::String
    asset_class::String
    locale::String

    # constructor -
    PolygonTickerTypesEndpointModel() = new()
end


text:
function _applicable()

end

# Returns whether it is applicable and the matched wilds
# a truth can match with an expr in many ways
function _isapplicable(expr::Node, truth::Node, wilds::Array{Node})::Tuple{Bool,Dict{Node,Set{Node}}}
    result = true
    matched = Dict(wild => Set{Node}() for wild in wilds)  # wehre the wild needs to take on this value and this values

    # println(class(truth))
    # println(quick_range(truth))
    # println(expr)
    # println(quick_range(expr))
    # println()

    # if (class(truth) <: ACall) && (class(args(truth)[1]) <: Wild) && (quick_range(expr) <: quick_range(args(truth)[1]))
    #     inner_result, inner_matched = _isapplicable(expr, )
    #     push!(matched[args(truth)[1]], expr)
    # else
    if (class(truth) <: Wild) && (quick_range(expr) <: quick_range(truth))
        result = true
        push!(matched[truth], expr)
    elseif expr.class != truth.class
        result = false
    elseif expr.data != truth.data
        result = false
    elseif length(expr.args) != length(truth.args)
        result = false
    else
        for i in 1:length(expr.args)
            inner_result, inner_matched = _isapplicable(expr.args[i], truth.args[i], wilds)
            for wild in wilds
                matched[wild] = union(matched[wild], inner_matched[wild])
            end
            if !inner_result
                result = false
            end
        end
    end

    for wild in wilds
        if length(matched[wild]) > 1
            result = false
        end
    end

    return result, matched
end

preorder_traversal(x::Node) = vcat([x], [preorder_traversal(arg) for arg in args(x)]...)
function count_ops(x::Node)
    if class(x) <: AAtom
        return 0
    else
        return length(args(x)) + sum([count_ops(arg) for arg in args(x)])
    end
end

function replace(expr::Node, old::Node, new::Node)::Node
    if expr == old
        return new
    else
        new_args = Tuple(replace(arg, old, new) for arg in args(expr))
        return Node(data(expr), class(expr), new_args)
    end
end

function replace(expr::Node, replacements::Dict{Node,Node})::Node
    result = expr
    for (key, value) in replacements
        result = replace(result, key, value)
    end
    return result
end

function get_wilds(expr::Node)::Array{Node}
    result = Node[]

    nodes = preorder_traversal(expr)
    for node in nodes
        if (class(node) <: ACall) && (class(args(node)[1]) <: Wild)
            push!(result, node)
        elseif class(node) <: AWild
            push!(result, node)
        end
    end

    return unique(result)
end

function isapplicable(expr::Node, truth::Node)
    wilds = get_wilds(truth)
    result, matched = _isapplicable(expr, truth, wilds)

    new_matched = Dict{Node,Node}()
    if result
        for wild in wilds
            new_matched[wild] = pop!(matched[wild])
        end
    end

    return result, new_matched
end

function rewrites(expr::Node, truth_table=truths)::Array{Node}
    result = Node[]
    for truth in truth_table
        app, matches = isapplicable(expr, truth.left)
        if app
            subs_right = replace(truth.right, matches)
            push!(result, subs_right)
        end
    end

    # Repeat recursivley on the arguments
    for (i, child) in enumerate(expr.args)

        child_rewrites = rewrites(child, truth_table)
        for j in 1:length(child_rewrites)
            # args = copy(expr.args)
            args = collect(expr.args)
            args[i] = child_rewrites[j]
            child_rewrites[j] = Node(expr.data, expr.class, Tuple(args))
        end
        append!(result, child_rewrites)
    end

    return unique(result)
end

function result(x, y=TRUE, maxdepth=2, truth_table=truths)  # if x ⟹ y it will return true, false or missing
    depth = 0
    values = Node[x]
    old_values = values  # checked trees
    while !(y in values) && !(~(y) in values) && (depth < maxdepth)
        depth += 1
        new_values = []
        for value in values
            append!(new_values, rewrites(value, truth_table))
        end
        values = setdiff(new_values, old_values)
        append!(old_values, values)
        # values = values[map(count_ops, values) .<= maxops]
        println(length(values))
    end

    if y in values
        return TRUE
    elseif ~(y) in values
        return FALSE
    else
        return UNKNOWN
    end
end

function rewrites_wp(expr::Node, node_id::Array{Int64}=Int64[], truth_table=truths)::Array{Tuple{Node,Int64,Array{Int64}}}
    result::Array{Tuple{Node,Int64,Array{Int64, M} where M},N} where N = []
    for (truth_id, truth) in enumerate(truth_table)
        app, matches = isapplicable(expr, truth.left)
        if app
            subs_right = replace(truth.right, matches)
            push!(result, (subs_right, truth_id, node_id))
        end
    end

    # Repeat recursivley on the arguments
    for (i, child) in enumerate(expr.args)

        child_rewrites = rewrites_wp(child, vcat(node_id, [i]), truth_table)
        for j in 1:length(child_rewrites)
            args = collect(expr.args)
            args[i] = child_rewrites[j][1]
            child_rewrites[j] = (Node(expr.data, expr.class, Tuple(args)), child_rewrites[j][2], child_rewrites[j][3])
        end
        append!(result, child_rewrites)
    end

    return unique(result)
end

function prove(x, y=TRUE, maxdepth=2, truth_table=truths)
    depth = 0
    values = Proof[Proof([(x, 1, [])])]
    old_values = values  # checked trees
    last_line = [steps(p)[end] for p in values]
    while !(y in last_line) && !(~(y) in last_line) && (depth < maxdepth)
        depth += 1
        new_values = []
        for value in values
            append!(new_values, extensions(value, truth_table))
        end
        # values = [Set(new_values)...]
        values = setdiff(new_values, old_values)
        append!(old_values, values)
        last_line = [steps(p)[end] for p in values]
        println(length(values))
        for value in last_line
            if length(get_wilds(value)) > 0
                println(value)
            end
        end
    end
    last_line = [steps(p)[end] for p in values]
    if y in last_line
        proof = values[[steps(p)[end]==y for p in values]][1]
        return proof
    elseif ~(y) in last_line
        proof = values[[steps(p)[end]==~(y) for p in values]][1]
        return proof
    else
        return UNKNOWN
    end
end

function result_wp(x, y=TRUE, maxdepth=2)
    return english(prove(x, y, maxdepth))
end


text:
## TODO: pending https://github.com/JuliaLang/julia/issues/29420
# this one is suggested in the issue, but it looks like time_t and tm are two different things?
# const Ctime_t = Base.Libc.TmStruct

const Ctm = Base.Libc.TmStruct
const Ctime_t = UInt
const Cclock_t = UInt

const int8_t = Int8
const uint8_t = UInt8
const int16_t = Int16
const uint16_t = UInt16
const int32_t = Int32
const uint32_t = UInt32
const uintptr_t = Csize_t


text:
using Devectorize
using Base.Test

#################################################
#
#   devec_transform
#
#################################################

n = 4
d = {:x => rand(n),
     :y => pi * [1:n],
     :z => 1 / [1:n]}

@devec_transform d  xd = x .* y + z

@test isequal(d[:xd], d[:x] .* d[:y] + d[:z])

@devec_transform(d,
                 xd => x .* y + z,
                 yd => x .* x,
                 zd => blend(y .> 7, z, -z))

@test isequal(d[:yd], d[:x] .* d[:x])


text:
using SyntaxTree, Reduce
using Test
@force using Reduce.Algebra

# write your own tests here
@test showerror(stdout,ReduceError("A Portable General-Purpose Computer Algebra System")) == nothing
@test rcall(:((1+pi)^2)) == convert(Expr,RExpr(rcall("(1+pi)**2")))
@test try; "1/0" |> rcall; false; catch; true; end
@test Reduce.Reset() == nothing
@test display(RExpr("(x+i)^3")) == nothing; print('\n')
@test Reduce._syme(Reduce.r_to_jl) |> typeof == String
@test R"x+2" == R"2+x-1+1"
@test :((x+1+π)^2; int(1/(1+x^3),x)) |> RExpr |> Reduce.parse |> typeof == Expr
@test !Base.process_exited(Reduce.rs)
@test string(R"x+1") |> typeof == String
@test RExpr(:x) == R"x"
@test RExpr(:x)*R"x" == R"x^2"
@test convert(RExpr,R"x").str == convert(Array{String,1},R"x")
@test load_package([:rlfi]) == load_package(:rlfi,:rlfi)
@test show(stdout, R"") == nothing
@test show(stdout,"text/latex",R"int(sinh(e**i*z),z)") == nothing
@test Base.write(Reduce.rs, R"") |> typeof == Int; rcall(1)
@test (x = :(x^2+2x+1); rcall(x,off=[:factor]) == x)
@test rcall("x + 1","factor") == "x + 1"
@test Expr(:function,:fun,:(return begin; x = 700; y = x; end)) |> RExpr |> Reduce.parse |> typeof == Expr
@test Expr(:for,:(i=2:34),:(product(i))) |> rcall |> eval |> typeof == BigInt
@test try; Expr(:type,false,:x) |> RExpr; false; catch; true; end
@test try; :(@time f(x)) |> RExpr; false; catch; true; end
@test (x = Expr(:function,:fun,:(return y=a^3+3*a^2*b+3*a*b^2+b^3)); x==x |> Reduce.factor |> expand)
@test try; Expr(:for,:(i=2:34),:(product(i))) |> RExpr |> Reduce.parse; false; catch; true; end
@test R"begin; 1:2; end" |> Reduce.parse |> RExpr |> string == "1:2 "
@test latex(:(x+1)) |> typeof == String
@test length(:(x+y)) |> typeof == Int
@test log(:(ℯ^x)) == :x
@test nextprime(100) == 101
@test ceiling(1.2) == 2
@test impart(:(1+2*im)) == 2
@test impart(2+1.7im) == 17/10
@test bernoulli(2) == 1/6
@test Reduce.parsegen(:parsetest,:expr) |> typeof == Expr
@test Reduce.parsegen(:calctest,:args) |> typeof == Expr
@test Reduce.parsegen(:switchtest,:switch) |> typeof == Expr
@test Reduce.parsegen(:unarytest,:unary) |> typeof == Expr
Sys.islinux() && @test Reduce.RSymReplace("!#03a9; *x**2 + !#03a9;") |> typeof == String
@test :((x+im+π)^2; int(1/(1+x^3),x)) |> RExpr |> rcall |> Reduce.parse |> typeof == Expr
@test :(int(sin(im*x+pi)^2-1,x)) |> rcall |> typeof == Expr
@test int(:(x^2+y),:x) |> RExpr == int("x^2+y","x") |> RExpr
@test R"/(2,begin 2; +(7,4); return +(4,*(2,7))+9 end)" |> Reduce.parse |> typeof == Expr
@test df(Expr(:function,:fun,:(return begin; zn = z^2+c; nz = z^3-1; end))|>RExpr,:z) |> typeof == RExpr
@test :([1 2; 3 4]) |> RExpr |> Reduce.parse |> RExpr == [1 2; 3 4] |> RExpr
println()
#@test Reduce.repl_init(Base.active_repl)==nothing

@test nextprime("3") == "5"
@test expand("(x-2)^2") |> RExpr == R"(x-2)^2"
@test nat("x+1") |> RExpr == "\nx + 1\n" |> RExpr
@test macroexpand(@__MODULE__(),@factor(:(x^2+2x+1))) == :((x+1)^2)
@test :x^2 == :(x^2)
@test NaN//NaN |> isnan
@test join(split(R"x+1;x+2"))|> string == "x+1;\nx+2"
@test sub(:x=>7,:x+7) == sub([:x=>7,:z=>21],:z-:x)
@test SyntaxTree.sub(Float64,prod((:x-:n)^:n,:n,1,7)|>horner) |> typeof == Expr
@test squash(Expr(:function,:(fun(x)),:(z=3;z+=:x))).args[2] == squash(:(y=:x;y+=3))
@test squash(:(sqrt(x)^2)) == :x
@test Expr(:block,:(x+1)) |> RExpr == R"1+x"
@test limit((1-1/:n)^-:n,:n,Inf) == ℯ
@test log(exp(:pi)) == π
@test 2//Inf == 0
@test Inf//2 == Inf

@test (rcall("x"); Algebra.ws() == :x)
@test (operator(:x); operator(:x); clear(:x); true)
@test try det([:x :y]) catch; true end
@test join([R"1",R"1"]) == R"1;1"
@test list([R"1",R"x"]) == list((1,:x))
@test Reduce.lister(:x) == R"x"
@test !latex(false)
@test (@rounded @factor x^2-2x+1) == :((x-1)^2)
@test (@rounded @off_factor @rcall x^2) == :(x^2)
@test det([:a :b; :c :d]) == :(a*d-b*c)
@test tp([:a;:b]) == [:a :b;]
@test transpose(:x) == :x
@test adjoint(:x) == :(repart(x)-impart(x)*im)
operator(:cbrt)
@test (rlet(:(cbrt(~x))=>:(x^(1/3))); true)
@test (rlet([:(cbrt(~x))=>:(x^(1/3))]); true)
@test (rlet(Dict(:(cbrt(~x))=>:(x^(1/3)))); true)
@test 1+R"x" == :(x+1)
@test inv([:a :b; :c :d]) |> typeof <: Array
@test inv(1) == 1.0
@test [1 2; 3 4]\[1,2] == [0.0,0.5]
@test 1.0//1.0 == 1
@test [:a :b; :c :d]/2 |> typeof <: Array
@test [:a :b; :c :d]+1 == 1+[:a :b; :c :d]
@test [:x] + 1 == 1 + [:x]
@test [:x,:y]' + 1 == 1 + [:x,:y]'
@test solve(:(x-1),:x) == solve((:(x-1),),:x)
@test (order(nothing); korder(nothing); true)
@test (Reduce.@subtype FakeReal <: Real; FakeReal(:(x+1)) + FakeReal(:y) == FakeReal(:(x+1+y)))


text:
using LinearAlgebra
using Fermi.DIIS

# Functions specific to DF-CTF
#
############ PREPROCESSING ################
"""
    Fermi.CoupledCluster.RCCSD{T}(Alg::CTF)

Compute a RCCSD wave function using the Compiled time factorization algorithm (CTF)
"""
function RCCSD{T}(guess::RCCSD{Tb},Alg::DFCTF) where { T <: AbstractFloat,
                                                    Tb <: AbstractFloat }
    refwfn = Fermi.HartreeFock.RHF()

    drop_occ = Fermi.CurrentOptions["drop_occ"]
    drop_vir = Fermi.CurrentOptions["drop_vir"]

    ints = refwfn.ints
    Fermi.Integrals.aux_ri!(ints)
    RCCSD{T}(refwfn, guess, ints, Alg)
end

function RCCSD{T}(refwfn::RHF, guess::RCCSD{Tb}, ints::IntegralHelper{Tc}, Alg::DFCTF) where { T <: AbstractFloat,
                                                                                Tb <: AbstractFloat,
                                                                                Tc <: AbstractFloat }
    d = [i - a for i = diag(ints["FOO"]), a = diag(ints["FVV"])]
    D = [i + j - a - b for i = diag(ints["FOO"]), j = diag(ints["FOO"]), a = diag(ints["FVV"]), b = diag(ints["FVV"])]
    newT1 = ints["FOV"]./d
    Bov = ints["BOV"]
    @tensor oovv[i,j,a,b] := Bov[Q,i,a]*Bov[Q,j,b]
    newT2 = oovv ./ D
    RCCSD{T}(refwfn, ints, newT1, newT2, Alg)
end
############ END PREPROCESSING ################
#
#
#
############# KERNEL FUNCTIONS ###################

function print_alg(Alg::DFCTF)
    @output "\n    • Computing CCSD with the DF-CTF algorithm .\n\n"
end

function compute_integrals(ints,Alg::DFCTF)
    Fermi.Integrals.aux_ri!(ints)
    @output "Aux basis: {}\n" ints.bname["aux"]
    ints["BOV"]
    ints["BOO"]
    ints["BVV"]
end

"""
    Fermi.Coupled Cluster.update_energy(T1::Array{T, 2}, T2::Array{T, 4}, f::Array{T,2}, Voovv::Array{T, 4}) where T <: AbstractFloat

Compute CC energy from amplitudes and integrals.
"""
function update_energy(T1::Array{T, 2}, T2::Array{T, 4}, f::Array{T,2}, ints, alg::DFCTF) where { T <: AbstractFloat }

    Bov = ints["BOV"]
    @tensoropt (k=>x, l=>x, c=>100x, d=>100x)  begin
        CC_energy = 2.0*f[k,c]*T1[k,c]
        B[l,c,k,d] := -1.0*T1[l,c]*T1[k,d]
        B[l,c,k,d] += -1.0*T2[l,k,c,d]
        B[l,c,k,d] += 2.0*T2[k,l,c,d]
        CC_energy += B[l,c,k,d]*Bov[Q,k,c]*Bov[Q,l,d]#Voovv[k,l,c,d]
        CC_energy += 2.0*T1[l,c]*T1[k,d]*Bov[Q,l,c]*Bov[Q,k,d]#Voovv[l,k,c,d]
    end

    return CC_energy
end
function compute_oovv(ints,alg::DFCTF)
    Bov = ints["BOV"]
    @tensor oovv[i,j,a,b] := Bov[Q,i,a]*Bov[Q,j,b]
    oovv
end

function update_T1(T1::Array{T,2}, T2::Array{T,4}, newT1::Array{T,2}, foo, fov, fvv, ints::IntegralHelper, alg::DFCTF) where { T <: AbstractFloat }
    Bov = ints["BOV"]
    Boo = ints["BOO"]
    Bvv = ints["BVV"]
    @tensor begin
        Voooo[i,j,k,l] := Boo[Q,i,k]*Boo[Q,j,l]
        Vooov[i,j,k,a] := Boo[Q,i,k]*Bov[Q,j,a]
        Voovv[i,j,a,b] := Bov[Q,i,a]*Bov[Q,j,b]
        Vovov[i,a,j,b] := Boo[Q,i,j]*Bvv[Q,a,b]
    end
    @tensoropt (i=>x, j=>x, k=>x, l=>x, a=>10x, b=>10x, c=>10x, d=>10x,Q=>50x) begin
        newT1[i,a] += fov[i,a]
        newT1[i,a] -= foo[i,k]*T1[k,a]
        newT1[i,a] += fvv[c,a]*T1[i,c]
        newT1[i,a] -= fov[k,c]*T1[i,c]*T1[k,a]
        newT1[i,a] += 2.0*fov[k,c]*T2[i,k,a,c]
        newT1[i,a] -= fov[k,c]*T2[k,i,a,c]
        newT1[i,a] -= T1[k,c]*Vovov[i,c,k,a]
        newT1[i,a] += 2.0*T1[k,c]*Voovv[k,i,c,a]
        newT1[i,a] -= T2[k,i,c,d]*Bov[Q,k,d]*Bvv[Q,a,c]
        newT1[i,a] += 2.0*T2[i,k,c,d]*Bov[Q,k,d]*Bvv[Q,a,c]
        newT1[i,a] += -2.0*T2[k,l,a,c]*Vooov[k,l,i,c]
        newT1[i,a] += T2[l,k,a,c]*Vooov[k,l,i,c]
        newT1[i,a] += -2.0*T1[k,c]*T1[l,a]*Vooov[l,k,i,c]
        newT1[i,a] -= T1[k,c]*T1[i,d]*Bov[Q,k,d]*Bvv[Q,a,c]
        newT1[i,a] += 2.0*T1[k,c]*T1[i,d]*Bov[Q,k,c]*Bvv[Q,a,d]
        newT1[i,a] += T1[k,c]*T1[l,a]*Vooov[k,l,i,c]
        newT1[i,a] += -2.0*T1[k,c]*T2[i,l,a,d]*Voovv[l,k,c,d]
        newT1[i,a] += -2.0*T1[k,c]*T2[l,i,a,d]*Voovv[k,l,c,d]
        newT1[i,a] += T1[k,c]*T2[l,i,a,d]*Voovv[l,k,c,d]
        newT1[i,a] += -2.0*T1[i,c]*T2[l,k,a,d]*Voovv[l,k,c,d]
        newT1[i,a] += T1[i,c]*T2[l,k,a,d]*Voovv[k,l,c,d]
        newT1[i,a] += -2.0*T1[l,a]*T2[i,k,d,c]*Voovv[k,l,c,d]
        newT1[i,a] += T1[l,a]*T2[i,k,c,d]*Voovv[k,l,c,d]
        newT1[i,a] += T1[k,c]*T1[i,d]*T1[l,a]*Voovv[l,k,c,d]
        newT1[i,a] += -2.0*T1[k,c]*T1[i,d]*T1[l,a]*Voovv[k,l,c,d]
        newT1[i,a] += 4.0*T1[k,c]*T2[i,l,a,d]*Voovv[k,l,c,d]
    end
end

function update_T2(T1::Array{T,2},T2::Array{T,4},newT2::Array{T,4},foo,fov,fvv,ints::IntegralHelper,alg::DFCTF) where T <: AbstractFloat
    Bov = ints["BOV"]
    Boo = ints["BOO"]
    Bvv = ints["BVV"]
    @tensor begin
        Voooo[i,j,k,l] := Boo[Q,i,k]*Boo[Q,j,l]
        Vooov[i,j,k,a] := Boo[Q,i,k]*Bov[Q,j,a]
        Voovv[i,j,a,b] := Bov[Q,i,a]*Bov[Q,j,b]
        Vovov[i,a,j,b] := Boo[Q,i,j]*Bvv[Q,a,b]
    end
    dfsz,nocc,nvir = size(Bov)
    @tensor _T[i,j,c,d] := T1[i,c]*T1[j,d] + T2[i,j,c,d]
    T2_inter = zeros(T,dfsz,nvir,nvir)
    T2_slice = zeros(T,nvir,nvir)
    Tslice = zeros(T,nvir,nvir)
    for i=1:nocc, j=1:nocc
        T2_inter .= 0
        T2_slice .= 0
        @views Tslice .= _T[i,j,:,:]
        Fermi.contract!(T2_inter,Tslice,Bvv,"Qad","cd","Qca")
        Fermi.contract!(T2_slice,T2_inter,Bvv,"ab","Qad","Qdb")
        newT2[i,j,:,:] += T2_slice
    end

    @tensoropt (i=>x, j=>x, k=>x, l=>x, a=>10x, b=>10x, c=>10x, d=>10x, Q=>50x) begin
        newT2[i,j,a,b] += Voovv[i,j,a,b]
        newT2[i,j,a,b] += T1[k,a]*T1[l,b]*Voooo[i,j,k,l]
        newT2[i,j,a,b] += T2[k,l,a,b]*Voooo[i,j,k,l]
        newT2[i,j,a,b] -= T1[i,c]*T1[j,d]*T1[k,a]*Bov[Q,k,c]*Bvv[Q,b,d]
        newT2[i,j,a,b] -= T1[i,c]*T1[j,d]*T1[k,b]*Bov[Q,k,d]*Bvv[Q,a,c]
        newT2[i,j,a,b] += T1[i,c]*T1[k,a]*T1[l,b]*Vooov[l,k,j,c]
        newT2[i,j,a,b] += T1[j,c]*T1[k,a]*T1[l,b]*Vooov[k,l,i,c]
        newT2[i,j,a,b] += T2[k,l,a,c]*T2[i,j,d,b]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += -2.0*T2[i,k,a,c]*T2[l,j,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += -2.0*T2[l,k,a,c]*T2[i,j,d,b]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T2[k,i,a,c]*T2[l,j,d,b]*Voovv[l,k,c,d]
        newT2[i,j,a,b] += T2[i,k,a,c]*T2[l,j,b,d]*Voovv[l,k,c,d]
        newT2[i,j,a,b] += -2.0*T2[i,k,a,c]*T2[j,l,b,d]*Voovv[l,k,c,d]
        newT2[i,j,a,b] += T2[k,i,a,c]*T2[l,j,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += -2.0*T2[k,i,a,c]*T2[j,l,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T2[i,j,a,c]*T2[l,k,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += -2.0*T2[i,j,a,c]*T2[k,l,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T2[k,j,a,c]*T2[i,l,d,b]*Voovv[l,k,c,d]
        newT2[i,j,a,b] += 4.0*T2[i,k,a,c]*T2[j,l,b,d]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T2[i,j,d,c]*T2[l,k,a,b]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T1[i,c]*T1[j,d]*T1[k,a]*T1[l,b]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += T1[i,c]*T1[j,d]*T2[l,k,a,b]*Voovv[l,k,c,d]
        newT2[i,j,a,b] += T1[k,a]*T1[l,b]*T2[i,j,d,c]*Voovv[l,k,c,d]
        P_OoVv[i,j,a,b] := -1.0*foo[i,k]*T2[k,j,a,b]
        P_OoVv[i,j,a,b] += fvv[c,a]*T2[i,j,c,b]
        P_OoVv[i,j,a,b] += -1.0*T1[k,b]*Vooov[j,i,k,a]
        P_OoVv[i,j,a,b] += T1[j,c]*Bov[Q,i,a]*Bvv[Q,c,b]
        P_OoVv[i,j,a,b] += -1.0*fov[k,c]*T1[i,c]*T2[k,j,a,b]
        P_OoVv[i,j,a,b] += -1.0*fov[k,c]*T1[k,a]*T2[i,j,c,b]
        P_OoVv[i,j,a,b] += -1.0*T2[k,i,a,c]*Voovv[k,j,c,b]
        P_OoVv[i,j,a,b] += -1.0*T1[i,c]*T1[k,a]*Voovv[k,j,c,b]
        P_OoVv[i,j,a,b] += -1.0*T1[i,c]*T1[k,b]*Vovov[j,c,k,a]
        P_OoVv[i,j,a,b] += 2.0*T2[i,k,a,c]*Voovv[k,j,c,b]
        P_OoVv[i,j,a,b] += -1.0*T2[i,k,a,c]*Vovov[j,c,k,b]
        P_OoVv[i,j,a,b] += -1.0*T2[k,j,a,c]*Vovov[i,c,k,b]
        P_OoVv[i,j,a,b] += -2.0*T1[l,b]*T2[i,k,a,c]*Vooov[l,k,j,c]
        P_OoVv[i,j,a,b] += T1[l,b]*T2[k,i,a,c]*Vooov[l,k,j,c]
        P_OoVv[i,j,a,b] += -1.0*T1[j,c]*T2[i,k,d,b]*Bov[Q,k,c]*Bvv[Q,a,d]
        P_OoVv[i,j,a,b] += -1.0*T1[j,c]*T2[k,i,a,d]*Bov[Q,k,d]*Bvv[Q,b,c]
        P_OoVv[i,j,a,b] += -1.0*T1[j,c]*T2[i,k,a,d]*Bov[Q,k,c]*Bvv[Q,b,d]
        P_OoVv[i,j,a,b] += 2.0*T1[k,c]*T2[i,j,a,d]*Bov[Q,k,c]*Bvv[Q,b,d]
        P_OoVv[i,j,a,b] += T1[j,c]*T2[l,k,a,b]*Vooov[l,k,i,c]
        P_OoVv[i,j,a,b] += T1[l,b]*T2[i,k,a,c]*Vooov[k,l,j,c]
        P_OoVv[i,j,a,b] += -1.0*T1[k,a]*T2[i,j,d,c]*Bov[Q,k,d]*Bvv[Q,b,c]
        P_OoVv[i,j,a,b] += T1[k,a]*T2[i,l,c,b]*Vooov[l,k,j,c]
        P_OoVv[i,j,a,b] += 2.0*T1[j,c]*T2[i,k,a,d]*Bov[Q,k,d]*Bvv[Q,b,c]
        P_OoVv[i,j,a,b] -= 1.0*T1[k,c]*T2[i,j,a,d]*Bov[Q,k,d]*Bvv[Q,b,c]
        P_OoVv[i,j,a,b] += T1[k,c]*T2[i,l,a,b]*Vooov[k,l,j,c]
        P_OoVv[i,j,a,b] += -2.0*T1[k,c]*T2[i,l,a,b]*Vooov[l,k,j,c]
        P_OoVv[i,j,a,b] += T2[j,k,c,d]*T2[i,l,a,b]*Voovv[k,l,c,d]
        P_OoVv[i,j,a,b] += -2.0*T1[k,c]*T1[j,d]*T2[i,l,a,b]*Voovv[k,l,c,d]
        P_OoVv[i,j,a,b] += T1[k,c]*T1[j,d]*T2[i,l,a,b]*Voovv[l,k,c,d]
        P_OoVv[i,j,a,b] += -2.0*T1[k,c]*T1[l,a]*T2[i,j,d,b]*Voovv[k,l,c,d]
        P_OoVv[i,j,a,b] += T1[k,c]*T1[l,a]*T2[i,j,d,b]*Voovv[l,k,c,d]
        P_OoVv[i,j,a,b] += T1[i,c]*T1[k,a]*T2[l,j,b,d]*Voovv[k,l,c,d]
        P_OoVv[i,j,a,b] += -2.0*T1[i,c]*T1[k,a]*T2[j,l,b,d]*Voovv[k,l,c,d]
        P_OoVv[i,j,a,b] += T1[i,c]*T1[k,a]*T2[l,j,d,b]*Voovv[l,k,c,d]
        P_OoVv[i,j,a,b] += T1[i,c]*T1[l,b]*T2[k,j,a,d]*Voovv[k,l,c,d]
        P_OoVv[i,j,a,b] += -2.0*T2[i,k,d,c]*T2[l,j,a,b]*Voovv[k,l,c,d]
        newT2[i,j,a,b] += P_OoVv[i,j,a,b] + P_OoVv[j,i,b,a]
    end
end
#
#
############# END KERNEL FUNCTIONS ###################


text:
module AntOptim

using JuMP

include("Parser.jl")
include("Model/Model.jl")

include("Ants.jl")

end


text:
promptname(level, name) = "$level|$name> "
function RunDebugger(frame, repl = nothing, terminal = nothing; initial_continue=false)
    if repl === nothing
        if !isdefined(Base, :active_repl)
            error("Debugger.jl needs to be run in a Julia REPL")
        end
        repl = Base.active_repl
    end
    if !isa(repl, REPL.LineEditREPL)
        error("Debugger.jl requires a LineEditREPL type of REPL")
    end

    if terminal === nothing
        terminal = Base.active_repl.t
    end
    state = DebuggerState(; frame=frame, repl=repl, terminal=terminal)

    # Setup debug panel
    normal_prefix = Sys.iswindows() ? "\e[33m" : "\e[38;5;166m"
    compiled_prefix = "\e[96m"
    panel = LineEdit.Prompt(promptname(state.level, "debug");
        prompt_prefix = () -> state.mode == Compiled() ? compiled_prefix : normal_prefix,
        prompt_suffix = Base.text_colors[:normal],
        on_enter = s->true)

    panel.hist = REPL.REPLHistoryProvider(Dict{Symbol,Any}(:debug => panel))
    REPL.history_reset_state(panel.hist)

    search_prompt, skeymap = LineEdit.setup_search_keymap(panel.hist)
    search_prompt.complete = REPL.LatexCompletions()

    state.main_mode = panel

    panel.on_done = (s,buf,ok)->begin
        line = String(take!(buf))
        old_level = state.level
        if !ok || strip(line) == "q"
            LineEdit.transition(s, :abort)
            LineEdit.reset_state(s)
            return false
        end
        if length(panel.hist.history) == 0
            printstyled(stderr, "no previous command executed\n"; color=Base.error_color())
            return false
        end
        if isempty(strip(line))
            command = panel.hist.history[end]
        else
            command = strip(line)
        end
        do_print_status = true
        cmd1 = split(command,' ')[1]
        do_print_status = try
            execute_command(state, Val{Symbol(cmd1)}(), command)
        catch err
            # This will only show the stacktrae up to the current frame because
            # currently, the unwinding in JuliaInterpreter unlinks the frames to
            # where the error is thrown

            # Buffer error printing
            io = IOContext(IOBuffer(), Base.pipe_writer(terminal))
            Base.display_error(io, err, JuliaInterpreter.leaf(state.frame))
            print(Base.pipe_writer(terminal), String(take!(io.io)))
            # Comment below out if you are debugging the Debugger
            #Base.display_error(Base.pipe_writer(terminal), err, catch_backtrace())
            LineEdit.transition(s, :abort)
            LineEdit.reset_state(s)
           return false
        end
        if old_level != state.level
            panel.prompt = promptname(state.level, "debug")
        end
        LineEdit.reset_state(s)
        if state.frame === nothing
            LineEdit.transition(s, :abort)
            LineEdit.reset_state(s)
            return false
        end
        if do_print_status
            print_status(Base.pipe_writer(terminal), active_frame(state); force_lowered = state.lowered_status)
        end
        return true
    end

    repl_switch = Dict{Any,Any}(
        '`' => function (s,args...)
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                prompt = julia_prompt(state)
                buf = copy(LineEdit.buffer(s))
                LineEdit.transition(s, prompt) do
                    LineEdit.state(s, prompt).input_buffer = buf
                end
            else
                LineEdit.edit_insert(s, '`')
            end
        end,
        'C' => function (s, args...)
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                toggle_mode(state)
                write(state.terminal, '\r')
                LineEdit.write_prompt(state.terminal, panel)
            else
                LineEdit.edit_insert(s, "C")
            end
        end,
        'L' => function (s, args...)
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                toggle_lowered(state)
                println(Base.pipe_writer(terminal))
                print_status(Base.pipe_writer(terminal), active_frame(state); force_lowered=state.lowered_status)
                LineEdit.write_prompt(state.terminal, panel)
            else
                LineEdit.edit_insert(s, "L")
            end
        end,
        '+' => function (s, args...)
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                NUM_SOURCE_LINES_UP_DOWN[] += 1
                println(Base.pipe_writer(terminal))
                print_status(Base.pipe_writer(terminal), active_frame(state); force_lowered=state.lowered_status)
                LineEdit.write_prompt(state.terminal, panel)
            else
                LineEdit.edit_insert(s, "+")
            end
        end,
        '-' => function (s, args...)
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                NUM_SOURCE_LINES_UP_DOWN[] = max(1, NUM_SOURCE_LINES_UP_DOWN[] - 1)
                println(Base.pipe_writer(terminal))
                print_status(Base.pipe_writer(terminal), active_frame(state); force_lowered=state.lowered_status)
                LineEdit.write_prompt(state.terminal, panel)
            else
                LineEdit.edit_insert(s, "-")
            end
        end
    )

    state.standard_keymap = Dict{Any,Any}[skeymap, LineEdit.history_keymap, LineEdit.default_keymap, LineEdit.escape_defaults]
    panel.keymap_dict = LineEdit.keymap([repl_switch;state.standard_keymap])

    if initial_continue
        try
            execute_command(state, Val(:c), "c")
        catch err
            # Buffer error printing
            io = IOContext(IOBuffer(), Base.pipe_writer(terminal))
            Base.display_error(io, err, JuliaInterpreter.leaf(state.frame))
            print(Base.pipe_writer(terminal), String(take!(io.io)))
            return
        end
        state.frame === nothing && return state.overall_result
    end
    print_status(Base.pipe_writer(terminal), active_frame(state); force_lowered=state.lowered_status)
    REPL.run_interface(terminal, LineEdit.ModalInterface([panel,search_prompt]))

    return state.overall_result
end

function julia_prompt(state::DebuggerState)
    # Return early if this has already been called on the state
    isassigned(state.julia_prompt) && return state.julia_prompt[]

    julia_prompt = LineEdit.Prompt(() -> promptname(state.level, "julia");
        # Copy colors from the prompt object
        prompt_prefix = state.repl.prompt_color,
        prompt_suffix = (state.repl.envcolors ? Base.input_color : state.repl.input_color),
        complete = DebugCompletionProvider(state),
        on_enter = REPL.return_callback)
    julia_prompt.hist = state.main_mode.hist
    julia_prompt.hist.mode_mapping[:julia] = julia_prompt

    julia_prompt.on_done = (s,buf,ok)->begin
        if !ok
            LineEdit.transition(s, :abort)
            return false
        end
        xbuf = copy(buf)
        command = String(take!(buf))
        @static if VERSION >= v"1.2.0-DEV.253"
            response = _eval_code(active_frame(state), command)
            REPL.print_response(state.repl, response, true, true)
        else
            ok, result = _eval_code(active_frame(state), command)
            REPL.print_response(state.repl, ok ? result : result[1], ok ? nothing : result[2], true, true)
        end
        println(state.terminal)
        LineEdit.reset_state(s)
    end
    julia_prompt.keymap_dict = LineEdit.keymap([REPL.mode_keymap(state.main_mode); state.standard_keymap])
    state.julia_prompt[] = julia_prompt
    return julia_prompt
end

@static if VERSION >= v"1.2.0-DEV.253"
    function _eval_code(frame::Frame, code::AbstractString)
        try
            return eval_code(frame, code), false
        catch
            return Base.catch_stack(), true
        end
    end
else
    function _eval_code(frame::Frame, code::AbstractString)
        try
            return true, eval_code(frame, code)
        catch err
            return false, (err, catch_backtrace())
        end
    end
end

function eval_code(frame::Frame, command::AbstractString)
    expr = Base.parse_input_line(command)
    isexpr(expr, :toplevel) && (expr = expr.args[end])

    if isexpr(expr, :toplevel)
      expr = Expr(:block, expr.args...)
    end
    # see https://github.com/JuliaLang/julia/issues/31255 for the Symbol("") check
    vars = filter(v -> v.name != Symbol(""), JuliaInterpreter.locals(frame))
    res = gensym()
    eval_expr = Expr(:let,
        Expr(:block, map(x->Expr(:(=), x...), [(v.name, maybe_quote(v.value)) for v in vars])...),
        Expr(:block,
            Expr(:(=), res, expr),
            Expr(:tuple, res, Expr(:tuple, [v.name for v in vars]...))
        ))
    eval_res, res = Core.eval(moduleof(frame), eval_expr)
    j = 1
    for (i, v) in enumerate(vars)
        if v.isparam
            frame.framedata.sparams[j] = res[i]
            j += 1
        else
            frame.framedata.locals[frame.framedata.last_reference[v.name]] = Some{Any}(res[i])
        end
    end
    eval_res
end

# Completions

mutable struct DebugCompletionProvider <: REPL.CompletionProvider
    state::DebuggerState
end

function LineEdit.complete_line(c::DebugCompletionProvider, s)
    partial = REPL.beforecursor(s.input_buffer)
    full = LineEdit.input_string(s)
    ret, range, should_complete = completions(c, full, lastindex(partial))
    return unique!(map(REPLCompletions.completion_text, ret)), partial[range], should_complete
end

function completions(c::DebugCompletionProvider, full, partial)
    mod = moduleof(c.state.frame)
    ret, range, should_complete = REPLCompletions.completions(full, partial, mod)

    # TODO Add local variable completions
    return ret, range, should_complete
end


text:
using Catalyst, DiffEqBase, ModelingToolkit, Test

using ModelingToolkit: value

@parameters t k1 k2
@variables S(t) I(t) R(t)
rxs = [Reaction(k1, [S,I], [I], [1,1], [2]),
       Reaction(k2, [I], [R]) ]
@named rs = ReactionSystem(rxs, t, [S,I,R], [k1,k2])

specset = Set([value(S) =>1, value(I) => 2, value(R) => 3])
@test issetequal(specset, speciesmap(rs))

pset = Set([value(k1) => 1, value(k2) => 2])
@test issetequal(pset, paramsmap(rs))

rxs2 = [Reaction(k2, [I], [R], [1], [1]),
        Reaction(k1, [S,I], [I], [1,1], [2])]
@named rs2 = ReactionSystem(rxs2, t, [R,I,S], [k2,k1])
@test rs == rs2

rs3 = make_empty_network()
@parameters k3 k4
@variables D(t)
addspecies!(rs3, S)
addspecies!(rs3, D)
addparam!(rs3, k3)
addparam!(rs3, k4)
@test issetequal(species(rs3), [S, D])
@test issetequal(params(rs3), [k3, k4])
addreaction!(rs3, Reaction(k3, [S], [D]))
addreaction!(rs3, Reaction(k4, [S,I], [D]))
merge!(rs, rs3)
addspecies!(rs2, S)
addspecies!(rs2, D)
addparam!(rs2, k3)
addparam!(rs2, k4)
addreaction!(rs2, Reaction(k3, [S], [D]))
addreaction!(rs2, Reaction(k4, [S,I], [D]))
@test rs2 == rs

rxs = [Reaction(k1, [S,I], [I], [1,1], [2]),
       Reaction(k2, [I], [R]) ]
@named rs = ReactionSystem(rxs, t, [S,I,R], [k1,k2])
rs3 = make_empty_network()
addspecies!(rs3, S)
addspecies!(rs3, D)
addparam!(rs3, k3)
addparam!(rs3, k4)
addreaction!(rs3, Reaction(k3, [S], [D]))
addreaction!(rs3, Reaction(k4, [S,I], [D]))
rs4 = merge(rs, rs3)
@test rs2 == rs4

rxs = [Reaction(k1*S, [S,I], [I], [2,3], [2]),
       Reaction(k2*R, [I], [R]) ]
@named rs = ReactionSystem(rxs, t, [S,I,R], [k1,k2])
deps = dependents(rxs[2], rs)
@test isequal(deps, [R,I])
@test isequal(dependents(rxs[1], rs), dependants(rxs[1], rs))
addspecies!(rs, S)
@test numspecies(rs) == 3
addspecies!(rs, S, disablechecks=true)
@test numspecies(rs) == 4
addparam!(rs, k1)
@test numparams(rs) == 2
addparam!(rs, k1, disablechecks=true)
@test numparams(rs) == 3

rnmat = @reaction_network begin
       α, S + 2I --> 2I
       β, 3I --> 2R + S
   end α β

smat = [1 2 0;
        0 3 0]
pmat = [0 2 0;
        1 0 2]
@test all(smat .== substoichmat(rnmat))
@test all(pmat .== prodstoichmat(rnmat))

############## testing newly added intermediate complexes reaction networks##############
rns  = Vector{ReactionSystem}(undef,6)
# mass-action non-catalytic
rns[1] = @reaction_network begin
    k₁, 2A --> B
    k₂, A --> C
    k₃, C --> D
    k₄, B + D --> E
end k₁ k₂ k₃ k₄
Z = [2 0 1 0 0 0 0;
     0 1 0 0 0 1 0;
     0 0 0 1 0 0 0;
     0 0 0 0 1 1 0;
     0 0 0 0 0 0 1]
B = [-1 0 0 0;
      1 0 0 0;
      0 -1 0 0;
      0 1 -1 0;
      0 0 1 0;
      0 0 0 -1;
      0 0 0 1]
Δ = [-1 0 0 0; 0 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 0; 0 0 0 -1; 0 0 0 0]
rcs,B2 = reactioncomplexes(rns[1])
@test B == B2
@test Z == complexstoichmat(rns[1])
@test Δ == complexoutgoingmat(rns[1])

# mass-action rober
rns[2] = @reaction_network begin
    k₁, A --> B
    k₂, B + B --> C + B
    k₃, B + C --> A + C
end k₁ k₂ k₃
Z = [1 0 0 0 1;
      0 1 2 1 0;
      0 0 0 1 1]
B = [-1 0 0;
      1 0 0;
      0 -1 0;
      0 1 -1;
      0 0 1]
Δ = [-1 0 0; 0 0 0; 0 -1 0; 0 0 -1; 0 0 0]
rcs,B2 = reactioncomplexes(rns[2])
@test B == B2
@test Z == complexstoichmat(rns[2])
@test Δ == complexoutgoingmat(rns[2])

#  some rational functions as rates
rns[3] = @reaction_network begin
    k₁ , ∅ --> X₁
    ( k₂/(1 + X₁*X₂ + X₃*X₄ ), k₃/(1 + X₁*X₂ + X₃*X₄ )), 2X₁ + X₂ ↔ 3X₃ + X₄
    k₄, X₄ --> ∅
end k₁ k₂ k₃ k₄
Z = [0 1 2 0 0;
      0 0 1 0 0;
      0 0 0 3 0;
      0 0 0 1 1]
B = [-1 0 0 1;
      1 0 0 0;
      0 -1 1 0;
      0 1 -1 0;
      0 0 0 -1]
Δ = [-1 0 0 0; 0 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 -1]
rcs,B2 = reactioncomplexes(rns[3])
@test B == B2
@test Z == complexstoichmat(rns[3])
@test Δ == complexoutgoingmat(rns[3])

# repressilator
rns[4]  = @reaction_network begin
   hillr(P₃,α,K,n), ∅ --> m₁
   hillr(P₁,α,K,n), ∅ --> m₂
   hillr(P₂,α,K,n), ∅ --> m₃
   (δ,γ), m₁ ↔ ∅
   (δ,γ), m₂ ↔ ∅
   (δ,γ), m₃ ↔ ∅
   β, m₁ --> m₁ + P₁
   β, m₂ --> m₂ + P₂
   β, m₃ --> m₃ + P₃
   μ, P₁ --> ∅
   μ, P₂ --> ∅
   μ, P₃ --> ∅
end α K n δ γ β μ;
Z = [0 1 0 0 1 0 0 0 0 0;
      0 0 1 0 0 1 0 0 0 0;
      0 0 0 1 0 0 1 0 0 0;
      0 0 0 0 1 0 0 1 0 0;
      0 0 0 0 0 1 0 0 1 0;
      0 0 0 0 0 0 1 0 0 1]
B = [-1 -1 -1 1 -1 1 -1 1 -1 0 0 0 1 1 1;
         1 0 0 -1 1 0 0 0 0 -1 0 0 0 0 0;
         0 1 0 0 0 -1 1 0 0 0 -1 0 0 0 0;
       0 0 1 0 0 0 0 -1 1 0 0 -1 0 0 0;
       0 0 0 0 0 0 0 0 0 1 0 0 0 0 0;
       0 0 0 0 0 0 0 0 0 0 1 0 0 0 0;
       0 0 0 0 0 0 0 0 0 0 0 1 0 0 0;
       0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0;
       0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0;
       0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1]
Δ = [-1 -1 -1 0 -1 0 -1 0 -1 0 0 0 0 0 0; 0 0 0 -1 0 0 0 0 0 -1 0 0 0 0 0; 0 0 0 0 0 -1 0 0 0 0 -1 0 0 0 0;
       0 0 0 0 0 0 0 -1 0 0 0 -1 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
       0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1]
rcs,B2 = reactioncomplexes(rns[4])
@test B == B2
@test Z == complexstoichmat(rns[4])
@test Δ == complexoutgoingmat(rns[4])

#brusselator
rns[5] = @reaction_network begin
   A, ∅ → X
   1, 2X + Y → 3X
   B, X → Y
   1, X → ∅
end A B
Z = [0 1 2 3 0;
      0 0 1 0 1]
B = [-1 0 0 1;
      1 0 -1 -1;
      0 -1 0 0;
      0 1 0 0;
      0 0 1 0]
Δ = [-1 0 0 0; 0 0 -1 -1; 0 -1 0 0; 0 0 0 0; 0 0 0 0]
rcs,B2 = reactioncomplexes(rns[5])
@test B == B2
@test Z == complexstoichmat(rns[5])
@test Δ == complexoutgoingmat(rns[5])

# some rational functions as rates
rns[6] = @reaction_network begin
    (k₁, k₋₁), X₁ + X₂ <--> X₃ + 2X₄
    (k₂/(1 + X₄*X₅ + X₆*X₇), k₋₂/(1 + X₄*X₅ + X₆*X₇)), 3X₄ + X₅ <--> X₆ + X₇
    (k₃/(1 + X₇ + X₈ + X₉ + X₁₀) , k₋₃/(1 + X₇ + X₈ + X₉ + X₁₀ )), 5X₇ + X₈ <--> X₉ + X₁₀
end k₁ k₋₁ k₂ k₋₂ k₃ k₋₃
Z = [1 0 0 0 0 0;
      1 0 0 0 0 0;
      0 1 0 0 0 0;
      0 2 3 0 0 0;
      0 0 1 0 0 0;
      0 0 0 1 0 0;
      0 0 0 1 5 0;
      0 0 0 0 1 0;
      0 0 0 0 0 1;
      0 0 0 0 0 1]
B = [-1 1 0 0 0 0;
      1 -1 0 0 0 0;
      0 0 -1 1 0 0;
      0 0 1 -1 0 0;
      0 0 0 0 -1 1;
      0 0 0 0 1 -1]
Δ = [-1 0 0 0 0 0; 0 -1 0 0 0 0; 0 0 -1 0 0 0; 0 0 0 -1 0 0; 0 0 0 0 -1 0; 0 0 0 0 0 -1]
rcs,B2 = reactioncomplexes(rns[6])
@test all(B .== B2)
@test Z == complexstoichmat(rns[6]; rcs=rcs)
@test Δ == complexoutgoingmat(rns[6], B=B2)

reaction_networks_standard = Vector{ReactionSystem}(undef,10)
reaction_networks_hill = Vector{ReactionSystem}(undef,10)
reaction_networks_real = Vector{ReactionSystem}(undef,3)
### some more networks to test whether Z*B == netstoichmat(rn)' or not. ###
reaction_networks_standard[1] = @reaction_network begin
    (p1,p2,p3), ∅ → (X1,X2,X3)
    (k1,k2), X2 ⟷ X1 + 2X3
    (k3,k4), X1 ⟷ X3
    (d1,d2,d3), (X1,X2,X3) → ∅
end p1 p2 p3 k1 k2 k3 k4 d1 d2 d3

reaction_networks_standard[2] = @reaction_network begin
    mmr(X2,v1,K1), ∅ → X1
    mm(X1,v2,K2), ∅ → X2
    d, X1+X2 → ∅
end v1 K1 v2 K2 d

reaction_networks_standard[3] = @reaction_network begin
    mm(X2,v1,K1), ∅ → X1
    mm(X3,v2,K2), ∅ → X2
    (k1,k2), X1 ⟷ X3
    (k3,k4), X3 + X2 ⟷ X4 + X1
    d, (X1,X2,X3,X4) → ∅
end v1 K1 v2 K2 k1 k2 k3 k4 d

reaction_networks_standard[4] = @reaction_network begin
    mmr(X4,v1,K1), ∅ → X1
    mmr(X1,v2,K2), ∅ → X2
    mmr(X2,v3,K3), ∅ → X3
    mmr(X3,v4,K4), ∅ → X4
    (d1,d2,d3,d4), (X1,X2,X3,X4) → ∅
end v1 K1 v2 K2 v3 K3 v4 K4 d1 d2 d3 d4

reaction_networks_standard[5] = @reaction_network begin
    p, ∅ → X1
    (k1,k2), X1 ⟷ X2
    (k3,k4), X2 ⟷ X3
    (k5,k6), X3 ⟷ X4
    d, X4 → ∅
end p k1 k2 k3 k4 k5 k6 d

reaction_networks_standard[6] = @reaction_network begin
    (p1,p2), ∅ → (X1,X2)
    (k1,k2), 2X1 ⟷ X3
    (k3,k4), X2 + X3 ⟷ 4X4
    (k5,k6), X4 + X1 ⟷ 2X3
    d, (X1,X2,X3,X4) → ∅
end p1 p2 k1 k2 k3 k4 k5 k6 d

reaction_networks_standard[7] = @reaction_network begin
    (p1,p2,p3), ∅ → (X1,X2,X3)
    (k1,k2), X1 + 2X2 ⟷ X4
    (mm(X3,v1,K1),k3), X4 ⟷ X5
    (d1,d2,d3,d4,d5), (X1,X2,X3,X4,X5) → ∅
end p1 p2 p3 k1 k2 k3 v1 K1 d1 d2 d3 d4 d5

reaction_networks_standard[8] = @reaction_network begin
    p, ∅ → 2X1
    k1, X1 → X2
    (k2, k3), X2 → X3
    d, X3 → ∅
end p k1 k2 k3 d

reaction_networks_standard[9] = @reaction_network begin
    (p1,p2,p3), ∅ ⟶ (X1,X2,X3)
    (d1,d2,d3), (X1,X2,X3) ⟶ ∅
    (k1,k2), X1 + X2 ⟷ X3
    (k3,k4), X3 ⟷ X4
    (k5,k6), X4 ⟷ X1 + X2
end p1 p2 p3 d1 d2 d3 k1 k2 k3 k4 k5 k6

reaction_networks_standard[10] = @reaction_network begin
    p, ∅ ⟶ X1
    (k1,k2), X1 → X2
    (k3,k4), X2 → X3
    (k5,k6), X3 → X4
    (k7,k8), X4 → X5
    d, X5  ⟶ ∅
end p k1 k2 k3 k4 k5 k6 k7 k8 d

### Network with hill functions ###.

reaction_networks_hill[1] = @reaction_network begin
    hillr(X2,v1,K1,n1),   ∅ → X1
    hillr(X1,v2,K2,n2),   ∅ → X2
    (d1,d2),          (X1,X2) → ∅
end v1 v2 K1 K2 n1 n2 d1 d2

reaction_networks_hill[2] = @reaction_network begin
    hillr(X3,v1,K1,n1),   ∅ → X1
    hillr(X1,v2,K2,n2),   ∅ → X2
    hillr(X2,v3,K3,n3),   ∅ → X3
    (d1,d2,d3),          (X1,X2,X3) → ∅
end v1 v2 v3 K1 K2 K3 n1 n2 n3 d1 d2 d3

reaction_networks_hill[3] = @reaction_network begin
    hillr(X2,v1,K1,n1), ∅ → X1
    hill(X1,v2,K2,n2), ∅ → X2
    d, X1+X2 → ∅
end v1 K1 n1 v2 K2 n2 d

reaction_networks_hill[4] = @reaction_network begin
    hillr(X2,v1,K1,n1)*hillr(X3,v1,K1,n1), ∅ → X1
    hillr(X1,v2,K2,n2)*hillr(X3,v2,K2,n2), ∅ → X2
    hillr(X1,v3,K3,n3)*hillr(X2,v3,K3,n3), ∅ → X3
    (d1,d2,d3), (X1,X2,X3)  ⟶ ∅
end v1 K1 n1 v2 K2 n2 v3 K3 n3 d1 d2 d3

reaction_networks_hill[5] = @reaction_network begin
    hillr(X2,v1,K1,n1)*hill(X4,v1,K1,n1), ∅ → X1
    hill(X5,v2,K2,n2), ∅ → X2
    hill(X3,v3,K3,n3), ∅ → X3
    hillr(X1,v4,K4,n4), ∅ → X4
    hill(X2,v5,K5,n5), ∅ → X5
    (k1,k2), X2 ⟷ X1 + 2X4
    (k3,k4), X4 ⟷ X3
    (k5,k6), 3X5 + X1 ⟷ X2
    (d1,d2,d3,d4,d5), (X1,X2,X3,X4,X5)  ⟶ ∅
end v1 K1 n1 v2 K2 n2 v3 K3 n3 v4 K4 n4 v5 K5 n5 k1 k2 k3 k4 k5 k6 d1 d2 d3 d4 d5

reaction_networks_hill[6] = @reaction_network begin
    v/10+hill(X1,v,K,n), ∅ → X1
    d, X1 → ∅
end v K n d

reaction_networks_hill[7] = @reaction_network begin
    v/10 + hill(X1,v,K,n), ∅ → X1 + X2
    (k1,k2), X1 + X2 ↔ X3
    k3, X3 → X1
    d, (X1,X2,X3) → ∅
end v K n k1 k2 k3 d

reaction_networks_hill[8] = @reaction_network begin
    hill(X2,v1,K1,n1), ∅ → X1
    hillr(X1,v2,K2,n2)*hill(X3,v3,K3,n3), ∅ → X2
    hill(X2,v4,K4,n4), ∅ → X3
    (d1,d2,d3), (X1,X2,X3) → ∅
end v1 K1 n1 v2 K2 n2 v3 K3 n3 v4 K4 n4 d1 d2 d3

reaction_networks_hill[9] = @reaction_network begin
    hill(X1,v1,K1,n1)*hillr(X1,v2,K2,n2), ∅ → X1
    d, X1 → ∅
end v1 K1 n1 v2 K2 n2 d

reaction_networks_hill[10] = @reaction_network begin
    hill(X2,v1,K1,n1), ∅ → X1
    hillr(X4,v2,K2,n2), ∅ → X2
    (k1,k2), 2X1 + X2 ⟷ X3
    (k3,k4), 2X2 + X3 ⟷ X4
    (k5,k6), X1 + X2 + X3 + X4 ⟷ X5 + X6
    (d1,d2), (X5,X6) → ∅
end v1 K1 n1 v2 K2 n2 k1 k2 k3 k4 k5 k6 d1 d2

reaction_networks_real[1] = @reaction_network begin
    v0 + hill(σ,v,K,n), ∅ → (σ+A)
    deg, (σ,A,Aσ) → ∅
    (kB,kD), A + σ ↔ Aσ
    S*kC, Aσ → σ
end v0 v K n kD kB kC deg S;

# A cell cycle model
reaction_networks_real[2] = @reaction_network begin
  k1, 0 --> Y
  k2p, Y --> 0
  k2pp*P, Y --> 0
  (k3p+k3pp*A)/(J3+Po), Po-->P
  (k4*m)/(J4+P), Y + P --> Y + Po
end k1 k2p k2pp k3p k3pp A J3 k4 m J4
#@add_constraint cc_network P+Po=1
#reaction_networks_real[3] = cc_network

# A bistable switch
reaction_networks_real[3] = @reaction_network begin
    d,    (X,Y) → ∅
    hillr(Y,v1,K1,n1), ∅ → X
    hillr(X,v2,K2,n2), ∅ → Y
end d v1 K1 n1 v2 K2 n2

myrn = [reaction_networks_standard;reaction_networks_hill;reaction_networks_real]
for i in 1:length(myrn)
    local rcs,B = reactioncomplexes(myrn[i])
    local Z = complexstoichmat(myrn[i])
    @test Z*B == (netstoichmat(myrn[i])')
end