text:
module TestUndirectedWiringDiagrams
using Test

using Catlab.CategoricalAlgebra.FinSets: FinFunction
using Catlab.WiringDiagrams.UndirectedWiringDiagrams

const UWD = UndirectedWiringDiagram

# Imperative interface
######################

# Untyped wiring diagrams.
d = UndirectedWiringDiagram(3)
@test nboxes(d) == 0
@test njunctions(d) == 0
add_box!(d, 2); add_box!(d, 2); add_box!(d, 2)
add_junctions!(d, 4)
@test box(d, 1) == 1
@test box(d, 1:4) == [1,1,2,2]
@test nboxes(d) == 3
@test njunctions(d) == 4
@test collect(boxes(d)) == collect(1:3)
@test collect(junctions(d)) == collect(1:4)
@test collect(ports(d)) == collect(1:6)
@test collect(ports(d, outer=true)) == collect(1:3)
@test ports(d, outer_box(d)) == [1,2,3]
@test ports.([d], boxes(d)) == [[1,2], [3,4], [5,6]]
set_junction!(d, [1,3,5], 1:3)
set_junction!(d, [2,4,6], 4)
set_junction!(d, 1:3, 1:3, outer=true)
@test junction(d, 1) == 1
@test junction(d, 1, outer=true) == 1
@test junction(d, [2,4,6]) == [4,4,4]
@test ports_with_junction(d, 1) == [1]
@test ports_with_junction(d, 1, outer=true) == [1]
@test ports_with_junction(d, 4) == [2,4,6]

# Typed wiring diagrams.
d = UndirectedWiringDiagram([:X,:Y,:Z])
add_box!(d, [:X,:W]); add_box!(d, [:Y,:W]); add_box!(d, [:Z,:W])
add_junctions!(d, [:X,:Y,:Z,:W])
@test port_type(d, 1) == :X
@test port_type(d, [1,3,5]) == [:X,:Y,:Z]
@test port_type(d, 3, outer=true) == :Z
@test junction_type(d, 3) == :Z
set_junction!(d, [1,3,5], 1:3)
set_junction!(d, [2,4,6], 4)
set_junction!(d, 1:3, 1:3, outer=true)
@test junction(d, 1:6) == [1,4,2,4,3,4]
@test junction(d, 1:3, outer=true) == [1,2,3]
@test_throws ErrorException set_junction!(d, 1, 2)

# Interface for ports that is relative to their boxes.
d_previous = d
d = UndirectedWiringDiagram([:X,:Y,:Z])
add_box!(d, [:X,:W]); add_box!(d, [:Y,:W]); add_box!(d, [:Z,:W])
@test port_type(d, (1,1)) == :X
@test port_type(d, (2,1)) == :Y
add_wires!(d, (i,1) => (outer_box(d),i) for i in 1:3)
add_wire!(d, (1,2) => (2,2))
add_wire!(d, (2,2) => (3,2))
@test d == d_previous

# Other constructors
#-------------------

# Singleton diagrams.
d = UndirectedWiringDiagram([:X,:Y,:Z])
add_box!(d, [:X,:Y,:Z])
add_junctions!(d, [:X,:Y,:Z])
set_junction!(d, 1:3); set_junction!(d, 1:3, outer=true)
@test singleton_diagram(UWD, [:X,:Y,:Z]) == d

# Diagrams from cospans.
@test cospan_diagram(UWD, FinFunction(1:3), FinFunction(1:3), [:X,:Y,:Z]) ==
  singleton_diagram(UWD, [:X,:Y,:Z])
@test cospan_diagram(UWD, FinFunction(1:3), FinFunction(1:3)) ==
  singleton_diagram(UWD, 3)

d = UndirectedWiringDiagram(1)
add_box!(d, 2); add_junctions!(d, 1)
set_junction!(d, [1,1]); set_junction!(d, 1, outer=true)
@test cospan_diagram(UWD, FinFunction([1,1]), FinFunction([1])) == d

# Diagrams from functions: junctions only.
d = UndirectedWiringDiagram(4)
add_junctions!(d, 3)
set_junction!(d, [1,2,2,3], outer=true)
@test junction_diagram(UWD, FinFunction([1,2,2,3])) == d

# Operadic interface
####################

f = UndirectedWiringDiagram(3)
add_box!(f, 2); add_box!(f, 2); add_box!(f, 2)
add_junctions!(f, 4)
set_junction!(f, [1,4,2,4,3,4])
set_junction!(f, 1:3, outer=true)

g1 = UndirectedWiringDiagram(2)
add_box!(g1, 1); add_box!(g1, 1)
add_junctions!(g1, 2)
set_junction!(g1, 1:2)
set_junction!(g1, 1:2, outer=true)
g2 = copy(g1)
set_junction!(g2, 2:-1:1)

h = UndirectedWiringDiagram(3)
add_box!(h, 2); add_box!(h, 1); add_box!(h, 1); add_box!(h, 2)
add_junctions!(h, 4)
set_junction!(h, [1,4,4,2,3,4])
set_junction!(h, 1:3, outer=true)
@test ocompose(f,2,g2) == h
@test ocompose(ocompose(f,1,g1),3,g2) == ocompose(ocompose(f,2,g2),1,g1)

h = UndirectedWiringDiagram(3)
for i in 1:6; add_box!(h, 1) end
add_junctions!(h, 4)
set_junction!(h, [1,4,4,2,3,4])
set_junction!(h, 1:3, outer=true)
@test ocompose(f, [g1,g2,g1]) == h

end


text:
"""
    clusterdata = organizeclusters(smd::SMD)

Organize pre-clusters into a vector indexing distinct clusters.

# Description
Pre-clusters in the input `smd`, as related by their shared integer value k_off
`smd.connectID`, are organized into a vector of matrices `clusterdata`.  Each
index of `clusterdata` corresponds to a distinct cluster.
"""
function organizeclusters(smd::SMD)
    # Isolate and sort some arrays from the smd structure.
    connectID = smd.connectID
    sortindices = sortperm(connectID)
    combineddata = [smd.x smd.y smd.x_se smd.y_se smd.framenum smd.datasetnum smd.connectID]
    combineddata = combineddata[sortindices, :]

    # Loop over preclusters and organize member localizations into a more
    # accessible format.
    npercluster = counts(connectID)
    npercluster = npercluster[npercluster .!= 0]
    ncumulative = [0; cumsum(npercluster)]
    uniqueIDs = unique(connectID)
    nclusters = length(uniqueIDs)
    clusterdata = Vector{Matrix{Float32}}(undef, nclusters)
    for nn = 1:nclusters
        clusterind = (1:npercluster[nn]) .+ ncumulative[nn]
        clusterdata[nn] = [combineddata[clusterind, :] sortindices[clusterind]]
    end

    return clusterdata
end

"""
    clusterdata = organizeclusters(smd)

Organize pre-clusters into a vector indexing distinct clusters.

# Description
Pre-clusters in the input `smd`, as related by their shared integer value k_off
`smd.connectID`, are organized into a vector of matrices `clusterdata`.  Each
index of `clusterdata` corresponds to a distinct cluster.
"""
function organizeclusters(smd)
    # Isolate and sort some arrays from the smd structure.
    connectID = smd.connectID
    sortindices = sortperm(connectID)
    combineddata = [smd.x smd.y smd.x_se smd.y_se smd.framenum smd.datasetnum smd.connectID]
    combineddata = combineddata[sortindices, :]

    # Loop over preclusters and organize member localizations into a more
    # accessible format.
    npercluster = counts(connectID)
    npercluster = npercluster[npercluster .!= 0]
    ncumulative = [0; cumsum(npercluster)]
    uniqueIDs = unique(connectID)
    nclusters = length(uniqueIDs)
    clusterdata = Vector{Matrix{Float32}}(undef, nclusters)
    for nn = 1:nclusters
        clusterind = (1:npercluster[nn]) .+ ncumulative[nn]
        clusterdata[nn] = [combineddata[clusterind, :] sortindices[clusterind]]
    end

    return clusterdata
end

ext: jl


text:
### Effective sample size (ESS)

ess(mcvariance::N, iidvariance::N, len::Integer) where {N<:Real} = len*iidvariance/mcvariance

ess(v::AbstractVector, vtype::Symbol, args...) = ess(mcvar(v, Val{vtype}, args...), mcvar(v, Val{:iid}), length(v))

ess(v::AbstractArray, vtype::Symbol, region, args...) = mapslices(x -> ess(vec(x), vtype, args...), v, region)

ess(s::VariableNState{Univariate}, vtype::Symbol, args...) = ess(s.value, vtype, args...)

ess(s::VariableNState{Multivariate}, vtype::Symbol, i::Integer, args...) = ess(s.value[i, :], vtype, args...)

ess(s::VariableNState{Multivariate}, vtype::Symbol, r::AbstractVector=1:s.size, args...) =
  eltype(s)[ess(s, vtype, i, args...) for i in r]

ess(v::AbstractVector, args...) = ess(v, :imse, args...)

ess(v::AbstractArray, region, args...) = ess(v, :imse, region, args...)

ess(s::VariableNState{Univariate}, args...) = ess(s.value, :imse, args...)

ess(s::VariableNState{Multivariate}, i::Integer, args...) = ess(s, :imse, i, args...)

ess(s::VariableNState{Multivariate}, r::AbstractVector=1:s.size, args...) = ess(s, :imse, r, args...)


text:
"""
```
blanchard_kahn(m::RiskAdjustedLinearization; deterministic::Bool = false, verbose::Symbol = :high)
```

checks the Blanchard-Kahn conditions for whether a first-order perturbation is saddle-path stable or not.

If `verbose` is `:low` or `:high`, a print statement will be shown if the Blanchard-Kahn conditions are satisfied.
"""
function blanchard_kahn(m::RiskAdjustedLinearization; deterministic::Bool = false, verbose::Symbol = :high)

    li = linearized_system(m)

    Nz = m.Nz
    Ny = m.Ny
    N_zy = m.Nz + m.Ny
    ztype = eltype(m.z)
    AA = Matrix{ztype}(undef, N_zy, N_zy)
    BB = similar(AA)

    # Populate AA
    AA[1:Ny, 1:Nz]                 = li[:Γ₅]
    AA[1:Ny, (Nz + 1):end]         = li[:Γ₆]
    AA[(Ny + 1):end, 1:Nz]         = Matrix{ztype}(I, m.Nz, m.Nz) # faster but makes allocations, unlike Diagonal(Ones{ztype}(Nz))
    AA[(Ny + 1):end, (Nz + 1):end] = Zeros{ztype}(m.Nz, m.Ny)

    # Populate BB
    BB[1:Ny, 1:Nz]                 = deterministic ? -li[:Γ₃] : -(li[:Γ₃] + li[:JV])
    BB[1:Ny, (Nz + 1):end]         = -li[:Γ₄]
    BB[(Ny + 1):end, 1:Nz]         =  li[:Γ₁]
    BB[(Ny + 1):end, (Nz + 1):end] =  li[:Γ₂]

    if count(abs.(eigen(AA, BB).values) .> 1) != m.Nz
        if deterministic
            throw(BlanchardKahnError("First-order perturbation around deterministic steady state is not saddle-path stable"))
        else
            throw(BlanchardKahnError("First-order perturbation around stochastic steady state is not saddle-path stable"))
        end
    else
        if verbose in [:low, :high]
            if deterministic
                println("Blanchard-Kahn conditions for a unique locally bounded deterministic " *
                        "steady-state perturbation are satisfied")
            else
                println("Blanchard-Kahn conditions for a unique locally bounded stochastic " *
                        "steady-state perturbation are satisfied")
            end
        end

        return true
    end
end

mutable struct BlanchardKahnError <: Exception
    msg::String
end
BlanchardKahnError() = BlanchardKahnError("First-order perturbation is not saddle-path stable")
Base.showerror(io::IO, ex::BlanchardKahnError) = print(io, ex.msg)


text:
# Temporary solution type has superset of all results
struct StateSpaceSolution{T1,T2,T3,T4,T5}
    z::T1 # observables, if relevant
    u::T2 # hidden state, or mean of prior if filtering/estimating
    W::T3 # shocks, if not filtering
    P::T4 # prior variances if filtering, would prefer actually having priors instead
    logpdf::T5  #likelihood of observables
end

## Some algorithm types for dispatching
# LTI = Linear Time Invariant
# QTI = Quadratic Time Invariant?  Not sure term, but mean 2nd order quadratic form
abstract type AbstractStateSpaceAlgorithm end
struct GeneralTimeVarying <: AbstractStateSpaceAlgorithm end
struct LTI <: AbstractStateSpaceAlgorithm end
struct QTI <: AbstractStateSpaceAlgorithm end

# Hopefully temporary algorithm types.  Allows us to optimize a likelihood only solution
# then later we can ensure there is no overhead
struct LTILikelihood <: AbstractStateSpaceAlgorithm end
struct QTILikelihood <: AbstractStateSpaceAlgorithm end

#############################################
## Fully general solve
# Calculates the logpdf if observables passed in.
## Can tweak and eventually check lack of overhead
function solve(
    f,
    g,
    u0,
    tspan,
    p = nothing,
    alg = GeneralTimeVarying();
    h = (u, p, t) -> u,
    D = nothing,
    noise = nothing,
    observables = nothing,
    save_everystep = true,
    save_posteriors = true,
    calculate_logpdf = true,
    simulate_observation_noise = false,
)
    return _solve(
        alg,
        noise,
        observables,
        f,
        g,
        h,
        D,
        u0,
        tspan,
        p,
        save_everystep,
        save_posteriors,
        calculate_logpdf,
        simulate_observation_noise,
    )
end

# Likelhood remains 0 if no obserables.  Need to verify no overhead
maybe_logpdf(observables::Nothing, D, z, i) = 0.0
maybe_logpdf(observables, D, z, i) = logpdf(D, observables[i] .- z)

function _solve(
    alg::GeneralTimeVarying,
    noise,
    observables,
    f,
    g,
    h,
    D,
    u0,
    tspan,
    p,
    save_everystep,
    save_posteriors,
    calculate_logpdf,
    simulate_observation_noise,
)
    T = tspan[2]
    @assert tspan[1] == 0
    @assert length(noise) == T

    z0 = h(u0, p, 0)
    u = Zygote.Buffer(Vector{typeof(u0)}(undef, T + 1))
    z = Zygote.Buffer(Vector{typeof(z0)}(undef, T + 1))
    u[1] = u0
    z[1] = z0
    loglik = 0.0  # remains 0 if no observables
    for t = 2:T+1
        u[t] = f(u[t-1], p, t - 1) .+ g(u[t-1], p, t - 1) * noise[t-1]
        z[t] = h(u[t], p, t)
        loglik += maybe_logpdf(observables, D, z[t], t - 1)  # z_0 doesn't enter likelihood
    end
    return StateSpaceSolution(copy(z), copy(u), noise, nothing, loglik)
end

############################################
# QTI Specializations
# This is doing state-space pruning, so the u0 will be the size of the underlying state rather than the doubling of it used in the general observation equations
solve(
    A_0,
    A_1,
    A_2,
    B,
    C_0,
    C_1,
    C_2,
    D,
    u0,
    tspan,
    alg = QTI();
    noise = nothing,
    observables = nothing,
) = _solve(alg, noise, observables, A_0, A_1, A_2, B, C_0, C_1, C_2, D, u0, tspan)

# With a symmetric A_2 tensor, there should be a specializaiton for the rrule of the quadratic form
#=
function quad(A_2, x)
    n = size(A_2, 1)
    return map( j -> x' * A_2[j, :, :] * x, 1:n)
end
=#
quad(A::AbstractArray{<:Number,3}, x) =
    @matmul c[l] := sum(j) (@matmul [l, j] := sum(k) A[l, j, k] * x[k]) * x[j]

function ChainRulesCore.rrule(::typeof(quad), A::AbstractArray{<:Number,3}, x)
    res = quad(A, x)
    function quad_pb(Δres)
        ΔA = similar(A)
        Δx = zeros(length(x))
        tmp = x * x'
        n = size(A, 1)
        for i in 1:n
            ΔA[i, :, :] .= tmp .* Δres[i]
            Δx += (A[i, :, :] + A[i, :, :]') * x .* Δres[i]
        end
        return NoTangent(), ΔA, Δx
    end
    return res, quad_pb
end

# QTI with noise and observables is a joint likelihood
function _solve(alg::QTI, noise, observables, A_0, A_1, A_2, B, C_0, C_1, C_2, D, u0, tspan)
    T = tspan[2]
    @assert tspan[1] == 0
    @assert length(noise) == T

    z0 = C_0 + C_1 * u0 + quad(C_2, u0)  # the same u0 used for initial observation
    u_f = Zygote.Buffer(Vector{typeof(u0)}(undef, T + 1))
    u = Zygote.Buffer(Vector{typeof(u0)}(undef, T + 1))
    z = Zygote.Buffer(Vector{typeof(z0)}(undef, T + 1))
    u[1] = u0
    u_f[1] = u0
    z[1] = z0
    loglik = 0.0  # remains 0 if no observables
    for t = 2:T+1
        u_f[t] = A_1 * u_f[t-1] .+ B * noise[t-1]
        u[t] = A_0 + A_1 * u[t-1] + quad(A_2, u_f[t-1]) .+ B * noise[t-1]
        z[t] = C_0 + C_1 * u[t] + quad(C_2, u_f[t])
        err = observables[t-1] - z[t]
        loglik += logpdf(D, err)  # z_0 doesn't enter likelihood
    end
    return StateSpaceSolution(copy(z), copy(u), noise, nothing, loglik)
end

# QTI with noise and observables is a joint likelihood
function _solve(
    alg::QTILikelihood,
    noise,
    observables,
    A_0,
    A_1,
    A_2,
    B,
    C_0,
    C_1,
    C_2,
    D,
    u0,
    tspan,
)
    T = tspan[2]
    @assert tspan[1] == 0
    @assert length(noise) == T

    z0 = C_0 + C_1 * u0 + quad(C_2, u0)  # the same u0 used for initial observation
    u_f = Zygote.Buffer(Vector{typeof(u0)}(undef, T + 1))
    u = Zygote.Buffer(Vector{typeof(u0)}(undef, T + 1))
    z = Zygote.Buffer(Vector{typeof(z0)}(undef, T + 1))
    u[1] = u0
    u_f[1] = u0
    z[1] = z0
    loglik = 0.0  # remains 0 if no observables
    for t = 2:T+1
        u_f[t] = A_1 * u_f[t-1] .+ B * noise[t-1]
        u[t] = A_0 + A_1 * u[t-1] + quad(A_2, u_f[t-1]) .+ B * noise[t-1]
        z[t] = C_0 + C_1 * u[t] + quad(C_2, u_f[t])
        err = observables[t-1] - z[t]
        loglik += logpdf(D, err)  # z_0 doesn't enter likelihood
    end
    return StateSpaceSolution(copy(z), copy(u), noise, nothing, loglik)
end

# QTI with noise and observables is a joint likelihood
function _solve(
    alg::QTI,
    noise,
    observables::Nothing,
    A_0,
    A_1,
    A_2,
    B,
    C_0,
    C_1,
    C_2,
    D,
    u0,
    tspan,
)
    T = tspan[2]
    @assert tspan[1] == 0
    @assert length(noise) == T

    z0 = C_0 + C_1 * u0 + quad(C_2, u0)  # the same u0 used for initial observation
    u_f = Zygote.Buffer(Vector{typeof(u0)}(undef, T + 1))
    u = Zygote.Buffer(Vector{typeof(u0)}(undef, T + 1))
    z = Zygote.Buffer(Vector{typeof(z0)}(undef, T + 1))
    u[1] = u0
    u_f[1] = u0
    z[1] = z0
    loglik = 0.0  # remains 0 if no observables
    for t = 2:T+1
        u_f[t] = A_1 * u_f[t-1] .+ B * noise[t-1]
        u[t] = A_0 + A_1 * u[t-1] + quad(A_2, u_f[t-1]) .+ B * noise[t-1]
        z[t] = C_0 + C_1 * u[t] + quad(C_2, u_f[t])
    end
    return StateSpaceSolution(copy(z), copy(u), noise, nothing, loglik)
end

## General function mapping for the Perturbation Soultions
# ideally, we wouldn't need the specializations for the QTI/etc if these are just as efficient

# First order
dssm_evolution(x, sol::FirstOrderPerturbationSolution, t) = sol.A * x
dssm_volatility(x, sol::FirstOrderPerturbationSolution, t) = sol.B
dssm_observation(x, sol::FirstOrderPerturbationSolution, t) = sol.C * x

# The 2nd order are trickier because of state-space pruning
function dssm_evolution(x, sol::SecondOrderPerturbationSolution, t)
    # assume x is stacked as [x_f, x]
    # Try a view?
    x_f_new = sol.A_1 * x[1:sol.n_x]
    x_new = sol.A_0 + sol.A_1 * x[(sol.n_x+1):end] + quad(sol.A_2, x[1:sol.n_x])
    return [x_f_new; x_new]
end

dssm_volatility(x, sol::SecondOrderPerturbationSolution, t) = [sol.B; sol.B]

function dssm_observation(x, sol::SecondOrderPerturbationSolution, t)
    # assume x is stacked as [x_f, x]
    # Try a view?
    return sol.C_0 + sol.C_1 * x[(sol.n_x+1):end] + quad(sol.C_2, x[1:sol.n_x])
end

## Mapping to the Pertubation solutions

# Default algorithm is the likelihood-only for now
solve(sol::FirstOrderPerturbationSolution, u0, tspan, alg = LTILikelihood(); kwargs...) =
    solve(sol.A, sol.B, sol.C, sol.D, u0, tspan, alg; kwargs...)

solve(sol::SecondOrderPerturbationSolution, u0, tspan, alg = QTILikelihood(); kwargs...) =
    solve(
        sol.A_0,
        sol.A_1,
        sol.A_2,
        sol.B,
        sol.C_0,
        sol.C_1,
        sol.C_2,
        sol.D,
        u0,
        tspan,
        alg;
        kwargs...,
    )

############################################
# Utility functions
function save_model_results(file_prefix, path, sol::FirstOrderPerturbationSolution, ϵ, z)
    try
        # For first-order ones, save A, B, C, D
        CSV.write(
            joinpath(path, "$(file_prefix)_A.csv"),
            DataFrame(sol.A),
            writeheader = false,
        )
        CSV.write(
            joinpath(path, "$(file_prefix)_B.csv"),
            DataFrame(sol.B),
            writeheader = false,
        )
        CSV.write(
            joinpath(path, "$(file_prefix)_C.csv"),
            DataFrame(sol.C),
            writeheader = false,
        )
        CSV.write(
            joinpath(path, "$(file_prefix)_D.csv"),
            DataFrame(column_name = sol.D.σ),
            writeheader = false,
        ) # Extract the diagonal
        # In addition, save ϵ and z with different file names
        CSV.write(
            joinpath(path, "$(file_prefix)_noise.csv"),
            DataFrame(transpose(hcat(ϵ...))),
            writeheader = false,
        )
        CSV.write(
            joinpath(path, "$(file_prefix)_observables.csv"),
            DataFrame(transpose(hcat(z...))),
            writeheader = false,
        )
        # Save the ergodic distribution for first order
        CSV.write(
            joinpath(path, "$(file_prefix)_ergodic.csv"),
            DataFrame(sol.x_ergodic.C.L * sol.x_ergodic.C.U),
            writeheader = false,
        )
    catch err
        println("Errors occur when saving the CSV files for first-order")
    end
end

function save_model_results(file_prefix, path, sol::SecondOrderPerturbationSolution, ϵ, z)
    try
        # For second-order ones, save A_0, A_1, A_2, B, C_0, C_1, C_2, D
        # Notice that the B and D objects are the same as those in the first-order case
        CSV.write(
            joinpath(path, "$(file_prefix)_A_0.csv"),
            DataFrame(column_name = sol.A_0),
            writeheader = false,
        ) # Vector
        CSV.write(
            joinpath(path, "$(file_prefix)_A_1.csv"),
            DataFrame(sol.A_1),
            writeheader = false,
        )
        CSV.write(
            joinpath(path, "$(file_prefix)_A_2.csv"),
            DataFrame(column_name = vec(sol.A_2)),
            writeheader = false,
        ) # Save as a vector, reshape it when loading
        CSV.write(
            joinpath(path, "$(file_prefix)_B.csv"),
            DataFrame(sol.B),
            writeheader = false,
        )
        CSV.write(
            joinpath(path, "$(file_prefix)_C_0.csv"),
            DataFrame(column_name = sol.C_0),
            writeheader = false,
        ) # Vector
        CSV.write(
            joinpath(path, "$(file_prefix)_C_1.csv"),
            DataFrame(sol.C_1),
            writeheader = false,
        )
        CSV.write(
            joinpath(path, "$(file_prefix)_C_2.csv"),
            DataFrame(column_name = vec(sol.C_2)),
            writeheader = false,
        ) # Save as a vector, reshape it when loading
        CSV.write(
            joinpath(path, "$(file_prefix)_D.csv"),
            DataFrame(column_name = sol.D.σ),
            writeheader = false,
        ) # Extract the diagonal
        # In addition, save ϵ and z with different file names
        CSV.write(
            joinpath(path, "$(file_prefix)_noise.csv"),
            DataFrame(transpose(hcat(ϵ...))),
            writeheader = false,
        )
        CSV.write(
            joinpath(path, "$(file_prefix)_observables.csv"),
            DataFrame(transpose(hcat(z...))),
            writeheader = false,
        )
    catch err
        println("Errors occur when saving the CSV files for second-order")
    end
end


text:
"""
"""
function is_one_to_one(x::AbstractDict)::Bool
    if (length(keys(x)) == length(unique(keys(x)))) &&
            (length(values(x)) == length(unique(values(x))))
        return true
    else
        return false
    end
end

"""
"""
function inverse(x::AbstractDict)::Dict
    if !is_one_to_one(x)
        error(
            string(
                "Input directory is not one-to-one."
                )
            )
    end

    keys_array = fix_type(collect(keys(x)))
    values_array = fix_type(collect(values(x)))

    key_type=eltype(keys_array)
    value_type=eltype(values_array)

    result = Dict{value_type, key_type}()

    for i = 1:length(keys_array)
        k = keys_array[i]
        v = values_array[i[]]
        result[v] = k
    end

    result_typefixed = fix_type(result)

    return result_typefixed
end


text:
"""
    smoothstep(x, x_step, dx, a=zero(x), b=one(x))

smooth step function:
- = a for x <= x_step - dx
- = b for x >= x_step + dx
- = (a-b)/2 for x == x_step
and smooth in between.
"""
function smoothstep(x, x_step, dx, a=zero(x), b=one(x))
    @assert dx > 0
    edge0 = x_step - dx
    x <= edge0 && return(a)
    edge1 = x_step + dx
    x >= edge1 && return(b)
    y0 = _smoothstep(x,edge0,edge1)
    y0 * (b-a) + a
end

function _smoothstep(x, edge0, edge1)
    #https://en.wikipedia.org/wiki/Smoothstep
    # Scale, bias and saturate x to 0..1 range
    # assume that
    @assert edge0 < x < edge1
    #x = clamp((x - edge0) / (edge1 - edge0), 0.0, 1.0);
    xs = (x - edge0) / (edge1 - edge0)
    # Evaluate polynomial
    xs * xs * (3 - 2 * xs)
  end

max_stars_repo_name: bgctw/MTKHelpers

max_issues_repo_name: bgctw/MTKHelpers

max_forks_repo_name: bgctw/MTKHelpers


text:
import Base: convert
import Base: ==

@compat abstract type InferenceType end
@compat abstract type PackedInferenceType end

@compat abstract type FunctorInferenceType <: Function end

abstract type InferenceVariable end
abstract type ConvolutionObject <: Function end

# been replaced by Functor types, but may be reused for non-numerical cases
@compat abstract type Pairwise <: InferenceType end
@compat abstract type Singleton <: InferenceType end

@compat abstract type FunctorSingleton <: FunctorInferenceType end
# @compat abstract type FunctorPartialSingleton <: FunctorInferenceType end
@compat abstract type FunctorSingletonNH <: FunctorSingleton end

@compat abstract type FunctorPairwise <: FunctorInferenceType end
@compat abstract type FunctorPairwiseMinimize <: FunctorInferenceType end
@compat abstract type FunctorPairwiseNH <: FunctorPairwise end
# @compat abstract type FunctorPairwiseNHMinimize <: FunctorPairwiseMinimize end # TODO

struct ContinuousScalar <: InferenceVariable
  dims::Int
  labels::Vector{String}
  ContinuousScalar() = new(1, String["";])
end
struct ContinuousMultivariate <:InferenceVariable
  dims::Int
  labels::Vector{String}
  ContinuousMultivariate() = new()
  ContinuousMultivariate(x) = new(x, String["";])
end

const FGG = Graphs.GenericIncidenceList{Graphs.ExVertex,Graphs.Edge{Graphs.ExVertex},Array{Graphs.ExVertex,1},Array{Array{Graphs.Edge{Graphs.ExVertex},1},1}}
const FGGdict = Graphs.GenericIncidenceList{Graphs.ExVertex,Graphs.Edge{Graphs.ExVertex},Dict{Int,Graphs.ExVertex},Dict{Int,Array{Graphs.Edge{Graphs.ExVertex},1}}}

# Condensed representation of KernelDensityEstimate, by saving points and bandwidth
mutable struct EasyMessage
  pts::Array{Float64,2}
  bws::Array{Float64,1}
  EasyMessage() = new()
  EasyMessage(a::Array{Float64,2}, b::Array{Float64,1}) = new(a,b)
  EasyMessage(p::BallTreeDensity) = new(getPoints(p), getBW(p)[:,1])
end

mutable struct FactorGraph
  g::FGGdict
  bn
  IDs::Dict{Symbol,Int}
  fIDs::Dict{Symbol,Int}
  id::Int
  nodeIDs::Array{Int,1} # TODO -- ordering seems improved to use adj permutation -- pending merge JuliaArchive/Graphs.jl/#225
  factorIDs::Array{Int,1}
  bnverts::Dict{Int,Graphs.ExVertex} # TODO -- not sure if this is still used, remove
  bnid::Int # TODO -- not sure if this is still used
  dimID::Int
  cg
  cgIDs::Dict{Int,Int} # cgIDs[exvid] = neoid
  sessionname::String
  robotname::String
  username::String
  registeredModuleFunctions::VoidUnion{Dict{Symbol, Function}}
  reference::VoidUnion{Dict{Symbol, Tuple{Symbol, Vector{Float64}}}}
  stateless::Bool
  FactorGraph() = new()
  FactorGraph(
    x1,
    x2,
    x3,
    x4,
    x5,
    x6,
    x7,
    x8,
    x9,
    x10,
    x11,
    x12,
    x13,
    x14,
    x15,
    x16,
    x17
   ) = new(
    x1,
    x2,
    x3,
    x4,
    x5,
    x6,
    x7,
    x8,
    x9,
    x10,
    x11,
    x12,
    x13,
    x14,
    x15,
    x16,
    x17,
    false )
end

"""
    $(SIGNATURES)

Construct an empty FactorGraph object with the minimum amount of information / memory populated.
"""
function emptyFactorGraph(;reference::VoidUnion{Dict{Symbol, Tuple{Symbol, Vector{Float64}}}}=nothing)
    fg = FactorGraph(Graphs.incdict(Graphs.ExVertex,is_directed=false),
                     Graphs.incdict(Graphs.ExVertex,is_directed=true),
                    #  Dict{Int,Graphs.ExVertex}(),
                    #  Dict{Int,Graphs.ExVertex}(),
                     Dict{Symbol,Int}(),
                     Dict{Symbol,Int}(),
                     0,
                     [],
                     [],
                     Dict{Int,Graphs.ExVertex}(),
                     0,
                     0,
                     nothing,
                     Dict{Int,Int}(),
                     "",
                     "",
                     "",
                     Dict{Symbol, Function}(:IncrementalInference=>IncrementalInference.getSample), # TODO likely to be removed
                     reference  ) #evalPotential
    return fg
end

mutable struct VariableNodeData
  initval::Array{Float64,2}
  initstdev::Array{Float64,2}
  val::Array{Float64,2}
  bw::Array{Float64,2}
  BayesNetOutVertIDs::Array{Int,1}
  dimIDs::Array{Int,1}
  dims::Int
  eliminated::Bool
  BayesNetVertID::Int
  separator::Array{Int,1}
  groundtruth::VoidUnion{ Dict{ Tuple{Symbol, Vector{Float64}} } } # not packed yet
  softtype
  initialized::Bool
  VariableNodeData() = new()
  function VariableNodeData(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11)
    warn("Deprecated use of VariableNodeData(11 param), use 13 parameters instead")
    new(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11, nothing, true) # TODO ensure this is initialized true is working for most cases
  end
  VariableNodeData(x1::Array{Float64,2},
                   x2::Array{Float64,2},
                   x3::Array{Float64,2},
                   x4::Array{Float64,2},
                   x5::Vector{Int},
                   x6::Vector{Int},
                   x7::Int,
                   x8::Bool,
                   x9::Int,
                   x10::Vector{Int},
                   x11::VoidUnion{ Dict{ Tuple{Symbol, Vector{Float64}} } },
                   x12,
                   x13::Bool) =
    new(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13)
end

mutable struct FactorMetadata
  factoruserdata
  variableuserdata::Union{Vector, Tuple}
  variablesmalldata::Union{Vector, Tuple}
  solvefor::Union{Symbol, Void}
  variablelist::Union{Void, Vector{Symbol}}
  dbg::Bool
  FactorMetadata() = new() # [], []
  FactorMetadata(x1, x2::Union{Vector,Tuple},x3) = new(x1, x2, x3, nothing, nothing, false)
  FactorMetadata(x1, x2::Union{Vector,Tuple},x3,x4::Symbol) = new(x1, x2, x3, x4, nothing, false)
  FactorMetadata(x1, x2::Union{Vector,Tuple},x3,x4::Symbol,x5::Vector{Symbol};dbg::Bool=false) = new(x1, x2, x3, x4, x5, dbg)
end

struct SingleThreaded
end
struct MultiThreaded
end

mutable struct ConvPerThread
  thrid_::Int
  particleidx::Int # the actual particle being solved at this moment
  factormetadata::FactorMetadata # additional data passed to user function -- optionally used by user function
  activehypo::Union{UnitRange{Int},Vector{Int}} # subsection indices to select which params should be used for this hypothesis evaluation
  p::Vector{Int} # a permutation vector for low-dimension solves (FunctorPairwise only)
  perturb::Vector{Float64} # slight numerical perturbation for degenerate solver cases such as division by zero
  X::Array{Float64,2}
  Y::Vector{Float64}
  res::Vector{Float64}
  ConvPerThread() = new()
end

function ConvPerThread(X::Array{Float64,2},
                       zDim::Int;
                       factormetadata::FactorMetadata=FactorMetadata(),
                       particleidx::Int=1,
                       activehypo= 1:length(params),
                       p=collect(1:size(X,1)),
                       perturb=zeros(zDim),
                       Y=zeros(size(X,1)),
                       res=zeros(0)  )
  #
  cpt = ConvPerThread()
  cpt.thrid_ = 0
  cpt.X = X
  cpt.factormetadata = factormetadata
  cpt.particleidx = particleidx
  cpt.activehypo = activehypo
  cpt.p = p
  cpt.perturb = perturb
  cpt.Y = Y
  cpt.res = res
  return cpt
end

mutable struct CommonConvWrapper{T} <: ConvolutionObject where {T<:FunctorInferenceType}
  ### Values consistent across all threads during approx convolution
  usrfnc!::T # user factor / function
  # general setup
  xDim::Int
  zDim::Int
  # special case settings
  specialzDim::Bool # is there a special zDim requirement -- defined by user
  partial::Bool # is this a partial constraint -- defined by user
  # multi hypothesis settings
  hypotheses::Union{Void, Distributions.Categorical} # categorical to select which hypothesis is being considered during convolugtion operation
  # values specific to one complete convolution operation
  params::Vector{Array{Float64,2}} # parameters passed to each hypothesis evaluation event on user function
  varidx::Int # which index is being solved for in params?
  measurement::Tuple # user defined measurement values for each approxConv operation
  threadmodel::Union{Type{SingleThreaded}, Type{MultiThreaded}}

  ### particular convolution computation values per particle idx (varies by thread)
  cpt::Vector{ConvPerThread}
  # varidx::Int # which index is being solved for in params?
  # factormetadata::FactorMetadata # additional data passed to user function -- optionally used by user function
  # activehypo::Union{UnitRange{Int},Vector{Int}} # subsection indices to select which params should be used for this hypothesis evaluation
  # particleidx::Int # the actual particle being solved at this moment
  # p::Vector{Int} # a permutation vector for low-dimension solves (FunctorPairwise only)
  # perturb::Vector{Float64} # slight numerical perturbation for degenerate solver cases such as division by zero
  # X::Array{Float64,2}
  # Y::Vector{Float64}
  # res::Vector{Float64}

  CommonConvWrapper{T}() where {T<:FunctorInferenceType} = new{T}()
end

function CommonConvWrapper(fnc::T,
                           X::Array{Float64,2},
                           zDim::Int,
                           params::Vector{Array{Float64,2}};
                           factormetadata::FactorMetadata=FactorMetadata(),
                           specialzDim::Bool=false,
                           partial::Bool=false,
                           hypotheses=nothing,
                           activehypo= 1:length(params),
                           varidx::Int=1,
                           measurement::Tuple=(zeros(0,1),),
                           particleidx::Int=1,
                           p=collect(1:size(X,1)),
                           perturb=zeros(zDim),
                           Y=zeros(size(X,1)),
                           xDim=size(X,1),
                           res=zeros(0),
                           threadmodel=MultiThreaded  ) where {T<:FunctorInferenceType}
  #
  ccw = CommonConvWrapper{T}()

  ccw.usrfnc! = fnc
  ccw.xDim = xDim
  ccw.zDim = zDim
  ccw.specialzDim = specialzDim
  ccw.partial = partial
  ccw.hypotheses = hypotheses
  ccw.params = params
  ccw.varidx = varidx
  ccw.threadmodel = threadmodel

  # thread specific elements
  ccw.cpt = Vector{ConvPerThread}(Threads.nthreads())
  for i in 1:Threads.nthreads()
    ccw.cpt[i] = ConvPerThread(X, zDim,
                    factormetadata=factormetadata,
                    particleidx=particleidx,
                    activehypo=activehypo,
                    p=p,
                    perturb=perturb,
                    Y=Y,
                    res=res )
  end
  # ccw.varidx = varidx
  # ccw.measurement = measurement
  # ccw.activehypo = activehypo
  # ccw.factormetadata = factormetadata
  # ccw.particleidx = particleidx
  # ccw.p = p
  # ccw.perturb = perturb
  # ccw.X = X
  # ccw.Y = Y
  # ccw.res = res

  return ccw
end

mutable struct GenericFunctionNodeData{T, S}
  fncargvID::Array{Int,1}
  eliminated::Bool
  potentialused::Bool
  edgeIDs::Array{Int,1}
  frommodule::S #Union{Symbol, AbstractString}
  fnc::T
  multihypo::String # likely to moved when GenericWrapParam is refactored
  GenericFunctionNodeData{T, S}() where {T, S} = new{T,S}()
  GenericFunctionNodeData{T, S}(x1, x2, x3, x4, x5::S, x6::T, x7::String="") where {T, S} = new{T,S}(x1, x2, x3, x4, x5, x6, x7)
  GenericFunctionNodeData(x1, x2, x3, x4, x5::S, x6::T, x7::String="") where {T, S} = new{T,S}(x1, x2, x3, x4, x5, x6, x7)
  # GenericFunctionNodeData(x1, x2, x3, x4, x5::S, x6::T, x7::String) where {T, S} = new{T,S}(x1, x2, x3, x4, x5, x6, x7)
end

# where {T <: Union{InferenceType, FunctorInferenceType}}
const FunctionNodeData{T} = GenericFunctionNodeData{T, Symbol}
FunctionNodeData(x1, x2, x3, x4, x5::Symbol, x6::T, x7::String="") where {T <: Union{FunctorInferenceType, ConvolutionObject}}= GenericFunctionNodeData{T, Symbol}(x1, x2, x3, x4, x5, x6, x7)

# where {T <: PackedInferenceType}
const PackedFunctionNodeData{T} = GenericFunctionNodeData{T, <: AbstractString}
PackedFunctionNodeData(x1, x2, x3, x4, x5::S, x6::T, x7::String="") where {T <: PackedInferenceType, S <: AbstractString} = GenericFunctionNodeData(x1, x2, x3, x4, x5, x6, x7)

###

function addGraphsVert!(fgl::FactorGraph,
            exvert::Graphs.ExVertex;
            labels::Vector{<:AbstractString}=String[])
  #
  Graphs.add_vertex!(fgl.g, exvert)
end

function getVertNode(fgl::FactorGraph, id::Int; nt::Symbol=:var, bigData::Bool=false)
  return fgl.g.vertices[id] # check equivalence between fgl.v/f[i] and fgl.g.vertices[i]
  # return nt == :var ? fgl.v[id] : fgl.f[id]
end
function getVertNode(fgl::FactorGraph, lbl::Symbol; nt::Symbol=:var, bigData::Bool=false)
  return getVertNode(fgl, (nt == :var ? fgl.IDs[lbl] : fgl.fIDs[lbl]), nt=nt, bigData=bigData)
end
getVertNode{T <: AbstractString}(fgl::FactorGraph, lbl::T; nt::Symbol=:var, bigData::Bool=false) = getVertNode(fgl, Symbol(lbl), nt=nt, bigData=bigData)

# excessive function, needs refactoring
function updateFullVertData!(fgl::FactorGraph,
    nv::Graphs.ExVertex;
    updateMAPest::Bool=false )
  #

  # not required, since we using reference -- placeholder function CloudGraphs interface
  # getVertNode(fgl, nv.index).attributes["data"] = nv.attributes["data"]
  nothing
end

function makeAddEdge!(fgl::FactorGraph, v1::Graphs.ExVertex, v2::Graphs.ExVertex; saveedgeID::Bool=true)
  edge = Graphs.make_edge(fgl.g, v1, v2)
  Graphs.add_edge!(fgl.g, edge)
  if saveedgeID push!(getData(v2).edgeIDs,edge.index) end #.attributes["data"]
  edge
end

function graphsOutNeighbors(fgl::FactorGraph, vert::Graphs.ExVertex; ready::Int=1,backendset::Int=1, needdata::Bool=false)
  Graphs.out_neighbors(vert, fgl.g)
end
function graphsOutNeighbors(fgl::FactorGraph, exVertId::Int; ready::Int=1,backendset::Int=1, needdata::Bool=false)
  graphsOutNeighbors(fgl.g, getVert(fgl,exVertId), ready=ready, backendset=backendset, needdata=needdata)
end

function graphsGetEdge(fgl::FactorGraph, id::Int)
  nothing
end

function graphsDeleteVertex!(fgl::FactorGraph, vert::Graphs.ExVertex)
  warn("graphsDeleteVertex! -- not deleting Graphs.jl vertex id=$(vert.index)")
  nothing
end

#


text:
export RAFFOutput

"""

This type defines the output file for the RAFF algorithm.

    RAFFOutput(status::Int, solution::Vector{Float64}, iter::Int,
               p::Int, f::Float64, nf::Int, nj::Int, outliers::Vector{Int})

where

  - `status`: is 1 if converged and 0 if not
  - `solution`: vector with the parameters of the model
  - `iter`: number of iterations up to convergence
  - `p`: number of trusted points
  - `f`: the residual value
  - `nf`: number of function evaluations
  - `nj`: number of Jacobian evaluations
  - `outliers`: the possible outliers detected by the method, for the
    given `p`

    RAFFOutput()

Creates a null version of output, equivalent to `RAFFOutput(0, [], -1, 0, Inf, -1, -1, [])`

    RAFFOuput(p::Int)
    RAFFOuput(sol::Vector{Float64}, p::Int)

Creates a null version of output for the given `p` and a null version
with the given solution, respectively.

"""
struct RAFFOutput

    status :: Int
    solution :: Vector{Float64}
    iter :: Int
    p :: Int
    f :: Float64
    nf :: Int
    nj :: Int
    outliers :: Vector{Int}

end

RAFFOutput(p) = RAFFOutput(0, [], -1, p, Inf, -1, -1, [])

RAFFOutput(sol, p) = RAFFOutput(0, sol, -1, p, Inf, -1, -1, [])

# Deprecated compatibility function.
RAFFOutput(status, solution, iter, p, f, outliers) = begin

    Base.depwarn("The call to `RAFFOutput` without `nf` and `nj` will be deprecated in future versions. Use the ful version `RAFFOutput(status, solution, iter, p, f, nf, nj, outliers)` instead.", :raff)

    RAFFOutput(status, solution, iter, p, f, -1, -1, outliers)

end

RAFFOutput() = RAFFOutput(0)

# Overload the == for RAFFOutput

import Base.==

==(a::RAFFOutput, b::RAFFOutput) = (
    (a.status == b.status) &&
    (a.solution == b.solution) &&
    (a.iter == b.iter) &&
    (a.p == b.p) &&
    (a.f == b.f) &&
    (a.outliers == b.outliers) &&
    (a.nf == b.nf) &&
    (a.nj == b.nj)
)

import Base.show

function show(io::IO, ro::RAFFOutput)
    print(io,"** RAFFOutput ** \n")
    print(io,"Status (.status) = $(ro.status) \n")
    print(io,"Solution (.solution) = $(ro.solution) \n")
    print(io,"Number of iterations (.iter) = $(ro.iter) \n")
    print(io,"Number of trust points (.p) = $(ro.p) \n")
    print(io,"Objective function value (.f) = $(ro.f) \n")
    print(io,"Number of function evaluations (.nf) = $(ro.nf)\n")
    print(io,"Number of Jacobian evaluations (.nj) = $(ro.nj)\n")
    print(io,"Index of outliers (.outliers) = $(ro.outliers)\n")
end


text:
using Glutils
using Test

@testset "Glutils.jl" begin
    include("../examples/runexamples.jl")
end # testset


text:
module Pythagorus

include("World.jl")

export World
export explore_concepts

include("Concept.jl")

export Concept
export find_examples_of,make_specialization

end # module


text:
module GeoPhyInv

# load all necessary packages
using Misfits
using TimerOutputs
using LinearMaps
using Ipopt
using Optim, LineSearches
using DistributedArrays
using Calculus
using ProgressMeter
using Distributed
using DistributedArrays
using SharedArrays
using Printf
using DataFrames
using SparseArrays
using Interpolations
using OrderedCollections
using CSV
using Statistics
using LinearAlgebra
using Random
using ImageFiltering
using NamedArrays
using DSP
using Test
using AxisArrays
using Distributions
using StatsBase
using RecipesBase
using FFTW

# create a timer object, used throughout this package, see TimerOutputs.jl
global const to = TimerOutput();

# include modules (note: due to dependencies, order is important!)
include("Interpolation/Interpolation.jl")
include("Utils/Utils.jl")
include("Operators.jl")
include("Smooth.jl")
include("IO.jl")
include("media/medium.jl")
include("srcwav/wavelets.jl")

# need to define supersource, source and receiver structs and export them (necessary for multiple dispatch)
struct Srcs
	n::Int
end
Srcs()=Srcs(0)
struct SSrcs
	n::Int
end
SSrcs()=SSrcs(0)
struct Recs
	n::Int
end
Recs()=Recs(0)

include("ageom/core.jl")

include("database/core.jl")

include("srcwav/core.jl")

include("Coupling.jl")
include("data/core.jl")

# Pressure and velocity fields (used for multiple dispatch)
struct P end
struct Vx end
struct Vz end

#include("Data/Data.jl")
include("Born/Born.jl")
include("fdtd/fdtd.jl")

include("fwi/fwi.jl")

include("Poisson/Poisson.jl")
include("plots.jl")

# export stuff from GeoPhyInv
export Data
export SrcWav
export update!, Medium
export ricker, ormsby
export Srcs, Recs, SSrcs
export AGeom, AGeomss
export update!, SeisForwExpt, SeisInvExpt, Fdtd, FdtdBorn, LS, LS_prior, Migr, Migr_FD

# export the Expt for Poisson
const PoissonExpt=GeoPhyInv.Poisson.ParamExpt
export PoissonExpt
mod!(a::PoissonExpt,b,c)=GeoPhyInv.Poisson.mod!(a,b,c)
mod!(a::PoissonExpt,b)=GeoPhyInv.Poisson.mod!(a,b)
mod!(a::PoissonExpt)=GeoPhyInv.Poisson.mod!(a)
operator_Born(a::PoissonExpt,b)=GeoPhyInv.Poisson.operator_Born(a,b)

end # module


text:
"""
    ICDO3

An ICD-O-3 cancer morphology code.
"""
struct ICDO3 <: AbstractICD10
  data::NTuple{6,UInt8}
end

function validICDO3input(str)
  icdO3fmt = r"^M[[:digit:]]{4}(/){0,1}[[:digit:]]{1}$"
  return occursin(icdO3fmt, str)
end

"""
    ICDO3(str::String, validateinput=false)

Create an ICDO3 code from a string.

ICDO3 codes consist of the letter `M` followed by 4 digits, a slash `/`, and a final behaviour digit.

If `validateinput=true` then the input is checked using a regex. Note this only checks the input format; it does not
check if the code actually exists.
"""
function ICDO3(str::String, validateinput = false)
  if validateinput
    validICDO3input(str) || throw(
      DomainError(
        str,
        "ICD-O-3 codes should have format `MNNNN/N` where `M` is the upper case letter 'M', `N` is a decimal digit.",
        ),
        )
  end
  return ICDO3((NTuple{5,UInt8}.(str[1:5])..., UInt8(str[end])))
end

## ICD-O-3 Base functions #####

function Base.print(io::IO, icdO3::ICDO3, punct = ICDOPTS[:punct])
  if punct
    write(io, Ref(icdO3.data[1:5]), "/", Ref(icdO3.data[6]))
  else
    write(io, Ref(icdO3))
  end
end

Base.show(io::IO, ::MIME"text/plain", icdO3::ICDO3) = print(io, icdO3, ICDOPTS[:punct])

Base.convert(::Type{T}, x::AbstractString) where {T<:ICDO3} = T(x)


text:
module PayloadUtils

import DataFrames

"""
    convert_array_to_dataframe(array::Array{Any,1})

Convert Array{Any,1} to DataFrame
"""
function convert_array_to_dataframe(array::Array{Any,1}) :: DataFrames.DataFrame
    df::DataFrames.DataFrame = DataFrames.DataFrame()

    namelist::Array = Symbol.(array[1])

    for (i, name) in enumerate(namelist)
        df[!, name] = [array[j][i] for j in 2:length(array)]
    end

    return df
end

"""
    convert_dataframe_to_array(df::DataFrames.DataFrame)

Convert DataFrame to Dict{Array}, including column names via first line of array
"""
function convert_dataframe_to_array(df::DataFrames.DataFrame) :: Array{Any,1}
    array::Array{Any, 1} = []

    push!(array, names(df))
    for row in eachrow(df)
        push!(array, [row[name] for name in names(df)])
    end

    return array
end

end


text:
"""
This example uses the
[Housing](https://archive.ics.uci.edu/ml/datasets/Housing) dataset
from the UCI Machine Learning Repository to demonstrate a linear
regression model. The dataset has housing related information for 506
neighborhoods in Boston from 1978. Each neighborhood has 14
attributes, the goal is to use the first 13, such as average number of
rooms per house, or distance to employment centers, to predict the
14’th attribute: median dollar value of the houses.

You can run the demo using `julia housing.jl`.  Use `julia housing.jl
--help` for a list of options.  The dataset will be automatically
downloaded and randomly split into training and test sets.  The
quadratic loss for the training and test sets will be printed at every
epoch and optimized parameters will be returned.

"""
module Housing
using Knet,ArgParse

function main(args=ARGS)
    s = ArgParseSettings()
    s.description="housing.jl (c) Deniz Yuret, 2016. Linear regression model for the Housing dataset from the UCI Machine Learning
Repository."
    s.exc_handler=ArgParse.debug_handler
    @add_arg_table s begin
        ("--seed"; arg_type=Int; default=-1; help="random number seed: use a nonnegative int for repeatable results")
        ("--epochs"; arg_type=Int; default=20; help="number of epochs for training")
        ("--lr"; arg_type=Float64; default=0.1; help="learning rate")
        ("--test"; arg_type=Float64; default=0.0; help="ratio of data to split for testing")
        ("--atype"; default=(gpu()>=0 ? "KnetArray" : "Array"); help="array type: Array for cpu, KnetArray for gpu")
        ("--fast"; action=:store_true; help="skip loss printing for faster run")
        ("--gcheck"; arg_type=Int; default=0; help="check N random gradients")
    end
    println(s.description)
    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)
    println("opts=",[(k,v) for (k,v) in o]...)
    o[:seed] > 0 && srand(o[:seed])
    atype = eval(parse(o[:atype]))
    w = Any[convert(atype, 0.1*randn(1,13)), 0]
    (xtrn,ytrn,xtst,ytst) = map(x->convert(atype,x), loaddata(o[:test]))
    println((:epoch,0,:trn,loss(w,xtrn,ytrn),:tst,loss(w,xtst,ytst)))
    if o[:fast]
        @time train(w, xtrn, ytrn; lr=o[:lr], epochs=o[:epochs])
        println((:epoch,o[:epochs],:trn,loss(w,xtrn,ytrn),:tst,loss(w,xtst,ytst)))
    else
        @time for epoch=1:o[:epochs]
            train(w, xtrn, ytrn; lr=o[:lr], epochs=1)
            println((:epoch,epoch,:trn,loss(w,xtrn,ytrn),:tst,loss(w,xtst,ytst)))
            if o[:gcheck] > 0
                gradcheck(loss, w, xtst, ytst; gcheck=o[:gcheck])
            end
        end
    end
    return w
end

predict(w,x)=(w[1]*x.+w[2])

loss(w,x,y)=(sum(abs2(y-predict(w,x))) / size(x,2))

lossgradient = grad(loss)

function train(w, x, y; lr=.1, epochs=20)
    for epoch=1:epochs
        g = lossgradient(w, x, y)
        for i in 1:length(w)
            w[i] -= lr * g[i]
        end
    end
    return w
end

function loaddata(test=0.0)
    url = "https://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data"
    file=joinpath(Knet.datapath,"housing.data")
    if !isfile(file)
        info("Downloading $url to $file")
        download(url, file)
    end
    data = readdlm(file)'
    @show size(data) # (14,506)
    x = data[1:13,:]
    y = data[14:14,:]
    x = (x .- mean(x,2)) ./ std(x,2) # Data normalization
    if test == 0
        xtrn = xtst = x
        ytrn = ytst = y
    else
        r = randperm(size(x,2))          # trn/tst split
        n = round(Int, (1-test) * size(x,2))
        xtrn=x[:,r[1:n]]
        ytrn=y[:,r[1:n]]
        xtst=x[:,r[n+1:end]]
        ytst=y[:,r[n+1:end]]
    end
    (xtrn, ytrn, xtst, ytst)
end

# This allows both non-interactive (shell command) and interactive calls like:
# $ julia housing.jl --epochs 10
# julia> Housing.main("--epochs 10")
!isinteractive() && !isdefined(Core.Main,:load_only) && main(ARGS)

end # module Housing

# SAMPLE RUN 65f57ff+ Wed Sep 14 10:02:30 EEST 2016
#
# housing.jl (c) Deniz Yuret, 2016. Linear regression model for the Housing dataset from the UCI Machine Learning
# Repository.
# opts=(:seed,-1)(:epochs,20)(:lr,0.1)(:atype,"KnetArray")(:gcheck,0)(:fast,true)
# size(data) = (14,506)
# (:epoch,0,:trn,577.8060720494859,:tst,656.5939685291422)
#   0.018364 seconds (7.15 k allocations: 347.266 KB)
# (:epoch,20,:trn,20.908099852265146,:tst,33.9047224887537)


text:
# Copyright 2017 Massachusetts Institute of Technology. See LICENSE file for details.

__precompile__(true)
module Defer
import Base.GC

export push_scope!, pop_scope!, scope, scope_nogc, @scope, defer, @defer, defer_call, @!

const scopes = Any[Any[]]
const ExceptionWrapper = Union{Some{Any}, Nothing}

function push_scope!()
  push!(scopes, Any[])
  length(scopes)
end

function pop_scope!(i::Int, e::ExceptionWrapper=nothing)
  if i == length(scopes)
    pop_scope!(e)
  else
    @assert i > 0
    warn("Popping scope $(length(scopes)), expected $i")
    if length(scopes) < i
      if e != nothing
        throw(something(e))
      end
    else
      while length(scopes) > i
        pop_scope!()
      end
      pop_scope!(e)
    end
  end
end

function pop_scope!(e::ExceptionWrapper=nothing)
  exceptions = e==nothing ? Any[] : Any[something(e)]
  this_scope = pop!(scopes)
  if isempty(scopes)
    push!(scopes, Any[])
  end
  for fin in reverse(this_scope)
    try
      fin()
    catch e
      push!(exceptions, e)
    end
  end
  empty!(this_scope)
  if !isempty(exceptions)
    if length(exceptions) == 1
      throw(exceptions[1])
    else
      e = CompositeException()
      append!(e.exceptions, exceptions)
      throw(e)
    end
  end
  nothing
end

function scope(f)
  sc = push_scope!()
  ex = Ref{Union{Some{Any},Nothing}}(nothing)
  try
    f()
  catch e
    ex[] = Some(e)
  finally
    pop_scope!(sc, Some(e))
  end
end

function scope_nogc(f)
  gc_enabled = GC.enable(false)
  sc = push_scope!()
  ex = Ref{Union{Some{Any},Nothing}}(nothing)
  try
    f()
  catch e
    ex[] = Some(e)
  finally
    pop_scope!(sc, ex[])
    if gc_enabled
      GC.enable(true)
    end
  end
end

_scope(code::Expr) = quote
  sc = push_scope!()
  ex = Ref{Union{Some{Any},Nothing}}(nothing)
  try
    $code
  catch e
    ex[] = Some(e)
  finally
    pop_scope!(sc, ex[])
  end
end

macro scope(code)
  if !isa(code, Expr) return esc(code) end
  if code.head == :let
    flet(x::Expr) = :($(esc(x.args[1])) = defer_call($(esc(x.args[2]))))
    newcode = Expr(:let, code.args[1], esc(code.args[2]))
    if code.args[1].head == :(=)
      newcode.args[1] = flet(code.args[1])
    else
      for i in 1:length(code.args[1].args)
        @assert code.args[1].args[i].head == :(=)
        newcode.args[1].args[i] = flet(code.args[1].args[i])
      end
    end
    return _scope(newcode)
  elseif code.head == :(=)
    return Expr(:(=), esc(code.args[1]), _scope(esc(code.args[2])))
  elseif code.head == :function
    return Expr(:function, esc(code.args[1]), _scope(esc(code.args[2])))
   else
    return _scope(esc(code))
  end
end

function defer(fin, sc::Integer=0)
  push!(scopes[sc > 0 ? sc : length(scopes) + sc], fin)
  nothing
end

function defer_call(x, f::Function=close, sc::Integer=0)
  defer(()->f(x), sc)
  x
end

macro defer(code)
  quote
    defer(()->$(esc(code)))
  end
end

macro !(value)
  quote
    defer_call($(esc(value)))
  end
end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess
using JuliaFEM.Testing

@testset "test that interface transfers constant field without error" begin
    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_3d.med"
    mesh = aster_read_mesh(meshfile)

    upper = Problem(Heat, "upper", 1)
    upper.elements = create_elements(mesh, "UPPER")
    update!(upper, "temperature thermal conductivity", 1.0)
    lower = Problem(Heat, "lower", 1)
    lower.elements = create_elements(mesh, "LOWER")
    update!(lower, "temperature thermal conductivity", 1.0)

    bc_upper = Problem(Dirichlet, "upper boundary", 1, "temperature")
    bc_upper.elements = create_elements(mesh, "UPPER_TOP")
    update!(bc_upper, "temperature 1", 0.0)

    bc_lower = Problem(Dirichlet, "lower boundary", 1, "temperature")
    bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM")
    update!(bc_lower, "temperature 1", 1.0)

    interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
    interface_slave_elements = create_elements(mesh, "LOWER_TOP")
    interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
    update!(interface_slave_elements, "master elements", interface_master_elements)
    interface.elements = [interface_master_elements; interface_slave_elements]

    solver = LinearSolver(upper, lower, bc_upper, bc_lower, interface)
    solver()

    node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0)
    T = [t[1] for t in temperature]
    minT = minimum(T)
    maxT = maximum(T)
    info("minT = $minT, maxT = $maxT")
    @test isapprox(minT, 0.5)
    @test isapprox(maxT, 0.5)
end


text:
# Atmos RHS

function ∑tendencies_atm!(dY, Y, (parameters, T_sfc), t)
    # Heat diffusion:
    # ∂_t ρ =  ∇ (μ ∇ ρ - w ρ)
    # ∂_t ρθ =  ∇ (μ ∇ ρθ - w ρθ)
    # ∂_t w =  ∇ (μ ∇ w - w w)

    # where
    # ∂_t T = n \cdot F   at z = zmin_atm
    # ∂_t T = 0           at z = zmax_atm
    # We also use this model to accumulate fluxes
    # ∂_t ϕ_bottom = n \cdot F

    # μ = FT(0.0001) # diffusion coefficient

    # T = u.x[1]

    # F_sfc = - calculate_flux( T_sfc[1], parent(T)[1] )

    # # set BCs
    # bcs_bottom = Operators.SetValue(F_sfc) # struct w bottom BCs
    # bcs_top = Operators.SetValue(FT(280.0))

    # gradc2f = Operators.GradientC2F(top = bcs_top) # gradient struct w BCs
    # gradf2c = Operators.GradientF2C(bottom = bcs_bottom)
    # @. du.x[1] = gradf2c( μ * gradc2f(T))
    # du.x[2] .= - F_sfc[1]

    UnPack.@unpack Cd, f, ν, ug, vg, C_p, MSLP, R_d, R_m, C_v, grav = parameters

    (Yc, Yf, F_sfc) = Y.x
    (dYc, dYf, dF_sfc) = dY.x
    UnPack.@unpack ρ, u, v, ρθ = Yc
    UnPack.@unpack w = Yf
    dρ = dYc.ρ
    du = dYc.u
    dv = dYc.v
    dρθ = dYc.ρθ
    dw = dYf.w

    # auxiliary calculations
    u_1 = parent(u)[1]
    v_1 = parent(v)[1]
    u_wind = sqrt(u_1^2 + v_1^2)

    # density (centers)
    gradc2f = Operators.GradientC2F()
    gradf2c = Operators.GradientF2C(bottom = Operators.SetValue(0.0), top = Operators.SetValue(0.0))

    If = Operators.InterpolateC2F(bottom = Operators.Extrapolate(), top = Operators.Extrapolate())
    @. dρ = gradf2c(-w * If(ρ)) # Eq. 4.11

    # potential temperature (centers)
    gradc2f = Operators.GradientC2F()
    gradf2c = Operators.GradientF2C(bottom = Operators.SetValue(0.0), top = Operators.SetValue(0.0)) # Eq. 4.20, 4.21

    @. dρθ = gradf2c(-w * If(ρθ) + ν * gradc2f(ρθ / ρ)) # Eq. 4.12

    # u velocity (centers)
    gradc2f = Operators.GradientC2F(top = Operators.SetValue(ug)) # Eq. 4.18
    gradf2c = Operators.GradientF2C(bottom = Operators.SetValue(Cd * u_wind * u_1)) # Eq. 4.16

    A = Operators.AdvectionC2C(bottom = Operators.SetValue(0.0), top = Operators.SetValue(0.0))
    @. du = gradf2c(ν * gradc2f(u)) + f * (v - vg) - A(w, u) # Eq. 4.8

    # v velocity (centers)
    gradc2f = Operators.GradientC2F(top = Operators.SetValue(vg)) # Eq. 4.18
    gradf2c = Operators.GradientF2C(bottom = Operators.SetValue(Cd * u_wind * v_1)) # Eq. 4.16

    A = Operators.AdvectionC2C(bottom = Operators.SetValue(0.0), top = Operators.SetValue(0.0))
    @. dv = gradf2c(ν * gradc2f(v)) - f * (u - ug) - A(w, v) # Eq. 4.9

    # w velocity (faces)
    gradc2f = Operators.GradientC2F()
    gradf2c = Operators.GradientF2C(bottom = Operators.SetValue(0.0), top = Operators.SetValue(0.0))

    B = Operators.SetBoundaryOperator(bottom = Operators.SetValue(0.0), top = Operators.SetValue(0.0))
    If = Operators.InterpolateC2F(bottom = Operators.Extrapolate(), top = Operators.Extrapolate())
    Π(ρθ) = C_p .* (R_d .* ρθ ./ MSLP) .^ (R_m ./ C_v)
    #dw = ρθ .* 0.0
    @. dw = B(-(If(ρθ / ρ) * gradc2f(Π(ρθ))) - grav + gradc2f(ν * gradf2c(w)) - w * If(gradf2c(w))) # Eq. 4.10

    return dY

end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
@testset "Basics" begin

    @testset "constructors" begin
        A = round.(Int, max.(0.0, 5 * rand(6, 10) .- 1))
        f = Polynomial(rand(10), A)
        @test typeof(f) <: Polynomial{Float64}

        @test_throws AssertionError Polynomial(rand(9), A)

        @polyvar x y
        f2 = Polynomial(2x^2+4y^2+3x*y+1)
        @test exponents(f2) == [0 2 1 0; 0 0 1 2]
        @test nvariables(f2) == 2
        @test coefficients(f2) == [1, 2, 3, 4]
        @test coefficienttype(f2) == Int64
        @test permutation(f2) == [4, 1, 2, 3]
        f2_2 = Polynomial(2x^2+4y^2+3x*y+1)
        @test f2 == f2_2

        @polyvar x y a b
        f = Polynomial(2x+b*y+a, parameters=[a, b])
        @test f isa Polynomial
        @test nvariables(f) == 2
        @test parameters(f) == [:a, :b]

        @test_throws ErrorException Polynomial(2x+b*y+a, variables=[x], parameters=[a, b])
    end

    @testset "evaluation" begin
        @polyvar x y
        f2 = 2x^2+4y^2+3x*y^4+1
        g = Polynomial(f2)
        w = rand(2)

        @test abs(SP.evaluate(g, w) - f2(x => w[1], y => w[2])) < 1e-13

        f2_x = MP.differentiate(2x^2+4y^2+3x*y^4+1, x)
        f2_y = MP.differentiate(2x^2+4y^2+3x*y^4+1, y)
        @test norm(SP.gradient(g, w) - [f2_x(x => w[1], y => w[2]), f2_y(x => w[1], y => w[2])]) < 1e-13

        u = zeros(2)
        SP.gradient!(u, g, w)
        @test u ≈ SP.gradient(g, w)

        @test all(SP.evaluate_and_gradient(g, w) .≈ (g(w), u))

        @test all(SP.evaluate_and_gradient!(u, g, w) .≈ g(w))
        @test u ≈ SP.gradient(g, w)
    end

    @testset "Evaluation with parameters" begin
        @polyvar x y z a b
        f = 2x^2+4a*y^2+(a^2+b^4+a+b+1)*3x*y^4+1+a*z+b-2x*b^6
        g = Polynomial(f, parameters=[a, b])

        w = [5, 3, -4]
        p = [-2, 7]
        f_wp = f([x, y, z] => w, [a, b] => p)
        ∇_wp = map(fi -> fi([x, y, z] => w, [a, b] => p), MP.differentiate(f, [x, y, z]))
        ∇p_wp = map(fi -> fi([x, y, z] => w, [a, b] => p), MP.differentiate(f, [a, b]))

        @test g(w, p) == f_wp
        @test gradient(g, w,p) == ∇_wp
        @test evaluate_and_gradient(g, w, p) == (f_wp, ∇_wp)
        @test differentiate_parameters(g, w, p) == ∇p_wp

        u = zeros(Int, 3)
        gradient!(u, g, w, p)
        @test u == ∇_wp
        u .= 0
        @test evaluate_and_gradient!(u, g, w, p) == f_wp
        @test u == ∇_wp

        v = zeros(Int, 2)
        differentiate_parameters!(v, g, w, p)
        @test v == ∇p_wp

        f = Polynomial(a + b + x + y, parameters=[a, b])
        @test differentiate_parameters(f, [0, 0], [1, 1]) == [1, 1]

        f = Polynomial(x*y + x + y, parameters=[a, b])
        @test differentiate_parameters(f, [0, 0], [1, 1]) == [0, 0]
    end

    @testset "show" begin
        @polyvar x[0:9]
        @test string(Polynomial(x[1])) == "x₀"
        @test string(Polynomial(sum((-1)^i * x[i]^i for i=1:length(x)))) ==
            "-x₀ + x₁² - x₂³ + x₃⁴ - x₄⁵ + x₅⁶ - x₆⁷ + x₇⁸ - x₈⁹ + x₉¹⁰"
        @test string(Polynomial(2im*x[1] - x[2])) == "(0 + 2im)x₀ - x₁"
        @test sprint(show, Polynomial((2+0im)*x[1] - x[2])) == "2x₀ - x₁"
        @test_nowarn string(typeof(Polynomial(x[1]^2+x[4]^2)))

        @polyvar x a b
        @test sprint(show, Polynomial(x+a^2+b, parameters=[a, b])) == "a² + b + x"
        @test sprint(show, Polynomial(zero(x))) == "0"
    end

    @testset "helpers" begin
        x = rand()
        z = rand(Complex{Float64})
        for k = 4:15
            @test abs(SP.pow(x, k) - x^k) < 1e-13
            @test abs(SP.pow(z, k) - z^k) < 1e-13
        end
    end

    @testset "system constructor" begin
        @polyvar x y
        f1 = x^2+y^2
        f2 = 2x^2+4y^2+3x*y^4+1
        g1 = Polynomial(f1)
        g2 = Polynomial(f2)
        @test_deprecated SP.system(g1, g2)
        @test PolynomialSystem(g1, g2) isa PolynomialSystem{2, 2}
        @test PolynomialSystem(g1, g2, g2) isa PolynomialSystem{3, 2}
        @test PolynomialSystem([f1, f2, y, x]) isa PolynomialSystem{4, 2}
        @test length(SP.PolynomialSystem([f1, f2, y, x])) == 4

        @test_nowarn sprint(show, PolynomialSystem([f1, f2, y, x]))
    end

    @testset "System evaluation and Jacobian" begin
        @polyvar x y
        f1 = x^2+y^2
        f2 = 2x^2+4y^2+3x*y^4+1
        g1 = Polynomial(f1)
        g2 = Polynomial(f2)

        G = PolynomialSystem(g1, g2)

        w = rand(2)
        @test [evaluate(g1, w), evaluate(g2, w)] ≈ evaluate(G, w)
        @test [g1(w), g2(w)] ≈ evaluate(G, w)

        w = SVector{2}(w)
        @test evaluate(G, w) isa SVector{2}
        @test [evaluate(g1, w), evaluate(g2, w)] ≈ evaluate(G, w)

        @test jacobian(G, [2, 3]) == [4 6; 251 672]
        U = zeros(Int, 2, 2)
        @test jacobian!(U, G, [2, 3]) == [4 6; 251 672] == U
        u = zeros(Int, 2)
        U .= 0
        evaluate_and_jacobian!(u, U, G, [2, 3])
        @test U == jacobian(G, [2, 3])
        @test u == G([2, 3])

        @polyvar x y a b

        F = PolynomialSystem(x^2+x*y+a+a*x, x+y+b, parameters=[a, b])

        @test F([0,0], [1, 2]) == [1, 2]
        @test jacobian(F, [0, 0], [3, 2]) == [3 0; 1 1]

        @test parameters(F) == [:a, :b]
        @test nparameters(F) == 2
        @test variables(F) == [:x, :y]
        @test nvariables(F) == 2
        @test F isa PolynomialSystem{2, 2, 2}

        @polyvar x y z a b
        f = 2x^2+4a*y^2+(a^2+b^4+a+b+1)*3x*y^4+1+a*z+b-2x*b^6
        g = (x^2+z)*(x+b*y+a)

        F = PolynomialSystem(f, g, parameters=[a, b])
        differentiate_parameters(F, [2, 3, -3], [4, -2]) == [4407 -14297; 1 3]
        differentiate_parameters(F, [2, 3, -3], @SVector([4, -2])) == [4407 -14297; 1 3]
        @test differentiate_parameters(F, [2, 3, -3], @SVector([4, -2])) isa SMatrix{2, 2, Int, 4}
        U = zeros(Int, 2, 2)
        differentiate_parameters!(U, F, [2, 3, -3], [4, -2])
        @test U == [4407 -14297; 1 3]
    end

    @testset "Hessian" begin
        @polyvar x y z

        w = [2, 3, 5]

        f = x^4*y+x*y*z+y+x+z
        g = Polynomial(f)
        grad_g = Polynomial.(MP.differentiate(f, [x, y, z]))
        hess_w = vcat(map(grad_g) do g_i
            gradient(g_i, w)'
        end...)

        @test hessian(g, w) == hess_w
        @test gradient_and_hessian(g, w) == (gradient(g, w), hess_w)

        U = zeros(Int, 3, 3)
        @test hessian!(U, g, w) == hess_w

        p = [x^4*y+x*y*z+y+x+z, (x^3+y^2)*z]
        P = PolynomialSystem(p)

        u = zeros(2, 3)
        U = zeros(2, 3, 3)

        hessian!(U, P, w)
        @test U[1,:,:] == hess_w
        jacobian_and_hessian!(u, U, P, w)
        @test u[1,:] == gradient(g, w)
        @test U[1,:,:] == hess_w
    end

    @testset "foreach" begin
        @polyvar x y
        g = [Polynomial(x^2+y^2), Polynomial(2x^2+4y^2+3x*y^4+1)]
        G = PolynomialSystem(g...)
        i = 1
        foreach(G) do gi
            @test exponents(gi) == exponents(g[i])
            i += 1
        end
    end

    @testset "scale coefficients" begin
        @polyvar x y
        f = Polynomial(x^2+y^2)
        scale_coefficients!(f, 2)
        @test coefficients(f) == [2, 2]

        g1 = Polynomial(x^2+y^2)
        g2 = Polynomial(2x^2+4y^2+3x*y^4+1)
        G = PolynomialSystem(g1, g2)
        w = rand(2)
        x1 = evaluate(G, w)
        scale_coefficients!(G, [-2, 3])
        x2 = evaluate(G, w)
        @test x2 ≈ (-2, 3) .* x1
    end
end


text:
"""
    require(package::Symbol; [fun_name]::String="", [explanation]::String="")

Helper method to check for optional packages and print an error message.

### Input

- `package`     -- symbol (the package name)
- `fun_name`    -- (optional; default: `""`) name of the function that requires
                   the package
- `explanation` -- (optional; default: `""`) additional explanation in the error
                   message

### Output

If the package is loaded, this function has no effect.
Otherwise it prints an error message.

### Algorithm

This function uses `@assert` and hence loses its ability to print an error
message if assertions are deactivated.
"""
function require(package::Symbol; fun_name::String="", explanation::String="")
    check = isdefined(@__MODULE__, package)
    @assert check "package '$package' not loaded" *
        (fun_name == "" ? "" :
            " (it is required for executing `$fun_name`" *
            (explanation == "" ? "" : " " * explanation) * ")")
end

"""
    require(packages::Vector{Symbol}; [fun_name]::String="",
                                      [explanation]::String="")

Helper method to check for optional packages and print an error message.

### Input

- `packages`    -- list of symbols (the package names)
- `fun_name`    -- (optional; default: `""`) name of the function that requires
                   the package
- `explanation` -- (optional; default: `""`) additional explanation in the error
                   message

### Output

If at least one of the packages is loaded, this function has no effect.
Otherwise it prints an error message.

### Algorithm

This function uses `@assert` and hence loses its ability to print an error
message if assertions are deactivated.
"""
function require(packages::Vector{Symbol}; fun_name::String="",
                                           explanation::String="")
    check = any(isdefined(@__MODULE__, package) for package in packages)
    @assert check "no package from '$packages' loaded" *
        (fun_name == "" ? "" :
            " (one of them is required for executing `$fun_name`" *
            (explanation == "" ? "" : " " * explanation) * ")")
end


text:
"""
    f, p = logsmooth(X, inBin, nbin=8, n=3)

Smooth the fft, `X`, and convert it to dB. Use `nbin` bins from 0 to 3×`inBin`,
thereafter increase bin sizes by a factor of 1.1, staying less than 2¹⁰.
For the `n` sets of bins `inBin`+[0:2], 2×`inBin`+[0:2], ... , `n`×`inBin`+[0:2],
don't do averaging. This way, the noise BW and the scaling of the tone and its
harmonics are unchanged. Unfortunately, harmonics above the nth appear smaller
than they really are because their energy is averaged over many bins.
"""
function logsmooth(X, inBin, nbin=8, n=3)
    N = length(X)
    N2 = N÷2
    f1 = rem(inBin-1, nbin) + 1

    startbin = [f1:nbin:inBin-1; inBin:inBin+2]
    for i = 1:n
        startbin = [startbin; (startbin[end]+1):nbin:(i*inBin-1); i*inBin .+ (0:2)]
    end
    m = startbin[length(startbin)] + nbin
    while m < N2
        startbin = [startbin; m]
        nbin = min(nbin*1.1, 2^10)
        m = round(Int, m .+ nbin)
    end
    stopbin = [startbin[2:length(startbin)] .- 1; N2]

    f = ((startbin+stopbin)/2 .- 1)/N
    p = zeros(size(f))
    for i = 1:length(f)
        p[i] = dbp(norm(X[startbin[i]:stopbin[i]])^2 / (stopbin[i] - startbin[i] + 1))
    end

    return f, p
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
# Xpress environment and other supporting facilities
"""
    Env()

Initialize Xpress environment.
Need to do ir once and only once.
"""
function Env()
    ret = @xprs_ccall(init, Cint, (Ptr{Cchar},), C_NULL)
    if ret != 0
        if ret == :number
            error("Invalid Xpress license")
        else
            error("Failed to create environment (error $ret).")
        end
    end
    # finalizer(env, free_env)  ## temporary disable: which tends to sometimes caused warnings
end

function is_valid()
    #env.ptr_env == 1
    error("is_valid does not apply to xpress models")
end

"""
    free_env()

Finalize Xpress environment.
"""
function free_env()
    ret = @xprs_ccall(free, Cint, ())
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id: HairpinIncr.mso								$
#*------------------------------------------------------------------
type Hairpin_PressureDrop
	Hairpin_PressureDrop()=begin
		N=outers.N
		Npoints=outers.Npoints
		new(
			length(Dict{Symbol,Any}(
				:Brief=>"Hydraulic Diameter of Pipe for Pressure Drop",
				:Lower=>1e-6
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of zones",
				:Default=>2
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of incremental points",
				:Default=>3
			)),
			fill(pressure(Dict{Symbol,Any}(
				:Brief=>"Incremental Local  Pressure",
				:Default=>1,
				:Lower=>1e-10,
				:Upper=>2e4,
				:DisplayUnit=>"kPa"
			)),(Npoints)),
			fill(press_delta(Dict{Symbol,Any}(
				:Brief=>"Incremental Pressure Drop for friction",
				:Default=>1e-3,
				:Lower=>0,
				:DisplayUnit=>"kPa",
				:Symbol=>"\\Delta P_{fric}"
			)),(Npoints)),
			fill(fricfactor(Dict{Symbol,Any}(
				:Brief=>"Incremental Friction Factor",
				:Default=>0.05,
				:Lower=>1e-10,
				:Upper=>2000
			)),(N)),
			fill(positive(Dict{Symbol,Any}(
				:Brief=>"Incremental Reynolds Number",
				:Default=>100,
				:Lower=>1
			)),(N)),
			[:Dh,:N,:Npoints,],
			[:Plocal,:Pd_fric,:fi,:Re,]
		)
	end
	Dh::length
	N::DanaInteger
	Npoints::DanaInteger
	Plocal::Array{pressure}
	Pd_fric::Array{press_delta}
	fi::Array{fricfactor}
	Re::Array{positive}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Hairpin_PressureDrop
function atributes(in::Hairpin_PressureDrop,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="to be documented"
	fields[:Info]="to be documented"
	drive!(fields,_)
	return fields
end
Hairpin_PressureDrop(_::Dict{Symbol,Any})=begin
	newModel=Hairpin_PressureDrop()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the ideal MHD equations
gamma = 5/3
equations = IdealGlmMhdEquations1D(gamma)

"""
    initial_condition_ryujones_shock_tube(x, t, equations::IdealGlmMhdEquations1D)

Ryu and Jones shock tube test case for one dimensional ideal MHD equations. Contains
fast shocks, slow shocks, and rational discontinuities that propagate on either side
of the contact discontinuity. Exercises the scheme to capture all 7 types of waves
present in the one dimensional MHD equations. It is the second test from Section 4 of
- Ryu and Jones (1995)
  Numerical Magnetohydrodynamics in Astrophysics: Algorithm and Tests
  for One-Dimensional Flow
  [DOI: 10.1086/175437](https://doi.org/10.1086/175437)
!!! note
    This paper has a typo in the initial conditions. Their variable `E` should be `p`.
"""
function initial_condition_ryujones_shock_tube(x, t, equations::IdealGlmMhdEquations1D)
  # domain must be set to [0, 1], γ = 5/3, final time = 0.2
  rho = x[1] <= 0.5 ? 1.08 : 1.0
  v1 = x[1] <= 0.5 ? 1.2 : 0.0
  v2 = x[1] <= 0.5 ? 0.01 : 0.0
  v3 = x[1] <= 0.5 ? 0.5 : 0.0
  p = x[1] <= 0.5 ? 0.95 : 1.0
  inv_sqrt4pi = 1.0 / sqrt(4 * pi)
  B1 = 2 * inv_sqrt4pi
  B2 = x[1] <= 0.5 ? 3.6 * inv_sqrt4pi : 4.0 * inv_sqrt4pi
  B3 = B1

  return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations)
end
initial_condition = initial_condition_ryujones_shock_tube

boundary_conditions = BoundaryConditionDirichlet(initial_condition)

surface_flux = flux_hll
volume_flux  = flux_hindenlang_gassner
basis = LobattoLegendreBasis(3)

indicator_sc = IndicatorHennemannGassner(equations, basis,
                                         alpha_max=0.5,
                                         alpha_min=0.001,
                                         alpha_smooth=true,
                                         variable=Trixi.density)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)
solver = DGSEM(basis, surface_flux, volume_integral)

coordinates_min = 0.0
coordinates_max = 1.0
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=7,
                n_cells_max=10_000,
                periodicity=false)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 0.2)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=0.8)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
# # Internal wave example
#
# In this example, we initialize an internal wave packet in two-dimensions
# and watch it propagate.

using Oceananigans, Oceananigans.Grids, Plots, Printf

# ## Numerical, domain, and internal wave parameters
#
# First, we pick a resolution and domain size,

Nx = 128 # resolution
Lx = 2π  # domain extent
nothing # hide

# We set up an internal wave with the pressure field
#
# $ p(x, y, z, t) = a(x, z) \, \cos(kx + mz - ω t) $.
#
# where $m$ is the vertical wavenumber, $k$ is the horizontal wavenumber,
# $ω$ is the wave frequncy, and $a(x, z)$ is a Gaussian envelope.

## Non-dimensional internal wave parameters
m = 16      # vertical wavenumber
k = 8       # horizontal wavenumber
N = 1       # buoyancy frequency
f = 0.2     # inertial frequency
nothing # hide

# ## A Gaussian wavepacket
#
# Next, we set up an initial condition corresponding to a propagating
# wave packet with a Gaussian envelope. The internal wave dispersion relation yields

ω² = (N^2 * k^2 + f^2 * m^2) / (k^2 + m^2)

## and thus
ω = sqrt(ω²)
nothing # hide

# The internal wave polarization relations follow from the linearized
# Boussinesq equations,

U = k * ω   / (ω^2 - f^2)
V = k * f   / (ω^2 - f^2)
W = m * ω   / (ω^2 - N^2)
B = m * N^2 / (ω^2 - N^2)
nothing # hide

# Finally, we set-up a small-amplitude, Gaussian envelope for the wave packet

## Some Gaussian parameters
A, x₀, z₀, δ = 1e-9, Lx/2, -Lx/2, Lx/15

## A Gaussian envelope
a(x, z) = A * exp( -( (x - x₀)^2 + (z - z₀)^2 ) / 2δ^2 )
nothing # hide

# Create initial condition functions
u₀(x, y, z) = a(x, z) * U * cos(k*x + m*z)
v₀(x, y, z) = a(x, z) * V * sin(k*x + m*z)
w₀(x, y, z) = a(x, z) * W * cos(k*x + m*z)
b₀(x, y, z) = a(x, z) * B * sin(k*x + m*z) + N^2 * z
nothing # hide

# We are now ready to instantiate our model on a uniform grid.
# We give the model a constant rotation rate with background vorticity `f`,
# use temperature as a buoyancy tracer, and use a small constant viscosity
# and diffusivity to stabilize the model.

model = IncompressibleModel(    grid = RegularCartesianGrid(size=(Nx, 1, Nx), extent=(Lx, Lx, Lx)),
                             closure = IsotropicDiffusivity(ν=1e-6, κ=1e-6),
                            coriolis = FPlane(f=f),
                             tracers = :b,
                            buoyancy = BuoyancyTracer())
nothing # hide

# We initialize the velocity and buoyancy fields
# with our internal wave initial condition.

set!(model, u=u₀, v=v₀, w=w₀, b=b₀)

# ## A wave packet on the loose
#
# Finally, we release the packet and watch it go!

simulation = Simulation(model, Δt = 0.001 * 2π/ω, stop_iteration = 0,
                        iteration_interval = 20)

anim = @animate for i=0:100
    x, z = xnodes(Cell, model.grid)[:], znodes(Face, model.grid)[:]
    w = model.velocities.w

    contourf(x, z, w.data[1:Nx, 1, 1:Nx+1]',
                   title = @sprintf("ωt = %.2f", ω * model.clock.time),
                  levels = range(-1e-8, stop=1e-8, length=10),
                   clims = (-1e-8, 1e-8),
                  xlabel = "x",
                  ylabel = "z",
                   xlims = (0, Lx),
                   ylims = (-Lx, 0),
               linewidth = 0,
                   color = :balance,
                  legend = false,
             aspectratio = :equal)

    simulation.stop_iteration += 20
    run!(simulation)
end

mp4(anim, "internal_wave.mp4", fps = 15) # hide


text:
using Plots
using Colors
pyplot()
include("N_plots.jl")

# Model and Objective
N = 51
model, obj = Dynamics.dubinscar_parallelpark
obj_uncon = UnconstrainedObjective(obj)

opts = SolverOptions()
opts.verbose = false
opts.cost_tolerance = 1e-6
opts.cost_tolerance_intermediate = 1e-5
opts.constraint_tolerance = 1e-4
opts.resolve_feasible = false
opts.outer_loop_update_type = :default
opts.penalty_scaling = 200
opts.penalty_initial = .1
opts.R_infeasible = 20
dircol_options = Dict("tol"=>opts.cost_tolerance,"constr_viol_tol"=>opts.constraint_tolerance)

#### UNCONSTRAINED #####
solver = Solver(model, obj_uncon, N=N, opts=opts)
n,m,N = get_sizes(solver)
U0 = ones(m,N)
X0 = line_trajectory(solver)
X0_rollout = rollout(solver,U0)

res_i, stats_i = solve(solver,U0)
stats_i["iterations"]
stats_i["runtime"]
evals(solver,:f) / stats_i["iterations"]

res_p, stats_p = solve_dircol(solver,X0_rollout,U0)
stats_p["iterations"]
evals(solver,:f)/stats_p["iterations"]
stats_p["runtime"]

Ns = [101,201,301,401]
group = "parallelpark/unconstrained"
run_step_size_comparison(model, obj_uncon, U0, group, Ns, opts=opts, integrations=[:rk3,:ipopt],benchmark=true)
plot_stat("runtime",group,legend=:bottomright,["rk3","ipopt"],title="Unconstrained Parallel Park")
plot_stat("iterations",group,legend=:bottom,["rk3","ipopt"],title="Unconstrained Parallel Park")
plot_stat("error",group,yscale=:log10,legend=:right,["rk3","ipopt"],title="Unconstrained Parallel Park")
plot_stat("std",group,yscale=:log10,legend=:right,["rk3","ipopt"],title="Unconstrained Parallel Park")

#### CONSTRAINED ####
solver = Solver(model, obj, N=N, opts=opts)
n,m,N = get_sizes(solver)
U0 = ones(m,N)
X0 = line_trajectory(solver)
X0_rollout = rollout(solver,U0)

res_i, stats_i = solve(solver,U0)
res_p, stats_p = solve_dircol(solver,X0_rollout,U0)

stats_i["iterations"]
stats_i["runtime"]
evals(solver,:f) / stats_i["iterations"]

stats_p["iterations"]
evals(solver,:f)/stats_p["iterations"]
stats_p["runtime"]

Ns = [101,201,301,401]
group = "parallelpark/constrained"
run_step_size_comparison(model, obj, U0, group, Ns, opts=opts, integrations=[:rk3,:ipopt],benchmark=true)
plot_stat("runtime",group,legend=:bottomright,["rk3","ipopt",""],title="Constrained Parallel Park",color=[:blue :darkorange2])
savefig(joinpath(IMAGE_DIR,"ppark_runtime.eps"))
plot_stat("iterations",group,legend=:bottom,["rk3","ipopt"],title="Constrained Parallel Park",color=[:blue :darkorange2])
savefig(joinpath(IMAGE_DIR,"ppark_iterations.eps"))
plot_stat("error",group,yscale=:log10,legend=:right,["rk3","ipopt"],title="Constrained Parallel Park")
Plots.eps(joinpath(IMAGE_DIR,"ppark_runtime"))

# Constraint vs time (Newton tail)
# using TrajectoryOptimization: gen_usrfun_newton, NewtonVars, gen_newton_functions, newton_projection
# t_start = time_ns()
# V_ = newton_projection(solver,res_i,eps=1e-8,verbose=false)
# res_ = ConstrainedVectorResults(solver,V_.Z.X,V_.Z.U)
# backwardpass!(res_,solver)
# rollout!(res_,solver,0.0)
# max_violation(res_)
# t_newton = float(time_ns()-t_start)/1e9
#
# t_i = 0.469548  # from running @btime
# t_p = 0.659831  # from running @btime
# t_2 = 11.325    # from running @btime
# time_i = collect(range(0,stop=t_i,length=stats_i["iterations"]))
# time_p = range(0,stop=t_p,length=stats_p["iterations"])
# p = plot(time_p,stats_p["c_max"][2:end],yscale=:log10,label="Ipopt",color=:blue,width=2,
#     markershape=:circle,markerstrokecolor=:blue)
# c_max = [stats_i["c_max"]; max_violation(res_)]
# push!(time_i,t_i+t_newton)
# plot!(time_i,c_max,label="ALTRO*",color=:green,width=2,markershape=:circle,markerstrokecolor=:green,
#     xlabel="runtime (s)",ylabel="max constraint violation")
# plot_vertical_lines!(p,[t_i])
# time_2 = range(0,step=time_i[2]-time_i[1],length=stats_2["iterations"])
# plot!(time_2,stats_2["c_max"],label="ALTRO",width=2,color=:darkorange2,
#     markershape=:circle,markerstrokecolor=:darkorange2,xlim=[0,1.5],size=(500,250))
# savefig(p,joinpath(IMAGE_DIR,"ppark_newton.eps"))

# Constraint vs Iteration (Newton tail)
using TrajectoryOptimization: gen_usrfun_newton, NewtonVars, gen_newton_functions, newton_projection
t_start = time_ns()
solver.opts.iterations = 25
solver.opts.constraint_tolerance = 1e-2
res_i, stats_i = solve(solver,U0)
max_violation(res_i)

plot(stats_i["c_max"],yscale=:log10)

V_ = newton_projection(solver,res_i,eps=1e-8,verbose=false)
res_ = ConstrainedVectorResults(solver,V_.Z.X,V_.Z.U)
backwardpass!(res_,solver)
rollout!(res_,solver,0.0)
max_violation(res_)
t_newton = float(time_ns()-t_start)/1e9

t_i = 0.469548  # from running @btime
t_p = 0.659831  # from running @btime
t_2 = 11.325    # from running @btime
time_i = collect(range(0,stop=t_i,length=stats_i["iterations"]))
time_p = range(0,stop=t_p,length=stats_p["iterations"])
p = plot(time_p,stats_p["c_max"][2:end],yscale=:log10,label="Ipopt",color=:blue,width=2,
    markershape=:circle,markerstrokecolor=:blue)
c_max = [stats_i["c_max"]; max_violation(res_)]
push!(time_i,t_i+t_newton)
plot!(time_i,c_max,label="ALTRO*",color=:green,width=2,markershape=:circle,markerstrokecolor=:green,
    xlabel="runtime (s)",ylabel="max constraint violation")
plot_vertical_lines!(p,[t_i])
time_2 = range(0,step=time_i[2]-time_i[1],length=stats_2["iterations"])
plot!(time_2,stats_2["c_max"],label="ALTRO",width=2,color=:darkorange2,
    markershape=:circle,markerstrokecolor=:darkorange2,xlim=[0,1.5],size=(500,250))
savefig(p,joinpath(IMAGE_DIR,"ppark_newton.eps"))

# #### INFEASIBLE ####
# opts = SolverOptions()
# opts.verbose = false
# opts.cost_tolerance = 1e-6
# opts.cost_tolerance_intermediate = 1e-3
# opts.cost_tolerance_infeasible = 1e-4
# opts.constraint_tolerance = 1e-4
# opts.resolve_feasible = false
# opts.outer_loop_update_type = :default
# opts.use_nesterov = true
# opts.penalty_scaling = 200
# opts.penalty_initial = 10
# opts.R_infeasible = 10
# opts.square_root = true
# opts.constraint_decrease_ratio = 0.25
# opts.penalty_update_frequency = 2
#
# solver = Solver(model, obj, N=101, opts=opts)
# n,m,N = get_sizes(solver)
# U0 = ones(m,N)
# X0 = line_trajectory(solver)
# X0_rollout = rollout(solver,U0)
#
# solver.opts.verbose = false
# solver.opts.resolve_feasible = true
# solver.opts.cost_tolerance_infeasible = 1e-5
# @time res_i, stats_i = solve(solver,X0,U0)
# stats_i["iterations"]
# res_i.U[1]
# stats_i["runtime"]
# stats_i["iterations (infeasible)"]
# @time res_p, stats_p = solve_dircol(solver,X0_rollout,U0)
# # res_s, stats_s = solve_dircol(solver,X0_rollout,U0,nlp=:snopt)
#
# # import TrajectoryOptimization: _solve, get_feasible_trajectory
# # results = copy(res_i)
# # solver.state.infeasible = true
# # results_feasible = get_feasible_trajectory(results, solver)
# # results_feasible.λ_prev[1]
#
# # solver.state.infeasible
# # res, stats = _solve(solver,to_array(results_feasible.U),prevResults=results_feasible); stats["iterations"]
# # res, stats = _solve(solver,to_array(results_feasible.U)); stats["iterations"]
# # res, stats = _solve(solver,to_array(results_feasible.U),λ=results_feasible.λ); stats["iterations"]
# # res, stats = _solve(solver,to_array(results_feasible.U),λ=results_feasible.λ,μ=results_feasible.μ); stats["iterations"]
# #
#
# constraint_plot(solver,X0,U0)
#
# Ns = [101,201,301,401]
# group = "parallelpark/infeasible"
# run_step_size_comparison(model, obj, U0, group, Ns, opts=opts, integrations=[:rk3,:ipopt],benchmark=true, infeasible=true)
# plot_stat("runtime",group,legend=:bottomright,["rk3","ipopt"],title="Constrained Parallel Park (infeasible)",ylim=[0,1.5])
# plot_stat("iterations",group,legend=:bottom,["rk3","ipopt"],title="Constrained Parallel Park (infeasible)")
# plot_stat("error",group,yscale=:log10,legend=:right,["rk3","ipopt"],title="Constrained Parallel Park (infeasible)")
#
#
#
# Combined Plot
group = "parallelpark/constrained"
Ns, data = load_data("runtime","ipopt","parallelpark/constrained")
Ns, err = load_data("std","ipopt","parallelpark/constrained")
p1 = plot(Ns,data,yerr=err,label="DIRCOL",color=:blue,marker=:circle,markerstrokecolor=:blue,ylabel="runtime",markersize=6)
Ns, data = load_data("runtime","rk3","parallelpark/constrained")
Ns, err = load_data("std","rk3","parallelpark/constrained")
plot!(Ns,data,yerr=err,label="ALTRO",color=:darkorange2,marker=:circle,markerstrokecolor=:darkorange2,markersize=6,ylim=(0,1.9))
Ns, data = load_data("runtime","rk3","parallelpark/infeasible")
Ns, err = load_data("std","rk3","parallelpark/infeasible")
plot!(Ns,data,yerr=err,label="ALTRO (inf)",color=:darkorange2,style=:dash,
    marker=:utriangle,markerstrokecolor=:darkorange2,markersize=8,
    title="Constrained",titlefontsize=10)

Ns, err = load_data("std","rk3","parallelpark/infeasible")

Ns, data = load_data("runtime",["rk3","ipopt"],"parallelpark/unconstrained")
Ns, err = load_data("std",["rk3","ipopt"],"parallelpark/unconstrained")
p2 = plot(Ns,data[2],yerr=err[2],color=:blue,style=:dot,label="DIRCOL", width=1.5,
    marker=:square,markerstrokecolor=:blue,legend=:topleft,markersize=4,ylim=ylims(p1))
plot!(Ns,data[1],yerr=err[1],color=:darkorange2,style=:dot,label="ALTRO",width=1.5,
    marker=:square,markerstrokecolor=:darkorange2,markersize=4,
    title="Unconstrained",titlefontsize=10)
plot(p1,p2,layout=(1,2),size=(500,300),dpi=400,xlabel="Knot points",ylabel="Runtime")

savefig(joinpath(IMAGE_DIR,"ppark_runtime.png"))


text:
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(path::String)::Int
    base_exp = map(x -> parse.(Int, split(x, ",")), split(read(path, String), "\n"))
    max_exponent = 0
    result = 0
    for i ∈ 1:length(base_exp)
        base, exponent = base_exp[i]
        exponent = log10(base) * exponent
        if exponent > max_exponent
            max_exponent = exponent
            result = i
        end
    end
    return result
end

compute("problems/0099/p099_base_exp.txt")

@benchmark compute("problems/0099/p099_base_exp.txt")

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
function safe_rm(file::String)
  try
    rm(file)
  catch err
    @warn(string("Can't remove ", file, ": throws error ", err))
  end
  return nothing
end


text:
const CITE_STR_JL =
"This work was done in \\cite[Julia v$VERSION]{Julia-2017} "*
"and made use of the following packages: "*
"ArrayInterface.jl\\cite[ArrayInterface]{DifferentialEquations.jl-2017}, "*
"LabelledArrays.jl\\cite[LabelledArrays]{DifferentialEquations.jl-2017}, "*
"QuadGK.jl\\cite[QuadGK]{quadgk}, "*
"RecursiveArrayTools.jl\\cite[RecursiveArrayTools]{DifferentialEquations.jl-2017} and "*
"Symbolics.jl\\cite[Symbolics]{gowda2021high}."

const CITE_STR_JL_AUTO =
"This work was done in \\cite[Julia v$VERSION]{Julia-2017} "*
"and made use of the following packages: "*
"ArrayInterface.jl\\cite[ArrayInterface]{DifferentialEquations.jl-2017}, "*
"LabelledArrays.jl\\cite[LabelledArrays]{DifferentialEquations.jl-2017}, "*
"QuadGK.jl\\cite[QuadGK]{quadgk}, "*
"RecursiveArrayTools.jl\\cite[RecursiveArrayTools]{DifferentialEquations.jl-2017} and "*
"Symbolics.jl\\autocite[Symbolics]{gowda2021high}."

const CITE_STR_JL_SINGLE =
"This work was done in \\cite[Julia v$VERSION]{Julia-2017} "*
"and made use of the Symbolics.jl\\cite[Symbolics]{gowda2021high} package."

const CITE_STR_JL_SINGLE_TT =
"This work was done in \\cite[Julia v$VERSION]{Julia-2017} "*
"and made use of the \\texttt{Symbolics.jl}\\cite[Symbolics]{gowda2021high} package."

const CITE_STR =
"This work was done in \\cite[Julia v$VERSION]{Julia-2017} "*
"and made use of the following packages: "*
"ArrayInterface\\cite[ArrayInterface]{DifferentialEquations.jl-2017}, "*
"LabelledArrays\\cite[LabelledArrays]{DifferentialEquations.jl-2017}, "*
"QuadGK\\cite[QuadGK]{quadgk}, "*
"RecursiveArrayTools\\cite[RecursiveArrayTools]{DifferentialEquations.jl-2017} and "*
"Symbolics\\cite[Symbolics]{gowda2021high}."

const CITE_BADGE =
"This work was done in \\cite[Julia v$VERSION]{Julia-2017} "*
"and made use of the Plots.jl\\cite[Plots]{https://doi.org/10.5281/zenodo.4725317} "*
"and WriteVTK.jl\\cite[WriteVTK]{https://doi.org/10.5281/zenodo.5068937} packages."


text:
using TuringCLIExample
TuringCLIExample.main()

max_stars_repo_name: essiequoi/TuringCLIExample

max_issues_repo_name: essiequoi/TuringCLIExample

max_forks_repo_name: essiequoi/TuringCLIExample


text:
using BinaryBuilder

name = "blaspp"
version = v"2021.04.02" # +1 patch version for switching back to openblas

# Collection of sources required to build PETSc. Avoid using the git repository, it will
# require building SOWING which fails in all non-linux platforms.
sources = [
    GitSource("https://bitbucket.org/icl/blaspp.git", "314bafceead689a35aab826e03aa76bf329cfb0e")
]

# Bash recipe for building across all platforms

# Needs to add -Dcapi eventually once it's added to the cmake build system. Note yet available under CMAKAE toolchain.
script = raw"""
cd blaspp
mkdir build && cd build
cmake \
  -DCMAKE_INSTALL_PREFIX=${prefix} \
  -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
  -DCMAKE_BUILD_TYPE="Release" \
  -Drun_result="0" \
  -Drun_result__TRYRUN_OUTPUT="ok" \
  -Dgpu_backend=none \
  -Dbuild_tests=no \
  ..
make -j${nproc}
make install
"""

# We attempt to build for all defined platforms
platforms = expand_cxxstring_abis(supported_platforms(;experimental=true))
products = [
    LibraryProduct("libblaspp", :libblaspp)
]

dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("OpenBLAS32_jll"),
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"7")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
function find_in_path(name::AbstractString; constants::Bool=true)
    isabspath(name) && return name
    isfile(name) && return abspath(name)
    base = name
    if endswith(name, ".jl")
        base = name[1:end-3]
    else
        name = string(base,".jl")
        isfile(name) && return abspath(name)
    end
    for prefix in [Pkg.dir(); LOAD_PATH]
        path = joinpath(prefix, name)
        if constants
            file_name = name
        else
            file_name = string(name[1:end-3], "_reload", ".jl")
            if !isfile(prefix, base, "src", file_name)
                #warn("Reload version of $name not found")
                file_name = name
            end
        end
        isfile(path) && return abspath(path)
        path = joinpath(prefix, base, "src", file_name)
        isfile(path) && return abspath(path)
        path = joinpath(prefix, name, "src", file_name)
        isfile(path) && return abspath(path)
    end
    return nothing
end

function standarize(filename)
    if endswith(filename, ".jl")
        filename
    else
        string(filename, ".jl")
    end
end

function find_file(filename; base_file=nothing, constants=true)
    path = find_in_path(filename, constants=constants)
    if path == nothing && base_file!=nothing
        base_file = Base.find_in_path(base_file)
        path = joinpath(dirname(base_file), filename)
        if !isfile(path)
            path = nothing
        end
    end
    return path
end

function reload_mtime(filename)
    path = find_file(filename)
    m = mtime(path)
    if m == 0.0
        warn("Could not find edit time of $filename")
    end
    m
end

function parse_file(filename; kwargs...)
    path = find_file(filename; kwargs...)
    handle = open(path)
    source = string("begin\n", readall(handle), "\n end")
    close(handle)
    parsed = parse(source)
    return parsed
end


text:
## initialize module
using Pkg;Pkg.activate(".")
using VortexDistributions

## other stuff and a change
using Test, Plots

## set up mock data
N = 10
xi = rand(N)
yi = rand(N)

## make cluster using vortices and tree
vort = PointVortex.(xi,yi,one.(xi))

# default constructor
c0 = Cluster(vort,spanning_tree(xi,yi))

## callable type
c1 = Cluster(vort)
@test c1.vortices |> length == N
@test c0.vortices[1] == c1.vortices[1]
@test c0.tree[1] == c1.tree[1]

## plot cluster
p1 = plot()
plot_cluster!(p1,c1) # check that the plot method still works

## do some random tests of robustness
N = 50
xi,yi = randn(N),randn(N)
tree = spanning_tree(xi,yi)
c1 = Cluster(PointVortex.(xi,yi,one.(xi)),tree)
p1 = plot(grid=false,axis=false)
plot_cluster!(p1,c1)

xi,yi = randn(N),randn(N)
xi .+= 5
tree = spanning_tree(xi,yi)
c1 = Cluster(PointVortex.(xi,yi,-one.(xi)),tree)
plot_cluster!(p1,c1)

## vector of clusters of different sizes
xi,yi = randn(2N),randn(2N)
xi .+= 10
tree = spanning_tree(xi,yi)
c2 = Cluster(PointVortex.(xi,yi,one.(xi)),tree)
cvec = [c0;c1;c2]

xi,yi = randn(3N),randn(3N)
xi .+= 20
tree = spanning_tree(xi,yi)
c3 = Cluster(PointVortex.(xi,yi,-one.(xi)),tree)
push!(cvec,c3)

## cvec
Base.size(c::Cluster) = length(c.vortices)
cvec[3] |> size
size(cvec)

## save
#savefig(p1,"rca.pdf")

## grow_plus_clusters


text:
#=
Load the PF00011 family from PFAM.
=#

import Pfam
msa = Pfam.MSA("PF00011");  ## describe the requested MSA
df = Pfam.load(msa); ## load as a DataFrame
@show df[1:50,:]; ## show some rows


text:
@testset "service" begin
    @testset "detach_indexer" begin
        repo = Repository(
            [
                Indexer("0x00", 10.0, [Allocation("Qmaaa", 10.0, 0)]),
                Indexer("0x01", 20.0, [Allocation("Qmaaa", 10.0, 0)]),
            ],
            [SubgraphDeployment("1x00", "Qmaaa", 10.0)],
        )

        # Should detach the specified indexer
        id = "0x00"
        indexer, frepo = detach_indexer(repo, id)
        @test frepo.indexers[1].id == "0x01"
        @test length(frepo.indexers) == 1
        @test indexer.id == "0x00"

        # Should throw an exception when `id` not in repo
        id = "0x10"
        @test_throws UnknownIndexerError detach_indexer(repo, id)
    end

    @testset "signal" begin
        subgraphs = [SubgraphDeployment("1x00", "Qmaaa", 10.0)]

        # Should get the signal of the one defined subgraph
        @test signal.(subgraphs) == [10.0]
    end

    @testset "ipfshash" begin
        subgraphs = [SubgraphDeployment("1x00", "Qmaaa", 10.0)]

        # Should get the ipfshash of the one defined subgraph
        @test ipfshash.(subgraphs) == ["Qmaaa"]
    end

    @testset "allocation" begin
        indexers = [
            Indexer("0x00", 10.0, [Allocation("Qmaaa", 10.0, 0)]),
            Indexer("0x01", 20.0, [Allocation("Qmaaa", 20.0, 0)]),
        ]

        # Should get the allocations of indexers
        @test allocation.(indexers) ==
            [[Allocation("Qmaaa", 10.0, 0)], [Allocation("Qmaaa", 20.0, 0)]]
    end

    @testset "allocated_stake" begin
        allocs = [Allocation("Qmaaa", 10.0, 0), Allocation("Qmaaa", 20.0, 0)]

        # Should get the allocations of indexers
        @test allocated_stake.(allocs) == [10.0, 20.0]
    end

    @testset "stakes" begin
        repo = Repository(
            [
                Indexer(
                    "0x00",
                    10.0,
                    [Allocation("Qmaaa", 10.0, 0), Allocation("Qmbbb", 10.0, 0)],
                ),
                Indexer("0x01", 20.0, [Allocation("Qmaaa", 10.0, 0)]),
            ],
            [SubgraphDeployment("1x00", "Qmaaa", 10.0)],
        )
        # Should get the sum of stake of the one defined subgraph
        @test stakes(repo) == [20.0]

        repo = Repository(
            [
                Indexer("0x00", 10.0, [Allocation("Qmaaa", 10.0, 0)]),
                Indexer("0x01", 20.0, [Allocation("Qmccc", 8.0, 0)]),
            ],
            [
                SubgraphDeployment("1x00", "Qmaaa", 10.0),
                SubgraphDeployment("1x01", "Qmbbb", 5.0),
                SubgraphDeployment("1x01", "Qmccc", 5.0),
            ],
        )

        # Should get proper sum of stake according to the subgraph ids
        @test stakes(repo) == [10.0, 0.0, 8.0]
    end

    @testset "optimize" begin
        repo = Repository(
            [
                Indexer(
                    "0x01", 10.0, [Allocation("Qmaaa", 2.0, 0), Allocation("Qmbbb", 8.0, 0)]
                ),
            ],
            [
                SubgraphDeployment("1x00", "Qmaaa", 10.0),
                SubgraphDeployment("1x01", "Qmbbb", 5.0),
            ],
        )
        indexer = Indexer("0x00", 5.0, Allocation[])
        ω = optimize(indexer, repo)
        @test isapprox(ω, [4.2, 0.8], atol=0.1)
    end
end


text:
# Anonymous functions have no `new` constructor. Here we generate
# one for them based on the types of args passed to FunctionConstructor

struct FunctionConstructor{F} end

@generated function (fc::FunctionConstructor{F})(args...) where F
    T = getfield(parentmodule(F), nameof(F))
    # We assume all gensym names are anonymous functions
    _isgensym(nameof(F)) || return :($T(args...))
    # Define `new` for rebuilt function type that matches args
    exp = Expr(:new, Expr(:curly, T, args...))
    for i in 1:length(args)
        push!(exp.args, :(args[$i]))
    end
    return exp
end

function ConstructionBase.constructorof(f::Type{F}) where F <: Function
    FunctionConstructor{F}()
end

_isgensym(s::Symbol) = occursin("#", string(s))