text:
"""
    $(TYPEDEF)

Pickup a block using a robot arm modeled after the [Modular Prosthetic Limb]
(https://www.jhuapl.edu/Content/techdigest/pdf/V30-N03/30-3-Johannes.pdf)
developed by the Applied Physics Laboratory, The Johns Hopkins University.

# Spaces

* **State: (106, )**
* **Action: (36, )**
* **Observation: (19, )**
"""
struct ArmHandPickup{Sim,OSpace} <: AbstractMuJoCoEnvironment
    sim::Sim
    observationspace::OSpace
    ball::Int
    palm::Int
    thumb::Int
    index::Int
    middle::Int
    ring::Int
    pinky::Int
    goal::SVector{3,Float64}

    function ArmHandPickup(sim::MJSim)
        m = sim.m
        mn, dn = sim.mn, sim.dn

        observationspace = MultiShape(
            d_thumb = ScalarShape(Float64),
            d_index = ScalarShape(Float64),
            d_middle = ScalarShape(Float64),
            d_ring = ScalarShape(Float64),
            d_pinky = ScalarShape(Float64),
            a_thumb = ScalarShape(Float64),
            a_index = ScalarShape(Float64),
            a_middle = ScalarShape(Float64),
            a_ring = ScalarShape(Float64),
            a_pinky = ScalarShape(Float64),
            a_close = ScalarShape(Float64),
            handball = ScalarShape(Float64),
            ballgoal = ScalarShape(Float64),
            ball = VectorShape(Float64, 3),
            palm = VectorShape(Float64, 3),
        )

        ball = jl_name2id(sim.m, MuJoCo.MJCore.mjOBJ_SITE, "ball")
        palm = jl_name2id(sim.m, MuJoCo.MJCore.mjOBJ_SITE, "palm")

        thumb = jl_name2id(sim.m, MuJoCo.MJCore.mjOBJ_SITE, "thumb_IMU")
        index = jl_name2id(sim.m, MuJoCo.MJCore.mjOBJ_SITE, "index_IMU")
        middle = jl_name2id(sim.m, MuJoCo.MJCore.mjOBJ_SITE, "middle_IMU")
        ring = jl_name2id(sim.m, MuJoCo.MJCore.mjOBJ_SITE, "ring_IMU")
        pinky = jl_name2id(sim.m, MuJoCo.MJCore.mjOBJ_SITE, "pinky_IMU")

        goal = SA_F64[0.0, 0.2, 0.5]

        new{typeof(sim),typeof(observationspace)}(
            sim,
            observationspace,
            ball,
            palm,
            thumb,
            index,
            middle,
            ring,
            pinky,
            goal,
        )
    end
end

ArmHandPickup() = first(tconstruct(ArmHandPickup, 1))

function LyceumBase.tconstruct(::Type{ArmHandPickup}, n::Integer)
    modelpath = joinpath(@__DIR__, "armhand.xml")
    Tuple(ArmHandPickup(s) for s in tconstruct(MJSim, n, modelpath, skip = 3))
end

@inline getsim(env::ArmHandPickup) = env.sim

@inline LyceumBase.observationspace(env::ArmHandPickup) = env.observationspace

@propagate_inbounds function LyceumBase.getobservation!(obs, env::ArmHandPickup)
    @boundscheck checkaxes(observationspace(env), obs)

    m, d = env.sim.m, env.sim.d
    sx = d.site_xpos
    dmin = 0.5

    ball = SPoint3D(sx, env.ball)
    palm = SPoint3D(sx, env.palm)
    thumb = SPoint3D(sx, env.thumb)
    index = SPoint3D(sx, env.index)
    middle = SPoint3D(sx, env.middle)
    ring = SPoint3D(sx, env.ring)
    pinky = SPoint3D(sx, env.pinky)
    goal = ball - env.goal

    @uviews obs @inbounds begin
        shaped = observationspace(env)(obs)

        shaped.ball .= ball
        shaped.palm .= palm .- ball

        shaped.handball = _sitedist(palm, ball, dmin)
        shaped.ballgoal = _sitedist(ball, goal, dmin)

        shaped.d_thumb = _sitedist(thumb, ball, dmin)
        shaped.d_index = _sitedist(index, ball, dmin)
        shaped.d_middle = _sitedist(middle, ball, dmin)
        shaped.d_ring = _sitedist(ring, ball, dmin)
        shaped.d_pinky = _sitedist(pinky, ball, dmin)
        shaped.a_thumb = cosine_dist(goal, thumb)
        shaped.a_index = cosine_dist(goal, index)
        shaped.a_middle = cosine_dist(goal, middle)
        shaped.a_ring = cosine_dist(goal, ring)
        shaped.a_pinky = cosine_dist(goal, pinky)
        shaped.a_close = cosine_dist(middle, thumb)
    end

    obs
end

@propagate_inbounds function LyceumBase.getreward(state, action, obs, env::ArmHandPickup)
    @boundscheck checkaxes(observationspace(env), obs)

    os = observationspace(env)(obs)
    handball = os.handball / 0.5
    ballgoal = os.ballgoal / 0.5

    reward = -handball
    if handball < 0.06
        reward = 2.0 - 2 * ballgoal
        reward -= 0.1 * (os.d_thumb + os.d_index + os.d_middle + os.d_ring + os.d_pinky)
    end
    reward
end

@propagate_inbounds function LyceumBase.randreset!(rng::Random.AbstractRNG, env::ArmHandPickup)
    fastreset_nofwd!(env.sim)
    env.sim.d.qpos[1] = rand(rng, Uniform(-0.15, 0.15))
    env.sim.d.qpos[2] = rand(rng, Uniform(-0.1, 0.1))
    forward!(env.sim)
    env
end

@inline _sitedist(s1, s2, dmin) = min(euclidean(s1, s2), dmin)


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# No regularization
# ------------------------------------------------------------
"""
    NoRegularization

Empty functor object for running an L-shaped algorithm without regularization.

"""
struct NoRegularization <: AbstractRegularization end

function initialize_regularization!(::AbstractLShaped, ::NoRegularization)
    return nothing
end

function restore_regularized_master!(::AbstractLShaped, ::NoRegularization)
    return nothing
end

function filter_variables!(::NoRegularization, ::Any)
    return nothing
end

function filter_constraints!(::NoRegularization, ::Any)
    return nothing
end

function log_regularization!(::AbstractLShaped, ::NoRegularization)
    return nothing
end

function log_regularization!(::AbstractLShaped, ::Integer, ::NoRegularization)
    return nothing
end

function take_step!(::AbstractLShaped, ::NoRegularization)
    return nothing
end

function decision(lshaped::AbstractLShaped, ::NoRegularization)
    return lshaped.x
end

function objective_value(lshaped::AbstractLShaped, ::NoRegularization)
    return lshaped.data.Q
end

function gap(lshaped::AbstractLShaped, ::NoRegularization)
    @unpack Q,θ = lshaped.data
    return abs(θ-Q)/(abs(Q)+1e-10)
end

# API
# ------------------------------------------------------------
"""
    DontRegularize

Factory object for [`NoRegularization`](@ref). Passed by default to `regularize` in `LShaped.Optimizer`.

"""
struct DontRegularize <: AbstractRegularizer end

function (::DontRegularize)(::DecisionMap, ::AbstractVector)
    return NoRegularization()
end

function str(::DontRegularize)
    return "L-shaped"
end

ext: jl


text:
export RandomizedNewton

"""
Randomized Newton
-> for symmetric positive definite matrix

Section 3.6 in Gower, R. M., & Richtárik, P. (2015).
Randomized iterative methods for linear systems.
SIAM Journal on Matrix Analysis and Applications, 36(4), 1660-1690.
"""
function RandomizedNewton(
  stp::AbstractStopping;
  is_zero_start::Bool = false,
  r::Int = 15,
  rand_r::Bool = false,
  kwargs...,
)
  A, b = get_matrix(stp.pb), get_vector(stp.pb)
  m, n = size(A)
  state = stp.current_state
  if is_zero_start
    state.res .= b
  else
    state.res .= b .- A * state.x
  end
  OK = start!(stp, no_opt_check = true)

  while !OK
    r = rand_r ? Int(floor(rand() * n) + 1) : r
    sub = zeros(Int64, r)
    sample!(1:n, sub) #x is a subset of [1,...,n] of size r

    Ai = A[sub, :]
    resk = state.res[sub] #resk = Ai*xk - b[sub]

    if 0 ∉ resk
      state.x .-= Ai \ resk
    end

    state.res .= b .- A * state.x
    OK = cheap_stop!(stp)
  end

  return stp
end


text:
function coreLaplace(logposterior::Function,
                 gradlogposterior::Function,
              hessianlogposterior::Function,
                               w0::Array{Float64,1};
                       iterations = iterations,
                        optimiser = optimiser,
                       show_every = show_every)

    @printf("Running laplace with D=%d for %d iterations\n", length(w0), iterations)

    negative_logposterior = x-> -1.0 * logposterior(x)

    negative_logposterior_gradient!(storage, w) = copyto!(storage, -1.0 * gradlogposterior(w))

    # find mode μ, we are *minimising*
    μ = Optim.optimize(negative_logposterior, negative_logposterior_gradient!, w0, optimiser,
            Optim.Options(iterations=iterations, show_trace=true, show_every=show_every)).minimizer

    # Covariance of Gaussian posterior at mode μ is the inverse of Hessian of negative log-likelihood
    Σ = get_covariance_at_mode(μ, x->-1*hessianlogposterior(x))

    # return Gaussian parameters
    return MvNormal(μ, Σ)

end

######################################################################
function get_covariance_at_mode(mode, hessianlogposterior)
######################################################################

    # find Hessian at mode
    H = hessianlogposterior(mode)

    # Covariance of Gaussian posterior at mode is the inverse of Hessian
    Σ = try

        ((H+H')/2)\I

    catch err
        if isa(err, LinearAlgebra.LAPACKException)
            local U,S,V = svd(H)
            V*Diagonal(1.0./max.(S, 1e-6))*U'
        else
            throw(err)
        end
    end

    # Enforce symmetry
    Σ = (Σ + Σ') * 0.5 + 1e-8*I

    return Σ

end


text:
import ChainRulesCore: rrule, @non_differentiable, NoTangent
using Yao, Yao.AD

function rrule(::typeof(apply), reg::ArrayReg, block::AbstractBlock)
    out = apply(reg, block)
    out, function (outδ)
        (in, inδ), paramsδ = apply_back((copy(out), outδ), block)
        return (NoTangent(), inδ, paramsδ)
    end
end

function rrule(::typeof(dispatch), block::AbstractBlock, params)
    out = dispatch(block, params)
    out, function (outδ)
        (NoTangent(), NoTangent(), outδ)
    end
end

function rrule(::typeof(expect), op::AbstractBlock, reg::AbstractRegister{B}) where {B}
    out = expect(op, reg)
    out, function (outδ)
        greg = Yao.AD.expect_g(op, reg)
        for b=1:B
            viewbatch(greg, b).state .*= 2*outδ[b]
        end
        return (NoTangent(), NoTangent(), greg)
    end
end

function rrule(::Type{Matrix}, block::AbstractBlock)
    out = Matrix(block)
    out, function (outδ)
        paramsδ = mat_back(block, outδ)
        return (NoTangent(), paramsδ)
    end
end

function rrule(::Type{ArrayReg{B}}, raw::AbstractArray) where B
    ArrayReg{B}(raw), adjy->(NoTangent(), reshape(adjy.state, size(raw)))
end

function rrule(::Type{ArrayReg{B}}, raw::ArrayReg) where B
    ArrayReg{B}(raw), adjy->(NoTangent(), adjy)
end

function rrule(::Type{ArrayReg}, raw::AbstractArray)
    ArrayReg(raw), adjy->(NoTangent(), reshape(adjy.state, size(raw)))
end

function rrule(::typeof(copy), reg::ArrayReg) where B
    copy(reg), adjy->(NoTangent(), adjy)
end

_totype(::Type{T}, x::AbstractArray{T}) where T = x
_totype(::Type{T}, x::AbstractArray{T2}) where {T,T2} = convert.(T, x)
rrule(::typeof(state), reg::ArrayReg{B,T}) where {B,T} = state(reg), adjy->(NoTangent(), ArrayReg(_totype(T, adjy)))
rrule(::typeof(statevec), reg::ArrayReg{B,T}) where {B,T} = statevec(reg), adjy->(NoTangent(), ArrayReg(_totype(T, adjy)))
rrule(::typeof(state), reg::AdjointArrayReg{B,T}) where {B,T} = state(reg), adjy->(NoTangent(), ArrayReg(_totype(T, adjy)')')
rrule(::typeof(statevec), reg::AdjointArrayReg{B,T}) where {B,T} = statevec(reg), adjy->(NoTangent(), ArrayReg(_totype(adjy)')')
rrule(::typeof(parent), reg::AdjointArrayReg) = parent(reg), adjy->(NoTangent(), adjy')
rrule(::typeof(Base.adjoint), reg::ArrayReg) = Base.adjoint(reg), adjy->(NoTangent(), parent(adjy))
@non_differentiable Yao.nparameters(::Any)
@non_differentiable Yao.zero_state(args...)
@non_differentiable Yao.rand_state(args...)
@non_differentiable Yao.uniform_state(args...)
@non_differentiable Yao.product_state(args...)


text:
###########
# simple project activate
using Pkg
path_project=join(split(@__DIR__, "/")[1:end-1],"/")
println(path_project)
Pkg.activate(path_project)
###########
include("../src/asm.jl")
include("../src/branching_network.jl")
include("../src/utils.jl")
include("../src/bivariate_gauss.jl")

using Statistics
using LinearAlgebra
using RollingFunctions
using Printf
using StatsBase

"""

measure mean and variance of correlated but short timeseries
"""
function produce_data_sampling_bias_illustration(;
        logTs = 0:0.2:5,
        tau = 1e1,
        trials::Int = Int(2e4),
    )
    rng = MersenneTwister(1000)

    # correlation
    rho = exp(-1/tau)

    result = Dict{String, Any}()
    result["T"] = zeros(length(logTs))
    result["mean"] = zeros(trials, length(logTs))
    result["var"]  = zeros(trials, length(logTs))
    for (i, logT) in enumerate(logTs)
        T = round(Int,10^logT)
        result["T"][i] = T
        println(T)
        println("logT = ", log10(T), " ~ ", logT)
        @showprogress 1 for j in 1:trials
            data = bivariate_gauss(rng, rho, T)
            result["mean"][j,i] = mean(data)
            result["var"][j,i]  = var(data, mean=result["mean"][j,i])
        end
    end
    return result
end

"""
    produce_data_avalanches_sandpile_illustration()

produce illustration of subsampling of 2D system as projection of activity to
1D (y axis) as a function of time (x-axis)
"""
function produce_data_avalanches_sandpile_illustration(;
        L::Int=Int(2^5), # linear extension of lattice
        Ls=(L,L),   # 2D lattice
        hc::Int=4,
        p_dis::Float64=1e-4, # dissipation probability per toppling (probably system size dependent)
        seed=1000,  # seed for random number generator
        periodic=true
    )
    rng = MersenneTwister(seed)

    lattice = LightGraphs.grid(Ls, periodic=periodic)
    hs = rand(rng, 0:hc-1, Ls)

    subsample = Dict{String, Any}("l" => Int[2^3])
    subsample["random"] = [Vector{Int}(undef,l*l) for l in subsample["l"]]
    subsample["window"] = [Vector{Int}(undef,l*l) for l in subsample["l"]]
    # fill the subsampling schemes
    for (i,l) in enumerate(subsample["l"])
        subsample["random"][i] .= sample(rng, vertices(lattice), length(subsample["random"][i]), replace=false)
        range =  floor(Int,L/2)-floor(Int,l/2)+1 : floor(Int,L/2)-floor(Int,l/2)+l
        subsample["window"][i] .= vcat(LinearIndices(size(hs))[range,range]...)
    end

    num_valid_avalanches = 0
    data = Dict()
    data["mean_h"] = Float64[]

    min_size_sub  = Int(10)
    min_size_full = Int(30)
    max_size_full = Int(1e3)
    println("find a suitable avalanche")
    center = LinearIndices(size(hs))[floor(Int, L/2),floor(Int, L/2)]
    for num_tries in 1:Int(1e5)
        # get list of lists of toppled sites (one list of toppled sites per time step)
        topple_sites_over_time = perturb_and_relax(rng, hs, lattice, hc, p_dis, start_site=center)

        println("size total = ", length(topple_sites_over_time))
        if (length(topple_sites_over_time) > min_size_full) && (length(topple_sites_over_time) < max_size_full)
            # get avalanche sizes for subsamples
            size_rand = 0
            size_wind = 0
            for sites in topple_sites_over_time
                size_rand += sum(in.(sites, Ref(subsample["random"][1])))
                size_wind += sum(in.(sites, Ref(subsample["window"][1])))
            end
            println("size rand = ", size_rand)
            println("size wind = ", size_wind)
            if (size_rand > min_size_sub) && (size_wind > min_size_sub)
                #want 2D array of (toppled sites projected, time)
                T = length(topple_sites_over_time)
                result = Dict{String, Any}()
                # subsamples
                result["subsample/rand"] = subsample["random"][1]
                result["subsample/wind"] = subsample["window"][1]
                #log
                result["states/full"] = zeros(Ls...,T)
                #time
                result["time/full"] = zeros((L,T))
                result["time/rand"] = zeros((L,T))
                result["time/wind"] = zeros((L,T))
                #space
                result["space/full"] = zeros(Ls)
                result["space/rand"] = zeros(Ls)
                result["space/wind"] = zeros(Ls)
                #layer
                result["layer/full"] = zeros(Ls)
                result["layer/rand"] = zeros(Ls)
                result["layer/wind"] = zeros(Ls)
                result["layer/full"][result["subsample/rand"]] .= 1
                result["layer/wind"][result["subsample/wind"]] .= 1
                for (t,sites) in enumerate(topple_sites_over_time)
                    sites_full = sites
                    sites_rand = sites[in.(sites, Ref(subsample["random"][1]))]
                    sites_wind = sites[in.(sites, Ref(subsample["window"][1]))]
                    result["states/full"][CartesianIndices(Ls)[sites_full],t] .= 1
                    result["space/full"][sites_full] .+= 1
                    result["space/rand"][sites_rand] .+= 1
                    result["space/wind"][sites_wind] .+= 1
                    result["time/full"][:,t] .= projection_from_sites(sites_full, Ls)
                    result["time/rand"][:,t] .= projection_from_sites(sites_rand, Ls)
                    result["time/wind"][:,t] .= projection_from_sites(sites_wind, Ls)
                end
                return result
            end
        end
    end
end

function projection_from_sites(active_sites, Ls)
    state = zeros(Ls)
    state[active_sites] .= 1
    return vec(sum(state, dims=1))
end

"""
    produce_data_avalanches_sandpile()

runs the default simulation of an abelian sandpile model with cyclic (periodic)
boundary conditions. Parameter can be specified if desired.

Model is Abelian sandpile plus periodic boundary conditions and the chance
that toppling does not propagate grains (to make it subcritical)

L=4;
Ls=(L,L);
seed=1000;
hc=4;

"""
function produce_data_avalanches_sandpile(;
        L::Int=Int(2^7), # linear extension of lattice
        Ls=(L,L),   # 2D lattice
        hc::Int=4,
        p_dis::Float64=1e-4, # dissipation probability per toppling (probably system size dependent)
        seed=1000,  # seed for random number generator
        num_avalanches::Int = Int(1e6),
        periodic=true
    )
    rng = MersenneTwister(seed)

    lattice = grid(Ls, periodic=periodic)
    hs = rand(rng, 0:hc-1, Ls)

    subsample = Dict{String, Any}("l" => Int[2^2, 2^4, 2^6])
    subsample["random"] = [Vector{Int}(undef,l*l) for l in subsample["l"]]
    subsample["window"] = [Vector{Int}(undef,l*l) for l in subsample["l"]]
    # fill the subsampling schemes
    for (i,l) in enumerate(subsample["l"])
        subsample["random"][i] .= sample(rng, vertices(lattice), length(subsample["random"][i]), replace=false)
        range =  floor(Int,L/2)-floor(Int,l/2)+1 : floor(Int,L/2)-floor(Int,l/2)+l
        subsample["window"][i] .= vcat(LinearIndices(size(hs))[range,range]...)
    end

    num_valid_avalanches = 0
    data = Dict()
    data["mean_h"] = Float64[]

    data["duration"] = Dict{String,Any}()
    data["duration"]["full"] = zeros(Int64, num_avalanches)

    data["size"] = Dict{String,Any}()
    data["size"]["full"] = zeros(Int64, num_avalanches)
    data["size"]["random"] = [zeros(Int64, num_avalanches) for i in 1:length(subsample["random"])]
    data["size"]["window"] = [zeros(Int64, num_avalanches) for i in 1:length(subsample["window"])]

    response_steps = 1000
    data["response"] = Dict{String, Any}()
    data["response"]["full"] = zeros(Int64, num_avalanches, response_steps)
    data["response"]["random"] = [zeros(Int64, num_avalanches, response_steps) for i in 1:length(subsample["random"])]
    data["response"]["window"] = [zeros(Int64, num_avalanches, response_steps) for i in 1:length(subsample["window"])]

    progress = Progress(num_avalanches, 1, "sampling: ", offset=0)
    while num_valid_avalanches < num_avalanches
        push!(data["mean_h"], mean(hs))
        # get list of lists of toppled sites (one list of toppled sites per time step)
        topple_sites_over_time = perturb_and_relax(rng, hs, lattice, hc, p_dis)
        ProgressMeter.update!(progress, num_valid_avalanches)
        if length(topple_sites_over_time) > 0
            num_valid_avalanches += 1

            # full sampling
            data["duration"]["full"][num_valid_avalanches] = length(topple_sites_over_time)
            data["size"]["full"][num_valid_avalanches] = sum(length.(topple_sites_over_time))

            for (t,sites) in enumerate(topple_sites_over_time)
                if t <= response_steps
                    data["response"]["full"][num_valid_avalanches,t] = length(sites)
                end
                # random spatial subsampling
                for (i, test_sample) in enumerate(subsample["random"])
                    a_t = sum(in.(sites, Ref(test_sample)))
                    data["size"]["random"][i][num_valid_avalanches] += a_t
                    if t <= response_steps
                        data["response"]["random"][i][num_valid_avalanches,t] = a_t
                    end
                end

                # field-of-view spatial subsampling
                for (i, test_sample) in enumerate(subsample["window"])
                    a_t = sum(in.(sites, Ref(test_sample)))
                    data["size"]["window"][i][num_valid_avalanches] += a_t
                    if t <= response_steps
                        data["response"]["window"][i][num_valid_avalanches,t] = a_t
                    end
                end
            end
        end
    end

    return subsample, data
end
#function produce_data_perturb_relax_bn(;
#        ref=4,
#        L::Int=Int(2^7), # linear extension of lattice
#        Ls=(L,L),          # 1D lattice
#        p_rec::Float64=0.644, #
#        seed=1000,  # seed for random number generator
#        num_avalanches::Int = Int(1e5),
#        periodic=true
#    )
#    rng = MersenneTwister(seed)
#
#    lattice = grid(Ls, periodic=periodic)
#    state = zeros(Bool, Ls)
#    state_ref = zeros(Int, Ls)
#    state_new = zeros(Bool, Ls)
#
#    subsample = Dict{String, Any}("l" => Int[2^2, 2^4, 2^6])
#    subsample["random"] = [Vector{Int}(undef,l*l) for l in subsample["l"]]
#    subsample["window"] = [Vector{Int}(undef,l*l) for l in subsample["l"]]
#    # fill the subsampling schemes
#    for (i,l) in enumerate(subsample["l"])
#        subsample["random"][i] .= sample(rng, vertices(lattice), length(subsample["random"][i]), replace=false)
#        range =  floor(Int,L/2)-floor(Int,l/2)+1 : floor(Int,L/2)-floor(Int,l/2)+l
#        subsample["window"][i] .= vcat(LinearIndices(size(state))[range, range]...)
#    end
#
#    num_valid_avalanches = 0
#    data = Dict()
#    data["duration"] = Dict{String,Any}()
#    data["duration"]["full"] = zeros(Int64, num_avalanches)
#
#    data["size"] = Dict{String,Any}()
#    data["size"]["full"] = zeros(Int64, num_avalanches)
#    data["size"]["random"] = [zeros(Int64, num_avalanches) for i in 1:length(subsample["random"])]
#    data["size"]["window"] = [zeros(Int64, num_avalanches) for i in 1:length(subsample["window"])]
#
#    progress = Progress(num_avalanches, 1, "sampling: ", offset=0)
#    for num_valid_avalanches in 1:num_avalanches
#        next!(progress)
#
#        # single perturbation
#        state_ref .= 0
#        state[rand(rng, vertices(lattice))] = true
#        state_ref[state] .= ref
#        activity = sum(state)
#        while !(activity == 0)
#            # full sampling
#            data["duration"]["full"][num_valid_avalanches] += 1
#            data["size"]["full"][num_valid_avalanches] += activity
#            # random spatial subsampling
#            for (i, test_sample) in enumerate(subsample["random"])
#                a_t = sum(state[test_sample])
#                data["size"]["random"][i][num_valid_avalanches] += a_t
#            end
#
#            # field-of-view spatial subsampling
#            for (i, test_sample) in enumerate(subsample["window"])
#                a_t = sum(state[test_sample])
#                data["size"]["window"][i][num_valid_avalanches] += a_t
#            end
#            # next step
#            step!(rng, state, state_new, state_ref, ref, lattice, p_rec)
#            activity = sum(state)
#        end
#    end
#
#    return data
#end

"""
    produce_data_box_scaling_branching_network(m, h)

simulation of branching network on 2D lattice with nearest neighbor interaction.
Returns estimates of the connected correlation function for box scaling analysis.
"""
function produce_data_box_scaling_branching_network(m, h::Float64;
        L::Int=Int(2^10), # linear extension of lattice
        Ls=(L,L),        # 2D lattice
        seed=1000,  # seed for random number generator
        steps_total::Int = Int(1e5),
        steps_equil::Int = Int(4e4),
        dstep_correlation::Int = 100,
    )
    string(l) = @sprintf("l=%d",l)

    dt=1
    p_ext = 1-exp(-h*dt)

    # initialize random number generator
    rng = MersenneTwister(seed)

    # construct 2D lattice with additional neighbors along the diagonal
    lattice = lattice_nearest_neighbor(Ls[1],Ls[2]);
    p_rec = m/mean(degree(lattice));
    state = zeros(Bool, Ls);
    state_new = zeros(Bool, Ls);
    #state_ref = zeros(Int, Ls);
    #val_ref = 2;

    # setup the spatial subsampling schemes
    data = Dict{String, Any}()
    data["l"] = Int[L/2, L/2^2, L/2^3, L/2^4, L/2^5, L/2^6]
    for l in data["l"]
        data["random/"*string(l)] = sample(rng, vertices(lattice), l*l, replace=false)
        range =  floor(Int,L/2)-floor(Int,l/2)+1 : floor(Int,L/2)-floor(Int,l/2)+l
        data["window/"*string(l)] =  LinearIndices(size(state))[range,range]
    end

    # prepare output dictionary
    # number of active sites over time
    data["activity/full"] = zeros(Int64, steps_total)
    for l in data["l"]
        data["activity/random/"*string(l)] = zeros(Int64, steps_total)
        data["activity/window/"*string(l)] = zeros(Int64, steps_total)
    end
    # time average of spatial correlation function
    data["correlation/full_lag"] = collect(0:floor(Int,L/2))
    data["correlation/full"] = zeros(Float64, length(data["correlation/full_lag"]))
    for l in data["l"]
        l_range = collect(0:floor(Int,l/2))
        data["correlation/window_lag/"*string(l)] = l_range
        data["correlation/window/"*string(l)] = zeros(Float64, length(l_range))
    end

    # actual simulation loop
    progress = Progress(steps_total, 1, "sampling: ", offset=0)
    number_estimates_correlation = 0
    for s in 1:steps_total
        next!(progress)
        # external input independently added to state_new
        for i in 1:length(state)
            if rand(rng) < p_ext
                state_new[i] = true
            end
        end
        # recurrent update step
        step!(rng, state, state_new, lattice, p_rec)
        #step!(rng, state, state_new, lattice, p_rec, state_ref, val_ref)

        # activity
        data["activity/full"][s] = sum(state)
        for l in data["l"]
            # random spatial subsampling
            state_sub = state[data["random/"*string(l)]]
            data["activity/random/"*string(l)][s] = sum(state_sub)
            # field-of-view spatial subsampling
            state_sub = state[data["window/"*string(l)]]
            data["activity/window/"*string(l)][s] = sum(state_sub)
        end

        # start measurements after initial equilibration
        # (needs to be checked afterwards with the time trance of active sites)
        if s > steps_equil
            # correlation estimates
            if s%dstep_correlation ==0
                number_estimates_correlation += 1
                state_cor = state
                data["correlation/full"] .+= correlation_x(state_cor,data["correlation/full_lag"])
                for l in data["l"]
                    state_cor = state[data["window/"*string(l)]]
                    data["correlation/window/"*string(l)] .+= correlation_x(state_cor, data["correlation/window_lag/"*string(l)])
                end
            end
        end
    end
    # normalization of correlation function estimate
    data["correlation/full"] ./= number_estimates_correlation
    for l in data["l"]
        data["correlation/window/"*string(l)] ./= number_estimates_correlation
    end

    return data
end


text:
# License for this file: MIT (expat)
# Copyright 2020-2021, DLR Institute of System Dynamics and Control
# Author: Martin Otter, DLR-SR

#=
Algorithm to exactly process linear Integer equations and hereby removing underdeterminism, overdeterminism,
and state constraints that are not structurally visible, as well as eliminating simple equations.

The following is exported:

- [`simplifyLinearIntegerEquations`](@ref)
- [`printSimplifiedLinearIntegerEquations`](@ref)

and the following utility functions that can be applied on the return
argument `vProperty` of [`simplifyLinearIntegerEquations`](@ref):

- `isNotEliminated(vProperty, v)` - if variable v is not eliminated.
- `isEliminated(vProperty, v)` - if variable v is eliminated.
- `isZero(vProperty,v)` - if eliminated variable v is zero.
- `isAlias(vProperty,v)` - if eliminated variable v is an alias variable v = v_alias
- `isNegAlias(vProperty,v)` - if eliminated variable v is a negative alias variable v = -v_alias
- `alias(vProperty,v)` - alias variable v_alias of eliminated variable v (v = v_alias).
- `negAlias(vProperty,v)` - negated alias variable v_alias of eliminated variable v (v = -v_alias).

# Main developer

[Martin Otter](https://rmc.dlr.de/sr/en/staff/martin.otter/),
[DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)
=#

import OrderedCollections

export simplifyLinearIntegerEquations!
export printSimplifiedLinearIntegerEquations

export isNotEliminated, isEliminated, isZero, isAlias, isNegAlias, alias, negAlias

# Inquire properties of eliminated variables (vEliminated is returned from simplifyLinearIntegerEquations(..)
const IS_PRESENT = typemin(Int)   # Variable is not eliminated

isNotEliminated(vProperty::Vector{Int}, v::Int) =  vProperty[v] == IS_PRESENT
isEliminated(   vProperty::Vector{Int}, v::Int) =  vProperty[v] != IS_PRESENT
isZero(         vProperty::Vector{Int}, v::Int) =  vProperty[v] == 0
isAlias(        vProperty::Vector{Int}, v::Int) =  vProperty[v] > 0
isNegAlias(     vProperty::Vector{Int}, v::Int) =  vProperty[v] < 0 && vProperty[v] != IS_PRESENT
alias(          vProperty::Vector{Int}, v::Int) =  vProperty[v]
negAlias(       vProperty::Vector{Int}, v::Int) = -vProperty[v]

"""
    getFirstActiveIndex(Gint, vActive, i::Int)

Return the first index j of equation i, such that vActive[ Gint[i][j] ] == true.
If there is no such index, return zero.
"""
function getFirstActiveIndex(Gint, vActive::Vector{Bool}, i::Int)::Int
    found = false
    j = 0
    for (k,vk) in enumerate(Gint[i])
        if vActive[vk]
            j = k
            break
        end
    end
    return j
end

"""
    swapEquations!(Gint, eInt, GcInt, GcInt2, i, j)

Swap equations i and j
"""
function swapEquations!(Gint, eInt::Vector{Int}, GcInt, GcInt2, i::Int, j::Int)::Nothing
    Gint_i   = Gint[i]
    GcInt_i  = GcInt[i]
    GcInt2_i = GcInt2[i]
    eInt_i   = eInt[i]

    Gint[i]   = Gint[j]
    GcInt[i]  = GcInt[j]
    GcInt2[i] = GcInt2[j]
    eInt[i]   = eInt[j]

    Gint[j]   = Gint_i
    GcInt[j]  = GcInt_i
    GcInt2[j] = GcInt2_i
    eInt[j]   = eInt_i
    return nothing
end

"""
    cj = find_v(eq_i, eqc_i, vj)

Return coefficient of variable vj in eq_i or zero, if vj not in eq_i.
"""
function find_v(eq_i, eqc_i, vj)::Int
    for (j, vi) in enumerate(eq_i)
        if vi == vj
            return eqc_i[j]
        end
    end
    return 0
end

mutable struct Buffer
    Gint_i::Vector{Int}
    GcInt2_i::Vector{Int}

    Buffer() = new(Int[], Int[])
end

"""
    eliminateVariable!(Gint, GcInt2, k, k_perm, i, pivot, vj, oldPivot)

Eliminate variable vj from equation i and store the result as equation i.

k_perm is the sorting order of equation k (return argument of sortperm(..)).

If equation i does not contain variable vj, the equation is not changed.

If equation i contains variable vj, then multiply equation i with pivot and
equation k with the coefficient of variable vj in equation i.
Subtract the latter from the former equation and divide by oldPivot
(it is guaranteed that the division has no reminder).

If Gc[i,j] would be the coefficient of variable j in equation i, then the
following operation is carried out:

pivot_i = Gc[i,vj]
for j = 1:nv  # nv: number of variables
    Gc[i,j] = div(pivot * Gc[i,j] - pivot_i * Gc[k,j], oldPivot)  # guaranteed no remainder
end

If one of the coeffients of equation i, so GcInt2[i], becomes zero, then this coefficient and the corresponding
variable is removed from the equation. It is therefore guaranteed that after the operation,
no coefficient of the new equation i is zero, so no element of GcInt2[i][:] is zero
and no element of Gint[i][:] is vj.
"""
function eliminateVariable!(Gint, GcInt2, k::Int, i::Int, pivot::Int, vj::Int, oldPivot::Int, buffer::Buffer)::Nothing
    # Return if equation i does not contain variable vj
    vj_present=false
    c_vj = 0   # Coefficient of vj
    for (index, value) in enumerate(Gint[i])
        if value == vj
            vj_present = true
            c_vj = GcInt2[i][index]
            deleteat!(Gint[i] ,index)
            deleteat!(GcInt2[i],index)
            break
        end
    end
    if !vj_present
        return nothing
    end

    #println("... Eliminate variable $vj from equation $i = ", Gint[i])
    #println("    Gint[$k] = ", Gint[k], ", GcInt2[$k] = ", GcInt2[k], ", pivot = $pivot, c_vj = $c_vj")

    # Union of relevant variables
    eq_i   = Gint[i]
    eq_k   = Gint[k]
    eqc_i  = GcInt2[i]
    eqc_k  = GcInt2[k]
    v_all  = union(eq_i, eq_k)
    #println("v_all = $v_all")

    # Eliminate variable vj from equation i
    empty!(buffer.Gint_i)
    empty!(buffer.GcInt2_i)
    for v in v_all
        ci = find_v(eq_i, eqc_i, v)
        ck = find_v(eq_k, eqc_k, v)
        ci_new = div(pivot*ci - c_vj*ck, oldPivot)
        if ci_new != 0
            push!(buffer.Gint_i , v)
            push!(buffer.GcInt2_i, ci_new)
        end
    end
    temp1     = Gint[i]
    temp2     = GcInt2[i]
    Gint[i]   = buffer.Gint_i
    GcInt2[i] = buffer.GcInt2_i
    buffer.Gint_i   = temp1
    buffer.GcInt2_i = temp2
    return nothing
end

"""
    rk = upperTrapezoidal!(Gint, eInt, GcInt2, vActive, ieq)

Transform equations ieq..end to upper trapezoidal form and return the
rank `rk` (`ieq <= rk <= length(Gint)`).
"""
function upperTrapezoidal!(Gint, eInt::Vector{Int}, GcInt, GcInt2, vActive::Vector{Bool}, ieq::Int)::Int
    j        = 0
    e        = 0
    vj       = 0
    pivot    = 0
    oldPivot = 1
    neq      = length(Gint)
    buffer   = Buffer()
    for k = ieq:neq
        # Search for a pivot in equations k:neq
        j = 0
        e = 0

        # Inspect only equations with one variable
        for k2 = k:neq
            if length(Gint[k2]) == 1
                j = getFirstActiveIndex(Gint, vActive, k2)
                if j > 0
                    e = k2
                    break
                end
            end
        end

        if j == 0
            # Inspect only equations with two variables
            for k2 = k:neq
                if length(Gint[k2]) == 2
                    j = getFirstActiveIndex(Gint, vActive, k2)
                    if j > 0
                        e = k2
                        break
                    end
                end
            end
        end

        if j == 0
            # Inspect all equations
            for k2 = k:neq
                j = getFirstActiveIndex(Gint, vActive, k2)
                if j > 0
                    e = k2
                    break
                end
            end
        end

        if j > 0
            # Extract infos
            pivot = GcInt2[e][j]
            vj    = Gint[e][j]
            deleteat!(Gint[e] , j)  # Remove variable vj in equation e (is added later to the front)
            deleteat!(GcInt2[e], j)

            # Swap equations k and e
            if e != k
                swapEquations!(Gint, eInt, GcInt, GcInt2, k, e)
            end
        else
            # no pivot found: equations k:neq have no active variables (rank = k-1)
            return k-1
        end

        # Eliminate variable vj from equations k+1:neq
        for i = k+1:neq
            eliminateVariable!(Gint, GcInt2, k, i, pivot, vj, oldPivot, buffer)
        end
        oldPivot = pivot

        # Add variable vj in equation k at the front
        pushfirst!(Gint[k] , vj)
        pushfirst!(GcInt2[k], pivot)

        #println("... k = $k")
        #println("        Gint  = ", Gint)
        #println("        GcInt2 = ", GcInt2)
    end
    return neq
end

"""
    AvarRev = revertAssociation(Avar::Vector{Int})

Revert the association Vector `Avar[i] = j`, such that
`AvarRev[j] = i` (`Avar[i]=0` is allowed and is ignored).
"""
function revertAssociation(Avar::Vector{Int})::Vector{Int}
    AvarRev  = fill(0, length(Avar))
    for (i, value) in enumerate(Avar)
        if value != 0
            AvarRev[ value ] = i
        end
    end
    return AvarRev
end

needsVariableElimination(vProperty, eq_k::Vector{Int}) = findfirst(v->isEliminated(vProperty,v), eq_k) != nothing

"""
    equationRemoved = simplifyOneEquation!(eq_k, eqc_k, AvarRev, vEliminated, vProperty)

Simplify one equation has much as possible. Return true, if equation is removed, otherwise return false.
"""
function simplifyOneEquation!(eq_k, eqc_k, AvarRev, vEliminated, vProperty)::Bool
    while length(eq_k) > 1 && needsVariableElimination(vProperty, eq_k[2:end])
        # Equation has more as one variable and at least one of the variables is zero, alias or negative alias variable
        for j = 2:length(eq_k)
            v_j = eq_k[j]

            # Eliminate variable if possible
            if isNotEliminated(vProperty, v_j)
                continue

            elseif isZero(vProperty, v_j)
                deleteat!(eq_k , j)
                deleteat!(eqc_k, j)
                break

            else  # alias or negAlias
                if isAlias(vProperty, v_j)
                    v_add  = alias(vProperty, v_j)
                    vc_add = eqc_k[j]
                else
                    v_add  = negAlias(vProperty, v_j)
                    vc_add = -eqc_k[j]
                end

                # Check whether v_add appears in equation
                isPresent = false
                for i = 2:length(eq_k)
                    if i != j && eq_k[i] == v_add
                        isPresent = true
                        eqc_k[i] += vc_add
                        if eqc_k[i] == 0
                            if i < j
                                deleteat!(eq_k , [i,j])
                                deleteat!(eqc_k, [i,j])
                            else
                                deleteat!(eq_k , [j,i])
                                deleteat!(eqc_k, [j,i])
                            end
                        else
                            deleteat!(eq_k , j)
                            deleteat!(eqc_k, j)
                        end
                        break
                    end
                end
                if isPresent
                    break
                else
                    eq_k[j]  = v_add
                    eqc_k[j] = vc_add
                end
            end
        end
    end

    # Check if equation can be removed
    vk = eq_k[1]
    if AvarRev[vk] == 0
        # vk is not a derivative of a variable -> it can be removed
        if length(eq_k) == 1  # equation k is a function of one variable
            # Variable is zero -> remove equation
            push!(vEliminated, vk)
            vProperty[vk] = 0
            empty!(eq_k)
            empty!(eqc_k)
            return true

        elseif length(eq_k) == 2 && abs(eqc_k[1]) == abs(eqc_k[2]) # equation k is a function of alias variables
            if eqc_k[1] > 0 && eqc_k[2] < 0  ||
               eqc_k[1] < 0 && eqc_k[2] > 0
                # Alias variable -> remove equation
                push!(vEliminated, vk)
                vProperty[vk] = eq_k[2]
                empty!(eq_k)
                empty!(eqc_k)
            else
                # Negative alias variable -> remove equation
                push!(vEliminated, vk)
                vProperty[vk] = -eq_k[2]
                empty!(eq_k)
                empty!(eqc_k)
            end
            return true
        end
    end
    return false
end

function printLinearIntegerEquations(Gint, eInt, GcInt, var_name::Function; rk=(0,0,0))::Int
    ne = 0
    for i = 1:length(Gint)
        e = Gint[i]

        if length(e) > 0
            ne += 1
            # Construct equation
            vc     = GcInt[i][1]
            prefix = (vc == 1 ? "" : (vc == -1 ? "-" : string(vc) ))
            str    = "0 = " * prefix * var_name(e[1])
            for j = 2:length(e)
                v  = e[j]
                vc = GcInt[i][j]
                str = str *  (vc > 0 ? " + " : " - ")

                if abs(vc) != 1
                    str = str * string(abs(vc)) * "*"
                end
                str = str * var_name(v)
            end
            println(lpad(string(eInt[i]), 8), ": ", str)
        end
        if i == rk[1]
            println("    ---------- rk1")
        elseif i == rk[2]
            println("    ---------- rk2")
        elseif i == rk[3]
            println("    ---------- rk3")
        end
    end
    return ne
end

"""
    (vEliminated, vProperty, nvArbitrary, redundantEquations) =
        simplifyLinearIntegerEquations!(G, eInt, GcInt, Avar)

Remove singularities of the **linear Integer equations** of a DAE system and simplify these equations as much as possible.

The following **singularities** are fixed by this function:

- Equations that are redundant are removed.

- State constraint that are not structurally visible are transformed
  to a structurally visible form so that structural index reduction algorithms,
  such as the Pantelides algorithm, can handle these state constraints.

- Variables that can have an arbitrary value and do not appear in the remaining set of equations
  (so must be solved from the linear Integer equations) are set to zero.
  The calling function should print a warning message in such a case (if `nvArbitrary > 0`).

The following **simplifications** are performed recursively (`c` is an arbitrary Integer literal):

- An equation `c*v1 = 0` is removed and `v1` is replaced by zero in all occurrences
  of the linear Integer equations and all expressions with `v1` are simplified.

- An equation `c*v2 + c*v3 = 0` or `c*v2 - c*v3 = 0` is removed and
  `v2` is replaced by `v3` or by `-v3` in all occurrences of the linear Integer equations
  and all expressions with `v2` and `v3` are simplified.

Note:

- Input arguments `G, eInt, GcInt` are changed by a call of this function and
  represent the transformed linear Integer equations.
  The Abstract Syntax Tree (AST) of all these equations must be
  replaced by new AST equations defined by the returned arguments.

- The number of operations in the transformed linear Integer equations is
  guaranteed to be not larger as the original equations - with exception
  of the structurally visible state constraints, where the number of operations
  could be larger.

- Potential states (variables appearing differentiated) and derivatives
  of potential states are **not** eliminated by this function.

# Input arguments

- `G`: Bi-partite graph/incidence matrix of all equations. Typically: `G::Vector{Vector{Int}}`
       On entry, `G` is the original graph. On exit, the linear Integer equations of `G` are
       typically changed.

- `eInt::Vector{Int}`: `G[eInt]` are the linear Integer equations of `G`.
       On exit, `eInt` is reordered. If `length(G[eInt[i]]) = 0`, then the corresponding equation is eliminated.

- `GcInt`: `GcInt[i]` is the vector of Integer coefficients that are associated with the variables
           of `G[eInt[i]]`. Typically: `GcInt::Vector{Vector{Int}}`.
           On exit, `GcInt` is reordered according to `eInt` and typically most of the coefficients have been changed.

- `Avar::Vector{Int}`: Defines the derivatives of the variables:
           `A[i] = if der(v_i) == v_k then k else 0`. This vector is not changed by `simplifyLinearIntegerEquations!`.

# Output arguments

- `vEliminated::Vector{Int}`: Variables that are eliminated.

- `vProperty::Vector{Int}`:  Defines the properties of the eliminated variables.
  These properties can be inquired with the following exported functions:

  o `isNotEliminated(vProperty, v)` - if variable v is not eliminated.
  o `isEliminated(vProperty, v)` - if variable v is eliminated.
  o `isZero(vProperty,v)` - if eliminated variable v is zero.
  o `isAlias(vProperty,v)` - if eliminated variable v is an alias variable v = v_alias
  o `isNegAlias(vProperty,v)` - if eliminated variable v is a negative alias variable v = -v_alias
  o `alias(vProperty,v)` - alias variable v_alias of eliminated variable v (v = v_alias).
  o `negAlias(vProperty,v)` - negated alias variable v_alias of eliminated variable v (v = -v_alias).

- `nvArbitrary::Int`: Variables `vEliminated[1:nvArbitrary]` are variables that can be arbitrarily set and that
              have been set to zero.

- `redundantEquations::Vector{Int}`: G[redundantEquations] are redundant equations that have been removed.

# Algorithm

The algorithm to remove the singularities is sketched in the paper:

- Otter, Elmqvist (2017):
  [Transformation of Differential Algebraic Array Equations to Index One Form](http://www.ep.liu.se/ecp/132/064/ecp17132565.pdf),
  section 5. Modelica'2017 Conference.

An error in this algorithm was fixed, the algorithm was improved to handle large equation systems
and to simplify equations as much as possible.

# Main developer

[Martin Otter](https://rmc.dlr.de/sr/en/staff/martin.otter/),
[DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)

"""
function simplifyLinearIntegerEquations!(G, eInt::Vector{Int}, GcInt, Avar::Vector{Int}; log::Bool=false, var_name::Function = v -> "???")
    nv = length(Avar)

    # Revert derivative association vector Avar
    AvarRev = revertAssociation(Avar)

    # Construct vActive1 and vActive2 for first and second transformation to upper trapezoidal form
    # (vActiveX[v] = false, if variable shall be ignored when transforming to upper trapezoidal form)
    #     vActive1[v] = true, if variable is not a potential state (does not appear differentiated)
    #                         and variable v must be solved from the linear Integer equations
    #     vActive2[v] = true, if variable is not a potential state (does not appear differentiated)
    vActive2 = [v==0 for v in Avar]
    vActive1 = copy(vActive2)
    for e in setdiff(collect(1:length(G)), eInt) # Equations that are not linear Integer equations
        for v in G[e]
            vActive1[v] = false
        end
    end

    # Save all variables with vActive1[v] = true in vShouldBeSolved
    vShouldBeSolved = findall(vActive1)

    # Construct the bi-partite graph of the linear Integer equations
    # (a deepcopy of the relevant part of G)
    Gint = deepcopy(G[eInt])

    if log
        println("\n+++ Remove singularities")
        println("    Linear Integer equations:")
        printLinearIntegerEquations(Gint, eInt, GcInt, var_name)

        println("    Unknown variables:")
        unknowns = collect(OrderedCollections.OrderedSet([v for e in Gint for v in e]))
        for v in unknowns
            print(lpad(string(v),8), ": ", var_name(v))
            if v in vShouldBeSolved
                print("    (to be solved by equations)\n")
            elseif Avar[v] > 0
                print("    (potential state)\n")
            else
                print("\n")
            end
        end
    end

    # Construct a deepcopy of GcInt (this copy is modified by the transformation to upper trapezoidal form)
    GcInt2 = deepcopy(GcInt)

    # First transformation to upper trapezoidal form:
    #    Diagonal entries: Variables v that must be solved from the linear Integer equations
    #println("\n... before upperTrapezoidal!:")
    #println("     Gint  = $Gint")
    #println("     GcInt2 = $GcInt2")
    rk1 = upperTrapezoidal!(Gint, eInt, GcInt, GcInt2, vActive1, 1)

        # Eliminate variables that must be solved from the linear Integer equations
        # but are not diagonal entries (= variables can be arbitrarily set)
        vSolved     = [Gint[i][1] for i = 1:rk1]
        vEliminated = setdiff(vShouldBeSolved, vSolved)
        vProperty   = fill(IS_PRESENT, nv)
        for v in vEliminated
            vProperty[v] = 0
        end
        nvArbitrary = length(vEliminated)

        if log
            println("\n    After first transformation to trapezoidal form (eliminate variables that must be solved):")
            printLinearIntegerEquations(Gint, eInt, GcInt2, var_name, rk=(rk1,0,0))
        end

    # Second transformation to upper trapezoidal form: Ignore potential states
    rk2 = upperTrapezoidal!(Gint, eInt, GcInt, GcInt2, vActive2, rk1+1)

    if log
        println("\n    After second transformation to trapezoidal form (ignore potential states):")
        printLinearIntegerEquations(Gint, eInt, GcInt2, var_name, rk=(rk1,rk2,0))
    end

    # Third transformation to upper trapezoidal form: All remaining variables are potential states
    fill!(vActive2, true)
    rk3 = upperTrapezoidal!(Gint, eInt, GcInt, GcInt2, vActive2, rk2+1)
    #println("\n... after upperTrapezoidal!:")
    #println("     Gint  = $Gint")
    #println("     GcInt2 = $GcInt2")
    #println("     rk1 = $rk1, rk2 = $rk2, rk3 = $rk3")

    if log
        println("\n    After third transformation to trapezoidal form (eliminate potential states):")
        printLinearIntegerEquations(Gint, eInt, GcInt2, var_name, rk=(rk1, rk2, rk3))
    end

    # Simplify equations from equation rk2 upto equation 1
    for k = rk2:-1:1
        simplifyOneEquation!(Gint[k], GcInt2[k], AvarRev, vEliminated, vProperty)
    end

    if log
        println("\n    After alias elimination:")
        printLinearIntegerEquations(Gint, eInt, GcInt2, var_name, rk=(rk1, rk2, rk3))
    end

    #println("\n... after equation simplification:")
    #println("     Gint  = $Gint")
    #println("     GcInt2 = $GcInt2")
    #println("     vEliminated = ", vEliminated)
    #println("     vProperty[vEliminated] = ",  vProperty[vEliminated])

    # Update GcInt (use GcInt2[i] if it has not more unknowns as the corresponding GcInt[i] equation)
    equationsRemoved = false
    for i = 1:rk2
        if length(GcInt2[i]) <= length(GcInt[i])
            # Use transformed equation
            GcInt[i] = GcInt2[i]
        else
            # Use original equation
            Gint[i] = G[ eInt[i] ]
            equationRemoved = simplifyOneEquation!(Gint[i], GcInt[i], AvarRev, vEliminated, vProperty)
            equationsRemoved = equationsRemoved || equationRemoved
        end
    end

    if equationsRemoved
        # Simplify equations, until no equation is removed anymore
        equationRemoved = false
        while true
            for i = 1:rk2
                if length(Gint[i]) > 0
                    equationRemoved = simplifyOneEquation!(Gint[i], GcInt[i], AvarRev, vEliminated, vProperty)
                    if equationRemoved
                        break
                    end
                end
            end
            if !equationRemoved
                break
            end
        end
    end

    # For constraint equations and for removed equations, use transformed equations
    for i = rk2+1:length(eInt)
        GcInt[i] = GcInt2[i]
    end

    if log
        println("\n    Final, simplified equations:")
        printLinearIntegerEquations(Gint, eInt, GcInt2, var_name, rk=(rk1, rk2, rk3))
    end

    # Update G
    for (i,e) in enumerate(eInt)
        G[e] = Gint[i]
    end

    redundantEquations = rk3 < length(eInt) ? eInt[rk3+1:end] : Int[]

    return (vEliminated, vProperty, nvArbitrary, redundantEquations)
end

"""
    printSimplifiedLinearIntegerEquations(G, eInt, GcInt, vEliminated, vProperty,
        nvArbitrary, redundantEquations, var_name::Function; printTest=false)

Print result of [`simplifyLinearIntegerEquations!`](@ref).

Function `var_name(v)` returns the name of variable `v` as String.

If `printTest=true`, statements are printed that can be included in a Testset.
"""
function printSimplifiedLinearIntegerEquations(G, eInt, GcInt, vEliminated, vProperty, nvArbitrary, redundantEquations, var_name::Function; printTest=false)::Nothing
    if nvArbitrary > 0
        println("\n    Variables that can be arbitrarily set and have been set to zero:")
        for v in vEliminated[1:nvArbitrary]
            println(lpad(string(v), 8), ": ", var_name(v), " = 0")
        end
    end

    if length(vEliminated) - nvArbitrary > 0
        println("\n    Variables that have been eliminated:")
        for v in vEliminated[nvArbitrary+1:end]
            print(lpad(string(v), 8), ": ", var_name(v), " = ")
            if isZero(vProperty, v)
                println("0")
            elseif isAlias(vProperty, v)
                println(var_name(alias(vProperty,v)))
            else
                println("-", var_name(negAlias(vProperty,v)))
            end
        end
    end

    if length(redundantEquations) > 0
        println("\n    Redundant equations that have been removed:")
        for e in redundantEquations
            println("      ", e)
        end
    end

    println("\n    Remaining transformed linear Integer equations:")
    ne = printLinearIntegerEquations(G[eInt], eInt, GcInt, var_name)
    if ne == 0
        println("      none (all linear Integer equations are removed)")
    end
    println()

    if printTest
        println("@test nvArbitrary == ", nvArbitrary)
        println("@test vEliminated == ", vEliminated)
        println("@test vProperty[vEliminated] == ", vProperty[vEliminated])
        println("@test redundantEquations == ", redundantEquations)
        println("@test eInt  == ", eInt)
        println("@test G[eInt] == ", G[eInt])
        println("@test GcInt == ", GcInt)
    end

    return nothing
end


text:
function groupsort_indexer(x::AbstractVector, ngroups::Integer)
    # translated from Wes McKinney's groupsort_indexer in pandas (file: src/groupby.pyx).

    # count group sizes, location 0 for NA
    n = length(x)
    # counts = x.pool
    counts = fill(0, ngroups + 1)
    for i = 1:n
        counts[x[i] + 1] += 1
    end

    # mark the start of each contiguous group of like-indexed data
    where = fill(1, ngroups + 1)
    for i = 2:ngroups+1
        where[i] = where[i - 1] + counts[i - 1]
    end

    # this is our indexer
    result = fill(0, n)
    for i = 1:n
        label = x[i] + 1
        result[where[label]] = i
        where[label] += 1
    end
    result, where, counts
end

groupsort_indexer(pv::PooledDataVector) = groupsort_indexer(pv.refs, length(pv.pool))


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libpcap_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libpcap")
JLLWrappers.@generate_main_file("libpcap", UUID("f78539f5-9ba6-5ded-ae79-6d41194c59b7"))
end  # module libpcap_jll


text:
abstract type TypeOfUnknown end

struct RandomStateType <: TypeOfUnknown end

struct TupleType <: TypeOfUnknown
    elements :: Array{TypeOfUnknown, 1}
end

struct ArrayType <: TypeOfUnknown
    shape :: Nullable{Array{Nullable{Int64}, 1}}
    eltype :: Nullable{DataType}
end

struct OpaqueType <: TypeOfUnknown
end

array_typeof(arr::Number) = ArrayType([], eltype(arr))
array_typeof(arr::AbstractArray) = ArrayType(collect(size(arr)), eltype(arr))

abstract type ArrayExpr end

struct Unknown <: ArrayExpr
    typ :: TypeOfUnknown

    # Technically, this object does not need a name, but:
    # 1. It makes debugging easier
    # 2. Julia automatically caches objects, so, say, two `Unknown(OpaqueType())`
    #    created seaparetly would have the same object ID, which would ruin the whole
    #    expression graph processing.
    name :: Symbol

    Unknown(typ; tag=nothing) = new(typ, tag === nothing ? gensym() : gensym(tag))
end

struct LazyArray <: ArrayExpr
    value :: Any
end

struct CallNode <: ArrayExpr
    func :: Symbol
    args :: Array{ArrayExpr, 1}
    typ :: TypeOfUnknown
end

map_accum(func, state, aexpr) = func(state, aexpr)

function map_accum(func, state, aexpr::CallNode)
    new_args = []
    for arg in aexpr.args
        state, new_arg = func(state, arg)
        push!(new_args, new_arg)
    end
    state, CallNode(aexpr.func, new_args, aexpr.typ)
end

expr_type(aexpr::Unknown) = aexpr.typ
expr_type(aexpr::LazyArray) = array_typeof(aexpr.value)
expr_type(aexpr::CallNode) = aexpr.typ

is_leaf(aexpr::ArrayExpr) = true
is_leaf(aexpr::CallNode) = false

is_node(aexpr::ArrayExpr) = false
is_node(aexpr::CallNode) = true

children(aexpr::ArrayExpr) = []
children(aexpr::CallNode) = aexpr.args


text:
using PracticePackage
using Test
using TestSetExtensions


text:
# # Tutorial
#
# ## Installing everything
#
# To install the default CSS framework [Bulma](https://bulma.io/), simply type
# ```julia
# Pkg.add("InteractBulma");
# ```

# in the REPL.
#
# The basic behavior is as follows: Interact provides a series of widgets, each widgets has a primary observable that can be obtained with `observe(widget)` and adding listeners to that observable one can provide behavior. Let's see this in practice.
#
# ## Displaying a widget
using InteractBulma, WebIO
ui = button()
display(ui)

# Note that `display` works in a [Jupyter notebook](https://github.com/JuliaLang/IJulia.jl) or in [Atom/Juno IDE](https://github.com/JunoLab/Juno.jl).
# InteractBase can also be deployed in Jupyter Lab, but that requires installing an extension first:
cd(Pkg.dir("WebIO", "assets"))
;jupyter labextension install webio
;jupyter labextension enable webio/jupyterlab_entry
# To deploy the app as a standalone Electron window, one would use [Blink.jl](https://github.com/JunoLab/Blink.jl):
using Blink
w = Window()
body!(w, ui);
# The app can also be served in a webpage via [Mux.jl](https://github.com/JuliaWeb/Mux.jl):
using Mux
webio_serve(page("/", req -> ui), rand(8000:9000)) # serve on a random port
#
# ## Adding behavior
# For now this button doesn't do anything. This can be changed by adding callbacks to its primary observable:
o = observe(ui)
# Each observable holds a value and its value can be inspected with the `[]` syntax:
o[]
# In the case of the button, the observable represents the number of times it has been clicked: click on it and check the value again.
#
# To add some behavior to the widget we can use the `on` construct. `on` takes two arguments, a function and an observable. As soon as the observable is changed, the function is called with the latest value.
on(println, o)
# If you click again on the button you will see it printing the number of times it has been clicked so far.
#
# *Tip*: anonymous function are very useful in this programming paradigm. For example, if you want the button to say "Hello!" when pressed, you should use:
on(n -> println("Hello!"), o)
#
# *Tip n. 2*: using the `[]` syntax you can also set the value of the observable:
o[] = 33;
# To learn more about Observables, check out their documentation [here](https://juliagizmos.github.io/Observables.jl/latest/).
# ## What widgets are there?
#
# Once you have grasped this paradigm, you can play with any of the many widgets available:
filepicker() |> display # observable is the path of selected file
textbox("Write here") |> display # observable is the text typed in by the user
autocomplete(["Mary", "Jane", "Jack"]) |> display # as above, but you can autocomplete words
checkbox(label = "Check me!") |> display # observable is a boolean describing whether it's ticked
toggle(label = "I have read and agreed") |> display # same as a checkbox but styled differently
slider(1:100, label = "To what extent?", value = 33) |> display # Observable is the number selected

# As well as the option widgets, that allow to choose among options:

dropdown(["a", "b", "c"]) |> display # Observable is option selected
togglebuttons(["a", "b", "c"]) |> display # Observable is option selected
radiobuttons(["a", "b", "c"]) |> display # Observable is option selected

# The option widgets can also take as input a dictionary (ordered dictionary is preferrable, to avoid items getting scrambled), in which case the label displays the key while the observable stores the value:
using DataStructures
s = dropdown(OrderedDict("a" => "Value 1", "b" => "Value 2"))
display(s)
#-
observe(s)[]
#
# ## A simpler approach for simpler cases
#
# While the approach sketched above works for all sorts of situations, there is a specific marcro to simplify it in some specific case. If you want to update some result (maybe a plot) as a function of some parameters (discrete or continuous) simply write `@manipulate` before the `for` loop. Discrete parameters will be replaced by `togglebuttons` and continuous parameters by `slider`: the result will be updated as soon as you click on a button or move the slider:
#
width, height = 700, 300
colors = ["black", "gray", "silver", "maroon", "red", "olive", "yellow", "green", "lime", "teal", "aqua", "navy", "blue", "purple", "fuchsia"]
color(i) = colors[i%length(colors)+1]
ui = @manipulate for nsamples in 1:200,
        sample_step in slider(0.01:0.01:1.0, value=0.1, label="sample step"),
        phase in slider(0:0.1:2pi, value=0.0, label="phase"),
        radii in 0.1:0.1:60
    cxs_unscaled = [i*sample_step + phase for i in 1:nsamples]
    cys = sin.(cxs_unscaled) .* height/3 .+ height/2
    cxs = cxs_unscaled .* width/4pi
    dom"svg:svg[width=$width, height=$height]"(
        (dom"svg:circle[cx=$(cxs[i]), cy=$(cys[i]), r=$radii, fill=$(color(i))]"()
            for i in 1:nsamples)...
    )
end
# or, if you want a plot with some variables taking discrete values:
using Plots, DataStructures

x = y = 0:0.1:30

freqs = OrderedDict(zip(["pi/4", "π/2", "3π/4", "π"], [π/4, π/2, 3π/4, π]))

mp = @manipulate for freq1 in freqs, freq2 in slider(0.01:0.1:4π; label="freq2")
    y = @. sin(freq1*x) * sin(freq2*x)
    plot(x, y)
end

# ## Layout
#
# To create a full blown web-app, you should learn the layout tools that the CSS framework you are using provides. See for example the [columns](https://bulma.io/documentation/columns/) and [layout](https://bulma.io/documentation/layout/) section of the Bulma docs. You can use [WebIO](https://github.com/JuliaGizmos/WebIO.jl) to create from Julia the HTML required to create these layouts.
#
# However, this can be overwhelming at first (especially for users with no prior experience in web design). A simpler solution is [CSSUtil](https://github.com/JuliaGizmos/CSSUtil.jl), a package that provides some tools to create simple layouts.
using CSSUtil
loadbutton = filepicker()
hellobutton = button("Hello!")
goodbyebutton = button("Good bye!")
ui = vbox( # put things one on top of the other
    loadbutton,
    hbox( # put things one next to the other
        pad(1em, hellobutton), # to allow some white space around the widget
        pad(1em, goodbyebutton),
    )
)
display(ui)
#
# ## Update widgets as function of other widgets
#
# Sometimes the full structure of the GUI is not known in advance. For example, let's imagine we want to load a DataFrame and create a button per column. Not to make it completely trivial, as soon as a button is pressed, we want to plot a histogram of the corresponding column.
#
# *Important note*: this app needs to run in Blink, as the browser doesn't allow us to get access to the local path of a file.
#
# We start by adding a `filepicker` to choose the file, and only once we have a file we want to update the GUI. this can be done as follows:
loadbutton = filepicker()
columnbuttons = Observable{Any}(dom"div"())
# `columnbuttons` is the `div` object that will contain all the relevant buttons. it is an `Observable` as we want its value to change over time.
# To add behavior, we can use `map!`:
using CSV, DataFrames
data = Observable{Any}(DataFrame)
map!(CSV.read, data, observe(loadbutton))
#
# Now as soon as a file is uploaded, the `Observable` `data` gets updated with the correct value. Now, as soon as `data` is updated, we want to update our buttons.
using CSSUtil
function makebuttons(df)
    buttons = button.(names(df))
    dom"div"(hbox(buttons))
end

map!(makebuttons, columnbuttons, data)
# Note that `data` is already an `Observable`, so there's no need to do `observe(data)`, `observe` can only be applied on a widget.
# We are almost done, we only need to add a callback to the buttons. The cleanest way is to do it during button initialization, meaning during our `makebuttons` step:
using Plots
plt = Observable{Any}(plot()) # the container for our plot
function makebuttons(df)
    buttons = button.(string.(names(df)))
    for (btn, name) in zip(buttons, names(df))
        map!(t -> histogram(df[name]), plt, observe(btn))
    end
    dom"div"(hbox(buttons))
end
#
# To put it all together:
using CSV, DataFrames, InteractBulma, WebIO, Observables, Plots, CSSUtil
loadbutton = filepicker()
columnbuttons = Observable{Any}(dom"div"())
data = Observable{Any}(DataFrame)
plt = Observable{Any}(plot())
map!(CSV.read, data, observe(loadbutton))

function makebuttons(df)
    buttons = button.(string.(names(df)))
    for (btn, name) in zip(buttons, names(df))
        map!(t -> histogram(df[name]), plt, observe(btn))
    end
    dom"div"(hbox(buttons))
end

map!(makebuttons, columnbuttons, data)

ui = dom"div"(loadbutton, columnbuttons, plt)
#
# And now to serve it in Blink:
using Blink
w = Window()
body!(w, ui)


text:
## extended methods for computing covariance and scatter matrix

# auxiliary functions

function _symmetrize!(a::DenseMatrix)
    m, n = size(a)
    m == n || error("a must be a square matrix.")
    for j = 1:n
        @inbounds for i = j+1:n
            vl = a[i,j]
            vr = a[j,i]
            a[i,j] = a[j,i] = middle(vl, vr)
        end
    end
    return a
end

function _scalevars(x::DenseMatrix, s::DenseVector, dims::Int)
    dims == 1 ? Diagonal(s) * x :
    dims == 2 ? x * Diagonal(s) :
    error("dims should be either 1 or 2.")
end

## scatter matrix

_unscaled_covzm(x::DenseMatrix, dims::Colon)   = unscaled_covzm(x)
_unscaled_covzm(x::DenseMatrix, dims::Integer) = unscaled_covzm(x, dims)

_unscaled_covzm(x::DenseMatrix, wv::AbstractWeights, dims::Colon)   =
    _symmetrize!(unscaled_covzm(x, _scalevars(x, values(wv))))
_unscaled_covzm(x::DenseMatrix, wv::AbstractWeights, dims::Integer) =
    _symmetrize!(unscaled_covzm(x, _scalevars(x, values(wv), dims), dims))

"""
    scattermat(X, [wv::AbstractWeights]; mean=nothing, dims=1)

Compute the scatter matrix, which is an unnormalized covariance matrix.
A weighting vector `wv` can be specified to weight
the estimate.

# Arguments
* `mean=nothing`: a known mean value. `nothing` indicates that the mean is
  unknown, and the function will compute the mean. Specifying `mean=0` indicates
  that the data are centered and hence there's no need to subtract the mean.
* `dims=1`: the dimension along which the variables are organized.
  When `dims = 1`, the variables are considered columns with observations in rows;
  when `dims = 2`, variables are in rows with observations in columns.
"""
function scattermat end

"""
    cov(X, w::AbstractWeights, vardim=1; mean=nothing,  corrected=false)

Compute the weighted covariance matrix. Similar to `var` and `std` the biased covariance
matrix (`corrected=false`) is computed by multiplying `scattermat(X, w)` by
``\\frac{1}{\\sum{w}}`` to normalize. However, the unbiased covariance matrix
(`corrected=true`) is dependent on the type of weights used:
* `AnalyticWeights`: ``\\frac{1}{\\sum w - \\sum {w^2} / \\sum w}``
* `FrequencyWeights`: ``\\frac{1}{\\sum{w} - 1}``
* `ProbabilityWeights`: ``\\frac{n}{(n - 1) \\sum w}`` where ``n`` equals `count(!iszero, w)`
* `Weights`: `ArgumentError` (bias correction not supported)
"""
cov

"""
    mean_and_cov(x, [wv::AbstractWeights,] vardim=1; corrected=false) -> (mean, cov)

Return the mean and covariance matrix as a tuple. A weighting
vector `wv` can be specified. `vardim` that designates whether
the variables are columns in the matrix (`1`) or rows (`2`).
Finally, bias correction is applied to the covariance calculation if
`corrected=true`. See [`cov`](@ref) documentation for more details.
"""
function mean_and_cov end

scattermat(x::DenseMatrix; mean=nothing, dims::Int=1) =
    _scattermatm(x, mean, dims)
_scattermatm(x::DenseMatrix, ::Nothing, dims::Int) =
    _unscaled_covzm(x .- mean(x, dims=dims), dims)
_scattermatm(x::DenseMatrix, mean, dims::Int=1) =
    _unscaled_covzm(x .- mean, dims)

scattermat(x::DenseMatrix, wv::AbstractWeights; mean=nothing, dims::Int=1) =
    _scattermatm(x, wv, mean, dims)
_scattermatm(x::DenseMatrix, wv::AbstractWeights, ::Nothing, dims::Int) =
    _unscaled_covzm(x .- mean(x, wv, dims=dims), wv, dims)
_scattermatm(x::DenseMatrix, wv::AbstractWeights, mean, dims::Int) =
    _unscaled_covzm(x .- mean, wv, dims)

## weighted cov
covm(x::DenseMatrix, mean, w::AbstractWeights, dims::Int=1;
     corrected::DepBool=nothing) =
    rmul!(scattermat(x, w, mean=mean, dims=dims), varcorrection(w, depcheck(:covm, corrected)))

cov(x::DenseMatrix, w::AbstractWeights, dims::Int=1; corrected::DepBool=nothing) =
    covm(x, mean(x, w, dims=dims), w, dims; corrected=depcheck(:cov, corrected))

function corm(x::DenseMatrix, mean, w::AbstractWeights, vardim::Int=1)
    c = covm(x, mean, w, vardim; corrected=false)
    s = stdm(x, w, mean, vardim; corrected=false)
    cov2cor!(c, s)
end

"""
    cor(X, w::AbstractWeights, dims=1)

Compute the Pearson correlation matrix of `X` along the dimension
`dims` with a weighting `w` .
"""
cor(x::DenseMatrix, w::AbstractWeights, dims::Int=1) =
    corm(x, mean(x, w, dims=dims), w, dims)

function mean_and_cov(x::DenseMatrix, dims::Int=1; corrected::Bool=true)
    m = mean(x, dims=dims)
    return m, covm(x, m, dims, corrected=corrected)
end
function mean_and_cov(x::DenseMatrix, wv::AbstractWeights, dims::Int=1;
                      corrected::DepBool=nothing)
    m = mean(x, wv, dims=dims)
    return m, cov(x, wv, dims; corrected=depcheck(:mean_and_cov, corrected))
end

"""
    cov2cor(C, s)

Compute the correlation matrix from the covariance matrix `C` and a vector of standard
deviations `s`. Use `StatsBase.cov2cor!` for an in-place version.
"""
cov2cor(C::AbstractMatrix, s::AbstractArray) = cov2cor!(copy(C), s)

"""
    cor2cov(C, s)

Compute the covariance matrix from the correlation matrix `C` and a vector of standard
deviations `s`. Use `StatsBase.cor2cov!` for an in-place version.
"""
cor2cov(C::AbstractMatrix, s::AbstractArray) = cor2cov!(copy(C), s)

"""
    cor2cov!(C, s)

Converts the correlation matrix `C` to a covariance matrix in-place using a vector of
standard deviations `s`.
"""
function cor2cov!(C::AbstractMatrix, s::AbstractArray)
    n = length(s)
    size(C) == (n, n) || throw(DimensionMismatch("inconsistent dimensions"))
    for i in CartesianIndices(size(C))
        @inbounds C[i] *= s[i[1]] * s[i[2]]
    end
    return C
end


text:
using Test
using MDCDL
using Random
using LinearAlgebra

@testset "Orthonormal matrices" begin
    Random.seed!(883290234387)

    @testset "factorization" begin
        for p = 2:32
            for n = 1:5
                # generate random orthonormal matrix
                A = qr(rand(p, p)).Q |> Matrix

                θ, μ = MDCDL.mat2rotations(A)

                @test length(θ) == fld(p * (p-1), 2)
                @test length(μ) == p
                @test isreal(θ)
                @test isreal(μ)
                @test all(abs.(μ) .≈ 1)
            end
        end
    end

    @testset "generation" begin
        for p = 2:32
            for n = 1:5
                θ = pi * (rand(fld(p * (p-1), 2)) .- 0.5)
                μ = rand([-1.0, 1.0], p)

                A = MDCDL.rotations2mat(θ, μ)

                @test all(size(A) .== (p, p))
                @test isreal(A)
                @test A' * A ≈ I
            end
        end
    end

    @testset "reconstruction" begin
        for p = 2:32
            for n = 1:5
                A = qr(rand(p, p)).Q |> Matrix

                θ, μ = MDCDL.mat2rotations(A)
                rA = MDCDL.rotations2mat(θ, μ)

                @test rA ≈ A
            end
        end
    end

    @testset "non-destructivity" begin
        for p = 2:32
            for n = 1:5
                A = qr(rand(p, p)).Q |> Matrix

                expctdA = copy(A)
                θ, μ = MDCDL.mat2rotations(A)
                @test A == expctdA

                expctdθ = copy(θ)
                expctdμ = copy(μ)
                rA = MDCDL.rotations2mat(θ, μ)
                @test θ == expctdθ
                @test μ == expctdμ
            end
        end
    end
end

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
# Digital signal plots

function PlotDigital(x::AbstractArray{T,1}, y::AbstractArray{T,1}; label::String = "",
                     count::Integer = min(length(x),length(y)), offset::Integer = 0,
                     stride::Integer = 1) where {T}

    if eltype(x) == Float64
        LibCImPlot.PlotDigitaldoublePtrdoublePtr(label, x, y, Cint(count), Cint(offset),
                                                 Cint(stride * sizeof(Float64)))
    elseif eltype(x) == Float32
        LibCImPlot.PlotDigitalFloatPtrFloatPtr(label, x, y, Cint(count), Cint(offset),
                                               Cint(stride * sizeof(Float32)))
    else
        x = convert.(Float32, x)
        y = convert.(Float32, y)
        LibCImPlot.PlotDigitalFloatPtrFloatPtr(label, x, y, Cint(count), Cint(offset),
                                               Cint(stride * sizeof(Float32)))
    end
end

function PlotDigital(x::AbstractArray{T1,1}, y::AbstractArray{T2,1}; kwargs...) where {T1,T2}
    (x, y) = promote(x,y)
    PlotDigital(x,y; kwargs...)
end


text:
benchmarks = [
    Dict(
        :name => "Modified LV for testing",
        :ode => @ODEmodel(
            x1'(t) = (a + b) * x1(t) - c * x1(t) * x2(t),
            x2'(t) = -a * b * x2(t) + d * x1(t) * x2(t),
            y1(t) = x1(t)
        )
    ),
    Dict(
        :name => "SIWR original",
        :ode => @ODEmodel(
            S'(t) = mu - bi * S(t) * I(t) - bw * S(t) * W(t) - mu * S(t) + a * R(t),
            I'(t) = bw * S(t) * W(t) + bi * S(t) * I(t) - (gam + mu) * I(t),
            W'(t) = xi * (I(t) - W(t)),
            R'(t) = gam * I(t) - (mu + a) * R(t),
            y(t) = k * I(t)
        ),
        :skip => false
    ),
    Dict(
        :name => "SIWR with extra output",
        :ode => @ODEmodel(
            S'(t) = mu - bi * S(t) * I(t) - bw * S(t) * W(t) - mu * S(t) + a * R(t),
            I'(t) = bw * S(t) * W(t) + bi * S(t) * I(t) - (gam + mu) * I(t),
            W'(t) = xi * (I(t) - W(t)),
            R'(t) = gam * I(t) - (mu + a) * R(t),
            y(t) = k * I(t),
            y2(t) = S(t) + I(t) + R(t)
        ),
        :skip => false
    ),
    Dict(
        :name => "Pharm",
        :ode => @ODEmodel(
            x0'(t) = a1 * (x1(t) - x0(t)) - (ka * n * x0(t)) / (kc * ka + kc * x2(t) + ka * x0(t)),
            x1'(t) = a2 * (x0(t) - x1(t)),
            x2'(t) = b1 * (x3(t) - x2(t)) - (kc * n * x2(t)) / (kc * ka + kc * x2(t) + ka * x0(t)),
            x3'(t) = b2 * (x2(t) - x3(t)),
            y1(t) = x0(t)
        ),
        :skip => false
    ),
    Dict(
        :name => "SEAIJRC Covid model",
        :ode => @ODEmodel(
            S'(t) = -b * S(t) * (I(t) + J(t) + q * A(t)) * Ninv(t),
            E'(t) = b * S(t) * (I(t) + J(t) + q * A(t)) * Ninv(t) - k * E(t),
            A'(t) = k * (1 - r) * E(t) - g1 * A(t),
            I'(t) = k * r * E(t) - (alpha + g1) * I(t),
            J'(t) = alpha * I(t) - g2 * J(t),
            C'(t) = alpha * I(t),
            Ninv'(t) = 0,
            y(t) = C(t),
            y2(t) = Ninv(t)
        ),
        :skip => false
    ),
    Dict(
        :name => "MAPK model (5 outputs)",
        :ode => @ODEmodel(
            KS00'(t) = -a00 * K(t) * S00(t) + b00 * KS00(t) + gamma0100 * FS01(t) + gamma1000 * FS10(t) + gamma1100 * FS11(t),
            KS01'(t) = -a01 * K(t) * S01(t) + b01 * KS01(t) + c0001 * KS00(t) - alpha01 * F(t) * S01(t) + beta01 * FS01(t) + gamma1101 * FS11(t),
            KS10'(t) = -a10 * K(t) * S10(t) + b10 * KS10(t) + c0010 * KS00(t) - alpha10 * F(t) * S10(t) + beta10 * FS10(t) + gamma1110 * FS11(t),
            FS01'(t) = -alpha11 * F(t) * S11(t) + beta11 * FS11(t) + c0111 * KS01(t) + c1011 * KS10(t) + c0011 * KS00(t),
            FS10'(t) = a00 * K(t) * S00(t) - (b00 + c0001 + c0010 + c0011) * KS00(t),
            FS11'(t) = a01 * K(t) * S01(t) - (b01 + c0111) * KS01(t),
            K'(t) = a10 * K(t) * S10(t) - (b10 + c1011) * KS10(t),
            F'(t) = alpha01 * F(t) * S01(t) - (beta01 + gamma0100) * FS01(t),
            S00'(t) = alpha10 * F(t) * S10(t) - (beta10 + gamma1000) * FS10(t),
            S01'(t) = alpha11 * F(t) * S11(t) - (beta11 + gamma1101 + gamma1110 + gamma1100) * FS11(t),
            S10'(t) = -a00 * K(t) * S00(t) + (b00 + c0001 + c0010 + c0011) * KS00(t) - a01 * K(t) * S01(t) + (b01 + c0111) * KS01(t) - a10 * K(t) * S10(t) + (b10 + c1011) * KS10(t),
            S11'(t) = -alpha01 * F(t) * S01(t) + (beta01 + gamma0100) * FS01(t) - alpha10 * F(t) * S10(t) + (beta10 + gamma1000) * FS10(t) - alpha11 * F(t) * S11(t) + (beta11 + gamma1101 + gamma1110 + gamma1100) * FS11(t),
            y1(t) = F(t),
            y2(t) = S00(t),
            y3(t) = S01(t),
            y4(t) = S10(t),
            y5(t) = S11(t)
        )
    ),
    Dict(
        :name => "MAPK model (5 outputs bis)",
        :ode => @ODEmodel(
            KS00'(t) = -a00 * K(t) * S00(t) + b00 * KS00(t) + gamma0100 * FS01(t) + gamma1000 * FS10(t) + gamma1100 * FS11(t),
            KS01'(t) = -a01 * K(t) * S01(t) + b01 * KS01(t) + c0001 * KS00(t) - alpha01 * F(t) * S01(t) + beta01 * FS01(t) + gamma1101 * FS11(t),
            KS10'(t) = -a10 * K(t) * S10(t) + b10 * KS10(t) + c0010 * KS00(t) - alpha10 * F(t) * S10(t) + beta10 * FS10(t) + gamma1110 * FS11(t),
            FS01'(t) = -alpha11 * F(t) * S11(t) + beta11 * FS11(t) + c0111 * KS01(t) + c1011 * KS10(t) + c0011 * KS00(t),
            FS10'(t) = a00 * K(t) * S00(t) - (b00 + c0001 + c0010 + c0011) * KS00(t),
            FS11'(t) = a01 * K(t) * S01(t) - (b01 + c0111) * KS01(t),
            K'(t) = a10 * K(t) * S10(t) - (b10 + c1011) * KS10(t),
            F'(t) = alpha01 * F(t) * S01(t) - (beta01 + gamma0100) * FS01(t),
            S00'(t) = alpha10 * F(t) * S10(t) - (beta10 + gamma1000) * FS10(t),
            S01'(t) = alpha11 * F(t) * S11(t) - (beta11 + gamma1101 + gamma1110 + gamma1100) * FS11(t),
            S10'(t) = -a00 * K(t) * S00(t) + (b00 + c0001 + c0010 + c0011) * KS00(t) - a01 * K(t) * S01(t) + (b01 + c0111) * KS01(t) - a10 * K(t) * S10(t) + (b10 + c1011) * KS10(t),
            S11'(t) = -alpha01 * F(t) * S01(t) + (beta01 + gamma0100) * FS01(t) - alpha10 * F(t) * S10(t) + (beta10 + gamma1000) * FS10(t) - alpha11 * F(t) * S11(t) + (beta11 + gamma1101 + gamma1110 + gamma1100) * FS11(t),
            y0(t) = K(t),
            y1(t) = F(t),
            y2(t) = S00(t),
            y3(t) = S01(t) + S10(t),
            y4(t) = S11(t)
        ),
        :skip => false
    ),
    Dict(
        :name => "MAPK model (6 outputs)",
        :ode => @ODEmodel(
            KS00'(t) = -a00 * K(t) * S00(t) + b00 * KS00(t) + gamma0100 * FS01(t) + gamma1000 * FS10(t) + gamma1100 * FS11(t),
            KS01'(t) = -a01 * K(t) * S01(t) + b01 * KS01(t) + c0001 * KS00(t) - alpha01 * F(t) * S01(t) + beta01 * FS01(t) + gamma1101 * FS11(t),
            KS10'(t) = -a10 * K(t) * S10(t) + b10 * KS10(t) + c0010 * KS00(t) - alpha10 * F(t) * S10(t) + beta10 * FS10(t) + gamma1110 * FS11(t),
            FS01'(t) = -alpha11 * F(t) * S11(t) + beta11 * FS11(t) + c0111 * KS01(t) + c1011 * KS10(t) + c0011 * KS00(t),
            FS10'(t) = a00 * K(t) * S00(t) - (b00 + c0001 + c0010 + c0011) * KS00(t),
            FS11'(t) = a01 * K(t) * S01(t) - (b01 + c0111) * KS01(t),
            K'(t) = a10 * K(t) * S10(t) - (b10 + c1011) * KS10(t),
            F'(t) = alpha01 * F(t) * S01(t) - (beta01 + gamma0100) * FS01(t),
            S00'(t) = alpha10 * F(t) * S10(t) - (beta10 + gamma1000) * FS10(t),
            S01'(t) = alpha11 * F(t) * S11(t) - (beta11 + gamma1101 + gamma1110 + gamma1100) * FS11(t),
            S10'(t) = -a00 * K(t) * S00(t) + (b00 + c0001 + c0010 + c0011) * KS00(t) - a01 * K(t) * S01(t) + (b01 + c0111) * KS01(t) - a10 * K(t) * S10(t) + (b10 + c1011) * KS10(t),
            S11'(t) = -alpha01 * F(t) * S01(t) + (beta01 + gamma0100) * FS01(t) - alpha10 * F(t) * S10(t) + (beta10 + gamma1000) * FS10(t) - alpha11 * F(t) * S11(t) + (beta11 + gamma1101 + gamma1110 + gamma1100) * FS11(t),
            y0(t) = K(t),
            y1(t) = F(t),
            y2(t) = S00(t),
            y3(t) = S01(t),
            y4(t) = S10(t),
            y5(t) = S11(t)
        )
    ),
    Dict(
        :name => "Goodwin oscillator",
        :ode => @ODEmodel(
            x1'(t) = -b * x1(t) + 1 / (c + x4(t)),
            x2'(t) = alpha * x1(t) - beta * x2(t),
            x3'(t) = gama * x2(t) - delta * x3(t),
            x4'(t) = sigma * x4(t) * (gama * x2(t) - delta * x3(t)) / x3(t),
            y(t) = x1(t)
        )
    ),
    Dict(
        :name => "HIV",
        :ode => @ODEmodel(
            x'(t) = lm - d * x(t) - beta * x(t) * v(t),
            y'(t) = beta * x(t) * v(t) - a * y(t),
            v'(t) = k * y(t) - u * v(t),
            w'(t) = c * x(t) * y(t) * w(t) - c * q * y(t) * w(t) - b * w(t),
            z'(t) = c * q * y(t) * w(t) - h * z(t),
            y1(t) = w(t),
            y2(t) = z(t)
        )
    ),
    Dict(
        :name => "SIRS forced",
        :ode => @ODEmodel(
            s'(t) = mu - mu * s(t) - b0 * (1 + b1 * x1(t)) * i(t) * s(t) + g * r(t),
            i'(t) = b0 * (1 + b1 * x1(t)) * i(t) * s(t) - (nu + mu) * i(t),
            r'(t) = nu * i(t) - (mu + g) * r(t),
            x1'(t) = -M * x2(t),
            x2'(t) = M * x1(t),
            y1(t) = i(t),
            y2(t) = r(t)
        )
    ),
    Dict(
        :name => "Akt pathway",
        :ode => @ODEmodel(
            EGFR'(t) = EGFR_turnover * pro_EGFR(t) + EGF_EGFR(t) * reaction_1_k2 - EGFR(t) * EGFR_turnover - EGF_EGFR(t) * reaction_1_k1,
            pEGFR'(t) = EGF_EGFR(t) * reaction_9_k1 - pEGFR(t) * reaction_4_k1 + pEGFR_Akt(t) * reaction_2_k2 + pEGFR_Akt(t) * reaction_3_k1 - Akt(t) * pEGFR(t) * reaction_2_k1,
            pEGFR_Akt'(t) = Akt(t) * pEGFR(t) * reaction_2_k1 - pEGFR_Akt(t) * reaction_3_k1 - pEGFR_Akt(t) * reaction_2_k2,
            Akt'(t) = pAkt(t) * reaction_7_k1 + pEGFR_Akt(t) * reaction_2_k2 - Akt(t) * pEGFR(t) * reaction_2_k1,
            pAkt'(t) = pAkt_S6(t) * reaction_5_k2 - pAkt(t) * reaction_7_k1 + pAkt_S6(t) * reaction_6_k1 + pEGFR_Akt(t) * reaction_3_k1 - S6(t) * pAkt(t) * reaction_5_k1,
            S6'(t) = pAkt_S6(t) * reaction_5_k2 + pS6(t) * reaction_8_k1 - S6(t) * pAkt(t) * reaction_5_k1,
            pAkt_S6'(t) = S6(t) * pAkt(t) * reaction_5_k1 - pAkt_S6(t) * reaction_6_k1 - pAkt_S6(t) * reaction_5_k2,
            pS6'(t) = pAkt_S6(t) * reaction_6_k1 - pS6(t) * reaction_8_k1,
            EGF_EGFR'(t) = EGF_EGFR(t) * reaction_1_k1 - EGF_EGFR(t) * reaction_9_k1 - EGF_EGFR(t) * reaction_1_k2,
            y1(t) = a1 * (pEGFR(t) + pEGFR_Akt(t)),
            y2(t) = a2 * (pAkt(t) + pAkt_S6(t)),
            y3(t) = a3 * pS6(t)
        )
    ),
    Dict(
        :name => "CD8 T cell differentiation",
        :ode => @ODEmodel(
            N'(t) = - N(t) * mu_N - N(t) * P(t) * delta_NE,
            E'(t) = N(t) * P(t) * delta_NE - E(t)^2 * mu_EE - E(t) * delta_EL + E(t) * P(t) * rho_E,
            S'(t) = S(t) * delta_EL - S(t) * delta_LM - S(t)^2 * mu_LL - E(t) * S(t) * mu_LE,
            M'(t) = S(t) * delta_LM - mu_M * M(t),
            P'(t) = P(t)^2 * rho_P - P(t) * mu_P - E(t) * P(t) * mu_PE - S(t) * P(t) * mu_PL,
            y1(t) = N(t),
            y2(t) = E(t) + S(t),
            y3(t) = M(t)
        )
    ),
    Dict(
        :name => "Chemical reaction network",
        :ode => @ODEmodel(
            x1'(t) = -k1 * x1(t) * x2(t) + k2 * x4(t) + k4 * x6(t),
            x2'(t) = -k1 * x1(t) * x2(t) + k2 * x4(t) + k3 * x4(t),
            x3'(t) = k3 * x4(t) + k5 * x6(t) - k6 * x3(t) * x5(t),
            x4'(t) = k1 * x1(t) * x2(t) - k2 * x4(t) - k3 * x4(t),
            x5'(t) = k4 * x6(t) + k5 * x6(t) - k6 * x3(t) * x5(t),
            x6'(t) = -k4 * x6(t) - k5 * x6(t) + k6 * x3(t) * x5(t),
            y1(t) = x3(t),
            y2(t) = x2(t)
        )
    )
]

# the NFkB example

ode = @ODEmodel(
    x1'(t) = k_prod  - k_deg * x1(t) - k1 * x1(t) * u(t),
    x2'(t) = -k3 * x2(t) - k_deg * x2(t) - a2 * x2(t) * x10(t) + t1 * x4(t) - a3 * x2(t) * x13(t) + t2 * x5(t) + (k1 * x1(t) - k2 * x2(t) * x8(t)) * u(t),
    x3'(t) = k3 * x2(t) - k_deg * x3(t) + k2 * x2(t) * x8(t) * u(t),
    x4'(t) = a2 * x2(t) * x10(t) - t1 * x4(t),
    x5'(t) = a3 *  x2(t) * x13(t) - t2 *  x5(t),
    x6'(t) = c6a * x13(t) - a1 * x6(t) * x10(t) + t2 * x5(t) - i1 * x6(t),
    x7'(t) = i1 * kv * x6(t) - a1 * x11(t) *  x7(t),
    x8'(t) = c4 * x9(t) - c5 * x8(t),
    x9'(t) = c2 - c1 * x7(t) - c3 * x9(t),
    x10'(t) = -a2 * x2(t) * x10(t) - a1 * x10(t) * x6(t) + c4a * x12(t) - c5a * x10(t) - i1a * x10(t) + e1a * x11(t),
    x11'(t) = -a1 * x11(t) * x7(t) + i1a * kv * x10(t) - e1a * kv * x11(t),
    x12'(t) = c2a + c1a * x7(t) - c3a * x12(t),
    x13'(t) = a1 * x10(t) * x6(t) - c6a * x13(t) - a3 * x2(t) * x13(t) + e2a * x14(t),
    x14'(t) = a1 * x11(t) * x7(t) - e2a * kv * x14(t),
    x15'(t) = c2c + c1c * x7(t) - c3c * x15(t),
    y1(t) = x7(t),
    y2(t) = x10(t) + x13(t),
    y3(t) = x9(t),
    y4(t) = x1(t) + x2(t) + x3(t),
    y5(t) = x2(t),
    y6(t) = x12(t)
)

QQ = StructuralIdentifiability.Nemo.QQ

ode = set_parameter_values(ode, Dict(
    a1 => QQ(1, 2),
    a2 => QQ(1, 5),
    a3 => QQ(1),
    c1a => QQ(5, 10^(7)),
    c2a => QQ(0),
    c5a => QQ(1, 10^(4)),
    c6a => QQ(2, 10^(5)),
    c1 => QQ(5, 10^(7)),
    c2 => QQ(0),
    c3 => QQ(4, 10^(4)),
    c4 => QQ(1, 2),
    kv => QQ(5),
    e1a => QQ(5, 10^(4)),
    c1c => QQ(5, 10^(7)),
    c2c => QQ(0),
    c3c => QQ(4, 10^(4))
))

push!(
    benchmarks,
    Dict(
        :name => "NFkB",
        :ode => ode,
        :skip => true
    )
)


text:
module TestUtils

using AMGX, Test, Defer

@testset "Utils" begin
    @testset "api_version" begin
        v = AMGX.api_version()
        @test v isa VersionNumber
    end

    @testset "build_info" begin
        version, date, time = AMGX.build_info()
        @test version isa String
        @test date isa String
        @test time isa String
    end

    @testset "versioninfo" begin
        str = sprint(AMGX.versioninfo)
        @test occursin("AMGX version", str)
        @test occursin("Built on", str)
        @test occursin("API version", str)
    end

    @testset "error_string begin" begin
        @test AMGX.error_string(AMGX.API.AMGX_RC_BAD_PARAMETERS) == "Incorrect parameters for amgx call."
        @test AMGX.error_string(AMGX.API.AMGX_RC_OK) == "No error."
        @test AMGX.error_string(AMGX.API.AMGX_RC_INTERNAL) == "Internal error."
    end

    @scope @testset "print callback" begin
        result_print = ""
        AMGX.register_print_callback(x -> (result_print = x; nothing))
        c = @! AMGX.Config("")
        @test !isempty(result_print)
    end

    @testset "pin / unpin" begin
        v = rand(10);
        AMGX.pin_memory(v)
        AMGX.unpin_memory(v)
    end

    @testset "signal handler" begin
        AMGX.install_signal_handler()
        AMGX.reset_signal_handler()
    end
end

end # module


text:
# This doesn't work. Just needed to see if it did

file = open("interp_fun.txt")

txt  = readall(file)
println(txt)
f(hello,interp_me) = txt

print(f("poop", "pooppoop pooop"))

max_stars_repo_name: sjkelly/JuliaPlayground

max_issues_repo_name: sjkelly/JuliaPlayground

max_forks_repo_name: sjkelly/JuliaPlayground