text:
typealias SCIP_Retcode Int8

# Return codes of SCIP functions. We check these for errros.
SCIP_OKAY               = int8(  1)
SCIP_ERROR              = int8(  0)
SCIP_NOMEMORY           = int8( -1)
SCIP_READERROR          = int8( -2)
SCIP_WRITEERROR         = int8( -3)
SCIP_NOFILE             = int8( -4)
SCIP_FILECREATEERROR    = int8( -5)
SCIP_LPERROR            = int8( -6)
SCIP_NOPROBLEM          = int8( -7)
SCIP_INVALIDCALL        = int8( -8)
SCIP_INVALIDDATA        = int8( -9)
SCIP_INVALIDRESULT      = int8(-10)
SCIP_PLUGINNOTFOUND     = int8(-11)
SCIP_PARAMETERUNKNOWN   = int8(-12)
SCIP_PARAMETERWRONGTYPE = int8(-13)
SCIP_PARAMETERWRONGVAL  = int8(-14)
SCIP_KEYALREADYEXISTING = int8(-15)
SCIP_MAXDEPTHLEVEL      = int8(-16)
SCIP_BRANCHERROR        = int8(-17)

SCIP_RETCODE = Dict{SCIP_Retcode, String}()
SCIP_RETCODE[SCIP_OKAY]               = "normal termination"
SCIP_RETCODE[SCIP_ERROR]              = "unspecified error"
SCIP_RETCODE[SCIP_NOMEMORY]           = "insufficient memory error"
SCIP_RETCODE[SCIP_READERROR]          = "read error"
SCIP_RETCODE[SCIP_WRITEERROR]         = "write error"
SCIP_RETCODE[SCIP_NOFILE]             = "file not found error"
SCIP_RETCODE[SCIP_FILECREATEERROR]    = "cannot create file"
SCIP_RETCODE[SCIP_LPERROR]            = "error in LP solver"
SCIP_RETCODE[SCIP_NOPROBLEM]          = "no problem exists"
SCIP_RETCODE[SCIP_INVALIDCALL]        = "method cannot be called at this time in the solution process"
SCIP_RETCODE[SCIP_INVALIDDATA]        = "error in input data"
SCIP_RETCODE[SCIP_INVALIDRESULT]      = "method returned an invalid result code"
SCIP_RETCODE[SCIP_PLUGINNOTFOUND]     = "a required plugin was not found"
SCIP_RETCODE[SCIP_PARAMETERUNKNOWN]   = "the parameter with the given name was not found"
SCIP_RETCODE[SCIP_PARAMETERWRONGTYPE] = "the parameter is not of the expected type"
SCIP_RETCODE[SCIP_PARAMETERWRONGVAL]  = "the value is invalid for the given parameter"
SCIP_RETCODE[SCIP_KEYALREADYEXISTING] = "the given key is already existing in table"
SCIP_RETCODE[SCIP_MAXDEPTHLEVEL]      = "maximal branching depth level exceeded"
SCIP_RETCODE[SCIP_BRANCHERROR]        = "no branching could be created"


text:
using Random, LinearAlgebra

# TODO: Remove imports
import RecursiveArrayTools.copyat_or_push!

"""
    mcwf_h(tspan, rho0, Hnh, J; <keyword arguments>)

Calculate MCWF trajectory where the Hamiltonian is given in hermitian form.

For more information see: [`mcwf`](@ref)
"""
function mcwf_h(tspan, psi0::T, H::AbstractOperator{B,B}, J::Vector;
        seed=rand(UInt), rates::DecayRates=nothing,
        fout=nothing, Jdagger::Vector=dagger.(J),
        tmp::T=copy(psi0),
        display_beforeevent=false, display_afterevent=false,
        kwargs...) where {B<:Basis,T<:Ket{B}}
    check_mcwf(psi0, H, J, Jdagger, rates)
    f(t::Float64, psi::T, dpsi::T) = dmcwf_h(psi, H, J, Jdagger, dpsi, tmp, rates)
    j(rng, t::Float64, psi::T, psi_new::T) = jump(rng, t, psi, J, psi_new, rates)
    integrate_mcwf(f, j, tspan, psi0, seed, fout;
        display_beforeevent=display_beforeevent,
        display_afterevent=display_afterevent,
        kwargs...)
end

"""
    mcwf_nh(tspan, rho0, Hnh, J; <keyword arguments>)

Calculate MCWF trajectory where the Hamiltonian is given in non-hermitian form.

```math
H_{nh} = H - \\frac{i}{2} \\sum_k J^†_k J_k
```

For more information see: [`mcwf`](@ref)
"""
function mcwf_nh(tspan, psi0::T, Hnh::AbstractOperator{B,B}, J::Vector;
        seed=rand(UInt), fout=nothing,
        display_beforeevent=false, display_afterevent=false,
        kwargs...) where {B<:Basis,T<:Ket{B}}
    check_mcwf(psi0, Hnh, J, J, nothing)
    f(t::Float64, psi::T, dpsi::T) = dmcwf_nh(psi, Hnh, dpsi)
    j(rng, t::Float64, psi::T, psi_new::T) = jump(rng, t, psi, J, psi_new, nothing)
    integrate_mcwf(f, j, tspan, psi0, seed, fout;
        display_beforeevent=display_beforeevent,
        display_afterevent=display_afterevent,
        kwargs...)
end

"""
    mcwf(tspan, psi0, H, J; <keyword arguments>)

Integrate the master equation using the MCWF method.

There are two implementations for integrating the non-hermitian
schroedinger equation:

* [`mcwf_h`](@ref): Usual formulation with Hamiltonian + jump operators
separately.
* [`mcwf_nh`](@ref): Variant with non-hermitian Hamiltonian.

The `mcwf` function takes a normal Hamiltonian, calculates the
non-hermitian Hamiltonian and then calls [`mcwf_nh`](@ref) which is
slightly faster.

# Arguments
* `tspan`: Vector specifying the points of time for which output should
be displayed.
* `psi0`: Initial state vector.
* `H`: Arbitrary Operator specifying the Hamiltonian.
* `J`: Vector containing all jump operators which can be of any arbitrary
operator type.
* `seed=rand()`: Seed used for the random number generator.
* `rates=ones()`: Vector of decay rates.
* `fout`: If given, this function `fout(t, psi)` is called every time an
output should be displayed. ATTENTION: The state `psi` is neither
normalized nor permanent! It is still in use by the ode solve
and therefore must not be changed.
* `Jdagger=dagger.(J)`: Vector containing the hermitian conjugates of the jump
operators. If they are not given they are calculated automatically.
* `display_beforeevent=false`: `fout` is called before every jump.
* `display_afterevent=false`: `fout` is called after every jump.
* `display_jumps=false`: If set to true, an additional list of times and indices
is returned. These correspond to the times at which a jump occured and the index
of the jump operators with which the jump occured, respectively.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function mcwf(tspan, psi0::T, H::AbstractOperator{B,B}, J::Vector;
        seed=rand(UInt), rates::DecayRates=nothing,
        fout=nothing, Jdagger::Vector=dagger.(J),
        display_beforeevent=false, display_afterevent=false,
        kwargs...) where {B<:Basis,T<:Ket{B}}
    isreducible = check_mcwf(psi0, H, J, Jdagger, rates)
    if !isreducible
        tmp = copy(psi0)
        dmcwf_h_(t::Float64, psi::T, dpsi::T) = dmcwf_h(psi, H, J, Jdagger, dpsi, tmp, rates)
        j_h(rng, t::Float64, psi::T, psi_new::T) = jump(rng, t, psi, J, psi_new, rates)
        integrate_mcwf(dmcwf_h_, j_h, tspan, psi0, seed,
            fout;
            display_beforeevent=display_beforeevent,
            display_afterevent=display_afterevent,
            kwargs...)
    else
        Hnh = copy(H)
        if typeof(rates) == Nothing
            for i=1:length(J)
                Hnh -= 0.5im*Jdagger[i]*J[i]
            end
        else
            for i=1:length(J)
                Hnh -= 0.5im*rates[i]*Jdagger[i]*J[i]
            end
        end
        dmcwf_nh_(t::Float64, psi::T, dpsi::T) = dmcwf_nh(psi, Hnh, dpsi)
        j_nh(rng, t::Float64, psi::T, psi_new::T) = jump(rng, t, psi, J, psi_new, rates)
        integrate_mcwf(dmcwf_nh_, j_nh, tspan, psi0, seed,
            fout;
            display_beforeevent=display_beforeevent,
            display_afterevent=display_afterevent,
            kwargs...)
    end
end

"""
    mcwf_dynamic(tspan, psi0, f; <keyword arguments>)

Integrate the master equation using the MCWF method with dynamic
Hamiltonian and Jump operators.

The `mcwf` function takes a normal Hamiltonian, calculates the
non-hermitian Hamiltonian and then calls [`mcwf_nh`](@ref) which is
slightly faster.

# Arguments
* `tspan`: Vector specifying the points of time for which output should
be displayed.
* `psi0`: Initial state vector.
* `f`: Function `f(t, psi) -> (H, J, Jdagger)` or `f(t, psi) -> (H, J, Jdagger, rates)`
    that returns the time-dependent Hamiltonian and Jump operators.
* `seed=rand()`: Seed used for the random number generator.
* `rates=ones()`: Vector of decay rates.
* `fout`: If given, this function `fout(t, psi)` is called every time an
output should be displayed. ATTENTION: The state `psi` is neither
normalized nor permanent! It is still in use by the ode solve
and therefore must not be changed.
* `display_beforeevent=false`: `fout` is called before every jump.
* `display_afterevent=false`: `fout` is called after every jump.
* `display_jumps=false`: If set to true, an additional list of times and indices
is returned. These correspond to the times at which a jump occured and the index
of the jump operators with which the jump occured, respectively.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function mcwf_dynamic(tspan, psi0::T, f::Function;
    seed=rand(UInt), rates::DecayRates=nothing,
    fout=nothing, display_beforeevent=false, display_afterevent=false,
    kwargs...) where {T<:Ket}
    tmp = copy(psi0)
    dmcwf_(t::Float64, psi::T, dpsi::T) = dmcwf_h_dynamic(t, psi, f, rates, dpsi, tmp)
    j_(rng, t::Float64, psi::T, psi_new::T) = jump_dynamic(rng, t, psi, f, psi_new, rates)
    integrate_mcwf(dmcwf_, j_, tspan, psi0, seed,
        fout;
        display_beforeevent=display_beforeevent,
        display_afterevent=display_afterevent,
        kwargs...)
end

"""
    mcwf_nh_dynamic(tspan, rho0, f; <keyword arguments>)

Calculate MCWF trajectory where the dynamic Hamiltonian is given in non-hermitian form.

For more information see: [`mcwf_dynamic`](@ref)
"""
function mcwf_nh_dynamic(tspan, psi0::T, f::Function;
    seed=rand(UInt), rates::DecayRates=nothing,
    fout=nothing, display_beforeevent=false, display_afterevent=false,
    kwargs...) where T<:Ket
    dmcwf_(t::Float64, psi::T, dpsi::T) = dmcwf_nh_dynamic(t, psi, f, dpsi)
    j_(rng, t::Float64, psi::T, psi_new::T) = jump_dynamic(rng, t, psi, f, psi_new, rates)
    integrate_mcwf(dmcwf_, j_, tspan, psi0, seed,
        fout;
        display_beforeevent=display_beforeevent,
        display_afterevent=display_afterevent,
        kwargs...)
end

function dmcwf_h_dynamic(t::Float64, psi::T, f::Function, rates::DecayRates,
                    dpsi::T, tmp::T) where T<:Ket
    result = f(t, psi)
    QO_CHECKS[] && @assert 3 <= length(result) <= 4
    if length(result) == 3
        H, J, Jdagger = result
        rates_ = rates
    else
        H, J, Jdagger, rates_ = result
    end
    QO_CHECKS[] && check_mcwf(psi, H, J, Jdagger, rates_)
    dmcwf_h(psi, H, J, Jdagger, dpsi, tmp, rates_)
end

function dmcwf_nh_dynamic(t::Float64, psi::T, f::Function, dpsi::T) where T<:Ket
    result = f(t, psi)
    QO_CHECKS[] && @assert 3 <= length(result) <= 4
    H, J, Jdagger = result[1:3]
    QO_CHECKS[] && check_mcwf(psi, H, J, Jdagger, nothing)
    dmcwf_nh(psi, H, dpsi)
end

function jump_dynamic(rng, t::Float64, psi::T, f::Function, psi_new::T, rates::DecayRates) where T<:Ket
    result = f(t, psi)
    QO_CHECKS[] && @assert 3 <= length(result) <= 4
    J = result[2]
    if length(result) == 3
        rates_ = rates
    else
        rates_ = result[4]
    end
    jump(rng, t, psi, J, psi_new, rates_)
end

"""
    integrate_mcwf(dmcwf, jumpfun, tspan, psi0, seed; fout, kwargs...)

Integrate a single Monte Carlo wave function trajectory.

# Arguments
* `dmcwf`: A function `f(t, psi, dpsi)` that calculates the time-derivative of
        `psi` at time `t` and stores the result in `dpsi`.
* `jumpfun`: A function `f(rng, t, psi, dpsi)` that uses the random number
        generator `rng` to determine if a jump is performed and stores the
        result in `dpsi`.
* `tspan`: Vector specifying the points of time for which output should
        be displayed.
* `psi0`: Initial state vector.
* `seed`: Seed used for the random number generator.
* `fout`: If given, this function `fout(t, psi)` is called every time an
        output should be displayed. ATTENTION: The state `psi` is neither
        normalized nor permanent! It is still in use by the ode solver
        and therefore must not be changed.
* `kwargs`: Further arguments are passed on to the ode solver.
"""
function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan,
                        psi0::T, seed, fout::Function;
                        display_beforeevent=false, display_afterevent=false,
                        display_jumps=false,
                        save_everystep=false, callback=nothing,
                        alg=OrdinaryDiffEq.DP5(),
                        kwargs...) where T

    # Display before or after events
    function save_func!(affect!,integrator)
        affect!.saveiter += 1
        copyat_or_push!(affect!.saved_values.t, affect!.saveiter, integrator.t)
        copyat_or_push!(affect!.saved_values.saveval, affect!.saveiter,
            affect!.save_func(integrator.u, integrator.t, integrator),Val{false})
        return nothing
    end
    save_before! = display_beforeevent ? save_func! : (affect!,integrator)->nothing
    save_after! = display_afterevent ? save_func! : (affect!,integrator)->nothing

    # Display jump operator index and times
    jump_t = Float64[]
    jump_index = Int[]
    save_t_index = if display_jumps
        function(t,i)
            push!(jump_t,t)
            push!(jump_index,i)
            return nothing
        end
    else
        (t,i)->nothing
    end

    function fout_(x::Vector{ComplexF64}, t::Float64, integrator)
        recast!(x, state)
        fout(t, state)
    end

    state = copy(psi0)
    dstate = copy(psi0)
    out_type = pure_inference(fout, Tuple{eltype(tspan),typeof(state)})
    out = DiffEqCallbacks.SavedValues(Float64,out_type)
    scb = DiffEqCallbacks.SavingCallback(fout_,out,saveat=tspan,
                                         save_everystep=save_everystep,
                                         save_start = false)

    cb = jump_callback(jumpfun, seed, scb, save_before!, save_after!, save_t_index, psi0)
    full_cb = OrdinaryDiffEq.CallbackSet(callback,cb,scb)

    function df_(dx::D, x::D, p, t) where D<:Vector{ComplexF64}
        recast!(x, state)
        recast!(dx, dstate)
        dmcwf(t, state, dstate)
        recast!(dstate, dx)
    end

    prob = OrdinaryDiffEq.ODEProblem{true}(df_, as_vector(psi0),(tspan[1],tspan[end]))

    sol = OrdinaryDiffEq.solve(
                prob,
                alg;
                reltol = 1.0e-6,
                abstol = 1.0e-8,
                save_everystep = false, save_start = false,
                save_end = false,
                callback=full_cb, kwargs...)

    if display_jumps
        return out.t, out.saveval, jump_t, jump_index
    else
        return out.t, out.saveval
    end
end

function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan,
                        psi0::T, seed, fout::Nothing;
                        kwargs...) where T
    function fout_(t::Float64, x::T)
        return normalize(x)
    end
    integrate_mcwf(dmcwf, jumpfun, tspan, psi0, seed, fout_; kwargs...)
end

function jump_callback(jumpfun::Function, seed, scb, save_before!::Function,
                        save_after!::Function, save_t_index::Function, psi0::Ket)

    tmp = copy(psi0)
    psi_tmp = copy(psi0)

    rng = MersenneTwister(convert(UInt, seed))
    jumpnorm = Ref(rand(rng))
    djumpnorm(x::Vector{ComplexF64}, t::Float64, integrator) = norm(x)^2 - (1-jumpnorm[])

    function dojump(integrator)
        x = integrator.u
        t = integrator.t

        affect! = scb.affect!
        save_before!(affect!,integrator)
        recast!(x, psi_tmp)
        i = jumpfun(rng, t, psi_tmp, tmp)
        x .= tmp.data
        save_after!(affect!,integrator)
        save_t_index(t,i)

        jumpnorm[] = rand(rng)
        return nothing
    end

    return OrdinaryDiffEq.ContinuousCallback(djumpnorm,dojump,
                     save_positions = (false,false))
end
as_vector(psi::StateVector) = psi.data

"""
    jump(rng, t, psi, J, psi_new)

Default jump function.

# Arguments
* `rng:` Random number generator
* `t`: Point of time where the jump is performed.
* `psi`: State vector before the jump.
* `J`: List of jump operators.
* `psi_new`: Result of jump.
"""
function jump(rng, t::Float64, psi::T, J::Vector, psi_new::T, rates::Nothing) where T<:Ket
    if length(J)==1
        QuantumOpticsBase.gemv!(complex(1.), J[1], psi, complex(0.), psi_new)
        psi_new.data ./= norm(psi_new)
        i=1
    else
        probs = zeros(Float64, length(J))
        for i=1:length(J)
            QuantumOpticsBase.gemv!(complex(1.), J[i], psi, complex(0.), psi_new)
            probs[i] = dot(psi_new.data, psi_new.data)
        end
        cumprobs = cumsum(probs./sum(probs))
        r = rand(rng)
        i = findfirst(cumprobs.>r)
        QuantumOpticsBase.gemv!(complex(1.)/sqrt(probs[i]), J[i], psi, complex(0.), psi_new)
    end
    return i
end

function jump(rng, t::Float64, psi::T, J::Vector, psi_new::T, rates::Vector{Float64}) where T<:Ket
    if length(J)==1
        QuantumOpticsBase.gemv!(complex(sqrt(rates[1])), J[1], psi, complex(0.), psi_new)
        psi_new.data ./= norm(psi_new)
        i=1
    else
        probs = zeros(Float64, length(J))
        for i=1:length(J)
            QuantumOpticsBase.gemv!(complex(sqrt(rates[i])), J[i], psi, complex(0.), psi_new)
            probs[i] = dot(psi_new.data, psi_new.data)
        end
        cumprobs = cumsum(probs./sum(probs))
        r = rand(rng)
        i = findfirst(cumprobs.>r)
        QuantumOpticsBase.gemv!(complex(sqrt(rates[i]/probs[i])), J[i], psi, complex(0.), psi_new)
    end
    return i
end

"""
Evaluate non-hermitian Schroedinger equation.

The non-hermitian Hamiltonian is given in two parts - the hermitian part H and
the jump operators J.
"""
function dmcwf_h(psi::T, H::AbstractOperator{B,B},
                 J::Vector, Jdagger::Vector, dpsi::T, tmp::T, rates::Nothing) where {B<:Basis,T<:Ket{B}}
    QuantumOpticsBase.gemv!(complex(0,-1.), H, psi, complex(0.), dpsi)
    for i=1:length(J)
        QuantumOpticsBase.gemv!(complex(1.), J[i], psi, complex(0.), tmp)
        QuantumOpticsBase.gemv!(-complex(0.5,0.), Jdagger[i], tmp, complex(1.), dpsi)
    end
    return dpsi
end

function dmcwf_h(psi::T, H::AbstractOperator{B,B},
                 J::Vector, Jdagger::Vector, dpsi::T, tmp::T, rates::Vector{Float64}) where {B<:Basis,T<:Ket{B}}
    QuantumOpticsBase.gemv!(complex(0,-1.), H, psi, complex(0.), dpsi)
    for i=1:length(J)
        QuantumOpticsBase.gemv!(complex(rates[i]), J[i], psi, complex(0.), tmp)
        QuantumOpticsBase.gemv!(-complex(0.5,0.), Jdagger[i], tmp, complex(1.), dpsi)
    end
    return dpsi
end

"""
Evaluate non-hermitian Schroedinger equation.

The given Hamiltonian is already the non-hermitian version.
"""
function dmcwf_nh(psi::T, Hnh::AbstractOperator{B,B}, dpsi::T) where {B<:Basis,T<:Ket{B}}
    QuantumOpticsBase.gemv!(complex(0,-1.), Hnh, psi, complex(0.), dpsi)
    return dpsi
end

"""
    check_mcwf(psi0, H, J, Jdagger, rates)

Check input of mcwf.
"""
function check_mcwf(psi0::Ket{B}, H::AbstractOperator{B,B}, J::Vector, Jdagger::Vector, rates::DecayRates) where B<:Basis
    # TODO: replace type checks by dispatch; make types of J known
    isreducible = true
    if !(isa(H, DenseOperator) || isa(H, SparseOperator))
        isreducible = false
    end
    for j=J
        @assert isa(j, AbstractOperator{B,B})
        if !(isa(j, DenseOperator) || isa(j, SparseOperator))
            isreducible = false
        end
    end
    for j=Jdagger
        @assert isa(j, AbstractOperator{B,B})
        if !(isa(j, DenseOperator) || isa(j, SparseOperator))
            isreducible = false
        end
    end
    @assert length(J) == length(Jdagger)
    if typeof(rates) == Matrix{Float64}
        throw(ArgumentError("Matrix of decay rates not supported for MCWF!
            Use diagonaljumps(rates, J) to calculate new rates and jump operators."))
    elseif typeof(rates) == Vector{Float64}
        @assert length(rates) == length(J)
    end
    isreducible
end

"""
    diagonaljumps(rates, J)

Diagonalize jump operators.

The given matrix `rates` of decay rates is diagonalized and the
corresponding set of jump operators is calculated.

# Arguments
* `rates`: Matrix of decay rates.
* `J`: Vector of jump operators.
"""
function diagonaljumps(rates::Matrix{Float64}, J::Vector{T}) where {B<:Basis,T<:AbstractOperator{B,B}}
    @assert length(J) == size(rates)[1] == size(rates)[2]
    d, v = eigen(rates)
    d, [sum([v[j, i]*J[j] for j=1:length(d)]) for i=1:length(d)]
end

function diagonaljumps(rates::Matrix{Float64}, J::Vector{T}) where {B<:Basis,T<:Union{LazySum{B,B},LazyTensor{B,B},LazyProduct{B,B}}}
    @assert length(J) == size(rates)[1] == size(rates)[2]
    d, v = eigen(rates)
    d, [LazySum([v[j, i]*J[j] for j=1:length(d)]...) for i=1:length(d)]
end

ext: jl


text:
@testset "6.3 Hyperbolic tangent" begin
    include("6.3.1 (c+d x)^m (a+b tanh)^n.jl")
    include("6.3.2 Hyperbolic tangent functions.jl")
    include("6.3.7 (d hyper)^m (a+b (c tanh)^n)^p.jl")
end

max_stars_repo_path: test/6 Hyperbolic functions/6.3 Hyperbolic tangent/.jl


text:
### ============== ### ============ ###
##    Extended Inertial Spin Model   ##
##    (SHARED ARRAYS)                ##
##    Martin Zumaya Hernandez        ##
##    EXAMPLE SIMULATION SCRIPT      ##
### ============== ### ============ ###

### ============ INCLUDE PACKAGES ============ ###

@everywhere using CollectiveDynamics.InertialSpin

### =============== ### =============== ###
###   DEFINITION OF INITIAL PARAMETERS  ###
### =============== ### =============== ###

# N   = 128
# η   = 1.0
# τ   = 6
# rep = 1
# T   = 8*exp10(-5)
# δ   = 0.35

n   = parse(Int, ARGS[1]) # Number of particles
k   = parse(Float64, ARGS[2]) # non local interactions per particle
Tf  = parse(Int, ARGS[3]) # number of iterations
rep = parse(Int, ARGS[4]) # repetition number

# η    = parse(Float64, ARGS[2]) # Dissipation term
# T    = parse(Float64, ARGS[3]) # temperature, noise

@eval @everywhere N    = $n
@eval @everywhere n_nl = $k

@everywhere χ   = 1.25
@everywhere η   = 1.5
@everywhere J   = 0.8
@everywhere T   = 0.01
@everywhere n_t = 6
@everywhere ρ   = 0.3
@everywhere v0  = 0.1

@everywhere pars = InertialParameters(χ, J, η, v0*sqrt(J/χ), ρ, v0, N, T, n_t, n_nl)

# σ = sqrt((2pars.d) * η * T) # noise std deviation ( square root of variance )
@everywhere σ = sqrt((2pars.d) * η * T * pars.dt) # noise std deviation ( square root of variance )
@everywhere L = cbrt(N / pars.ρ) # 3D

# @everywhere κ_dist = Poisson(n_nl)

### =============== ### =============== ###
### SET UP SYSTEM AND OUTPUT STRUCTURE  ###
### =============== ### =============== ###

Rij = SharedArray{Float64}(N, N)

k_nl = SharedArray{Int64}(N) # particles positions

pos  = SharedArray{Float64}(3N) # particles positions
vel  = SharedArray{Float64}(3N) # array of particles' velocities
v_t  = SharedArray{Float64}(3N) # short-range topological interactions
v_nl = SharedArray{Float64}(3N) # long-range topological interactions
spin = SharedArray{Float64}(3N) # long-range topological interactions

for i in 1:length(pos)
    pos[i]  = 2*rand()*L - L
    vel[i]  = 2*rand() - 1
    spin[i] = 2*rand() - 1
    v_t[i]  = 0.0
    v_nl[i] = 0.0
end

# for i in 1:3:length(vel)
for i in 0:N-1
    k_nl[i+1] = 0
    # norm = sqrt(vel[i]^2 + vel[i+1]^2 + vel[i+2]^2)
    # vel[i]   /= norm
    # vel[i+1] /= norm
    # vel[i+2] /= norm
    p_vel = v0 * normalize([vel[3i+1],vel[3i+2],vel[3i+3]])

    vel[3i+1] = p_vel[1]
    vel[3i+2] = p_vel[2]
    vel[3i+3] = p_vel[3]

    p_spin = normalize([spin[3i+1],spin[3i+2],spin[3i+3]])

    temp_spin = normalize(cross(p_spin, p_vel))
    # temp_spin = normalize(cross([spin[i],spin[i+1],spin[i+2]],[vel[i],vel[i+1],vel[i+2]]))

    spin[3i+1] = temp_spin[1]
    spin[3i+2] = temp_spin[2]
    spin[3i+3] = temp_spin[3]
end

output_path = set_output_data_structure_lr("EXTENDED_INERTIAL_SPIN", N, η, T, n_nl)

pos_file  = open(output_path * "/pos_$(rep).dat", "w+")
vel_file  = open(output_path * "/vel_$(rep).dat", "w+")
spin_file = open(output_path * "/spin_$(rep).dat", "w+")

# write initial conditions
println("//////// ", 1)
write(pos_file,  pos)
write(vel_file,  vel)
write(spin_file, spin)

### ============== ### ============== ###
###          SYSTEM EVOLUTION         ###
### ============== ### ============== ###

times = [convert(Int, exp10(i)) for i in 0:Tf]

for i in 1:(length(times) - 1)

    for t in (times[i]+1):times[i+1]

        evolve_extended_system_sh(pos, vel, v_t, v_nl, spin, k_nl, Rij, pars, σ)

        if t % times[i] == 0 || t % div(times[i], exp10(1)) == 0
            println("//////// ", t)
            write(pos_file,  pos)
            write(vel_file,  vel)
            write(spin_file, spin)
        end
    end

end

### ============== ### ============== ### ============== ###

close(pos_file)
close(vel_file)
close(spin_file)

rmprocs(workers())

println("Done all")

### ============== ### ============== ### ============== ###


text:
using DynamicalSystemsBase
using Test, StaticArrays
using DynamicalSystemsBase: DDS, DS
using LinearAlgebra

println("\nTesting discrete system evolution...")

@testset "Logistic Map" begin
  d1 = Systems.logistic(0.1)
  d2 = DDS(d1.f, 0.1, d1.p)
  d3 = DDS(d1.f, big(0.1), d1.p, d1.jacobian)

  @testset "trajectory" begin
    ts1 = trajectory(d1, 100)
    ts3 = trajectory(d3, 100)
    @test ts1[10] ≈ ts3[10]
    @test eltype(ts3) == BigFloat
  end
  @testset "Derivatives" begin
    f1 = jacobian(d1)
    f2 = jacobian(d2)
    f3 = jacobian(d3)
    @test isapprox(f1, f2;rtol = 1e-12)
    @test isapprox(f1, f3;rtol = 1e-12)
    @test typeof(f3) == BigFloat
  end
end

@testset "Henon map" begin

  d1 = Systems.henon()
  d2 = DDS(d1.f, get_state(d1), d1.p)
  d3 = DDS(d1.f, big.(get_state(d1)), d1.p, d1.jacobian)

  @testset "trajectory" begin
    ts1 = trajectory(d1, 100)
    ts3 = trajectory(d3, 100)
    @test ts1[10] ≈ ts3[10]
    @test eltype(ts3) == BigFloat
  end
  @testset "jacobian" begin
    f1 = jacobian(d1)
    f2 = jacobian(d2)
    f3 = jacobian(d3)
    @test isapprox(f1, f2;rtol = 1e-12)
    @test isapprox(f1, f3;rtol = 1e-12)
    @test eltype(f3) == BigFloat
  end
end

@testset "Coupled standard maps" begin
    M = 5; ks = 0.5ones(M); Γ = 0.05;
    ds = Systems.coupledstandardmaps(M, 0.1rand(2M); ks=ks, Γ = Γ)

    u0 = copy(get_state(ds))
    data = trajectory(ds, 100)
    @test u0 != data[end]

    @test jacobian(ds) != jacobian(ds, rand(2M))

    data2 = trajectory(ds, 100; save_idxs = 1:2)
    @test size(data2) == (101, 2)
    data3 = trajectory(ds, 100; save_idxs = 2)
    @test size(data3) == (101, 1)
end

using SparseArrays
@testset "Sparse Matrix Jacobian" begin

ds = Systems.coupledstandardmaps(10)
J = ds.J

sJ = sparse(J)

sparseds = DiscreteDynamicalSystem(ds.f, ds.u0, ds.p, ds.jacobian, sJ)

j = jacobian(sparseds)
@test typeof(j) <: SparseMatrixCSC

tinteg = tangent_integrator(sparseds, 4)

prev = deepcopy(get_deviations(tinteg))
step!(tinteg)

curr = deepcopy(get_deviations(tinteg))
@test prev != curr
end


text:
# ---
# title: 292. Nim Game
# id: problem292
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Brainteaser, Minimax
# link: <https://leetcode.com/problems/nim-game/description/>
# hidden: true
# ---
#
# You are playing the following Nim Game with your friend:
#
#   * Initially, there is a heap of stones on the table.
#   * You and your friend will alternate taking turns, and **you go first**.
#   * On each turn, the person whose turn it is will remove 1 to 3 stones from the heap.
#   * The one who removes the last stone is the winner.
#
# Given `n`, the number of stones in the heap, return `true` _if you can win the
# game assuming both you and your friend play optimally, otherwise
# return_`false`.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 4
#     Output: false
#     Explanation: These are the possible outcomes:
#     1. You remove 1 stone. Your friend removes 3 stones, including the last stone. Your friend wins.
#     2. You remove 2 stones. Your friend removes 2 stones, including the last stone. Your friend wins.
#     3. You remove 3 stones. Your friend removes the last stone. Your friend wins.
#     In all outcomes, your friend wins.
#
#
# **Example 2:**
#
#
#
#     Input: n = 1
#     Output: true
#
#
# **Example 3:**
#
#
#
#     Input: n = 2
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 231 - 1`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Compat
using Compat: Random, mapreduce
import Compat.IteratorSize, Base.length, Compat.IteratorEltype, Base.eltype

using Phylo.API

abstract type AbstractTreeIterator{T <: AbstractTree} end

function IteratorSize(::Type{AbstractTreeIterator})
    return HasLength()
end

function IteratorEltype(::Type{AbstractTreeIterator})
    return HasEltype()
end

abstract type AbstractNodeIterator{T <: AbstractTree} <: AbstractTreeIterator{T} end

function length(ni::It) where It <: AbstractNodeIterator
    return ni.filterfn === nothing ? length(_getnodes(ni.tree)) :
        mapreduce(val -> ni.filterfn(_extractnode(ni.tree, val)) ? 1 : 0,
                  +, ni; init = 0)
end

abstract type AbstractBranchIterator{T <: AbstractTree} <: AbstractTreeIterator{T} end

function length(bi::It) where It <: AbstractBranchIterator
    return bi.filterfn === nothing ? length(_getbranches(bi.tree)) :
        mapreduce(val -> bi.filterfn(_extractbranch(bi.tree, val)) ? 1 : 0,
                  +, bi; init = 0)
end

struct NodeIterator{T <: AbstractTree} <: AbstractNodeIterator{T}
    tree::T
    filterfn::Union{Function, Nothing}
end
"""
    nodeiter(tree::AbstractTree)

Returns an iterator over the nodes of any tree.
"""
nodeiter(tree::T) where T <: AbstractTree =
    NodeIterator{T}(tree, nothing)

"""
    nodefilter(filterfn::Function, tree::AbstractTree)

Returns an iterator over the nodes of any tree, where the
`AbstractNode` is filtered by the function `filterfn`.
"""
nodefilter(filterfn::Function, tree::T) where T <: AbstractTree =
    NodeIterator{T}(tree, filterfn)

eltype(ni::NodeIterator{T}) where T <: AbstractTree = nodetype(ni.tree)

struct NodeNameIterator{T <: AbstractTree} <: AbstractNodeIterator{T}
    tree::T
    filterfn::Union{Function, Nothing}
end

"""
    nodenameiter(tree::AbstractTree)

Returns an iterator over the names of the nodes of any tree.
"""
nodenameiter(tree::T) where T <: AbstractTree =
    NodeNameIterator{T}(tree, nothing)

"""
    nodenamefilter(filterfn::Function, tree::AbstractTree)

Returns an iterator over the nodenames of any tree, where the
`AbstractNode` itself is filtered by the function `filterfn`.
"""
nodenamefilter(filterfn::Function, tree::T) where T <: AbstractTree =
    NodeNameIterator{T}(tree, filterfn)

eltype(ni::NodeNameIterator{T}) where T <: AbstractTree = nodenametype(ni.tree)

struct BranchIterator{T <: AbstractTree} <: AbstractBranchIterator{T}
    tree::T
    filterfn::Union{Function, Nothing}
end

"""
    branchiter(tree::AbstractTree)

Returns an iterator over the branches of any tree.
"""
branchiter(tree::T) where T <: AbstractTree =
    BranchIterator{T}(tree, nothing)

"""
    branchfilter(filterfn::Function, tree::AbstractTree)

Returns an iterator over the branches of any tree, where the
`AbstractBranch` is filtered by the function `filterfn`.
"""
branchfilter(filterfn::Function, tree::T) where T <: AbstractTree =
    BranchIterator{T}(tree, filterfn)

eltype(bi::BranchIterator{T}) where T <: AbstractTree = branchtype(bi.tree)

struct BranchNameIterator{T <: AbstractTree} <: AbstractBranchIterator{T}
    tree::T
    filterfn::Union{Function, Nothing}
end

"""
    branchnameiter(tree::AbstractTree)

Returns an iterator over the names of branches of any tree.
"""
branchnameiter(tree::T) where T <: AbstractTree =
    BranchNameIterator{T}(tree, nothing)

"""
    branchnamefilter(filterfn::Function, tree::AbstractTree)

Returns an iterator over the names of the branches of any tree, where
the `AbstractBranch` is filtered by the function `filterfn`.
"""
branchnamefilter(filterfn::Function, tree::T) where T <: AbstractTree =
    BranchNameIterator{T}(tree, filterfn)

eltype(bi::BranchNameIterator{T}) where T <: AbstractTree = branchnametype(bi.tree)

@static if VERSION >= v"0.7.0"
import Base.iterate
function iterate(ni::NodeIterator, state = nothing)
    nodes = _getnodes(ni.tree)
    if state === nothing
        result = iterate(nodes)
    else
        result = iterate(nodes, state)
    end

    result === nothing && return nothing

    if ni.filterfn === nothing
        return _extractnode(ni.tree, result[1]), result[2]
    end

    val, state = result
    node = _extractnode(ni.tree, val)
    while !ni.filterfn(node)
        result = iterate(nodes, state)
        result === nothing && return nothing
        val, state = result
        node = _extractnode(ni.tree, val)
    end

    return node, state
end

function iterate(ni::NodeNameIterator, state = nothing)
    nodes = _getnodes(ni.tree)
    if state === nothing
        result = iterate(nodes)
    else
        result = iterate(nodes, state)
    end

    result === nothing && return nothing

    if ni.filterfn === nothing
        return _extractnodename(ni.tree, result[1]), result[2]
    end

    val, state = result
    node = _extractnode(ni.tree, val)
    while !ni.filterfn(node)
        result = iterate(nodes, state)
        result === nothing && return nothing
        val, state = result
        node = _extractnode(ni.tree, val)
    end

    name = _extractnodename(ni.tree, val)
    return name, state
end

function iterate(bi::BranchIterator, state = nothing)
    branches = _getbranches(bi.tree)
    if state === nothing
        result = iterate(branches)
    else
        result = iterate(branches, state)
    end

    result === nothing && return nothing

    if bi.filterfn === nothing
        return _extractbranch(bi.tree, result[1]), result[2]
    end

    val, state = result
    branch = _extractbranch(bi.tree, val)
    while !bi.filterfn(branch)
        result = iterate(branches, state)
        result === nothing && return nothing
        val, state = result
        branch = _extractbranch(bi.tree, val)
    end

    return branch, state
end

function iterate(bi::BranchNameIterator, state = nothing)
    branches = _getbranches(bi.tree)
    if state === nothing
        result = iterate(branches)
    else
        result = iterate(branches, state)
    end

    result === nothing && return nothing

    if bi.filterfn === nothing
        return _extractbranchname(bi.tree, result[1]), result[2]
    end

    val, state = result
    branch = _extractbranch(bi.tree, val)
    while !bi.filterfn(branch)
        result = iterate(branches, state)
        result === nothing && return nothing
        val, state = result
        branch = _extractbranch(bi.tree, val)
    end
    name = _extractbranchname(bi.tree, val)

    return name, state
end

else
import Base.start, Base.next, Base.done
function start(ni::It) where It <: AbstractNodeIterator
    nodes = _getnodes(ni.tree)
    state = start(nodes)

    if ni.filterfn === nothing || done(nodes, state)
        return state
    end

    val, after = next(nodes, state)
    while !ni.filterfn(_extractnode(ni.tree, val))
        state = after
        if done(nodes, state)
            return state
        end
        val, after = next(nodes, state)
    end

    return state
end

function done(ni::It, state) where It <: AbstractNodeIterator
    return done(_getnodes(ni.tree), state)
end

function next(ni::NodeIterator, state)
    nodes = _getnodes(ni.tree)
    val, state = next(nodes, state)
    node = _extractnode(ni.tree, val)

    if ni.filterfn === nothing || done(ni, state)
        return node, state
    end

    val, after = next(nodes, state)
    while !ni.filterfn(_extractnode(ni.tree, val))
        state = after
        if done(nodes, state)
            return node, state
        end
        val, after = next(nodes, state)
    end

    return node, state
end

function next(ni::NodeNameIterator, state)
    nodes = getnodes(ni.tree)
    val, state = next(nodes, state)
    node = _extractnode(ni.tree, val)
    name = _extractnodename(ni.tree, val)

    if ni.filterfn === nothing || done(ni, state)
        return name, state
    end

    val, after = next(nodes, state)
    while !ni.filterfn(_extractnode(ni.tree, val))
        state = after
        if done(nodes, state)
            return name, state
        end
        val, after = next(nodes, state)
    end

    return name, state
end

function start(bi::It) where It <: AbstractBranchIterator
    branches = _getbranches(bi.tree)
    state = start(branches)

    if bi.filterfn === nothing || done(branches, state)
        return state
    end

    val, after = next(branches, state)
    while !bi.filterfn(_extractbranch(bi.tree, val))
        state = after
        if done(branches, state)
            return state
        end
        val, after = next(branches, state)
    end

    return state
end

function done(bi::It, state) where It <: AbstractBranchIterator
    return done(_getbranches(bi.tree), state)
end

function next(bi::BranchIterator, state)
    branches = _getbranches(bi.tree)
    val, state = next(branches, state)
    branch = _extractbranch(bi.tree, val)

    if bi.filterfn === nothing || done(bi, state)
        return branch, state
    end

    val, after = next(branches, state)
    while !bi.filterfn(_extractbranch(bi.tree, val))
        state = after
        if done(branches, state)
            return branch, state
        end
        val, after = next(branches, state)
    end

    return branch, state
end

function next(bi::BranchNameIterator, state)
    branches = _getbranches(bi.tree)
    val, state = next(branches, state)
    branch = _extractbranch(bi.tree, val)
    name = _extractbranchname(bi.tree, val)

    if bi.filterfn === nothing || done(bi, state)
        return name, state
    end

    val, after = next(branches, state)
    while !bi.filterfn(_extractbranch(bi.tree, val))
        state = after
        if done(branches, state)
            return name, state
        end
        val, after = next(branches, state)
    end

    return name, state
end

end


text:
@testset "007" begin
    @test tdpe.e007(6) == 13
end

max_stars_repo_name: mithodin/ttpe

max_issues_repo_name: mithodin/ttpe

max_forks_repo_name: mithodin/ttpe


text:
@testset "Multigrid.Multigrid" begin
    @testset "testGMG.jl" begin
      include("testGMG.jl")
    end
    @testset "testGMGRAPforPoisson.jl" begin
       include("testGMGRAPforPoisson.jl")
    end
	@testset "testSAforDivSigGrad.jl" begin
		include("testSAforDivSigGrad.jl")
    end
    @testset "testLinSolveMGWrapper.jl" begin
        include("testLinSolveMGWrapper.jl")
    end
	@testset "testLinSolveAMGWrapper.jl" begin
        include("testLinSolveAMGWrapper.jl")
    end
	@testset "testGMGRAPforElasticity.jl" begin
       include("testGMGRAPforElasticity.jl")
    end
	@testset "testGMGforElasticity.jl" begin
       include("testGMGforElasticity.jl")
    end
	@testset "testGMGRAPforElasticityVanka.jl" begin
        include("testGMGRAPforElasticityVanka.jl")
    end
	@testset "testGMGforElasticityVanka.jl" begin
       include("testGMGforElasticityVanka.jl")
    end
	@testset "testHybridKaczmarz.jl" begin
       include("testHybridKaczmarz.jl")
    end
end


text:
using Test
using StaticArrays, IntervalSets
import ClimaCore.DataLayouts: IJFH
import ClimaCore:
    Fields, slab, Domains, Topologies, Meshes, Operators, Spaces, Geometry

using LinearAlgebra: norm
using ForwardDiff

function spectral_space_2D(; n1 = 1, n2 = 1, Nij = 4)
    domain = Domains.RectangleDomain(
        Geometry.XPoint(-3.0) .. Geometry.XPoint(5.0),
        Geometry.YPoint(-2.0) .. Geometry.YPoint(8.0),
        x1periodic = false,
        x2periodic = false,
        x1boundary = (:east, :west),
        x2boundary = (:south, :north),
    )
    mesh = Meshes.RectilinearMesh(domain, n1, n2)
    grid_topology = Topologies.Topology2D(mesh)

    quad = Spaces.Quadratures.GLL{Nij}()
    space = Spaces.SpectralElementSpace2D(grid_topology, quad)
    return space
end

@testset "1×1 2D domain space" begin
    Nij = 4
    n1 = n2 = 1
    space = spectral_space_2D(n1 = n1, n2 = n2, Nij = Nij)

    field =
        Fields.Field(IJFH{ComplexF64, Nij}(ones(Nij, Nij, 2, n1 * n2)), space)

    @test sum(field) ≈ Complex(1.0, 1.0) * 8.0 * 10.0 rtol = 10eps()
    @test sum(x -> 3.0, field) ≈ 3 * 8.0 * 10.0 rtol = 10eps()
    @test norm(field) ≈ sqrt(2.0 * 8.0 * 10.0) rtol = 10eps()
    @test norm(field, 1) ≈ norm(Complex(1.0, 1.0)) * 8.0 * 10.0 rtol = 10eps()
    @test norm(field, Inf) ≈ norm(Complex(1.0, 1.0)) rtol = 10eps()

    @test extrema(real, field) == (1.0, 1.0)

    @test Operators.matrix_interpolate(field, 4) ≈
          [Complex(1.0, 1.0) for i in 1:(4 * n1), j in 1:(4 * n2)]

    field_sin = map(x -> sin((x.x) / 2), Fields.coordinate_field(space))
    M = Operators.matrix_interpolate(field_sin, 20)
    @test size(M) == (20, 20)  # 20 x 20 for a 1 element field

    real_field = field.re

    # test broadcasting
    res = field .+ 1
    @test parent(Fields.field_values(res)) == Float64[
        f == 1 ? 2 : 1 for i in 1:Nij, j in 1:Nij, f in 1:2, h in 1:(n1 * n2)
    ]

    res = field.re .+ 1
    @test parent(Fields.field_values(res)) ==
          Float64[2 for i in 1:Nij, j in 1:Nij, f in 1:1, h in 1:(n1 * n2)]

end

@testset "Broadcasting interception for tuple-valued fields" begin
    n1 = n2 = 1
    Nij = 4
    space = spectral_space_2D(n1 = n1, n2 = n2, Nij = Nij)

    nt_field = Fields.Field(
        IJFH{NamedTuple{(:a, :b), Tuple{Float64, Float64}}, Nij}(
            ones(Nij, Nij, 2, n1 * n2),
        ),
        space,
    )
    nt_sum = sum(nt_field)
    @test nt_sum isa NamedTuple{(:a, :b), Tuple{Float64, Float64}}
    @test nt_sum.a ≈ 8.0 * 10.0 rtol = 10eps()
    @test nt_sum.b ≈ 8.0 * 10.0 rtol = 10eps()
    @test norm(nt_field) ≈ sqrt(2.0 * 8.0 * 10.0) rtol = 10eps()

    # test scalar asignment
    nt_field.a .= 0.0
    @test sum(nt_field.a) == 0.0
end

@testset "Special case handling for broadcased norm to pass through space local geometry" begin
    space = spectral_space_2D()
    u = Geometry.Covariant12Vector.(ones(space), ones(space))
    @test norm.(u) ≈ hypot(4 / 8 / 2, 4 / 10 / 2) .* ones(space)
end

@testset "FieldVector" begin
    space = spectral_space_2D()
    u = Geometry.Covariant12Vector.(ones(space), ones(space))
    x = Fields.coordinate_field(space)
    y = [1.0, 2.0, 3.0]
    z = 1.0
    Y = Fields.FieldVector(u = u, k = (x = x, y = y, z = z))

    @test propertynames(Y) == (:u, :k)
    @test propertynames(Y.k) == (:x, :y, :z)
    @test Y.u === u
    @test Y.k.x === x
    @test Y.k.y === y
    @test Y.k.z === z

    Y1 = 2 .* Y
    @test parent(Y1.u) == 2 .* parent(u)
    @test parent(Y1.k.x) == 2 .* parent(x)
    @test Y1.k.y == 2 .* y
    @test Y1.k.z === 2 * z

    Y1 .= Y1 .+ 2 .* Y
    @test parent(Y1.u) == 4 .* parent(u)
    @test parent(Y1.k.x) == 4 .* parent(x)
    @test Y1.k.y == 4 .* y
    @test Y1.k.z === 4 * z

    Y.k.z = 3.0
    @test Y.k.z === 3.0
end

@testset "FieldVector basetype replacement and deepcopy" begin
    domain_z = Domains.IntervalDomain(
        Geometry.ZPoint(-1.0) .. Geometry.ZPoint(1.0),
        periodic = true,
    )
    mesh_z = Meshes.IntervalMesh(domain_z; nelems = 10)
    topology_z = Topologies.IntervalTopology(mesh_z)

    domain_x = Domains.IntervalDomain(
        Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0),
        periodic = true,
    )
    mesh_x = Meshes.IntervalMesh(domain_x; nelems = 10)
    topology_x = Topologies.IntervalTopology(mesh_x)

    domain_xy = Domains.RectangleDomain(
        Geometry.XPoint(-1.0) .. Geometry.XPoint(1.0),
        Geometry.YPoint(-1.0) .. Geometry.YPoint(1.0),
        x1periodic = true,
        x2periodic = true,
    )
    mesh_xy = Meshes.RectilinearMesh(domain_xy, 10, 10)
    topology_xy = Topologies.Topology2D(mesh_xy)

    quad = Spaces.Quadratures.GLL{4}()

    space_vf = Spaces.CenterFiniteDifferenceSpace(topology_z)
    space_ifh = Spaces.SpectralElementSpace1D(topology_x, quad)
    space_ijfh = Spaces.SpectralElementSpace2D(topology_xy, quad)
    space_vifh = Spaces.ExtrudedFiniteDifferenceSpace(space_ifh, space_vf)
    space_vijfh = Spaces.ExtrudedFiniteDifferenceSpace(space_ijfh, space_vf)

    space2field(space) = map(
        coord -> (coord, Geometry.Covariant12Vector(1.0, 2.0)),
        Fields.coordinate_field(space),
    )

    Y = Fields.FieldVector(
        field_vf = space2field(space_vf),
        field_if = slab(space2field(space_ifh), 1),
        field_ifh = space2field(space_ifh),
        field_ijf = slab(space2field(space_ijfh), 1, 1),
        field_ijfh = space2field(space_ijfh),
        field_vifh = space2field(space_vifh),
        field_vijfh = space2field(space_vijfh),
        array = [1.0, 2.0, 3.0],
        scalar = 1.0,
    )

    Yf = ForwardDiff.Dual{Nothing}.(Y, 1.0)
    Yf .= Yf .^ 2 .+ Y
    @test all(ForwardDiff.value.(Yf) .== Y .^ 2 .+ Y)
    @test all(ForwardDiff.partials.(Yf, 1) .== 2 .* Y)

    dual_field = Yf.field_vf
    dual_field_original_basetype = similar(Y.field_vf, eltype(dual_field))
    @test eltype(dual_field_original_basetype) === eltype(dual_field)
    @test eltype(parent(dual_field_original_basetype)) === Float64
    @test eltype(parent(dual_field)) === ForwardDiff.Dual{Nothing, Float64, 1}

    object_that_contains_Yf = (; Yf)
    @test axes(deepcopy(Yf).field_vf) === space_vf
    @test axes(deepcopy(object_that_contains_Yf).Yf.field_vf) === space_vf
end


text:
module ModelTypes
export ModelType, LES, SCM

abstract type ModelType end
struct LES <: ModelType end
struct SCM <: ModelType end

end


text:
@doc """
Generate some data using
y = mat*x + n
where x = vector of interest
mat = some transformation matrix
n = noise
This assumes p(n) ~ iid Normal(0, sigma^2)
""" ->
function generate_data(x::Array{Float64, 1}, mat::Array{Float64, 2}, sigma=1.0)
	@assert length(x) == size(mat)[2]
	return mat*x + sigma*randn(size(mat)[1])
end

max_stars_repo_name: eggplantbren/MeanImage

max_issues_repo_name: eggplantbren/MeanImage

max_forks_repo_name: eggplantbren/MeanImage


text:
@testset "Stacks" begin

    @testset "Constructors" begin
        s = Stack{Int}()
        @test isa(s, Stack{Int})

        s = Stack{Int}(5)
        @test isa(s, Stack{Int})

        # deprecated
        # s = Stack(Int)
        # s = Stack(Int, 5)
    end

    @testset "Core Functionality" begin
        s = Stack{Int}(5)
        n = 100

        @test isa(s, Stack{Int})
        @test length(s) == 0
        @test isempty(s)
        @test_throws ArgumentError top(s)
        @test_throws ArgumentError pop!(s)

        for i = 1 : n
            push!(s, i)
            @test top(s) == i
            @test !isempty(s)
            @test length(s) == i
        end

        for i = 1 : n
            x = pop!(s)
            @test x == n - i + 1
            if i < n
                @test top(s) == n - i
            else
                @test_throws ArgumentError top(s)
            end
            @test isempty(s) == (i == n)
            @test length(s) == n - i
        end
    end

    @testset "empty!" begin
        s = Stack{Int}(1)
        push!(s, 10)
        @test length(empty!(s)) == 0
    end

    @testset "iter should return a LIFO collection" begin
        stk = Stack{Int}(10)
        #an array to check iteration sequence against
        arr = Int64[]

        n = 100

        for i = 1:n
            push!(stk, i)
            push!(arr, i)
        end

        @testset "iterator" begin
            index = length(arr)
            for i in stk
                @test(arr[index] == i)
                index -= 1
            end
        end

        @testset "reverse iterator" begin
            index = 1
            for i in reverse_iter(stk)
                @test(arr[index] == i)
                index += 1
            end
        end

        @test arr == [i for i in reverse_iter(stk)]
        @test reverse(arr) == [i for i in stk]
    end

end # @testset Stack


text:
module ScriptIt
export startScript, toggle_unique_key_backend

using Juno
using Atom
using Khepri
using JuliaInterpreter
using Core

isControlFlow = true #true = control flow, false = last call
is_autocad = false #to be used for the unique_key function in order to check reverseHighlight

function toggle_unique_key_backend()
    global is_autocad=!is_autocad
    @info "Unique key for: $(is_autocad ? "autocad backend" : "all backends except autocad")."
end

function script_it_start()
    @info "Welcome to ScriptIt! If you need help press ctrl+alt+h."
end

function extractCode(filePath)
    beforeAll = time()
    #backend(autocad)
    f=  with(traceability, false) do
            all_shapes()
        end
    s=""
    adToolString = "using Khepri"
    backendString = "delete_all_shapes()\nclear_trace!()"
    s = string(adToolString,"\n",backendString, "\n")
    temp_range_to_shape_Dictionary = Dict()
    for el in f
        elementString = string(meta_program(el))
        s = string(s,elementString,"\n")
    end
    afterAll = time()
    s
end

function toggleControlFlow()
    global isControlFlow=!isControlFlow
    @info "Highlight mode: $(isControlFlow ? "control flow" : "last call")."
end

function changeDetected()
    @info "Free-Refactoring modifications detected. Re-running code."
end

function toggleTraceability(toggle)
    @info "Traceability mode: $(toggle ? "On" : "Off")."
end

function toggleFreeForAllRefactoring(toggle)
    @info "Free For All Refactoring mode: $(toggle ? "On" : "Off")."
end

function highlightShape(linesAndTabPath)
    currentTabPath = Symbol(linesAndTabPath[1])
    lines = linesAndTabPath[2]
    shapesToHighlight = Shape[]
    if (lines != "")
        lineList = split(lines,",")
        for line in lineList
            if (line!="")
                lineToHighlight = parse(Int, line)
                append!(shapesToHighlight, source_shapes(currentTabPath, lineToHighlight))
            end
        end
    end
    highlight_shapes(shapesToHighlight)
end

function reverseHighlight()
    highlight_shapes(Shape[])
    @info("Select every shape you want and then press the enter key.");
    selected_shapes = select_shapes()
    actualCodeLinesHighlighted = []
    fileAndLinesToHighlight = Dict()
    stringRanges = ""
    for shape in selected_shapes
        refValVal = unique_key(shape) #TODO see new method with AML
        for dShape in keys(Khepri.shape_to_file_locations)
            dictShape = unique_key(dShape)
            if(refValVal == dictShape && length(Khepri.shape_to_file_locations[dShape])>0)
                if(isControlFlow)
                    actualCodeLinesHighlighted = vcat(actualCodeLinesHighlighted,[Khepri.shape_to_file_locations[dShape]])
                else
                    actualCodeLinesHighlighted = vcat(actualCodeLinesHighlighted,first(Khepri.shape_to_file_locations[dShape]))
                end
            end
        end
    end
    actualCodeLinesHighlighted = unique(actualCodeLinesHighlighted)
    fileAndLinesToHighlight = convertToDictAndToList(actualCodeLinesHighlighted)
    #println("files and lines to highlight: ", fileAndLinesToHighlight)
    return fileAndLinesToHighlight
end

function convertToDictAndToList(list)
    tempDict=Dict()
    file = nothing
    for el in list
        lineJoin = []
        if (isa(el, Tuple)) #This if is necessary to allow the reverseHighlight of last function call mode
            file = string(el[1])
            line = el[2]
            lineJoin = vcat(lineJoin, line)
        else
            for lineIter in el
                file = string(lineIter[1])
                line = lineIter[2]
                lineJoin = vcat(lineJoin, line)
            end
        end
        if(haskey(tempDict,file))
            tempDict[file] = vcat(tempDict[file],[lineJoin])
        else
            tempDict[file] = [lineJoin]
        end
    end
    finalList = []
    for key in keys(tempDict)
        finalList = vcat(finalList, [[key, unique(tempDict[key])]])
    end
    return finalList
end

function unique_key(shape) #TODO improve this thing with AML
    if is_autocad
        Khepri.ACADGetHandleFromShape(connection(current_backend()), Khepri.ref(shape).value)
    else
        shape
    end
end

function isEmpty()
    if (all_shapes() == [])
        return true
    else
        return false
    end
end

function startScript()
    Atom.handle(extractCode, "extractCode")
    Atom.handle(highlightShape, "highlightShape")
    Atom.handle(toggleControlFlow, "toggleControlFlow")
    Atom.handle(changeDetected, "changeDetected")
    Atom.handle(toggleTraceability, "toggleTraceability")
    Atom.handle(toggleFreeForAllRefactoring, "toggleFreeForAllRefactoring")
    Atom.handle(reverseHighlight, "reverseHighlight")
    Atom.handle(isEmpty, "isEmpty")
    Atom.handle(script_it_start, "script_it_start")
end

function __init__()
    @eval Main using Khepri
    traceability(true)
    excluded_modules([Base, Khepri, Atom, Atom.CodeTools, Base.CoreLogging,ScriptIt])
    ScriptIt.startScript()
end

end # module

max_stars_repo_name: Daniel-Fermoselle/Script-It-Julia

max_issues_repo_name: Daniel-Fermoselle/Script-It-Julia

max_forks_repo_name: Daniel-Fermoselle/Script-It-Julia


text:
# This script experiments with AD-ing a composition of functions
# that involves a transform from a property to a number via `val`.
#
# At some point there was a concern that we need to fully rrule the
# outer function as well, but this test here shows this is not necessary.
# Even if the rrules return Vectors of invariants (or similar), then
# Zygote with convert them, using the `ProjectTo`, which is implemented
# in properties.jl:
#   (::ProjectTo{T})(φ::Invariant{T}) where {T} = val(φ)
#
# Another key point to note here (not really related to val I think)
# is that using Zygote inside the rrules means that SVectors are
# converted to SizedVectors; sticking with SVectors requires manual
# implementation of the rrules. But it is not entirely clear to me why.

##

using ACE, StaticArrays, Zygote, ChainRules, LinearAlgebra
using ACE: Invariant, val

import ChainRules: rrule, NoTangent, ZeroTangent

##

function inner(x)
   xi1 = Invariant.(x[1])
   xi2 = Invariant.(x[2])
   return SVector(xi1[1] * xi1[2], xi2[1] + xi2[2])
end

function rrule_inner(dp, x, T=Invariant)
   xi1 = T.(x[1])
   xi2 = T.(x[2])
   o = T(1.0)
   g = [ SA[ dp[1] * xi1[2], dp[1] * xi2[1] ],
         SA[ dp[2] * o, dp[2] * o ] ]
   @show eltype(g)
   return NoTangent(), g
end

rrule(::typeof(inner), x) =
         inner(x), dp -> rrule_inner(dp, x)

function rrule2_inner(dq, dp, x)
   @assert dq[1] isa ZeroTangent
   @assert length(dq) == 2
   dq2 = dq[2]
   @assert dq2 isa AbstractVector{<: SVector}
   # now we need
   #  ∂(dq[2] ⋅ g) / ∂dq and  ∂(dq[2] ⋅ g) / ∂x

   # here I'm hacking a bit so I can use Zygote to take the second
   # derivative for me but then convert back to invariants to make sure
   # the "units" are correct. Note that `g` and hence `dq[2] ⋅ g)` should have
   # both be invariants!!
   f(dp_, x_) = ACE.contract(dq2, rrule_inner(dp_, x_, identity)[2])
   gg = Zygote.gradient(f, dp, x)
   gg_dp = SVector{2, Invariant{Float64}}(gg[1])
   gg_x = SVector{2, Invariant{Float64}}.(gg[2])
   @show eltype(gg_dp)
   @show eltype(gg_x)
   return NoTangent(), gg_dp, gg_x
end

rrule(::typeof(rrule_inner), dp, x) =
         rrule_inner(dp, x), dq -> rrule2_inner(dq, dp, x)

##

outer(m) = sum( val.(m).^2 )

toy(x) = outer(inner(x))

##

x = randn(SVector{2, Float64}, 2)
toy(x)
Zygote.gradient(toy, x)[1]

##
toy2(x) = sum(ACE.normsq, Zygote.gradient(toy, x)[1])
toy2(x)

Zygote.refresh()

Zygote.gradient(toy2, x)


text:
function POMDPs.updater(p::POMCPOWRAVEPlanner)
    rng = MersenneTwister(rand(p.solver.rng, UInt32))
    return BootstrapFilter(p.problem, 10*p.solver.tree_queries, rng)
end


text:
# Copyright (c) 2020 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export forward!

"""
`forward!(alpha,c,Nt,a,A,y)`

In-place version of [`forward`](@ref).
"""
function forward!(alpha::AbstractArray{T},c::Vector{T},
                  Nt::Integer,
                  a::AbstractVector{T},
                  A::AbstractMatrix{T},
                  y::AbstractArray{T}) where {T <: AbstractFloat}
  Ns = size(A,1)
  fill!(alpha,zero(T))
  fill!(c,zero(T))
  for j = 1:Ns
    alpha[1,j] = a[j]*y[1,j]
    c[1] += alpha[1,j]
  end

  # normalize
  for j = 1:Ns
    alpha[1,j] /= c[1]
  end

  for t = 2:Nt
    for j = 1:Ns
      for i = 1:Ns
        aij = A[i,j]
        if aij != T(0)
          alpha[t,j] += alpha[t-1,i]*aij*y[t,j]
        end
      end
      c[t] += alpha[t,j]
    end
    for j = 1:Ns
      alpha[t,j] /= c[t]
    end
  end
  return alpha, c
end

# forward with constrained path
function forward!(alpha::AbstractArray{T},c::Vector{T},
                  Nt::Integer,
                  t2IJ::Vector{Pair{Vector{D},Vector{D}}},
                  A::AbstractMatrix{T},
                  y::AbstractArray{T}) where {D<:Integer, T <: AbstractFloat}
  Ns = size(A,1)

  fill!(alpha,zero(T))
  fill!(c,zero(T))
  idx0 = unique(t2IJ[1][1])
  for i in idx0
    alpha[1,i] = y[1,i]/length(idx0)
    c[1] += alpha[1,i]
  end

  # normalize
  for j = 1:Ns
    alpha[1,j] /= c[1]
  end

  for t = 2:Nt
    I,J = t2IJ[t-1]
    for k in eachindex(I)
      i,j = I[k],J[k]
      alpha[t,j] += alpha[t-1,i]*A[i,j]*y[t,j]
    end
    for j = 1:Ns
      c[t] += alpha[t,j]
    end
    for j = 1:Ns
      alpha[t,j] /= c[t]
    end
  end
  return alpha, c
end

export forward

"""
  `alpha, c = forward(Nt,a,A,y)`

Computes the scaled forward probabilities (``\\bar{\\alpha}``) of an Hidden Markov Model (HMM) of `Ns` states with transition matrix `A` (must be a `Ns` × `Ns` matrix), initial probabilities `a` (must be a `Ns`-long vector) and observation likelihoods `y` (must be of size `Nt2` × `Ns` with `Nt2` ≥ `Nt`).

Returns:
* `alpha` a matrix of size `(Nt,Ns)` containing the forward probabilities
* `c` a `Nt`-long vector containing the normalization coefficients which can be used to compute the backward probabilities (see `backward`) and the log maximum likelihood (`sum(log.(c))`).

`alpha, c = forward(Nt,t2tr,A,y)`

Here `t2tr` can be an `Nt-1` vector of `Pair{Vector{D},Vector{D}}`. For example given `Nt=3`:
```julia
t2tr = [[1,1]=>[1,2],[1,2]=>[3,3],[3]=>[4]]
```
indicates that at time `t=1` transitions from state `1` to `1` and from `1` to `2` are allowed.
At time `t=2` transitions from state `1` to `3` and from `2` to `3` are allowed.
At time `t=3` only the transition from state `3` to `4` is allowed.
In practice these indices can be used to construct the (possibly sparse) time-dependent transition matrices described in the [theory section](@ref theory).
"""
function forward(Nt::Integer,a,A::AbstractArray{T},y) where {T <: AbstractFloat}
  Ns = size(y,2)
  alpha = zeros(T,Nt,Ns)
  c = zeros(T,Nt)
  forward!(alpha,c,Nt,a,A,y)
end

# y is a tensor Nt2 × Ns × Nb
function forward(Nts::Array{D},
                 as::AbstractMatrix{T},
                 As::AbstractArray{T,3},
                 ys::AbstractArray{T,3}) where {D <: Integer,
                                       T <: AbstractFloat}
  Ns = size(As,1)
  Nb = size(ys,3) # batch number
  @assert Nb == length(Nts) == size(as,2) == size(As,3)
  alphas = [zeros(T,Nt,Ns) for Nt in Nts]
  cs = [zeros(T,Nt) for Nt in Nts]
  for i = 1:Nb
    forward!(alphas[i],cs[i],Nts[i],view(as,:,i),view(As,:,:,i),view(ys,:,:,i))
  end
  return alphas, cs
end

# y is a tensor Nt2 × Ns × Nb, constraint path
function forward(Nts::Array{D},
                 as::AbstractVector,
                 A::AbstractMatrix{T},
                 y::Array{T,3}) where {D <: Integer,
                                       T <: AbstractFloat}
  Ns = size(A,1)
  Nb = size(y,3) # batch number
  @assert Nb == length(Nts) == length(as)
  alphas = [zeros(T,Nt,Ns) for Nt in Nts]
  cs = [zeros(T,Nt) for Nt in Nts]
  for i = 1:Nb
    forward!(alphas[i],cs[i],Nts[i],as[i],A,view(y,:,:,i))
  end
  return alphas, cs
end


text:
### Thermal Generation Formulations

abstract type AbstractThermalFormulation <: AbstractDeviceFormulation end

abstract type AbstractThermalDispatchForm <: AbstractThermalFormulation end

struct ThermalUnitCommitment <: AbstractThermalFormulation end

struct ThermalDispatch <: AbstractThermalDispatchForm end

struct ThermalRampLimited <: AbstractThermalDispatchForm end

struct ThermalDispatchNoMin <: AbstractThermalDispatchForm end

# Variables for Thermal Generation

"""
This function add the variables for power generation output to the model
"""
function activepower_variables(ps_m::CanonicalModel,
                               devices::Array{T,1},
                               time_range::UnitRange{Int64}) where {T <: PSY.ThermalGen}

    add_variable(ps_m, devices, time_range, :Pth, false, :var_active)

    return

end

"""
This function add the variables for power generation output to the model
"""
function reactivepower_variables(ps_m::CanonicalModel,
                                 devices::Array{T,1},
                                 time_range::UnitRange{Int64}) where {T <: PSY.ThermalGen}

    add_variable(ps_m, devices, time_range, :Qth, false, :var_reactive)

    return

end

"""
This function add the variables for power generation commitment to the model
"""
function commitment_variables(ps_m::CanonicalModel,
                              devices::Array{T,1},
                              time_range::UnitRange{Int64}) where {T <: PSY.ThermalGen}

    add_variable(ps_m, devices, time_range, :on_th, true)
    add_variable(ps_m, devices, time_range, :start_th, true)
    add_variable(ps_m, devices, time_range, :stop_th, true)

    return

end

"""
This function adds the active power limits of generators when there are no CommitmentVariables
"""
function activepower_constraints(ps_m::CanonicalModel,
                                 devices::Array{T,1},
                                 device_formulation::Type{D},
                                 system_formulation::Type{S},
                                 time_range::UnitRange{Int64}) where {T <: PSY.ThermalGen,
                                                                      D <: AbstractThermalDispatchForm,
                                                                      S <: PM.AbstractPowerFormulation}

    range_data = [(g.name, g.tech.activepowerlimits) for g in devices]

    device_range(ps_m, range_data, time_range, :thermal_active_range, :Pth)

    return

end

"""
This function adds the active power limits of generators when there are CommitmentVariables
"""
function activepower_constraints(ps_m::CanonicalModel,
                                 devices::Array{T,1},
                                 device_formulation::Type{D},
                                 system_formulation::Type{S},
                                 time_range::UnitRange{Int64}) where {T <: PSY.ThermalGen,
                                                                      D <: AbstractThermalFormulation,
                                                                      S <: PM.AbstractPowerFormulation}

    range_data = [(g.name, g.tech.activepowerlimits) for g in devices]

    device_semicontinuousrange(ps_m, range_data, time_range, :thermal_active_range, :Pth, :on_th)

    return

end

"""
This function adds the reactive  power limits of generators when there are CommitmentVariables
"""
function reactivepower_constraints(ps_m::CanonicalModel,
                                   devices::Array{T,1},
                                   device_formulation::Type{D},
                                   system_formulation::Type{S},
                                   time_range::UnitRange{Int64}) where {T <: PSY.ThermalGen,
                                                                        D <: AbstractThermalDispatchForm,
                                                                        S <: PM.AbstractPowerFormulation}

    range_data = [(g.name, g.tech.reactivepowerlimits) for g in devices]

    device_range(ps_m, range_data , time_range, :thermal_reactive_range, :Qth)

    return

end

"""
This function adds the reactive power limits of generators when there CommitmentVariables
"""
function reactivepower_constraints(ps_m::CanonicalModel,
                                   devices::Array{T,1},
                                   device_formulation::Type{D},
                                   system_formulation::Type{S},
                                   time_range::UnitRange{Int64}) where {T <: PSY.ThermalGen,
                                                                        D <: AbstractThermalFormulation,
                                                                        S <: PM.AbstractPowerFormulation}

    range_data = [(g.name, g.tech.reactivepowerlimits) for g in devices]

    device_semicontinuousrange(ps_m, range_data , time_range, :thermal_reactive_range, :Qth, :on_th)

    return

end

"""
This function adds the active power limits of generators when there are no CommitmentVariables
"""
function activepower_constraints(ps_m::CanonicalModel,
                                 devices::Array{T,1},
                                 device_formulation::Type{PSI.ThermalDispatchNoMin},
                                 system_formulation::Type{S},
                                 time_range::UnitRange{Int64}) where {T <: PSY.ThermalGen,
                                                                      S <: PM.AbstractPowerFormulation}

    range_data = [(g.name, (min = 0.0, max=g.tech.activepowerlimits.max)) for g in devices]

    device_range(ps_m, range_data, time_range, :thermal_active_range, :Pth)

    return

end

"""
This function adds the reactive  power limits of generators when there are CommitmentVariables
"""
function reactivepower_constraints(ps_m::CanonicalModel,
                                   devices::Array{T,1},
                                   device_formulation::Type{PSI.ThermalDispatchNoMin},
                                   system_formulation::Type{S},
                                   time_range::UnitRange{Int64}) where {T <: PSY.ThermalGen,
                                                                        S <: PM.AbstractPowerFormulation}

    range_data = [(g.name, (min = 0.0, max=g.tech.reactivepowerlimits.max)) for g in devices]

    device_range(ps_m, range_data , time_range, :thermal_reactive_range, :Qth)

    return

end

### Constraints for Thermal Generation without commitment variables ####

"""
This function adds the Commitment Status constraint when there are CommitmentVariables
"""
function commitment_constraints(ps_m::CanonicalModel,
                                devices::Array{T,1},
                                device_formulation::Type{D},
                                system_formulation::Type{S},
                                time_range::UnitRange{Int64},
                                parameters::Bool) where {T <: PSY.ThermalGen,
                                                                     D <: AbstractThermalFormulation,
                                                                     S <: PM.AbstractPowerFormulation}

    if !(:thermal_status in keys(ps_m.initial_conditions))
        @info("Initial status conditions not provided. This can lead to unwanted results")
        status_init(ps_m, devices, parameters)
    end

    device_commitment(ps_m,
                      ps_m.initial_conditions[:thermal_status],
                      time_range, :commitment_th,
                      (:start_th, :stop_th, :on_th))

    return

end

# ramping constraints

function _get_data_for_rocc(devices::Array{T,1}) where {T <: PSY.ThermalGen}

    set_name = Vector{String}(undef, length(devices))
    ramp_params = Vector{UpDown}(undef, length(devices))
    minmax_params = Vector{MinMax}(undef, length(devices))

    idx = eachindex(devices)
    i, state = iterate(idx)
    for g in devices
        if !isnothing(g.tech.ramplimits)
            set_name[i] = g.name
            ramp_params[i] = g.tech.ramplimits
            minmax_params[i] = g.tech.activepowerlimits
            y = iterate(idx, state)
            y === nothing && (i += 1; break)
            i, state = y
        end
    end

    deleteat!(set_name, i:last(idx))
    deleteat!(ramp_params, i:last(idx))
    deleteat!(minmax_params, i:last(idx))

    return set_name, ramp_params, minmax_params

end

"""
This function adds the ramping limits of generators when there are CommitmentVariables
"""
function ramp_constraints(ps_m::CanonicalModel,
                          devices::Array{T,1},
                          device_formulation::Type{D},
                          system_formulation::Type{S},
                          time_range::UnitRange{Int64},
                          parameters::Bool) where {T <: PSY.ThermalGen,
                                                               D <: AbstractThermalFormulation,
                                                               S <: PM.AbstractPowerFormulation}

    data = _get_data_for_rocc(devices)

    if !isempty(data[2])
        if !(:thermal_output in keys(ps_m.initial_conditions))
            @info("Initial conditions for rate of change not provided. This can lead to unwanted results")
            output_init(ps_m, devices, parameters)
        end

        @assert length(data[2]) == length(ps_m.initial_conditions[:thermal_output])
        # Here goes the reactive power ramp limits
        device_mixedinteger_rateofchange(ps_m,
                                        data,
                                        ps_m.initial_conditions[:thermal_output],
                                        time_range,
                                        :ramp_thermal, (:Pth, :start_th, :stop_th))
    else
        @warn "Data doesn't contain generators with ramp limits, consider adjusting your formulation"
    end

    return

end

function ramp_constraints(ps_m::CanonicalModel,
                          devices::Array{T,1},
                          device_formulation::Type{D},
                          system_formulation::Type{S},
                          time_range::UnitRange{Int64},
                          parameters::Bool) where {T <: PSY.ThermalGen,
                                                               D <: AbstractThermalDispatchForm,
                                                               S <: PM.AbstractPowerFormulation}

    data = _get_data_for_rocc(devices)

    if !isempty(data[2])
        if !(:thermal_output in keys(ps_m.initial_conditions))
            @info("Initial conditions for rate of change not provided. This can lead to unwanted results")
            output_init(ps_m, devices, parameters)
        end
        # Here goes the reactive power ramp limits
        device_linear_rateofchange(ps_m,
                                   (data[1], data[2]),
                                   ps_m.initial_conditions[:thermal_output], time_range,
                                   :ramp_thermal, :Pth)
    else
        @warn "Data doesn't contain generators with ramp limits, consider adjusting your formulation"
    end

    return

end

"""
This function adds the ramping limits of generators when there are CommitmentVariables
"""
function ramp_constraints(ps_m::CanonicalModel,
                          devices::Array{T,1},
                          device_formulation::Type{D},
                          system_formulation::Type{S},
                          time_range::UnitRange{Int64},
                          parameters::Bool) where {T <: PSY.ThermalGen,
                                                               D <: AbstractThermalFormulation,
                                                               S <: PM.AbstractActivePowerFormulation}

    data = _get_data_for_rocc(devices)

    if !isempty(data[2])
        if !(:thermal_output in keys(ps_m.initial_conditions))
            @info("Initial conditions for rate of change not provided. This can lead to unwanted results")
            output_init(ps_m, devices, parameters)
        end

        @assert length(data[2]) == length(ps_m.initial_conditions[:thermal_output])

        device_mixedinteger_rateofchange(ps_m,
                                        data,
                                        ps_m.initial_conditions[:thermal_output],
                                        time_range,
                                        :ramp_thermal, (:Pth, :start_th, :stop_th))

    else
        @warn "Data doesn't contain generators with ramp limits, consider adjusting your formulation"
    end

    return

end

function ramp_constraints(ps_m::CanonicalModel,
                          devices::Array{T,1},
                          device_formulation::Type{D},
                          system_formulation::Type{S},
                          time_range::UnitRange{Int64},
                          parameters::Bool) where {T <: PSY.ThermalGen,
                                                               D <: AbstractThermalDispatchForm,
                                                               S <: PM.AbstractActivePowerFormulation}

    data = _get_data_for_rocc(devices)

    if !isempty(data[2])
        if !(:thermal_output in keys(ps_m.initial_conditions))
            @info("Initial conditions for rate of change not provided. This can lead to unwanted results")
            output_init(ps_m, devices, parameters)
        end
        device_linear_rateofchange(ps_m,
                                    (data[1], data[2]),
                                    ps_m.initial_conditions[:thermal_output], time_range,
                                    :ramp_thermal, :Pth)
    else
        @warn "Data doesn't contain generators with ramp limits, consider adjusting your formulation"
    end

    return

end

# time constraints

function _get_data_for_tdc(devices::Array{T,1}) where {T <: PSY.ThermalGen}

    set_name = Vector{String}(undef, length(devices))
    time_params = Vector{UpDown}(undef, length(devices))

    idx = eachindex(devices)
    i, state = iterate(idx)
    for g in devices
        if !isnothing(g.tech.timelimits)
            set_name[i] = g.name
            time_params[i] = g.tech.timelimits
            y = iterate(idx, state)
            y === nothing && (i += 1; break)
            i, state = y
        end
    end

    deleteat!(set_name, i:last(idx))
    deleteat!(time_params, i:last(idx))

    return set_name, time_params
end

function time_constraints(ps_m::CanonicalModel,
                          devices::Array{T,1},
                          device_formulation::Type{D},
                          system_formulation::Type{S},
                          time_range::UnitRange{Int64},
                          parameters::Bool = false) where {T <: PSY.ThermalGen,
                                                               D <: AbstractThermalFormulation,
                                                               S <: PM.AbstractPowerFormulation}

    duration_data = _get_data_for_tdc(devices)

    if !isempty(duration_data)
        if !(:thermal_duration_on in keys(ps_m.initial_conditions))
            @info("Initial conditions for time up/down not provided. This can lead to unwanted results")
            duration_init(ps_m, devices, parameters)
        end
        @warn("Currently Min-up/down time limits only work if time resoultion they are defined at matches the simulation's time resoultion")
        device_duration_retrospective(ps_m,
                                      duration_data,
                                      ps_m.initial_conditions[:thermal_duration_on],
                                      ps_m.initial_conditions[:thermal_duration_off],
                                      time_range,
                                      :duration, (:on_th, :start_th, :stop_th))
    else
        @warn "Data doesn't contain generators with time-up/down limits, consider adjusting your formulation"
    end

    return

end

# thermal generation cost

function cost_function(ps_m::CanonicalModel,
                       devices::Array{T,1},
                       device_formulation::Type{D},
                       system_formulation::Type{S}) where {T <: PSY.ThermalGen,
                                                           D <: AbstractThermalDispatchForm,
                                                           S <: PM.AbstractPowerFormulation}

    add_to_cost(ps_m, devices, :Pth, :variablecost)

    return

end

function cost_function(ps_m::CanonicalModel,
                       devices::Array{T,1},
                       device_formulation::Type{D},
                       system_formulation::Type{S}) where {T <: PSY.ThermalGen,
                                                           D <: AbstractThermalFormulation,
                                                           S <: PM.AbstractPowerFormulation}

    #Variable Cost component
    add_to_cost(ps_m, devices, :Pth, :variablecost)

    #Commitment Cost Components
    add_to_cost(ps_m, devices, :start_th, :startupcost)
    add_to_cost(ps_m, devices, :stop_th, :shutdncost)
    add_to_cost(ps_m, devices, :on_th, :fixedcost)

    return

end


text:
struct HttpClient
    endpoint:: String
end
export HttpClient

# https://coincheck.com/ja/documents/exchange/api#about
const default_http_client = HttpClient("https://coincheck.com")

function call_http_api(path; method = Methods.GET, credential = Nullable(), args = Dict())
    call_http_api(defaut_http_client, path, method = method, credential = credential, args = args)
end
function call_http_api(client:: HttpClient, path; method = Methods.GET, credential = Nullable(), args = Nullable())
    # URL
    url = HttpUtil.convert_to_url(client.endpoint, path, (method == Methods.GET) ? args : Nullable())
    # Body
    body = (method == Methods.GET || isnull(args)) ? Nullable() : JSON.json(args)

    headers = Nullable()

    if !isnull(credential)
        # nonce
        nonce = string(UInt64(Dates.time() * 1e6))

        message = nonce * url * (isnull(body) ? "" : body)
        signature = Nettle.hexdigest("sha256", credential.secret_key, message)

        headers = Dict{String, String}("ACCESS-KEY" => credential.access_key, "ACCESS-NONCE" => nonce, "ACCESS-SIGNATURE" => signature)
    end

    response = HttpUtil.make_http_request(method, url, headers = headers, body = body)
    return JSON.parse(response.body)
end
export call_http_api


text:
@params mutable struct ReuseAdaptiveSIMP{T} <: AbstractSIMP
    simp::SIMP{T}
    pstart::T
    pfinish::T
    reuse::Bool
    result::SIMPResult{T}
    ftol_cont
    entertol::T
    Δp::T
    maxiter::Int
    innerpolynomial::PolynomialFit{T}
    prev_f::T
    prev_dfdp::T
    prev_ftol::T
    usederivative::Bool
    cheqfilter::CheqFilter
    adapt_pstep::Bool
end

function ReuseAdaptiveSIMP(simp::SIMP{T}, ::Type{Val{filtering}}=Val{false};
    pstart = T(1),
    pfinish = T(5),
    ftol_cont = ExponentialContinuation(T(0.01), T(4log(0.7)), T(0), 100, T(1e-4)), #speed = T(0.1),
    entertol = T(1),
    Δp = T(0.1),
    reuse = true,
    maxiter = 100,
    adalength = 3,
    usederivative = true,
    ratio = 0.5,
    rmin = 1.0,
    adapt_pstep = true) where {T, filtering}

    ncells = getncells(simp.optimizer.obj.problem)
    result = NewSIMPResult(T, simp.optimizer, ncells)

    innerpolynomial = PolynomialFit{T}(2, ratio)
    # `adalength` is the number of steps backward we consider when trying to adapt the penalty steps
    cheqfilter = CheqFilter{filtering}(simp.optimizer.obj.solver, T(rmin))

    return ReuseAdaptiveSIMP(simp, pstart, pfinish, reuse, result, ftol_cont, #=speed,=# entertol, Δp, maxiter, innerpolynomial, T(NaN), T(NaN), T(0), usederivative, cheqfilter, adapt_pstep)
end

function (rasimp::ReuseAdaptiveSIMP)(x0=rasimp.simp.optimizer.obj.solver.vars)
    @unpack pstart, pfinish, Δp = rasimp
    obj = rasimp.simp.optimizer.obj

    # Make an initial inner solve with the pstart creating the MMA workspace in the process
    p = pstart
    workspace = innersolve!(rasimp, x0, p)
    # Record the first value in the polynomial fit objects
    # Make a second data point
    while p < pfinish && obj.fevals < rasimp.maxiter
        #obj.cheqfilter(obj.solver.vars)
        deltap = ceil(abs(rasimp.prev_f) * 4 * MMA.ftol(workspace.model) / dFdp(rasimp.simp.optimizer.obj), 0.1)
        p = p + deltap
        maxiter = workspace.model.maxiter[]
        workspace.converged = false
        innersolve!(rasimp, workspace, p, false)
        # Check that we got past the reuse iteration, if any. If not, keeps going.
        obj.fevals >= rasimp.maxiter && break
    end
    rasimp.result = rasimp.simp.result
    return rasimp.result
end

function innersolve!(asimp::AdaptiveSIMP, x0::AbstractArray{T}, p) where T
    prev_l = length(asimp.simp.topologies)
    prev_fevals = asimp.simp.optimizer.obj.fevals

    innerpolynomial = asimp.innerpolynomial

    # Set penalty as starting penalty
    setpenalty!(asimp.simp, p)
    # Turn off solution reuse
    asimp.simp.optimizer.obj.reuse = false

    # Solve the SIMP with reuse as false and max 3 function evaluations and eps tolerance, 2 iterations with trace recording the objective values

    @unpack model, s_init, s_decr, s_incr, optimizer, suboptimizer, dual_caps, obj, constr, x_abschange, x_converged, f_abschange, f_converged, g_residual, g_converged = asimp.simp.optimizer

    # Resetting the adaptation polynomial
    #reset!(innerpolynomial)

    # Does the first function evaluation
    # Number of function evaluations is the number of iterations plus 1
    workspace = MMA.Workspace(model, x0, optimizer, suboptimizer, s_init=s_init, s_incr=s_incr, s_decr=s_decr, dual_caps=dual_caps)
    workspace.model.maxiter[] = 0
    # Record the first value in the polynomial fit struct
    newvalue!(innerpolynomial, obj.fevals, workspace.f_x)

    # Initial workspace.iter is 0
    # Tolerance is pretty low
    MMA.ftol!(workspace.model, sqrt(eps(T)))
    MMA.grtol!(workspace.model, sqrt(eps(T)))
    MMA.xtol!(workspace.model, sqrt(eps(T)))
    ftol = sqrt(eps(T))
    if obj.fevals >= asimp.maxiter
        return workspace
    end
    mma_results = _innersolve!(asimp, workspace)
    update_result!(asimp.simp, mma_results, prev_l, prev_fevals)
    return workspace
end

function innersolve!(asimp::AdaptiveSIMP, workspace::MMA.Workspace, p, reuse = true)
    prev_l = length(asimp.simp.topologies)
    prev_fevals = asimp.simp.optimizer.obj.fevals

    setpenalty!(asimp.simp, p)
    # Turn on solution reuse
    obj = asimp.simp.optimizer.obj
    obj.reuse = reuse && asimp.reuse

    # Solve the SIMP with reuse as true and max 3 function evaluations with last tolerance, 3 iterations with trace recording the objective values and dCdp

    # Resetting the adaptation polynomial
    #reset!(asimp.innerpolynomial)

    # Number of function evaluations is the number of iterations (not plus one as before)
    workspace.model.maxiter[] += 1
    # Use the last tolerance value
    # This runs 1 dummy iteration since reuse is true
    @debug println("Before MMA")
    @debug @show asimp.simp.optimizer.obj.fevals
    @debug @show MMA.xtol(workspace.model)
    @debug @show MMA.ftol(workspace.model)
    @debug @show MMA.grtol(workspace.model)
    @debug @show workspace.iter
    @debug @show workspace.model.maxiter[]
    @debug @show asimp.simp.optimizer.obj.reuse
    mma_results = asimp.simp.optimizer(workspace)
    @debug @show workspace.iter
    @debug @show asimp.simp.optimizer.obj.fevals
    @debug @show workspace.x_converged
    @debug @show workspace.f_converged
    @debug @show workspace.gr_converged
    @debug println("After MMA")
    newvalue!(asimp.innerpolynomial, obj.fevals, workspace.f_x)
    if workspace.converged || obj.fevals >= asimp.maxiter
        asimp.prev_f = workspace.f_x
        asimp.prev_dfdp = dFdp(asimp.simp.optimizer.obj)
    else
        @debug println("_innersolve!")
        mma_results = _innersolve!(asimp, workspace, 1)
    end
    update_result!(asimp.simp, mma_results, prev_l, prev_fevals)
    return workspace
end

function _innersolve!(asimp::AdaptiveSIMP{T}, workspace::MMA.Workspace, offset=0) where T
    @unpack innerpolynomial, pstart, pfinish = asimp
    p = getpenalty(asimp.simp).p
    obj = asimp.simp.optimizer.obj
    xmin = obj.solver.xmin

    # This runs 2 iterations
    workspace.model.maxiter[] += 1
    mma_results = asimp.simp.optimizer(workspace)
    ftol = MMA.ftol(workspace.model)
    prev_ftol = T(NaN)
    @debug @show ftol
    #@show workspace.model.maxiter[]
    if !workspace.converged && obj.fevals < asimp.maxiter
        # Record the second value in the polynomial fit struct
        newvalue!(innerpolynomial, obj.fevals, workspace.f_x)
        #@assert order(innerpolynomial) == 1
        workspace.model.maxiter[] += 1
        mma_results = asimp.simp.optimizer(workspace)
        if !workspace.converged && obj.fevals < asimp.maxiter
            # Record the third value in the polynomial fit struct
            newvalue!(innerpolynomial, obj.fevals, workspace.f_x)

            # If not converged yet, time to do tolerance adaptation
            # Fit the quadratic
            innerpolycoeffs = solve!(innerpolynomial)
            @debug @show innerpolycoeffs
            # Get the second derivative of the quadratic
            zero_derivative = innerpolynomial(obj.fevals)
            first_derivative = innerpolynomial'(obj.fevals)
            second_derivative = innerpolynomial''(obj.fevals)
            @debug @show second_derivative
            # Set the gradient tolerance as the speed/innerpolynomial''
            prev_ftol = ftol
            progress = log(1 + 1/max(abs(second_derivative/zero_derivative) + abs(first_derivative/zero_derivative), sqrt(eps(T))))
            maxprogress = log(1 + 1/sqrt(eps(T)))
            #progress = p - pstart
            #maxprogress = pfinish - pstart
            #tol_cont = ExponentialContinuation(asimp.speed, T(4log(0.7)), T(0), 100, T(1e-4))
            ##Uncomment##grtol = max(asimp.speed*(1 - progress/maxprogress), sqrt(eps(T)))
            ftol = asimp.ftol_cont(progress/maxprogress)

            #grtol = speed*progress/maxprogress
            #grtol = max(speed*sqrt(abs(second_derivative))/abs(first_derivative), sqrt(eps(T)))
            @debug @show ftol
            #@show workspace.model.maxiter[]
            #MMA.grtol!(workspace.model, grtol)
            #MMA.xtol!(workspace.model, grtol)
            MMA.ftol!(workspace.model, ftol)
            while ftol != prev_ftol && !MMA.assess_convergence(workspace)[end] && obj.fevals < asimp.maxiter
                workspace.model.maxiter[] += 1
                prev_l = length(asimp.simp.topologies)
                prev_fevals = asimp.simp.optimizer.obj.fevals
                prev_f = workspace.f_x
                mma_results = asimp.simp.optimizer(workspace)
                # Should we check if converged here and break??
                # Record the workspace.model.maxiter value in the polynomial fit struct
                #if (obj.fevals != prev_fevals || workspace.f_x != prev_f)
                    newvalue!(innerpolynomial, obj.fevals, workspace.f_x)
                    # Solve for coefficients
                    innerpolycoeffs = solve!(innerpolynomial)
                    @debug @show innerpolycoeffs
                    # Find second derivative
                    zero_derivative =  innerpolynomial(obj.fevals)
                    first_derivative = innerpolynomial'(obj.fevals)
                    second_derivative = innerpolynomial''(obj.fevals)
                    @debug @show second_derivative
                    # Update the gradient tolerance
                    prev_ftol = ftol
                    progress = log(1 + 1/max(abs(second_derivative/zero_derivative) + abs(first_derivative/zero_derivative), sqrt(eps(T))))
                    maxprogress = log(1 + 1/sqrt(eps(T)))
                    #progress = p - pstart
                    #maxprogress = pfinish - pstart
                    ftol = asimp.ftol_cont(progress/maxprogress)
                    #grtol = max(speed*sqrt(abs(second_derivative))/abs(first_derivative), sqrt(eps(T)))
                    @debug @show ftol
                    #MMA.grtol!(workspace.model, grtol)
                    #MMA.xtol!(workspace.model, grtol)
                    MMA.ftol!(workspace.model, ftol)
                #else
                #    @show grtol
                #    break
                #end
            end
        end
    end
    asimp.prev_f = workspace.f_x
    dCdp = dFdp(asimp.simp.optimizer.obj)
    @debug @show dCdp
    asimp.prev_dfdp = dCdp
    asimp.prev_ftol = ftol

    # Last tolerance: grtol
    # Last objective: workspace.f_x
    @debug @show workspace.f_x

    return mma_results
end


text:
using CausalDiscovery
using Test

# @testset "CausalDiscovery.jl" begin
#   # Write your own tests here.
#   include("CISC.jl")
#   include("model.jl")
# end

@testset "Autumn" begin
  # include("Autumn/sexpr.jl")
  # include("Autumn/transform.jl")
  include("Autumn/particles.jl")
  include("Autumn/compile.jl")
end


text:
using Test
using Petri
import Petri: fluxes, odefunc
using ModelingToolkit
using LabelledArrays
using OrdinaryDiffEq
import OrdinaryDiffEq: solve

@testset "Generation of ODE formulas" begin
    @variables S,E,I,R, β,γ,μ
    sir = Petri.Model([S,I,R],[(S+I, 2I), (I,R)])
    seir = Petri.Model([S,E,I,R],[(S+I, E+I), (E,I), (I,R)])
    seirs = Petri.Model([S,E,I,R],[(S+I, E+I), (E,I), (I,R), (R,S)])

    sirf = fluxes(sir)
    @test length(sirf) == 3

    seirf = fluxes(seir)
    @test length(seirf) == 4

    seirsf = fluxes(seirs)
    @test length(seirsf) == 4
end

N(x) = sum(x)

@testset "Generation of ODE solutions" begin
    @variables S,E,I,R, β,γ,μ
    sir = Petri.Model([S,I,R],[(S+I, 2I), (I,R)])
    seir = Petri.Model([S,E,I,R],[(S+I, E+I), (E,I), (I,R)])
    seirs = Petri.Model([S,E,I,R],[(S+I, E+I), (E,I), (I,R), (R,S)])
    fex = odefunc(sir,:sir)
    f = eval(fex)
    u0 = @LArray [100.0, 1, 0] (:S, :I, :R)
    dusir = @LArray [0.0, 0.0, 0.0] (:S, :I, :R)
    p = @LArray [0.35, 0.05] (:μ, :β)
    f(dusir, u0, p, 0.0)
    @show dusir

    @testset "SIR Solving" begin
        u0 = @LArray [100.0, 1, 0] (:S, :I, :R)
        p = @LArray [0.35, 0.05] (:μ, :β)
        fex = odefunc(sir, :sir)
        f = eval(fex)
        prob = ODEProblem(f,u0,(0.0,365.0),p)
        sol = solve(prob,Tsit5())
        @test sol[end].S < 1
        @test sol[end].I < 1
        @test sol[end].R > 99
    end

    @testset "SEIR Solving" begin
        u0 = @LArray [100.0, 1, 0, 0] (:S, :E, :I, :R)
        p = @LArray [0.35, 0.05, 0.05] (:μ, :β, :γ)
        fex = odefunc(seir, :seir)
        f = eval(fex)
        prob = ODEProblem(f,u0,(0.0,365.0),p)
        sol = solve(prob,Tsit5())
        @test sol[end].S < 1
        @test sol[end].E < 1
        @test sol[end].I < 1
        @test sol[end].R > 99
    end

    @testset "SEIRS Solving" begin
        u0 = @LArray [100.0, 1, 0, 0] (:S, :E, :I, :R)
        p = @LArray [0.35, 0.05, 0.07, 0.3] (:μ, :β, :γ, :η)
        fex = odefunc(seirs, :seirs)
        f = eval(fex)
        prob = ODEProblem(f,u0,(0.0,365.0),p)
        sol = solve(prob,Tsit5())
        @test sol[end].S > 5
        @test sol[end].E > 5
        @test sol[end].I > 5
        @test sol[end].R > 5
        @test sol[end].S < 95
        @test sol[end].E < 95
        @test sol[end].I < 95
        @test sol[end].R < 95
    end
end


text:
seatcodes = readlines("./src/Day5/input.txt")
using Chain
seatbin = @chain seatcodes begin
    replace.(_,'R' => '1')
    replace.(_, 'L' => '0')
    replace.(_, 'F' => '0')
    replace.(_, 'B' => '1')
end

rowbin = first.(seatbin, 7)
colbin = last.(seatbin, 3)
rownum = parse.(Int, rowbin, base = 2)
colnum = parse.(Int, colbin, base = 2)

seatid = (8 .* rownum) .+ colnum
max_seat_id = maximum(seatid)

[x for x in 1:max_seat_id if ((x - 1) ∈ seatid) & ((x + 1) ∈ seatid) & (x ∉ seatid)]


text:
#!/bin/bash
# -*- mode: julia -*-
#=
exec julia --project=@. --color=yes --startup-file=no -e "include(popfirst!(ARGS))" "${BASH_SOURCE[0]}" "$@"
=#

using JuMP, Cbc

if length(ARGS) != 1
	print("usage: ./<script_name> <file supported by JuMP.read_from_file>")
end

m = JuMP.read_from_file(ARGS[1])
set_optimizer(m, Cbc.Optimizer)
optimize!(m)
@show JuMP.primal_status(m)
if JuMP.primal_status(m) == MOI.FEASIBLE_POINT
	@show JuMP.objective_value(m)
end


text:
module PDiscRefFEsTests

using Test
using Gridap

D = 2
order = 1

reffe = PDiscRefFE(Float64,D,order)

@test reffe.nfacedofs[end] == [1,2,3]

dofbasis(reffe)
polytope(reffe)
shfbasis(reffe)
nfacedofs(reffe)

end # module


text:
# All functions in this file end up calling `binarydither!` and don't return IndirectArrays.

"""
    dither!([out,] img, alg::AbstractDither, args...; kwargs...)

Dither image `img` using algorithm `alg`.

# Output
If `out` is specified, it will be changed in place. Otherwise `img` will be changed in place.
"""
dither!

"""
    dither([T::Type,] img, alg::AbstractDither, args...; kwargs...)

Dither image `img` using algorithm `alg`.

# Output
If no return type is specified, `dither` will default to the type of the input image.
"""
dither

##############
# Public API #
##############

# If `out` is specified, it will be changed in place...
function dither!(out::GenericImage, img::GenericImage, alg::AbstractDither; kwargs...)
    return _binarydither!(out, img, alg; kwargs...)
end

# ...otherwise `img` will be changed in place.
function dither!(img::GenericImage, alg::AbstractDither; kwargs...)
    tmp = copy(img)
    return _binarydither!(img, tmp, alg; kwargs...)
end

# Otherwise the return type can be chosen...
function dither(::Type{T}, img::GenericImage, alg::AbstractDither; kwargs...) where {T}
    out = similar(Array{T}, axes(img))
    return _binarydither!(out, img, alg; kwargs...)
end

# ...and defaults to the type of the input image.
function dither(img::GenericImage{T,N}, alg::AbstractDither; kwargs...) where {T<:Pixel,N}
    return dither(T, img, alg; kwargs...)
end

#############################
# Low-level algorithm calls #
#############################

# Dispatch to binary dithering on grayscale images
# when no color palette is provided
function _binarydither!(
    out::GenericGrayImage, img::GenericGrayImage, alg::AbstractDither; to_linear=false
)
    to_linear && (img = srgb2linear.(img))
    return binarydither!(alg, out, img)
end

# Dispatch to per-channel dithering on color images when no color palette is provided
function _binarydither!(
    out::GenericImage{T,2},
    img::GenericImage{T,2},
    alg::AbstractDither;
    to_linear=false,
    kwargs...,
) where {T<:Color{<:Real,3}}
    to_linear && (@warn "Skipping transformation `to_linear` when dithering color images.")

    cvout = channelview(out)
    cvimg = channelview(img)
    for c in axes(cvout, 1)
        # Note: the input `out` will be modified
        # since binarydither! modifies the view of the channelview of `out`.
        binarydither!(alg, view(cvout, c, :, :), view(cvimg, c, :, :))
    end
    return out
end


text:
function collision_cross_section_pisigma11(T, nch)

    # pisigma
    pisigma11 = zeros(nch, nch)

    # [A, B, C, D]
    A = 0.0
    B = 0.0
    C = 0.0
    D = 0.0

    # A review of reaction rates and thermodynamic and transport properties
    # for an 11-species air model for chemical and
    # thermal nosequilibrium calculations to 30,000 K

    # N2-N2
    A = 0.0
    B = -0.0112
    C = -0.1182
    D = 4.8464

    pisigma11[1,1] = pisigma(T, A, B, C, D)

    # N2-N
    A = 0.0
    B = -0.0194
    C = 0.0119
    D = 4.1055

    pisigma11[1,2] = pisigma(T, A, B, C, D)
    pisigma11[2,1] = 0.0

    # N-N
    A = 0.0
    B = -0.0033
    C = -0.0572
    D = 5.0452

    pisigma11[2,2] = pisigma(T, A, B, C, D)

    return pisigma11
end

function collision_cross_section_pisigma22(T, nch)

    # pisigma
    pisigma22 = zeros(nch, nch)

    # [A, B, C, D]
    A = 0.0
    B = 0.0
    C = 0.0
    D = 0.0

    # A review of reaction rates and thermodynamic and transport properties
    # for an 11-species air model for chemical and
    # thermal nosequilibrium calculations to 30,000 K

    # N2-N2
    A = 0.0
    B = -0.0203
    C = 0.0683
    D = 4.0900

    pisigma22[1,1] = pisigma(T, A, B, C, D)

    # N2-N
    A = 0.0
    B = -0.0190
    C = 0.0239
    D = 4.1782

    pisigma22[1,2] = pisigma(T, A, B, C, D)
    pisigma22[2,1] = 0.0

    # N-N
    A = 0.0
    B = -0.0118
    C = -0.0960
    D = 4.3252

    pisigma22[2,2] = pisigma(T, A, B, C, D)

    return pisigma22
end

function pisigma(T, A, B, C, D)
    ans = exp(D) * T^(A^2*log(T) + B*log(T) + C)
    return ans
end

function freaction_rate(Tfr, nreaction)
    kf = zeros(nreaction)

    # N2+N2 -> 2N+N2
    Cr = 7.0*10^(21)
    sr = -1.60
    tr = 113200

    kf[1] = reaction_rate_arrhcnius(Tfr, Cr, sr, tr)

    # N2+N  -> 2N+N
    Cr = 3.0*10^(22)
    sr = -1.60
    tr = 113200

    kf[2] = reaction_rate_arrhcnius(Tfr, Cr, sr, tr)

    return kf
end

function reaction_rate_arrhcnius(Tfr, Cr, sr, tr)
    k = Cr * Tfr^(sr) * exp(-tr/Tfr)
    return k
end

function k_eq(T, nsb, nreaction)
    keq = zeros(nreaction)

    # Equilibrium constants curve-fits
    # C. Park, “Nonequilibrium Hypersonic Aerothermodynamics”, Wiley, New York, 1990.

    # N2+N2 -> 2N+N2
    ar1 = [3.490700,
            2.072300,
            1.606000,
            1.535100,
            1.476600,
            1.476600]

    ar2 = [8.313300,
            1.389700,
            1.573200,
            1.606100,
            1.629100,
            1.629100]

    ar3 = [4.097800,
            2.061700,
            1.392300,
            1.299300,
            1.215300,
            1.215300]

    ar4 = [-1.272800,
            -1.182800,
            -1.153300,
            -1.149400,
            -1.145700,
            -1.145700]

    ar5 = [7.487000,
            1.510500,
            -4.543000,
            -6.980000,
            -9.444000,
            -9.444000]

    keq[1] = keq_curve_fitting(T, ar1[nsb], ar2[nsb], ar3[nsb], ar4[nsb], ar5[nsb])

    # N2+N  -> 2N+N
    # 上式と一緒

    keq[2] = keq_curve_fitting(T, ar1[nsb], ar2[nsb], ar3[nsb], ar4[nsb], ar5[nsb])

    return keq
end

function keq_curve_fitting(T, A1, A2, A3, A4, A5)
    Z = 10^4/T
    keq = exp(A1/Z + A2 + A3*log(Z) + A4*Z + A5*Z^2)
    keq = min(keq, 700.0)
    return keq
end

function breaction_rate(kf, keq, nreaction)
    kb = zeros(nreaction)
    for i in 1:nreaction
        kb[i] = kf[i] / keq[i]
    end
    return kb
end

function chemical_enthalpy(T)
    # 複数の温度モードを考慮した方がよい
    delta_hs_0 = [0.0, 3.364e7]  # N2,N[J/kg]
    mw  = [28e-3, 14e-3]          # [kg/mol]
    R   = 8.314                   # [J/K/mol]
    nch = length(delta_hs_0)

    num_molecule = 1             # 二原子分子の数
    theta_vib = [3.353]

    Rs = zeros(nch)              # 化学種気体定数
    for s in 1:nch
        Rs[s] = R / mw[s]
    end

    hs = zeros(nch)
    for s in 1:num_molecule
        etr   = 1.5 * Rs[s] * T
        erot  = Rs[s] * T
        evib  = Rs[s] * theta_vib[s] / (exp(theta_vib[s]/T)-1)
        hs[s] = etr + erot + evib + delta_hs_0[s] + Rs[s]*T
    end
    for s in num_molecule+1:nch
        etr   = 1.5 * Rs[s] * T
        hs[s] = etr + delta_hs_0[s] + Rs[s]*T
    end

    return hs
end