text:
"""
Copy of Cosmology.jl to remove Unitful dependence.

Original license:

The MIT License (MIT)

Copyright (c) 2013 Mike Nolta <mike@nolta.net>

Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the "Software"),
to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense,
and/or sell copies of the Software, and to permit persons to whom the
Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
DEALINGS IN THE SOFTWARE.
"""
module CosmologyTools

using Slomo.Constants: rho_crit_h2
export Ωm, ρcrit, ρm, default_cosmo

abstract type AbstractCosmology end
abstract type AbstractClosedCosmology <: AbstractCosmology end
abstract type AbstractFlatCosmology <: AbstractCosmology end
abstract type AbstractOpenCosmology <: AbstractCosmology end

struct FlatLCDM{T <: Real} <: AbstractFlatCosmology
    h::T
    Ω_Λ::T
    Ω_m::T
    Ω_r::T
end
FlatLCDM(h::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real) =
    FlatLCDM(promote(float(h), float(Ω_Λ), float(Ω_m), float(Ω_r))...)

a2E(c::FlatLCDM, a) = sqrt(c.Ω_r + c.Ω_m * a + c.Ω_Λ * a^4)

struct ClosedLCDM{T <: Real} <: AbstractClosedCosmology
    h::T
    Ω_k::T
    Ω_Λ::T
    Ω_m::T
    Ω_r::T
end
ClosedLCDM(h::Real, Ω_k::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real) =
    ClosedLCDM(promote(float(h), float(Ω_k), float(Ω_Λ), float(Ω_m),
                       float(Ω_r))...)

struct OpenLCDM{T <: Real} <: AbstractOpenCosmology
    h::T
    Ω_k::T
    Ω_Λ::T
    Ω_m::T
    Ω_r::T
end
OpenLCDM(h::Real, Ω_k::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real) =
    OpenLCDM(promote(float(h), float(Ω_k), float(Ω_Λ), float(Ω_m),
                     float(Ω_r))...)

function a2E(c::Union{ClosedLCDM,OpenLCDM}, a)
    a2 = a * a
    sqrt(c.Ω_r + c.Ω_m * a + (c.Ω_k + c.Ω_Λ * a2) * a2)
end

for c in ("Flat", "Open", "Closed")
    name = Symbol("$(c)WCDM")
    @eval begin
        struct $(name){T <: Real} <: $(Symbol("Abstract$(c)Cosmology"))
            h::T
            Ω_k::T
            Ω_Λ::T
            Ω_m::T
            Ω_r::T
            w0::T
            wa::T
        end
        function $(name)(h::Real, Ω_k::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real,
                         w0::Real, wa::Real)
            $(name)(promote(float(h), float(Ω_k), float(Ω_Λ), float(Ω_m),
                            float(Ω_r), float(w0), float(wa))...)
        end
    end
end

function WCDM(h::Real, Ω_k::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real, w0::Real, wa::Real)
    if Ω_k < 0
        ClosedWCDM(h, Ω_k, Ω_Λ, Ω_m, Ω_r, w0, wa)
    elseif Ω_k > 0
        OpenWCDM(h, Ω_k, Ω_Λ, Ω_m, Ω_r, w0, wa)
    else
        FlatWCDM(h, Ω_k, Ω_Λ, Ω_m, Ω_r, w0, wa)
    end
end

function a2E(c::Union{FlatWCDM,ClosedWCDM,OpenWCDM}, a)
    ade = exp((1 - 3 * (c.w0 + c.wa)) * log(a) + 3 * c.wa * (a - 1))
    sqrt(c.Ω_r + (c.Ω_m + c.Ω_k * a) * a + c.Ω_Λ * ade)
end

"""
    cosmology(;h = 0.69,
               Neff = 3.04,
               OmegaK = 0,
               OmegaM = 0.29,
               OmegaR = nothing,
               Tcmb = 2.7255,
               w0 = -1,
               wa = 0)
# Parameters
* `h` - Dimensionless Hubble constant
* `OmegaK` - Curvature density (Ω_k)
* `OmegaM` - Matter density (Ω_m)
* `OmegaR` - Radiation density (Ω_r)
* `Tcmb` - CMB temperature in Kelvin; used to compute Ω_γ
* `Neff` - Effective number of massless neutrino species; used to compute Ω_ν
* `w0` - CPL dark energy equation of state; `w = w0 + wa(1-a)`
* `wa` - CPL dark energy equation of state; `w = w0 + wa(1-a)`
# Examples
```jldoctest
julia> c = cosmology()
Cosmology.FlatLCDM{Float64}(0.69, 0.7099122024007928, 0.29, 8.77975992071536e-5)
julia> c = cosmology(OmegaK=0.1)
Cosmology.OpenLCDM{Float64}(0.69, 0.1, 0.6099122024007929, 0.29, 8.77975992071536e-5)
julia> c = cosmology(w0=-0.9, OmegaK=-0.1)
Cosmology.ClosedWCDM{Float64}(0.69, -0.1, 0.8099122024007929, 0.29, 8.77975992071536e-5, -0.9, 0.0)
```
"""
function cosmology(;h = 0.69,
                   Neff = 3.04,
                   OmegaK = 0,
                   OmegaM = 0.29,
                   OmegaR = nothing,
                   Tcmb = 2.7255,
                   w0 = -1,
                   wa = 0)

    if OmegaR === nothing
        OmegaG = 4.48131e-7 * Tcmb^4 / h^2
        OmegaN = Neff * OmegaG * (7 / 8) * (4 / 11)^(4 / 3)
        OmegaR = OmegaG + OmegaN
    end

    OmegaL = 1 - OmegaK - OmegaM - OmegaR

    if !(w0 == -1 && wa == 0)
        return WCDM(h, OmegaK, OmegaL, OmegaM, OmegaR, w0, wa)
    end

    if OmegaK < 0
        return ClosedLCDM(h, OmegaK, OmegaL, OmegaM, OmegaR)
        elseif OmegaK > 0
        return OpenLCDM(h, OmegaK, OmegaL, OmegaM, OmegaR)
        else
        return FlatLCDM(h, OmegaL, OmegaM, OmegaR)
    end
end

# hubble rate
scale_factor(z) = 1 / (1 + z)
E(c::AbstractCosmology, z) = (a = scale_factor(z); a2E(c, a) / a^2)

const planck18 = cosmology(h = 0.6766,
                           Neff = 3.046,
                           OmegaM = 0.3111)

const planck15 = cosmology(h = 0.6774,
                           Neff = 3.046,
                           OmegaM = 0.3089)

const planck13 = cosmology(h = 0.6777,
                           Neff = 3.046,
                           OmegaM = 0.3071)

const wmap9 = cosmology(h = 0.6932,
                        Neff = 3.046,
                        OmegaM = 0.2865)

default_cosmo = planck18

"""
Compute the matter density of the universe at redshift z, in units of the critical density.
"""
function Ωm(z; cosmo = default_cosmo)
    return cosmo.Ω_m * (1.0 + z) ^ 3 / E(cosmo, z) ^ 2
end

"""
Compute the critical density of the universe at redshift z, in Msun / kpc3
"""
function ρcrit(z; cosmo = default_cosmo)
    return rho_crit_h2 * cosmo.h^2 * E(cosmo, z) ^ 2
end

"""
Compute the mass density of the universe at redshift z, in Msun / kpc3
ρm(z) = Ωm(z) * ρcrit(z)
"""
function ρm(z; cosmo = default_cosmo)
    return cosmo.Ω_m * (1 + z)^3 * rho_crit_h2 * cosmo.h^2
end

end


text:
using MultivariatePolynomials
using JuMP
using PolyJuMP
using SumOfSquares
using DynamicPolynomials
using Mosek

# Create symbolic variables (not JuMP decision variables)
@polyvar x1 x2

# Create a Sum of Squares JuMP model with the Mosek solver
m = SOSModel(solver = MosekSolver())

# Create a JuMP decision variable for the lower bound
@variable m γ

# f(x) is the Goldstein-Price function
f1 = x1+x2+1
f2 = 19-14*x1+3*x1^2-14*x2+6*x1*x2+3*x2^2
f3 = 2*x1-3*x2
f4 = 18-32*x1+12*x1^2+48*x2-36*x1*x2+27*x2^2

f = (1+f1^2*f2)*(30+f3^2*f4)

# Constraints f(x) - γ to be sum of squares
@constraint m f >= γ

@objective m Max γ

status = solve(m)

# The lower bound found is 3
println(getobjectivevalue(m))

ext: jl


text:
#========== adjusting index ranges, runtime ==========#

# This is to get the range of j in A[2j], from axes(A,1):

function divrange(r::AbstractUnitRange, f::Integer)
    if f > 0
        # a = div(first(r), f, RoundUp) # onnly 1.4 it seems?
        a = cld(first(r), f)
        # z = div(last(r), f, RoundDown)
        z = fld(last(r), f)
    elseif f < 0
        # a = div(last(r), f, RoundUp)
        a = cld(last(r), f)
        # z = div(first(r), f, RoundDown)
        z = fld(first(r), f)
    else
        throw("can't scale indices by zero")
    end
    a:z
end

#=
divrange(1:10, 2) .* 2
divrange(0:10, 2) .* 2
divrange(1:11, 2) .* 2
divrange(1:10, 3) .* 3

divrange(1:10, -1) .* -1 |> sort
divrange(1:10, -2) .* -2 |> sort
divrange(0:10, -2) .* -2 |> sort
divrange(0:11, -2) .* -2 |> sort
=#

# Special case of A[-i]:

function minusrange(r::AbstractRange)
    -last(r):-first(r)
end

#=
minusrange(1:11) == divrange(1:11, -1)
minusrange(1:10) == divrange(1:10, -1)
=#

# This is to get the range of j in A[j+k], given axes(A,1) and the range of k

function subranges(r::AbstractUnitRange, s::AbstractRange)
    first(r)-minimum(s) : last(r)-maximum(s)
end

function addranges(r::AbstractUnitRange, s::AbstractRange)
    first(r)+maximum(s) : last(r)+minimum(s)
end

#=
issubset(subranges(1:10, 1:3) .+ 1, 1:10)
issubset(subranges(1:10, 1:3) .+ 3, 1:10)

issubset(addranges(1:10, 1:3) .- 1, 1:10)
issubset(addranges(1:10, 1:3) .- 3, 1:10)
=#

# # This is to get range of j in A[j÷2], from axes(A,1):

function mulrange(r::AbstractUnitRange, f::Integer)
    first(r)*f : last(r)*f
end

function dotdiv(r::AbstractUnitRange, f::Integer)
    first(r)÷f : last(r)÷f
end

#=
mulrange(1:10, 2) .÷ 2 |> unique
mulrange(0:10, 2) .÷ 2 |> unique
mulrange(1:11, 2) .÷ 2 |> unique

mulrange(1:10, 3) .÷ 3 |> unique
mulrange(-10:-1, 3) .÷ 3 |> unique
mulrange(-11:-1, 3) .÷ 3 |> unique
mulrange(-11:-2, 3) .÷ 3 |> unique
mulrange(-12:-3, 3) .÷ 3 |> unique
=#

# This is for A[I[j]] (where this range must be a subset of axes(A,1))
# and for A[I[j]+k] (where it enters into the calculation of k's range).

function extremerange(A)
    α, ω = minimum(A), maximum(A)
    α isa Integer && ω isa Integer || throw("expected integers!")
    α:ω
end

# This gives the range of j implied by A[i, pad(j,3)]

function padrange(r::AbstractUnitRange, lo::Integer, hi::Integer)
    first(r)-lo : last(r)+hi
end

#========== functions used by the macro ==========#

@nospecialize

# This is for the bounds check on A[I[j],k]:

function extremeview(ex::Expr)
    @assert ex.head == :ref
    A = ex.args[1]
    if any(is_const, ex.args[2:end])
        ind = map(i -> is_const(i) ? i : (:), ex.args[2:end])
        :(@view $A[$(ind...)])
    else
        A
    end
end

"""
    range_expr_walk(:(axes(A,1)), :(2i+1)) -> range, :i

Given the axis of `A`, and the expression inside `A[2i+1]`,
this returns an expression for the resulting range of index `i`.
Understands operations `+, -, *, ÷`.
(Don't really need `÷`, as this results in a non-`UnitRange`
which can't be a valid index.)

If the expression is from something like `A[2i+j]`, then it returns a tuple of ranges
and a tuple of symbols. The range for `:j` contains `:$(AXIS)i` and v-v.

If the expression is from `A[I[j]]` then it returns `(min:max, nothing)`,
and the caller should check `issubset(min:max, axes(A,1))`.
"""
function range_expr_walk(r, ex::Expr, con=[])
    ex.head == :kw && return range_expr_kw(r, ex)
    if ex.head == :ref # case of M[I[j], k] with r=axes(M,1)
        A = ex.args[1]
        A = extremeview(ex)
        push!(con, :(minimum($A) in $r && maximum($A) in $r || throw("not safe!"))) # not used??
        return (:($extremerange($A)),nothing)
    end
    ex.head == :call || throw("not sure what to do with $ex")
    if ex.args[1] == :pad && length(ex.args) == 3
        _, a, p = ex.args
        return range_expr_walk(:($padrange($r, $p, $p)), a)
    elseif ex.args[1] == :pad && length(ex.args) == 4
        _, a, lo, hi = ex.args
        return range_expr_walk(:($padrange($r, $lo, $hi)), a)
    elseif length(ex.args) == 2
        op, a = ex.args
        if op == :+
            return range_expr_walk(r, a)
        elseif op == :-
            return range_expr_walk(:($minusrange($r)), a)
        end
    elseif length(ex.args) == 3
        op, a, b = ex.args
        if op == :+
            is_const(a) && return range_expr_walk(:($r .- $a), b)
            is_const(b) && return range_expr_walk(:($r .- $b), a)
            # with neither constant, first go outwards from index j to expression b...
            ax_a = range_unwrap(a)
            ax_b = range_unwrap(b)
            #... then use that with given size(A,d) to constrain range of i, and v-v:
            range_a, i_a = range_expr_walk(:($subranges($r, $ax_b)), a)
            range_b, i_b = range_expr_walk(:($subranges($r, $ax_a)), b)
            return (range_a, range_b), (i_a, i_b)

        elseif op == :-
            is_const(a) && return range_expr_walk(:($minusrange($r .- $a)), b)
            is_const(b) && return range_expr_walk(:($r .+ $b), a)
            ax_a = range_unwrap(a)
            ax_b = range_unwrap(b)
            range_a, i_a = range_expr_walk(:($addranges($r, $ax_b)), a)
            range_b, i_b = range_expr_walk(:($minusrange($subranges($r, $ax_a))), b)
            return (range_a, range_b), (i_a, i_b)

        elseif op == :*
            is_const(a) && return range_expr_walk(:($divrange($r, $a)), b)
            is_const(b) && return range_expr_walk(:($divrange($r, $b)), a)
        elseif op == :÷
            is_const(b) && return range_expr_walk(:($mulrange($r, $b)), a)
        elseif op == :/
            throw("not sure what to do with $ex, perhaps you wanted ÷")
        end
    elseif length(ex.args) > 3
        op, a, b, c = ex.args[1:4]
        ds = ex.args[5:end]
        if op == :+
            is_const(a) && return range_expr_walk(:($r .- $a), :(+($b, $c, $(ds...))))
            is_const(b) && return range_expr_walk(:($r .- $b), :(+($a, $c, $(ds...))))
            is_const(c) && return range_expr_walk(:($r .- $c), :(+($a, $b, $(ds...))))
        end
    end
    throw("not sure what to do with $ex, sorry")
end

range_expr_walk(range, s::Symbol) = range, s
range_expr_walk(range, n::Integer) = range, nothing

is_const(::Int) = true
is_const(::Any) = false
is_const(s::Symbol) = s in [:(:), :begin, :end] # : for the purposes of saveconstraints setting :intersect
is_const(ex::Expr) = begin
    ex.head == :$ && return true # what's returned by range_expr_walk will still contain $
    if ex.head == :call && ex.args[1] in (:+, :-, :*, :÷)
        return all(is_const, ex.args[2:end])
    end
    false
end

"""
    range_expr_walk(:(axes(A,1)), :(i=j)) -> :(axes(A, :i)), :j

Special case for keyword indexing, `A[i=j, k=j+2]` comes here.
"""
function range_expr_kw(r::Expr, ex::Expr)
    @assert ex.head == :kw
    # @assert r.head == :call && r.args[1] == :axes
    r.args[3] = QuoteNode(ex.args[1])
    range_expr_walk(r, ex.args[2])
end

"""
    range_unwrap(:(2i+1)) -> :(2 .* AXIS_i .+ 1)

This goes in the opposite direction to `range_expr_walk`, and gives
the range of values taken by the expression, in terms of `Symbol($AXIS, i)`.
"""
range_unwrap(i::Symbol) = Symbol(AXIS, i)
range_unwrap(ex::Expr) = begin
    if ex.head == :ref # case of A[I[j]+k] comes here
        A = ex.args[1]
        return :($extremerange($A))
    end
    ex.head == :call || throw("don't know how to handle $ex")
    # if ex.args[1] == :pad or :clamp or :mod --- find test cases.
    if length(ex.args) == 2
        op, a = ex.args
        if op == :-
            return :($minusrange($(range_unwrap(a))))
        end
    elseif length(ex.args) == 3
        op, a, b = ex.args
        if op == :*
            a == -1 && return :($minusrange($(range_unwrap(b))))
            b == -1 && return :($minusrange($(range_unwrap(a))))
            is_const(a) && return :($a .* $(range_unwrap(b)))
            is_const(b) && return :($b .* $(range_unwrap(a)))
        elseif op == :+
            is_const(a) && return :($a .+ $(range_unwrap(b)))
            is_const(b) && return :($b .+ $(range_unwrap(a)))
        elseif op == :-
            is_const(a) && return :($a .- $(range_unwrap(b)))
            is_const(b) && return :($(range_unwrap(a)) .- $b)
        elseif op == :÷
            is_const(b) && return :($dotdiv($(range_unwrap(a)), $b))
        end
    end
    throw("don't know how to handle $ex, sorry")
end

"""
    range_expr_walk(nothing, :(i+2)) -> nothing, :j

Special case used for `A[mod(i+2)]` etc.
"""
range_expr_walk(r::Nothing, s::Symbol) = r, s
range_expr_walk(r::Nothing, n::Integer) = r, nothing
function range_expr_walk(r::Nothing, ex::Expr)
    is_const(ex) && return r, nothing
    ex.head == :ref && return r, nothing
    ex.head == :call || throw("not sure what to do with $ex")
    # if ex.args[1] in [:+, :-, :*, :÷]
        syms = []
        for x in ex.args[2:end]
            _, i = range_expr_walk(r, x)
            i isa Tuple && push!(syms, i...)
            i isa Symbol && push!(syms, i)
        end
        if isempty(syms)
            return r, nothing
        elseif length(syms) == 1
            return r, syms[1]
        else
            return r, Tuple(syms)
        end
    # end
end

@specialize

#========== the end ==========#


text:
# function for measurements and checkpointing
function measure(
    symmetry :: String,
    obs_file :: String,
    cp_file  :: String,
    Λ        :: Float64,
    dΛ       :: Float64,
    χ        :: Vector{Matrix{Float64}},
    χ_tol    :: NTuple{2, Float64},
    t        :: DateTime,
    t0       :: DateTime,
    r        :: Reduced_lattice,
    m        :: Mesh,
    a        :: Action,
    wt       :: Float64,
    ct       :: Float64
    )        :: Tuple{DateTime, Bool}

    # open files
    obs = h5open(obs_file, "cw")
    cp  = h5open(cp_file, "cw")

    # compute correlations
    χp = deepcopy(χ)
    compute_χ!(Λ, r, m, a, χ, χ_tol)

    # save correlations and self energy if respective dataset does not yet exist
    if haskey(obs, "χ/$(Λ)") == false
        save_self!(obs, Λ, m, a)
        save_χ!(obs, Λ, symmetry, m, χ)
    end

    # check for monotonicity of static part of dominant on-site correlation
    idx      = argmax(Float64[maximum(abs.(χ[i])) for i in eachindex(χ)])
    monotone = χ[idx][1, 1] / χp[idx][1, 1] >= 0.995

    # compute current run time (in hours)
    h0 = 1e-3 * (Dates.now() - t0).value / 3600.0

    # if more than half an hour is left to the wall time limit, use ct as timer heuristic for checkpointing
    if wt - h0 > 0.5
        # test if time limit for checkpoint (in hours) has been reached
        h = 1e-3 * (Dates.now() - t).value / 3600.0

        if h >= ct
            # generate checkpoint if it does not exist yet
            if haskey(cp, "a/$(Λ)") == false
                println();
                println("Generating checkpoint at cutoff Λ / |J| = $(Λ) ...")
                checkpoint!(cp, Λ, dΛ, m, a)
                println("Successfully generated checkpoint.")
                println();
            end

            # reset timer
            t = Dates.now()
        end
    # if less than half an hour is left to the wall time, generate checkpoints as if ct = 0. Lower bound to prevent cancellation during checkpoint writing
    elseif 0.1 < wt - h0 <= 0.5
        # generate checkpoint if it does not exist yet
        if haskey(cp, "a/$(Λ)") == false
            println(); println()
            println("Generating checkpoint at cutoff Λ / |J| = $(Λ) ...")
            checkpoint!(cp, Λ, dΛ, m, a)
            println("Successfully generated checkpoint.")
            println(); println()
        end
    end

    # close files
    close(obs)
    close(cp)

    return t, monotone
end

"""
    save_launcher!(
        path        :: String,
        f           :: String,
        name        :: String,
        size        :: Int64,
        model       :: String,
        symmetry    :: String,
        J           :: Vector{<:Any}
        ;
        A           :: Float64            = 0.0,
        S           :: Float64            = 0.5,
        β           :: Float64            = 1.0,
        euclidean   :: Bool               = false,
        num_σ       :: Int64              = 25,
        num_Ω       :: Int64              = 15,
        num_ν       :: Int64              = 10,
        num_χ       :: Int64              = 10,
        p_σ         :: NTuple{2, Float64} = (0.3, 1.0),
        p_Ωs        :: NTuple{5, Float64} = (0.3, 0.05, 0.10, 0.1, 50.0),
        p_νs        :: NTuple{5, Float64} = (0.3, 0.05, 0.10, 0.1, 50.0),
        p_Ωt        :: NTuple{5, Float64} = (0.3, 0.15, 0.20, 0.3, 50.0),
        p_νt        :: NTuple{5, Float64} = (0.3, 0.15, 0.20, 0.5, 50.0),
        p_χ         :: NTuple{5, Float64} = (0.3, 0.05, 0.10, 0.1, 50.0),
        lins        :: NTuple{5, Float64} = (5.0, 4.0, 8.0, 6.0, 4.0),
        bounds      :: NTuple{5, Float64} = (0.5, 50.0, 250.0, 150.0, 50.0),
        max_iter    :: Int64              = 10,
        min_eval    :: Int64              = 10,
        max_eval    :: Int64              = 100,
        Σ_tol       :: NTuple{2, Float64} = (1e-8, 1e-5),
        Γ_tol       :: NTuple{2, Float64} = (1e-8, 1e-5),
        χ_tol       :: NTuple{2, Float64} = (1e-8, 1e-5),
        parquet_tol :: NTuple{2, Float64} = (1e-8, 1e-5),
        ODE_tol     :: NTuple{2, Float64} = (1e-8, 1e-2),
        loops       :: Int64              = 1,
        parquet     :: Bool               = false,
        Σ_corr      :: Bool               = true,
        initial     :: Float64            = 50.0,
        final       :: Float64            = 0.05,
        bmin        :: Float64            = 1e-4,
        bmax        :: Float64            = 0.2,
        overwrite   :: Bool               = true,
        cps         :: Vector{Float64}    = Float64[],
        wt          :: Float64            = 23.5,
        ct          :: Float64            = 4.0
        )           :: Nothing

Generate executable Julia file `path` which sets up and runs the FRG solver.
For more information on the different solver parameters see documentation of `launch!`.
"""
function save_launcher!(
    path        :: String,
    f           :: String,
    name        :: String,
    size        :: Int64,
    model       :: String,
    symmetry    :: String,
    J           :: Vector{<:Any}
    ;
    A           :: Float64            = 0.0,
    S           :: Float64            = 0.5,
    β           :: Float64            = 1.0,
    euclidean   :: Bool               = false,
    num_σ       :: Int64              = 25,
    num_Ω       :: Int64              = 15,
    num_ν       :: Int64              = 10,
    num_χ       :: Int64              = 10,
    p_σ         :: NTuple{2, Float64} = (0.3, 1.0),
    p_Ωs        :: NTuple{5, Float64} = (0.3, 0.05, 0.10, 0.1, 50.0),
    p_νs        :: NTuple{5, Float64} = (0.3, 0.05, 0.10, 0.1, 50.0),
    p_Ωt        :: NTuple{5, Float64} = (0.3, 0.15, 0.20, 0.3, 50.0),
    p_νt        :: NTuple{5, Float64} = (0.3, 0.15, 0.20, 0.5, 50.0),
    p_χ         :: NTuple{5, Float64} = (0.3, 0.05, 0.10, 0.1, 50.0),
    lins        :: NTuple{5, Float64} = (5.0, 4.0, 8.0, 6.0, 4.0),
    bounds      :: NTuple{5, Float64} = (0.5, 50.0, 250.0, 150.0, 50.0),
    max_iter    :: Int64              = 10,
    min_eval    :: Int64              = 10,
    max_eval    :: Int64              = 100,
    Σ_tol       :: NTuple{2, Float64} = (1e-8, 1e-5),
    Γ_tol       :: NTuple{2, Float64} = (1e-8, 1e-5),
    χ_tol       :: NTuple{2, Float64} = (1e-8, 1e-5),
    parquet_tol :: NTuple{2, Float64} = (1e-8, 1e-5),
    ODE_tol     :: NTuple{2, Float64} = (1e-8, 1e-2),
    loops       :: Int64              = 1,
    parquet     :: Bool               = false,
    Σ_corr      :: Bool               = true,
    initial     :: Float64            = 50.0,
    final       :: Float64            = 0.05,
    bmin        :: Float64            = 1e-4,
    bmax        :: Float64            = 0.2,
    overwrite   :: Bool               = true,
    cps         :: Vector{Float64}    = Float64[],
    wt          :: Float64            = 23.5,
    ct          :: Float64            = 4.0
    )           :: Nothing

    # convert J for type safety
    J = Array{Array{Float64,1},1}([[x...] for x in J])

    open(path, "w") do file
        # load source code
        write(file, "using PFFRGSolver \n \n")

        # setup for launcher function
        write(file, """launch!("$(f)",
                    "$(name)",
                    $(size),
                    "$(model)",
                    "$(symmetry)",
                    $(J),
                    A           = $(A),
                    S           = $(S),
                    β           = $(β),
                    euclidean   = $(euclidean),
                    num_σ       = $(num_σ),
                    num_Ω       = $(num_Ω),
                    num_ν       = $(num_ν),
                    num_χ       = $(num_χ),
                    p_σ         = $(p_σ),
                    p_Ωs        = $(p_Ωs),
                    p_νs        = $(p_νs),
                    p_Ωt        = $(p_Ωt),
                    p_νt        = $(p_νt),
                    p_χ         = $(p_χ),
                    lins        = $(lins),
                    bounds      = $(bounds),
                    max_iter    = $(max_iter),
                    min_eval    = $(min_eval),
                    max_eval    = $(max_eval),
                    Σ_tol       = $(Σ_tol),
                    Γ_tol       = $(Γ_tol),
                    χ_tol       = $(χ_tol),
                    parquet_tol = $(parquet_tol),
                    ODE_tol     = $(ODE_tol),
                    loops       = $(loops),
                    parquet     = $(parquet),
                    Σ_corr      = $(Σ_corr),
                    initial     = $(initial),
                    final       = $(final),
                    bmin        = $(bmin),
                    bmax        = $(bmax),
                    overwrite   = $(overwrite),
                    cps         = $(cps),
                    wt          = $(wt),
                    ct          = $(ct))""")
    end

    return nothing
end

"""
    make_job!(
        path        :: String,
        dir         :: String,
        input       :: String,
        exe         :: String,
        sbatch_args :: Dict{String, String}
        )           :: Nothing

Generate a SLURM job file `path` to run the FRG solver on a cluster node. `dir` is the job working directory.
`input` is the launcher file generated by `save_launcher!`. `exe` is the path to the Julia executable.
`sbatch_args` is used to set SLURM parameters, e.g. `sbatch_args = Dict(["account" => "my_account", "mem" => "8gb"])`.
"""
function make_job!(
    path        :: String,
    dir         :: String,
    input       :: String,
    exe         :: String,
    sbatch_args :: Dict{String, String}
    )           :: Nothing

    # assert that input is a valid Julia script
    @assert endswith(input, ".jl") "Input must be *.jl file."

    # make local copy to prevent global modification of sbatch_args
    args = copy(sbatch_args)

    # set thread affinity, if not done already
    if haskey(args, "export")
        if occursin("JULIA_EXCLUSIVE", args["export"]) == false
            args["export"] *= ",JULIA_EXCLUSIVE=1"
        end
    else
        push!(args, "export" => "ALL,JULIA_EXCLUSIVE=1")
    end

    # set working directory, if not done already
    if haskey(args, "chdir")
        @warn "Overwriting working directory passed via SBATCH dict ..."
        args["chdir"] = dir
    else
        push!(args, "chdir" => dir)
    end

    # set output file, if not done already
    if haskey(args, "output") == false
        output = split(input, ".jl")[1] * ".out"
        push!(args, "output" => output)
    end

    open(path, "w") do file
        # set SLURM parameters
        write(file, "#!/bin/bash \n")

        for arg in keys(args)
            write(file, "#SBATCH --$(arg)=$(args[arg]) \n")
        end

        write(file, "\n")

        # start calculation
        write(file, "$(exe) -O3 -t \$SLURM_CPUS_PER_TASK $(input)")
    end

    return nothing
end

"""
    make_repository!(
        dir         :: String,
        exe         :: String,
        sbatch_args :: Dict{String, String}
        )           :: Nothing

Generate file structure for several runs of the FRG solver with presumably different parameters.
Assumes that the target folder `dir` contains only launcher files generated by `save_launcher!`.
`exe` is the path to the Julia executable. For each launcher file in `dir` a separate folder and job file are created.
`sbatch_args` is used to set SLURM parameters, e.g. `sbatch_args = Dict(["account" => "my_account", "mem" => "8gb"])`.
"""
function make_repository!(
    dir         :: String,
    exe         :: String,
    sbatch_args :: Dict{String, String}
    )           :: Nothing

    # init folder for saving finished calculations
    fin_dir = joinpath(dir, "finished")

    if isdir(fin_dir) == false
        mkdir(fin_dir)
    end

    # for each *.jl file, init a new folder, move the *.jl file into it and create a job file
    for file in readdir(dir)
        if endswith(file, ".jl")
            # buffer paths
            subdir = joinpath(dir, split(file, ".jl")[1])
            input  = joinpath(subdir, file)
            path   = joinpath(subdir, split(file, ".jl")[1] * ".job")

            # create subdir and job file
            mkdir(subdir)
            mv(joinpath(dir, file), input)
            make_job!(path, subdir, file, exe, sbatch_args)
        end
    end

    return nothing
end

"""
    collect_repository!(
        dir :: String
        )   :: Nothing

Collect final results in file structure generated by `make_repository!`.
Finished calculations are moved to the finished folder.
"""
function collect_repository!(
    dir :: String
    )   :: Nothing

    println("Collecting results from repository ...")

    # check that finished folder exists
    @assert isdir(joinpath(dir, "finished")) "Folder $(joinpath(dir, "finished")) does not exist."

    # for each folder move *_obs, *_cp and *.out, then remove their parent dir. If calculation has not finished set overwrite = false in *.jl file
    for file in readdir(dir)
        if isdir(joinpath(dir, file)) && file != "finished"
            # get file list of subdir
            subdir   = joinpath(dir, file)
            subfiles = readdir(subdir)

            # check if output files exist
            obs_filter = filter(x -> endswith(x, "_obs"), subfiles)

            if length(obs_filter) != 1
                @warn "Could not find unique *_obs file in $(subdir), skipping ..."
                continue
            end

            cp_filter = filter(x -> endswith(x, "_cp"), subfiles)

            if length(cp_filter) != 1
                @warn "Could not find unique *_cp file in $(subdir), skipping ..."
                continue
            end

            out_filter = filter(x -> endswith(x, ".out"), subfiles)

            if length(out_filter) != 1
                @warn "Could not find unique *.out file in $(subdir), skipping ..."
                continue
            end

            # buffer names of output files
            obs_name = obs_filter[1]
            cp_name  = cp_filter[1]
            out_name = out_filter[1]

            # buffer paths of output files
            obs_file = joinpath(subdir, obs_name)
            cp_file  = joinpath(subdir, cp_name)
            out_file = joinpath(subdir, out_name)

            # check if calculation is finished
            obs_data = h5open(obs_file, "r")
            cp_data  = h5open(cp_file, "r")

            if haskey(obs_data, "finished") && haskey(cp_data, "finished")
                # close files
                close(obs_data)
                close(cp_data)

                # move files to finished folder
                mv(obs_file, joinpath(joinpath(dir, "finished"), obs_name))
                mv(cp_file,  joinpath(joinpath(dir, "finished"), cp_name))
                mv(out_file, joinpath(joinpath(dir, "finished"), out_name))

                # remove parent dir
                rm(subdir, recursive = true)
            else
                # close files
                close(obs_data)
                close(cp_data)

                # load the launcher
                launcher_file = joinpath(subdir, file * ".jl")
                stream        = open(launcher_file, "r")
                launcher      = read(stream, String)

                if occursin("overwrite   = true", launcher)
                    # replace overwrite flag and overwrite stream
                    launcher = replace(launcher, "overwrite   = true" => "overwrite   = false")
                    stream   = open(launcher_file, "w")
                    write(stream, launcher)

                    # close the stream
                    close(stream)
                elseif occursin("overwrite   = false", launcher)
                    # close the stream
                    close(stream)
                else
                    # close the stream
                    close(stream)

                    # print error message
                    println("Parameter file $(launcher_file) seems to be broken.")
                end
            end
        end
    end

    println("Done. Results collected in $(joinpath(dir, "finished")).")

    return nothing
end

# load launchers for parquet equations and FRG
include("parquet.jl")
include("launcher_1l.jl")
include("launcher_2l.jl")
include("launcher_ml.jl")

"""
    launch!(
        f           :: String,
        name        :: String,
        size        :: Int64,
        model       :: String,
        symmetry    :: String,
        J           :: Vector{<:Any}
        ;
        A           :: Float64            = 0.0,
        S           :: Float64            = 0.5,
        β           :: Float64            = 1.0,
        euclidean   :: Bool               = false,
        num_σ       :: Int64              = 25,
        num_Ω       :: Int64              = 15,
        num_ν       :: Int64              = 10,
        num_χ       :: Int64              = 10,
        p_σ         :: NTuple{2, Float64} = (0.3, 1.0),
        p_Ωs        :: NTuple{5, Float64} = (0.3, 0.05, 0.10, 0.1, 50.0),
        p_νs        :: NTuple{5, Float64} = (0.3, 0.05, 0.10, 0.1, 50.0),
        p_Ωt        :: NTuple{5, Float64} = (0.3, 0.15, 0.20, 0.3, 50.0),
        p_νt        :: NTuple{5, Float64} = (0.3, 0.15, 0.20, 0.5, 50.0),
        p_χ         :: NTuple{5, Float64} = (0.3, 0.05, 0.10, 0.1, 50.0),
        lins        :: NTuple{5, Float64} = (5.0, 4.0, 8.0, 6.0, 4.0),
        bounds      :: NTuple{5, Float64} = (0.5, 50.0, 250.0, 150.0, 50.0),
        max_iter    :: Int64              = 10,
        min_eval    :: Int64              = 10,
        max_eval    :: Int64              = 100,
        Σ_tol       :: NTuple{2, Float64} = (1e-8, 1e-5),
        Γ_tol       :: NTuple{2, Float64} = (1e-8, 1e-5),
        χ_tol       :: NTuple{2, Float64} = (1e-8, 1e-5),
        parquet_tol :: NTuple{2, Float64} = (1e-8, 1e-5),
        ODE_tol     :: NTuple{2, Float64} = (1e-8, 1e-2),
        loops       :: Int64              = 1,
        parquet     :: Bool               = false,
        Σ_corr      :: Bool               = true,
        initial     :: Float64            = 50.0,
        final       :: Float64            = 0.05,
        bmin        :: Float64            = 1e-4,
        bmax        :: Float64            = 0.2,
        overwrite   :: Bool               = true,
        cps         :: Vector{Float64}    = Float64[],
        wt          :: Float64            = 23.5,
        ct          :: Float64            = 4.0
        )           :: Nothing

Runs the FRG solver. A detailed explanation of the solver parameters is given below:
* `f`           : name of the output files. The solver will generate two files (`f * "_obs"` and `f * "_cp"`) containing observables and checkpoints respectively.
* `name`        : name of the lattice
* `size`        : size of the lattice. Correlations are truncated beyond this range.
* `model`       : name of the spin model. Defines coupling structure.
* `symmetry`    : symmetry of the spin model. Used to reduce computational complexity.
* `J`           : coupling vector of the spin model. J is normalized together with A during initialization of the solver.
* `A`           : on-site repulsion term. A is normalized together with J during initialization of the solver.
* `S`           : total spin quantum number (only relevant for pure Heisenberg models)
* `β`           : damping factor for fixed point iterations of parquet equations
* `euclidean`   : flag to build lattice by Euclidean (aka real space) instead of bond distance
* `num_σ`       : number of non-zero, positive frequencies for the self energy
* `num_Ω`       : number of non-zero, positive frequencies for the bosonic axis of the two-particle irreducible channels
* `num_ν`       : number of non-zero, positive frequencies for the fermionic axis of the two-particle irreducible channels
* `num_χ`       : number of non-zero, positive frequencies for the correlations
* `p_σ`         : parameters for updating self energy mesh between ODE steps \n
                  p_σ[1] gives the percentage of linearly spaced frequencies
                  p_σ[2] sets the linear extent relative to the position of the quasi-particle peak
* `p_Ω / p_ν`   : parameters for updating bosonic / fermionic s (u) and t channel frequency meshes between ODE steps \n
                  p_Γ[1] gives the percentage of linearly spaced frequencies
                  p_Γ[2] (p_Γ[3]) sets the lower (upper) bound for the accepted relative deviation between the values at the origin and the first finite frequency
                  p_Γ[4] sets the lower bound for the linear spacing in units of the cutoff Λ
                  p_Γ[5] sets the upper bound for the linear extent in units of the cutoff Λ
* `p_χ`         : parameters for updating correlation mesh between ODE steps \n
                  p_χ[1] gives the percentage of linearly spaced frequencies
                  p_χ[2] (p_χ[3]) sets the lower (upper) bound for the accepted relative deviation between the values at the origin and the first finite frequency
                  p_χ[4] sets the lower bound for the linear spacing in units of the cutoff Λ
                  p_χ[5] sets the upper bound for the linear extent in units of the cutoff Λ
* `lins`        : parameters for controlling the scaling of frequency meshes before adaptive scanning is utilized \n
                  lins[1] gives the scale, in units of |J|, beyond which adaptive meshes are used
                  lins[2] gives the linear extent, in units of the cutoff Λ, for the self energy
                  lins[3] gives the linear extent, in units of the cutoff Λ, for the bosonic axis of the two-particle irreducible channels
                  lins[4] gives the linear extent, in units of the cutoff Λ, for the fermionic axis of the two-particle irreducible channels
                  lins[5] gives the linear extent, in units of the cutoff Λ, for the correlations
* `bounds`      : parameters for controlling the upper mesh bounds \n
                  bounds[1] gives, in units of |J|, the stopping scale beyond which no further contraction of the meshes is performed
                  bounds[2] gives, in units of the cutoff Λ, the upper bound for the self energy
                  bounds[3] gives, in units of the cutoff Λ, the upper bound for the bosonic axis of the two-particle irreducible channels
                  bounds[4] gives, in units of the cutoff Λ, the upper bound for the fermionic axis of the two-particle irreducible channels
                  bounds[5] gives, in units of the cutoff Λ, the upper bound for the correlations
* `max_iter`    : maximum number of parquet iterations
* `min_eval`    : minimum initial number of subdivisions for vertex quadrature. eval is min_eval for parquet iterations.
* `max_eval`    : maximum initial number of subdivisions for vertex quadrature. eval is ramped up from min_eval to max_eval as a function of the cutoff Λ (for Λ < |J|).
* `Σ_tol`       : absolute and relative error tolerance for self energy quadrature
* `Γ_tol`       : absolute and relative error tolerance for vertex quadrature
* `χ_tol`       : absolute and relative error tolerance for correlation quadrature
* `parquet_tol` : absolute and relative error tolerance for convergence of parquet iterations
* `ODE_tol`     : absolute and relative error tolerance for Bogacki-Shampine solver
* `loops`       : number of loops to be calculated
* `parquet`     : flag to enable parquet iterations. If `false`, initial condition is chosen as bare vertex.
* `Σ_corr`      : flag to enable self energy corrections. Has no effect for 'loops <= 2'.
* `initial`     : start value of the cutoff in units of |J|
* `final`       : final value of the cutoff in units of |J|. If `final = initial` and `parquet = true` only a solution of the parquet equations is computed.
* `bmin`        : minimum step size of the ODE solver in units of |J|
* `bmax`        : maximum step size of the ODE solver in units of Λ
* `overwrite`   : flag to indicate whether a new calculation should be started. If false, checks if `f * "_obs"` and `f * "_cp"` exist and continues calculation from available checkpoint with lowest cutoff.
* `cps`         : list of intermediate cutoffs in units of |J|, where a checkpoint with full vertex data shall be generated
* `wt`          : wall time (in hours) for the calculation. Should be set according to cluster configurations. If run remote, set `wt = Inf` to avoid data loss. \n
                  WARNING: For run times longer than wt, no checkpoints are created.
* `ct`          : minimum time (in hours) between subsequent checkpoints
"""
function launch!(
    f           :: String,
    name        :: String,
    size        :: Int64,
    model       :: String,
    symmetry    :: String,
    J           :: Vector{<:Any}
    ;
    A           :: Float64            = 0.0,
    S           :: Float64            = 0.5,
    β           :: Float64            = 1.0,
    euclidean   :: Bool               = false,
    num_σ       :: Int64              = 25,
    num_Ω       :: Int64              = 15,
    num_ν       :: Int64              = 10,
    num_χ       :: Int64              = 10,
    p_σ         :: NTuple{2, Float64} = (0.3, 1.0),
    p_Ωs        :: NTuple{5, Float64} = (0.3, 0.05, 0.10, 0.1, 50.0),
    p_νs        :: NTuple{5, Float64} = (0.3, 0.05, 0.10, 0.1, 50.0),
    p_Ωt        :: NTuple{5, Float64} = (0.3, 0.15, 0.20, 0.3, 50.0),
    p_νt        :: NTuple{5, Float64} = (0.3, 0.15, 0.20, 0.5, 50.0),
    p_χ         :: NTuple{5, Float64} = (0.3, 0.05, 0.10, 0.1, 50.0),
    lins        :: NTuple{5, Float64} = (5.0, 4.0, 8.0, 6.0, 4.0),
    bounds      :: NTuple{5, Float64} = (0.5, 50.0, 250.0, 150.0, 50.0),
    max_iter    :: Int64              = 10,
    min_eval    :: Int64              = 10,
    max_eval    :: Int64              = 100,
    Σ_tol       :: NTuple{2, Float64} = (1e-8, 1e-5),
    Γ_tol       :: NTuple{2, Float64} = (1e-8, 1e-5),
    χ_tol       :: NTuple{2, Float64} = (1e-8, 1e-5),
    parquet_tol :: NTuple{2, Float64} = (1e-8, 1e-5),
    ODE_tol     :: NTuple{2, Float64} = (1e-8, 1e-2),
    loops       :: Int64              = 1,
    parquet     :: Bool               = false,
    Σ_corr      :: Bool               = true,
    initial     :: Float64            = 50.0,
    final       :: Float64            = 0.05,
    bmin        :: Float64            = 1e-4,
    bmax        :: Float64            = 0.2,
    overwrite   :: Bool               = true,
    cps         :: Vector{Float64}    = Float64[],
    wt          :: Float64            = 23.5,
    ct          :: Float64            = 4.0
    )           :: Nothing

    # init timers for checkpointing
    t  = Dates.now()
    t0 = Dates.now()

    println()
    println("################################################################################")
    println("Initializing solver ...")
    println(); println()

    # check if symmetry parameter is valid
    symmetries = String["su2", "u1-dm"]
    @assert in(symmetry, symmetries) "Symmetry $(symmetry) unknown. Valid arguments are su2 and u1-dm."

    # init names for observables and checkpoints file
    obs_file = f * "_obs"
    cp_file  = f * "_cp"

    # test if a new calculation should be started
    if overwrite
        println("overwrite = true, starting from scratch ...")

        # delete existing observables
        if isfile(obs_file)
            rm(obs_file)
        end

        # delete existing checkpoints
        if isfile(cp_file)
            rm(cp_file)
        end

        # open new files
        obs = h5open(obs_file, "cw")
        cp  = h5open(cp_file, "cw")

        # convert J for type safety
        J = Vector{Vector{Float64}}([[x...] for x in J])

        # normalize couplings and level repulsion
        JAtemp = normalize([J, [[A]]])
        J      = JAtemp[1]
        A      = JAtemp[2][1][1]

        # build lattice and save to files
        println();
        l = get_lattice(name, size, euclidean = euclidean)

        println();
        r = get_reduced_lattice(model, J, l)

        save!(obs, r)
        save!(cp, r)

        # close files
        close(obs)
        close(cp)

        # build meshes
        σ  = get_mesh(lins[2] * initial, bounds[2] * max(initial, bounds[1]), num_σ, p_σ[1])
        Ωs = get_mesh(lins[3] * initial, bounds[3] * max(initial, bounds[1]), num_Ω, p_Ωs[1])
        νs = get_mesh(lins[4] * initial, bounds[4] * max(initial, bounds[1]), num_ν, p_νs[1])
        Ωt = get_mesh(lins[3] * initial, bounds[3] * max(initial, bounds[1]), num_Ω, p_Ωt[1])
        νt = get_mesh(lins[4] * initial, bounds[4] * max(initial, bounds[1]), num_ν, p_νt[1])
        χ  = get_mesh(lins[5] * initial, bounds[5] * max(initial, bounds[1]), num_χ, p_χ[1])
        m  = Mesh(num_σ + 1, num_Ω + 1, num_ν + 1, num_χ + 1, σ, Ωs, νs, Ωt, νt, χ)

        # build action
        a = get_action_empty(symmetry, r, m, S = S)
        init_action!(l, r, a)
        set_repulsion!(A, a)

        # initialize by parquet iterations
        if parquet
            println(); println()
            println("Warming up with some parquet iterations ...")
            flush(stdout)
            launch_parquet!(obs_file, cp_file, symmetry, l, r, m, a, initial, bmax * initial, β, max_iter, min_eval, Σ_tol, Γ_tol, χ_tol, parquet_tol, S = S)
            println("Done. Action is initialized with parquet solution.")
        end

        println(); println()
        println("Solver is ready.")
        println("################################################################################")
        println()

        # start calculation
        println("Renormalization group flow with ℓ = $(loops) ...")
        flush(stdout)

        if loops == 1
            launch_1l!(obs_file, cp_file, symmetry, l, r, m, a, p_σ, p_Ωs, p_νs, p_Ωt, p_νt, p_χ, lins, bounds, initial, final, bmax * initial, bmin, bmax, min_eval, max_eval, Σ_tol, Γ_tol, χ_tol, ODE_tol, t, t0, cps, wt, ct, A = A, S = S)
        elseif loops == 2
            launch_2l!(obs_file, cp_file, symmetry, l, r, m, a, p_σ, p_Ωs, p_νs, p_Ωt, p_νt, p_χ, lins, bounds, initial, final, bmax * initial, bmin, bmax, min_eval, max_eval, Σ_tol, Γ_tol, χ_tol, ODE_tol, t, t0, cps, wt, ct, A = A, S = S)
        elseif loops >= 3
            launch_ml!(obs_file, cp_file, symmetry, l, r, m, a, p_σ, p_Ωs, p_νs, p_Ωt, p_νt, p_χ, lins, bounds, loops, Σ_corr, initial, final, bmax * initial, bmin, bmax, min_eval, max_eval, Σ_tol, Γ_tol, χ_tol, ODE_tol, t, t0, cps, wt, ct, A = A, S = S)
        end
    else
        println("overwrite = false, trying to load data ...")

        if isfile(obs_file) && isfile(cp_file)
            println()
            println("Found existing output files, checking status ...")

            # open files
            obs = h5open(obs_file, "cw")
            cp  = h5open(cp_file, "cw")

            if haskey(obs, "finished") && haskey(cp, "finished")
                # close files
                close(obs)
                close(cp)

                println(); println()
                println("Calculation has finished already.")
                println("################################################################################")
                flush(stdout)
            else
                println()
                println("Final Λ has not been reached, resuming calculation ...")

                # load data
                l, r              = read_lattice(cp)
                Λ, dΛ, m, a_inter = read_checkpoint(cp, 0.0)
                χ                 = read_χ_all(obs, Λ; verbose = false)[2]

                # update frequency mesh
                a = deepcopy(a_inter)
                m = resample_from_to(Λ, p_σ, p_Ωs, p_νs, p_Ωt, p_νt, p_χ, lins, bounds, m, a_inter, a, χ)

                # close files
                close(obs)
                close(cp)

                println(); println()
                println("Solver is ready.")
                println("################################################################################")
                println()

                # resume calculation
                println("Renormalization group flow with ℓ = $(loops) ...")
                flush(stdout)

                if loops == 1
                    launch_1l!(obs_file, cp_file, symmetry, l, r, m, a, p_σ, p_Ωs, p_νs, p_Ωt, p_νt, p_χ, lins, bounds, Λ, final, dΛ, bmin, bmax, min_eval, max_eval, Σ_tol, Γ_tol, χ_tol, ODE_tol, t, t0, cps, wt, ct, A = A, S = S)
                elseif loops == 2
                    launch_2l!(obs_file, cp_file, symmetry, l, r, m, a, p_σ, p_Ωs, p_νs, p_Ωt, p_νt, p_χ, lins, bounds, Λ, final, dΛ, bmin, bmax, min_eval, max_eval, Σ_tol, Γ_tol, χ_tol, ODE_tol, t, t0, cps, wt, ct, A = A, S = S)
                elseif loops >= 3
                    launch_ml!(obs_file, cp_file, symmetry, l, r, m, a, p_σ, p_Ωs, p_νs, p_Ωt, p_νt, p_χ, lins, bounds, loops, Σ_corr, Λ, final, dΛ, bmin, bmax, min_eval, max_eval, Σ_tol, Γ_tol, χ_tol, ODE_tol, t, t0, cps, wt, ct, A = A, S = S)
                end
            end
        else
            println(); println()
            println("Found no existing output files, terminating solver ...")
            println("################################################################################")
            flush(stdout)
        end
    end

    println()
    println("################################################################################")
    println("Solver terminated.")
    println("################################################################################")
    flush(stdout)

    return nothing
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   This files must call all the tests related to the submodules.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "SGP4" begin
    cd("./SGP4")
    include("./SGP4/sgp4.jl")
    cd("../")
end

@testset "SatelliteToolboxTLE" begin
    cd("./SatelliteToolboxTLE")
    include("./SatelliteToolboxTLE/tle.jl")
    cd("../")
end


text:
function ee_position!(position, robot, q::AbstractArray{T}) where {T}
    state = robot.statecache[T]
    set_configuration!(state, q)

    point = Point3D(robot.frame_ee, 0.0, 0.0, 0.0)
    position[:] = transform(state, point, root_frame(robot.mechanism)).v
end

function cb_eval_fc_con_ee(kc, cb, evalRequest, evalResult, userParams)
    x = evalRequest.x
    problem, robot = userParams

    # dimension of each mesh point
    nₓ = robot.n_q + robot.n_v + robot.n_τ

    length_c = 3

    for (i, k) = enumerate(sort(collect(keys(problem.ee_pos))))
        ind_qᵢ = range(1 + (k - 1) * nₓ, length=robot.n_q)  # indices of the decision variables
        ind_c = (1:length_c) .+ ((i - 1) * length_c)
        @views ee_position!(evalResult.c[ind_c], robot, x[ind_qᵢ])  # constraint evaluation at that point
    end

    return 0
end

function cb_eval_ga_con_ee(kc, cb, evalRequest, evalResult, userParams)
    x = evalRequest.x
    problem, robot = userParams

    # dimension of each mesh point
    nₓ = robot.n_q + robot.n_v + robot.n_τ

    for (i, k) = enumerate(sort(collect(keys(problem.ee_pos))))
        ind_qᵢ = range(1 + (k - 1) * nₓ, length=robot.n_q)  # indices of the decision variables
        problem.jacdata_ee_position(x[ind_qᵢ])              # jacobian evaluation at that point
        ind_jac = (1:problem.jacdata_ee_position.length_jac) .+ ((i - 1) * problem.jacdata_ee_position.length_jac)
        evalResult.jac[ind_jac] = nonzeros(problem.jacdata_ee_position.jac)  # pass jacobian to Knitro
    end

    return 0
end

function get_ee_path(problem::Problem, robot::Robot, x)
    positions = zeros(3, problem.num_knots)

    # dimension of each mesh point
    nₓ = robot.n_q + robot.n_v + robot.n_τ

    length_c = 3

    for i = (1:problem.num_knots) .- 1
        ind_c = (1:length_c) .+ (i * length_c)
        ind_qᵢ = range(1 + i * nₓ, length=robot.n_q)
        @views ee_position!(positions[ind_c], robot, x[ind_qᵢ])
    end

    return positions
end

function show_ee_path(vis::Visualizer, ee_positions::Array{Float64})
    points = collect(eachcol(ee_positions))
    material = LineBasicMaterial(color=colorant"yellow", linewidth=2.0)
    setobject!(vis["ee path"], Object(PointCloud(points), material, "Line"))
    nothing
end


text:
# Functions for computing MFM partition distribution coefficients V_n(t) and also p(k|t).
module MFM

using SpecialFunctions
lgamma_(x) = logabsgamma(x)[1]

logsumexp(a,b) = (m = max(a,b); m == -Inf ? -Inf : log(exp(a-m) + exp(b-m)) + m)

# Compute log_v[t] = log(V_n(t)) under the given MFM parameters, for t=1:upto.
function coefficients(log_pk,gamma,n,upto)
    tolerance = 1e-12
    log_v = zeros(upto)
    for t = 1:upto
        if t>n; log_v[t] = -Inf; continue; end
        a,c,k,p = 0.0, -Inf, 1, 0.0
        while abs(a-c) > tolerance || p < 1.0-tolerance  # Note: The first condition is false when a = c = -Inf
            if k >= t
                a = c
                b = lgamma_(k+1)-lgamma_(k-t+1) - lgamma_(k*gamma+n)+lgamma_(k*gamma) + log_pk(k)
                c = logsumexp(a,b)
            end
            p += exp(log_pk(k))
            k = k+1
        end
        log_v[t] = c
    end
    return log_v
end

# Compute p[k,t]=p(k|t) under the given MFM parameters, for k=1:upto_k, t=1:upto_t.
function p_kt(log_pk,gamma,n,upto_k,upto_t)
    @assert(upto_t<=n,"p(k|t) is undefined for t>n.")
    log_v = coefficients(log_pk,gamma,n,upto_t)
    p = zeros(upto_k,upto_t)
    for k = 1:upto_k, t = 1:min(k,upto_t)
        b = lgamma_(k+1)-lgamma_(k-t+1) - lgamma_(k*gamma+n)+lgamma_(k*gamma) + log_pk(k)
        p[k,t] = exp(b - log_v[t])
    end
    return p
end

end # module MFM


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((3, true, false),),
    ((3, false, true),),
    ((3, false, true), :SOCExpPSD),
    ]
insts["fast"] = [
    ((10, true, false),),
    ((10, false, true),),
    ((10, false, true), :SOCExpPSD),
    ((50, true, false),),
    ((50, false, true),),
    ((50, false, true), :SOCExpPSD),
    ((400, true, false),),
    ((400, false, true),),
    ((400, true, false),),
    ((400, false, true),),
    ((400, false, true), :SOCExpPSD),
    ]
insts["various"] = [
    ((50, false, true), :SOCExpPSD),
    ((1000, true, false),),
    ((1000, false, true),),
    ((2000, true, false),),
    ((2000, false, true),),
    ((4000, true, false),),
    ((4000, false, true),),
    ((8000, true, false),),
    ((8000, false, true),),
    ]
return (PortfolioJuMP, insts)


text:
push!(LOAD_PATH,"app/resources/insurancecontracts")
import InsuranceContractsController, StippleController
using Genie, Genie.Router, Genie.Requests, Genie.Renderer.Html, Genie.Renderer.Json, Stipple, StippleController

route("/") do
  serve_static_file("welcome.html")
end

route("/contracts", InsuranceContractsController.insurancecontracts_view)

route("/history") do
  history_id = params(:history_id, "0")
  println("history_id= " * history_id)
  InsuranceContractsController.historyforest_view(parse(Int,history_id))
end

route("/csection") do
  history_id = params(:history_id, "0")
  version_id = params(:version_id, "0")
  println("history_id= " * history_id * " version_id= " * version_id)
  InsuranceContractsController.csection_view(parse(Int,history_id),parse(Int,version_id))
end

route("/jsonpayload", method = POST) do
  if haskey(jsonpayload(), "command")
      if jsonpayload()["command"] == "showsection"
          JSON.json(InsuranceContractsController.csection(3,3))
      else
          JSON.json(Dict{String,Any}("error" => "only showsection recognizable"))
      end
  else
      JSON.json(jsonpayload())
  end
end

initModel(3,3)

route("/reactive") do
  StippleController.render("initial")
end

route("/sub", method = POST) do
  StippleController.process("sub")
  redirect("/reactive")
end

Genie.up()
Genie.config.log_requests = false

max_stars_repo_name: michaelfliegner/BitempWebApp

max_issues_repo_name: michaelfliegner/BitempWebApp


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -1.3634297333714107 <= q <= -0.20981772159177886)

                     add_NL_constraint(m, :(1/(1 + exp(-(-0.25086622252953283 + -0.4158912166699471*1/(1 + exp(-(0.18862560740304835 + -0.15744636994779526*1/(1 + exp(-(0.5575683410422836 + 0.3104874388021095*1/(1 + exp(-(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4])))) + 0.12330243724213741*1/(1 + exp(-(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))))))) + 0.6894873557418655*1/(1 + exp(-(0.5818927234969289 + 0.4060546195732435*1/(1 + exp(-(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4])))) + 0.22410021012329473*1/(1 + exp(-(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4])))))))))) + 0.15556574286983027*1/(1 + exp(-(-0.6464546506721769 + -0.3052437329487665*1/(1 + exp(-(0.5575683410422836 + 0.3104874388021095*1/(1 + exp(-(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4])))) + 0.12330243724213741*1/(1 + exp(-(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))))))) + 0.4287563119683666*1/(1 + exp(-(0.5818927234969289 + 0.4060546195732435*1/(1 + exp(-(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4])))) + 0.22410021012329473*1/(1 + exp(-(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))))))))))))) + 1/(1 + exp(-(-0.9805026843034819 + -0.19393497198753096*1/(1 + exp(-(0.18862560740304835 + -0.15744636994779526*1/(1 + exp(-(0.5575683410422836 + 0.3104874388021095*1/(1 + exp(-(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4])))) + 0.12330243724213741*1/(1 + exp(-(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))))))) + 0.6894873557418655*1/(1 + exp(-(0.5818927234969289 + 0.4060546195732435*1/(1 + exp(-(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4])))) + 0.22410021012329473*1/(1 + exp(-(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4])))))))))) + -0.01642706539433636*1/(1 + exp(-(-0.6464546506721769 + -0.3052437329487665*1/(1 + exp(-(0.5575683410422836 + 0.3104874388021095*1/(1 + exp(-(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4])))) + 0.12330243724213741*1/(1 + exp(-(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))))))) + 0.4287563119683666*1/(1 + exp(-(0.5818927234969289 + 0.4060546195732435*1/(1 + exp(-(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4])))) + 0.22410021012329473*1/(1 + exp(-(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))))))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
"""
Define constant variables
"""

const USE_INVB = true

const INF = 1.0e+20
const EPS = 1e-6

const MAX_ITER = 30000000000
const MAX_HISTORY = 100
const PRINT_ITER_FREQ = 10

@enum(
    Status,
    NotSolved,
    Solve,
    Optimal,
    Unbounded,
    Infeasible,
    Feasible,
    Cycle,
)

const BASIS_BASIC = 1
const BASIS_AT_LOWER = 2
const BASIS_AT_UPPER = 3
const BASIS_FREE = 4

@enum(
    PivotRule,
    Bland,
    Steepest,
    Dantzig,
    Artificial,
)


text:
using Stheno, Optim, Plots, Zygote

using Stheno: derivative

const σ²_noise = 1e-3

# Helper function: Parameter handling.
function unpack(θ)
    σ² = exp(θ[1]) + 1e-3
    l = exp(θ[2]) + 1e-3
    return σ², l
end

# Helper function: construct the model and it's derivative.
function build_model(θ)
    σ², l = unpack(θ)
    k = σ² * stretch(EQ(), 1 / l)
    f = GP(k, GPC())
    df = derivative(f)
    return f, df
end

"""
    nlml(θ, x_f, y_f, x_df, y_df)

Compute the negative log marginal likelihood of the parameters θ given observations y_f of
f at x_f and observations y_df of df at x_df.
"""
function nlml(
    θ,
    x_f::AbstractVector{<:Real},
    y_f::AbstractVector{<:Real},
    x_df::AbstractVector{<:Real},
    y_df::AbstractVector{<:Real},
)
    f, df = build_model(θ)

    # Specify points in processes at which we'll make observations.
    fx = f(x_f, σ²_noise)
    dfx = df(x_df, σ²_noise)
    return -logpdf([fx, dfx], [y_f, y_df])
end

# Toy problem with synthetic data.
g = sin
dg = cos

x_g = range(-10.0, 10.0; length=15);
x_dg = range(-10.0, 10.0; length=20);

y_g = g.(x_g) .+ sqrt(σ²_noise) .* randn(length(x_g));
y_dg = dg.(x_dg) .+ sqrt(σ²_noise) .* randn(length(x_dg));

nlml(θ) = nlml(θ, x_g, y_g, x_dg, y_dg)

# Optimise the hyper-parameters using a gradient-free method to start with.
θ0 = randn(2);
results = Optim.optimize(nlml, θ -> only(Zygote.gradient(nlml, θ)), θ0, BFGS(); inplace=false)
θ_opt  = results.minimizer;

# Get process at optimal parameters, and produce the posteriors.
f, df = build_model(θ_opt);
f_post, df_post = (f, df) | (f(x_g, 1e-3) ← y_g, df(x_dg, 1e-3) ← y_dg);

# Make posterior predictions and visualise.
x_pr = range(-15.0, 15.0; length=250);

plt = plot();
scatter!(plt, x_g, y_g; color=:red, label="y_g");
scatter!(plt, x_dg, y_dg; color=:blue, label="y_dg");
plot!(plt, x_pr, g.(x_pr); color=:red, label="g");
plot!(plt, x_pr, dg.(x_pr); color=:blue, label="dg");
plot!(plt, f_post(x_pr, 1e-3); samples=10, label="f_post");
plot!(plt, df_post(x_pr, 1e-3); samples=10, label="df_post");
display(plt);


text:
#
# Mimics the setup in demo.sh for Stanford's C implementation.
# Roughly ~3x slower than the C model as of 09/28/15
#
# The dataset is from char-rnn.
#
# To Download go to
using Glove
using DataStructures

S = Float64
println("Using type: ", S)

file = "input.txt"
tokens = split(read(file, String))
println("Counted ", length(tokens), " tokens")
ct = counter(tokens)
println("Counted ", length(ct.map), " unique tokens")
filter!((e) -> e.second >= 5, ct.map)
println("Counted ", length(ct.map), " filtered tokens")

filtered_tokens = collect(keys(ct))
table = LookupTable(filtered_tokens)

println("Creating Co-occurences")

println("before jit...")
@time codict = weightedsums(S, table, tokens, window=15)
println("after jit ...")
@time codict = weightedsums(S, table, tokens, window=15)

# make the covec
n = length(codict)
covec = CooccurenceVector{S}(undef, n)
cokeys = collect(keys(codict))
covals = collect(values(codict))
@inbounds for i=1:n
  covec[i] = (cokeys[i][1], cokeys[i][2], covals[i])
end

vocabsize = length(table)
vecsize = 50
model = Model(S, vocabsize, vecsize)

println("before jit ...")
@time adagrad!(S, model, covec, 1, xmax=10)
println("after jit ...")
@time adagrad!(S, model, covec, 1, xmax=10)


text:
using Gadfly
srand(111)
plot(x=rand(7), y=rand(7))

using Gadfly
srand(111)
plot(x=rand(113), Geom.histogram(bincount=10))

using Winston

srand(111)

pl = Winston.plot(cumsum(rand(100) .-0.5), "g", cumsum(rand(100) .-0.5), "b")

display(pl)

max_stars_repo_path: Section 4/Julia+Visualization.jl

max_issues_repo_path: Section 4/Julia+Visualization.jl

max_forks_repo_path: Section 4/Julia+Visualization.jl


text:
function initial_conditions_steady_state(rdpp::ReactionDiffusionPricePaths)
    φ = [(xᵢ - rdpp.p₀) * rdpp.source_term.λ /
         (2 * rdpp.D * rdpp.source_term.μ) *
         exp((-rdpp.source_term.μ * rdpp.L^2) / 4) +
         (sqrt(pi) * rdpp.source_term.λ *
          erf(sqrt(rdpp.source_term.μ) * (rdpp.p₀ - xᵢ))) /
         (4 * rdpp.D * rdpp.source_term.μ^(3 / 2)) for xᵢ in rdpp.x]
    return φ
end

function sample_mid_price_path(rdpp, Δt, price_path)
    mid_prices = zeros(Float64, rdpp.T + 1)
    mid_prices[1] = price_path[1]
    for t = 1:rdpp.T
        close_ind = floor(Int, t / Δt) + 1
        mid_prices[t+1] = price_path[close_ind]
    end
    return mid_prices
end

function extract_mid_price(rdpp, lob_density)
    mid_price_ind = 2
    while (lob_density[mid_price_ind] > 0) | (lob_density[mid_price_ind+1]>lob_density[mid_price_ind])
        mid_price_ind += 1
        if mid_price_ind > rdpp.M
            mid_price_ind += 1
            break
        end
    end
    y1 = lob_density[mid_price_ind-1]
    y2 = lob_density[mid_price_ind]
    x1 = rdpp.x[mid_price_ind-1]
    mid_price = round(-(y1 * rdpp.Δx)/(y2 - y1) + x1, digits = 2)
    return mid_price
end

function dtrw_solver(rdpp::ReactionDiffusionPricePaths)
    Δx = rdpp.L / rdpp.M
    Δt = (Δx^2) / (2.0 * rdpp.D)
    time_steps = floor(Int, rdpp.T / rdpp.Δt)

    φ = ones(Float64, rdpp.M + 1, time_steps + 1)
    φ[:, 1] = initial_conditions_steady_state(rdpp)

    p = ones(Float64, time_steps + 1)
    p[1] = rdpp.p₀
    ϵ = rand(Normal(0.0, 1.0), time_steps)
    P⁺s = ones(Float64, time_steps)
    P⁻s = ones(Float64, time_steps)

    @inbounds for n = 1:time_steps
        Vₜ = sign(ϵ[n]) * min(abs(rdpp.σ * ϵ[n]), rdpp.Δx / rdpp.Δt)
        V = (-Vₜ .* rdpp.x) ./ (2.0 * rdpp.D)

        P⁺ = 1 / (1 + exp(-rdpp.β * Vₜ * rdpp.Δx / rdpp.D))
        P⁻ = 1 - P⁺

        P⁺s[n] = P⁺
        P⁻s[n] = P⁻

        φ₋₁ = (1-Vₜ*rdpp.Δx/rdpp.D) * φ[1, n]
        φₘ₊₁ = (1+Vₜ*rdpp.Δx/rdpp.D) * φ[end, n]

        φ[1, n+1] = P⁺ * φ₋₁ + P⁻ * φ[2, n] + rdpp.source_term(rdpp.x[1], p[n])

        φ[end, n+1] = P⁻ * φₘ₊₁ + P⁺ * φ[end-1, n] +
                      rdpp.source_term(rdpp.x[end], p[n])

        φ[2:end-1, n+1] = P⁺ * φ[1:end-2, n] + P⁻ * φ[3:end, n] +
                          [rdpp.source_term(xᵢ, p[n]) for xᵢ in rdpp.x[2:end-1]]

        p[n+1] = extract_mid_price(rdpp, φ[:, n+1])
    end
    mid_price_bars_close = sample_mid_price_path(rdpp, rdpp.Δt, p)
    return φ, p, mid_price_bars_close, P⁺s, P⁻s
end


text:
function stationary_algebraic_simple(parameters, settings)
    @unpack μ, υ, θ, r, ζ = parameters
    @unpack T = settings
    g = μ + (1 - (θ - 1)*ζ*(r(T) - μ))/((θ - 1)^2 * ζ) + υ^2/2 * (θ*(θ*(θ - 1)*(r(T) - μ - υ^2/2)*ζ - 2) + 1)/((θ - 1)*((θ - 1)*(r(T) - μ - υ^2/2)*ζ - 1)); # (9)
    ν = (μ-g)/υ^2 + √(((g-μ)/υ^2)^2 + (r(T) - g)/(υ^2/2)); # (11)
    v(z) = 1/(r(T) - μ - υ^2/2) * (1 + 1/ν * exp(-(ν + 1)*z)); # (10). Note that these values are in the transformed space.
    return (g = g, ν = ν, v = v)
end

# Minimizes value-matching condition (31)
function stationary_numerical_simple(parameters, settings)
    @unpack z_ex, T = settings
    x0 = settings.stationary_x0(parameters, settings)
    @unpack μ, υ, θ, r, ζ, π = parameters
    @assert z_ex[1] == 0.0 && issorted(z_ex)  # validate grid
    z = z_ex[2:end-1]
    r_T = r(T)

    ω = ω_weights(z_ex, θ, 1) # (26), see utils/quadrature.jl

    bc = (Mixed(ξ = 1), Mixed(ξ = 1)) # boundary conditions for differential operators
    L_1_minus = L₁₋bc(z_ex, bc) # use backward difference as the drift is negative
    L_2 = L₂bc(z_ex, bc)

    Ξ₁ = 1/(1 - 1*(z[1] - 0.0)) # (24)

    function stationary_numerical_given_g(in)
        g = in[1]
        A = (r_T - μ - υ^2/2)*I - (μ + υ^2 - g)*L_1_minus - υ^2/2 * L_2 # (19)
        v = A \ π.(T, z) # discretized system of ODE for v, where v'(T) = 0 (30)
        return Ξ₁*v[1] - dot(ω, v) + ζ # value matching condition (31)
    end

    sol = nlsolve(stationary_numerical_given_g, x0, inplace = false)
    g_T = sol.zero[1]
    @assert converged(sol) "The solver didn't converge; exiting now."
    @assert (μ + υ^2/2 - g_T < 0) # Negative drift condition (20)
    @assert r_T > g_T "r_T > g_T is required for convergence" # From (19)

    # Use the g_T to recreate L_T and v_T
    A_T = (r_T - μ - υ^2/2)*I - (μ + υ^2 - g_T)*L_1_minus - υ^2/2 * L_2 # (19)
    v_T = A_T \ π.(T, z) # (30)
    return (g = g_T, v = v_T, residual = stationary_numerical_given_g(g_T))
end


text:
register(DataDep(
    "pudl-raw-ferc1",
    """
    Dataset: PUDL Raw FERC Form 1
    Website: https://doi.org/10.5281/zenodo.4127044
    Author: Gosnell, Christina, Selvans, Zane
    Date of Publication: October 24, 2020
    License:

    Raw Federal Energy Regulatory Commission (FERC) Form 1 data, archived from https://www.ferc.gov/industries-data/electric/general-information/electric-industry-forms/form-1-electric-utility-annual

    Form No. 1 financial and operating report from the Federal Energy Regulatory Commission (FERC).

    $PUDL_DESCRIPTION
    """,
    [
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/datapackage.json",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-1994.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-1995.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-1996.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-1997.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-1998.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-1999.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2000.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2001.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2002.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2003.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2004.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2005.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2006.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2007.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2008.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2009.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2010.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2011.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2012.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2013.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2014.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2015.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2016.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2017.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2018.zip",
        "https://zenodo.org/api/files/c991f4cc-fdc7-43f7-807c-61428fa1b369/ferc1-2019.zip",
    ],
    "98f2b2a90ea917b11504015e141c1e172c5673de580009a8a8cf494dede479a3",
))


text:
using Test

m = page_model()
include("../src/components/Discontinuity.jl")

add_comp!(m, Discontinuity)

set_param!(m, :Discontinuity, :rt_g_globaltemperature, [0.75,0.77,0.99,1.27,1.62,1.99,3.07,3.90,5.10,6.03])
set_param!(m, :Discontinuity, :y_year_0, 2008.)
set_param!(m, :Discontinuity, :y_year, [2009, 2010, 2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200])
set_param!(m, :Discontinuity, :rgdp_per_cap_NonMarketRemainGDP, readpagedata(m, "test/validationdata/rgdp_per_cap_NonMarketRemainGDP.csv"))
set_param!(m, :Discontinuity, :rcons_per_cap_NonMarketRemainConsumption, readpagedata(m, "test/validationdata/rcons_per_cap_NonMarketRemainConsumption.csv"))
set_param!(m, :Discontinuity, :isatg_saturationmodification, 28.333333333333336)
##running Model
run(m)

@test !isnan(m[:Discontinuity, :irefeqdis_eqdiscimpact][8])
@test !isnan(m[:Discontinuity, :igdpeqdis_eqdiscimpact][10,8])
@test !isnan(m[:Discontinuity, :igdp_realizeddiscimpact][10,8])
@test !isnan(m[:Discontinuity, :occurdis_occurrencedummy][10])
@test !isnan(m[:Discontinuity, :expfdis_discdecay][10])
@test !isnan(m[:Discontinuity, :idis_lossfromdisc][10])
@test !isnan(m[:Discontinuity, :isat_satdiscimpact][10,8])
@test !isnan(m[:Discontinuity, :isat_per_cap_DiscImpactperCapinclSaturation][10,8])
@test !isnan(m[:Discontinuity, :rcons_per_cap_DiscRemainConsumption][10])

#validating - comparison spreadsheet has discontinuity occuring in 2200
#keep running model until m[:Discontinuity,:occurdis_occurrencedummy] shows discontiuity occuring in 2200
output=m[:Discontinuity,:rcons_per_cap_DiscRemainConsumption]
validation=readpagedata(m,"test/validationdata/rcons_per_cap_DiscRemainConsumption.csv")

@test output ≈ validation rtol=1e-2