text:
using DiffEqFlux, OrdinaryDiffEq, Flux, Plots
using Zygote
using DifferentialEquations
using ForwardDiff
using LinearAlgebra
using BSON: @save
using BSON: @load
LinearAlgebra.BLAS.set_num_threads(16)

function rober(du,u,p,t)
    y₁,y₂,y₃ = u
    k₁,k₂,k₃ = p
    du[1] = -k₁*y₁+k₃*y₂*y₃
    du[2] =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
    du[3] =  k₂*y₂^2
    nothing
end

function generate_data(u0, pz, tspan, solver)
    ode_prob = ODEProblem(rober,u0,tspan,pz)
    ode_sol = solve(ode_prob, solver)
    data = Array(ode_sol)
    return data, ode_sol.t
end

function predict_n_ode(k)
    Array(solve(prob, Rosenbrock23(), u0=u0,p=k, saveat=tsteps))
end

function loss_n_ode(k)
    #loss for stiff-net
      pred = predict_n_ode(k)
      loss = sum(abs2,ode_data .- pred)
      return loss
end

function grab_∇(par)
    #Grabs gradients of NN loss w.r.t params
    ForwardDiff.gradient(loss_n_ode, par)
end

function dθ_dt!(du,u,k,t)
    #ODE used to evolve parmaeters based on the negative gradient definition (evolving in the direction of steepest descent)
    du .= -1 .* grab_∇(u)
end

#################################################################################################################################################################################
#Generate the data
pz = [0.04,3e7,1e4]
u0 = [1.0,0.0,0.0]
tspan = (1.0e-5,1.0e5)
solver = Rosenbrock23(autodiff=false,diff_type=Val{:forward})
ode_data, tsteps = generate_data(u0, pz, tspan, solver)
plot(tsteps, ode_data')

#Create the NN
dudt2 = Chain(Dense(3,3, tanh),
             Dense(3, 3, tanh),
             Dense(3, 3, tanh),
             Dense(3,3))

@load "stiff-net.bson" dudt2

p,re = Flux.destructure(dudt2) # use this p as the initial condition!
dudt(u,p,t) = re(p)(u) # need to restrcture for backprop!
prob = ODEProblem(dudt,u0,tspan)
#Test predicion functions
predict_n_ode(p)
#Test loss functions, check for good initialisation (loss must be <100)
loss_n_ode(p)
#Save well behaved initialisation
#@save "stiff-net.bson" dudt2

#Train with solver
tspan_θ = (0.0,0.125)
prob_θ = ODEProblem(dθ_dt!,p,tspan_θ)
integrator = init(prob_θ,Rosenbrock23(autodiff=false))
loss_arr = []
params_arr = []
for (u,t) in tuples(integrator)
    loss = loss_n_ode(u)
    push!(loss_arr, loss)
    push!(params_arr, u)
    println(loss)
end

plot(loss_arr)

#Plotting first concentration
plotly()
plot(tsteps, predict_n_ode(param_his[:, end])[1, :])
scatter!(tsteps, ode_data[1, :], title = "First concentration", label = ["NN-pred" "Ground truth"])

#Plotting second concentration
plot(tsteps, predict_n_ode(param_his[:, end])[2, :])
scatter!(tsteps, ode_data[2, :], title = "Second concentration", label = ["NN-pred" "Ground truth"])

#Plotting third concentration
plot(tsteps, predict_n_ode(param_his[:, end])[3, :])
scatter!(tsteps, ode_data[3, :], title = "Third concentration", label = ["NN-pred" "Ground truth"])

#Plotting losss
loss = [loss_n_ode(param_his[:,i]) for i in 1:size(param_his)[2]]
plot(loss, title = "Loss")


text:
function spg_get_ir_reciprocal_mesh(mesh::Array{Int64, 1},
                            is_shift::Array{Int64, 1},
                            is_time_reversal::Int64,
                            lattice::Array{Float64, 2},
                            positions::Array{Float64, 2},
                            types::Array{Int64, 1},
                            num_atom::Int64,
                            symprec::Float64=1e-5)
    #
    mesh = Base.cconvert( Array{Int32,1}, mesh)
    is_shift = Base.cconvert( Array{Int32,1}, is_shift )
    types = Base.cconvert( Array{Int32,1}, types)
    num_atom = Base.cconvert( Int32, num_atom )
    is_time_reversal = Base.cconvert( Int32, is_time_reversal )

    # Prepare for output
    nkpts = prod(mesh)
    grid_address = zeros(Int32,3,nkpts)
    ir_mapping_table = zeros(Int32,nkpts)

    num_ir =
    ccall((:spg_get_ir_reciprocal_mesh, libsymspg), Int32,
          (Ptr{Int32}, Ptr{Int32}, Ptr{Int32},
           Ptr{Int32}, Int32, Ptr{Float64}, Ptr{Float64},
           Ptr{Int32}, Int32, Float64),
           grid_address, ir_mapping_table, mesh, is_shift, is_time_reversal,
           lattice, positions, types, num_atom, symprec)

    return Base.cconvert(Int64, num_ir),
           Base.cconvert(Array{Int64,2}, grid_address),
           Base.cconvert(Array{Int64,1}, ir_mapping_table)
end

ext: jl


text:
"""
    parameters(t, A, n0, m_max[, p_max[, tol]]) -> NamedTuple{(:m, :s)}

Compute Taylor series parameters needed for `exp(t*A) * B`.

This is Code Fragment 3.1 from [^AlMohyHigham2011].
"""
function parameters(
    t, A, n0, m_max, p_max=p_from_m(m_max), tol=eps(float(real(Base.promote_eltype(t, A))))
)
    tnorm = abs(t)
    iszero(tnorm) && return (m=0, s=1)
    Anorm = opnormest1(A)
    iszero(Anorm) && return (m=0, s=1)
    tAnorm = tnorm * Anorm
    T = float(real(Base.promote_eltype(t, A)))
    θ = coefficients(T(tol))
    ℓ = 2 # §3: “where the positive integer ℓ is a parameter (typically set to 1 or 2)”
    if tAnorm * (n0 * m_max) ≤ θ[m_max] * (2 * ℓ * p_max * (p_max + 3)) # (3.13) is satisfied
        m_opt = 0
        Cm_opt = T(Inf)
        # work around argmin not taking a function
        for m in 1:m_max
            Cm = m * ceil(Int, tAnorm / θ[m])
            if Cm < Cm_opt
                m_opt = m
                Cm_opt = Cm
            end
        end
        s = ceil(Int, Cm_opt / m_opt)
    else
        # TODO: replace powers of A here and below with opnormest(pow, A, 1)
        # see https://github.com/JuliaLang/julia/pull/39058
        Apow = A * A
        d = tnorm * sqrt(opnormest1(A))
        m_opt = m_max + 1
        Cm_opt = T(Inf)
        for p in 2:p_max # Compute minimum in (3.11)
            Apow *= A
            # (3.7)
            d, d_old = tnorm * opnormest1(Apow)^(1//(p + 1)), d
            α = max(d, d_old)
            m_min = p * (p - 1) - 1
            # work around argmin not taking a function
            for m in m_min:m_max
                Cm = m * ceil(Int, α / θ[m])
                if Cm < Cm_opt || (Cm == Cm_opt && m < m_opt)
                    m_opt = m
                    Cm_opt = Cm
                end
            end
        end
        s = max(ceil(Int, Cm_opt / m_opt), 1)
    end
    return (m=m_opt, s=s)
end
# work around opnorm(A, 1) and (A^2)*A having very slow defaults for these arrays
# https://github.com/sethaxen/ExponentialAction.jl/issues/3
function parameters(
    t,
    A::Union{Bidiagonal,Tridiagonal},
    n0,
    m_max,
    p_max=p_from_m(m_max),
    tol=eps(float(real(Base.promote_eltype(t, A)))),
)
    return parameters(t, sparse(A), n0, m_max, p_max, tol)
end

# avoid differentiating through parameters with ChainRules-compatible ADs
ChainRulesCore.@non_differentiable parameters(t, A, n0, m_max, p_max, tol)

# solution to p(p-1) ≤ m + 1
p_from_m(m) = oftype(m, fld(1 + sqrt(5 + 4 * m), 2))


text:
## Coordinates
#==========================================================================================#

wing_bounds(lead, trail) = permutedims([ lead trail ])

chop_leading_edge(obj  :: HalfWing, span_num; y_flip = false) = chop_coordinates(leading_edge(obj, y_flip),  span_num)
chop_trailing_edge(obj :: HalfWing, span_num; y_flip = false) = chop_coordinates(trailing_edge(obj, y_flip), span_num)

chop_leading_edge(obj :: Wing, span_num :: Integer)  = chop_coordinates([ leading_edge(left(obj), true)[1:end-1]; leading_edge(right(obj)) ], span_num)
chop_trailing_edge(obj :: Wing, span_num :: Integer) = chop_coordinates([ trailing_edge(left(obj), true)[1:end-1]; trailing_edge(right(obj)) ], span_num)

coordinates(wing :: HalfWing, y_flip = false) = let (lead, trail) = wing_bounds(wing, y_flip); affine_transformation(wing).(wing_bounds(lead, trail)) end
coordinates(wing :: Wing) = let (lead, trail) = wing_bounds(wing); affine_transformation(wing).(wing_bounds(lead, trail)) end

chord_coordinates(wing :: HalfWing, span_num :: Vector{<: Integer}, chord_num :: Integer; spacings = symmetric_spacing(wing), flip = false) = chop_wing(coordinates(wing, flip), span_num, chord_num; span_spacing = spacings, flip = flip)

function camber_coordinates(wing :: HalfWing, span_num :: Vector{<: Integer}, chord_num :: Integer; spacings = symmetric_spacing(wing), flip = false)
    leading_xyz  = leading_edge(wing, flip)
    scaled_foils = @. wing.chords * (camber_coordinates ∘ camber_thickness)(wing.foils, chord_num)
    affine_transformation(wing).(chop_spanwise_sections(scaled_foils, twists(wing), leading_xyz, span_num, spacings, flip))
end

function surface_coordinates(wing :: HalfWing, span_num :: Vector{<: Integer}, chord_num :: Integer; spacings = symmetric_spacing(wing), flip = false)
    leading_xyz  = leading_edge(wing, flip)
    scaled_foils = @. wing.chords * (extend_yz ∘ coordinates ∘ cosine_interpolation)(wing.foils, chord_num)
    affine_transformation(wing).(chop_spanwise_sections(scaled_foils, twists(wing), leading_xyz, span_num, spacings, flip))
end

function number_of_spanwise_panels(wing :: HalfWing, span_num :: Integer)
    # Compute contribution of each section to total span length
    weights = spans(wing) / span(wing)

    weights[findall(<(0.2), weights)] .*= 3

    # Heuristic (aka hAx0rZ) check to ensure small sections also get some panel love
    # weights = ifelse(any(<(0.2), weights), fill(1. / length(spans(wing)), length(spans(wing))), weights)

    # Generate spanwise panel distribution
    ceil.(Int, span_num .* weights)
end

function number_of_spanwise_panels(wing :: HalfWing, span_num :: Vector{<: Integer})
    @assert (length ∘ spans)(wing) > 1 "Provide a positive integer of spanwise panels for 1 wing section."
    span_num
end

# Spacing
symmetric_spacing(wing :: HalfWing) = [ Sine(); fill(Cosine(), (length ∘ spans)(wing) - 1) ]

## Wing variants
#==========================================================================================#

function chord_coordinates(wing :: Wing, span_num :: Vector{<: Integer}, chord_num :: Integer; spacings = symmetric_spacing(wing))
    left_coord  = chord_coordinates(left(wing), reverse(span_num), chord_num; spacings = reverse(spacings), flip = true)
    right_coord = chord_coordinates(right(wing), span_num, chord_num; spacings = spacings)

    [ left_coord[:,1:end-1] right_coord ]
end

function camber_coordinates(wing :: Wing, span_num :: Vector{<: Integer}, chord_num :: Integer; spacings = symmetric_spacing(wing))
    left_coord  = camber_coordinates(left(wing), reverse(span_num), chord_num; spacings = reverse(spacings), flip = true)
    right_coord = camber_coordinates(right(wing), span_num, chord_num; spacings = spacings)

    [ left_coord[:,1:end-1] right_coord ]
end

function surface_coordinates(wing :: Wing, span_num :: Vector{<: Integer}, chord_num :: Integer; spacings = symmetric_spacing(wing))
    left_coord  = surface_coordinates(left(wing), reverse(span_num), chord_num; spacings = reverse(spacings), flip = true)
    right_coord = surface_coordinates(right(wing), span_num, chord_num; spacings = spacings)

    [ left_coord[:,1:end-1] right_coord ]
end

number_of_spanwise_panels(wing :: Wing, span_num :: Integer) = number_of_spanwise_panels(right(wing), span_num ÷ 2)

symmetric_spacing(wing :: Wing)     = [ Sine(); fill(Cosine(), (length ∘ spans ∘ right)(wing) - 1) ]

function number_of_spanwise_panels(wing :: Wing, span_num :: Vector{<: Integer})
    @assert (length ∘ spans ∘ right)(wing) > 1 "Provide a positive integer of spanwise panels for 1 wing section."
    span_num .÷ 2
end

# Coordinates
chord_coordinates(wing :: Wing, span_num :: Integer, chord_num :: Integer; spacings = symmetric_spacing(wing)) = chord_coordinates(wing, number_of_spanwise_panels(wing, span_num), chord_num; spacings = spacings)

camber_coordinates(wing :: Wing, span_num :: Integer, chord_num :: Integer; spacings = symmetric_spacing(wing)) = camber_coordinates(wing, number_of_spanwise_panels(wing, span_num), chord_num; spacings = spacings)

surface_coordinates(wing :: Wing, span_num :: Integer, chord_num :: Integer; spacings = symmetric_spacing(wing)) = surface_coordinates(wing, number_of_spanwise_panels(wing, span_num), chord_num; spacings = spacings)

## Panelling
#==========================================================================================#

mesh_chords(wing :: HalfWing, span_num :: Vector{<: Integer}, chord_num :: Integer; spacings = symmetric_spacing(wing), flip = false) = make_panels(chord_coordinates(wing, span_num, chord_num; spacings = spacings, flip = flip))

function mesh_chords(wing :: Wing, span_num, chord_num; spacings = symmetric_spacing(wing))
    left_panels  = mesh_chords(left(wing), reverse(span_num), chord_num; spacings = reverse(spacings), flip = true)
    right_panels = mesh_chords(right(wing), span_num, chord_num; spacings = spacings)

    [ left_panels right_panels ]
end

mesh_wing(wing :: HalfWing, span_num :: Vector{<: Integer}, chord_num :: Integer; spacings = symmetric_spacing(wing), flip = false) = make_panels(surface_coordinates(wing, span_num, chord_num, spacings = spacings, flip = flip))

function mesh_wing(wing :: Wing, span_num, chord_num; spacings = symmetric_spacing(wing))
    left_panels  = mesh_wing(left(wing), reverse(span_num), chord_num; spacings = reverse(spacings), flip = true)
    right_panels = mesh_wing(right(wing), span_num, chord_num; spacings = spacings)

    [ left_panels right_panels ]
end

mesh_cambers(wing :: HalfWing, span_num :: Vector{<: Integer}, chord_num :: Integer; spacings = symmetric_spacing(wing), flip = false) = make_panels(camber_coordinates(wing, span_num, chord_num; spacings = spacings, flip = flip))

function mesh_cambers(wing :: Wing, span_num, chord_num; spacings = symmetric_spacing(wing))
    left_panels  = mesh_cambers(left(wing), reverse(span_num), chord_num; spacings = reverse(spacings), flip = true)
    right_panels = mesh_cambers(right(wing), span_num, chord_num; spacings = spacings)

    [ left_panels right_panels ]
end

function make_panels(wing :: AbstractWing, span_num :: Vector{<: Integer}, chord_num :: Integer; spacings = symmetric_spacing(wing))
    horseshoe_panels = mesh_chords(wing, span_num, chord_num; spacings = spacings)
    camber_panels    = mesh_cambers(wing, span_num, chord_num; spacings = spacings)
    horseshoe_panels, normal_vector.(camber_panels)
end

make_panels(wing :: AbstractWing, span_num :: Integer, chord_num :: Integer; spacings = symmetric_spacing(wing)) = make_panels(wing, [span_num], chord_num; spacings = spacings)

panel_wing(comp :: AbstractWing, span_panels :: Union{Integer, Vector{<: Integer}}, chord_panels :: Integer; spacing = symmetric_spacing(comp)) = make_panels(comp, span_panels, chord_panels, spacings = spacing)

## Meshing type for convenience
#==========================================================================================#

struct WingMesh{M <: AbstractWing, N <: Integer, P, Q, T} <: AbstractWing
    surface       :: M
    num_span      :: Vector{N}
    num_chord     :: N
    chord_spacing :: P
    span_spacing  :: Q
    chord_mesh    :: Matrix{T}
    camber_mesh   :: Matrix{T}
end

"""
    WingMesh(
             surface :: AbstractWing,
             n_span :: Vector{Integer}, n_chord :: Integer;
             span_spacing :: AbstractSpacing = symmetric_spacing(surface)
            )

Define a container to generate meshes and panels for a given `AbstractWing` with a specified distribution of number of spanwise panels, and a number of chordwise panels.

Optionally a combination of `AbstractSpacing` types (`Sine(), Cosine(), Uniform()`) can be provided to the **named argument** `span_spacing`, either as a singleton or as a vector with length equal to the number of spanwise sections. By default, the combination is `[Sine(), Cosine(), ..., Cosine()]`.
"""
function WingMesh(surface :: M, n_span :: AbstractVector{N}, n_chord :: N; chord_spacing :: P = Cosine(), span_spacing :: Q = symmetric_spacing(surface)) where {M <: AbstractWing, N <: Integer, P <: AbstractSpacing, Q <: Union{AbstractSpacing, Vector{<:AbstractSpacing}}}
    check_definition(surface, n_span)
    chord_mesh  = chord_coordinates(surface, n_span, n_chord; spacings = span_spacing)
    camber_mesh = camber_coordinates(surface, n_span, n_chord; spacings = span_spacing)
    T = promote_type(eltype(chord_mesh), eltype(camber_mesh))
    WingMesh{M,N,P,Q,T}(surface, n_span, n_chord, chord_spacing, span_spacing, chord_mesh, camber_mesh)
end

WingMesh(surface, n_span :: Integer, n_chord :: Integer; chord_spacing = Cosine(), span_spacing = symmetric_spacing(surface)) = WingMesh(surface, number_of_spanwise_panels(surface, n_span), n_chord; chord_spacing = chord_spacing, span_spacing = span_spacing)

check_definition(surf :: HalfWing, n_span) = @assert length(n_span) == length(surf.spans) "The spanwise number vector's length must be the same as the number of sections of the surface."
check_definition(surf :: Wing, n_span) = @assert length(n_span) == length(surf.right.spans) == length(surf.left.spans) "The spanwise number vector's length must be the same as the number of sections of the surface."

##
"""
    chord_coordinates(wing :: WingMesh, n_span = wing.num_span, n_chord = wing.num_chord)

Generate the chord coordinates of a `WingMesh` with default spanwise ``n_s`` and chordwise ``n_c`` panel distributions from the mesh.
"""
chord_coordinates(wing :: WingMesh, n_span = wing.num_span, n_chord = wing.num_chord) = chord_coordinates(wing.surface, n_span, n_chord)

"""
    camber_coordinates(wing :: WingMesh, n_span = wing.num_span, n_chord = wing.num_chord)

Generate the camber coordinates of a `WingMesh` with default spanwise ``n_s`` and chordwise ``n_c`` panel distributions from the mesh.
"""
camber_coordinates(wing :: WingMesh, n_span = wing.num_span, n_chord = wing.num_chord) = camber_coordinates(wing.surface, n_span, n_chord)

"""
    surface_coordinates(wing :: WingMesh, n_span = wing.num_span, n_chord = wing.num_chord)

Generate the surface coordinates of a `WingMesh` with default spanwise ``n_s`` and chordwise ``n_c`` panel distributions from the mesh.
"""
surface_coordinates(wing :: WingMesh, n_span = wing.num_span, n_chord = wing.num_chord) = surface_coordinates(wing.surface, n_span, n_chord)

"""
    surface_panels(wing_mesh :: WingMesh,
                   n_s = wing_mesh.num_span,
                   n_c = length(first(foils(wing_mesh.surface))).x)

Generate the surface panel distribution from a `WingMesh` with the default spanwise ``n_s`` panel distribution from the mesh and the chordwise panel ``n_c`` distribution from the airfoil.

In case of strange results, provide a higher number of chordwise panels to represent the airfoils more accurately
"""
surface_panels(wing :: WingMesh, n_span = wing.num_span, n_chord = length(first(foils(wing.surface)).x))  = (make_panels ∘ surface_coordinates)(wing, n_span, n_chord)

"""
    chord_panels(wing_mesh :: WingMesh)

Generate the chord panel distribution from a `WingMesh`.
"""
chord_panels(wing :: WingMesh) = make_panels(wing.chord_mesh)

"""
    camber_panels(wing_mesh :: WingMesh)

Generate the camber panel distribution from a `WingMesh`.
"""
camber_panels(wing :: WingMesh) = make_panels(wing.camber_mesh)

"""
    wetted_area(wing_mesh :: WingMesh,
                n_s = wing_mesh.num_span,
                n_c = length(first(foils(wing_mesh.surface))).x)

Determine the wetted area ``S_{wet}`` of a `WingMesh` by calculating the total area of the surface panels.
"""
wetted_area(wing :: WingMesh, n_span = wing.num_span, n_chord = length(first(foils(wing.surface)).x)) = wetted_area(surface_panels(wing, n_span, n_chord))

"""
    wetted_area_ratio(wing_mesh :: WingMesh,
                      n_s = wing_mesh.num_span,
                      n_c = length(first(foils(wing_mesh.surface))).x)

Determine the wetted area ratio ``S_{wet}/S`` of a `WingMesh` by calculating the ratio of the total area of the surface panels to the projected area of the `Wing`.
"""
wetted_area_ratio(wing :: WingMesh, n_span = wing.num_span, n_chord = length(first(foils(wing.surface)).x)) = wetted_area(wing, n_span, n_chord) / projected_area(wing.surface)

function Base.show(io :: IO, mesh :: WingMesh)
    n_c, n_s = size(mesh.chord_mesh) .- 1
    println(io, "WingMesh —")
    println(io, "Spanwise panels: ", n_s)
    println(io, "Chordwise panels: ", n_c)
    println(io, "Spanwise spacing: ", mesh.span_spacing)
    println(io, "Chordwise spacing: ", mesh.chord_spacing)

    nothing
end


text:
module Cbc

if haskey(ENV,"JULIA_CBC_LIBRARY_PATH") || VERSION < v"1.3"
    deps_file = joinpath(dirname(@__DIR__), "deps", "deps.jl")
    if isfile(deps_file)
        using Libdl
        include(deps_file)
    else
        error("Cbc not properly installed. Please run import `Pkg; Pkg.build(\"Cbc\")`.")
    end
else
    import Cbc_jll: libcbcsolver
end

using CEnum

include("gen/ctypes.jl")
include("gen/libcbc_common.jl")
include("gen/libcbc_api.jl")

const _CBC_VERSION = VersionNumber(unsafe_string(Cbc_getVersion()))

if !(v"2.10.0" <= _CBC_VERSION <= v"2.10.5")
    error(
        "You have installed version $_CBC_VERSION of Cbc, which is not " *
        "supported by Cbc.jl. If the version change was breaking, changes " *
        "will need to be made to the Julia code. Please open an issue at " *
        "https://github.com/jump-dev/Cbc.jl."
    )
end

include("MOI_wrapper/MOI_wrapper.jl")

# TODO(odow): remove at Cbc.jl v1.0.0.
function CbcSolver(args...; kwargs...)
    error(
        "`CbcSolver` is no longer supported. If you are using JuMP, upgrade " *
        "to the latest version and use `Cbc.Optimizer` instead. If you are " *
        "using MathProgBase (e.g., via `lingprog`), you will need to upgrade " *
        "to MathOptInterface (https://github.com/jump-dev/MathOptInterface.jl)."
    )
end
export CbcSolver

end


text:
#==============================================================================
    "more flexible than other languages for object oriented programming"

    tested:     PASS
    author:     discourse.julialang.org/u/synxroform
 =============================================================================#

mutable struct MyStruct
    id::Int64
    name::String
end

function display(this::MyStruct)
    println("=> id:", this.id, " name:", this.name)
end

function merge(this::MyStruct, other::MyStruct)
    this.id += other.id
    this.name *= other.name
    return this
end

function Base.getproperty(this::MyStruct, name::Symbol)
    try
        f = eval(name)
        return (args...) -> f(this, args...)
    catch UndefVarError
        return Base.getfield(this, name)
    end
end

m1 = MyStruct(123, "hello_")
m2 = MyStruct(456, "world")

m1.merge(m2)
m1.display()

#=  RESULTS
MyStruct(123, "hello_")
MyStruct(456, "world")
MyStruct(579, "hello_world")
=> id:579 name:hello_world
=#

max_stars_repo_name: mahiki/julia-is-great

max_issues_repo_name: mahiki/julia-is-great

max_forks_repo_name: mahiki/julia-is-great


text:
"""
```
make_decomp_mbs(m_new, m_old, input_type, cond_new, cond_old, class;
    individual_shocks = false)
```

Construct and return the `MeansBands` (for shockdec, trend, dettrend, hist, and
forecast) necessary to call the plotting function `shockdec` in
`plot_forecast_decomposition`.
"""
function make_decomp_mbs(m_new::M, m_old::M, input_type::Symbol,
                         cond_new::Symbol, cond_old::Symbol,
                         class::Symbol; individual_shocks::Bool = false, forecast_string_old = "", forecast_string_new = "") where M<:AbstractDSGEModel
    # Read in means
    input_file = get_decomp_mean_file(m_new, m_old, input_type, cond_new, cond_old, class, forecast_string_new = forecast_string_new, forecast_string_old = forecast_string_old)
    decomps = JLD2.jldopen(input_file, "r") do file
        read(file, "decomps")
    end

    # Common metadata
    comps = [:data, :news, :para]
    dates = decomps[collect(keys(decomps))[1]][!,:date]
    vars  = collect(keys(get_dict(m_new, class)))

    metadata = Dict{Symbol, Any}()
    metadata[:para]            = input_type
    metadata[:indices]         = get_dict(m_new, class)
    metadata[:class]           = class
    metadata[:date_inds]       = OrderedDict(date => i for (i, date) in enumerate(dates))
    metadata[:forecast_string] = ""
    metadata[:cond_type]       = cond_new

    # Shock decomposition
    shockdec_mb = MeansBands(Dict(metadata), DataFrame(date = dates), Dict{Symbol, DataFrame}())
    shockdec_mb.metadata[:product] = :shockdec
    if individual_shocks
        for var in vars
            shocks = setdiff(names(decomps[var]), vcat([:date, :dettrend, :total], comps))
            if var == vars[1]
                shockdec_mb.metadata[:shock_indices] = OrderedDict(shock => i for (i, shock) in enumerate(shocks))
            end
            for shock in shocks
                varshock = Symbol(var, "__", shock)
                shockdec_mb.means[varshock] = decomps[var][!,shock]
                shockdec_mb.bands[varshock] = DataFrame(date = dates)
            end
        end
    else
        shockdec_mb.metadata[:shock_indices] = OrderedDict(comp => i for (i, comp) in enumerate(comps))
        for var in vars
            for comp in comps
                varcomp = Symbol(var, "__", comp)
                shockdec_mb.means[!,varcomp] = decomps[var][!,comp]
                shockdec_mb.bands[varcomp] = DataFrame(date = dates)
            end
        end
    end

    # Trend
    trend_mb = MeansBands(Dict(metadata), DataFrame(date = dates), Dict{Symbol, DataFrame}())
    trend_mb.metadata[:product] = :trend
    for var in vars
        trend_mb.means[!,var] = zeros(length(dates))
        trend_mb.bands[var] = DataFrame(date = dates)
    end

    # Deterministic trend
    dettrend_mb = MeansBands(Dict(metadata), DataFrame(date = dates), Dict{Symbol, DataFrame}())
    dettrend_mb.metadata[:product] = :dettrend
    for var in vars
        dettrend_mb.means[!,var] = individual_shocks ? decomps[var][:dettrend] : zeros(length(dates))
        dettrend_mb.bands[var] = DataFrame(date = dates)
    end

    # History
    hist_inds  = dates .<= date_mainsample_end(m_new)
    hist_dates = dates[hist_inds]
    if isempty(hist_dates)
        hist_mb = MeansBands()
    else
        hist_mb = MeansBands(Dict(metadata), DataFrame(date = hist_dates), Dict{Symbol, DataFrame}())
        hist_mb.metadata[:product]   = :hist
        hist_mb.metadata[:date_inds] = OrderedDict(date => i for (i, date) in enumerate(hist_dates))
        for var in vars
            hist_mb.means[!,var] = if individual_shocks
                decomps[var][hist_inds, :data] + decomps[var][hist_inds, :news]
            else
                decomps[var][hist_inds, :total]
            end
            hist_mb.bands[var] = DataFrame(date = hist_dates)
        end
    end

    # Forecast
    fcast_inds  = dates .> date_mainsample_end(m_new)
    fcast_dates = dates[fcast_inds]
    if isempty(fcast_dates)
        fcast_mb = MeansBands()
    else
        fcast_mb = MeansBands(Dict(metadata), DataFrame(date = fcast_dates), Dict{Symbol, DataFrame}())
        fcast_mb.metadata[:product]   = :forecast
        fcast_mb.metadata[:date_inds] = OrderedDict(date => i for (i, date) in enumerate(fcast_dates))
        for var in vars
            fcast_mb.means[!,var] = if individual_shocks
                decomps[var][fcast_inds, :data] + decomps[var][fcast_inds, :news]
            else
                decomps[var][fcast_inds, :total]
            end
            fcast_mb.bands[var] = DataFrame(date = fcast_dates)
        end
    end

    return shockdec_mb, trend_mb, dettrend_mb, hist_mb, fcast_mb
end

"""
```
plot_forecast_decomposition(m_new, m_old, var, class, input_type,
    cond_new, cond_old; titles = [], individual_shocks = false,
    groups = shock_groupings(m_new), verbose = :low, kwargs...)

plot_forecast_decomposition(m_new, m_old, vars, class, input_type,
    cond_new, cond_old; titles = [], individual_shocks = false,
    groups = shock_groupings(m_new), verbose = :low, kwargs...)
```

Plot forecast decomposition (looks like a shock decomposition). If
`individual_shocks = false`, then the black and red lines give the total
difference and the bars give the data revision, news, and parameter
components. Otherwise, the black and red lines give the total data revision +
news component and the bars give the individual shock contributions.

The `groups` keyword argument is only used if `individual_shocks = true`.
"""
function plot_forecast_decomposition(m_new::M, m_old::M, var::Symbol, class::Symbol,
                                     input_type::Symbol, cond_new::Symbol, cond_old::Symbol;
                                     title::String = "", kwargs...) where M<:AbstractDSGEModel

    plot_forecast_decomposition(m_new, m_old, [var], class, input_type, cond_new, cond_old;
                                titles = isempty(title) ? String[] : [title], kwargs...)
end

function plot_forecast_decomposition(m_new::M, m_old::M, vars::Vector{Symbol}, class::Symbol,
                                     input_type::Symbol, cond_new::Symbol, cond_old::Symbol;
                                     titles::Vector{String} = String[],
                                     individual_shocks::Bool = false,
                                     groups::Vector{ShockGroup} = shock_groupings(m_new),
                                     plotroot::String = figurespath(m_new, "forecast"),
                                     verbose::Symbol = :low, forecast_string_new = "", forecast_string_old = "",
                                     kwargs...) where M<:AbstractDSGEModel
    # Create MeansBands
    mbs = make_decomp_mbs(m_new, m_old, input_type, cond_new, cond_old, class,
                          individual_shocks = individual_shocks, forecast_string_new = forecast_string_new, forecast_string_old = forecast_string_old)

    # Create shock grouping
    if !individual_shocks
        groups = [ShockGroup("data", [:data], colorant"#9DE0AD"), # sea foam green
                  ShockGroup("news", [:news], colorant"#45ADA8"), # turquoise
                  ShockGroup("para", [:para], colorant"#547980")] # blue gray
    end

    # Get titles if not provided
    if isempty(titles)
        detexify_title = typeof(Plots.backend()) == Plots.GRBackend
        titles = map(var -> describe_series(m_new, var, class, detexify = detexify_title), vars)
    end

    # Loop through variables
    plots = OrderedDict{Symbol, Plots.Plot}()
    for (var, title) in zip(vars, titles)
        # Call recipe
        plots[var] = shockdec(var, mbs..., groups;
                              hist_label = "Historical Diff", forecast_label = "Forecast Diff",
                              ylabel = series_ylabel(m_new, var, class),
                              title = title, kwargs...)

        if !isempty(plotroot)
            # Save plot
            base = Symbol(:decomp, individual_shocks ? :shocks : :total, "_", var)
            output_file = get_decomp_filename(m_new, m_old, input_type, cond_new, cond_old, base, Symbol(),
                                              pathfcn = figurespath, fileformat = plot_extension(), forecast_string_new = forecast_string_new, forecast_string_old = forecast_string_old)
            output_file = joinpath(plotroot, basename(output_file))

            save_plot(plots[var], output_file, verbose = verbose)
        end
    end
end


text:
using FileIO;
using DelimitedFiles;
using Statistics: mean, std;
using Random: randperm;

# =============================================================================
# datasets.jl -> Funciones útiles aplicables a un dataset:
#   - oneHotEncoding
#   - normalizar:
#       -calculateMinMaxNormalizationParameters
#       -normalizeMinMax!
#       -calculateZeroMeanNormalizationParameters
#       -normalizeZeroMean!
#   - holdOut (sobrecargada 2 y 3 params)
#   - validacion cruzada(añadir en un futuro) -> FALTAS TU
#   - classifyOutputs
# =============================================================================

#
# Esta función sirve para normalizar las salidas deseadas del dataset para un
#   problema de clasificacion, tanto como si es binario o multiclase.
#
# @arguments
#   feature: Array con los outputs deseados
#
# @return: outputs en one-hot-encoding
#
function oneHotEncoding(feature::Array{Any,1}, classes::Array{Any,1})
    # Primero se comprueba que todos los elementos del vector esten en el vector de clases (linea adaptada del final de la practica 4)
    @assert(all([in(value, classes) for value in feature]));
    numClasses = length(classes);

    @assert(numClasses>1)
    if (numClasses==2)
        # Si solo hay dos clases, se devuelve una matriz con una columna
        oneHot = Array{Bool,2}(undef, size(feature,1), 1);
        oneHot[:,1] .= (feature.==classes[1]);
    else
        # Si hay mas de dos clases se devuelve una matriz con una columna por clase
        oneHot = Array{Bool,2}(undef, size(feature,1), numClasses);
        for numClass = 1:numClasses
            oneHot[:,numClass] .= (feature.==classes[numClass]);
        end;
    end;
    return oneHot;
end;
# Extraemos las clases distintas con unique
oneHotEncoding(feature::Array{Any,1}) = oneHotEncoding(feature::Array{Any,1}, unique(feature));
#  En este caso, el propio vector ya está codificado
oneHotEncoding(feature::Array{Bool,1}) = feature;

# =============================================================================
# Funciones útiles para la normalización de un dataset
# =============================================================================

#
# Función auxiliar para calcular el máximo y mínimo de cada atributo.
#
# @arguments
#   dataset: los inputs.
#   dataInRows: si esta traspuesta o no.
#
# @return: array con el máximo y mínimo de cada atributo.
#
calculateMinMaxNormalizationParameters(dataset::Array{Float64,2}; dataInRows=true) =
    ( minimum(dataset, dims=(dataInRows ? 1 : 2)), maximum(dataset, dims=(dataInRows ? 1 : 2)) );

#
# Función AUXILIAR para normalizar los inputs entre máximo y mínimo.
#
# @arguments
#   dataset: los inputs.
#   normalizationParameters: array con el máximo y mínimo de cada atributo.
#   dataInRows: si esta traspuesta o no.
#
function normalizeMinMax!(dataset::Array{Float64,2}, normalizationParameters::Tuple{Array{Float64,2},Array{Float64,2}}; dataInRows=true)
    min = normalizationParameters[1];
    max = normalizationParameters[2];
    dataset .-= min;
    dataset ./= (max .- min);
    # Si hay algun atributo en el que todos los valores son iguales, se pone a 0
    if (dataInRows)
        dataset[:, vec(min.==max)] .= 0;
    else
        dataset[vec(min.==max), :] .= 0;
    end;
end;

#
# Función que normaliza los inputs entre máximo y mínimo (es a la que llamamos).
#
# @arguments
#   dataset: los inputs.
#   dataInRows: si la matriz está traspuesta o no.
#
normalizeMinMax!(dataset::Array{Float64,2}; dataInRows=true) = normalizeMinMax!(dataset, calculateMinMaxNormalizationParameters(dataset; dataInRows=dataInRows); dataInRows=dataInRows);

#
# Función auxiliar para calcular la media y la desviacion típica de cada atributo.
#
# @arguments
#   dataset: los inputs.
#   dataInRows: si esta traspuesta o no.
#
# @return: array con la media y la desviacion típica de cada atributo.
#
calculateZeroMeanNormalizationParameters(dataset::Array{Float64,2}; dataInRows=true) =
( mean(dataset, dims=(dataInRows ? 1 : 2)), std(dataset, dims=(dataInRows ? 1 : 2)) );

#
# Función AUXILIAR para normalizar los inputs con media cero.
#
# @arguments
#   dataset: los inputs.
#   normalizationParameters: array con la media y la desviacion típica de cada atributo.
#   dataInRows: si esta traspuesta o no.
#
function normalizeZeroMean!(dataset::Array{Float64,2}, normalizationParameters::Tuple{Array{Float64,2},Array{Float64,2}}; dataInRows=true)
    avg  = normalizationParameters[1];
    stnd = normalizationParameters[2];
    dataset .-= avg;
    dataset ./= stnd;
    # Si hay algun atributo en el que todos los valores son iguales, se pone a 0
    if (dataInRows)
        dataset[:, vec(stnd.==0)] .= 0;
    else
        dataset[vec(stnd.==0), :] .= 0;
    end
end;

#
# Función que normaliza los inputs con media 0 (es a la que llamamos).
#
# @arguments
#   dataset: los inputs.
#   dataInRows: si la matriz está traspuesta o no.
#
normalizeZeroMean!(dataset::Array{Float64,2}; dataInRows=true) = normalizeZeroMean!(dataset, calculateZeroMeanNormalizationParameters(dataset; dataInRows=dataInRows); dataInRows=dataInRows);

# =============================================================================
# Obtener los conjuntos discuntos adecuados (holdOut)
# =============================================================================

#
# Función que divide el conjunto de entrenamiento en entrenamiento y test.
#
# @arguments
#   N: número de patrones.
#   P: porcentaje de división test/entrenamiento.
#
# @return nuevos conjuntos de entrenamiento y test
#
function holdOut(N::Int, P::Float64)
    @assert ((P>=0.) & (P<=1.));
    indices = randperm(N);
    numTrainingInstances = Int(round(N*(1-P)));
    return (indices[1:numTrainingInstances], indices[numTrainingInstances+1:end]);
end

#
# Función que divide el conjunto de entrenamiento en entrenamiento, test y validación.
#
# @arguments
#   N: número de patrones.
#   Ptest: porcentaje de división test/entrenamiento.
#   Pval: porcentaje de división validación/entrenamiento.
#
# @return nuevos conjuntos de entrenamiento, test y validacion.
#
function holdOut(N::Int, Ptest::Float64, Pval::Float64)
    @assert ((Pval>=0.) & (Pval<=1.));
    @assert ((Ptest>=0.) & (Ptest<=1.));
    @assert ((Pval+Ptest)<=1.);
    # Primero separamos en entrenamiento+validation y test
    (trainingValidationIndices, testIndices) = holdOut(N, Ptest);
    # Después separamos el conjunto de entrenamiento+validation
    (trainingIndices, validationIndices) = holdOut(length(trainingValidationIndices), Pval*N/length(trainingValidationIndices))
    return (trainingValidationIndices[trainingIndices], trainingValidationIndices[validationIndices], testIndices);
end;

# =============================================================================
# classifyOutputs -> Pasar de una salida en float64 a bool
# =============================================================================

function classifyOutputs(outputs::Array{Float64,2}; dataInRows::Bool=true, threshold::Float64=0.5)
    numOutputs = size(outputs, dataInRows ? 2 : 1);
    @assert(numOutputs!=2)
    if numOutputs==1
        return convert(Array{Bool,2}, outputs.>=threshold);
    else
        # Miramos donde esta el valor mayor de cada instancia con la funcion findmax
        (_,indicesMaxEachInstance) = findmax(outputs, dims= dataInRows ? 2 : 1);
        # Creamos la matriz de valores booleanos con valores inicialmente a false y asignamos esos indices a true
        outputsBoolean = Array{Bool,2}(falses(size(outputs)));
        outputsBoolean[indicesMaxEachInstance] .= true;
        # Comprobamos que efectivamente cada patron solo este clasificado en una clase
        @assert(all(sum(outputsBoolean, dims=dataInRows ? 2 : 1).==1));
        return outputsBoolean;
    end;
end;


text:
function makealltrans(N,n,Ω,basis="Hermite")
dim=length(N)
if dim==1
  Nx = N[1]
  ωx = Ω[1]
  #n-field transforms for PGPE
  x,wx,Tx = nfieldtrans(Nx,n,ω=ωx,basis=basis)
  return x,wx,Tx
elseif dim==2
  Nx,Ny = N
  ωx,ωy = Ω
  #n-field transforms for PGPE
  x,wx,Tx = nfieldtrans(Nx,n,ω=ωx,basis=basis)
  y,wy,Ty = nfieldtrans(Ny,n,ω=ωy,basis=basis)
  return x,wx,Tx,y,wy,Ty
elseif dim==3
  Nx,Ny,Nz = N
  ωx,ωy,ωz = Ω
  #n-field transforms for PGPE
  x,wx,Tx = nfieldtrans(Nx,n,ω=ωx,basis=basis)
  y,wy,Ty = nfieldtrans(Ny,n,ω=ωy,basis=basis)
  z,wz,Tz = nfieldtrans(Nz,n,ω=ωz,basis=basis)
  return x,wx,Tx,y,wy,Ty,z,wz,Tz
end
end


text:
using JGUI
using Base.Test

### Simple usage
w = window()
gr = grid(w)
push!(w, gr)

b11 = button(gr, "11")
b22 = button(gr, "22")
b33 = button(gr, "33")

## basic array-like interface to add widget
gr[1,1] = b11
gr[2,2] = b22
gr[3,3] = b33

raise(w)

### alignment
w = window()
gr = grid(w)
push!(w, gr)

b11 = button(gr, "1,1")
b11[:alignment] = (:left, :top)
b11[:sizepolicy] = (:expand, :expand)

b22 = button(gr, "2,2")
b33 = button(gr, "3,3")
b33[:alignment] = (:right, :bottom)

gr[1,1] = b11
gr[2,2] = b22
gr[3,3] = b33

### can add through an array
using JGUI
w = window()
gr = grid(w)
push!(w, gr)

b11 = button(gr, "1,1")
b11[:alignment] = (:left, :top)
b11[:sizepolicy] = (:expand, :expand)

b22 = button(gr, "2,2")
b33 = button(gr, "3,3")
b33[:alignment] = (:right, :bottom)

gr[:,:] = [b11 b22; nothing b33]

raise(w)
### spanning multiple columns
w = window()
gr = grid(w)
push!(w, gr)

b11 = button(gr, "1,1")
b22 = button(gr, "2,2")
b33 = button(gr, "3,3")
[btn[:sizepolicy] = (:expand, :expand) for btn in [b11, b22, b33]]

gr[1:2, 1] = b11
gr[1,3] = b22
gr[2, 2:3] = b33
raise(w)
### configure weights -- abuse of notation

w = window()
gr = grid(w)
push!(w, gr)

b11 = button(gr, "1,1")
b22 = button(gr, "2,2")
b33 = button(gr, "3,3")
[btn[:sizepolicy] = (:expand, :expand) for btn in [b11, b22, b33]]

gr[1,:] = [b11 b22 b33]
column_stretch(gr, 1, 1)
column_stretch(gr, 2, 2)
column_stretch(gr, 3, 3)
raise(w)


text:
export isean

function getweigth(len, index)
    if len === 14 || len === 8
        return index % 2 === 0 ? 1 : 3
    end

    return index % 2 === 0 ? 3 : 1
end

function isean(str::AbstractString)::Bool
    eanReg = r"^(\d{8}|\d{13}|\d{14})$"

    if !contains(str, eanReg)
        return false
    end

    checkDigit = parse(Int, str[end])
    checksum = sum(map(((index, value), ) -> value * getweigth(length(str), index) , enumerate(parse.(Int8, split(str[1:end-1], "")))))
    remainder = 10 - (checksum % 10)
    expectedCheckDigit = remainder < 10 ? remainder : 0

    return checkDigit === expectedCheckDigit
end


text:
defaultsamplesregion() = ((60, 60, 60), Rect(Vec(-2,-2,-2), Vec(4,4,4)))
samples, region = defaultsamplesregion()

function gyroid()
    gyroid(v) = cos(v[1])*sin(v[2])+cos(v[2])*sin(v[3])+cos(v[3])*sin(v[1])
    gyroid_shell(v) = max(gyroid(v)-0.4,-gyroid(v)-0.4)

    # generate directly using GeometryBasics API
    # Rect specifies the sampling intervals
    gy_mesh = Mesh(gyroid_shell, Rect(Vec(0,0,0),Vec(pi*4,pi*4,pi*4)),
                    MarchingCubes(), samples=(50,50,50))

    return gy_mesh
end

spheredist(v) = norm(v) - 1.0

function boxdist(v, b)
    q = abs.(v) - b
    return norm(max.(q, 0)) + min.(maximum(q), 0)
end

function sphere()

    # generate directly using GeometryBasics API
    # Rect specifies the sampling intervals
    sphere_mesh = Mesh(spheredist, Rect(Vec(-1.2,-1.2,-1.2),Vec(2.4, 2.4, 2.4)),
                    MarchingCubes(), samples=(10,10,10))

    return sphere_mesh
end

D(f) = x -> ForwardDiff.derivative(f, x)

linecurve(t) = Point3{Float64}(1., 1., 1.) .* t
function curvedistgen(curve::Function, mint, maxt; eps=1e-13)
    function curvedist(p)
        disttocurve(t) = normsq(curve(t) - p)
        Ddisttocurve = D(disttocurve)
        return sqrt(minimum(disttocurve.(find_zeros(Ddisttocurve, mint, maxt))))
    end
    return curvedist
end

function arc_length_parameterization(γ::Function, tspan)
    normtangentγ(t) = norm(ForwardDiff.derivative(γ, t))

    f_t_to_s(s, p, t) = normtangentγ(t)
    s0 = 0
    t_to_s_prob = ODEProblem(f_t_to_s, s0, tspan)
    t_to_s = solve(t_to_s_prob, Tsit5(), reltol=1e-8, abstol=1e-8)

    maxs = t_to_s(1.0)

    f_s_to_t(t, p, s) = 1 / normtangentγ(t)
    t0 = tspan[1]
    sspan = (0.0, maxs)
    s_to_t_prob = ODEProblem(f_s_to_t, t0, sspan)
    s_to_t = solve(s_to_t_prob, Tsit5(), reltol=1e-8, abstol=1e-8)

    reparamγ(s) = γ(s_to_t(s))
    t_to_s, s_to_t, reparamγ, (0., maxs)
end

function curveofshapes(curve::AbstractArray{Tuple{Point3{Float64}, Float64}, 1}, shapefunc::Function)
    function innerdistance(v)
        distances = map(p->shapefunc((v - p[1]) .* (1 / p[2])), curve)
        return minimum(distances)
    end
    return innerdistance
end

function curvedistgenkdtree(γ::Function, mint, maxt; arc_length_sampling_freq=0.01)
    tspan = (mint, maxt)

    t_to_s, s_to_t, γs, sspan = arc_length_parameterization(γ, tspan)
    mins, maxs = sspan[1:2]

    γsamples = (mins:arc_length_sampling_freq:maxs)
    γsampled_points = map(γs, γsamples)

    γkdtree = KDTree(γsampled_points)

    function curvedist(p)
        #disttocurvet(t) = normsq(γ(t) - p)
        #Ddisttocurvet = D(disttocurvet)
        #Ddisttocurves = Ddisttocurvet ∘ s_to_t

        s_nearest_index, distance = knn(γkdtree, p, 1)
        #s_nearest = γsamples[s_nearest_index][1]
        #t_nearest = s_to_t(s_nearest)
        #s_search_range = (s_nearest - arc_length_sampling_freq, s_nearest + arc_length_sampling_freq)
        #@show s_search_range
        #try
            #t_min = find_zero(Ddisttocurvet, t_nearest; maxevals=100, atol=1e-4, rtol=1e-3)
            #return sqrt(disttocurvet(t_nearest))
        distance[1]
            #=
        catch e
            @show t_nearest, s_nearest, s_nearest_index
            @show p
            @show disttocurvet(t_nearest)
            @show e
            throw(e)
        end
        =#

    end
    curvedist, γsampled_points
end

struct Curve
    γ::Function
    mint::Float64
    maxt::Float64
end

function curvesdistgenkdtree(γs::Vector{Curve}; arc_length_sampling_freq=0.01)
    tspan = (mint, maxt)

    t_to_s, s_to_t, γs, sspan = arc_length_parameterization(γ, tspan)
    mins, maxs = sspan[1:2]

    γsamples = (mins:arc_length_sampling_freq:maxs)
    γsampled_points = map(γs, γsamples)

    γkdtree = KDTree(γsampled_points)

    function curvedist(p)
        #disttocurvet(t) = normsq(γ(t) - p)
        #Ddisttocurvet = D(disttocurvet)
        #Ddisttocurves = Ddisttocurvet ∘ s_to_t

        s_nearest_index, distance = knn(γkdtree, p, 1)
        #s_nearest = γsamples[s_nearest_index][1]
        #t_nearest = s_to_t(s_nearest)
        #s_search_range = (s_nearest - arc_length_sampling_freq, s_nearest + arc_length_sampling_freq)
        #@show s_search_range
        #try
            #t_min = find_zero(Ddisttocurvet, t_nearest; maxevals=100, atol=1e-4, rtol=1e-3)
            #return sqrt(disttocurvet(t_nearest))
        distance[1]
            #=
        catch e
            @show t_nearest, s_nearest, s_nearest_index
            @show p
            @show disttocurvet(t_nearest)
            @show e
            throw(e)
        end
        =#

    end
    curvedist, γsampled_points
end

function boundary_triangles(tetrahedra::Array{SimplexFace{4, Int32}, 1})
    tritonum = Dict{SVector{3, UInt32}, Int}()
    @show tritonum
    #tritonum[sv(1, 2, 3)] = 1
    @show tritonum

    tricombinations = [sv(1,2,3), sv(1,2,4), sv(1,3,4), sv(2,3,4)]
    @show tricombinations
    @show tricombinations[1]
    alltris = Array{NgonFace{3, UInt32}, 1}(undef, 0)
    for (i, tet) in enumerate(tetrahedra)
        for combination in tricombinations
            index = sortsv(tet[combination]...)
            if !haskey(tritonum, index)
                tritonum[index] = 1
            else
                tritonum[index] += 1
            end
            tri = NgonFace{3, UInt32}(index...)
            push!(alltris, tri)
        end
    end
    @show length(tritonum)

    boundarytris = Array{NgonFace{3, UInt32}, 1}(undef, 0)
    for (key, val) in pairs(tritonum)
        if val == one(UInt32)
            tri = NgonFace{3, UInt32}(key...)
            push!(boundarytris, tri)
        end
    end

    return boundarytris, alltris
end

function rectGBtoGT(rect::GeometryBasics.HyperRectangle{N, T}) where {N, T <: Number}
    return GeometryTypes.HyperRectangle{N, T}(origin(rect), widths(rect))
end

function rectGTtoGB(rect::GeometryTypes.HyperRectangle{N, T}) where {N, T <: Number}
    return GeometryBasics.HyperRectangle{N, T}(rect.origin, rect.widths)
end

function origin(rect::GeometryTypes.HyperRectangle{N, T}) where {N, T <: Number}
    return rect.origin
end

function widths(rect::GeometryTypes.HyperRectangle{N, T}) where {N, T <: Number}
    return rect.widths
end

function spheresdf(samples::NTuple{3, T}, region::Rect) where {T <: Integer}
    regionsamples = sampleregion(samples, region)
    return spheredist.(regionsamples)
end

function sampleregion(samples::NTuple{3, T}, region::Rect) where {T <: Integer}
    regionsamples = zeros(Point3{Float64}, samples[1], samples[2], samples[3])
    regionwidths = widths(region)
    regionorigin = origin(region)
    regionincrements = regionwidths ./ (samples .- 1)
    for i = 1:samples[1], j = 1:samples[2], k = 1:samples[3]
        samplepositions = ((i, j, k) .- 1) .* regionincrements .+ regionorigin
        regionsamples[i, j, k] = Point3{Float64}(samplepositions...)
    end
    return regionsamples
end

function facearray(mesh)
    mesh_faces = faces(mesh)
    m = length(mesh_faces)
    face_arr = Array{NgonFace{3, Int64}, 1}(undef, m)
    for i = 1:m
        face_arr[i] = Int64.(mesh_faces[i])
    end
    return face_arr
end

function singletri()
    verts = Array{Point{3, Float64}, 1}(undef, 3)
    pointtype = Point{3, Float64}
    verts[1] = pointtype([0.0, 0.0, 1.0])
    verts[2] = pointtype([0.0, 1.0, 0.0])
    verts[3] = pointtype([1.0, 0.0, 0.0])
    faces = Array{GLTriangleFace, 1}(undef, 1)
    faces[1] = GLTriangleFace(1, 2, 3)
    geom_mesh = Mesh(verts, faces)
    return geom_mesh
end

function connectedvertices(faces::AbstractArray{NgonFace{3, N}}) where {N <: Integer}
    Ntype = eltype(N)
    connectedvertexsets = Vector{Set{Ntype}}(undef, 0)
    #@show connectedvertexsets
    #@show length(connectedvertexsets)
    for i =  1:length(faces)
    #for i =  1:10
        if i % 100 == 0
            @show i
        end
        #@show i
        curset = Set{Ntype}([Ntype(faces[i][v]) for v = 1:3])
        #@show curset
        #@show typeof(curset)
        for v = 1:3
            for k = length(connectedvertexsets):-1:1
            #@show k
                vert = Ntype(faces[i][v])
                if vert in connectedvertexsets[k]
                    #@show typeof(curset)
                    #@show connectedvertexsets[k]
                    vertexset = popat!(connectedvertexsets, k)
                    #@show vertexset
                    curset = union!(vertexset, curset)
                    #@show typeof(curset)
                end
            end
        end
        #@show typeof(connectedvertexsets)
        push!(connectedvertexsets, curset)
    end
    return connectedvertexsets
end

function removesmallcomponents(mesh, connectedvertexsets; sizetokeep=1)
    sizessorted = sort(length.(connectedvertexsets))
    componentsizetokeep = sizessorted[sizetokeep]
    biggestcomponent = collect(filter(x -> length(x) == componentsizetokeep, connectedvertexsets)[1])
    oldtonew = Dict{UInt32, UInt32}()
    numnewverts = length(biggestcomponent)
    newverts = Array{Point3{Float64}, 1}(undef, numnewverts)
    oldverts = coordinates(mesh)
    #for i = 1:4
    for i = 1:numnewverts
        oldindex = biggestcomponent[i]
        oldtonew[oldindex] = i
        #@show oldverts[oldindex]
        newverts[i] = oldverts[oldindex]
    end
    @show length(newverts)

    oldfaces = faces(mesh)
    newfaces = Array{NgonFace{3, UInt32}, 1}(undef, 0)
    @show newfaces
    for i = 1:length(oldfaces)
    #for i = 1:100
        #@show i
        #@show oldfaces[i]
        newindices = Array{UInt32, 1}(undef, 0)
        for j = 1:3
            oldindex = UInt32(oldfaces[i][j])
            #@show oldindex
            if haskey(oldtonew, oldindex)
                push!(newindices, oldtonew[oldindex])
            end
        end
        #@show newindices
        if length(newindices) == 3
            newface = NgonFace{3, UInt32}(newindices...)
            push!(newfaces, newface)
        end
        #newindices = [oldtonew[] for j = 1:3]
        #@show newindices
        #newface = NgonFace{3, UInt32}([])

    end
    #@show newfaces
    newmesh = GeometryBasics.Mesh(newverts, newfaces)
    return newmesh
    #return biggestcomponent

end

function mandelbulbmesh(samples::NTuple{3, T}, region::Rect, numfractaliter::Integer) where {T <: Integer}
    mandelbulbdist = mandelbulbdistfunc(numfractaliter)
    return Mesh(mandelbulbdist, region, samples, MarchingTetrahedra(), pointtype=Point3{Float64})
end

function meshview(inmesh, color)
    return Makie.mesh(inmesh, color=color, show_axis=false)
end

function meshview(inmesh)
    return Makie.mesh(inmesh, color=[norm(v) for v in coordinates(inmesh)], show_axis=false)
end

function normsq(p::AbstractArray{T}) where {T <: Number}
    return sum(p .* p)
end

function inversesqrt(num::T) where {T <: Number}
    return 1.0 / sqrt(num)
end

function mandelbulbdistfunc(numfractaliter::Integer)
    function mandelbulbdist(p::Point{3, T}) where {T <: Real}
        px, py, pz = p[1], p[2], p[3]
        dz = 1.0
        w = Array{T, 1}(undef, 3)
        w .= 0

        stop_iter = 0

        w .+= p
        #@show typeof(w)
        m = normsq(w)
        #trap = vcat(abs.(w), m)

        for i = 1:numfractaliter

            w2 = w .* w
            w4 = w2 .* w2

            x, y, z = w[1], w[2], w[3]
            x2, y2, z2 = w2[1], w2[2], w2[3]
            x4, y4, z4 = w4[1], w4[2], w4[3]

            m2 = m * m
            m4 = m2 * m2

            dz *= 8.0 * sqrt(m4 * m2 * m)
            dz += 1.0

            k3 = x2 + z2
            k2 = inversesqrt(k3 * k3 * k3 * k3 * k3 * k3 * k3)
            k1 = x4 + y4 + z4 -
                6.0 * y2 * z2 -
                6.0 * x2 * y2 +
                2.0 * z2 * x2
            k4 = x2 - y2 + z2

            neww = [
                px +
                64.0 * x * y * z * (x2 - z2) * k4 *
                (x4 - 6.0 * x2 * z2 + z4) * k1 * k2,
                py +
                -16.0 * y2 * k3 * k4 * k4 + k1 * k1,
                pz +
                -8.0 * y * k4 * (x4 * x4 -
                28.0 * x4 * x2 * z2 +
                70.0 * x4 * z4 -
                28.0 * x2 * z2 * z4 +
                z4 * z4) * k1 * k2
                ]
            w[:] = neww[:]

            #trap[:, :, :] = min.(trap, vcat(abs.(w), m)) .* (1.0 .- stop_iter) .+ stop_iter .* trap
            m = normsq(w)
            if m > 256.0
                return -0.0005
            end
            #m[:, :, :] = min.(normsq(w), 256.0)
            #m[map(isnan, m)] .= 256.0
            #dz[map(isnan, dz)] .= 1.0

            #@show i
            #@show w
            #@show m
            #@show stop_iter
            #@show dz
            stop_iter = i

        end

        dist = 0.25 * log(m) * sqrt(m) / dz
        return dist
    end
    return mandelbulbdist
end

function init()
    blob2 = mandelbulbmesh(samples, region, 2)
    connecteds = connectedvertices(faces(blob2))
    #newblob2 = removesmallcomponents(blob2, connecteds)
    return blob2, connecteds #, newblob2
end

function loadbunny()
    bunnymeshpre = load("meshes/original/small_bunny.ply")
    vertspre, tris = coordinates(bunnymeshpre), faces(bunnymeshpre)
    transformmat = [1 0 0; 0 0 1; 0 1 0]
    verts = map(vert->Point3{Float32}(transformmat * vert), vertspre)

    return Mesh(verts, tris)
end

function bunnyeyepot(t::Float64)
    eyepot = Array(sparsevec(
        # eyeleft, eyeright, max size, booty, ear left, ear right
        [1311, 1015, 1430, 680, 568, 747],
         [Float64(3 * cos(t)), Float64(3 * cos(t)), Float64(0), Float64(-2), Float64(0.5), Float64(0.5)]))
    return eyepot
end

#begin
    #p1 = solvepoissonproblem(bs, bunnyeyepot(0.))
    #for i in 1:100
        #p1 = solveheatstep(bs, p1, 0.04)
        #meshview(bunny, p1)
    #end
#end

max_stars_repo_name: zoemcc/geometry

max_issues_repo_name: zoemcc/geometry

max_forks_repo_name: zoemcc/geometry


text:
@DeepShiba_Func mutable struct Log <: Func end

@DeepShiba_Func mutable struct Sin <: Func end

@DeepShiba_Func mutable struct Cos <: Func end

@DeepShiba_Func mutable struct Tan <: Func end

@DeepShiba_Func mutable struct Tanh <: Func end

_Log(x::Variable) = Log([x, ], nothing, 0)(x)
_Sin(x) = Sin([x, ],   nothing, 0)(x)
_Cos(x) = Cos([x, ], nothing, 0)(x)
_Tan(x) = Tan([x, ], nothing, 0)(x)
_Tanh(x) = Tanh([x, ], nothing, 0)(x)

forward(f::Log, x) = log(x)
forward(f::Sin, x) = sin(x)
forward(f::Cos, x) = cos(x)
forward(f::Tan, x) = tan(x)
forward(f::Tanh, x) = tanh(x)

function backward(f::Log)::AbstractArray{Variable}
    gys = get_gys(f.outputs)
    x = f.inputs[1]
    return [1 / x,] .* gys
end

function backward(f::Sin)::AbstractArray{Variable}
    gys = get_gys(f.outputs)
    x, = f.inputs
    return [cos(x),] .* gys
end

function backward(f::Cos)::AbstractArray{Variable}
    gys = get_gys(f.outputs)
    x, = f.inputs
    return [-sin(x),] .* gys
end

function backward(f::Tan)::AbstractArray{Variable}
    gys = get_gys(f.outputs)
    x, = f.inputs
    return [1 / cos(x)^2,] .* gys
end

function backward(f::Tanh)::AbstractArray{Variable}
    gys = get_gys(f.outputs)
    y = f.outputs[1]
    return [1 - y^2,] .* gys
end


text:
@testset "largestconncomp on dense bitmap" begin
    pic1 = falses(150,150);
    pic1[55:75, 55:75] = trues(21,21);
    pic1[100:125, 100:125] = trues(26,26);
    pic1[130:140,130:140] = trues(11,11);
    indmap = [Int[] for i in 1:150, j in 1:150];
    for i in 55:75, j in 55:75
        append!(indmap[i,j], [1,2,3])
    end
    for i in 100:125, j in 100:125
        append!(indmap[i,j], [-99,-98])
    end
    for i in 130:140, j in 130:140
        append!(indmap[i,j], [0])
    end

    t1 = RANSAC.largestconncomp(pic1, indmap, 1:2)
    @test t1 == repeat([-99, -98], 26*26)
    @test t1 == RANSAC.largestconncomp(pic1, indmap)
    @test t1 == RANSAC.largestconncomp(pic1, indmap, trues(3,3))
end

@testset "8 connectivity" begin
    # largest patch is only 8 connectivity connected
    pic2 = falses(150,150);
    pic2[55:73, 55:73] = trues(19,19);
    pic2[100:2:126, 100:2:126] = trues(14,14);
    pic2[101:2:127, 101:2:127] = trues(14,14);
    pic2[130:140,130:140] = trues(11,11);
    indmap = [Int[] for i in 1:150, j in 1:150];
    for i in 55:73, j in 55:73
        append!(indmap[i,j], [1,2,3])
    end
    for i in 100:2:126, j in 100:2:126
        append!(indmap[i,j], [-99,-98])
    end
    for i in 101:2:127, j in 101:2:127
        append!(indmap[i,j], [-99,-98])
    end
    for i in 130:140, j in 130:140
        append!(indmap[i,j], [0])
    end

    t2d = RANSAC.largestconncomp(pic2, indmap)
    t2def = RANSAC.largestconncomp(pic2, indmap, :default)
    t24 = RANSAC.largestconncomp(pic2, indmap, 1:2)
    @test t24 == repeat([1,2,3], 19*19)
    @test t2d == t24
    @test t2def == t24

    t28 = RANSAC.largestconncomp(pic2, indmap, trues(3,3))
    t28eight = RANSAC.largestconncomp(pic2, indmap, :eight)
    @test t28 == repeat([-99,-98], 14*14*2)
    @test t28 == t28eight
end


text:
function test_modes()
  @testset "--- Testing Sparse/Dense Jacobian ---" begin

    fname = "input_vals_vortex1.jl"
    println("\n\ntesting ", fname)
  #  fname = "input_vals_vortex5.jl"
    mesh, sbp, eqn, opts = run_solver(fname)

    @test  calcNorm(eqn, eqn.res_vec)  < 1e-9

    println("\n\n----- Testing jacobian vector product -----")
    # test jacobian vector product
    # load a new initial condition
    ICFunc = EulerEquationMod.ICDict["ICIsentropicVortex"]
    ICFunc(mesh, sbp, eqn, opts, eqn.q_vec)
    NonlinearSolvers.array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec)
#=  # unneeded after LinearOperator introduction
    # compute jacobian
    Tsol = eltype(eqn.q)
    Tres = eltype(eqn.res)
    jac = SparseMatrixCSC(mesh, typeof(real(eqn.res[1])))
    epsilon = 1e-20
    pert = complex(0, epsilon)
    newton_data = NonlinearSolvers.NewtonData{Tsol, Tres}(mesh, sbp, eqn, opts)
    NonlinearSolvers.calcJacobianSparse(mesh, sbp, eqn, opts, EulerEquationMod.evalResidual, Array{Float64}( 0,0,0), pert, jac)

    # to mat-vec product
    newton_data = NonlinearSolvers.NewtonData{Tsol, Tres}(mesh, sbp, eqn, opts)
    v = ones(mesh.numDof)  # vector to multiply jacobian against
    result1 = jac*v
    result2 = zeros(mesh.numDof)

    NonlinearSolvers.calcJacVecProd(mesh, sbp, eqn, opts, pert, physicsRhs, (EulerEquationMod.evalResidual,), v, result2)

    # check the two products are equal
    for i=1:mesh.numDof
      @test isapprox( result1[i], result2[i]) atol=1e-13
    end
=#
  #  results_all = hcat(result1, result2)
  #  println("results_all = \n", results_all)

  #  results_diff = result1 - result2
  #  diff_norm = calcNorm(eqn, results_diff)
  #  println("diff_norm = ", diff_norm)

   #=
    # test entropy variables
    fname = "input_vals_vortexa.jl"
    println("\n\ntesting ", fname)
    mesh, sbp, eqn, opts = solvePDE(fname)

    @test  calcNorm(eqn, eqn.res_vec)  < 1e-9
    =#

    fname = "input_vals_vortex2.jl"
    println("\n\ntesting ", fname)
    mesh, sbp, eqn, opts = solvePDE(fname)
    @test  calcNorm(eqn, eqn.res_vec)  < 1e-9

  #=
    # test entropy variables
    fname = "input_vals_vortex2a.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)
    @test  calcNorm(eqn, eqn.res_vec)  < 1e-9
  =#

    fname = "input_vals_vortex3dg.jl"
    println("\n\ntesting ", fname)
    mesh, sbp, eqn, opts = solvePDE(fname)
    @test  calcNorm(eqn, eqn.res_vec)  < 1e-9

  #=
    # test entropy variables
    fname = "input_vals_vortex3a.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)
    @test  calcNorm(eqn, eqn.res_vec)  < 1e-9
  =#

    fname = "input_vals_vortex4.jl"
    println("\n\ntesting ", fname)
    mesh, sbp, eqn, opts = solvePDE(fname)
    @test  calcNorm(eqn, eqn.res_vec)  < 1e-9

  #=
    # test entropy variables
    fname = "input_vals_vortex4a.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)
    @test  calcNorm(eqn, eqn.res_vec)  < 1e-9
  =#
  end

  return nothing
end

add_func1!(EulerTests, test_modes, [TAG_SHORTTEST])


text:
#=====================Test different annotator pipelines.======================#
#Dependencies
    using JavaCoreNLP
    include("./data/texts.jl")
    #using Base.Test

#UNCOMMENT ONE TO SELECT AN ANNOTATION PIPELINE.
    ##pipeline = StanfordCoreNLP("tokenize")                #ok
    ##pipeline = StanfordCoreNLP("quote")                   #Runs to completion.
    ##pipeline = StanfordCoreNLP("sentiment")               #Runs to completion.
    ##pipeline = StanfordCoreNLP("ssplit")                  #ok
    ##pipeline = StanfordCoreNLP("pos")                     #ok
    ##pipeline = StanfordCoreNLP("lemma")                   #ok
    pipeline = StanfordCoreNLP("parse")                   #ok
    ##pipeline = StanfordCoreNLP("ner")                     #ok
    ##pipeline = StanfordCoreNLP("depparse")                #ok
    ##pipeline = StanfordCoreNLP("coref")                   #ok
    ##pipeline = StanfordCoreNLP("relation")                  #ok
    ##pipeline = StanfordCoreNLP("mention")                 #ok
    ##pipeline = StanfordCoreNLP("ner", "coref", "openie")  #ok

##UNCOMMENT ONE TO SELECT A TEXT.
    text = hitchhikers_text
    ##text = obama_text
    ##text = google_text
    ##text = happy_text
    ##text = subset_text_1
    ##text = instance_text_1

#Main function call.
    demo(pipeline, text)

#Report time consumed by each annotator in the pipeline.
    println("\n", timinginformation(pipeline))

#Report done.
    println("\nComplete.")


text:
module layout

import Parameters: @with_kw

import Term: int, get_last_valid_str_idx
import ..renderables: RenderablesUnion, Renderable, AbstractRenderable
import ..measure: Measure
import ..segment: Segment
using ..box
import ..consoles: console_width, console_height

export Padding, vstack, hstack, pad
export Spacer, vLine, hLine

# ---------------------------------------------------------------------------- #
#                                    PADDING                                   #
# ---------------------------------------------------------------------------- #
"""
    Padding

Stores information about ammount of padding.
"""
@with_kw struct Padding
    left::Int
    right::Int
    top::Int
    bottom::Int
end

"""
    pad(text::AbstractString, target_width::Int, method::Symbol)::String

Pad a string to width: `target_width` by adding empty spaces strings " ".
Where the spaces are added depends on the justification `method` ∈ (:left, :center, :right).
"""
function pad(text::AbstractString, target_width::Int, method::Symbol)::String
    # get total padding size
    lw = Measure(text).w
    if lw >= target_width
        return text
    end

    npads = target_width - lw

    if method == :left
        return text * " "^npads
    elseif method == :right
        return " "^npads * text
    else
        if npads % 2 == 0
            nl = nr = Int(npads/2)
        else
            nl = (Int ∘ floor)(npads/2)
            nr = nl + 1
        end
        return " "^nl * text * " "^nr
    end
end

"""
    pad(text::AbstractString, left::Int=0, right::Int=0)::String

Pad a string by a fixed ammount.
"""
function pad(text::AbstractString, left::Int=0, right::Int=0)::String
    return " "^left * text * " "^right
end

function pad(segments::Vector{Segment}, left::Int=0, right::Int=0)::Vector{Segment}
    return map(
        (s)->Segment(" "^left * s.text * " "^right),
        segments
    )
end

# ---------------------------------------------------------------------------- #
#                                   STACKING                                   #
# ---------------------------------------------------------------------------- #
"""
    vstack(r1::RenderablesUnion, r2::RenderablesUnion)

Vertically stack two renderables to give a new renderable.
"""
function vstack(r1::RenderablesUnion, r2::RenderablesUnion)
    r1 = Renderable(r1)
    r2 = Renderable(r2)

    # get dimensions of final renderable
    w1 = length(r1.segments) > 1 ? max([s.measure.w for s in r1.segments]...) : r1.measure.w
    w2 = length(r2.segments) > 1 ? max([s.measure.w for s in r2.segments]...) : r2.measure.w
    if w1 > w2
        s1 = r1.segments

        s2 = map(
            (s)->Segment(s.text * " ".^(w1-s.measure.w)),
            r2.segments
        )
    elseif w1 < w2
        s1 = map(
            (s)->Segment(s.text * " ".^(w2-s.measure.w)),
            r1.segments
        )
        s2 = r2.segments
    else
        s1, s2 = r1.segments, r2.segments
    end

    # create segments stack
    segments::Vector{Segment} = vcat(s1, s2)
    measure = Measure(max(w1, w2), length(segments))
    return Renderable(segments, measure)
end

"""
    vstack(renderables...)

Vertically stack a variable number of renderables
"""
function vstack(renderables...)
    renderable = Renderable()

    for ren in renderables
        renderable = vstack(renderable, ren)
    end
    return renderable
end

"""
    hstack(r1::RenderablesUnion, r2::RenderablesUnion)

Horizontally stack two renderables to give a new renderable.
"""
function hstack(r1::RenderablesUnion, r2::RenderablesUnion)
    r1 = r1 isa AbstractRenderable ? r1 : Renderable(r1)
    r2 = r2 isa AbstractRenderable ? r2 : Renderable(r2)

    # get dimensions of final renderable
    h1 = r1.measure.h
    h2 = r2.measure.h
    # h1 = max([s.measure.w for s in r1.segments]...)
    # h2 = max([s.measure.w for s in r2.segments]...)
    Δh = abs(h2 - h1)

    # make sure both renderables have the same number of segments
    if h1 > h2
        s1 = r1.segments
        s2 = vcat(r2.segments, [Segment(" "^(r2.measure.w)) for i in 1:(Δh)])
    elseif h1 < h2
        s1 = vcat(r1.segments, [Segment(" "^(r1.measure.w)) for i in 1:(Δh)])
        s2 = r2.segments
    else
        s1, s2, = r1.segments, r2.segments
    end

    # combine segments
    segments = [Segment(s1.text * s2.text) for (s1, s2) in zip(s1, s2)]

    return Renderable(segments, Measure(r1.measure.w+r2.measure.w, max(r1.measure.h, r2.measure.h)))
end

"""
    hstack(renderables...)

Horizonatlly stack a variable number of renderables
"""
function hstack(renderables...)
    renderable = Renderable()

    for ren in renderables
        renderable = hstack(renderable, ren)
    end
    return renderable
end

# --------------------------------- operators -------------------------------- #

"""
    Operators for more coincise layout of renderables
"""

Base.:/(r1::RenderablesUnion, r2::RenderablesUnion) = vstack(r1, r2)

Base.:*(r1::AbstractRenderable, r2::AbstractRenderable) = hstack(r1, r2)
Base.:*(r1::AbstractString, r2::AbstractRenderable) = hstack(r1, r2)
Base.:*(r1::AbstractRenderable, r2::AbstractString) = hstack(r1, r2)

# ---------------------------------------------------------------------------- #
#                                LINES & SPACER                                #
# ---------------------------------------------------------------------------- #
abstract type AbstractLayoutElement <: AbstractRenderable end

# ---------------------------------- spacer ---------------------------------- #
"""
        Spacer

A box of empty text with given width and height.
"""
mutable struct Spacer <: AbstractLayoutElement
    segments::Vector{Segment}
    measure::Measure
end

function Spacer(width::Number, height::Number; char::Char = ' ')
    width = int(width)
    height = int(height)

    line = char^width
    segments = [Segment(line) for i in 1:height]
    return Spacer(segments, Measure(segments))
end

# ----------------------------------- vline ---------------------------------- #
"""
    vLine

A multi-line renderable with each line made of a | to create a vertical line
"""
mutable struct vLine <: AbstractLayoutElement
    segments::Vector{Segment}
    measure::Measure
    height::Int
end

"""
    vLine(height::Number, style::Union{String, Nothing}; box::Symbol=:ROUNDED)

Create a `vLine` given a height and, optionally, style information.
"""
function vLine(
    height::Number; style::Union{String,Nothing} = nothing, box::Symbol = :ROUNDED
)
    height = int(height)
    char = string(eval(box).head.left)
    segments = [Segment(char, style) for i in 1:height]
    return vLine(segments, Measure(segments), height)
end

"""
    vLine(; style::Union{String, Nothing}=nothing, box::Symbol=:ROUNDED)

Create a `vLine` as tall as the `stdout` console
"""
function vLine(; style::Union{String,Nothing} = nothing, box::Symbol = :ROUNDED)
    return vLine(console_height(); style = style, box = box)
end

"""
    hLine

A 1-line renderable made of repeated character from a Box.
"""
mutable struct hLine <: AbstractLayoutElement
    segments::Vector{Segment}
    measure::Measure
    width::Int
end

"""
    hLine(width::Number, style::Union{String, Nothing}; box::Symbol=:ROUNDED)

Create a styled `hLine` of given width.
"""
function hLine(
    width::Number; style::Union{String,Nothing} = nothing, box::Symbol = :ROUNDED
)
    width = int(width)
    char = eval(box).row.mid
    segments = [Segment(char^width, style)]
    return hLine(segments, Measure(segments), width)
end

"""
    hLine(width::Number, text::String; style::Union{String, Nothing}=nothing, box::Symbol=:ROUNDED)

Creates an hLine object with texte centered horizontally.
"""
function hLine(
    width::Number,
    text::String;
    style::Union{String,Nothing} = nothing,
    box::Symbol = :ROUNDED,
)
    box = eval(box)

    initial_line = box.top.mid^width

    cutval = int(ncodeunits(initial_line) / 2 - ncodeunits(text) - 5)
    cut_start = get_last_valid_str_idx(initial_line, cutval)

    pre = Segment(
        Segment(initial_line[1:cut_start], style) *
        "\e[0m" *
        " " *
        Segment(text, style) *
        " ",
    )

    post = Segment(box.top.mid^(length(initial_line) - pre.measure.w - 1), style)

    segments = [Segment(pre * (post), style)]
    return hLine(segments, Measure(segments), width)
end

"""
    hLine(; style::Union{String, Nothing}=nothing, box::Symbol=:ROUNDED)

Construct an `hLine` as wide as the `stdout`
"""
function hLine(; style::Union{String,Nothing} = nothing, box::Symbol = :ROUNDED)
    return hLine(console_width(); style = style, box = box)
end

"""
    hLine(text::AbstractString; style::Union{String, Nothing}=nothing, box::Symbol=:ROUNDED)

Construct an `hLine` as wide as the `stdout` with centered text.
"""
function hLine(
    text::AbstractString; style::Union{String,Nothing} = nothing, box::Symbol = :ROUNDED
)
    return hLine(console_width(), text; style = style, box = box)
end
end


text:
#export nonlinear

mutable struct nonlinear
    jac;       # The Jacobian matrix
    res;       # The residual vector
    var;     # The solution veltor
	nlvar;
    max_iter;         # max iteration
    atol;
	rtol;
	bilinear;
	linear;
    nonlinear(max_iter, atol, rtol) = new(0,0,0,0,max_iter, atol, rtol,0,0);
end

function init_nonlinear(nl,var,nlvar,bi,li);
	nl.var = var;
	nl.nlvar = nlvar;
	#print("v = \n");
	#@show nlvar[4].values;
	nl.bilinear = bi;
	nl.linear = li;
end

function eval_jac(nl, formjac, t, dt)
	 (nl.jac, b) = formjac(nl.var, nl.bilinear, nl.linear, t, dt);
end

function eval_res(nl, formfunc, t, dt)
	b = formfunc(nl.var, nl.linear, t, dt);
	nl.res = b;#A*nl.var.values - b;
	#@show A;
end

function newton(nl,formjac, formfunc, nlvar, t=0, dt=0)
	debug = true;
	eval_res(nl, formfunc, t, dt);
	init_res = norm(nl.res);
	if (init_res < nl.atol)
		if debug print("\ninitial residual = ", init_res, ", already converged\n"); end
		return;
	end
	if debug print("\ninitial residual = ", init_res, "\n"); end

	i = 0;
	while (i < nl.max_iter)
		eval_jac(nl, formjac, t, dt);
		delta = - nl.jac \ nl.res;

		#@show nl.jac
		#@show nl.res
		#@show(delta)
		#stop
		# place the values in the variable value arrays

		#print("\nlength(nl.var) = ", length(nl.var), "\n");
		#@show nl.var[1].values;
		#@show nl.var[2].values;

		if typeof(nl.var) <: Array
			# for vi=1:length(nl.var)
			# 	components = length(nl.var[vi].symvar);
			# 	for compi=1:components
			# 		nl.var[vi].values[compi,:] = nl.var[vi].values[compi,:]+delta[:];
			# 	end
			# end
			tmp = 0;
			totalcomponents = 0;
			for vi=1:length(nl.var)
				totalcomponents = totalcomponents + length(nl.nlvar[vi].symvar);
				#@show(totalcomponents)
			end
			for vi=1:length(nl.var)
				components = length(nl.var[vi].symvar);
				#@show(components)
				for compi=1:components
					nl.var[vi].values[compi,:] = delta[(compi+tmp):totalcomponents:end];
					nl.nlvar[vi].values[compi,:] = nl.nlvar[vi].values[compi,:]+delta[(compi+tmp):totalcomponents:end];
					tmp = tmp + 1;
				end
			end
		else
			#print("\nHERE\n");
			components = length(nl.var.symvar);
			#print("\ncomponents = ", components, "\n");
			#print("\nsymvar = ", nl.var.symvar, "\n");
			#print("\nvals = ", nl.var.symvar, "\n");
			for compi=1:components
				nl.var.values[compi,:] =  delta[compi:components:end];
				nl.nlvar.values[compi,:] =  nl.nlvar.values[compi,:]+delta[compi:components:end];
			end
		end

		nlvar = nl.nlvar;
		#@show(nlvar.values)
		i = i+1;
		eval_res(nl, formfunc, t, dt);
		curr_res = norm(nl.res);
		if debug print(i,"th iteration residual = ", curr_res, "\n");end
		if (curr_res < nl.atol || curr_res/init_res < nl.rtol)
			if debug print("\nsolution is converged in ", i, " iterations\n");end
			return;
		end
	end
	print("\nsolution is not converged\n");
end


text:
using sem, Feather, ModelingToolkit, Statistics, LinearAlgebra,
    Optim, SparseArrays, Test, Zygote, LineSearches, ForwardDiff

## Observed Data
growth_dat = Feather.read("test/comparisons/growth_dat.feather")
growth_dat_miss30 = Feather.read("test/comparisons/growth_dat_miss30.feather")
#growth_par = Feather.read("test/comparisons/growth_par.feather")

definition_dat = Feather.read("test/comparisons/definition_dat.feather")
definition_dat_unique =
    Feather.read("test/comparisons/definition_dat_unique.feather")
definition_par = Feather.read("test/comparisons/definition_par.feather")
definition_par_missing =
    Feather.read("test/comparisons/definition_par_missing.feather")
definition_par_missing_unique =
    Feather.read("test/comparisons/definition_par_missing_unique.feather")

semobserved = SemObsCommon(data = Matrix(growth_dat); meanstructure = true)

diff_fin = SemFiniteDiff(BFGS(), Optim.Options())

start_val = [1.0, 1, 1, 1, 0.05, 0.05, 0.0, 1.0, 1.0]

###
## Model definition
@ModelingToolkit.variables x[1:7], m[1:2], load_t[1:4]

S =[x[1]  0     0     0     0     0
    0     x[2]  0     0     0     0
    0     0     x[3]  0     0     0
    0     0     0     x[4]  0     0
    0     0     0     0     x[5]  x[7]
    0     0     0     0     x[7]  x[6]]

F =[1.0 0 0 0 0 0
    0 1 0 0 0 0
    0 0 1 0 0 0
    0 0 0 1 0 0]

A =[0  0  0  0  1.0  load_t[1]
    0  0  0  0  1  load_t[2]
    0  0  0  0  1  load_t[3]
    0  0  0  0  1  load_t[4]
    0  0  0  0  0  0
    0  0  0  0  0  0]

M = [0.0, 0, 0, 0, m[1:2]...]

S = sparse(S)

#F
F = sparse(F)

#A
A = sparse(A)

data_def = Matrix(definition_dat)

#abstract type Imply end

struct ImplySymbolicDefinition{
        F <: Any,
        A <: AbstractArray,
        S <: Array{Float64},
        F2 <: Any,
        A2 <: Union{Nothing, AbstractArray},
        I <: Int64,
        P <: Int64,
        R <: AbstractArray,
        D <: AbstractArray} <: Imply
    imp_fun::F
    imp_cov::A # Array of matrices
    start_val::S
    imp_fun_mean::F2
    imp_mean::A2 # Array of matrices of meanstructure
    n_obs::I
    n_patterns::P
    rows::R
    data_def::D
end

function ImplySymbolicDefinition(
    A::Spa1,
    S::Spa2,
    F::Spa3,
    M::Spa4,
    parameters,
    def_vars,
    start_val,
    data_def
        ) where {
        Spa1 <: SparseMatrixCSC,
        Spa2 <: SparseMatrixCSC,
        Spa3 <: SparseMatrixCSC,
        Spa4 <: Union{Nothing, AbstractArray}
        }

    n_obs = size(data_def, 1)
    n_def_vars = Float64(size(data_def, 2))

    patterns = [data_def[i, :] for i = 1:n_obs]
    remember = Vector{Vector{Float64}}()
    rows = [Vector{Int64}(undef, 0) for i = 1:n_obs]

    for i = 1:size(patterns, 1)
        unknown = true
        for j = 1:size(remember, 1)
            if patterns[i] == remember[j]
                push!(rows[j], i)
                unknown = false
            end
        end
        if unknown
            push!(remember, patterns[i])
            push!(rows[size(remember, 1)], i)
        end
    end

    rows = rows[1:length(remember)]
    n_patterns = size(rows, 1)

    pattern_n_obs = size.(rows, 1)

    #############################################
    #Model-implied covmat
    invia = sem.neumann_series(A)

    imp_cov_sym = F*invia*S*permutedims(invia)*permutedims(F)

    imp_cov_sym = Array(imp_cov_sym)
    imp_cov_sym = ModelingToolkit.simplify.(imp_cov_sym)

    imp_fun =
        eval(ModelingToolkit.build_function(
            imp_cov_sym,
            parameters,
            def_vars
        )[2])

    imp_cov = [zeros(size(F)[1], size(F)[1]) for i = 1:n_patterns]

    #Model implied mean
    imp_mean_sym = F*invia*M
    imp_mean_sym = Array(imp_mean_sym)
    imp_mean_sym = ModelingToolkit.simplify.(imp_mean_sym)

    imp_fun_mean =
        eval(ModelingToolkit.build_function(
            imp_mean_sym,
            parameters,
            def_vars
        )[2])

    imp_mean = [zeros(size(F)[1]) for i = 1:n_patterns]

    data_def = remember

    return ImplySymbolicDefinition(
        imp_fun,
        imp_cov,
        copy(start_val),
        imp_fun_mean,
        imp_mean,
        n_obs,
        n_patterns,
        rows,
        data_def
    )
end

function(imply::ImplySymbolicDefinition)(parameters)
    for i = 1:imply.n_patterns
        let (cov,
            mean,
            def_vars) =
                (imply.imp_cov[i],
                imply.imp_mean[i],
                imply.data_def[i])
            imply.imp_fun(cov, parameters, def_vars)
            imply.imp_fun_mean(mean, parameters, def_vars)
        end
    end
end

imply = ImplySymbolicDefinition(
    A,
    S,
    F,
    M,
    [x[1:7]..., m[1:2]...],
    load_t,
    start_val,
    data_def
    )

@benchmark imply($start_val)

struct SemDefinition{ #################### call it per person or sth????
        INV <: AbstractArray,
        C <: AbstractArray,
        L <: AbstractArray,
        M <: AbstractArray,
        I <: AbstractArray,
        IM <: Union{AbstractArray, Nothing},
        T <: AbstractArray,
        DP <: AbstractArray,
        K <: Union{AbstractArray, Nothing},
        R <: Union{AbstractArray, Nothing},
        U,
        V} <: LossFunction
    inverses::INV #preallocated inverses of imp_cov
    choleskys::C #preallocated choleskys
    logdets::L #logdets of implied covmats
    meandiff::M
    imp_inv::I
    imp_mean::IM
    mult::T
    data_perperson::DP
    keys::K
    rows::R
    objective::U
    grad::V
end

function SemDefinition(
        observed::O where {O <: SemObs},
        imply::I where {I <: ImplySymbolicDefinition},
        objective,
        grad)
    n_obs = Int64(observed.n_obs)
    n_man = Int64(observed.n_man)
    n_patterns = imply.n_patterns

    inverses = [zeros(n_man, n_man) for i = 1:n_patterns]
    choleskys = Array{Cholesky{Float64,Array{Float64,2}},1}(undef, n_patterns)
    logdets = zeros(n_patterns)

    meandiff = [zeros(n_man) for i = 1:n_obs]

    imp_inv = zeros(size(observed.data, 2), size(observed.data, 2))
    mult = similar.(inverses)

    data_perperson = [observed.data[i, :] for i = 1:n_obs]

    return SemDefinition(
    inverses,
    choleskys,
    logdets,
    meandiff,
    imp_inv,
    nothing,
    mult,
    data_perperson,
    nothing,
    nothing,
    copy(objective),
    copy(grad)
    )
end

function (semdef::SemDefinition)(par, model::Sem{O, I, L, D}) where
    {O <: SemObs, L <: Loss, I <: Imply, D <: SemFiniteDiff}

    if isnothing(model.imply.imp_mean)
        error("A model implied meanstructure is needed for Definition Variables")
    end

    for i = 1:size(semdef.choleskys, 1)
        semdef.choleskys[i] =
            cholesky!(Hermitian(model.imply.imp_cov[i]); check = false)
    end

    if any(.!isposdef.(semdef.choleskys))
        return Inf
    end

    semdef.logdets .= logdet.(semdef.choleskys)

    for i = 1:size(semdef.inverses, 1)
        semdef.inverses[i] .= LinearAlgebra.inv!(semdef.choleskys[i])
    end

    F = zero(eltype(par))

    for i = 1:Int64(model.imply.n_patterns)
        for j in model.imply.rows[i]
            let (imp_mean, meandiff, inverse, data, logdet) =
                (model.imply.imp_mean[i],
                semdef.meandiff[i],
                semdef.inverses[i],
                semdef.data_perperson[j],
                semdef.logdets[i])

                F += F_one_person(imp_mean, meandiff, inverse, data, logdet)

            end
        end
    end
    return F
end

loss = Loss([SemDefinition(semobserved, imply, 0.0, 0.0)])

model_fin = Sem(semobserved, imply, loss, diff_fin)

model_fin(start_val)

solution = sem_fit(model_fin)
@benchmark sem_fit(model_fin)

par_order = [collect(1:5); 7; 6; 8; 9]

all(
    abs.(solution.minimizer .- definition_par.Estimate[par_order]
        ) .< 0.05*abs.(definition_par.Estimate[par_order]))

## missings ############################################################

data_miss = Matrix(growth_dat_miss30)
semobserved = SemObsMissing(data_miss)

function SemDefinition(
        observed::O where {O <: SemObsMissing},
        imply::I where {I <: ImplySymbolicDefinition},
        objective,
        grad)

    keys = findrow.(1:Int64(observed.n_obs), [observed.rows])

    rows_nested = Vector{Vector{Vector{Int64}}}()
    keys_inner_vec = Vector{Vector{Int64}}()

    for i in 1:length(imply.rows)
        nest = Vector{Vector{Int64}}()
        keys_inner = Vector{Int64}()
        for j in imply.rows[i]
            unknown = true
            key = keys[j]
            for k in 1:size(keys_inner, 1)
                if keys_inner[k] == key
                    push!(nest[k], j)
                    unknown = false
                end
            end
            if unknown
                push!(keys_inner, key)
                push!(nest, [j])
            end
        end
        push!(rows_nested, nest)
        push!(keys_inner_vec, keys_inner)
    end

    data_perperson = [observed.data[i, :] for i = 1:Int64(observed.n_obs)]

    #rows_nested

    inverses = Vector{Vector{Array{Float64, 2}}}()
    for i = 1:size(keys_inner_vec, 1)
        inverses_defgroup = Vector{Array{Float64, 2}}()
        for j in keys_inner_vec[i]
            nvar = Int64(observed.pattern_nvar_obs[j])
            push!(inverses_defgroup, zeros(nvar, nvar))
        end
        push!(inverses, inverses_defgroup)
    end

    #inverses

    #choleskys = Array{Cholesky{Float64,Array{Float64,2}},1}(undef, length(inverses))
    choleskys = [
        Array{Cholesky{Float64,Array{Float64,2}},1}(
            undef, length(inverses[i])) for i in 1:length(inverses)]

    n_patterns = sum([size(inverses[i], 1) for i = 1:length(inverses)])
    logdets = [zeros(size(keys_inner_vec[i], 1)) for i = 1:size(keys_inner_vec, 1)]

    imp_mean = [zeros.(size.(inverses[i], 1)) for i = 1:size(inverses, 1)]
    meandiff = [zeros.(size.(inverses[i], 1)) for i = 1:size(inverses, 1)]

    imp_inv = similar(imply.imp_cov)
    mult = similar.(inverses)

    return SemDefinition(
    inverses,
    choleskys,
    logdets,
    meandiff,
    imp_inv,
    imp_mean,
    mult,
    data_perperson,
    keys_inner_vec,
    rows_nested,
    copy(objective),
    copy(grad)
    )
end

function (semdef::SemDefinition)(par, model::Sem{O, I, L, D}) where
    {O <: SemObsMissing, L <: Loss, I <: Imply, D <: SemFiniteDiff}

    if isnothing(model.imply.imp_mean)
        error("A model implied meanstructure is needed for Definition Variables")
    end

    ################################################
    #copyto!(semfiml.imp_inv, model.imply.imp_cov)
    #a = cholesky!(Hermitian(semfiml.imp_inv); check = false)

    #if !isposdef(a)
    #    F = Inf
    #else
    for i = 1:size(semdef.keys, 1)
        for j = 1:size(semdef.keys[i], 1)
            @views semdef.inverses[i][j] .=
                    model.imply.imp_cov[i][
                        model.observed.patterns[semdef.keys[i][j]],
                        model.observed.patterns[semdef.keys[i][j]]]
            semdef.choleskys[i][j] =
                cholesky!(Hermitian(semdef.inverses[i][j]); check = false)
            if !isposdef(semdef.choleskys[i][j]) return Inf end
        end
    end

    @views  for i = 1:size(semdef.keys, 1)
                for j = 1:size(semdef.keys[i], 1)
                    semdef.imp_mean[i][j] .=
                        model.imply.imp_mean[i][
                            model.observed.patterns[semdef.keys[i][j]]]
                end
            end

    for i = 1:size(semdef.keys, 1)
        for j = 1:size(semdef.keys[i], 1)
            semdef.logdets[i][j] = logdet(semdef.choleskys[i][j])
        end
    end

    #semml.imp_inv .= LinearAlgebra.inv!(a)
    for i = 1:size(semdef.keys, 1)
        for j = 1:size(semdef.keys[i], 1)
            semdef.inverses[i][j] .= LinearAlgebra.inv!(semdef.choleskys[i][j])
        end
    end

    F = zero(eltype(par))

    for i = 1:size(semdef.rows, 1)
        for j = 1:size(semdef.rows[i], 1)
            for k = 1:size(semdef.rows[i][j], 1)
                let (imp_mean, meandiff, inverse, data, logdet) =
                    (semdef.imp_mean[i][j],
                    semdef.meandiff[i][j],
                    semdef.inverses[i][j],
                    model.observed.data_perperson[semdef.rows[i][j][k]],
                    semdef.logdets[i][j])

                    F += F_one_person(imp_mean, meandiff, inverse, data, logdet)
                end
            end
        end
    end
    return F
end

function findrow(r, rows)
    for i in 1:length(rows)
        if r ∈ rows[i]
            return i
        end
    end
end

function remove_all_missing(data)
    keep = Vector{Int64}()
    for i = 1:size(data, 1)
        if any(.!ismissing.(data[i, :]))
            push!(keep, i)
        end
    end
    return data[keep, :], keep
end

data_def = data_def[remove_all_missing(data_miss)[2], :]

imply = ImplySymbolicDefinition(
    A,
    S,
    F,
    M,
    [x[1:7]..., m[1:2]...],
    load_t,
    start_val,
    data_def
    )

loss = Loss([SemDefinition(semobserved, imply, 0.0, 0.0)])

semdef = SemDefinition(semobserved, imply, 0.0, 0.0)

typeof(semdef.keys)

model_fin = Sem(semobserved, imply, loss, diff_fin)

model_fin(start_val)

solution = sem_fit(model_fin)

function F_one_person(imp_mean, meandiff, inverse, data, logdet)
    F = logdet
    @. meandiff = data - imp_mean
    F += meandiff'*inverse*meandiff
    return F
end

par_order = [collect(1:5); 7; 6; 8; 9]

all(
    abs.(solution.minimizer .- definition_par_missing.Estimate[par_order]
        ) .< 0.05*abs.(definition_par_missing.Estimate[par_order]))

############################################
### missings with unique defvar patterns ###
############################################

data_def = Matrix(definition_dat_unique)
data_def = data_def[remove_all_missing(data_miss)[2], :]

imply = ImplySymbolicDefinition(
    A,
    S,
    F,
    M,
    [x[1:7]..., m[1:2]...],
    load_t,
    start_val,
    data_def
    )

loss = Loss([SemDefinition(semobserved, imply, 0.0, 0.0)])

#semdef = SemDefinition(semobserved, imply, 0.0, 0.0)

#typeof(semdef.keys)

model_fin = Sem(semobserved, imply, loss, diff_fin)

model_fin(start_val)

solution = sem_fit(model_fin)

function F_one_person(imp_mean, meandiff, inverse, data, logdet)
    F = logdet
    @. meandiff = data - imp_mean
    F += meandiff'*inverse*meandiff
    return F
end

par_order = [collect(1:5); 7; 6; 8; 9]

all(
    abs.(solution.minimizer .- definition_par_missing.Estimate[par_order]
        ) .< 0.05*abs.(definition_par_missing.Estimate[par_order]))


text:
using BHAPtfem
using Test

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  @test media[:mud] == 11.35
end # cd


text:
"""
    nutrient_load(model::ABM, series<:NutrientSeries)

Updates lake nutrient concentration accorting to `series` type.

See [NutrientSeries](@ref) for more details.
"""
nutrient_load!(m::ABM, series::Constant) = nothing

function nutrient_load!(m::ABM, series::Dynamic)
    if m.lake.p.nutrients >= m.init_nutrients
        nutrients_from_sewage = 0.0
        for municipality in municipalities(m)
            b =
                m.max_sewage_water /
                (municipality.households - municipality.tolerance_level_affectors)
            a = -b * municipality.tolerance_level_affectors
            nutrients_from_sewage += a + b * municipality.affectors
        end
        if m.lake.p.nutrients + nutrients_from_sewage < m.init_nutrients
            # In this case sewage recycling drops the nutrient level below initial
            m.lake.p.nutrients += m.init_nutrients - m.lake.p.nutrients
        else
            m.lake.p.nutrients += nutrients_from_sewage
        end
    end
end

function nutrient_load!(m::ABM, series::Noise)
    if haskey(m.properties, :nutrient_stabilise) && m.nutrient_stabilise == m.year && !isempty(series.process.S₁.data)
        # Time to flip to stable
        # TODO: Generalise this
        # For now, the only time we need this is in a S2-T2 transition
        m.nutrient_series = Noise(GeometricBrownianMotionProcess(0.0, 0.05, m.year, m.lake.p.nutrients), 1.0, 2.5)
        series = m.nutrient_series
    end
    step_noise!(series, 1.0, m)
    m.lake.p.nutrients = series.process.curW
end

function step_noise!(noise::Noise, dt, model)
    N = noise.process
    N.dt = dt # Per year
    current = N.curW
    failure = 0
    while failure <= 1000
        setup_next_step!(N, nothing, nothing)
        if noise.min <= current + N.dW <= noise.max
            accept_step!(N, dt, nothing, nothing, false)
            break
        else
            reject_step!(N, dt, nothing, nothing)
            failure += 1
        end
    end
    if failure == 1000
        error("Failed to step nutrient noise. Attempt to decrease range, increase σ, or increase time for bridges.\n
               Debug. Test: $(haskey(model.properties, :test)), Year: $(model.year)\n
               Current: $(current), dW: $(N.dW), minmax: $(noise.min), $(noise.max)\n
               $(model.nutrient_series.process.W)")
    end
    return nothing
end

function nutrient_load!(m::ABM, series::TransientUp)
    if m.year > series.start_year
        if round(m.lake.p.nutrients + m.nutrient_change; sigdigits = 2) >=
           m.target_nutrients
            m.lake.p.nutrients = m.target_nutrients
            if series.post_target_series isa TransientDown
                m.target_nutrients = m.init_nutrients
            end
            m.nutrient_series = series.post_target_series
        else
            m.lake.p.nutrients += m.nutrient_change
        end
    end
end

function nutrient_load!(m::ABM, series::TransientDown)
    if m.year > series.start_year
        if round(m.lake.p.nutrients - m.nutrient_change; sigdigits = 2) <=
           m.target_nutrients
            m.lake.p.nutrients = m.target_nutrients
            if series.post_target_series isa TransientUp
                m.target_nutrients = m.init_nutrients
            end
            m.nutrient_series = series.post_target_series
        else
            m.lake.p.nutrients -= m.nutrient_change
        end
    end
end


text:
module ARGSITER_TESTS

using Reproduce, Test, FileIO

const ARGS_DICT = Dict(
    "a"=>[1,2,3,4],
    "b"=>["1", "2"],
    "c1+c2"=>[("1", 103.0), ("2", 2034.0)]
)

function setup_tests()
end

function reset_tests()
end

function test_args_iter()

        answer = [(1,  Dict(["c1"=>"1", "c2"=>103.0,  "a"=>1, "b"=>"1"]))
                  (2,  Dict(["c1"=>"1", "c2"=>103.0,  "a"=>2, "b"=>"1"]))
                  (3,  Dict(["c1"=>"1", "c2"=>103.0,  "a"=>3, "b"=>"1"]))
                  (4,  Dict(["c1"=>"1", "c2"=>103.0,  "a"=>4, "b"=>"1"]))
                  (5,  Dict(["c1"=>"2", "c2"=>2034.0, "a"=>1, "b"=>"1"]))
                  (6, Dict(["c1"=>"2", "c2"=>2034.0, "a"=>2, "b"=>"1"]))
                  (7, Dict(["c1"=>"2", "c2"=>2034.0, "a"=>3, "b"=>"1"]))
                  (8, Dict(["c1"=>"2", "c2"=>2034.0, "a"=>4, "b"=>"1"]))
                  (9,  Dict(["c1"=>"1", "c2"=>103.0,  "a"=>1, "b"=>"2"]))
                  (10,  Dict(["c1"=>"1", "c2"=>103.0,  "a"=>2, "b"=>"2"]))
                  (11,  Dict(["c1"=>"1", "c2"=>103.0,  "a"=>3, "b"=>"2"]))
                  (12,  Dict(["c1"=>"1", "c2"=>103.0,  "a"=>4, "b"=>"2"]))
                  (13, Dict(["c1"=>"2", "c2"=>2034.0, "a"=>1, "b"=>"2"]))
                  (14, Dict(["c1"=>"2", "c2"=>2034.0, "a"=>2, "b"=>"2"]))
                  (15, Dict(["c1"=>"2", "c2"=>2034.0, "a"=>3, "b"=>"2"]))
                  (16, Dict(["c1"=>"2", "c2"=>2034.0, "a"=>4, "b"=>"2"]))]

    arg_iter = ArgIterator(ARGS_DICT)
    return all([answer[arg[1]] == arg for arg in arg_iter])
end

function test_args_loop_dict()

    answer = [(1, Dict(["c1"=>1, "seed"=> 2, "s1"=>"hello"])),
              (2, Dict(["c1"=>1, "seed"=> 3, "s1"=>"hello"])),
              (3, Dict(["c1"=>1, "seed"=> 4, "s1"=>"hello"])),
              (4, Dict(["c1"=>1, "seed"=> 5, "s1"=>"hello"])),
              (5, Dict(["c1"=>1, "seed"=> 6, "s1"=>"hello"])),
              (6, Dict(["c1"=>2, "seed"=> 2, "s1"=>"hello"])),
              (7, Dict(["c1"=>2, "seed"=> 3, "s1"=>"hello"])),
              (8, Dict(["c1"=>2, "seed"=> 4, "s1"=>"hello"])),
              (9, Dict(["c1"=>2, "seed"=> 5, "s1"=>"hello"])),
              (10, Dict(["c1"=>2, "seed"=> 6, "s1"=>"hello"]))]

    arg_loop = ArgLooper([Dict(["c1"=>1]), Dict(["c1"=>2])], Dict(["s1"=>"hello"]), 2:6, "seed")
    return all([answer[arg[1]] == arg for arg in arg_loop])
end

macro tests()
    @testset "ArgIterator Tests" begin
        setup_tests()
        @test test_args_iter() == true
        @test test_args_loop_dict() == true
        reset_tests()
    end
end

end