text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antddropdown

"""
    'feffery'_antddropdown(;kwargs...)

An AntdDropdown component.

Keyword arguments:
- `id` (String; optional)
- `arrow` (Bool; optional)
- `buttonMode` (Bool; optional)
- `className` (String; optional)
- `clickedKey` (String; optional)
- `disabled` (Bool; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `menuItems` (optional): . menuItems has the following type: Array of lists containing elements 'title', 'href', 'target', 'disabled', 'icon', 'key', 'isDivider'.
Those elements have the following types:
  - `title` (String; optional)
  - `href` (String; optional)
  - `target` (String; optional)
  - `disabled` (Bool; optional)
  - `icon` (String; optional)
  - `key` (String; optional)
  - `isDivider` (Bool; optional)s
- `overlayClassName` (String; optional)
- `overlayStyle` (Dict; optional)
- `placement` (a value equal to: 'bottomLeft', 'bottomCenter', 'bottomRight', 'topLeft', 'topCenter', 'topRight'; optional)
- `style` (Dict; optional)
- `title` (String; optional)
- `trigger` (a value equal to: 'click', 'hover'; optional)
- `visible` (Bool; optional)
"""
function 'feffery'_antddropdown(; kwargs...)
        available_props = Symbol[:id, :arrow, :buttonMode, :className, :clickedKey, :disabled, :loading_state, :menuItems, :overlayClassName, :overlayStyle, :placement, :style, :title, :trigger, :visible]
        wild_props = Symbol[]
        return Component("'feffery'_antddropdown", "AntdDropdown", "feffery_antd_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: ReilYoung/feffery-antd-components

max_issues_repo_name: ReilYoung/feffery-antd-components

max_forks_repo_name: ReilYoung/feffery-antd-components


text:
if in("../src", LOAD_PATH)
    push!(LOAD_PATH,"../src")
end

using ccall
using Base.Test

# write your own tests here
@test 1 == 1

ext: jl


text:
module Learning

abstract type AbstractValueFunction end

mutable struct ValueFunction{T<:AbstractArray, F<:Function, F′<:Function}
    w::T
    σ::F
    σ′::F′
end

(vf::VF)(x) where {VF<:ValueFunction} = vf.σ(vf.w'x)
deriv(vf::VF)(x) where {VF<:ValueFunction} = vf.σ′(vf.w'x).*x

# mutable struct DeepValueFunction{M}
#     model::M
# end

abstract type AbstractLearningUpdate end

struct TD{T<:Number} <: AbstractLearningUpdate
    γ::T
end

function update!(model,
                 lu::TDLearning,
                 opt,
                 s_t::Array{<:AbstractFloat, 1},
                 a_t::I,
                 s_tp1::Array{<:AbstractFloat, 1},
                 r::AF,
                 terminal::Bool) where {AF<:AbstractFloat, I<:Integer}

    ps = params(model)
    γ = lu.γ.*(1 .- terminal)

    v_tp1 = model(s_tp1)

    target = (r .+ γ.*v_tp1)
    gs = Flux.gradient(ps) do
        v_t = model(s_t)
        return Flux.mse(target, v_t)
    end
    Flux.Optimise.update!(opt, ps, gs)
end

end

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
if Pkg.installed("Ledger") == nothing
    Pkg.clone("https://github.com/ozanarkancan/Ledger.jl.git")
end

for p in ("Knet","ArgParse")
    Pkg.installed(p) == nothing && Pkg.add(p)
end

using Knet, Ledger

include(Pkg.dir("Knet","data","housing.jl"))

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

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

lossgradient = grad(loss)

function train(w, x, y; lr=.1, epochs=20)
    for epoch=1:epochs
        g = lossgradient(w, x, y)
        update!(w, g; lr=lr)
    end
    return w
end

function experimentf(;seed=-1, epochs=20, lr=0.1, test=0.0, winit=0.1)
    srand(seed)
    w = map(Array{Float32}, [winit*randn(1,13), winit*randn(1,1) ])
    (xtrn,ytrn,xtst,ytst) = map(Array{Float32}, Main.housing(test))
    experiment_log = []

    push!(experiment_log, Dict("epoch" => 0, "train-mse" => loss(w, xtrn, ytrn), "test-mse" => loss(w, xtst, ytst)))
    println(experiment_log[end])

    for i=1:epochs
        train(w, xtrn, ytrn; lr=lr, epochs=1)
        push!(experiment_log, Dict("epoch" => i, "train-mse" => loss(w, xtrn, ytrn), "test-mse" => loss(w, xtst, ytst)))
        println(experiment_log[end])
        sleep(0.5)
    end
    return experiment_log, Dict("test-mse" => loss(w, xtst, ytst))
end

function main()
    engine = LedgerEngine("mongodb://localhost:27017")
    configf = Dict(:seed=>123, :epochs=>100, :lr=>0.1, :test=>0.0, :winit=>0.1)
    Ledger.run(engine, "housing", "linear-regression", experimentf, configf, ["epoch", "train-mse", "test-mse"])
end

PROGRAM_FILE=="knet_linreg.jl" && main()


text:
struct MCopula{d} <: Copula{d} end
MCopula(d) = MCopula{d}()
Distributions.cdf(::MCopula{d},u) where {d} = min(u)
function Distributions._rand!(rng::Distributions.AbstractRNG, ::MCopula{d}, x::AbstractVector{T}) where {d,T<:Real}
    x .= rand(rng)
end
function Base.rand(rng::Distributions.AbstractRNG,::MCopula{d}) where {d}
    repeat([rand(rng)],d)
end


text:
using Test
using CausalityTools
import DelayEmbeddings: Dataset
import DynamicalSystemsBase: DiscreteDynamicalSystem, ContinuousDynamicalSystem

include("methods/test_smeasure.jl")
include("methods/test_joint_distance_distribution.jl")
include("methods/test_predictive_asymmetry.jl")
include("methods/test_crossmapping.jl")

include("systems/discrete/test_discrete_systems.jl")
include("systems/continuous/test_continuous_systems.jl")

#include("integrations/test_uncertaindata_integration.jl")


text:
using RCall, CSV, DataFrames
include("ProblemGenerator.jl")
include("VirusCallBacks.jl")

function DataRangle(sol)
  timeLength = length(sol.t)
  df = DataFrame()
  df.x = repeat(repeat(1:N,N),timeLength)
  df.y = repeat(repeat(1:N,inner=N),timeLength)
  df.t = repeat(sol.t,inner=N^2)
  df.IFN = vec(sol[:,:,7,:])

  return df
end

function CellStateGrid(θ::ParContainer,t::Float64)
  grid = zeros(Int64,size(θ.cellsInfected))
  grid[θ.cellsInfected .<= t] .= 1
  grid[θ.cellsDead .<= t] .= 2

  return vec(grid)
end

#-----------ISD, Not Stochastic, Homo-----------
probISD_nS_Ho = ModelSetup(:ISD,:notStochastic,:Homo)
    solISD_nS_Ho = @time solve(probISD_nS_Ho,CVODE_BDF(linear_solver=:GMRES),saveat=0.1)

ISD_nS_Ho = DataRangle(solISD_nS_Ho)
CSV.write("/home/robertgregg/R/Scripts/ABMMovie/ISD_nS_Ho.csv", ISD_nS_Ho)

#-----------Virus, Not Stochastic, Homo-----------
probVirus_nS_Ho = ModelSetup(:Virus,:notStochastic,:Homo)
    solVirus_nS_Ho = @time solve(probVirus_nS_Ho,CVODE_BDF(linear_solver=:GMRES),saveat=0.1,callback=cb)

Virus_nS_Ho = DataRangle(solVirus_nS_Ho)
CSV.write("/home/robertgregg/R/Scripts/ABMMovie/Virus_nS_Ho.csv", Virus_nS_Ho)

#-----------ISD, Stochastic, Hetero-----------
probISD_S_He = ModelSetup(:ISD,:Stochastic,:Hetero)
    solISD_S_He = @time solve(probISD_S_He,CVODE_BDF(linear_solver=:GMRES),saveat=0.1)

ISD_S_He = DataRangle(solISD_S_He)
CSV.write("/home/robertgregg/R/Scripts/ABMMovie/ISD_S_He.csv", ISD_S_He)

#-----------Virus, Stochastic, Hetero-----------
probVirus_S_He = ModelSetup(:Virus,:Stochastic,:Hetero)
    solVirus_S_He = @time solve(probVirus_S_He,CVODE_BDF(linear_solver=:GMRES),saveat=0.1,callback=cb)

Virus_S_He = DataRangle(solVirus_S_He)
CSV.write("/home/robertgregg/R/Scripts/ABMMovie/Virus_S_He.csv", Virus_S_He)

#-----------Cell States----------
timeLength = length(solVirus_nS_Ho.t)
CellStateVirus_nS_Ho = DataFrame()
CellStateVirus_nS_Ho.x = repeat(repeat(1:N,N),timeLength)
CellStateVirus_nS_Ho.y = repeat(repeat(1:N,inner=N),timeLength)
CellStateVirus_nS_Ho.t = repeat(solVirus_S_He.t,inner=N^2)
CellStateVirus_nS_Ho.Cell = vcat([CellStateGrid(probVirus_nS_Ho.p,t) for t in solVirus_nS_Ho.t]...)
CSV.write("/home/robertgregg/R/Scripts/ABMMovie/CellStateVirus_nS_Ho.csv", CellStateVirus_nS_Ho)

CellStateVirus_S_He = DataFrame()
CellStateVirus_S_He.x = repeat(repeat(1:N,N),timeLength)
CellStateVirus_S_He.y = repeat(repeat(1:N,inner=N),timeLength)
CellStateVirus_S_He.t = repeat(solVirus_S_He.t,inner=N^2)
CellStateVirus_S_He.Cell = vcat([CellStateGrid(probVirus_S_He.p,t) for t in solVirus_S_He.t]...)
CSV.write("/home/robertgregg/R/Scripts/ABMMovie/CellStateVirus_S_He.csv", CellStateVirus_S_He)


text:
# # Mountain Runners
# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../runners.mp4" type="video/mp4">
# </video>
# ```
#
# Let's consider a race to the top of a mountain. Runners have been scattered about
# a map in some low lying areas and need to find the best path up to the peak.
#
# We'll use [`Pathfinding.AStar`](@ref) and a [`Pathfinding.PenaltyMap`](@ref) to simulate this.

# ## Setup
using Agents, Agents.Pathfinding
using Random
using FileIO # To load images you also need ImageMagick available to your project

@agent Runner GridAgent{2} begin end

# Our agent, as you can see, is very simple. Just an `id` and `pos`ition provided by
# [`@agent`](@ref). The rest of the dynamics of this example will be provided by the model.

function initialize(map_url; goal = (128, 409), seed = 88)
    ## Load an image file and convert it do a simple representation of height
    heightmap = floor.(Int, convert.(Float64, load(download(map_url))) * 255)
    ## The space of the model can be obtained directly from the image.
    ## Our example file is (400, 500).
    space = GridSpace(size(heightmap); periodic = false)
    ## The pathfinder. We use the `MaxDistance` metric since we want the runners
    ## to look for the easiest path to run, not just the most direct.
    pathfinder = AStar(space; cost_metric = PenaltyMap(heightmap, MaxDistance{2}()))
    model = ABM(
        Runner,
        space;
        rng = MersenneTwister(seed),
        properties = Dict(:goal => goal, :pathfinder => pathfinder)
    )
    for _ in 1:10
        ## Place runners in the low-lying space in the map.
        runner = add_agent!((rand(model.rng, 100:350), rand(model.rng, 50:200)), model)
        ## Everyone wants to get to the same place.
        set_target!(runner, goal, model.pathfinder)
    end
    return model
end

# The example heightmap we use here is a small region of countryside in Sweden, obtained
# with the [Tangram heightmapper](https://github.com/tangrams/heightmapper).

# ## Dynamics
# With the pathfinder in place, and all our runners having a goal position set, stepping
# is now trivial.

agent_step!(agent, model) = move_along_route!(agent, model, model.pathfinder)

# ## Let's Race
# %% #src
# Plotting is simple enough. We just need to use the [`InteractiveDynamics.abm_plot`](@ref)
# for our runners, and display the heightmap for our reference. A better interface to do
# this is currently a work in progress.
using InteractiveDynamics
using CairoMakie
CairoMakie.activate!() # hide

# We load the sample heightmap
map_url =
    "https://raw.githubusercontent.com/JuliaDynamics/" *
    "JuliaDynamics/master/videos/agents/runners_heightmap.jpg"
model = initialize(map_url)

# and plot
static_preplot!(ax, model) = scatter!(ax, model.goal; color = (:red, 50), marker = 'x')

abm_video(
    "runners.mp4",
    model,
    agent_step!;
    resolution = (700, 700),
    frames = 410,
    framerate = 45,
    ac = :black,
    as = 8,
    scatterkwargs = (strokecolor = :white, strokewidth = 2),
    heatarray = model -> penaltymap(model.pathfinder),
    heatkwargs = (colormap = :terrain,),
    static_preplot!
)
nothing # hide

# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../runners.mp4" type="video/mp4">
# </video>
# ```


text:
#=
Base.@kwdef mutable struct SIPOptimizer <: MOI.AbstractOptimizer
    objective_sense::MOI.OptimizationSense = MOI.FEASIBLE_SENSE
    time_limit::Float64 = 1000.0
    is_silent::Bool = false
    x::Vector{Float64} = Float64[]
    obj_bound::Float64 = 0.0
    obj_value::Float64 = 0.0
    termination_status::MOI.TerminationStatusCode = MOI.OPTIMIZE_NOT_CALLED
    result_status_code::MOI.ResultStatusCode = MOI.UNKNOWN_RESULT_STATUS
    run_time::Float64 = 0.0
end

function MOI.empty!(m::SIPOptimizer)
end
function MOI.is_empty(m::Optimizer)
end

function MOI.set(m::SIPOptimizer, ::MOI.Silent, value)
    m.is_silent = value
    return nothing
end

function MOI.set(m::SIPOptimizer, ::MOI.TimeLimitSec, value::Nothing)
    m.time_limit = Inf
    return nothing
end
function MOI.set(m::SIPOptimizer, ::MOI.TimeLimitSec, value::Float64)
    m.time_limit = value
    return nothing
end

MOI.set(m::SIPOptimizer, p::MOI.RawParameter, value) = error("The SIPOptimizer does not have any raw parameters.")

function MOI.get(m::SIPOptimizer, ::MOI.ListOfVariableIndices)
end
function MOI.get(m::SIPOptimizer, ::MOI.NumberOfVariables)
end

function MOI.get(m::SIPOptimizer, ::MOI.ObjectiveValue)

end
function MOI.get(m::SIPOptimizer, ::MOI.ObjectiveBound)

end

MOI.get(m::SIPOptimizer, p::MOI.RawParameter) = error("The SIPOptimizer does not have any raw parameters.")

function MOI.get(m::SIPOptimizer, ::MOI.RelativeGap)
    abs(objective_bound() - objective_value())/objective_value()
end
MOI.get(m::SIPOptimizer, ::MOI.SolverName) = "EAGO: SIPOptimizer"
MOI.get(m::SIPOptimizer, ::MOI.TerminationStatus) = m._termination_status_code
MOI.get(m::SIPOptimizer, ::MOI.PrimalStatus) = m.result_status_code
MOI.get(m::SIPOptimizer, ::MOI.SolveTime) = m.run_time
MOI.get(m::SIPOptimizer, ::MOI.ResultCount) = (m.result_status_code === MOI.FEASIBLE_POINT) ? 1 : 0
MOI.get(m::SIPOptimizer, ::MOI.TimeLimitSec) = m.time_limit

MOI.supports(::SIPOptimizer, ::MOI.TimeLimitSec) = true
MOI.supports(::SIPOptimizer, ::MOI.ObjectiveSense) = true
function MOI.set(m::SIPOptimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    m.objective_sense = sense
    return nothing
end

function MOI.optimize!(m::SIPOptimizer)
end
=#

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module AlgoSetmod

using Penalty
using DDirection
using LineSearch

"""
Liste les choix algorithmiques utilisés dans la résolution du MPCC :
* choix de la fonction de pénalité
* choix de la direction de descente
* choix de la recherche linéaire
"""

type AlgoSet
 penalty::Function
 direction::Function
 linesearch::Function

 # dans ALAS le critère d'arrêt sur le gradient du Lagrangien est donné par :
 # dual_feasible/scaling_dual<=prec
 scaling_dual::Function #dépend des multiplicateurs, la précision(r,s,t),
			#la précision demandé et le gradient du lagrangien
 unconstrained_stopping::Function
 crho_update::Function
end

#Constructeur par défaut
function AlgoSet()

  #scaling_dual=(usg,ush,uxl,uxu,ucl,ucu,lg,lh,lphi,precrst,prec,rho,dualfeas)->max(norm([usg;ush;uxl;uxu;ucl;ucu;lg;lh;lphi]),1)
  scaling_dual=(usg,ush,uxl,uxu,ucl,ucu,lg,lh,lphi,precrst,prec,rho,dualfeas)->1
  #précision du problème pénalisé "sans contraintes"
  #unconstrained_stopping=(prec,rho)->max(maximum(1./rho),prec)
  unconstrained_stopping=(prec,rho)->prec

  crho_update=(feas,rho)->min(feas*maximum(rho),100)

 return AlgoSet(Penalty.Quadratic,DDirection.NwtdirectionSpectral,LineSearch.ArmijoWolfe,scaling_dual,unconstrained_stopping,crho_update)
end

function AlgoSet(penalty::Function,direction::Function, linesearch::Function)

  scaling_dual=(usg,ush,uxl,uxu,ucl,ucu,lg,lh,lphi,precrst,prec,rho,dualfeas)->max(norm([usg;ush;uxl;uxu;ucl;ucu;lg;lh;lphi]),1)
  unconstrained_stopping=(prec,rho)->max(max.(1./rho),prec)
  crho_update=(feas,rho)->min(feas*maximum(rho),100)

 return AlgoSet(penalty,direction, linesearch,scaling_dual,unconstrained_stopping,crho_update)
end

#end of module
end


text:
# Whistle blower - complain about inappropriate content

"""
Number of whistle blower reports for a message.
- Key: `message_id`
- Value: count of whistle blower reports
"""
const WHISTLE_REPORTS_COUNT = Dict{UInt64,UInt64}()

"""
Which messages a whistle blower has reported so far.
- Key: `user_id`
- Value: Set of `message_id`
"""
const WHISTLE_REPORTS = Dict{UInt64,Set{UInt64}}()

"""
How many reports needs to be received before the message is deleted.
"""
const WHISTLE_MIN_REPORTS_FOR_REMOVAL = 3

"""
Emoji used to trigger a whistle blower report
"""
const WHISTLE_EMOJI = "🚷"

const WHISTLE_BLOWER_THANK_YOU_MESSAGE = replace_newlines("""
    **Thank you for the report.** :thumbsup:
    The message with potentially inappropriate content will be
    removed automatically once we receive enough reports.
    Please direct any question to the Moderator team.
    Your reaction emoji has also been instantly removed from
    the original message in order to keep your report anonymous
    from the public.
    """)

"""
    handler(c::Client, e::MessageReactionAdd, ::Val{:whistle})

Handle whistle blower reports. The complete flow is as follows:

1. Reporter places a special emoji on the potentially inappropriate message
2. This handler function is called.
3. Increment the counter for whistle blower reports on the reported message
4. If there are enough reports, delete the reported message
5. Otherwise, just remove the emoji.
6. Send a DM to the report about the action that's taken.
7. Log an audit event.
"""
function handler(c::Client, e::MessageReactionAdd, ::Val{:whistle})
    if e.emoji.name === WHISTLE_EMOJI
        @info "MessageReactionAdd (whistle)" e.user_id e.channel_id e.message_id e.emoji

        # Find previous reports by this user
        prior_user_reports = get!(WHISTLE_REPORTS, e.user_id, Set{UInt64}())

        # Increment counter
        count = get!(WHISTLE_REPORTS_COUNT, e.message_id, 0)
        if e.message_id ∉ prior_user_reports  # do not over-count
            count += 1
        end
        WHISTLE_REPORTS_COUNT[e.message_id] = count

        if count >= WHISTLE_MIN_REPORTS_FOR_REMOVAL
            @info "Max reports reached" e.user_id e.channel_id e.message_id
            # delete the message
            m = Message(; id=e.message_id, channel_id=e.channel_id)
            @discord delete(c, m)
        else
            # remove the emoji to protect reporter's identity
            m = Message(; id=e.message_id, channel_id=e.channel_id)
            u = User(; id=e.user_id)
            @discord delete(c, e.emoji, m, u)
        end

        # thank reporter
        dm_channel = fetchval(create_dm(c; recipient_id=e.user_id))
        uuid = string(uuid4())[(end - 11):end]
        reference = " (Reference ID: $uuid)"
        @discord create(
            c, Message, dm_channel; content=WHISTLE_BLOWER_THANK_YOU_MESSAGE * reference
        )

        # remember this report
        push!(prior_user_reports, e.message_id)

        # log an audit event
        user = @discord retrieve(c, User, e.user_id)
        audit(
            "whistle",
            user.id,
            user.username,
            user.discriminator,
            "message_id=$(e.message_id) uuid=$uuid count=$count",
        )
    end
end


text:
"""
Modia module for executing a model including code generation and calling DAE solver.

* Original developer: Toivo Henningsson, Lund
* Developer: Hilding Elmqvist, Mogram AB
* Copyright (c) 2016-2018: Hilding Elmqvist, Toivo Henningsson, Martin Otter
* License: MIT (expat)

"""
module Execution

export simulate_ida
export setOptions

using Base.Meta: quot, isexpr
using DataStructures.OrderedDict
@static if ! (VERSION < v"0.7.0-DEV.2005")
  using LinearAlgebra
  using SparseArrays
end

import ..Instantiation: Symbolic, Der, Instance, AbstractDict, VariableDict, Variable, Nothing, time_symbol, simulationModel_symbol, vars_of, check_start, GetField, This, time_global, simulationModel_global, eqs_of, get_start, get_dims, model_name_of, operator_table, prettyPrint
import ..Instantiation: Variability, constant, parameter, discrete, continuous

import ModiaMath #0.7
using Unitful
using ..ModiaLogging

const PrintJSONsolved = false

const showCode = false        # Show the code for initialization and residual calculations
const logComputations = false # Insert logging of variable values in the code
const callF = false
const showJacobian = false

global logTiming               # Show timing for each major task, simulate twice to see effect of compilation time.
global storeEliminated
global handleImpulses

function setOptions(options)
  global storeEliminated = true
  if haskey(options, :storeEliminated)
    global storeEliminated = options[:storeEliminated]
    @show storeEliminated
  end
  global handleImpulses = false
  if haskey(options, :handleImpulses)
    global handleImpulses = options[:handleImpulses]
    @show handleImpulses
  end
  global logTiming = false
  if haskey(options, :logTiming)
    global logTiming = options[:logTiming]
    @show logTiming
  end
end

# ----------------------------------

function split_variables(src::VariableDict)
    params, vars = VariableDict(), VariableDict()
    for (name, var) in src
        if !isa(var, Instance) && (!isa(var, Variable) || var.variability <= parameter)
          params[name] = var
        else
          vars[name] = var
        end
    end
    (params, vars)
end

# Substitution of variable names for vector elements, e.g. z => x[3]
const Subs = Dict{Symbolic,Any}

subs(s::Subs, ex, complete::Bool) = get(s, ex, ex)
function subs(s::Subs, ex::Symbolic, complete::Bool)
    if haskey(s, ex);  s[ex]
    elseif complete;   ModiaLogging.increaseLogCategory(:NoSubstitution); error("No substitution provided for: ", ex)
    else               ex
    end
end

function subs(s::Subs, ex::Expr, complete::Bool)
    if isexpr(ex, :quote)
        ex
    else
        Expr(ex.head, [subs(s, arg, complete) for arg in ex.args]...)
    end
end

subs(s::Subs, ex) = subs(s, ex, false)

get_value(v::Variable) = v.value
get_value(x) = x

code_reshape(ex, dims::Int...) = :( $(quot(reshape))($ex, $dims) )

code_state_read(x::Symbol, k::Int) = :($x[$k])
code_state_read(x::Symbol, k::Int, size::Int) = :($x[$k:$(k+size-1)])
code_state_read(x::Symbol, k::Int, dims::Int...) = code_reshape(code_state_read(x, k, prod(dims)), dims...)

get_name(g::GetField) = (@assert g.base == This(); g.name)

mark_solved_eq(ex, solved::Bool) = error("Only support = and := equations, got ", ex)
function mark_solved_eq(ex::Expr, solved::Bool)
    if isexpr(ex, [:(=), :(:=)], 2);  Expr(solved ? :(:=) : :(=), ex.args...)
    else                              error("Only support = and := equations, got ", ex)
    end
end
mark_solved_equations(eqs::Vector, solved::Vector{Bool}) = [mark_solved_eq(eq, sol) for (eq, sol) in zip(eqs, solved)]

find_diffstates!(d::Set{Symbolic}, ex) = nothing
find_diffstates!(d::Set{Symbolic}, der::Der) = (push!(d, der.base); nothing)
function find_diffstates!(d::Set{Symbolic}, ex::Expr)
    if !isexpr(ex, :quote)
        for arg in ex.args
            find_diffstates!(d, arg)
        end
    end
    nothing
end

function code_eliminated_func(fname, unpack, eliminated_computations, vars, x::Symbol, der_x::Symbol)
    @gensym xs der_xs results ts k
    fname = Symbol(fname)

    alloc_eliminated = []
    push_eliminated = []
    for computation in eliminated_computations
        name = (computation::Expr).args[1]::Symbol
        res_name = gensym(string(name, "_results"))

        var = vars[name]
        T, dims = var.typ, get_dims(var)
#        if T == A; T = Float64; end
#        if T != Any || T != Float64; T = Any; end  # Hack to handle units
        if typeof(T) <: Unitful.Unitlike; T = Float64; end # To handle units
#        @show T, string(var)

        if !isempty(dims);  T = Array{T,length(dims)};  end

        push!(alloc_eliminated, :($res_name = $results[$(string(name))] = Vector{$(quot(T))}(0)) )
        push!(push_eliminated,  :($(quot(push!))($res_name, $name)) )
    end
#    @show eliminated_computations
#    @show push_eliminated

    eliminated_code = quote
        function $(fname)($results, $ts, $xs, $der_xs)
            $(alloc_eliminated...)
            for ($k, $time_symbol) in enumerate($ts)
                $x = $(quot(vec))($xs[$k, :])
                $der_x = $(quot(vec))($der_xs[$k, :])
                $(unpack...)
                $(eliminated_computations...)
                $(push_eliminated...)
            end
            nothing
        end
    end
    if showCode
      @show eliminated_code
    end
    eliminated_code
end

der_name_of(name::Symbol) = Symbol("der(", name, ")")

function store_eliminated!(results::Vector{Vector}, k::Int, result)
    dest = results[k]
    if !isa(result, eltype(dest))
        T = promote_type(eltype(dest), typeof(result))
        dest = results[k] = Vector{T}(dest)
    end
    push!(dest, result)
end

function substituteExpr(ex, equations, s)
  if true # Check expression
    for eq in equations
      if isexpr(eq, :(:=), 2)
  #      lhs_name, rhs = get_name(eq.args[1]), subs(s, eq.args[2], true)
        if ex == eq.args[1]
          ex = eq.args[2]
          break
        end
      end
    end
  end
  cond = subs(s, ex, true)
  evalCond = eval(cond)
  return evalCond
end

global oldCond = false
global F_Dict = Dict()

struct Eval{F}
  f::F
end
(m::Eval)(x...) = eval(m.f)(x...)

@static if VERSION < v"0.7.0-DEV.2005"
  const letArgs = 1
else
  const letArgs = 2
end

function prepare_ida(instance::Instance, first_F_args, initial_bindings::AbstractDict{Symbol,Any}; store_eliminated = false, need_eliminated_f = false)

    proceed = zeros(Bool,1)
    global F_Dict

    params, vars = split_variables(vars_of(instance))
    for (name, var) in vars
        check_start(var, name)
    end

    s = Subs()
    for (k,(name,var)) in enumerate(vars)
        s[GetField(This(), name)] = name
        s[Der(GetField(This(), name))] = der_name_of(name)
    end
    for (name,var) in params
        s[GetField(This(), name)] = get_value(var)
    end

    s[time_global] = time_symbol
    s[simulationModel_global] = simulationModel_symbol

    computations = []
#    push!(computations, :(if ModiaMath.ModiaToModiaMath.isPreInitial($simulationModel_symbol)
#       $simulationModel_symbol.store = Synchronous.Store()
#    end))
    if logComputations
      push!(computations, :(if ! $proceed[1]; println("\nCOMPUTATIONS") end))
#      push!(computations, :(@show $simulationModel_symbol))
    end
    eliminated_computations = []

    residuals = Symbol[]
    eliminated = Symbol[]
    states = copy(vars)

    function getType(vars, v)
      vars[v].typ
    end

    modeConditions = []
    for eq in eqs_of(instance)
        if isexpr(eq, :(:=), 2)
            lhs_name, rhs = get_name(eq.args[1]), subs(s, eq.args[2], true)

            @assert haskey(states, lhs_name)
            delete!(states, lhs_name)
            T = getType(vars, lhs_name)
            ex = :($lhs_name = $rhs)
#            ex = :($lhs_name::$T = $rhs)
            push!(computations, ex)
            if logComputations
              push!(computations, :(if ! $proceed[1]; @show $lhs_name end) )
            end
            push!(eliminated_computations, ex)
            push!(eliminated, lhs_name)
        elseif isexpr(eq, :(=), 2)
            eq = subs(s, eq, true)
            lhs, rhs = eq.args
            residual = gensym("residual")
            push!(residuals, residual)
            push!(computations, :($residual = $rhs - $lhs) )
#            push!(computations, :($residual = Unitful.ustrip($rhs - $lhs)) )  # need using Unitful in function
            if logComputations
              push!(computations, :(if ! $proceed[1]; @show $residual end) )
            end
        else # conditional equations
            time = initial_bindings[time_symbol]
            s[time_global] = time
            condExpr = eq.args[1]
            evalCond = substituteExpr(condExpr, eqs_of(instance), s)
            push!(modeConditions, evalCond)
            if time == 0.0
              println("Mode condition: ", prettyPrint(condExpr), " is $evalCond at time=$time")
              oldCond = evalCond
            end

            global oldCond
            if evalCond != oldCond
              println("Mode change: ", prettyPrint(condExpr), " became ", evalCond, " at time ", time)
            end
            oldCond = evalCond

            e = if evalCond; eq.args[2] else eq.args[3] end
            e = subs(s, e, true)

            eblock = e.args

            for ei in eblock
              if ei != nothing
                lhs, rhs = ei.args
                residual = gensym("residual")
                push!(residuals, residual)
                push!(computations, :($residual = $rhs - $lhs) )
                if logComputations
                  push!(computations, :(if ! $proceed[1]; @show $residual end) )
              end
            end
          end

#            error("Unsupported equation type: ", eq)
        end
    end

    d = Set{Symbolic}()
    for eq in eqs_of(instance);  find_diffstates!(d, eq);  end

    x0 = Float64[]
    diffstates = Bool[]
    for (name,var) in states
        is_diffstate = GetField(This(), name) in d
        s = get_start(var)
        if isa(s, AbstractArray)
            append!(x0, vec(s))
#            @show name, vec(s)
            append!(diffstates, fill(is_diffstate, length(s)))
        else
            push!(x0, s)
#            push!(x0, ustrip(s))
#            @show name, s
            push!(diffstates, is_diffstate)
        end
    end
#    @show x0 size(x0)
#    @show diffstates
#=
    println("State vector allocation:")
    i = 0
    for ((name, var), offset) in zip(states, state_offsets)
        der_name = der_name_of(name)
        print("$name = x[$offset]")
        i += 1
        if diffstates[i]
          println("   $der_name = der_x[$offset]")
        else
          println()
        end
    end
=#

    # Create mapping between states and state vector
    state_sizes = [prod(get_dims(var)) for var in values(states)]
    state_offsets = cumsum(vcat([1], state_sizes[1:end-1]))

    if state_sizes != Any[]
      state_size = sum(state_sizes)
    else
      state_size = 0
    end

    if state_size > 0
      ModiaLogging.increaseLogCategory(:DynamicModel)
    else
      ModiaLogging.increaseLogCategory(:StaticModel)
    end

    loglnModia("statesize = ", state_size)

#    if ! haskey(F_Dict, modeConditions)

      @gensym x der_x r

      unpack = []

      if logComputations
        push!(unpack, :(
          if ! $proceed[1];
            println("Press enter to continue, q to stop, p to proceed: "); l = readline(STDIN); if l != "" && l[1] == 'q'; error("quit") elseif l != "" && l[1] == 'p'; $proceed[1] = true end
          end))
        push!(unpack, :(if ! $proceed[1]; println("\nUnpack:") end))
      end

  #    println("State vector allocation:")
      initials = []
      i = 1
      for ((name, var), offset) in zip(states, state_offsets)
          der_name = der_name_of(name)

          if true # diffstates[i]  # Wonder about reason???
            push!(unpack, :(global $name = $(code_state_read(x, offset, get_dims(var)...))) )
            push!(unpack, :($der_name = $(code_state_read(der_x, offset, get_dims(var)...))) )
          else
            push!(unpack, :(global $name = $(code_state_read(der_x, offset, get_dims(var)...))) )
          end
          i += prod(get_dims(var))
          if logComputations
            push!(unpack, :(if ! $proceed[1]; @show $name, $der_name end) )
          end
          push!(initials, :($name = $(quot(get_start(var)))) )
  #        push!(initials, :(@show $name) )
#          push!(initials, :($der_name = 0 * $(quot(get_start(var))) / SIUnits.Second) )
          push!(initials, :($der_name = 0 * $(quot(get_start(var))) ))
      end
      for (name,value) in initial_bindings
          push!(initials, :($name = $(quot(value))) )
      end

      initial_ex = :(let; end)
      initial_body = (initial_ex.args[letArgs].args)::Vector{Any}
      @assert isempty(initial_body)

      append!(initial_body, initials)
      append!(initial_body, instance.initial_pre)
      append!(initial_body, computations)
      append!(initial_body, instance.initial_post)
      # todo: avoid ... here too?
      push!(  initial_body, :(Any[$(residuals...)], Any[$(eliminated...)]) )

      if showCode
        println("\nINITIALIZATION CODE")
        @show initial_ex
      end
      initial_residuals, initial_eliminated = eval(initial_ex)
      residual_dims = [size(res) for res in initial_residuals]
  #    @show residual_dims
      residual_sizes = [prod(dims) for dims in residual_dims]
  #    @show residual_sizes
      if residual_sizes != Any[]
        residual_size = sum(residual_sizes)
      else
        residual_size = 0
      end
      loglnModia("residual_size = ", residual_size)
      @assert state_size == residual_size

      residual_offsets = cumsum(vcat([1], residual_sizes[1:end-1]))
  #    @show residual_offsets

      pack = []
      for (residual, dims, offset) in zip(residuals, residual_dims, residual_offsets)
          rhs = residual
          if dims === ()
              lhs = code_state_read(r, offset)
          else
              sz = prod(dims)
              lhs = code_state_read(r, offset, sz)
              if length(dims) > 1;  rhs = code_reshape(rhs, sz);  end
          end
          push!(pack, :($lhs = $rhs) )
      end

      store_elim = :nothing
      if store_eliminated
          @gensym eliminated_results
          store_elims = [
              :( $(quot(store_eliminated!))($eliminated_results, $k, $name) )
              for (k,name) in enumerate(eliminated)
          ]
          store_elim = quote
              if ($(quot(is_output_point))($simulationModel_symbol))
                  $eliminated_results = $(quot(get_eliminated_results_object))($simulationModel_symbol)
                  $(store_elims...)
              end
          end
      end

      name = Symbol("F_", model_name_of(instance))
      F_code = :( function $(name)($(first_F_args...), _t, $x, $der_x, $r, _w)
#                  $r[0] = 0.0 # For the case of no differential equations
                  end )
      F_body = (F_code.args[2].args)::Vector{Any}
  #    @assert isempty(F_body)

      append!(F_body, unpack)
      push!(F_body, :($time_symbol = $simulationModel_symbol.simulationState.time)) ####might break Sundials (not DAE)
      append!(F_body, instance.F_pre)
      append!(F_body, computations)
      append!(F_body, instance.F_post)
      append!(F_body, pack)
  #    push!(  F_body, store_elim)
      push!(  F_body, nothing)
      F_code = quote
          let
              $(F_code)
          end
      end

      if showCode
        println("\nFUNCTION F CODE")
        @show F_code
      end

#      F = Eval(F_code)
      F = eval(Module(), F_code)
#=
      if modeConditions != []
        F_Dict[modeConditions] = (F, initial_eliminated)
      end
    else
      (F, initial_eliminated) = F_Dict[modeConditions]
    end
=#

    if need_eliminated_f
        eliminated_code = code_eliminated_func(
            string("eliminated_", model_name_of(instance), "!"),
            unpack, eliminated_computations, vars, x, der_x
        )
        eliminated_f = eval(Module(), eliminated_code)
    else
        eliminated_f = nothing
    end

    eliminated_Ts = OrderedDict{Symbol,Type}()
    for (name, value) in zip(eliminated, initial_eliminated)
        eliminated_Ts[name] = typeof(value)
    end

    der_x0 = zeros(size(x0)) # todo: allow other initial guesses for der_x?

    (F, eliminated_f, x0, der_x0, diffstates, params, states, state_sizes, state_offsets, eliminated_Ts)
end

# Temporary until we can store it in SimulationModel
global_elim_results = Vector{Vector}()

# Temporary definitions
import ModiaMath.ModiaToModiaMath.ModiaSimulationModel
#is_output_point(m::ModiaSimulationModel) = true # Should not always be true!
is_output_point(m::ModiaSimulationModel) = false # Should not always be true!
get_eliminated_results_object(m::ModiaSimulationModel) = global_elim_results

function extract_results_ida(x_res, der_x_res, states, state_offsets, params)
    results = Dict{AbstractString,Any}()
#    @show keys(states)
    for ((name,var), offset) in zip(states, state_offsets)
        inds = get_dims(var) == () ? offset : (offset:(offset+prod(get_dims(var))-1))
        results[string(name)] = x_res[:,inds]
        results[string("der(",name,")")] = der_x_res[:,inds]
    end
    for (name,par) in params
        results[string(name)] = get_value(par)
    end
    results
end

# -----------------------------

function simulate_ida(instance::Instance, t, args...; kwargs...)
    simulate_ida(instance, collect(Float64, t), args...; kwargs...)
end

function simulate_ida(instance::Instance, t::Vector{Float64},
                      jac::Union{SparseMatrixCSC{Bool,Int},Nothing};# =nothing;
                      log=false, relTol=1e-4, maxSparsity=0.1,
                      store_eliminated=storeEliminated)

    if PrintJSONsolved
      printJSONforSolvedEquations(instance)
    end

    initial_bindings = Dict{Symbol,Any}(time_symbol => t[1])
    initial_m = ModiaSimulationModel()

    initial_bindings[simulationModel_symbol] = initial_m

    prep = prepare_ida(instance, [simulationModel_symbol], initial_bindings, store_eliminated=storeEliminated, need_eliminated_f=storeEliminated)
    F, eliminated_f, x0, der_x0, diffstates, params, states, state_sizes, state_offsets, eliminated = prep

    eliminated_results = Vector[Vector{T}() for T in values(eliminated)]
    # temporary until we can store it in simulationModel
    global global_elim_results = eliminated_results

    if callF || handleImpulses || showJacobian
      callResidualFunction(F, callF, handleImpulses, showJacobian, x0, der_x0, diffstates, instance)
    end

    xNames = fill("[]",length(x0))
    ii = 0
    for (name,var) in states
      ii += 1
      name = string(name)
      xNames[state_offsets[ii]] = name
      if state_sizes[ii] > 1
        dims = get_dims(var)
        dimsArray = collect(dims)
        if length(dimsArray) == 1
          for j in 1:dimsArray[1]
            xNames[state_offsets[ii]+j-1] = name*"["*string(j)*"]"
          end
        else
          xNames[state_offsets[ii]] = name*"[]"
          xNames[state_offsets[ii]+state_sizes[ii]-1] = string(collect(dimsArray))
        end
      end
    end
#    @show xNames

    start = now()

    if length(x0) > 0
      if false
        m = ModiaSimulationModel(model_name_of(instance), F, x0, der_x0, jac;
                        xw_states=diffstates, maxSparsity=maxSparsity, nc=1, nz=initial_m.nz_preInitial,
                        xNames=xNames)
      else
        m = ModiaSimulationModel(string(model_name_of(instance)), F, x0;
                        maxSparsity=maxSparsity, nc=1, nz=initial_m.nz_preInitial, jac=jac, x_fixed=diffstates)
      end
      if logTiming
        print("\n  ModiaMath:           ")
        @time ModiaMath.ModiaToModiaMath.simulate(m, t; log=log, tolRel=relTol)
        # @show now()-start
        if logTiming
          print("  ModiaMath:           ")
          @time ModiaMath.ModiaToModiaMath.simulate(m, t; log=log, tolRel=relTol)
          # @show now()-start
        end
      else
        ModiaMath.ModiaToModiaMath.simulate(m, t; log=log, tolRel=relTol)
      end
      (t_res,x_res,der_x_res) = ModiaMath.ModiaToModiaMath.getRawResult(m)
      # @show now()-start
    else
      t_res = t
      x_res = t
      der_x_res = t
    end

    results = extract_results_ida(x_res, der_x_res, states, state_offsets, params)
#    @show keys(results)
    # @show now()-start
    if store_eliminated
       Base.invokelatest(eliminated_f, results, t_res, x_res, der_x_res)
       # @show now()-start
    else
        for (name, result) in zip(keys(eliminated), eliminated_results)
            results[string(name)] = result
        end
        # @show now()-start
    end
    results["time"] = t_res
    results
end

# -----------------------------

"""
Experimental function for calculating rank of system matrix and for impulse handling
"""
function callResidualFunction(F, callF, handleImpulses, showJacobian, x0, der_x0, diffstates, instance)
    n = length(x0)
    mF = ModiaSimulationModel()
    r0 = fill(0.0, n)
    w = []

    Base.invokelatest(F, mF, 0, x0, der_x0, r0, w)
#      @show x0 der_x0 r0

    r = fill(0.0, n)
    delta = 1E-10
    A = fill(0.0, (n, n))
    independent = 0

    for i in 1:n
      X0 = copy(x0)
      X0[i] += delta
      Base.invokelatest(F, mF, 0, X0, der_x0, r, w)
      delta_r = r-r0
      A[:,i] = delta_r/delta
      if rank(A[:, 1:i]) < independent + 1
        if showJacobian
          println("Dependent column: ", i)
          if all(delta_r .== 0)
            println("Zero disturbance in residuals")
          else
            @show delta_r
            x = A[:, 1:i-1] \ delta_r
            y = [if abs(xx) < 0.1; 0.0 else xx end for xx in x]
            nonzeros = [i for i in 1:length(y) if y[i] != 0.0]
#            @show x y nonzeros
          end
        end
      else
        independent += 1
      end
    end

    der_index = [i for i in 1:length(diffstates) if diffstates[i]]
    E = fill(0.0, (n, n))
    independent = 0
    for i in 1:n
      if i in der_index
        der_X0 = copy(der_x0)
        der_X0[i] += delta
        Base.invokelatest(F, mF, 0, x0, der_X0, r, w)
        delta_r = r-r0
        E[:,i] = delta_r/delta
        if rank(E[:, 1:i]) < independent + 1
          if showJacobian
            println("Dependent column: ", i)
            if all(delta_r .== 0)
              println("Zero disturbance in residuals")
            end
          end
        else
          independent += 1
        end
      end
    end

    if showJacobian
      @show A
      @show size(A)
      @show rank(A)
      @show der_index
      @show E
      @show size(E)
      @show rank(E)
      J = A + E/delta
      @show rank(J)
      if rank(J) != size(A,1)
        println("The implicit system matrix does not have full rank.")
      end
    end

    if handleImpulses
      oldx0 = copy(x0)
      X0 = copy(x0)
      der_X0 = copy(der_x0)
      repeat = true
      while repeat
        # Make Implicit Euler step to initialize x0
        J = A + E/delta
        # J*dx + r0 = 0
        dx = J \ -r0
        X0 += dx
        der_X0 = (X0-oldx0)/delta
        Base.invokelatest(F, mF, 0, X0, der_X0, r0, w)
        repeat = norm(dx) > 1E-3  # Change to r0 after updating der_x0?
      end

      params, vars = split_variables(vars_of(instance))
      names = collect(keys(vars))
#=
      for i in 1:length(x0)
        if abs(X0[i]-oldx0[i]) > 1E5
          println("Dirac impulse in variable ", names[i])
        end
      end
=#
      for i in der_index
        if abs(X0[i]-oldx0[i]) > 1E-7
          println("Dirac impulse caused a change in the state ", names[i], " from $(oldx0[i]) to $(X0[i])")
        end
      end

      oldx0 = copy(x0)
      repeat = true
      Base.invokelatest(F, mF, 0, x0, der_x0, r0, w)
      while repeat
        # Make Implicit Euler step to initialize x0
        J = A + E/delta
        # J*dx + r0 = 0
        dx = J \ -r0
        x0 += dx
        Base.invokelatest(F, mF, 0, x0, der_x0, r0, w)
        repeat = norm(dx) > 1E-3  # Change to r0 after updating der_x0?
      end
    end
end

# -------------------------

#= causes warning
import JSON
JSON.lower(t::This) = "this."
JSON.lower(s::SIUnits.SIUnit) = string(s)
=#

"""
Experiemental code for printing the AST of solved equations in JSON format
"""
function printJSONforSolvedEquations(instance)
  modelName = string(model_name_of(instance))
  fileName = Base.Filesystem.realpath( modelName*".json" )
  println("\nWriting model in JSON format to file \"", fileName, "\"")
  indent = 3
  open(fileName,"w") do file
     variables = JSON.json(instance.variables,indent)
     print(file, "{\n")
     print(file, "\"name\":\"", modelName, "\",\n")
     print(file, "\"variables\":\n"*variables)

     jsonSolved = JSON.json(makeJSON(instance.equations),indent)
     jsonSolved = replace(jsonSolved, "\\", "")
     print(file, ",\n")
     print(file, "\"equations\":\n", jsonSolved)
     print(file, "}\n")
  end

#      jsonSolved = JSON.json(makeJSON(instance))
#      jsonSolved = replace(jsonSolved, "\\", "")
#      println(jsonSolved)
end

#makeJSON(der::Der) = Symbol("der("*string(der.base.name)*")")
makeJSON(der::Der) = Symbol("der("*replace(string(der.base.name), ".", "_")*")")
#makeJSON(get::GetField) = get.name
makeJSON(get::GetField) = Symbol(replace(string(get.name), ".", "_")) # get.name # Handle dummy derivatives
makeJSON(s::Symbol) = string(s)

#makeJSON(ex) = ex
makeJSON(ex) = get(operator_table, string(ex), string(ex))
#makeJSON(s::Instance) = s
function makeJSON(ex::Array{Any})
  [makeJSON(e) for e in ex]
end

function makeJSON(ex::Expr)
  if isexpr(ex, :quote) || isexpr(ex, :line)
    nothing
  elseif isexpr(ex, :block)
    makeJSON(ex.args[2])
  else
    Expr(ex.head, [makeJSON(arg) for arg in ex.args]...)
  end
end

end


text:
#=
Creted: 2020-05-04
Author: Emil Dmitruk

Structures for storing matrices used at different stages of preprocessing for
topological analysis with Eirene library.
=#

# TODO draw a diagram of all structures
# ===
struct MethodsParams
	plot_filters::Bool
	plot_heatmaps::Bool
	plot_betti_figrues::Bool
	lower_dist_mat_resolution::Bool
	lower_ord_mat_resolution::Bool
	legend_on::Bool
	do_dsiplay::Bool
	save_gabor_params::Bool
	save_subelements::Bool
	save_figure::Bool

	function MethodsParams(;plot_filters = false,
								plot_heatmaps = true,
								plot_betti_figrues = true,
								lower_dist_mat_resolution = false,
								lower_ord_mat_resolution = false,
								legend_on = true,
								do_dsiplay = false,
								save_gabor_params = false,
								save_subelements = false,
								save_figure = false)

		new(plot_filters, plot_heatmaps, plot_betti_figrues,
			lower_dist_mat_resolution, lower_ord_mat_resolution,
			legend_on, do_dsiplay, save_gabor_params, save_subelements,
			save_figure)
	end
end

struct ImageTopologyParams
	total_bins::Int
	max_size_limiter::Int
	min_B_dim::Int
	max_B_dim::Int
	file_name::String
	img_path::String
	gruping::Bool
	sub_img_size::Int
	sub_sample_size::Int
	pooling_method::String
	gabor_set::Int
	overlap::Number
	gaussian_blurr::Number

	function ImageTopologyParams(;total_bins = 5, max_size_limiter = 200,
								min_B_dim = 1, max_B_dim = 4,
								file_name = "", img_path="img/",
	    						gruping = true, sub_img_size = 33, sub_sample_size=2,
								pooling_method = "avg_pooling", gabor_set = 4, overlap = 0.0,
								gaussian_blurr=0.25)

		new(total_bins, max_size_limiter, min_B_dim, max_B_dim,
			file_name, img_path, gruping, sub_img_size, sub_sample_size,
			pooling_method, gabor_set, overlap, gaussian_blurr)
	end
end

function get_params_description(par::ImageTopologyParams)
	if par.pooling_method == "gauss_pooling"
		met = "_$(par.pooling_method)$(ceil(Int, par.gaussian_blurr*100))"
	else
		met = "_$(par.pooling_method)"
	end

	return "file_$(split(par.file_name,'.')[1])"*
            "_subimgsize$(par.sub_img_size)"*
            "_maxB_$(par.max_B_dim)"*
			"_minB_$(par.min_B_dim)"*
			"_gaborset_$(par.gabor_set)"*
			"_poolingsize_$(par.sub_sample_size)"*
            "$(met)_"*
			"_overlap_$(Int(par.overlap*10))_"
end

function get_ord_mat_from_img(par::ImageTopologyParams, met_par::MethodsParams; get_distances=false)
	@info "Current img_size" par.sub_img_size
	@info "Using file: " par.file_name
	@debug "Used params: " par.total_bins, par.gabor_set
	size_limiter = par.max_size_limiter

	# =============================== Get image masks ======================
	masks = get_filter_set(par.gabor_set, par.sub_img_size; plot_filters = false,
										save_gabor_params = met_par.save_gabor_params)

	# =============================== Get image ================================
	file_n = split(par.file_name, ".")[1]
	loaded_img = load(par.img_path*par.file_name)
	img1_gray = Gray.(loaded_img)
	img_size = size(img1_gray)

	# ================================ Process Image =======================
	# get the correlation matrix accordingly to choosen method
	local_correlations = get_local_correlations("gabor", img1_gray,img_size,
										par.sub_img_size; masks = masks,
										overlap=par.overlap)

	# ======================== Compute pairwise correlation ================
	dist_mat = pairwise(Euclidean(), local_correlations, dims=2)

	# =============================== Ordered matrix =======================
	size_limiter = size(dist_mat,1)
	if size_limiter > par.max_size_limiter
		@warn "Restricting matrix size, because matrix is too big"
		size_limiter = par.max_size_limiter
	end

	# ===
	# Matrix gupping
	met_par.lower_dist_mat_resolution && group_distances!(dist_mat, par.total_bins)

	# ===
	# Matrix ordering
	ordered_matrix = get_ordered_matrix(dist_mat[1:size_limiter,1:size_limiter];
												assign_same_values=par.gruping)
	if met_par.lower_ord_mat_resolution
		ordered_matrix = lower_ordmat_resolution(ordered_matrix, par.total_bins)
	end
	if get_distances
		return dist_mat
	else
		return ordered_matrix
	end
end

# ===
struct TopologyMatrixSet
    file_name::String

	# 2 below are not necessary, if params are includede in this structure
	sub_sample_size::Int
	pooling_method::String

	# Matrices
	ordered_matrix::Array
	reordered_matrix
		# reordered_map_ref
	pooled_matrix
	pooled_reord_matrix
	renum_pooled_orig_matrix
	renum_pooled_reord_matrix
	reordered_renum_pooled_orig_matrix

	matrix_collection
	description_vector
	ranks_collection

	params::ImageTopologyParams

	function TopologyMatrixSet(input_matrix::Array, params::ImageTopologyParams
								; description_vector)
		# TODO add parameter which describes which methods should be used
		# @warn "Using constant in structure definition- TODO: change to variable"
	    # file_name = images_set[6]
		# sub_sample_size = 2
		# pooling_method = "avg_pooling"
		# ===
		file_name = params.file_name

		reordered_matrix, reordered_map_ref =
	 			order_max_vals_near_diagonal2(input_matrix; do_final_plot=false, do_all_plots = false);

		pooled_matrix   	= reorganize_matrix(input_matrix; subsamp_size=params.sub_sample_size, method=params.pooling_method, gauss_sigma=params.gaussian_blurr)
		pooled_reord_matrix = reorganize_matrix(reordered_matrix; subsamp_size=params.sub_sample_size, method=params.pooling_method, gauss_sigma=params.gaussian_blurr)
		# =
		# gaussian_blurr = g_blurr
		# used_kernel = Kernel.gaussian(gaussian_blurr)
		# pooled_matrix   	= ceil.(Int,imfilter(input_matrix, used_kernel))
		# pooled_reord_matrix = ceil.(Int,imfilter(reordered_matrix, used_kernel))
		# =

		renum_pooled_orig_matrix  = get_ordered_matrix(pooled_matrix; assign_same_values=true)
		renum_pooled_reord_matrix = get_ordered_matrix(pooled_reord_matrix; assign_same_values=true)

		reordered_renum_pooled_orig_matrix, reordered_renum_pooled_orig_matrix_ref =
	 			order_max_vals_near_diagonal2(renum_pooled_orig_matrix; do_final_plot=false, do_all_plots = false);
		matrix_collection = Array[]
		push!(matrix_collection, input_matrix)
		push!(matrix_collection, reordered_matrix)
		push!(matrix_collection, pooled_matrix)
		push!(matrix_collection, pooled_reord_matrix)
		push!(matrix_collection, renum_pooled_orig_matrix)
		push!(matrix_collection, renum_pooled_reord_matrix)
		push!(matrix_collection, reordered_renum_pooled_orig_matrix)
		ranks_collection = zeros(Int,size(matrix_collection)[1])
		for mat = 1: size(matrix_collection)[1]
			ranks_collection[mat] = rank(matrix_collection[mat])
		end

		new(file_name, params.sub_sample_size,  params.pooling_method, input_matrix,
			reordered_matrix,
			# reordered_map_ref,
			pooled_matrix, pooled_reord_matrix,
			renum_pooled_orig_matrix, renum_pooled_reord_matrix,
			reordered_renum_pooled_orig_matrix,
			matrix_collection, description_vector, ranks_collection, params)
	end
end

# ===
struct TopologyMatrixBettisSet
	min_B_dim::Int
	max_B_dim::Int
	bettis_collection

	function TopologyMatrixBettisSet(top_mat_set::TopologyMatrixSet;min_B_dim=1, max_B_dim=3)
		bettis_collection = Any[]
		for matrix = top_mat_set.matrix_collection
 		   # ===
 		   # Persistent homology
 		   eirene_geom = eirene(matrix,maxdim=top_mat_set.params.max_B_dim,model="vr")
 		   bett_geom = get_bettis(eirene_geom, top_mat_set.params.max_B_dim, min_dim = top_mat_set.params.min_B_dim)
		   push!(bettis_collection,bett_geom)
	   end
	   new(min_B_dim, max_B_dim, bettis_collection)
   end
end

# ===
struct MatrixHeatmap
	heat_map
	matrix_property::String

	function MatrixHeatmap(in_array, description)
		hmap_len = size(in_array)[1]
		ordered_map1 = plot_square_heatmap(in_array, 5, hmap_len; plt_title=description,)
		new(ordered_map1, description)
	end
end

# ===
struct TopologyMatrixHeatmapsSet
	heatmaps::Array{MatrixHeatmap}

	heatmap_plots_set

	function TopologyMatrixHeatmapsSet(topology_matrix::TopologyMatrixSet)
		heatmaps = [MatrixHeatmap(topology_matrix.ordered_matrix,"original"),
					MatrixHeatmap(topology_matrix.reordered_matrix,"reordered"),
					MatrixHeatmap(topology_matrix.pooled_matrix,"pooled_origi"),
					MatrixHeatmap(topology_matrix.pooled_reord_matrix,"pooled_reordered"),
					MatrixHeatmap(topology_matrix.renum_pooled_orig_matrix,"renum_pooled_orig"),
					MatrixHeatmap(topology_matrix.renum_pooled_reord_matrix,"renum_pooled_reord"),
					MatrixHeatmap(topology_matrix.reordered_renum_pooled_orig_matrix,"reordered_renum_pooled_original"),
					]
		heatmap_plots_set = Any[]
		for hmap in heatmaps
			push!(heatmap_plots_set,hmap.heat_map)
		end
		new(heatmaps,heatmap_plots_set)
	end
end

# ===
struct BettiPlot
	betti_plot

	function BettiPlot(in_array; min_B_dim=1)
		betti_plot = plot_bettis2(in_array, "", legend_on=false, min_dim=min_B_dim);
		xlabel!("Steps");

		new(betti_plot)
	end
end

# ===
struct TopologyMatrixBettisPlots
	betti_plots_set
	function TopologyMatrixBettisPlots(bettis_collection::TopologyMatrixBettisSet)
		total_bettis = size(bettis_collection.bettis_collection)[1]
		betti_plots_set = Any[]
		for bett = 1:total_bettis
			betti_plot = BettiPlot(bettis_collection.bettis_collection[bett])
			push!(betti_plots_set, betti_plot.betti_plot)
		end

		new(betti_plots_set)
	end
end


text:
#
# This file is part of module
#   ModiaMedia (ModiaMedia/src/ModiaMedia.jl)
#
using LinearAlgebra

const excludeEnthalpyOfFormation = true                         # If true, enthalpy of formation Hf is not included in specific enthalpy hc
const referenceEnthalpy          = ReferenceEnthalpy_UserDefined   # Choice of reference enthalpy
const h_offset                   = 0.0                          # User defined offset for reference enthalpy, if ReferenceEnthalpy = ReferenceEnthalpy_UserDefined

"""
    medium = SimpleIdealGasMedium(; mediumName     = nothing,
                                    reference_p    = 101325,
                                    reference_T    = 298.15,
                                    p_default      = 101325,
                                    T_default      = 293.15,
                                    fluidConstants = nothing,
                                    data           = nothing)

Generate a `SimpleIdealGasMedium <: PureSubstance` medium object.
"""

# what is the deal with the baseProperties symbol?
struct IdealMoistAirMedium <: CondensingGases
    infos::FluidInfos
    fluidConstants::SVector{2,IdealGasFluidConstants}
    fluidLimits::FluidLimits
    data::SVector{2,SingleGasNasaData}

    function IdealMoistAirMedium(; mediumName="MoistAir",
                                    substanceNames=["water", "air"],
                                    reference_p=101325,
                                    reference_T=298.15,
                                    reference_X=[0.01, 0.99],
                                    p_default=101325,
                                    T_default=293.15,
                                    X_default=[0.01, 0.99],
                                    fluidConstants=[mediumDict["H2O"].fluidConstants[1],
                                                    mediumDict["Air"].fluidConstants[1]],
                                    data=[mediumDict["H2O"].data,
                                          mediumDict["Air"].data])

        infos = FluidInfos(mediumName           = mediumName,
                           substanceNames       = substanceNames,
                           extraPropertiesNames = fill("",0),
                           ThermoStates         = IndependentVariables_pTX,
                           baseProperties       = :BaseProperties_SingleGasNasa,
                           singleState          = false,
                           reducedX             = true,
                           fixedX               = false,
                           reference_p          = reference_p,
                           reference_T          = reference_T,
                           reference_X          = reference_X,
                           p_default            = p_default,
                           T_default            = T_default,
                           X_default            = X_default,
                           h_default            = h_pTX(SVector{2,SingleGasNasaData}(data),p_default,T_default,X_default))

        # these limits should be fixed after the specific enthalpy can be calculated below 273.15K; the high end also has problems as well.
        fluidLimits = FluidLimits(TMIN = max(mediumDict["H2O"].fluidLimits.TMIN, mediumDict["Air"].fluidLimits.TMIN),
                                  TMAX = min(mediumDict["H2O"].fluidLimits.TMAX, mediumDict["Air"].fluidLimits.TMAX),
                                  HMIN = h_pTX(SVector{2,SingleGasNasaData}(data),p_default,280.0,X_default),
                                  HMAX = h_pTX(SVector{2,SingleGasNasaData}(data),p_default,600.0,X_default))

        new(infos, SVector{2,IdealGasFluidConstants}(fluidConstants), fluidLimits, SVector{2,SingleGasNasaData}(data))
    end
end

function saturationPressure(TSat)
    # this is only valid to 273.15 K (need to add the ice saturation pressure and spline both together)
    TCritical=647.096
    pCritical=22.064e6
    r1=(1 - TSat/TCritical)

    a = [-7.85951783,1.84408259,-11.7866497,22.6807411,-15.9618719,1.80122502]
    n = [1.0,1.5,3.0,3.5,4.0,7.5]
    pSat = exp(((a[1]*r1^n[1] + a[2]*r1^n[2] + a[3]*r1^n[3] + a[4]*r1^n[4]
     + a[5]*r1^n[5] + a[6]*r1^n[6])*TCritical)/TSat)*pCritical

     return pSat
end

# enthalpy of water only for liquid condition; need to splice with the enthalpy of ice to get number valid over whole range
enthalpyOfWater(T) = 4200.0

function h_pTX(data::SVector{2,SingleGasNasaData},p,T,X)
    # steam.MM/dryAir.MM
    k_mair = data[1].MM/data[2].MM

    p_steam_sat = saturationPressure(T)
    X_sat = min(p_steam_sat*k_mair/max(100*eps(), p-p_steam_sat)*(1-X[1]), 1.0)
    X_liquid = max(X[1]-X_sat, 0.0);
    X_steam = X[1]-X_liquid;
    X_air = 1-X[1];

    h = dot([h_T(data[1], T, excludeEnthalpyOfFormation, referenceEnthalpy, (46479.819 + 2501014.5));
             h_T(data[2], T, excludeEnthalpyOfFormation, referenceEnthalpy, 25104.684)],
            [X_steam; X_air])
        + enthalpyOfWater(T)*X_liquid

    return h
end

function u_pTX(data::SVector{2,SingleGasNasaData},p,T,X)
    # steam.MM/dryAir.MM
    k_mair = data[1].MM/data[2].MM

    p_steam_sat = saturationPressure(T)
    X_sat = min(p_steam_sat*k_mair/max(100*eps(), p-p_steam_sat)*(1-X[1]), 1.0)
    X_liquid = max(X[1]-X_sat, 0.0);
    X_steam = X[1]-X_liquid;
    X_air = 1-X[1];
    R_gas = data[2].R*X_air/(1-X_liquid) + data[1].R*X_steam/(1-X_liquid);

    u = h_pTX(data,p,T,X) - R_gas*T
    return u
end

p_dTX(data::SVector{2,SingleGasNasaData},d,T,X) = d*dot([data[1].R; data[2].R], X)*T

setState_pTX(m::IdealMoistAirMedium,p,T,X) = ThermodynamicState_pTX(p,T,X)
setState_dTX(m::IdealMoistAirMedium,d,T,X) = ThermodynamicState_pTX(p_dTX(m.data,d,T,X),T,X)

pressure(               m::IdealMoistAirMedium, state::ThermodynamicState_pTX)::Float64 = state.p
temperature(            m::IdealMoistAirMedium, state::ThermodynamicState_pTX)::Float64 = state.T
gasConstant(            m::IdealMoistAirMedium, state::ThermodynamicState_pTX)::Float64 = m.data[2].R*(1-state.X[1]) + m.data[1].R*state.X[1]
density(                m::IdealMoistAirMedium, state::ThermodynamicState_pTX)::Float64 = state.p/(gasConstant(m,state)*state.T)
specificEnthalpy(       m::IdealMoistAirMedium, state::ThermodynamicState_pTX)::Float64 = h_pTX(m.data,state.p,state.T,state.X)
specificInternalEnergy( m::IdealMoistAirMedium, state::ThermodynamicState_pTX)::Float64 = u_pTX(m.data,state.p,state.T,state.X)


text:
const LinsolveSubarray = SubArray{Float64,1,Vector{Float64},Tuple{UnitRange{Int64}},true}

function admm_step!(x::Vector{Float64},
	s::SplitVector{Float64},
	μ::Vector{Float64},
	ν::LinsolveSubarray,
	x_tl::LinsolveSubarray,
	s_tl::Vector{Float64},
	ls::Vector{Float64},
	sol::Vector{Float64},
	F,
	q::Vector{Float64},
	b::Vector{Float64},
	ρ::Vector{Float64},
	α::Float64,
	σ::Float64,
	m::Int64,
	n::Int64,
	set::CompositeConvexSet{Float64})
	# linear solve
	# Create right hand side for linear system
	# deconstructed solution vector is ls = [x_tl(n+1); ν(n+1)]
	# x_tl and ν are automatically updated, since they are views on sol
	@. ls[1:n] = σ * x - q
	@. ls[(n + 1):end] = b - s + μ / ρ
	sol .= F \ ls

	# Over relaxattion
	@. x = α * x_tl + (1.0 - α) * x
	@. s_tl = s - (ν + μ) / ρ
	@. s_tl = α * s_tl + (1.0 - α) * s
	@. s = s_tl + μ / ρ

	# Project onto cone
	p_time = @elapsed project!(s, set)

	# update dual variable μ
	@. μ = μ + ρ .* (s_tl - s)
	return p_time
end

# SOLVER ROUTINE
# -------------------------------------

"""
optimize!(model)

Attempts to solve the optimization problem defined in `COSMO.Model` object with the user settings defined in `COSMO.Settings`. Returns a `COSMO.Result` object.
"""
function optimize!(ws::COSMO.Workspace)
	solver_time_start = time()
	settings = ws.settings

	# create scaling variables
	# with scaling    -> uses mutable diagonal scaling matrices
	# without scaling -> uses identity matrices
	ws.sm = (settings.scaling > 0) ? ScaleMatrices(ws.p.model_size[1], ws.p.model_size[2]) : ScaleMatrices()

	# perform preprocessing steps (scaling, initial KKT factorization)
	ws.times.setup_time = @elapsed setup!(ws);
	ws.times.proj_time  = 0. #reset projection time

	# instantiate variables
	num_iter = 0
	status = :Unsolved
	cost = Inf
	r_prim = Inf
	r_dual = Inf

	# print information about settings to the screen
	settings.verbose && print_header(ws)
	time_limit_start = time()

	#preallocate arrays
	m, n = ws.p.model_size
	δx = zeros(n)
	δy =  zeros(m)
	s_tl = zeros(m) # i.e. sTilde

	ls = zeros(n + m)
	sol = zeros(n + m)
	x_tl = view(sol, 1:n) # i.e. xTilde
	ν = view(sol, (n + 1):(n + m))

	settings.verbose_timing && (iter_start = time())

	for iter = 1:settings.max_iter

		num_iter+= 1

		@. δx = ws.vars.x
		@. δy = ws.vars.μ

		ws.times.proj_time += admm_step!(
			ws.vars.x, ws.vars.s, ws.vars.μ, ν,
			x_tl, s_tl, ls, sol,
			ws.F, ws.p.q, ws.p.b, ws.ρvec,
			settings.alpha, settings.sigma,
			m, n, ws.p.C);

		# compute deltas for infeasibility detection
		@. δx = ws.vars.x - δx
		@. δy = -ws.vars.μ + δy

		# compute residuals (based on optimality conditions of the problem) to check for termination condition
		# compute them every {settings.check_termination} step
		mod(iter, settings.check_termination)  == 0 && ((r_prim, r_dual) = calculate_residuals(ws))

		# check convergence with residuals every {settings.checkIteration} steps
		if mod(iter, settings.check_termination) == 0
			# update cost
			cost = ws.sm.cinv[] * (1/2 * ws.vars.x' * ws.p.P * ws.vars.x + ws.p.q' * ws.vars.x)[1]

			if abs(cost) > 1e20
				status = :Unsolved
				break
			end

			# print iteration steps
			settings.verbose && print_iteration(settings, iter, cost, r_prim, r_dual)

			if has_converged(ws, r_prim, r_dual)
				status = :Solved
				break
			end
		end

		# check infeasibility conditions every {settings.checkInfeasibility} steps
		if mod(iter, settings.check_infeasibility) == 0
			if is_primal_infeasible(δy, ws)
				status = :Primal_infeasible
				cost = Inf
				break
			end

			if is_dual_infeasible(δx, ws)
				status = :Dual_infeasible
				cost = -Inf
				break
			end
		end

		# adapt rhoVec if enabled
		if settings.adaptive_rho && (mod(iter, settings.adaptive_rho_interval) == 0) && (settings.adaptive_rho_interval > 0)
			adapt_rho_vec!(ws)
		end

		if settings.time_limit !=0 &&  (time() - time_limit_start) > settings.time_limit
			status = :Time_limit_reached
			break
		end

	end #END-ADMM-MAIN-LOOP

	settings.verbose_timing && (ws.times.iter_time = (time() - iter_start))
	settings.verbose_timing && (ws.times.post_time = time())

	# calculate primal and dual residuals
	if num_iter == settings.max_iter
		r_prim, r_dual = calculate_residuals(ws)
		status = :Max_iter_reached
	end

	# reverse scaling for scaled feasible cases
	if settings.scaling != 0
		reverse_scaling!(ws)
		# FIXME: Another cost calculation is not necessary since cost value is not affected by scaling
		cost =  (1/2 * ws.vars.x' * ws.p.P * ws.vars.x + ws.p.q' * ws.vars.x)[1] #sm.cinv * not necessary anymore since reverseScaling
	end

	ws.times.solver_time = time() - solver_time_start
	settings.verbose_timing && (ws.times.post_time = time() - ws.times.post_time)
	# print solution to screen
	settings.verbose && print_result(status, num_iter, cost, ws.times.solver_time)

	# create result object
	res_info = ResultInfo(r_prim, r_dual)
	y = -ws.vars.μ
	return result = Result{Float64}(ws.vars.x, y, ws.vars.s.data, cost, num_iter, status, res_info, ws.times);

end


text:
using HierarchicalBetaBinomial
using Test

@testset "HierarchicalBetaBinomial.jl" begin
    # Write your tests here.
end


text:
"""
    Schur complement structure of the form:
       M = A B
           C D
       |A B| |x| = |u|
       |C D| |y| = |v|
    M ∈ R{n+m,n+m}
    A ∈ R{n,n}
    D ∈ R{m,m}
    The Schur complement is taken about A, the top-left block:
        SC = D - C * A^-1 * B
"""
mutable struct Schur{T,n,m,nn,nm,mm}
    A::SMatrix{n,n,T,nn}
    B::SMatrix{n,m,T,nm}
    C::SMatrix{m,n,T,nm}
    D::SMatrix{m,m,T,mm}
    Ai::SMatrix{n,n,T,nn} # Inverse of A
    CAi::SMatrix{m,n,T,nm} # C*A^{-1}
    CAiB::SMatrix{m,m,T,mm} # C*A^{-1}*B
    gs_data::SDMGSSolver{m,T} # Schur complement of the block A in M
    u::SVector{n,T}
    v::SVector{m,T}
    x::SVector{n,T}
    y::SVector{m,T}
end

"""
    Initialize a Schur Object from a matrix M. The matrix M is split in 4 blocks:
        - top left block = n x n
        - bottom right block = m x m
"""
function Schur(M::AbstractMatrix{T}; n::Int=1, m::Int=size(M)[1]-n) where {T}
    @assert all(size(M) .== n+m)
    A = SMatrix{n,n,T,n^2}(M[1:n,      1:n])
    B = SMatrix{n,m,T,n*m}(M[1:n,      n .+ (1:m)])
    C = SMatrix{m,n,T,n*m}(M[n .+ (1:m),1:n])
    D = SMatrix{m,m,T,m^2}(M[n .+ (1:m),n .+ (1:m)])
    Ai = inv(A)
    CAi = C*Ai
    CAiB = C*Ai*B
    gs_data = SDMGSSolver(m,T=T)
    factorize!(gs_data, D - CAiB)
    u = zeros(SVector{n,T})
    v = zeros(SVector{m,T})
    x = zeros(SVector{n,T})
    y = zeros(SVector{m,T})
    return Schur{T,n,m,n^2,n*m,m^2}(A,B,C,D,Ai,CAi,CAiB,gs_data,u,v,x,y)
end

"""
    Update the Schur complement structure to handle a change in the D matrix.
    It recomputes the Gram-Schmidt factorization of the Schur complement of the block A.
"""
function schur_factorize!(S::Schur{T,n,m,nn,nm,mm}, D::AbstractMatrix{T}) where {T,n,m,nn,nm,mm}
    # update D, gs_data and the shur complement
    # we may want to set S.D to D so that it's clean, even if it's not necessary
    # S.D = D
    CAiB = S.CAiB
    gs_data = S.gs_data
    factorize!(gs_data, D - CAiB)
    return nothing
end

"""
    Solve for [x;y] given [u,v] and the factorized matrix M.
"""
function schur_solve!(S::Schur{T,n,m,nn,nm,mm}, u::AbstractVector{T}, v::AbstractVector{T}) where {T,n,m,nn,nm,mm}
    B  = S.B
    Ai = S.Ai
    CAi = S.CAi
    gs_data = S.gs_data
    S.u = u
    S.v = v
    us = S.u # staticarrays
    vs = S.v # staticarrays
    # S.x = Ai*(us + B*(As*(C*(Ai*us) - vs)))
    # S.y = As*(- C*(Ai*us) + vs)

    qr_solve!(gs_data, CAi*us - vs)
    temp = gs_data.xs
    S.x = Ai*(us + B*temp)
    S.y = - temp
    return nothing
end


text:
using StatisticalRethinking, CmdStan
#gr(size=(600,800));

ProjDir = rel_path("..", "scripts", "03")
cd(ProjDir)

binomialstanmodel = "
// Inferring a Rate
data {
  int N;
  int<lower=0> k[N];
  int<lower=1> n[N];
}
parameters {
  real<lower=0,upper=1> theta;
  real<lower=0,upper=1> thetaprior;
}
model {
  // Prior Distribution for Rate Theta
  theta ~ beta(1, 1);
  thetaprior ~ beta(1, 1);

  // Observed Counts
  k ~ binomial(n, theta);
}
";

stanmodel = Stanmodel(name="binomial", monitors = ["theta"], model=binomialstanmodel,
  output_format=:mcmcchains);

N2 = 4^2
d = Binomial(9, 0.66)
n2 = Int.(9 * ones(Int, N2))
k2 = rand(d, N2);

binomialdata = Dict("N" => length(n2), "n" => n2, "k" => k2);

rc, chn, cnames = stan(stanmodel, binomialdata, ProjDir, diagnostics=false,
  CmdStanDir=CMDSTAN_HOME);

MCMCChains.describe(chn)

if rc == 0
  plot(chn)
end

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
import BlockArrays: BlockSizes, BlockIndex, globalrange, nblocks, global2blockindex, blockindex2global

#=
[1,1  1,2] | [1,3  1,4  1,5]
--------------------------
[2,1  2,2] | [2,3  2,4  2,5]
[3,1  3,2] | [3,3  3,4  3,5]
----------------------------
[4,1  4,2] | [4,3  4,4  4,5]
[5,1  5,2] | [5,3  5,4  5,5]
[6,1  6,2] | [6,3  6,4  6,5]
=#

block_size = BlockSizes([1,2,3], [2, 3])

@test nblocks(block_size) == (3,2)
@test nblocks(block_size, 1) == 3
@test nblocks(block_size, 2) == 2

@inferred globalrange(block_size, (1,1)) == (1:1, 1:2)
@inferred globalrange(block_size, (1,2)) == (1:1, 3:5)
@inferred globalrange(block_size, (2,1)) == (2:3, 1:2)
@inferred globalrange(block_size, (2,2)) == (2:3, 3:5)

@inferred global2blockindex(block_size, (3, 1)) == BlockIndex((2,1), (2,1))
@inferred global2blockindex(block_size, (1, 4)) == BlockIndex((1,2), (1,2))
@inferred global2blockindex(block_size, (4, 5)) == BlockIndex((3,2), (1,3))

@inferred blockindex2global(block_size, BlockIndex((2,1), (2,1))) == (3, 1)
@inferred blockindex2global(block_size, BlockIndex((1,2), (1,2))) == (1, 4)
@inferred blockindex2global(block_size, BlockIndex((3,2), (1,3))) == (4, 5)

buf = IOBuffer()
print(buf, block_size)
@test takebuf_string(buf) == "[1,2,3]×[2,3]"


text:
using Test
using QuantumLattices.Prerequisites.Combinatorics

@testset "Combinations" begin
    @test Combinations{0}("abc")|>collect == [()]
    @test Combinations{1}("abc")|>collect == [('a',), ('b',), ('c',)]
    @test Combinations{2}("abc")|>collect == [('a', 'b'), ('a', 'c'), ('b', 'c')]
    @test Combinations{3}("abc")|>collect == [('a', 'b', 'c')]
end

@testset "DulCombinations" begin
    @test DulCombinations{0}("abc")|>collect == [()]
    @test DulCombinations{1}("abc")|>collect == [('a',), ('b',), ('c',)]
    @test DulCombinations{2}("abc")|>collect == [('a', 'a'), ('a', 'b'), ('a', 'c'), ('b', 'b'), ('b', 'c'), ('c', 'c')]
    @test DulCombinations{3}("abc")|>collect == [('a', 'a', 'a'), ('a', 'a', 'b'), ('a', 'a', 'c'), ('a', 'b', 'b'), ('a', 'b', 'c'), ('a', 'c', 'c'), ('b', 'b', 'b'), ('b', 'b', 'c'), ('b', 'c', 'c'), ('c', 'c', 'c')]
end

@testset "Permutations" begin
    @test Permutations{0}("abc")|>collect == [()]
    @test Permutations{1}("abc")|>collect == [('a',), ('b',), ('c',)]
    @test Permutations{2}("abc")|>collect == [('a', 'b'), ('a', 'c'), ('b', 'a'), ('b', 'c'), ('c', 'a'), ('c', 'b')]
    @test Permutations{3}("abc")|>collect == [('a', 'b', 'c'), ('a', 'c', 'b'), ('b', 'a', 'c'), ('b', 'c', 'a'), ('c', 'a', 'b'), ('c', 'b', 'a')]
end

@testset "DulPermutations" begin
    @test DulPermutations{0}("abc")|>collect == [()]
    @test DulPermutations{1}("abc")|>collect == [('a',), ('b',), ('c',)]
    @test DulPermutations{2}("abc")|>collect == [('a', 'a'), ('a', 'b'), ('a', 'c'), ('b', 'a'), ('b', 'b'), ('b', 'c'), ('c', 'a'), ('c', 'b'), ('c', 'c')]
    @test DulPermutations{3}("abc")|>collect == [(i, j, k) for i∈"abc" for j∈"abc" for k∈"abc"]
end


text:
function extract_scope_wfm(scope, args...)
data = []
for ch in args
	x = (ch1, initial_x, x_increm) = GI.SCOPE.fetch_wfm(scope,ch)
	push!(data,x)
end
return data
# can destructure this array automatically when calling the function
# ch1,ch2,ch3,ch4 = extract_scope_wfm(scope1, 1,2,3,4)
end


text:
# Calcium buffering dynamics
include("common.jl")

# Calcium buffering parameters
@with_kw struct CABUFParams
    KM_CA
    ET
end

_β_ca(ca, KM, ET) = _mm((ca + KM)^2, KM * ET)
_β_ca(ca, p::CABUFParams) = _β_ca(ca, p.KM_CA, p.ET)


text:
# IntervalTree map operations

typealias IntervalMap{K, V} IntervalTree{K, IntervalValue{K, V}}

function Base.setindex!{K, V, B}(t::IntervalBTree{K, V, B}, value,
                                 key::Tuple{Any, Any})
    push!(t, V(key[1], key[2], value), true)
end

function Base.setindex!{K, V, B}(t::IntervalBTree{K, V, B}, value,
                                 key::Interval{K})
    push!(t, V(key.first, key.last, value), true)
end

function Base.setindex!{K, V, B}(t::IntervalBTree{K, V, B}, value,
                                 first, last)
    push!(t, V(first, last, value), true)
end

function Base.getindex{K, V, B}(t::IntervalBTree{K, V, B}, key::AbstractInterval{K})
    return _getindex(t.root, key)
end

function _getindex{K, V, B}(t::InternalNode{K, V, B}, key::AbstractInterval{K})
    i = findidx(t, key)
    if 1 <= length(t) - 1 && key >= t.keys[i]
        return _getindex(t.children[i+1], key)
    else
        return _getindex(t.children[i], key)
    end
end

function _getindex{K, V, B}(t::LeafNode{K, V, B}, key::AbstractInterval{K})
    i = findidx(t, key)
    if 1 <= i <= length(t) && first(t.entries[i]) == first(key) &&
        last(t.entries[i]) == last(key)
        return t.values[i]
    else
        error(KeyError((key.a, key.b)))
    end
end

function Base.get{K, V, B}(t::IntervalBTree{K, V, B}, key::Tuple{Any, Any}, default)
    return _get(t.root, Interval{K}(key[1], key[2]), default)
end

function Base.get{K, V, B}(t::IntervalBTree{K, V, B}, key::AbstractInterval{K}, default)
    return _get(t.root, key, default)
end

function _get{K, V, B}(t::InternalNode{K, V, B}, key::AbstractInterval{K}, default)
    i = findidx(t, key)
    if 1 <= length(t) - 1 && key >= t.keys[i]
        return _get(t.children[i+1], key, default)
    else
        return _get(t.children[i], key, default)
    end
end

function _get{K, V, B}(t::LeafNode{K, V, B}, key::AbstractInterval{K}, default)
    i = findidx(t, key)
    if 1 <= i <= length(t) &&
        first(t.entries[i]) == first(key) && last(t.entries[i]) == last(key)
        return t.entries[i]
    else
        return default
    end
end

function Base.get!{K, V, B}(t::IntervalBTree{K, V, B}, key::Tuple{Any, Any}, default)
    return push!(t, V(key[1], key[2], default), true, false)
end

function Base.get!{K, V, B}(t::IntervalBTree{K, V, B}, key::AbstractInterval{K}, default)
    return push!(t, V(key.first, key.last, default), true, false)
end

function Base.delete!{K, V, B}(t::IntervalBTree{K, V, B}, first, last)
    return deletefirst!(t, first, last)
end

function Base.delete!{K, V, B}(t::IntervalBTree{K, V, B}, key::Tuple{Any, Any})
    return deletefirst!(t, key)
end

function Base.delete!{K, V, B}(t::IntervalBTree{K, V, B}, key::Interval{K})
    return deletefirst!(t, key)
end

immutable IntervalKeyIterator{K, V, B}
    t::IntervalBTree{K, V, B}
end

function Base.keys(t::IntervalBTree)
    return IntervalKeyIterator(t)
end

function Base.start{K, V, B}(it::IntervalKeyIterator{K, V, B})
    # traverse to the first leaf node
    node = it.t.root
    while !isa(node, LeafNode{K, V, B})
        node = node.children[1]
    end

    return IntervalBTreeIteratorState(Nullable(node), 1)
end

function Base.next{K, V, B}(t::IntervalKeyIterator{K, V, B},
                            state::IntervalBTreeIteratorState{K, V, B})
    leaf = get(state.leaf)
    key = Interval{K}(first(leaf.entries[state.i]),
                      last(leaf.entries[state.i]))
    if state.i < length(leaf)
        state = IntervalBTreeIteratorState{K, V, B}(leaf, state.i + 1)
    else
        state = IntervalBTreeIteratorState{K, V, B}(leaf.right, 1)
    end
    return key, state
end

function Base.done{K, V, B}(t::IntervalKeyIterator{K, V, B},
                            state::IntervalBTreeIteratorState{K, V, B})
    return isnull(state.leaf) || isempty(get(state.leaf))
end

immutable IntervalValueIterator{K, V <: IntervalValue, B}
    t::IntervalBTree{K, V, B}
end

function Base.values{K, V <: IntervalValue, B}(t::IntervalBTree{K, V, B})
    return IntervalValueIterator(t)
end

function Base.start{K, V, B}(it::IntervalValueIterator{K, V, B})
    # traverse to the first leaf node
    node = it.t.root
    while !isa(node, LeafNode{K, V, B})
        node = node.children[1]
    end

    return IntervalBTreeIteratorState(Nullable(node), 1)
end

function Base.next{K, V, B}(t::IntervalValueIterator{K, V, B},
                            state::IntervalBTreeIteratorState{K, V, B})
    leaf = get(state.leaf)
    value = leaf.entries[state.i].value
    if state.i < length(leaf)
        state = IntervalBTreeIteratorState{K, V, B}(leaf, state.i + 1)
    else
        state = IntervalBTreeIteratorState{K, V, B}(leaf.right, 1)
    end
    return value, state
end

function Base.done{K, V, B}(t::IntervalValueIterator{K, V, B},
                            state::IntervalBTreeIteratorState{K, V, B})
    return isnull(state.leaf) || isempty(get(state.leaf))
end


text:
module SBTestSimpleSDMLayers

using SimpleSDMLayers
using SpatialBoundaries
using Test

precipitation = SimpleSDMPredictor(WorldClim, BioClim, 12; left=-80.0, right=-56.0, bottom=44.0, top=62.0)

W = wombling(precipitation)

@test isa(W, LatticeWomble)

Lt, Ld = SimpleSDMPredictor(W)
@test isa(Lt, SimpleSDMPredictor)

Rt, Rd = SimpleSDMResponse(W)
@test isa(Rt, SimpleSDMResponse)

@test !any(map(isnan, Rt.grid))
@test !any(map(isnan, Rd.grid))

end


text:
using PlotlyJS
using MAT

## ============================== DEP ================================
file = matread("/Users/martingurtner/Downloads/DEP_conv_larger_larger.mat")
t_elapsed_admm_mean = file["t_elapsed_admm_mean"]
t_elapsed_admm_std = file["t_elapsed_admm_std"]
t_elapsed_centralized_mean = file["t_elapsed_centralized_mean"]
t_elapsed_centralized_std = file["t_elapsed_centralized_std"]
N_agnts_vals = file["N_agnts_vals"]
actuators_per_agent_vals = file["actuators_per_agent_vals"]

##
x_grid = [x for x = N_agnts_vals for y = actuators_per_agent_vals]
y_grid = [y for x = N_agnts_vals for y = actuators_per_agent_vals]
data_admm = vec(transpose(t_elapsed_admm_mean))
data_centralzied = vec(transpose(t_elapsed_centralized_mean))
##
data = (t_elapsed_admm_mean - t_elapsed_centralized_mean)./t_elapsed_admm_mean
I_1 = 4
I_2 = 1
data = data[I_1:end,I_2:end]
x = N_agnts_vals[I_1:end]
y = actuators_per_agent_vals[I_2:end]
##

titlefontsize = 30;
labelfontsize = 25;
ticksfontsize = 20;
figsize_paramConv = (700,600);
figsize_exps = (figsize_paramConv[2],200);

cmax = maximum(data)
cmin = minimum(data)

c_zero = -cmin/(cmax - cmin)

plot_data = contour(;
    x=x,
    y=y,
    z=data,
    colorscale=[
        [0, :red],
        [c_zero, :white],
        [1, :blue]
        ],
    colorbar=attr(;
        title="(t₁-t₂)/t₁",titleside="right",
        tickfont=attr(;size=ticksfontsize),
        titlefont=attr(;size=labelfontsize),
    ),
    zauto=false,
    zmin=cmin,
    zmax=cmax,
)

p_DEP = Plot(plot_data,
    Layout(;
        # title="DEP",
        # titlefont=attr(;size=titlefontsize),
        width = figsize_paramConv[1],
        height = figsize_paramConv[2],
        margin=attr(;t=0,l=100,b=0,r=100),
        xaxis=attr(title="Number of agents",
                tickfont=attr(;size=ticksfontsize),
                titlefont=attr(;size=labelfontsize),
                automargin=true),
        yaxis=attr(title="Number of actuators per agent",
               tickfont=attr(;size=ticksfontsize),
                titlefont=attr(;size=labelfontsize),
                automargin=true),
        )
     )
savefig(p_DEP, "/Users/martingurtner/Downloads/DEP_distVsCentr.pdf");
## ============================== MAG ================================

file = matread("/Users/martingurtner/Downloads/MAG_conv_trim_M1.mat")
t_elapsed_admm_mean = file["t_elapsed_admm_mean"]
t_elapsed_admm_std = file["t_elapsed_admm_std"]
t_elapsed_centralized_mean = file["t_elapsed_centralized_mean"]
t_elapsed_centralized_std = file["t_elapsed_centralized_std"]
N_agnts_vals = file["N_agnts_vals"]
actuators_per_agent_vals = file["actuators_per_agent_vals"]

##
x_grid = [x for x = N_agnts_vals for y = actuators_per_agent_vals]
y_grid = [y for x = N_agnts_vals for y = actuators_per_agent_vals]
data_admm = vec(transpose(t_elapsed_admm_mean))
data_centralzied = vec(transpose(t_elapsed_centralized_mean))
##
data = (t_elapsed_admm_mean - t_elapsed_centralized_mean)./t_elapsed_admm_mean
I_1 = 1
I_1end = 7
I_2 = 1
data = data[I_1:I_1end,I_2:end]
x = N_agnts_vals[I_1:I_1end]
y = actuators_per_agent_vals[I_2:end]
##

titlefontsize = 30;
labelfontsize = 25;
ticksfontsize = 20;
figsize_paramConv = (700,600);
figsize_exps = (figsize_paramConv[2],200);

cmax = maximum(data)
cmin = minimum(data)

c_zero = -cmin/(cmax - cmin)

plot_data = contour(;
    x=x,
    y=y,
    z=data,
    colorscale=[
        [0, :red],
        [c_zero, :white],
        [1, :blue]
        ],
    colorbar=attr(;
        title="(t₁-t₂)/t₁",titleside="right",
        tickfont=attr(;size=ticksfontsize),
        titlefont=attr(;size=labelfontsize),
    ),
    zauto=false,
    zmin=cmin,
    zmax=cmax,
)

p_MAG = Plot(plot_data,
    Layout(;
        # title="DEP",
        # titlefont=attr(;size=titlefontsize),
        width = figsize_paramConv[1],
        height = figsize_paramConv[2],
        margin=attr(;t=0,l=100,b=0,r=100),
        xaxis=attr(title="Number of agents",
                tickfont=attr(;size=ticksfontsize),
                titlefont=attr(;size=labelfontsize),
                automargin=true),
        yaxis=attr(title="Number of actuators per agent",
               tickfont=attr(;size=ticksfontsize),
                titlefont=attr(;size=labelfontsize),
                automargin=true),
        )
     )
savefig(p_MAG, "/Users/martingurtner/Downloads/MAG_distVsCentr.pdf");

## ============================== ACU ================================

file = matread("/Users/martingurtner/Downloads/ACU_conv_trim_3.mat")
# t_elapsed_admm = file["t_elapsed_admm"]
t_elapsed_admm_mean = file["t_elapsed_admm_mean"]
t_elapsed_admm_std = file["t_elapsed_admm_std"]
# t_elapsed_centralized = file["t_elapsed_centralized"]
t_elapsed_centralized_mean = file["t_elapsed_centralized_mean"]
t_elapsed_centralized_std = file["t_elapsed_centralized_std"]
N_agnts_vals = file["N_agnts_vals"]
actuators_per_agent_vals = file["actuators_per_agent_vals"]

# data_mean = data -> [mean(trim(data[i,j,:], prop=0.2)) for i in 1:size(data)[1], j in 1:size(data)[2]]
# t_elapsed_admm_mean = data_mean(t_elapsed_admm)
# t_elapsed_centralized_mean = data_mean(t_elapsed_centralized)

##
x_grid = [x for x = N_agnts_vals for y = actuators_per_agent_vals]
y_grid = [y for x = N_agnts_vals for y = actuators_per_agent_vals]
data_admm = vec(transpose(t_elapsed_admm_mean))
data_centralzied = vec(transpose(t_elapsed_centralized_mean))
##
data = (t_elapsed_admm_mean - t_elapsed_centralized_mean)./t_elapsed_centralized_mean
# I_1 = 3
# I_1_step = 2
# I_1_end = 23
# I_2 = 1
# I_2_end = 9

I_1 = 1
I_1_step = 1
I_1_end = length(N_agnts_vals)
I_2 = 1
I_2_end = length(actuators_per_agent_vals)

data = data[I_1:I_1_step:I_1_end,I_2:I_2_end];
x = N_agnts_vals[I_1:I_1_step:I_1_end];
y = actuators_per_agent_vals[I_2:I_2_end];
##

titlefontsize = 30;
labelfontsize = 25;
ticksfontsize = 20;
figsize_paramConv = (700,600);
figsize_exps = (figsize_paramConv[2],200);

cmax = maximum(data)
cmin = minimum(data)

c_zero = -cmin/(cmax - cmin)

plot_data = contour(;
    x=x,
    y=y,
    z=data,
    line_smoothing=1,
    colorscale=[
        [0, :red],
        [c_zero, :white],
        [1, :blue]
        ],
    colorbar=attr(;
        title="(t₁-t₂)/t₁",titleside="right",
        tickfont=attr(;size=ticksfontsize),
        titlefont=attr(;size=labelfontsize),
    ),
    zauto=false,
    zmin=cmin,
    zmax=cmax,
)

p_ACU = Plot(plot_data,
    Layout(;
        # title="DEP",
        # titlefont=attr(;size=titlefontsize),
        width = figsize_paramConv[1],
        height = figsize_paramConv[2],
        margin=attr(;t=0,l=100,b=0,r=100),
        xaxis=attr(title="Number of agents",
                tickfont=attr(;size=ticksfontsize),
                titlefont=attr(;size=labelfontsize),
                automargin=true),
        yaxis=attr(title="Number of actuators per agent",
               tickfont=attr(;size=ticksfontsize),
                titlefont=attr(;size=labelfontsize),
                automargin=true),
        )
     )

savefig(p_ACU, "/Users/martingurtner/Downloads/ACU_distVsCentr.pdf");

##
# using Plots
# using MAT
# gr()

# file = matread("/Users/martingurtner/Downloads/ACU_conv_M1.mat")
# t_elapsed_admm_mean = file["t_elapsed_admm_mean"]
# t_elapsed_admm_std = file["t_elapsed_admm_std"]
# t_elapsed_centralized_mean = file["t_elapsed_centralized_mean"]
# t_elapsed_centralized_std = file["t_elapsed_centralized_std"]
# N_agnts_vals = file["N_agnts_vals"]
# actuators_per_agent_vals = file["actuators_per_agent_vals"]
# ##
# # plot(N_agnts_vals, [t_elapsed_admm_mean[:,9], t_elapsed_centralized_mean[:,9]])
# # plot!(N_agnts_vals, t_elapsed_centralized_mean[:,9])
# plot(N_agnts_vals, t_elapsed_admm_mean,
# ribbon=(t_elapsed_admm_std,t_elapsed_admm_std),
# xlabel="Number of agents", ylabel="Time [s]", label="ADMM")
# plot!(N_agnts_vals, t_elapsed_centralized_mean,
#     ribbon=(t_elapsed_centralized_std,t_elapsed_centralized_std),
#     label="Centralized")


text:
###############################################################################
#
# nanmean, nanmedian, nanstd, nanmax, nanmin
#
###############################################################################
"""

    nanmax(x::Array)

Return the maximum of array ommiting the NaN, given
- `x` Array of numbers, can be NaN

---
Example
```julia
xs = [1, 2, 4, NaN];
nmax = nanmax(xs);
```
"""
function nanmax(x::Array)
    _x = filter(!isnan, x);

    if length(_x) == 0 return NaN end;

    return maximum( _x )
end

"""

    nanmean(x::Array)

Return the mean of array by ommiting the NaN, given
- `x` Array of numbers, can be NaN

---
Example
```julia
xs = [1, 2, 4, NaN];
nmean = nanmean(xs);
```
"""
function nanmean(x::Array)
    _x = filter(!isnan, x);

    if length(_x) == 0 return NaN end;

    return mean( _x )
end

"""

    nanmedian(x::Array)

Return the median of array by ommiting the NaN, given
- `x` Array of numbers, can be NaN

---
Example
```julia
xs = [1, 2, 4, NaN];
nmed = nanmedian(xs);
```
"""
function nanmedian(x::Array)
    _x = filter(!isnan, x);

    if length(_x) == 0 return NaN end;

    return median( _x )
end

"""

    nanmin(x::Array)

Return the maximum of array ommiting the NaN, given
- `x` Array of numbers, can be NaN

---
Example
```julia
xs = [1, 2, 4, NaN];
nmin = nanmin(xs);
```
"""
function nanmin(x::Array)
    _x = filter(!isnan, x);

    if length(_x) == 0 return NaN end;

    return minimum( _x )
end

"""

    nanpercentile(x::Array, p::Number)

Return the percentile by excluding the NaN of given
- `x` Array of data
- `p` Percentile

---
Example
```julia
xs = rand(100);
pth = nanpercentile(rand(100), 50);
xs[1:10] .= NaN;
pth = nanpercentile(rand(100), 50);
```
"""
function nanpercentile(x::Array, p::Number)
    @assert 0 <= p <= 100

    _x = filter(!isnan, x);

    if length(_x) == 0 return NaN end;

    return percentile( _x, p )
end

"""

    nanstd(x::Array)

Return the std of array by ommiting the NaN, given
- `x` Array of numbers, can be NaN

---
Example
```julia
xs = [1, 2, 4, NaN];
nstd = nanstd(xs);
```
"""
function nanstd(x::Array)
    _x = filter(!isnan, x);

    if length(_x) == 0 return NaN end;

    return std( _x )
end

###############################################################################
#
# Error measures
#
###############################################################################
"""

    mae(y::Array, pred::Array)

Return the mean absolute error by ommiting the NaN, given
- `y` Array of numbers, can be NaN
- `pred` Array of predictions, can be NaN

---
Example
```julia
ys = [1, 2, 4, NaN];
preds = [2, 4, 1, NaN];
nmae = mae(ys, preds);
```
"""
function mae(y::Array, pred::Array)
    return nanmean( abs.(y .- pred) )
end

"""

    mape(y::Array, pred::Array)

Return the mean absolute percentage error by ommiting the NaN, given
- `y` Array of numbers, can be NaN
- `pred` Array of predictions, can be NaN

---
Example
```julia
ys = [1, 2, 4, NaN];
preds = [2, 4, 1, NaN];
nmape = mape(ys, preds);
```
"""
function mape(y::Array, pred::Array)
    _mean = abs( nanmean(y) );
    _diff = abs.(y .- pred) ./ _mean .* 100;

    return nanmean( _diff )
end

"""

    mase(y::Array, pred::Array)

Return the mean absolute standardized error by ommiting the NaN, given
- `y` Array of numbers, can be NaN
- `pred` Array of predictions, can be NaN

---
Example
```julia
ys = [1, 2, 4, NaN];
preds = [2, 4, 1, NaN];
nmase = mase(ys, preds);
```
"""
function mase(y::Array, pred::Array)
    _nstd = nanstd(y);
    _diff = abs.(y .- pred) ./ _nstd .* 100;

    return nanmean( _diff )
end

"""

    rmse(y::Array, pred::Array)

Return the root mean square error by ommiting the NaN, given
- `y` Array of numbers, can be NaN
- `pred` Array of predictions, can be NaN

---
Example
```julia
ys = [1, 2, 4, NaN];
preds = [2, 4, 1, NaN];
rmse = rmse(ys, preds);
```
"""
function rmse(y::Array, pred::Array)
    return sqrt( nanmean( (y .- pred) .^ 2 ) )
end


text:
@testset " TimeSeries " begin

    using Random

    nt, nv = 10, 3
    xt = TimeSeries(nt, nv)

    @test length(xt.t) == nt
    @test typeof(xt.t) == Array{Float64,1}

    t = collect(0:10.0)
    u = [rand(nv) for i in eachindex(t)]
    yo = TimeSeries(t, u)

    @test typeof(yo.t) == Array{Float64,1}
    @test typeof(yo.u) == Array{Array{Float64,1},1}

end


text:
using LinearAlgebra, Distributions, JLD2, DataDeps
using OceanBoundaryLayerParameterEstimation

using ParameterEstimocean
using ParameterEstimocean.Parameters: closure_with_parameters
using ParameterEstimocean.Observations: set!

using Oceananigans
using Oceananigans.Units
using Oceananigans.OutputWriters: TimeInterval, JLD2OutputWriter
using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity, RiBasedVerticalDiffusivity
using Oceananigans: run!
using ParameterEstimocean.PseudoSteppingSchemes

regenerate_synthetic_observations = true
architecture = CPU()
# Δt = 5.0

# parameter_names = (:CᵂwΔ,  :Cᵂu★, :Cᴰ,
#                    :Cˢc,   :Cˢu,  :Cˢe,
#                    :Cᵇc,   :Cᵇu,  :Cᵇe,
#                    :Cᴷc⁻,  :Cᴷu⁻, :Cᴷe⁻,
#                    :Cᴷcʳ,  :Cᴷuʳ, :Cᴷeʳ,
#                    :CᴷRiᶜ, :CᴷRiʷ)

# parameter_set = ParameterSet{CATKEVerticalDiffusivity}(Set(parameter_names),
#                              nullify = Set([:Cᴬu, :Cᴬc, :Cᴬe]))

# closure = closure_with_parameters(CATKEVerticalDiffusivity(Float64;), parameter_set.settings)

# true_parameters = (CᵂwΔ = 4.46,    Cᵂu★ = 4.56,   Cᴰ = 2.91,
#                    Cˢc = 0.426,    Cˢu = 0.628,   Cˢe = 0.711,
#                    Cᵇc = 0.0723,   Cᵇu = 0.596,   Cᵇe = 0.637,
#                    Cᴷc⁻ = 0.343,   Cᴷu⁻ = 0.343,  Cᴷe⁻ = 1.42,
#                    Cᴷcʳ = -0.891,  Cᴷuʳ = -0.721, Cᴷeʳ = 1.50,
#                    CᴷRiᶜ = 2.03,   CᴷRiʷ = 0.101)

# true_closure = closure_with_parameters(closure, true_parameters)

# directory = "lesbrary_catke_perfect_model_2_days"
# isdir(directory) || mkpath(directory)

Δt = 10minutes

parameter_set = RiBasedParameterSet

closure = closure_with_parameters(RiBasedVerticalDiffusivity(), parameter_set.settings)
# closure = RiBasedVerticalDiffusivity()

true_parameters = parameter_set.settings
true_closure = closure

directory = "lesbrary_ri_based_perfect_model_2_days"
isdir(directory) || mkpath(directory)

###
### Generate data at high resolution using 6-day LESbrary as a template for ICs and BCs but CATKE "perfect model" parameters.
###

function run_synthetic_single_column_simulation!(observation, prefix;
                                                architecture = CPU(),
                                                closure = ConvectiveAdjustmentVerticalDiffusivity(),
                                                Δt = 10.0,
                                                tracers = (:b, :e)
                                            )

    grid = RectilinearGrid(architecture,
                            size = observation.grid.Nz,
                            topology = (Flat, Flat, Bounded),
                            z = (-observation.grid.Lz, 0))

    p = observation.metadata.parameters
    Qᵘ = p.boundary_condition_u_top
    Qᵇ = p.boundary_condition_θ_top
    N² = p.N²_deep

    u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))
    b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ), bottom = GradientBoundaryCondition(N²))

    model = HydrostaticFreeSurfaceModel(; grid, tracers, closure,
                                          buoyancy = BuoyancyTracer(),
                                          boundary_conditions = (u=u_bcs, b=b_bcs),
                                          coriolis = FPlane(f = observation.metadata.coriolis.f))

    set!(model, observation, 1)

    init_with_parameters(file, model) = file["parameters"] = (; Qᵇ, Qᵘ, Δt, N², tracers)

    simulation = Simulation(model; Δt, stop_time=observation.times[end])

    model = simulation.model

    simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers),
                                                          schedule = TimeInterval(0.25days),
                                                          array_type = Array{Float64},
                                                          prefix = prefix,
                                                          with_halos = true,
                                                          init = init_with_parameters,
                                                          force = true)

    run!(simulation)
end

fields_by_case = Dict(
   "free_convection" => (:b, :e),
   "weak_wind_strong_cooling" => (:b, :u, :v, :e),
   "strong_wind_weak_cooling" => (:b, :u, :v, :e),
   "strong_wind" => (:b, :u, :v, :e),
   "strong_wind_no_rotation" => (:b, :u, :e)
)

regenerate_synthetic_observations && begin

    transformation = (b = Transformation(normalization=ZScore()),
                    u = Transformation(normalization=ZScore()),
                    v = Transformation(normalization=ZScore()),
                    e = Transformation(normalization=RescaledZScore(1e-1)))

    # lesbrary_observations = SixDaySuite(; transformation=transformation, times=nothing, Nz=128)

    field_names = (:b, :u, :v, :e)

    two_day_suite_path(case) = "two_day_suite_2m/$(case)_instantaneous_statistics.jld2"

    for (case, forward_map_names) in zip(keys(fields_by_case), values(fields_by_case))

        lesbrary_data_path = @datadep_str two_day_suite_path(case)

        lesbrary_observation = SyntheticObservations(lesbrary_data_path; transformation,
                                                times=nothing,
                                                field_names=(:b, :u, :v, :e),
                                                forward_map_names)

        prefix = joinpath(directory, case)
        run_synthetic_single_column_simulation!(lesbrary_observation, prefix;
                                                architecture,
                                                closure=true_closure,
                                                Δt)
    end
end

###
### Build coarse-grained SyntheticObservations from data.
###

Nz = 64
Nensemble = 30

training_times = [0.5days, 0.75days, 1.0days]
validation_times = [1.0days, 1.25days, 1.5days]
testing_times = [1.5days, 1.75days, 2.0days]

path_fn(case) = joinpath(directory, case) * ".jld2"

begin
    file = jldopen(path_fn("free_convection"))
    @show iterations = parse.(Int, keys(file["timeseries/t"]))
    @show times = [file["timeseries/t/$i"] for i in iterations]
    close(file)
end

build_prior(name) = ScaledLogitNormal(bounds=bounds(name, parameter_set))
free_parameters = FreeParameters(named_tuple_map(parameter_set.names, build_prior))

function inverse_problem(Nensemble, times)
    observations = SyntheticObservationsBatch(path_fn, transformation, times, Nz; datadep=false, architecture)

    simulation = ensemble_column_model_simulation(observations;
                                                  Nensemble,
                                                  architecture,
                                                  tracers = (:b, :e),
                                                  closure)

    simulation.Δt = Δt

    inverse_problem = InverseProblem(observations, simulation, free_parameters; output_map = ConcatenatedOutputMap())
    return inverse_problem
end

###
### Generate preliminary inverse problems.
###

training = inverse_problem(Nensemble, training_times)
validation = inverse_problem(Nensemble, validation_times)
testing = inverse_problem(Nensemble, testing_times)

resampler = Resampler(resample_failure_fraction=0.5, acceptable_failure_fraction=1.0)
eki = EnsembleKalmanInversion(training; resampler)

stepping_scheme = ConstantConvergence(convergence_ratio = 0.9)
params = iterate!(eki; iterations = 10, pseudo_stepping = stepping_scheme)

###
### Visualize
###

visualize!(training, true_parameters;
                    field_names = [:u, :v, :b, :e],
                    directory,
                    filename = "true_parameter_training_realizations.png"
                    )

visualize!(training, params;
                    field_names = [:u, :v, :b, :e],
                    directory,
                    filename = "perfect_model_visual_calibrated.png"
                    )

plot_parameter_convergence!(eki, directory; true_parameters, n_columns=3)
plot_pairwise_ensembles!(eki, directory, true_parameters)
plot_error_convergence!(eki, directory, true_parameters)