text:
include("../src/Util_multi_RJ.jl")
using ReversibleJumpMultipleParticle

function run_rjsmc(N, lbs_set, ubs_set, obj_funcs, ϵ₀, Δϵ_tol, α, Δϵ_window, output;
                   save_interval=0, target=0.234 , parallel=false, N_sample=100, γscale=1., nuscale=1.)
    mkdir(output)
    lap_start = time()
    lap_end = 0
    N_models = length(obj_funcs)
    stepsize_flags = zeros(N_models)
    model_indecies = rand(collect(1:N_models), N)
    particles = generate_init_particles(N, lbs_set, ubs_set, obj_funcs, model_indecies)

    γs = ones(N_models)*γscale
    nus = ones(N_models)*nuscale
    ϵ_old = ϵ₀
    ϵ_new = 0.0
    cnt = 0
    terminate = false
    ESS=N
    Δϵ_mean = zeros(Δϵ_window)
    logfilename = joinpath(output, "log.csv")

    while terminate==false
        lap_end = time()
        lap = lap_end - lap_start
        lap_start = lap_end
        open(logfilename, "a") do fp
            writecsv(fp, Array([cnt, ϵ_old, ESS, lap])')
        end
        cnt += 1

        ϵ_new = calcPA(particles, α)
        Δϵ = ϵ_old - ϵ_new
        Δϵ_mean[2:end] .= Δϵ_mean[1:end-1]
        Δϵ_mean[1] = Δϵ
        if mean(Δϵ_mean) < Δϵ_tol && cnt > (Δϵ_window+100)
            terminate=true
        end

        # resampling
        candidates = deepcopy(particles)
        while ESS < N
            candidates = vcat(candidates, resample(particles, N-ESS))
            ESS = length(candidates)
        end

        F_old = get_obj_func_vals(candidates)
        candidates = reset_weights(candidates)
        model_indecies = get_model_indecies(candidates)
        accept_flags = zeros(length(candidates))
        if parallel == true
            candidates, accept_flags = MCstep_all_parallel_rj(candidates, ϵ_old)
        else
            candidates, accept_flags = MCstep_all_rj(candidates, ϵ_old)
        end

        model_indecies = get_model_indecies(candidates)
        particles = update_weights(candidates, F_old, ϵ_old, ϵ_new)
        particles = drop_zero_weights(particles)
        particles = normalize_weights(particles)

        acceptance = length(particles) / length(candidates)
        MH_acceptances = zeros(length(obj_funcs))
        MH_acceptances .= NaN
        for id in unique(model_indecies)
            flag_temp = accept_flags[model_indecies .== id]
            MH_acceptances[id] = sum(flag_temp)/length(flag_temp)
        end
        open(joinpath(output, "acceptance.csv"), "a") do fp
            writecsv(fp, Array([MH_acceptances' acceptance]))
        end

        model_indecies = get_model_indecies(particles)
        for id in unique(model_indecies)
            if MH_acceptances[id] < target
                stepsize_flags[id] = 1
            end
        end

        for id in unique(model_indecies)
            particles_subset = get_subset(particles, id)
            if stepsize_flags[id] == 0
                particles_subset, γs[id], nus[id] = update_stepsize_inRKHS_all(particles_subset, γs[id], nus[id], target,
                                                                               MH_acceptances[id], parallel=parallel,
                                                                               N=N_sample, adaptive=false)
                particles[model_indecies .== id] = particles_subset
            else
                particles_subset, γs[id], nus[id] = update_stepsize_inRKHS_all(particles_subset, γs[id], nus[id], target,
                                                                               MH_acceptances[id], parallel=parallel,
                                                                               N=N_sample, adaptive=true)
                particles[model_indecies .== id] = particles_subset
            end
        end

        proportion = get_proportion_of_models(particles, length(obj_funcs))
        open(joinpath(output, "proportion.csv"), "a") do fp
            writecsv(fp, vcat(ϵ_new, proportion)')
        end

        if length(particles) == 0
            println("particles were vanished at ϵ=$ϵ_new")
            break
        end

        if save_interval==0
            nothing
        elseif mod(cnt, save_interval) == 0 && terminate == false
            model_indecies = get_model_indecies(particles)
            for id in unique(model_indecies)
                particles_subset = get_subset(particles, id)
                est_params = (lbs_set[id].*(ubs_set[id]./lbs_set[id]).^(hcat(get_values(particles_subset)...)./10))'
                open(joinpath(output, @sprintf("est_param_%d_%05d.csv", id, cnt)), "w") do fp
                    writecsv(fp, est_params)
                end
            end

        elseif terminate == true
            ESS = length(particles)
            lap_end = time()
            lap = lap_end - lap_start
            open(logfilename, "a") do fp
                writecsv(fp, Array([cnt, ϵ_old, ESS, lap])')
            end

            model_indecies = get_model_indecies(particles)
            for id in unique(model_indecies)
                particles_subset = get_subset(particles, id)
                est_params = (lbs_set[id].*(ubs_set[id]./lbs_set[id]).^(hcat(get_values(particles_subset)...)./10))'
                open(joinpath(output, @sprintf("est_param_%d_%05d.csv", id, cnt)), "w") do fp
                    writecsv(fp, est_params)
                end
            end

            proportion = get_proportion_of_models(particles, length(obj_funcs))
            open(joinpath(output, "proportion.csv"), "a") do fp
                writecsv(fp, vcat(ϵ_new, proportion)')
            end
            break
        end

        ϵ_old = ϵ_new
        ESS = length(particles)
    end

end


text:
function check_model_form(model)

  if typeof(model) <:Blanchard_Kahn_Form
    assessment = check_blanchard_kahn_form(model)
  elseif typeof(model) <:Klein_Form
    assessment = check_klein_form(model)
  elseif typeof(model) <:Binder_Pesaran_Form
    assessment = check_binder_pesaran_form(model)
  elseif typeof(model) <:Sims_Form
    assessment = check_sims_form(model)
  elseif typeof(model) <:Gomme_Klein_Form
    assessment = check_gomme_klein_form(model)
  elseif typeof(model) <:Lombardo_Sutherland_Form
    assessment = check_lombardo_sutherland_form(model)
  end

  return assessment

end

function check_blanchard_kahn_form(model::Blanchard_Kahn_Form)

  assessment = "Pass"

  if model.nx < 1
    assessment = "Fail: Too few predetermined variables"
  elseif model.ny < 1
    assessment = "Fail: Too few non-predetermined variables"
  elseif model.nx+model.ny != size(model.a,1)
    assessment = "Fail: A is not conformable with nx and ny"
  elseif size(model.a,1) != size(model.a,2)
    assessment = "Fail: A is not square"
  elseif size(model.a,1) != size(model.c,1)
    assessment = "Fail: A and C must have the same number of rows"
  elseif size(model.c,2) > model.nx
    assessment = "Fail: There are more innovations than predetermined variables"
  elseif size(model.sigma,1) != size(model.sigma,2)
    assessment = "Fail: SIGMA must be square"
  elseif size(model.c,2) != size(model.sigma,1)
    assessment = "Fail: C and SIGMA are not conformable"
  end

  return assessment

end

function check_klein_form(model::Klein_Form)

  assessment = "Pass"

  if model.nx < 1
    assessment = "Fail: Too few predetermined variables"
  elseif model.ny < 1
    assessment = "Fail: Too few non-predetermined variables"
  elseif model.nx+model.ny != size(model.a,1)
    assessment = "Fail: A is not conformable with nx and ny"
  elseif size(model.a,1) != size(model.a,2)
    assessment = "Fail: A is not square"
  elseif size(model.b,1) != size(model.b,2)
    assessment = "Fail: B is not square"
  elseif size(model.a) != size(model.b)
    assessment = "Fail: A and B are not conformable"
  elseif size(model.a,1) != size(model.c,1)
    assessment = "Fail: A and C must have the same number of rows"
  elseif size(model.c,2) > model.nx
    assessment = "Fail: There are more innovations than predetermined variables"
    elseif size(model.sigma,1) != size(model.sigma,2)
    assessment = "Fail: SIGMA must be square"
  elseif size(model.c,2) != size(model.sigma,1)
    assessment = "Fail: C and SIGMA are not conformable"
  end

  return assessment

end

function check_binder_pesaran_form(model::Binder_Pesaran_Form)

  assessment = "Pass"

  if size(model.a,1) != size(model.a,2)
    assessment = "Fail: A is not square"
  elseif size(model.a) != size(model.a1)
    assessment = "Fail: A and A1 are not conformable"
  elseif size(model.a) != size(model.b)
    assessment = "Fail: A and B are not conformable"
  elseif size(model.a,1) != size(model.c,1)
    assessment = "Fail: A and C must have the same number of rows"
  elseif size(model.sigma,1) != size(model.sigma,2)
    assessment = "Fail: SIGMA must be square"
  elseif size(model.c,2) != size(model.sigma,1)
    assessment = "Fail: C and SIGMA are not conformable"
  end

  return assessment

end

function check_sims_form(model::Sims_Form)

  assessment = "Pass"

  if size(model.gamma0,1) != size(model.gamma0,2)
    assessment = "Fail: GAMMA0 is not square"
  elseif size(model.gamma1,1) != size(model.gamma1,2)
    assessment = "Fail: GAMMA1 is not square"
  elseif size(model.gamma0) != size(model.gamma1)
    assessment = "Fail: GAMMA0 and GAMMA1 are not conformable"
  elseif size(model.gamma0,1) != size(model.c,1)
    assessment = "Fail: GAMMA0 and C must have the same number of rows"
  elseif size(model.c,2) != 1
    assessment = "Fail: C must be a vector"
  elseif size(model.gamma0,1) != size(model.psi,1)
    assessment = "Fail: GAMMA0 and PSI must have the same number of rows"
  elseif size(model.gamma0,1) > size(model.pi,1)
    assessment = "Fail: GAMMA0 and PI must have the same number of rows"
  elseif size(model.sigma,1) != size(model.sigma,2)
    assessment = "Fail: SIGMA must be square"
  elseif size(model.psi,2) != size(model.sigma,1)
    assessment = "Fail: PSI and SIGMA are not conformable"
  end

  return assessment

end

function check_gomme_klein_form(model::Gomme_Klein_Form)

  assessment = "Pass"

  if model.nx < 1
    assessment = "Fail: Too few predetermined variables"
  elseif model.ny < 1
    assessment = "Fail: Too few non-predetermined variables"
  elseif model.nx+model.ny != size(model.derivs,1)
    assessment = "Fail: Need first derivative for all model equations"
  elseif 2*size(model.derivs,1) != size(model.derivs,2)
    assessment = "Fail: The dimensions of the first-derivative matrix are incorrect"
  elseif size(model.derivs,2) != size(model.hessians,2)
    assessment = "Fail: The derivative matrix and the hessian matrix are not conformable"
  elseif size(model.hessians,1) != (model.nx+model.ny)*size(model.hessians,2)
    assessment = "Fail: The dimensions of the hessian matrix are incorrect"
  elseif size(model.eta,1) != model.nx
    assessment = "Fail: Eta has too few rows"
  elseif size(model.eta,2) != size(model.sigma,1)
    assessment = "Fail: Eta has too few columns"
  elseif size(model.sigma,1) != size(model.sigma,2)
    assessment = "Fail: SIGMA must be square"
  end

  return assessment

end

function check_lombardo_sutherland_form(model::Lombardo_Sutherland_Form)

  assessment = "Pass"

  if model.nx < 1
    assessment = "Fail: Too few predetermined variables"
  elseif model.ny < 1
    assessment = "Fail: Too few non-predetermined variables"
  elseif model.nx+model.ny != size(model.derivs,1)
    assessment = "Fail: Need first derivative for all model equations"
  elseif 2*size(model.derivs,1) != size(model.derivs,2)
    assessment = "Fail: The dimensions of the first-derivative matrix are incorrect"
  elseif size(model.derivs,2) != size(model.hessians,2)
    assessment = "Fail: The derivative matrix and the hessian matrix are not conformable"
  elseif size(model.hessians,1) != (model.nx+model.ny)*size(model.hessians,2)
    assessment = "Fail: The dimensions of the hessian matrix are incorrect"
  elseif size(model.eta,1) != model.nx
    assessment = "Fail: Eta has too few rows"
  elseif size(model.eta,2) != size(model.sigma,1)
    assessment = "Fail: Eta has too few columns"
  elseif size(model.sigma,1) != size(model.sigma,2)
    assessment = "Fail: SIGMA must be square"
  end

  return assessment

end

ext: jl


text:
function minvalue(values::Array{T,1})::T where T <: Number
    foldl( (lhs,rhs) -> muladd(10,lhs,rhs), unique(sort(values)) )
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.Variables (ModiaMath/Variables/_module.jl)
#

# Define ModelVariables struct and functions that holds all information about the variables in a simulation model

isArray(v::ModiaMath.AbstractVariable)     = typeof(v.value) <: AbstractArray
isScalar(v::ModiaMath.AbstractVariable)    = !(typeof(v.value) <: AbstractArray)
isReal(v::ModiaMath.AbstractVariable)      = typeof(v) <: ModiaMath.AbstractRealVariable
valueLength(v::ModiaMath.AbstractVariable) = (isArray(v) ? length(v.value) : 1)
isUsedInAnalysis(v::ModiaMath.AbstractVariable, analysis::AnalysisType) = Int(v.analysis) <= Int(analysis)

#------------------------------- Extract variables from a hierarchical component -------------------------------
function get_ModelVariables_aux!(model, modelType, var, dict, analysis::AnalysisType)
    dict[model] = true   # Mark that model is going to be inspected
    #println("... inspect modelType = ", modelType)

    for i = 1:fieldcount(modelType)
        field = getfield(model, fieldname(modelType, i))
        ftype = typeof(field)
        #println("... fieldname = ", fieldname(modelType,i), ", ftype = ", ftype)

        if ftype <: ModiaMath.AbstractVariable
            if isUsedInAnalysis(field, analysis)
                # Use field since defined for the desired analysis
                push!(var, field)
            end
        elseif ftype <: ModiaMath.AbstractComponentWithVariables && !haskey(dict, field)  # Only inspect further, if model/field was not yet inspected
            get_ModelVariables_aux!(field, ftype, var, dict, analysis)
        end
    end
end

"""
    var = get_ModelVariables(model) - Return all variables defined in model component

The function returns a vector var that contains all AbstractVariables defined in (the hierarchical) model
"""
function get_ModelVariables(model::Any, analysis::AnalysisType)
    var = ModiaMath.AbstractVariable[]
    modelType = typeof(model)

    dict = IdDict{Any,Bool}()

    if modelType <: ModiaMath.AbstractComponentWithVariables
        get_ModelVariables_aux!(model, modelType, var, dict, analysis)
    end
    #for v in var
    #   println(v._internal.name)
    #end
    return var
end

#------------------------------- Construct a variable table for printing/debugging -------------------------------

numericType(v, analysis) = (v.numericType == XD_EXP &&  analysis != ModiaMath.DynamicAnalysis &&
                              Int(v.derivative.analysis) > Int(analysis) ) ? XA : v.numericType

const staticArrays    = "StaticArrays."
const lenStaticArrays = length(staticArrays)

function shortenedTypeof(value)
    vtype = string(typeof(value))
    if length(vtype) >= lenStaticArrays  && vtype[1:lenStaticArrays ] == staticArrays
        vtype = vtype[lenStaticArrays + 1:end]
    end
    return Symbol(vtype)
end

vecIndex(v) = (isScalar(v) || v.ivar == 0) ? v.ivar : (v.ivar:v.ivar + length(v.value) - 1)

function get_variableTable(var::Vector{ModiaMath.AbstractVariable})
    v_table = DataFrames.DataFrame(name=Symbol[], ValueType=Symbol[], unit=String[],
                numericType=NumericType[], vec=Symbol[], vecIndex=Any[], resultIndex=Any[],
                der=Symbol[], causality=Causality[], min=Any[], max=Any[], nominal=Any[], fixed=Bool[],
                start=Any[], info=Symbol[])

    for v in var
        if isReal(v)
            unit        = v.unit
            numericType = v.numericType
            vec         = NumericTypeToVector[Int(v.numericType)]
            der         = typeof(v.derivative) == Nothing ? Symbol("") : instanceName(v.derivative)
            vmin        = v.min
            vmax        = v.max
            vnominal    = v.nominal
        else
            unit        = Symbol("")
            numericType = NoNumericType
            vec         = Symbol("")
            der         = Symbol("")
            vmin        = Symbol("")
            vmax        = Symbol("")
            vnominal    = Symbol("")
        end
        resultIndex = (isScalar(v) || v.iresult == 0) ? v.iresult : (v.iresult:v.iresult + length(v.value) - 1)
        push!(v_table, [instanceName(v), shortenedTypeof(v.value), unit, numericType, vec,
                      vecIndex(v), resultIndex, der, v.causality, vmin, vmax, vnominal, v.fixed, v.start, Symbol(v.info)])
    end

    return v_table
end

function indexToString(name, A, linearIndex)
    index = CartesianIndices(A)[linearIndex]
    s     = string(name, "[")
    for i in 1:length(index)
        s = string(s, index[i])

        if i == length(index)
            s = string(s, "]")
        else
            s = string(s, ",")
        end
    end
    return s
end

function add_xName!(v::RealVariable, vnumType, xNames, ix_beg)
    if vnumType == XD_EXP || vnumType == XD_IMP || vnumType == XA
        name = string(instanceName(v))
    elseif vnumType == LAMBDA
        name = "integral(" * string(instanceName(v)) * ")"
    elseif vnumType == MUE
        name = "integral(" * string(instanceName(v)) * ")"
    else
        error("... should not occur")
    end

    if isScalar(v)
        xNames[ix_beg] = string(name)
    else
        for i = 1:length(v.value)
            xNames[ix_beg + i - 1] = indexToString(name, v.value, i)
        end
    end
end

function pushNames!(v::RealVariable, vnumType, xNames, residueNames)
    if vnumType == XD_EXP || vnumType == XD_IMP || vnumType == XA || vnumType == MUE
        push!(xNames, string(instanceName(v)))
    elseif vnumType == LAMBDA
        push!(xNames, "integral(" * string(instanceName(v)) * ")")
    elseif vnumType == FD_IMP || vnumType == FC
        push!(residueNames, string(instanceName(v)))
    end

    if vnumType == XD_EXP
        if typeof(v.derivative) == Nothing
            error(instanceName(v), " is defined as XD_EXP, but no variable is defined as its derivative.")
        end
        push!(residueNames, string(instanceName(v.derivative)) * " - derx[.]")
    end
end

"""
    vars = ModiaMath.ModelVariables(model::ModiaMath.AbstractComponentWithVariables;
                                    analysis::AnalysisType = ModiaMath.DynamicAnalysis)

Return a struct that contains all the variables in `model` in a form so that
fast copying from the integrator interface to the variables can be performed,
and vice versa.
"""
mutable struct ModelVariables
   # Variables used by simulator
    var::Vector{ModiaMath.AbstractVariable}

    # Dimensions of x and derx vector
    nx::Int              # nx = nx_exp + nx_imp + nx_alg + nx_lambda + nx_mue
    nx_exp::Int
    nx_imp::Int
    nx_alg::Int
    nx_lambda::Int
    nx_mue::Int

    # Dimensions of residue vector
    nfd::Int
    nfd_imp::Int    # nfd_exp = nx_exp
    nfc::Int

    # Dimensions of auxiliary variables
    nwr::Int
    nwc::Int

    # Info to copy variables to x/derx/residue/result-values
    x_var::Vector{RealVariable}            # [x_exp_var     , x_imp_var , x_alg_var]
    derx_var::Vector{RealVariable}         # [der(x_imp_var), lambda_var, mue_var]
    residue_var::Vector{RealVariable}      # [fd_imp_var    , fc_var]
    result_var::Vector{ModiaMath.AbstractVariable}   # [time, x_var, derx_var, wr_var, wc_var]
    result_names::Vector{String}           # Names of the (scalarized) result-vector elements
    x_names::Vector{String}                # Names of the (scalarized) x-vector elements

    # Dimensions of the parts of the variable vectors (equivalent to the dimensions of the value vectors like "x" above)
    nx_exp_var::Int
    nx_imp_var::Int
    nx_alg_var::Int
    nx_lambda_var::Int
    nx_mue_var::Int
    nfd_imp_var::Int
    nfc_var::Int
    nwr_var::Int
    nwc_var::Int

    # Dummy
    dummyDifferentialEquation::Bool   # = true, if only dummy differential equation

    function ModelVariables(model::ModiaMath.AbstractComponentWithVariables; analysis::AnalysisType=ModiaMath.DynamicAnalysis)
        # Extract variables from model
        var = get_ModelVariables(model, analysis)
        @assert(length(var) > 0)

        # Determine dimensions of the various parts of the vectors
        ndim = fill(0, 12)
        nvar = fill(0, 12)

        for v in var
            numType = Int(numericType(v, analysis))
            ndim[numType] += valueLength(v)
            nvar[numType] += 1
        end

        (nx_exp, nx_imp, nx_alg, nx_lambda, nx_mue, nderx_exp, nderx_imp, nfd_imp, nfc, nwr, nwc, dummy1) = ndim
        (nx_exp_var, nx_imp_var, nx_alg_var, nx_lambda_var, nx_mue_var, nderx_exp_var, nderx_imp_var, nfd_imp_var, nfc_var, nwr_var, nwc_var, dummy2) = nvar
        @assert(nx_exp >= nderx_exp)
        @assert(nx_imp == nderx_imp)
        #println("nx_exp = ", nx_exp, ", nx_imp = ", nx_imp, ", nx_alg = ", nx_alg, ", nx_lambda = ", nx_lambda, ", nx_mue = ", nx_mue)
        nx = nx_exp + nx_imp + nx_alg + nx_lambda + nx_mue

        if nx == 0
            # add dummy differential equation der(x) = -x; x(0)=0
            dummyDifferentialEquation = true
            dummy_x    = RealScalar("_dummy_x"   ; fixed=true, numericType=XD_EXP)
            dummy_derx = RealScalar("_dummy_derx";             numericType=DER_XD_EXP, integral=dummy_x)

            pushfirst!(var, dummy_derx)
            pushfirst!(var, dummy_x)

            nx        = 1
            nx_exp    = 1
            nderx_exp = 1
            nx_exp_var    = 1
            nderx_exp_var = 1
        else
            dummyDifferentialEquation = false
        end
        time = RealScalar("time"; unit="s", causality=Independent, numericType=TIME)
        pushfirst!(var, time)

        # println("\n... v_table = ", get_variableTable(var))

        # Check dimensions
        if nx_exp + nfd_imp + nfc != nx
            xNames       = String[]
            residueNames = String[]
            for i in eachindex(var)
                v = var[i]
                vnumType = numericType(v, analysis)
                pushNames!(v, vnumType, xNames, residueNames)
            end

            error("The number of x-variables (= ", nx, ") is not identical to the number of equations (= ", nx_exp + nfd_imp + nfc, "):\n",
               "x-variables = ", xNames, "\n",
               "residues = ", residueNames)
        end

        # Allocate variable/index vectors to copy x-, derx-values to variables and variables to residues
        x_var     = Array{ModiaMath.AbstractRealVariable}(undef, nx_exp_var + nx_imp_var + nx_alg_var)   # = [x_exp_var, x_imp_var, x_alg_var]
        ix_exp    = 1
        ix_imp    = nx_exp    + 1
        ix_alg    = ix_imp    + nx_imp
        ix_lambda = ix_alg    + nx_alg
        ix_mue    = ix_lambda + nx_lambda
        ix_exp_var    = 1
        ix_imp_var    = nx_exp_var    + 1
        ix_alg_var    = ix_imp_var    + nx_imp_var
        ix_lambda_var = ix_alg_var    + nx_alg_var
        ix_mue_var    = ix_lambda_var + nx_lambda_var

        derx_var     = Array{ModiaMath.AbstractRealVariable}(undef, nx_imp_var + nx_lambda_var + nx_mue_var)  # = [der(x_imp_var), lambda_var, mue_var]
        iderx_imp    = 1
        iderx_lambda = nx_imp + 1
        iderx_mue    = nx_imp + nx_lambda + 1
        iderx_imp_var    = 1
        iderx_lambda_var = nx_imp_var + 1
        iderx_mue_var    = nx_imp_var + nx_lambda_var + 1

        residue_var  = Array{ModiaMath.AbstractRealVariable}(undef, nfd_imp_var + nfc_var)   # = [fd_imp_var, fc_var]
        ifd_imp = 1
        ifc     = nfd_imp + 1
        ifd_imp_var = 1
        ifc_var     = nfd_imp_var + 1

        result_var   = Array{ModiaMath.AbstractVariable}(undef, 1 + length(x_var) + nderx_exp_var + length(derx_var) + nwr_var + nwc_var - (dummyDifferentialEquation ? 2 : 0))    # [time, x_var, derx_var, wr_var, wc_var]
        result_names = Array{String}(undef, 1 + nx_exp + nx_imp + nx_alg + nderx_exp + nx_imp + nx_lambda + nx_mue + nwr + nwc - (dummyDifferentialEquation ? 2 : 0))
        result_names[1] = "time"
        iresult     = 1
        iresult_var = 1

        x_names      = Array{String}(undef, nx)

        #println("... length(result_var) = ", length(result_var), ", length(result_names) = ", length(result_names), ", nx_exp = ", nx_exp)

        for i in eachindex(var)
            v      = var[i]
            vnumType = numericType(v, analysis)
            if vnumType == FD_IMP
                v.ivar                   = nx_exp + ifd_imp
                residue_var[ifd_imp_var] = v
                ifd_imp                 += valueLength(v)
                ifd_imp_var             += 1
            elseif vnumType == FC
                v.ivar               = nx_exp + ifc
                residue_var[ifc_var] = v
                ifc                 += valueLength(v)
                ifc_var             += 1
            else
                if !(dummyDifferentialEquation && (vnumType == XD_EXP || vnumType == DER_XD_EXP))
                    # Determine information about result vectors
                    var_name = string(instanceName(v))
                    if typeof(v.value) == Float64
                        result_names[iresult] = var_name
                        v.iresult = iresult
                        iresult  += 1
                    else
                        v_value = v.value
                        for i in eachindex(v_value)
                            result_names[iresult + i - 1] = indexToString(var_name, v_value, i)
                        end
                        v.iresult = iresult
                        iresult  += length(v_value)
                    end
                    result_var[iresult_var] = v
                    iresult_var += 1
                end

                if vnumType == XD_EXP
                    add_xName!(v, vnumType, x_names, ix_exp)
                    v.ivar             = ix_exp
                    x_var[ix_exp_var]  = v
                    ix_exp            += valueLength(v)
                    ix_exp_var        += 1
                elseif vnumType == XD_IMP
                    add_xName!(v, vnumType, x_names, ix_imp)
                    v.ivar            = ix_imp
                    x_var[ix_imp_var] = v
                    ix_imp           += valueLength(v)
                    ix_imp_var       += 1
                elseif vnumType == XA
                    add_xName!(v, vnumType, x_names, ix_alg)
                    v.ivar            = ix_alg
                    x_var[ix_alg_var] = v
                    ix_alg           += valueLength(v)
                    ix_alg_var       += 1
                elseif vnumType == LAMBDA
                    add_xName!(v, vnumType, x_names, ix_lambda)
                    v.ivar                     = ix_lambda
                    derx_var[iderx_lambda_var] = v
                    ix_lambda        += valueLength(v)
                    iderx_lambda     += valueLength(v)
                    ix_lambda_var    += 1
                    iderx_lambda_var += 1
                elseif vnumType == MUE
                    add_xName!(v, vnumType, x_names, ix_mue)
                    v.ivar                  = ix_mue
                    derx_var[iderx_mue_var] = v
                    ix_mue        += valueLength(v)
                    iderx_mue     += valueLength(v)
                    ix_mue_var    += 1
                    iderx_mue_var += 1
                end
            end
        end

        for v in var
            if v.numericType == DER_XD_IMP
                v.ivar = (v.integral).ivar
                derx_var[iderx_imp_var] = v
                iderx_imp_var += 1
            elseif v.numericType == DER_XD_EXP
                v.ivar = (v.integral).ivar
            end
        end

        new(var, nx, nx_exp, nx_imp, nx_alg, nx_lambda, nx_mue, nx_exp + nfd_imp, nfd_imp, nfc, nwr, nwc,
          x_var, derx_var, residue_var, result_var, result_names, x_names,
          nx_exp_var, nx_imp_var, nx_alg_var, nx_lambda_var, nx_mue_var, nx_exp_var + nfd_imp_var, nfc_var, nwr_var, nwc_var,
          dummyDifferentialEquation)
    end
end

"""
   table = ModiaMath.get_xTable(vars::ModelVariables)

Function returns a DataFrames tables of all the variables stored in x-vector in `vars`.
"""
function get_xTable(m::ModelVariables)
    x_table = DataFrames.DataFrame(x=Symbol[], name=Symbol[], fixed=Bool[], start=Union{Float64,AbstractArray}[])

    for v in m.x_var
        push!(x_table, [Symbol("x[", vecIndex(v), "]"), instanceName(v), v.fixed, v.start])
    end

    # include integral(lambda) and integral(mue)
    for i in m.nx_imp_var+1:length(m.derx_var)
        v = m.derx_var[i]
        push!(x_table, [Symbol("x[", vecIndex(v), "]"), Symbol("integral(",instanceName(v),")"), false, 0.0])
    end
    return x_table
end

"""
   table = ModiaMath.get_copyToVariableTable(vars::ModelVariables)

Function returns a DataFrames tables of all the variables in `vars`
that are copied from x/derx to the variables.
"""
function get_copyToVariableTable(m::ModelVariables)
    copyToVariable_table = DataFrames.DataFrame(source=Symbol[], target=Symbol[])

    for v in m.x_var
        push!(copyToVariable_table, [Symbol("x[", vecIndex(v), "]"), instanceName(v)])
    end

    for v in m.derx_var
        push!(copyToVariable_table, [Symbol("derx[", vecIndex(v), "]"), instanceName(v)])
    end
    return copyToVariable_table
end

"""
   table = ModiaMath.get_copyToResidueTable(vars::ModelVariables)

Function returns a DataFrames tables of all the variables in `vars`
that are copied from variables to the residue vector.
"""
function get_copyToResidueTable(m::ModelVariables)
    r_table = DataFrames.DataFrame(source=Symbol[], target=Symbol[])
    x_var = m.x_var
    for i = 1:m.nx_exp_var
        v     = x_var[i]
        der_v = v.derivative
        push!(r_table, [Symbol("derx[", vecIndex(v), "] - ", instanceName(der_v)), Symbol("residue[", vecIndex(v), "]") ])
    end

    for v in m.residue_var
        push!(r_table, [instanceName(v), Symbol("residue[", vecIndex(v), "]") ])
    end
    return r_table
end

"""
   table = ModiaMath.get_copyToResultTable(vars::ModelVariables)

Function returns a DataFrames tables of all the variables in `vars`
that are copied from variables to the result.
"""
function get_copyToResultTable(m::ModelVariables)
    r_table = DataFrames.DataFrame(source=Symbol[], target=Symbol[], start=Union{Float64,AbstractArray}[])

    for v in m.result_var
        resultIndex = isScalar(v) ? v.iresult : (v.iresult:v.iresult + length(v.value) - 1)
        push!(r_table, [instanceName(v), Symbol("result[", resultIndex, "]"), v.start])
    end
    return r_table
end

"""
   table = print_ModelVariables(obj)

Print all the variables in `obj` in form of DataFrames tables.
`obj` can be of type ModiaMath.ModelVariables or ModiaMath.SimulationModel.
"""
function print_ModelVariables(m::ModelVariables)
    variabletable = ModiaMath.get_variableTable(m.var)
    x_table = ModiaMath.get_xTable(m)
    copyToVariableTable = ModiaMath.get_copyToVariableTable(m)
    copyToResidueTable = ModiaMath.get_copyToResidueTable(m)
    copyToResultTable = ModiaMath.get_copyToResultTable(m)

    print("\n\nvariables: ");                show(variabletable      , splitcols=true, summary=false)
    print("\n\n\nx vector: ");               show(x_table            , splitcols=true, summary=false)
    print("\n\n\ncopy to variables: ");      show(copyToVariableTable, splitcols=true, summary=false)
    print("\n\n\ncopy to residue vector: "); show(copyToResidueTable , splitcols=true, summary=false)
    print("\n\n\ncopy to results: ");        show(copyToResultTable  , splitcols=true, summary=false)
    print("\n")
end

"""
    copy_x_and_derx_to_variables!(time::Float64, x::Vector{Float64},
                                  derx::Vector{Float64}, vars::ModelVariables)

Copy `x`and `derx`of the integrator interface to the model variables `vars`.
"""
function copy_x_and_derx_to_variables!(time::Float64, x::Vector{Float64}, derx::Vector{Float64}, m::ModelVariables)
    @assert(length(x)    == m.nx)
    @assert(length(derx) == m.nx)

    m.var[1].value = time
    for v in m.x_var
        #println("... typeof(", fullName(v), ") = ", typeof(v), ", isimmutable(v) = ", isimmutable(v))
        if isScalar(v)
            v.value = x[ v.ivar ]
        elseif isimmutable(v.value)
            # v is an immutable array (e.g. SVector{3,Float64})
            v.value = x[ v.ivar:v.ivar + length(v.value) - 1 ]
        else
            vv = v.value
            for j in 1:length(vv)
                vv[j] = x[ v.ivar + j - 1 ]
            end
        end
    end

    for v in m.derx_var
        if isScalar(v)
            v.value = derx[ v.ivar ]
        elseif isimmutable(v.value)
            # v is an immutable array (e.g. SVector{3,Float64})
            v.value = derx[ v.ivar:v.ivar + length(v.value) - 1 ]
        else
            vv = v.value
            for j in 1:length(vv)
                vv[j] = derx[ v.ivar + j - 1 ]
            end
        end
    end

    return nothing
end

"""
    copy_start_to_x!(vars::ModelVariables, x::Vector{Float64}, x_fixed::Vector{Bool}, x_nominal::Vector{Float64})

Copy `start`, `fixed` and `nominal`values of variables `vars` to `x`, `x_fixed`, and `x_nominal` vectors.
"""
function copy_start_to_x!(m::ModelVariables, x::Vector{Float64}, x_fixed::Vector{Bool}, x_nominal::Vector{Float64})
    @assert(length(x)         == m.nx)
    @assert(length(x_fixed)   == m.nx)
    @assert(length(x_nominal) == m.nx)

    for v in m.x_var
        ibeg = v.ivar
        if isScalar(v)
            x[ibeg]         = v.start
            x_fixed[ibeg]   = v.fixed
            x_nominal[ibeg] = v.nominal
        else
            vv = v.start
            for j in 1:length(vv)
                k = ibeg + j - 1
                x[ k ]       = vv[j]
                x_fixed[   k ] = v.fixed
                x_nominal[ k ] = v.nominal
            end
        end
    end
    return nothing
end
copy_start_to_x!(m::ModelVariables, x::Vector{Float64}, x_fixed::Vector{Bool}) = copy_start_to_x!(m, x, x_fixed, fill(1.0,m.nx))

"""
    copy_variables_to_residue!(vars::ModelVariables, x::Vector{Float64},
                               derx::Vector{Float64}, residue::Vector{Float64})

Copy variables `vars` to `residue` vector and include the inherent residue equations
of `XD_EXP` variables (`residue = der(v) - derx`).
"""
function copy_variables_to_residue!(m::ModelVariables, x::Vector{Float64}, derx::Vector{Float64}, residue::Vector{Float64})
    @assert(length(x)       == m.nx)
    @assert(length(derx)    == m.nx)
    @assert(length(residue) == m.nx)

    # If dummy equations present, compute derivative of dummy equation
    if m.dummyDifferentialEquation
        # der(x) = -x
        m.var[3].value = -m.var[2].value
    end

    # Generate residues of explicitly solvable derivatives
    x_var = m.x_var
    for i = 1:m.nx_exp_var
        v = x_var[i]
        if isScalar(v)
            residue[v.ivar] = derx[v.ivar] - v.derivative.value
        else
            der_vv = v.derivative.value
            for j in 1:length(der_vv)
                k = v.ivar + j - 1
                residue[ k ] = derx[ k ] - der_vv[j]
            end
        end
    end

    # Copy residue variables of fd_imp and fc to residue vector
    for v in m.residue_var
        if isScalar(v)
            residue[ v.ivar ] = v.value
        else
            vv = v.value
            for j in 1:length(vv)
                residue[ v.ivar + j - 1 ] = vv[j]
            end
        end
    end

    return nothing
end

function get_variableValueTable(m::ModelVariables)
    v_table = DataFrames.DataFrame(simulator=Symbol[], variable=Symbol[], value=Any[])

    for v in m.x_var
        push!(v_table, [Symbol("x[", vecIndex(v), "]"), instanceName(v), v.value])
    end

    for v in m.derx_var
        push!(v_table, [Symbol("derx[", vecIndex(v), "]"), instanceName(v), v.value])
    end
    return v_table
end

function get_residueValueTable(m::ModelVariables, derx_integrator::Vector{Float64})
    r_table = DataFrames.DataFrame(simulator=Symbol[], variable=Symbol[], value=Any[])
    x_var = m.x_var
    for i = 1:m.nx_exp_var
        v     = x_var[i]
        der_v = v.derivative
        value = derx_integrator[ vecIndex(v) ] - der_v.value
        push!(r_table, [Symbol("residue[", vecIndex(v), "]"), Symbol("derx[", vecIndex(v), "] - ", instanceName(der_v)), value ])
    end

    for v in m.residue_var
        push!(r_table, [Symbol("residue[", vecIndex(v), "]"), instanceName(v), v.value ])
    end
    return r_table
end


text:
export Affine

struct Affine{T <: Real}
    A::Matrix{T}
    b::Vector{T}
end

using LinearAlgebra
function project!(s::Affine{T}, v::Vector{T}, y::Vector{T}) where {T <: Real}
    v = y .- LinearAlgebra.pinv(s.A) * (s.A * y .- s.b)
end

function project(s::Affine{T}, y::Vector{T}) where {T <: Real}
    project!(s, similar(y), y)
end


text:
cd(@__DIR__)
using Pkg; Pkg.activate("."); Pkg.instantiate()

using Statistics, Printf
using Plots

using Oceananigans
using Oceananigans.Diagnostics
using Oceananigans.OutputWriters

# Model parameters
# N = 32  # Number of grid points in each dimension.
Nx = Nz = 128
Ny = 2
L = 1   # Physical length of each dimension.

κ = 1e-4  # Diffusivity
ν = κ     # Assuming Prandtl number Pr = 1

end_time = 10

# Set up model
model = Model(
           architecture = CPU(),
             float_type = Float64,
                   grid = RegularCartesianGrid(size=(Nx, Ny, Nz), x=(-L/2, L/2), y=(-L/2, L/2), z=(-L/2, L/2)),
                tracers = (:b,),
               coriolis = nothing,
               buoyancy = BuoyancyTracer(),
                closure = ConstantIsotropicDiffusivity(ν=ν, κ=κ),
    boundary_conditions = HorizontallyPeriodicSolutionBCs()
)

# Setting initial conditions
ε(σ) = σ * randn()  # Noise

# b₀(x, y, z) = 1 * (z < 0) + -1 * (z > 0) + ε(1e-2)

@inline interface(x, y, z) = 0.05sin(6π*x)
@inline function b₀(x, y, z)
    z > interface(x, y, z) && return -1
    return 1
end

set!(model, b=b₀)

# Add JLD2 output writer for 3D fields.
fields = Dict(
    :u => model -> Array(interior(model.velocities.u)),
    :v => model -> Array(interior(model.velocities.v)),
    :w => model -> Array(interior(model.velocities.w)),
    :b => model -> Array(interior(model.tracers.b))
)

model.output_writers[:fields] =
    JLD2OutputWriter(model, fields, dir=".", prefix="rayleigh_taylor_instability_3d_fields",
                     interval=1, force=true, verbose=true)

# Add JLD2 output writer for horizontal averages.
 u̅ = HorizontalAverage(model.velocities.u; return_type=Array)
 v̅ = HorizontalAverage(model.velocities.v; return_type=Array)
 w̅ = HorizontalAverage(model.velocities.w; return_type=Array)
 b̅ = HorizontalAverage(model.tracers.b; return_type=Array)

 profiles = Dict(
    :u => model -> u̅(model)[2:end-1],
    :v => model -> v̅(model)[2:end-1],
    :w => model -> w̅(model)[2:end-1],
    :b => model -> b̅(model)[2:end-1]
)

model.output_writers[:horizontal_averages] =
    JLD2OutputWriter(model, profiles, dir=".", prefix="rayleigh_taylor_instability_3d_horizontal_averages",
                     interval=0.1, force=true, verbose=true)

# Set up adaptive time stepping.
wizard = TimeStepWizard(cfl=0.2, Δt=1e-6, max_change=1.2, max_Δt=5e-3)

# Set up CFL diagnostics.
cfl = AdvectiveCFL(wizard)

# Number of time steps to perform at a time before printing a progress
# statement and updating the adaptive time step.
Ni = 20

function coarse_grain(data, resolution)
    @assert length(data) % resolution == 0
    s = length(data) / resolution

    data_cs = zeros(resolution)
    for i in 1:resolution
        t = data[Int((i-1)*s+1):Int(i*s)]
        data_cs[i] = mean(t)
    end

    return data_cs
end

function plot_buoyancy(model)
    t_str = @sprintf("t = %.2f", model.clock.time)
    xC, zC = model.grid.xC, model.grid.zC

    b_horizontal_average = b̅(model)[2:end-1]
    b_profile = plot(b_horizontal_average, zC, title=t_str, label="simulation", xlabel="buoyancy", ylabel="z",
                     xlims=(-1, 1), ylims=(-L/2, L/2))

    coarse_resolution = cr = 16
    b_cr = coarse_grain(b_horizontal_average, cr)
    zC_cr = coarse_grain(collect(zC), cr)

    plot!(b_cr, zC_cr, label="coarsened")

    j½ = Int(model.grid.Ny/2)
    b_xz_data = interior(model.tracers.b)[:, j½, :]'
    b_slice = heatmap(xC, zC, b_xz_data, title="buoyancy slice", color=:balance,
                      xlims=(-L/2, L/2), ylims=(-L/2, L/2), clims=(-1, 1), aspect_ratio=:equal)

    display(plot(b_profile, b_slice, layout=(1, 2), show=true))
end

@time while model.clock.time < end_time
    walltime = @elapsed time_step!(model; Nt=Ni, Δt=wizard.Δt)

    # Calculate simulation progress in %.
    progress = 100 * (model.clock.time / end_time)

    # Calculate maximum velocities.
    umax = maximum(abs, model.velocities.u.data.parent)
    vmax = maximum(abs, model.velocities.v.data.parent)
    wmax = maximum(abs, model.velocities.w.data.parent)

    # Calculate a new adaptive time step.
    update_Δt!(wizard, model)

    # Print progress statement.
    i, t = model.clock.iteration, model.clock.time
    @printf("[%05.2f%%] i: %d, t: %.2e, umax: (%.3e, %.3e, %.3e), CFL: %.4e, next Δt: %.2e, ⟨wall time⟩: %s\n",
            progress, i, t, umax, vmax, wmax, cfl(model), wizard.Δt, prettytime(walltime / Ni))

    # plot_buoyancy(model)
end

# mp4(anim, "rayleigh_taylor_instability_3d.mp4", fps=15)

@time begin
    # Set up model
    model = Model(
               architecture = CPU(),
                 float_type = Float64,
                       grid = RegularCartesianGrid(size=(Nx, Ny, Nz), x=(-L/2, L/2), y=(-L/2, L/2), z=(-L/2, L/2)),
                    tracers = (:b,),
                   coriolis = nothing,
                   buoyancy = BuoyancyTracer(),
                    closure = ConstantIsotropicDiffusivity(ν=ν, κ=κ),
        boundary_conditions = HorizontallyPeriodicSolutionBCs()
    )
    set!(model, b=b₀)
    # Set up adaptive time stepping.
    wizard = TimeStepWizard(cfl=0.2, Δt=1e-6, max_change=1.2, max_Δt=5e-3)

    # Set up CFL diagnostics.
    cfl = AdvectiveCFL(wizard)

    while model.clock.time < end_time
        walltime = @elapsed time_step!(model; Nt=Ni, Δt=wizard.Δt)

        # Calculate simulation progress in %.
        progress = 100 * (model.clock.time / end_time)

        # Calculate maximum velocities.
        umax = maximum(abs, model.velocities.u.data.parent)
        vmax = maximum(abs, model.velocities.v.data.parent)
        wmax = maximum(abs, model.velocities.w.data.parent)

        # Calculate a new adaptive time step.
        update_Δt!(wizard, model)

        # Print progress statement.
        i, t = model.clock.iteration, model.clock.time
        @printf("[%05.2f%%] i: %d, t: %.2e, umax: (%.3e, %.3e, %.3e), CFL: %.4e, next Δt: %.2e, ⟨wall time⟩: %s\n",
                progress, i, t, umax, vmax, wmax, cfl(model), wizard.Δt, prettytime(walltime / Ni))

        # plot_buoyancy(model)
    end
end


text:
# # Schema Examination
# In this example we build schema of documents with complex structure and show how can we filter it and perform transformations.
# We start by adding libraries we want to use

#md # !!! tip
#md #     This example is also available as a Jupyter notebook, feel free to run it yourself:
#md #     [`schema_examination.ipynb`](@__NBVIEWER_ROOT_URL__/examples/schema_examination.ipynb)

#nb # We start by installing JsonGrinder and few other packages we need for the example.
#nb # Julia Ecosystem follows philosophy of many small single-purpose composable packages
#nb # which may be different from e.g. python where we usually use fewer larger packages.
#nb using Pkg
#nb pkg"add JsonGrinder#master Flux Mill#master MLDataPattern JSON HierarchicalUtils StatsBase OrderedCollections"

using JsonGrinder, Flux, Mill, MLDataPattern, JSON, HierarchicalUtils, StatsBase, OrderedCollections
using JsonGrinder: DictEntry, Entry

data_dir = "data/documents" #src
data_dir = "../../../data/documents" #nb
data_dir = "data/documents" #md
data_dir = "data/documents" #jl

# This is how some of the documents look like:
open(JSON.parse, first(readdir(data_dir, join=true)))

# We load files in data/documents and parse them
sch = JsonGrinder.schema(readdir(data_dir, join=true), x->open(JSON.parse, x))
# The default printing method restricts depth and width of the printed schema.
# We can see the whole schema using the `printtree` function from [HierarchicalUtils](https://github.com/CTUAvastLab/HierarchicalUtils.jl).
# The htrunc and vtrunc kwargs tell us maximum number of keys and max depth that will be rendered, respectively.
printtree(sch, htrunc=20, vtrunc=20)

# We suggest default extractor.
extractor = suggestextractor(sch)

# We show the almost whole extractor. Feel free to remove the htrunc and vtrunc kwargs if you want to
# see it whole.
printtree(extractor, htrunc=20, vtrunc=20)

# We see that there are some dictionaries with lots of keys, so let's examine the schema more.
#
# Mill.jl [treats Dictionaries as a cartesian product of their embeddings](https://ctuavastlab.github.io/Mill.jl/stable/manual/nodes/#[ProductNode](@ref)s-and-[ProductModel](@ref)s)
# which does make sense in case when there is consistent number of keys, and keys themselves don't carry semantic meaning.
# Looking at the schema, we can hypothesize many different keys, which occur very scarcely in data, carry semantic information.
#
# We want to examine how many unique keys are there in the schema in order to handle them differently and train also on key names in such case.
# So let's take a look at histogram of number of children per Dictionary.
#
# Function [list_lens](https://ctuavastlab.github.io/Mill.jl/stable/api/utilities/#Mill.list_lens) ¨
# from [Mill.jl](https://github.com/CTUAvastLab/Mill.jl) lets us iterate over all nodes in our tree structure
# in a way we know their position in the schema.
StatsBase.countmap([length(get(sch, i).childs) for i in list_lens(sch) if get(sch, i) isa DictEntry]) |> sort

# We see that 1 dict has 103 unique children, 1 dict has 13 unique children,
# 91 dicts have 9 unique children, 59 dicts don't have any children etc.
#
# We can take a more detailed look at Dicts with > 5 children.
#
# The following code prints paths to all Dictionaries in the schema and number of their children if they have more than 5 children.
# In total there is lots of diction
for i in list_lens(sch)
    e = get(sch, i)
    if e isa DictEntry && length(e.childs) > 5
        @info i length(e.childs)
    end
end

# The dictionaries with most unique children are following ones:
# ```
# ┌ Info: (@lens _.childs[:ref_entries])
# └   length(e.childs) = 13
# ┌ Info: (@lens _.childs[:bib_entries])
# └   length(e.childs) = 103
# ```
# because this is where keys have semantic meaning.
# JsonGrinder contains [ExtractKeyAsField](@ref) extractor, which treats
# dictionaries with large number of keys as array of pairs (key, value)
# which leads to more reasonable model.
#
# There is a default value, but we want to set it ourselves to 13 to cover
# both cases we see in out data. This can be performed by creating new extractor
# like this
extractor = suggestextractor(sch, (; key_as_field=13))

# When we look at the larger part of extractor
printtree(extractor, htrunc=20, vtrunc=20)
# we now see represenation of `bib_entries` and `ref_entries` is
# more reasonable now.
#
# So we can say this extractor looks much better.
#
# But still, some values are very sparse,
# let's print all parts of schema where each value is observed only once
for i in list_lens(sch)
    e = get(sch, i)
    if e isa Entry && maximum(values(e.counts)) == 1
        @info i
    end
end

#  we can see lots of leaves under `bib_entries`, which is cased by uniqueness of keys here
# but apart from that, we can see other interesting fields
# ```
# [ Info: (@lens _.childs[:metadata].childs[:authors].items.childs[:middle].items)
# [ Info: (@lens _.childs[:metadata].childs[:authors].items.childs[:last])
# [ Info: (@lens _.childs[:metadata].childs[:authors].items.childs[:affiliation].childs[:location].childs[:region])
# [ Info: (@lens _.childs[:paper_id])
# [ Info: (@lens _.childs[:body_text].items.childs[:text])
# [ Info: (@lens _.childs[:body_text].items.childs[:ref_spans].items.childs[:start])
# [ Info: (@lens _.childs[:body_text].items.childs[:ref_spans].items.childs[:end])
# [ Info: (@lens _.childs[:back_matter].items.childs[:text])
# [ Info: (@lens _.childs[:back_matter].items.childs[:cite_spans].items.childs[:ref_id])
# [ Info: (@lens _.childs[:back_matter].items.childs[:cite_spans].items.childs[:start])
# [ Info: (@lens _.childs[:back_matter].items.childs[:cite_spans].items.childs[:text])
# [ Info: (@lens _.childs[:back_matter].items.childs[:cite_spans].items.childs[:end])
# [ Info: (@lens _.childs[:back_matter].items.childs[:ref_spans].items.childs[:start])
# [ Info: (@lens _.childs[:back_matter].items.childs[:ref_spans].items.childs[:text])
# [ Info: (@lens _.childs[:back_matter].items.childs[:ref_spans].items.childs[:end])
# ```
#
# Let's remove some of them from the extractor so we don't train on them.
delete!(extractor.dict, :paper_id)
delete!(extractor.dict[:metadata].dict[:authors].item.dict, :last)
delete!(extractor.dict[:metadata].dict[:authors].item.dict, :middle)

# Now the extractor looks even better!
printtree(extractor, htrunc=20, vtrunc=20)

# This concludes example about examining schema and modifying extractor accordingly.


text:
using JultraDark
using Test
using MPI

resol = 16

output_dir = "output"
output_times = 0.0:0.01:0.1

output_config = OutputConfig(output_dir, output_times; box=false)
options = Config.SimulationConfig(10, t->1)

# Grids

grids = Grids(1.0, 16)
grids.ψx .= grids.dist ./ 1e9 # Set ψx to something non-zero
@test simulate(grids, options, output_config) == nothing

# PencilGrids

if ~MPI.Initialized()
    MPI.Init()
end

comm = MPI.COMM_WORLD
# Disable output on all but one process.
const DEV_NULL = @static Sys.iswindows() ? "nul" : "/dev/null"
rank = MPI.Comm_rank(comm)
rank == 0 || redirect_stdout(open(DEV_NULL, "w"))

grids = PencilGrids(1.0, 16)
grids.ψx .= grids.dist ./ 1e9 # Set ψx to something non-zero
@test simulate(grids, options, output_config) == nothing


text:
# # Rossby Haurwitz solutions
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add Oceananigans, JLD2, Plots"
# ```

# ## A spherical domain
#
# We use a one-dimensional domain of geophysical proportions,

using Oceananigans
using Oceananigans.Utils: prettytime
using Oceananigans.Units
using Oceananigans.Operators
using Printf
using Oceananigans.Diagnostics: accurate_cell_advection_timescale
using Oceananigans.Advection: VelocityStencil, VorticityStencil, EnstrophyConservingScheme
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom
# using GLMakie

#  λ for latitude and ϕ for latitude is
using Oceananigans.Advection: EnergyConservingScheme, EnstrophyConservingScheme
using Oceananigans.Coriolis: HydrostaticSphericalCoriolis

# include("visualization.jl")

# ## Building a `HydrostaticFreeSurfaceModel`
#
# We use `grid` and `coriolis` to build a simple `HydrostaticFreeSurfaceModel`,

function run_rossby_haurwitz(; architecture = CPU(),
                               Nx = 90,
                               Ny = 30,
                               coriolis_scheme = EnstrophyConservingScheme(),
                               advection_scheme = VectorInvariant(),
                               prefix = "vector_invariant")

    h₀ = 8e3

    coriolis = HydrostaticSphericalCoriolis(scheme=coriolis_scheme)

    grid = LatitudeLongitudeGrid(architecture, size = (Nx, Ny, 5),
                                 longitude = (-180, 180),
                                 latitude = (-80, 80),
                                 z = (-h₀, 0),
                                 halo = (4, 4, 4),
                                 precompute_metrics = true)

    grid = ImmersedBoundaryGrid(grid, GridFittedBottom((x, y) -> -h₀-1))

    free_surface = ImplicitFreeSurface(solver_method=:HeptadiagonalIterativeSolver, gravitational_acceleration=900)

    model = HydrostaticFreeSurfaceModel(; grid, free_surface,
                                        tracers = (),
                                        momentum_advection = advection_scheme,
                                        buoyancy = nothing,
                                        coriolis = coriolis,
                                        closure = nothing)

    R = model.grid.radius   # [m]
    ω = 0.0 # 7.848e-6            # [s⁻¹]
    K = 7.848e-6            # [s⁻¹]
    n = 4                   # dimensionless
    g = model.free_surface.gravitational_acceleration          # [m/s²]
    Ω = 7.292115e-5
    ϵ = 0.0 # perturbation veloctiy # [m/s]

    A(θ) = ω/2 * (2 * Ω + ω) * cos(θ)^2 + 1/4 * K^2 * cos(θ)^(2*n) * ((n+1) * cos(θ)^2 + (2 * n^2 - n - 2) - 2 * n^2 * sec(θ)^2 )
    B(θ) = 2 * K * (Ω + ω) * ((n+1) * (n+2))^(-1) * cos(θ)^(n) * ( n^2 + 2*n + 2 - (n+1)^2 * cos(θ)^2) # why not  (n+1)^2 sin(θ)^2 + 1
    C(θ)  = 1/4 * K^2 * cos(θ)^(2 * n) * ( (n+1) * cos(θ)^2 - (n+2))

    # here: θ ∈ [-π/2, π/2] is latitude, ϕ ∈ [0, 2π) is longitude
    u₁(θ, ϕ) =  R * ω * cos(θ) + R * K * cos(θ)^(n-1) * (n * sin(θ)^2 - cos(θ)^2) * cos(n*ϕ)
    v₁(θ, ϕ) = -n * K * R * cos(θ)^(n-1) * sin(θ) * sin(n*ϕ)
    h₁(θ, ϕ) =  0*h₀ + R^2/g * (  A(θ)  +  B(θ)  * cos(n * ϕ) + C(θ) * cos(2 * n * ϕ) )

    # Total initial conditions
    # previously: θ ∈ [-π/2, π/2] is latitude, ϕ ∈ [0, 2π) is longitude
    # oceanoganigans: ϕ ∈ [-90, 90], λ ∈ [-180, 180],
    rescale¹(λ) = (λ + 180)/ 360 * 2π # λ to θ
    rescale²(ϕ) = ϕ / 180 * π # θ to ϕ
    # arguments were u(θ, ϕ), λ |-> ϕ, θ |-> ϕ
    uᵢ(λ, ϕ, z) = u₁(rescale²(ϕ), rescale¹(λ))
    vᵢ(λ, ϕ, z) = v₁(rescale²(ϕ), rescale¹(λ))
    hᵢ(λ, ϕ)    = h₁(rescale²(ϕ), rescale¹(λ)) # (rescale¹(λ), rescale²(ϕ))

    u, v, w = model.velocities
    η = model.free_surface.η

    set!(u, uᵢ)
    set!(v, vᵢ)
    set!(η, hᵢ)

    # Time step restricted on the gravity wave speed. If using the implicit free surface method it is possible to increase it
    Δt =0.1*accurate_cell_advection_timescale(model)

    simulation = Simulation(model, Δt = Δt, stop_time = 50days)

    # wizard = TimeStepWizard(cfl=0.5, max_change=3.0, max_Δt=3minutes)

    progress(sim) = @printf("Iter: %d, time: %s, Δt: %s, max|u|: %.3f, max|η|: %.3f \n",
                            iteration(sim), prettytime(sim), prettytime(sim.Δt), maximum(abs, model.velocities.u), maximum(abs, model.free_surface.η))

    simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))
    # simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(100))

    u, v, w = model.velocities
    η=model.free_surface.η

    ζ_op = KernelFunctionOperation{Face, Face, Center}(ζ₃ᶠᶠᶜ, grid; computed_dependencies=(u, v))

    ζ = Field(ζ_op)

    output_fields = (; u = u, v = v, η = η, ζ = ζ)

    simulation.output_writers[:fields] = JLD2OutputWriter(model, output_fields,
                                                        schedule = TimeInterval(40Δt),
                                                        prefix = "rh_$(prefix)_Nx$(Nx)",
                                                        overwrite_existing = true)
    run!(simulation)

    return simulation.output_writers[:fields].filepath
end

filepath_w = run_rossby_haurwitz(architecture=GPU(), Nx=512, Ny=256, advection_scheme=WENO5(vector_invariant=VelocityStencil()), prefix = "WENOVectorInvariantVel")
filepath_w = run_rossby_haurwitz(architecture=GPU(), Nx=512, Ny=256, advection_scheme=WENO5(vector_invariant=VorticityStencil()), prefix = "WENOVectorInvariantVort")
filepath_w = run_rossby_haurwitz(architecture=GPU(), Nx=512, Ny=256, advection_scheme=VectorInvariant(), prefix = "VectorInvariant")


text:
ENV["PYTHON_JL_RUNTIME_PYTHON"] = Sys.which("python")
using SpikeNetOpt
import JLD

const ngt_spikes
const opt_vec
const extremum
const extremum_param
const ngt_spikes
const fitness

using Plots
###
const cell_type = "ADEXP"
global vecp = false
###
using SpikeNetOpt
SNO = SpikeNetOpt

(vmgtv, vmgtt, ngt_spikes, ground_spikes) = SNO.get_data()
println("Ground Truth")
plot(vmgtt[:], vmgtv[:]) |> display

lower, upper = get_adexp_ranges()
ALLEN_DURATION = 2000 * ms
ALLEN_DELAY = 1000 * ms

function loss(E, ngt_spikes, ground_spikes)
    spikes = raster_synchp(E)
    spikes = [s * ms for s in spikes]
    maxt = findmax(sort!(unique(vcat(spikes, ground_spikes))))[1]
    if size(spikes)[1] > 1
        t, S = SPIKE_distance_profile(spikes, ground_spikes; t0 = 0, tf = maxt)
        spkd = SpikeSynchrony.trapezoid_integral(t, S) / (t[end] - t[1]) # == SPIKE_distance(y1, y2)
    else
        spkd = 1.0
    end
    delta = abs(size(spikes)[1] - ngt_spikes)
    return spkd + delta

end

function zz_(param)
    if cell_type == "IZHI"
        pp = SNN.IZParameter(; a = param[1], b = param[2], c = param[3], d = param[4])
        E = SNN.IZ(; N = 1, param = pp)
    end
    if cell_type == "ADEXP"
        adparam = SNN.ADEXParameter(;
            a = param[1],
            b = param[2],
            cm = param[3],
            v0 = param[4],
            τ_m = param[5],
            τ_w = param[6],
            θ = param[7],
            delta_T = param[8],
            v_reset = param[9],
            spike_delta = param[10],
        )

        E = SNN.AD(; N = 1, param = adparam)
    end

    SNN.monitor(E, [:v])
    SNN.monitor(E, [:fire])
    ALLEN_DURATION = 2000 * ms
    ALLEN_DELAY = 1000 * ms
    current = current_search(cell_type, param, ngt_spikes)
    E.I = [current * nA]#[param[5]*nA]

    SNN.sim!(
        [E];
        dt = 1 * ms,
        delay = ALLEN_DELAY,
        stimulus_duration = ALLEN_DURATION,
        simulation_duration = ALLEN_DURATION + ALLEN_DELAY + 443ms,
    )
    vecplot(E, :v) |> display
    error = loss(E, ngt_spikes, ground_spikes)
    error
end

function z(x::AbstractVector)
    param = x
    params = SNN.IZParameter(; a = param[1], b = param[2], c = param[3], d = param[4])
    E = SNN.IZ(; N = 1, param = params)
    SNN.monitor(E, [:v])
    SNN.monitor(E, [:fire])
    ALLEN_DURATION = 2000 * ms
    ALLEN_DELAY = 1000 * ms
    E.I = [param[5] * nA]
    SNN.sim!(
        [E];
        dt = 1 * ms,
        delay = ALLEN_DELAY,
        stimulus_duration = ALLEN_DURATION,
        simulation_duration = ALLEN_DURATION + ALLEN_DELAY + 500ms,
    )
    error = loss(E, ngt_spikes)
    @show(error)
    error

end

function checkmodel(param)
    if cell_type == "IZHI"
        pp = SNN.IZParameter(; a = param[1], b = param[2], c = param[3], d = param[4])
        E = SNN.IZ(; N = 1, param = pp)
    end
    if cell_type == "ADEXP"

        adparam = SNN.ADEXParameter(;
            a = param[1],
            b = param[2],
            cm = param[3],
            v_rest = param[4],
            tau_m = param[5],
            tau_w = param[6],
            v_thresh = param[7],
            delta_T = param[8],
            v_reset = param[9],
            spike_height = param[10],
        )
    end
    E = SNN.AD(; N = 1, param = adparam)

    ALLEN_DURATION = 2000 * ms
    ALLEN_DELAY = 1000 * ms
    E.I = [current_search(param, ngt_spikes) * nA]

    SNN.monitor(E, [:v])
    SNN.sim!(
        [E];
        dt = 1 * ms,
        delay = ALLEN_DELAY,
        stimulus_duration = ALLEN_DURATION,
        simulation_duration = ALLEN_DURATION + ALLEN_DELAY + 443ms,
    )
    if vecp
        vec = SNN.vecplot(E, :v)
        vec |> display
        vec
    end

end

ɛ = 0.125
options = GA(
    populationSize = 20,
    ɛ = ɛ,
    selection = ranklinear(1.5),
    crossover = intermediate(1.0),
    mutation = uniform(1.0),
)
@time result = Evolutionary.optimize(
    zz_,
    lower,
    upper,
    initd,
    options,
    Evolutionary.Options(
        iterations = 20,
        successive_f_tol = 25,
        show_trace = true,
        store_trace = true,
    ),
)
fitness = minimum(result)
println("GA: ɛ=$ɛ) => F: $(minimum(result))")
extremum_param = Evolutionary.minimizer(result)
opt_vec = checkmodel(extremum_param)
println("probably jumbled extremum param")
println(extremum_param)
plot(opt_vec) |> display
plot(vmgtt[:], vmgtv[:]) |> display


text:
__precompile__()

module InteractBase

using WebIO, DataStructures, Observables, CSSUtil, Colors, Requires, JSExpr
import JSExpr: JSString
using Compat
using JSON
using Knockout
using Widgets
import Widgets:
    observe,
    AbstractWidget,
    div,
    Widget,
    widget,
    widgettype,
    ObservablePair,
    @map,
    @layout!,
    @output!,
    @display!,
    @widget,
    node,
    components,
    input

import Observables: throttle, _val

export filepicker, datepicker, timepicker, colorpicker, spinbox

export autocomplete, input, dropdown, checkbox, textbox, textarea, button, slider, toggle, togglecontent

export radiobuttons, togglebuttons, tabs, checkboxes, toggles

export latex, alert, confirm, highlight, accordion, tabulator, mask

export settheme!, resettheme!, gettheme, NativeHTML

export slap_design!

abstract type WidgetTheme; end
struct NativeHTML<:WidgetTheme; end

const font_awesome = joinpath(@__DIR__, "..", "assets", "all.js")
const prism_js = joinpath(@__DIR__, "..", "assets", "prism.js")
const prism_css = joinpath(@__DIR__, "..", "assets", "prism.css")
const highlight_css = joinpath(@__DIR__, "..", "assets", "highlight.css")
const nouislider_min_js = joinpath(@__DIR__, "..", "assets", "nouislider.min.js")
const nouislider_min_css = joinpath(@__DIR__, "..", "assets", "nouislider.min.css")
const style_css = joinpath(@__DIR__, "..", "assets", "style.css")

include("classes.jl")
include("backends.jl")
include("widget.jl")
include("widget_utils.jl")
include("input.jl")
include("slider.jl")
include("optioninput.jl")
include("defaults.jl")
include("manipulate.jl")
include("output.jl")
include("modifiers.jl")

include("providers/atom.jl")
include("providers/blink.jl")
include("providers/mux.jl")

end # module


text:
## conversion from the output of Mamba
function PosteriorDraws(chain::Mamba.ModelChains)
    all(search(name, '[') == 0 for name in chain.names) ||
        error("This function can't handle vector and matrix values (yet).")
    keys = map(Symbol, chain.names)
    vars = [vec(chain.value[:, var_index, :])
            for var_index in indices(chain.value,2)]
    PosteriorDraws(keys, vars)
end


text:
c1 = [-1; -4]
c2 = [-2; -3]

dim_x = length(c1)
dim_u = length(c2)

b = [-2; -3]

A1 = [1 -3;
     -1 -3]
A2 = [1 -2;
     -1 -1]

M = 1000;

# Loading the necessary packages
#-------------------------------
using JuMP
using GLPK
using LinearAlgebra
using Test

function build_subproblem()
    sub_problem_model = Model(GLPK.Optimizer)
    @variable(sub_problem_model, u[1:dim_u] >= 0)
    @constraint(sub_problem_model, constr_ref_subproblem[j = 1:size(A2, 2)], dot(A2[:, j], u) >= c2[j])
    return (sub_problem_model, u)
end

# Master Problem Description
# --------------------------

master_problem_model = Model(GLPK.Optimizer)
(sub_problem_model, u) = build_subproblem()

# Variable Definition
# ----------------------------------------------------------------
@variable(master_problem_model, 0 <= x[1:dim_x] <= 1e6, Int)
@variable(master_problem_model, t <= 1e6)

# Objective Setting
# -----------------
@objective(master_problem_model, Max, t)

print(master_problem_model)

iter_num = 0

function benders_lazy_constraint_callback(cb_data)
    global iter_num
    iter_num += 1
    println("Iteration number = ", iter_num)

    x_current = [callback_value(cb_data, x[j]) for j in eachindex(x)]
    fm_current = callback_value(cb_data, t)

    c_sub = b - A1 * x_current
    @objective(sub_problem_model, Min, dot(c1, x_current) + dot(c_sub, u))
    optimize!(sub_problem_model)

    print("\nThe current subproblem model is \n", sub_problem_model)

    t_status_sub = termination_status(sub_problem_model)
    p_status_sub = primal_status(sub_problem_model)

    fs_x_current = objective_value(sub_problem_model)

    u_current = value.(u)

    γ = dot(b, u_current)

    if p_status_sub == MOI.FEASIBLE_POINT && fs_x_current  ≈  fm_current # we are done
        @info("No additional constraint from the subproblem")
        # println("Optimal solution of the original problem found")
        # println("The optimal objective value t is ", fm_current)
        # println("The optimal x is ", x_current)
        # println("The optimal v is ", dual.(constr_ref_subproblem))
    end

    if p_status_sub == MOI.FEASIBLE_POINT && fs_x_current < fm_current
        println("\nThere is a suboptimal vertex, add the corresponding constraint")
        cv = A1' * u_current - c1
        new_optimality_cons = @build_constraint(t + dot(cv, x) <= γ)
        MOI.submit(master_problem_model, MOI.LazyConstraint(cb_data), new_optimality_cons)
    end

    if t_status_sub == MOI.INFEASIBLE_OR_UNBOUNDED
        println("\nThere is an  extreme ray, adding the corresponding constraint")
        ce = A1' * u_current
        new_feasibility_cons = @build_constraint(dot(ce, x) <= γ)
        MOI.submit(master_problem_model, MOI.LazyConstraint(cb_data), new_feasibility_cons)
    end
end

MOI.set(master_problem_model, MOI.LazyConstraintCallback(), benders_lazy_constraint_callback)

optimize!(master_problem_model)

t_status = termination_status(master_problem_model)
p_status = primal_status(master_problem_model)

if p_status == MOI.INFEASIBLE_POINT
    println("The problem is infeasible :-(")
end

if t_status == MOI.INFEASIBLE_OR_UNBOUNDED
    fm_current = M
    x_current = M * ones(dim_x)
end

if p_status == MOI.FEASIBLE_POINT
    fm_current = value(t)
    x_current = Float64[]
        for i in 1:dim_x
        push!(x_current, value(x[i]))
    end
end

println("Status of the master problem is ", t_status,
        "\nwith fm_current = ", fm_current,
        "\nx_current = ", x_current)

@test value(t) ≈ -4


text: transD_GP.SkyTEM1DInversion.splitlineconvandlast(soundings, 10, 2; zstart=zstart, extendfrac=extendfrac, dz=dz, nlayers=nlayers, figsize=(15,7))


text:
using NLopt
using ForwardDiff

function a_next(c, h, states, agg_state, 𝔼V, params, hh::Renter)
  @unpack ρ, r = agg_state
  @unpack w, z = states
  a = w

  a_next = (1+r) * (a + z - c - ρ * h)
end

function objective0(c, h, states, agg_state, 𝔼V, params, hh::Renter)
  @unpack β = params

  a_next_ = a_next(c, h, states, agg_state, 𝔼V, params, hh)

  uu = c > 0 && h > 0 ? u(c,h) : 10_000 * min(c,h) - 10_000

  uu + β * 𝔼V(a_next_)

end

function constraint0(c, h, states, agg_state, 𝔼V, params, hh::Renter)
  @unpack r = agg_state

  a_next_ = a_next(c, h, states, agg_state, 𝔼V, params, hh)
  # - y_min / r <= a
  #(1+r) * m <= p * h * (1-δ) * θ
  - a_next_
end

function objective_nlopt(x::Vector, grad::Vector, args...)
  if length(grad) > 0
    ForwardDiff.gradient!(grad, XX -> objective0(XX[1], XX[2], args...), x)
  end
  objective0(x[1], x[2], args...)
end

function constraint_nlopt(x::Vector, grad::Vector, args...)
  if length(grad) > 0
    ForwardDiff.gradient!(grad, xx -> constraint0(xx[1], xx[2], args...), x)
  end
  constraint0(x[1], x[2], args...)
end

function Aiyagari.get_optimum(states, agg_state, 𝔼V, params, a_grid, hh::Renter)
  opt = Opt(:LD_MMA, 2)
  #opt = Opt(:LD_SLSQP, 2)
  lower_bounds!(opt, [eps(), eps()])
  upper_bounds!(opt, [+Inf, params.h_thres])

  xtol_rel!(opt, √eps())
  ftol_rel!(opt, eps())

  maxeval!(opt, 150)

  max_objective!(opt, (x,g) -> objective_nlopt(x, g, states, agg_state, 𝔼V, params, hh))
  inequality_constraint!(opt, (x,g) -> constraint_nlopt(x, g, states, agg_state, 𝔼V, params, hh), eps())

  guess = (states.w + states.z)/2
  (max_f, max_x, ret) = NLopt.optimize(opt, [guess, min(guess / agg_state.ρ, params.h_thres)])

  val = max_f
  c, h = max_x
  w_ = a_next(c, h, states, agg_state, 𝔼V, params, hh)

  conv = ret in [:FTOL_REACHED, :XTOL_REACHED, :SUCCESS, :LOCALLY_SOLVED]

  pol_full = (c=c, h=h, w_next=w_, ret=ret, conv=conv, count= opt.numevals)
  pol = pol_full.w_next, pol_full.h

  (pol=pol, pol_full=pol_full, val=val, conv=all(conv))
end


text:
# implementation of algorithms from plotly.py/packages/python/plotly/plotly/subplots.py

# Constants
# ---------
# Subplot types that are each individually positioned with a domain
#
# Each of these subplot types has a `domain` property with `x`/`y`
# properties.
# Note that this set does not contain `xaxis`/`yaxis` because these behave a
# little differently.

const _single_subplot_types = Set(["scene", "geo", "polar", "ternary", "mapbox"])
const _subplot_types = union(_single_subplot_types, Set(["xy", "domain"]))

# For most subplot types, a trace is associated with a particular subplot
# using a trace property with a name that matches the subplot type. For
# example, a `scatter3d.scene` property set to `'scene2'` associates a
# scatter3d trace with the second `scene` subplot in the figure.
#
# There are a few subplot types that don't follow this pattern, and instead
# the trace property is just named `subplot`.  For example setting
# the `scatterpolar.subplot` property to `polar3` associates the scatterpolar
# trace with the third polar subplot in the figure
const _subplot_prop_named_subplot = Set(["polar", "ternary", "mapbox"])

@with_kw struct SubplotRef
    subplot_kind::String
    layout_keys::Vector{Symbol}
    trace_kwargs::PlotlyAttribute
end

const GridRef = Matrix{Vector{SubplotRef}}

function _get_initial_max_subplot_ids()
    max_subplot_ids = Dict(subplot_kind => 0 for subplot_kind in _single_subplot_types)
    max_subplot_ids["xaxis"] = 0
    max_subplot_ids["yaxis"] = 0
    return max_subplot_ids
end

"""
**Parameters**

- `kind` Subplot type. One of
    - 'xy': 2D Cartesian subplot type for scatter, bar, etc.
    - 'scene': 3D Cartesian subplot for scatter3d, cone, etc.
    - 'polar': Polar subplot for scatterpolar, barpolar, etc.
    - 'ternary': Ternary subplot for scatterternary
    - 'mapbox': Mapbox subplot for scattermapbox
    - 'domain': Subplot type for traces that are individually positioned. pie, parcoords, parcats, etc.
    - Trace type: Put the name of the type of trace here and we will determine the appropriate kind of subplot

- `secondary_y`: If true, create a secondary y-axis positioned on the right side of the subplot. Only valid if kind="xy".
- `colspan`: number of subplot columns for this subplot to span.
- `rowspan`: number of subplot rows for this subplot to span.
- `l`: padding left of cell
- `r`: padding right of cell
- `t`: padding right of cell
- `b`: padding bottom of cell
"""
@with_kw struct Spec
    kind::String = "xy"
    secondary_y::Bool = false
    colspan::Int = 1
    rowspan::Int = 1
    l::Float64 = 0.0
    r::Float64 = 0.0
    b::Float64 = 0.0
    t::Float64 = 0.0

    @assert !secondary_y || (kind == "xy")
end

"""
**Parameters**

- `cell`: index of the subplot cell to overlay inset axes onto.
- `kind`: Subplot kind (see `Spec` docs)
- `l`: padding left of inset in fraction of cell width
- `w`: inset width in fraction of cell width ('to_end': to cell right edge)
- `b`: padding bottom of inset in fraction of cell height
- `h`: inset height in fraction of cell height ('to_end': to cell top edge)
"""
@with_kw struct Inset
    cell::Tuple{Int,Int} = (1, 1)
    kind::String = "xy"
    colspan::Int = 1
    rowspan::Int = 1
    l::Float64 = 0.0
    w::Union{String,Float64} = "to_end"
    b::Float64 = 0.0
    h::Union{String,Float64} = "to_end"
end

function _check_hv_spacing(dimsize::Int, spacing::Float64, name::String, dimvarname::String)
    spacing < 0 || spacing > 1 && error("$name spacing must be between 0 and 1")
    if dimsize <= 1
        return true
    end
    max_spacing = 1.0 / (dimsize - 1)
    if spacing > max_spacing
        error("$name spacing cannot be greater than 1/($dimvarname - 1) = $max_spacing")
    end
    return true
end

"""
**Parameters**

- `rows`: Number of rows in the subplot grid. Must be greater than zero.
- `cols`: Number of columns in the subplot grid. Must be greater than zero.
- `shared_xaxes`: Assign shared (linked) x-axes for 2D cartesian subplots
    - true or "columns": Share axes among subplots in the same column
    - "rows": Share axes among subplots in the same row
    - "all": Share axes across all subplots in the grid
- `shared_yaxes`: Assign shared (linked) y-axes for 2D cartesian subplots
    - "columns": Share axes among subplots in the same column
    - true or "rows": Share axes among subplots in the same row
    - "all": Share axes across all subplots in the grid.
- `start_cell`:  Choose the starting cell in the subplot grid used to set the domains_grid of the subplots.
    - "top-left": Subplots are numbered with (1, 1) in the top left corner
    - "bottom-left": Subplots are numbererd with (1, 1) in the bottom left corner
- `horizontal_spacing`: Space between subplot columns in normalized plot coordinates. Must be a float between 0 and 1. Applies to all columns (use "specs" subplot-dependents spacing)
- `vertical_spacing`: Space between subplot rows in normalized plot coordinates. Must be a float between 0 and 1. Applies to all rows (use "specs" subplot-dependents spacing)
- `subplot_titles`: Title of each subplot as a list in row-major ordering. Empty strings ("") can be included in the list if no subplot title is desired in that space so that the titles are properly indexed.
- `specs`:  Per subplot specifications of subplot type, row/column spanning, and spacing.
    - The number of rows in "specs" must be equal to "rows".
    - The number of columns in "specs"
    - Each item in the "specs" list corresponds to one subplot
      in a subplot grid. (N.B. The subplot grid has exactly "rows"
      times "cols" cells.)
    - Use missing for a blank a subplot cell (or to move past a col/row span).
    - Each item in "specs" is an instance of `Spec`. See docs for `Spec` for more information
    - Note: Use `horizontal_spacing` and `vertical_spacing` to adjust the spacing in between the subplots.
- `insets`: Inset specifications.  Insets are subplots that overlay grid subplots
    - Each item in "insets" is an instance of `Inset`. See docs for `Inset` for more info
- `column_widths`:  Array of length `cols` of the relative widths of each column of suplots. Values are normalized internally and used to distribute overall width of the figure (excluding padding) among the columns.
- `row_heights`: Array of length `rows` of the relative heights of each row of subplots. Values are normalized internally and used to distribute overall height of the figure (excluding padding) among the rows
- `column_titles`: list of length `cols` of titles to place above the top subplot in each column.
- `row_titles`: list of length `rows` of titles to place on the right side of each row of subplots.
- `x_title`: Title to place below the bottom row of subplots, centered horizontally
- `y_title`: Title to place to the left of the left column of subplots, centered vertically
"""
@with_kw struct Subplots
    rows::Int = 1
    cols::Int = 1
    shared_xaxes::_Maybe{Union{String,Bool}} = false
    shared_yaxes::_Maybe{Union{String,Bool}} = false
    start_cell::String = "top-left"
    subplot_titles::_Maybe{Matrix{<:_Maybe{String}}} = missing
    column_widths::_Maybe{Vector{Float64}} = missing
    row_heights::_Maybe{Vector{Float64}} = missing
    specs::Matrix{<:_Maybe{Spec}} = [Spec() for _ in 1:rows, __ in 1:cols]
    insets::_Maybe{Union{Bool,Vector{Inset}}} = missing
    column_titles::_Maybe{Vector{String}} = missing
    row_titles::_Maybe{Vector{String}} = missing
    x_title::_Maybe{String} = missing
    y_title::_Maybe{String} = missing
    grid_ref::GridRef = let
        gr = GridRef(undef, (rows, cols))
        if rows == 1 && cols == 1
            gr[1, 1] = [SubplotRef(subplot_kind="xy", layout_keys=[:xaxis, :yaxis], trace_kwargs=attr(xaxis="x", yaxis="y"))]
        end
        gr
    end

    # computed
    has_secondary_y::Bool = any(!ismissing(s) && getfield(s, :secondary_y) for s in specs)
    horizontal_spacing::Float64 = has_secondary_y ? 0.4 / cols : 0.2 / cols
    vertical_spacing::Float64 = ismissing(subplot_titles) ? 0.3 / rows : 0.5 / rows
    max_width::Float64 = has_secondary_y ? 0.94 : (ismissing(row_titles) ? 1.0 : 0.98)
    _widths::Vector{Float64} = let
        if ismissing(column_widths)
            fill((max_width - horizontal_spacing * (cols - 1)) / cols, cols)
        else
            w_sum = sum(column_widths)
            [(max_width - horizontal_spacing * (cols - 1)) * (w / w_sum) for w in column_widths]
        end
    end
    _heights::Vector{Float64} = let
        if ismissing(row_heights)
            fill((1.0 - vertical_spacing * (rows - 1)) / rows, rows)
        else
            h_sum = sum(row_heights)
        [(1.0 - vertical_spacing * (rows - 1)) * (h / h_sum) for h in row_heights]
        end
    end

    @assert start_cell in ["top-left", "bottom-left"]
    @assert ismissing(shared_xaxes) || (shared_xaxes in [true, false, "rows", "columns", "all"])
    @assert ismissing(shared_yaxes) || (shared_yaxes in [true, false, "rows", "columns", "all"])
    @assert _check_hv_spacing(cols, horizontal_spacing, "Horizontal", "cols")
    @assert _check_hv_spacing(rows, vertical_spacing, "Vertical", "rows")
    @assert length(_widths) == cols
    @assert length(_heights) == rows
    @assert ismissing(column_titles) || (length(column_titles) == cols)
    @assert ismissing(row_titles) || (length(row_titles) == rows)

end


text:
using Compat.Test, MathOptInterface
using LinQuadOptInterface

const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const LQOI = LinQuadOptInterface

@testset "LinQuadOptInterface" begin
    solver = LQOI.MockLinQuadOptimizer()
    # TODO(@joaquim): test with solve=true
    config = MOIT.TestConfig(solve=false)

    @testset "Unit Tests" begin
        MOIT.basic_constraint_tests(solver, config)
        MOIT.unittest(solver, config, [
            "solve_affine_interval",
            "solve_qp_edge_cases",
            "solve_qcp_edge_cases",
            "solve_affine_deletion_edge_cases"
        ])
    end

    @testset "Linear tests" begin
        MOIT.contlineartest(solver, config)
    end

    @testset "Quadratic tests" begin
        MOIT.contquadratictest(solver, config)
    end

    @testset "Linear Conic tests" begin
        MOIT.lintest(solver, config)
    end

    @testset "Integer Linear tests" begin
        MOIT.intlineartest(solver, config, ["int2"])
    end

    @testset "ModelLike tests" begin
        @testset "nametest" begin
            MOIT.nametest(LQOI.MockLinQuadOptimizer())
        end
        @testset "validtest" begin
            MOIT.validtest(LQOI.MockLinQuadOptimizer())
        end
        @testset "emptytest" begin
            MOIT.emptytest(LQOI.MockLinQuadOptimizer())
        end
        @testset "orderedindicestest" begin
            MOIT.orderedindicestest(LQOI.MockLinQuadOptimizer())
        end
        @testset "copytest" begin
            MOIT.copytest(LQOI.MockLinQuadOptimizer(),
                          LQOI.MockLinQuadOptimizer())
        end
    end
end

@testset "Issue #52" begin
    @testset "scalaraffine" begin
        model = LQOI.MockLinQuadOptimizer()
        x = MOI.addvariable!(model)
        f1 =  MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 1.5)
        s1 = MOI.LessThan(2.0)
        c = MOI.addconstraint!(model, f1, s1)
        f2 = MOI.get(model, MOI.ConstraintFunction(), c)
        @test f1 ≈ f2
        s2 = MOI.get(model, MOI.ConstraintSet(), c)
        @test typeof(s1) == typeof(s2)
        @test s1.upper == s2.upper
    end
    @testset "vectoraffine" begin
        model = LQOI.MockLinQuadOptimizer()
        x = MOI.addvariable!(model)
        f1 =  MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.([1], MOI.ScalarAffineTerm.([1.0], [x])), [1.5])
        s1 = MOI.Zeros(1)
        c = MOI.addconstraint!(model, f1, s1)
        f2 = MOI.get(model, MOI.ConstraintFunction(), c)
        @test f1 ≈ f2
        s2 = MOI.get(model, MOI.ConstraintSet(), c)
        @test typeof(s1) == typeof(s2)
        @test s1.dimension == s2.dimension
    end
end

@testset "Issue #54" begin
    @testset "scalar, one-sided" begin
        model = LQOI.MockLinQuadOptimizer()
        x = MOI.addvariable!(model)
        f =  MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 1.5)
        s = MOI.LessThan(2.0)
        c = MOI.addconstraint!(model, f, s)
        # Change the constraint set, and verify that we get the same set
        # when we retrieve it:
        s2 = MOI.LessThan(1.0)
        MOI.set!(model, MOI.ConstraintSet(), c, s2)
        s3 = MOI.get(model, MOI.ConstraintSet(), c)
        @test typeof(s2) == typeof(s3)
        @test s2.upper == s3.upper
    end
    @testset "scalar, interval" begin
        model = LQOI.MockLinQuadOptimizer()
        x = MOI.addvariable!(model)
        f =  MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0], [x]), 1.5)
        s = MOI.Interval(2.0, 3.0)
        c = MOI.addconstraint!(model, f, s)
        # Change the constraint set, and verify that we get the same set
        # when we retrieve it:
        s2 = MOI.Interval(1.0, 2.0)
        MOI.set!(model, MOI.ConstraintSet(), c, s2)
        s3 = MOI.get(model, MOI.ConstraintSet(), c)
        @test typeof(s2) == typeof(s3)
        @test s2.lower == s3.lower
        @test s2.upper == s3.upper
    end
end


text:
module StableApproxEPW

using LinearAlgebra, SparseArrays, SpecialFunctions
using FastGaussQuadrature, Roots
using QuasiMonteCarlo

include("adaptive-quad-rule.jl")
export ϵsupport, adaptive_G_quad

include("circular-waves.jl")
export b̃p, βp, bp, solution_surrogate

include("herglotz-densities.jl")
export ãp, αp, ap, wz, logwz

include("plane-waves.jl")
export gpw, approximation_set

include("parametric_sampling.jl")
export TruncKernel, kernel_diagonal
export probabilitydensityfunction, probabilitydensityfunction_weighted
export cumulativedensityfunction
export Sampler, weight, inversion
export random_sampling, uniform_sampling, sobol_sampling

include("regularizedSVD.jl")
export RegularizedSVDPseudoInverse, condition_number, solve_via_regularizedSVD

include("dirichlet_sampling.jl")
export samples_from_nodes, boundary_sampling_nodes
export number_of_boundary_sampling_nodes
export Dirichlet_sampling

end

max_stars_repo_name: EmileParolin/STrAW

max_issues_repo_name: EmileParolin/STrAW

max_forks_repo_name: EmileParolin/STrAW


text:
@testset "Structured Matrices" begin
    @testset "Diagonal" begin
        N = 3
        rrule_test(Diagonal, randn(N, N), (randn(N), randn(N)))
        D = Diagonal(randn( N))
        rrule_test(Diagonal, D, (randn(N), randn(N)))
        # Concrete type instead of UnionAll
        rrule_test(typeof(D), D, (randn(N), randn(N)))

        # TODO: replace this with a `rrule_test` once we have that working
        # see https://github.com/JuliaDiff/ChainRulesTestUtils.jl/issues/24
        res, pb = rrule(Diagonal, [1, 4])
        @test pb(10*res) == (NO_FIELDS, [10, 40])
        comp = Composite{typeof(res)}(; diag=10*res.diag)  # this is the structure of Diagonal
        @test pb(comp) == (NO_FIELDS, [10, 40])
    end

    @testset "::Diagonal * ::AbstractVector" begin
        N = 3
        rrule_test(
            *,
            randn(N),
            (Diagonal(randn(N)), Diagonal(randn(N))),
            (randn(N), randn(N)),
        )
    end
    @testset "diag" begin
        N = 7
        rrule_test(diag, randn(N), (randn(N, N), randn(N, N)))
        rrule_test(diag, randn(N), (Diagonal(randn(N)), randn(N, N)))
        rrule_test(diag, randn(N), (randn(N, N), Diagonal(randn(N))))
        rrule_test(diag, randn(N), (Diagonal(randn(N)), Diagonal(randn(N))))
        VERSION ≥ v"1.3" && @testset "k=$k" for k in (-1, 0, 2)
            M = N - abs(k)
            rrule_test(diag, randn(M), (randn(N, N), randn(N, N)), (k, nothing))
        end
    end
    @testset "diagm" begin
        @testset "without size" begin
            M, N = 7, 9
            s = (8, 8)
            a, ā = randn(M), randn(M)
            b, b̄ = randn(M), randn(M)
            c, c̄ = randn(M - 1), randn(M - 1)
            ȳ = randn(s)
            ps = (0 => a, 1 => b, 0 => c)
            y, back = rrule(diagm, ps...)
            @test y == diagm(ps...)
            ∂self, ∂pa, ∂pb, ∂pc = back(ȳ)
            @test ∂self === NO_FIELDS
            ∂a_fd, ∂b_fd, ∂c_fd = j′vp(_fdm, (a, b, c) -> diagm(0 => a, 1 => b, 0 => c), ȳ, a, b, c)
            for (p, ∂px, ∂x_fd) in zip(ps, (∂pa, ∂pb, ∂pc), (∂a_fd, ∂b_fd, ∂c_fd))
                ∂px = unthunk(∂px)
                @test ∂px isa Composite{typeof(p)}
                @test ∂px.first isa AbstractZero
                @test ∂px.second ≈ ∂x_fd
            end
        end
        VERSION ≥ v"1.3" && @testset "with size" begin
            M, N = 7, 9
            a, ā = randn(M), randn(M)
            b, b̄ = randn(M), randn(M)
            c, c̄ = randn(M - 1), randn(M - 1)
            ȳ = randn(M, N)
            ps = (0 => a, 1 => b, 0 => c)
            y, back = rrule(diagm, M, N, ps...)
            @test y == diagm(M, N, ps...)
            ∂self, ∂M, ∂N, ∂pa, ∂pb, ∂pc = back(ȳ)
            @test ∂self === NO_FIELDS
            @test ∂M === DoesNotExist()
            @test ∂N === DoesNotExist()
            ∂a_fd, ∂b_fd, ∂c_fd = j′vp(_fdm, (a, b, c) -> diagm(M, N, 0 => a, 1 => b, 0 => c), ȳ, a, b, c)
            for (p, ∂px, ∂x_fd) in zip(ps, (∂pa, ∂pb, ∂pc), (∂a_fd, ∂b_fd, ∂c_fd))
                ∂px = unthunk(∂px)
                @test ∂px isa Composite{typeof(p)}
                @test ∂px.first isa AbstractZero
                @test ∂px.second ≈ ∂x_fd
            end
        end
    end
    @testset "Symmetric" begin
        N = 3
        rrule_test(Symmetric, randn(N, N), (randn(N, N), randn(N, N)))
    end
    @testset "$f" for f in (Adjoint, adjoint, Transpose, transpose)
        n = 5
        m = 3
        rrule_test(f, randn(m, n), (randn(n, m), randn(n, m)))
        rrule_test(f, randn(1, n), (randn(n), randn(n)))
    end
    @testset "$T" for T in (UpperTriangular, LowerTriangular)
        n = 5
        rrule_test(T, T(randn(n, n)), (randn(n, n), randn(n, n)))
    end
    @testset "$Op" for Op in (triu, tril)
        n = 7
        rrule_test(Op, randn(n, n), (randn(n, n), randn(n, n)))
        @testset "k=$k" for k in -2:2
            rrule_test(Op, randn(n, n), (randn(n, n), randn(n, n)), (k, nothing))
        end
    end
end


text:
# # Illustration of Jacobi Fields
#
# This tutorial illustrates the usage of Jacobi Fields within
# `Manopt.jl`.
# For this tutorial you should be familiar with the basic terminology on a
# manifold like the exponential and logarithmic map as well as
# [`geodesic`](@ref)s.
#
# We first initialize the manifold
exportFolder = joinpath(@__DIR__,"..","..","docs","src","assets","images","tutorials") #src
using Manopt
# and we define some colors from [Paul Tol](https://personal.sron.nl/~pault/)
using Colors
black = RGBA{Float64}(colorant"#000000")
TolVibrantOrange = RGBA{Float64}(colorant"#EE7733")
TolVibrantCyan = RGBA{Float64}(colorant"#33BBEE")
TolVibrantTeal = RGBA{Float64}(colorant"#009988")
nothing #hide
# Assume we have two [`SnPoint`](@ref)s on the equator of the [`Sphere`](@ref)`(2)` $\mathcal M = \mathbb S^2$
M = Sphere(2)
x,y = [ SnPoint([1.,0.,0.]), SnPoint([0.,1.,0.])]
# their connecting [`geodesic`](@ref) (sampled at `100` points)
geodesicCurve = geodesic(M,x,y,100);
asyResolution = 2
nothing #hide
# looks as follows using [`renderAsymptote`](@ref) with the [`asyExportS2Signals`](@ref) export
renderAsymptote(exportFolder*"/jacobiGeodesic.asy",asyExportS2Signals; #src
    render = asyResolution, #src
    curves=[geodesicCurve], points = [ [x,y] ], #src
    colors=Dict(:curves => [black], :points => [TolVibrantOrange]), #src
    dotSize = 3.5, lineWidth = 0.75, cameraPosition = (1.,1.,.5) #src
)#src
#md #
#md # ```julia
#md # renderAsymptote("jacobiGeodesic.asy",asyExportS2Signals;
#md #     render = asyResolution,
#md #     curves=[geodesicCurve], points = [ [x,y] ],
#md #     colors=Dict(:curves => [black], :points => [TolVibrantOrange]),
#md #     dotSize = 3.5, lineWidth = 0.75, cameraPosition = (1.,1.,.5)
#md # )
#md # ```
#md #
#md # ![A geodesic connecting two points on the equator](../assets/images/tutorials/jacobiGeodesic.png)
#
# where $x$ is on the left. Then this tutorial solves the following task:
#
# Given a direction $\xi_x\in T_x\mathcal M$, for example the [`SnTVector`](@ref)
ξx = SnTVector([0.,0.4,0.5])
# we move the start point $x$ into, how does any point on the geodesic move?
#
# Or mathematically: Compute $D_x g(t; x,y)$ for some fixed $t\in[0,1]$
# and a given direction $\xi_x$.
# Of course two cases are quite easy: For $t=0$ we are in $x$ and how $x$ “moves”
# is already known, so $D_x g(0;x,y) = \xi$. On the other side, for $t=1$,
# $g(1; x,y) = y$ which is fixed, so $D_x g(1; x,y)$ is the zero tangent vector
# (in $T_y\mathcal M$).
#
# For all other cases we employ a [`jacobiField`](@ref), which is a (tangent)
# vector field along the [`geodesic`](@ref) given as follows: The _geodesic variation_
# $\Gamma_{g,\xi}(s,t)$ is defined for some $\varepsilon > 0$ as
#
# $\Gamma_{g,\xi}(s,t):=\exp{\gamma_{x,\xi}(s)}[t\log_{g(s;x,\xi)}y],\qquad s\in(-\varepsilon,\varepsilon),\ t\in[0,1].$
#
# Intuitively we make a small step $s$ into direction $\xi$ using the geodesic
# $g(\cdot; x,\xi)$ and from $z=g(s; x,\xi)$ we follow (in $t$) the geodesic
# $g(\cdot; z,y)$. The corresponding Jacobi field~\(J_{g,\xi}\)
# along~\(g(\cdot; x,y\) is given
#
# $J_{g,\xi}(t):=\frac{D}{\partial s}\Gamma_{g,\xi}(s,t)\Bigl\rvert_{s=0}$
#
# which is an ODE and we know the boundary conditions $J_{g,\xi}(0)=\xi$ and
# $J_{g,\xi}(t) = 0$. In symmetric spaces we can compute the solution, since the
# system of ODEs decouples, see for example [do Carmo](#doCarmo1992),
# Chapter 4.2. Within `Manopt.jl` this is implemented as
# [`jacobiField`](@ref)`(M,x,y,t,ξ[,β])`, where the optional parameter (function)
# `β` specifies, which Jacobi field we want to evaluate and the one used here is
# the default.
#
# We can hence evaluate that on the points on the geodesic at
T = [0:0.1:1.0...]
nothing #hide
# namely
Z = geodesic(M,x,y,T)
nothing #hide
# the geodesic moves as
ηx = jacobiField.(Ref(M), Ref(x), Ref(y), T, Ref(ξx) )
# which can also be called using [`DxGeo`](@ref).
# We can add to the image above by creating extended tangent vectors
# [`TVectorE`](@ref) the include their base points
Vx = TVectorE.(ηx,Z)
# and add that as one further set to the Asymptote export.
renderAsymptote(exportFolder*"/jacobiGeodesicDxGeo.asy",asyExportS2Signals; #src
    render = asyResolution, #src
    curves=[geodesicCurve], points = [ [x,y], Z], tVectors = [Vx], #src
    colors=Dict( #src
        :curves => [black], #src
        :points => [TolVibrantOrange,TolVibrantCyan], #src
        :tvectors => [TolVibrantCyan] #src
    ), #src
    dotSizes = [3.5,2.], lineWidth = 0.75, cameraPosition = (1.,1.,.5) #src
) #src
#md #
#md # ```julia
#md # renderAsymptote("jacobiGeodesicDxGeo.asy",asyExportS2Signals;
#md #     render = asyResolution,
#md #     curves=[geodesicCurve], points = [ [x,y], Z], tVectors = [Vx],
#md #     colors=Dict(
#md #         :curves => [black],
#md #         :points => [TolVibrantOrange,TolVibrantCyan],
#md #         :tvectors => [TolVibrantCyan]
#md #     ),
#md #     dotSizes = [3.5,2.], lineWidth = 0.75, cameraPosition = (1.,1.,.5)
#md # )
#md # ```
#
#md # ![A Jacobi field for $D_xg(t,x,y)[\eta]$](../assets/images/tutorials/jacobiGeodesicDxGeo.png)
#
# If we further move the end point, too, we can derive that Differential in direction
ξy = SnTVector([0.2,0.,-0.5])
ηy = DyGeo.(Ref(M),Ref(x),Ref(y),T,Ref(ξy))
Vy = TVectorE.(ηy,Z)
# and we can look at the total effect, where the [`TVectorE`](@ref)s even verify
# that only tangent vectors are added that have a common base point
Vb = Vx .+ Vy
renderAsymptote(exportFolder*"/jacobiGeodesicResult.asy",asyExportS2Signals; #src
    render = asyResolution, #src
    curves=[geodesicCurve], points = [ [x,y], Z], tVectors = [Vx,Vy,Vb], #src
    colors=Dict( #src
        :curves => [black], #src
        :points => [TolVibrantOrange,TolVibrantCyan], #src
        :tvectors => [TolVibrantCyan,TolVibrantCyan,TolVibrantTeal] #src
    ), #src
    dotSizes = [3.5,2.], lineWidth = 0.75, cameraPosition = (1.,1.,0.) #src
) #src
#md # ```julia
#md # renderAsymptote("jacobiGeodesicResult.asy",asyExportS2Signals;
#md #    render = asyResolution,
#md #    curves=[geodesicCurve], points = [ [x,y], Z], tVectors = [Vx,Vy,Vb],
#md #    colors=Dict(
#md #        :curves => [black],
#md #        :points => [TolVibrantOrange,TolVibrantCyan],
#md #        :tvectors => [TolVibrantCyan,TolVibrantCyan,TolVibrantTeal]
#md #   ),
#md #   dotSizes = [3.5,2.], lineWidth = 0.75, cameraPosition = (1.,1.,0.)
#md # )
#md # ```
#md #
#md # ![A Jacobi field for the effect of two differentials (blue) in sum (teal)](../assets/images/tutorials/jacobiGeodesicResult.png)
#
# ## Literature
#
# ```@raw html
# <ul><li id="doCarmo1992">[<a>doCarmo1992</a>] do Carmo, M. P.:
#    <emph>Riemannian Geometry</emph>, Mathematics: Theory & Applications,
#    Birkhäuser Basel, 1992, ISBN: 0-8176-3490-8</li>
# <li id="BergmannGousenbourger2018">[<a>BergmannGousenbourger2018</a>]
#   Bergmann, R.; Gousenbourger, P.-Y.: <emph>A variational model for data
#   fitting on manifolds by minimizing the acceleration of a Bézier curve</emph>,
#   Frontiers in Applied Mathematics and Statistics, 2018.
#   doi: <a href="https://dx.doi.org/10.3389/fams.2018.00059">10.3389/fams.2018.00059</a></li>
# </ul>
# ```


text:
module Aardvark

    export EMD
    export emd

    using Distances

    struct EMD <: Distances.Metric
        ground_distance::AbstractArray{<:AbstractFloat, 2}
        max_iter::Integer
        raise_error::Bool
        function EMD(ground_distance, max_iter = 100000, raise_error=true)
            new(ground_distance, max_iter, raise_error)
        end
    end

    function (dist::EMD)(x, y)
        # I do not know how to pass a slice to Cxx so for now we need to make sure
        # we are not getting slices passed here
        return emd(Array(x),
                   Array(y),
                   dist.ground_distance,
                   dist.max_iter)
    end

    module CppEMD
        using CxxWrap

        filename = joinpath(string(@__DIR__), "../lib/libEMD.so")
        @wrapmodule(filename)

        function __init__()
            @initcxx
        end
    end

    function emd(a::Array{<:AbstractFloat, 1},
                 b::Array{<:AbstractFloat, 1},
                 M::Array{<:AbstractFloat, 2},
                 max_iter::Integer=100000;
                 raise_error = true,
                 log = false)

        """
        raise_error: Whether to raise an error or to silently return -1.0

        Example of a problem without a solution
        a = [92  0.0130435
             93  0.00434783
             95  0.0347826
             96  0.0695652
             97  0.0652174
             98  0.0434783
             99  0.0434783
             100 0.926087]
        b = [92   0.0130435
             93   0.00434783
             95   0.0347826
             96   0.0623188
             97   0.0652174
             98   0.0869565
             99   0.0434783
             100  0.889855]
        a = sparsevec(a[:, 1], a[:, 2])
        b = sparsevec(b[:, 1], b[:, 2])
        ground_distance = hist.hist_to_distance_matrix(100)

        """

        # We can't pass views here
        # Cxx seems to take care of the C major format of the arrays though
        n1 = size(M)[1]
        n2 = size(M)[2]
        nmax = n1 + n2 - 1
        result_code = 0
        nG = 0

        cost = Ref(0.0)
        alpha = zeros(n1)
        beta = zeros(n2)
        G = zeros((1, 1))
        Gv = zeros()
        iG = zeros()
        jG = zeros()

        G = zeros((n1, n2))

        result_code = CppEMD.EMD_wrap(n1, n2, a, b, M, G, alpha, beta, cost, max_iter)
        if result_code == 0 && raise_error
            error("Problem not feasible.")
        elseif result_code == 0
            cost = Ref(-1.0)
        end
        if log
            return G, cost, alpha, beta, result_code
        else
            return cost[]
        end
end

end