text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#                            Chron1.0DistOnly.jl                                #
#                                                                               #
#     Illustrates the use of the Chron.jl package to estimate eruption          #
#  and/or deposition ages.                                                      #
#                                                                               #
#     This file uses code cells (denoted by "## ---"). To evaluate a cell in    #
#  the Julia REPL and move to the next cell, the default shortcut in Atom is    #
#  alt-shift-enter.                                                             #
#                                                                               #
#     You may have to adjust the path below which specifies the location of     #
#  the CSV data files for each sample, depending on what you want to run.       #
#                                                                               #
#   Last modified by C. Brenhin Keller 2020-05-25                               #
#                                                                               #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
## --- Load the Chron package - - - - - - - - - - - - - - - - - - - - - - - - -

    try
        using Chron
    catch
        using Pkg
        Pkg.add("Chron")
        using Chron
    end

    using Statistics, StatsBase, DelimitedFiles, SpecialFunctions
    using Plots; gr();

## --- Define sample properties - - - - - - - - - - - - - - - - - - - - - - - -

    # This example data is from Clyde et al. (2016) "Direct high-precision
    # U–Pb geochronology of the end-Cretaceous extinction and calibration of
    # Paleocene astronomical timescales" EPSL 452, 272–280.
    # doi: 10.1016/j.epsl.2016.07.041

    nSamples = 5 # The number of samples you have data for
    smpl = NewChronAgeData(nSamples)
    smpl.Name      =  ("KJ08-157", "KJ04-75", "KJ09-66", "KJ04-72", "KJ04-70",)
    smpl.Height   .=  [     -52.0,      44.0,      54.0,      82.0,      93.0,]
    smpl.Height_sigma .= [    3.0,       1.0,       3.0,       3.0,       3.0,]
    smpl.Age_Sidedness .= zeros(nSamples) # Sidedness (zeros by default: geochron constraints are two-sided). Use -1 for a maximum age and +1 for a minimum age, 0 for two-sided
    smpl.Path = "examples/DenverUPbExampleData/" # Where are the data files?
    smpl.inputSigmaLevel = 2 # i.e., are the data files 1-sigma or 2-sigma. Integer.
    smpl.Age_Unit = "Ma" # Unit of measurement for ages and errors in the data files
    smpl.Height_Unit = "cm" # Unit of measurement for Height and Height_sigma

    # IMPORTANT: smpl.Height must increase with increasing stratigraphic height
    # -- i.e., stratigraphically younger samples must be more positive. For this
    # reason, it is convenient to represent depths below surface as negative
    # numbers.

    # For each sample in smpl.Name, we must have a csv file at smpl.Path which
    # contains each individual mineral age and uncertainty. For instance,
    # examples/DenverUPbExampleData/KJ08-157.csv contains:
    #
    #   66.12,0.14
    #   66.115,0.048
    #   66.11,0.1
    #   66.11,0.17
    #   66.096,0.056
    #   66.088,0.081
    #   66.085,0.076
    #   66.073,0.084
    #   66.07,0.11
    #   66.055,0.043
    #   66.05,0.16
    #   65.97,0.12

## --- Bootstrap pre-eruptive distribution - - - - - - - - - - - - - - - - - - -

    # Bootstrap a KDE of the pre-eruptive (or pre-deposition) zircon distribution
    # shape from individual sample datafiles using a KDE of stacked sample data
    BootstrappedDistribution = BootstrapCrystDistributionKDE(smpl)
    h = plot(range(0,1,length=length(BootstrappedDistribution)), BootstrappedDistribution,
        label="Bootstrapped distribution", xlabel="Time (arbitrary units)", ylabel="Probability Density",
        fg_color_legend=:white)
    savefig(h, joinpath(smpl.Path,"BootstrappedDistribution.pdf"))
    display(h)
## --- Estimate the eruption age distributions for each sample  - - - - - - - -

    # Configure distribution model here
    distSteps = 5*10^5 # Number of steps to run in distribution MCMC
    distBurnin = distSteps÷2 # Number to discard

    # Choose the form of the prior distribution to use
    # Some pre-defined possiblilities include UniformDisribution,
    # TriangularDistribution, and MeltsVolcanicZirconDistribution
    # or you can define your own as with BootstrappedDistribution above
    dist = BootstrappedDistribution

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

    # Run MCMC to estimate saturation and eruption/deposition age distributions
    @time tMinDistMetropolis(smpl,distSteps,distBurnin,dist)

    # This will save rank-order and distribution plots, and print results to a
    # csv file -- you can find them in smpl.Path

## ---


text:
type wfinfo
    mediumtype::String
    component :: Array{String} # provide information of output data component
    mediumsize :: Array{Int64}
    ext :: Int64
    gridsize :: Array{Float64}
    timesamp :: Int64 # Total time sampling nT
    timepoint :: Array{Float64}
end

function readwfinfo(DataPath::String)
    fid = open(DataPath,"r")
    iflag = read(fid,Int64,1)
    if iflag[1] == 1
        mediumtype = "free surface"
    elseif iflag[1] == 2
        mediumtype = "unlimited"
    else error("Check wavefield file")
    end

    CpntDict = Dict(1=>"P",2=>"Vx",3=>"Vz",
                  4=>"Txx",5=>"Tzz",6=>"Txz")
    nCpnt = read(fid,Int64,1)
    NumCpnt = read(fid,Int64,nCpnt[1])
    component = []
    for i in NumCpnt
      if i == 23
          push!(component,CpntDict[2])
          push!(component,CpntDict[3])
      elseif i == 45
          push!(component,CpntDict[4])
          push!(component,CpntDict[5])
      else push!(component,CpntDict[i])
      end
    end
    mediumsize = read(fid,Int64,2)
    ext = read(fid,Int64,1)
    gridsize = read(fid,Float64,2)
    timesamp = read(fid,Int64,1)
    timepoint = read(fid,Float64,timesamp[1])

    return wfinfo(mediumtype, component, mediumsize,
    ext[1], gridsize, timesamp[1], timepoint)

end
function readwfdata(fid::IOStream,compnumber::Int64,timesamp::Int64,
    startpoint::Int64,length::Int64)
    seek(fid,8*(8+compnumber+timesamp+startpoint))
    requireddata = read(fid,Float64,length)
    return requireddata
end

function readwfdata(fid::IOStream,startpoint::Int64,length::Int64)
    seek(fid,8)
    compnumber = read(fid,Int64,1)
    seek(fid,8*(7+compnumber[1]))
    timesamp = read(fid,Int64,1)
    seek(fid,8*(8+compnumber[1]+timesamp[1]+startpoint))
    requireddata = read(fid,Float64,length)
    return(requireddata)
end

function readwfdata(fid::IOStream,length::Int64)
    requireddata = read(fid,Float64,length)
    return requireddata
end

function readwfdata(DataPath::String)
    fid = open(DataPath,"r")
    iflag = read(fid,Int64,1)
    nCpnt = read(fid,Int64,1)
    NumCpnt = read(fid,Int64,nCpnt[1])
    mediumsize = read(fid,Int64,2)
    ext = read(fid,Int64,1)
    gridsize = read(fid,Float64,2)
    timesamp = read(fid,Int64,1)
    timepoint = read(fid,Float64,timesamp[1])
    CpntDict = Dict(1=>mediumsize[1]*mediumsize[2],
                    2=>mediumsize[1]*(mediumsize[2]-1),
                    3=>(mediumsize[1]-1)*mediumsize[2],
                    4=>mediumsize[1]*mediumsize[2],
                    5=>mediumsize[1]*mediumsize[2],
                    6=>(mediumsize[1]-1)*(mediumsize[2]-1))
    totalCpnt = 0
    for Cpnt in NumCpnt
        if Cpnt == 23
            totalCpnt = totalCpnt + CpntDict[2] + CpntDict[3]
        elseif Cpnt == 45
            totalCpnt = totalCpnt + CpntDict[4] + CpntDict[5]
        else totalCpnt = totalCpnt + CpntDict[Cpnt]
        end
    end
    requireddata = read(fid,Float64,totalCpnt*timesamp[1])
    requireddata = reshape(requireddata, totalCpnt, timesamp[1])
    return requireddata
end

ext: jl


text:
module CatPlanarGraphs

include("planar_graph.jl")

end


text:
using SDC
using Test

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


text:
# ------------------------------------------------------------------------------------------
# # Intro to Flux.jl
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We have learned how machine learning allows us to classify data as apples or bananas with
# a single neuron. However, some of those details are pretty fiddly! Fortunately, Julia has
# a powerful package that does much of the heavy lifting for us, called
# [`Flux.jl`](https://fluxml.github.io/).
#
# *Using `Flux` will make classifying data and images much easier!*
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Using `Flux.jl`
#
# In the next notebook, we are going to see how Flux allows us to redo the calculations from
# the previous notebook in a simpler way. We can get started with `Flux.jl` via:
# ------------------------------------------------------------------------------------------

using Flux

# ------------------------------------------------------------------------------------------
# #### Helpful built-in functions
#
# When working we'll `Flux`, we'll make use of built-in functionality that we've had to
# create for ourselves in previous notebooks.
#
# For example, the sigmoid function, σ, that we have been using already lives within `Flux`:
# ------------------------------------------------------------------------------------------

?σ

# ------------------------------------------------------------------------------------------
# Importantly, `Flux` allows us to *automatically create neurons* with the **`Dense`**
# function. For example, in the last notebook, we were looking at a neuron with 2 inputs and
# 1 output:
#
#  <img src="data/single-neuron.png" alt="Drawing" style="width: 500px;"/>
#
#  We could create a neuron with two inputs and one output via
# ------------------------------------------------------------------------------------------

model = Dense(2, 1, σ)

# ------------------------------------------------------------------------------------------
# This `model` object comes with places to store weights and biases:
# ------------------------------------------------------------------------------------------

model.W

model.b

typeof(model.W)

# ------------------------------------------------------------------------------------------
# Unlike in previous notebooks, note that `W` is no longer a `Vector` (1D `Array`) and `b`
# is no longer a number! Both are now stored in so-called `TrackedArray`s and `W` is
# effectively being treated as a matrix with a single row. We'll see why in the next
# notebook.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Other helpful built-in functionality includes ways to automatically calculate gradients
# and also the cost function that we've used in previous notebooks -
#
# $$L(w, b) = \sum_i \left[y_i - f(x_i, w, b) \right]^2$$
#
# This is the "mean square error" function, which in `Flux` is named **`Flux.mse`**.
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/13. ML - Intro to Flux.jl.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
# Implement A_mul_B! , evalT and evalP for AbstractMatrixLike
abstract type AbstractMatrixLike end

function Base.:*(A::AbstractMatrixLike,x::AbstractVector)
    b = similar(x)
    return A_mul_B!(b,A,x)
end

struct MimeticGrad{T,GradType,VecType,BcondType,NBF,NF} <: AbstractMatrixLike
    g::Vector{GradType} # cell stuff
    brcf::Vector{VecType} # cell-to-face vector
    rcf::Vector{NTuple{2,VecType}} # cell-to-face vector
    bvkinv::Vector{T} # 1/(cell_volume*K)
    vkinv::Vector{NTuple{2,T}} # 1/(cell_volume*K)
    bcond::BcondType
    f2cloop::Face2CellLoop{NBF,NF,VecType}
end

function MimeticGrad(k,bcond,m::HomogeneousMesh)
    T = eltype(k)
    vectype = vec_type(m)
    gradtype = grad_type(bcond)
    f2cloop = m.f2cloops
    g = zeros(gradtype,length(m.cells))
    brcf, rcf, bvkinv, vkinv = get_stuff_needed_for_MimeticGrad(k,m)
    return MimeticGrad{T,gradtype,vectype,typeof(bcond),nbfaces(m),nfaces(m)}(g,brcf,rcf,bvkinv,vkinv,bcond,f2cloop)
end

function get_stuff_needed_for_MimeticGrad(k,m::HomogeneousMesh)
    nbf = nbfaces(m)
    nf = nfaces(m)

    brcf = zeros(vec_type(m),nbf)
    rcf = Array{NTuple{2,vec_type(m)}}(nf)
    T = eltype(k)
    bvkinv = zeros(T,nbf)
    vkinv = Array{NTuple{2,T}}(nf)

    nodes = m.nodes
    cells = m.cells
    bf2c = m.f2cloops.bf2c
    bfaces = m.bfaces
    bcv = m.f2cloops.bcv
    @inbounds for i in 1:nbf
        j = bf2c[i][1]
        ccenter = center(cells[j],nodes)
        brcf[i] = center(bfaces[i],nodes) - ccenter
        bvkinv[i] = bcv[i]/k[i]
    end

    f2c = m.f2cloops.f2c
    faces = m.faces
    cv = m.f2cloops.cv
    @inbounds for i in 1:nf
        j1,j2 = f2c[i]
        cc1 = center(cells[j1],nodes)
        cc2 = center(cells[j2],nodes)
        fc = center(faces[i],nodes)
        rcf[i] = (fc - cc1, fc - cc2)
        cv1,cv2 = cv[i]
        vkinv[i] = (cv1/k[j1], cv2/k[j2])
    end

    return brcf, rcf, bvkinv, vkinv
end

ibfaces(::Type{MimeticGrad{T,GT,VecType,NC,NBF,NF}}) where {T,GT,VecType,NC,NBF,NF} = Base.OneTo(NBF)
@inline ibfaces(a::MimeticGrad) = ibfaces(typeof(a))

ifaces(::Type{MimeticGrad{T,GT,VecType,NC,NBF,NF}}) where {T,GT,VecType,NC,NBF,NF} = Base.OneTo(NF)
@inline ifaces(a::MimeticGrad) = ifaces(typeof(a))

function evalT(A::MimeticGrad,out::FaceVector{T},inp::FaceVector{T}) where {T}
    bf2c = A.f2cloop.bf2c
    brcf = A.brcf
    g = A.g
    bvkinv = A.bvkinv
    fill!(g,zero(eltype(g)))
    bcond = A.bcond
    bt = A.f2cloop.bft
    @inbounds for i in ibfaces(A)
        j = bf2c[i][1]
        if typeof(bcond[bt[i]]) <: Neumman
            inp[:b,i] = zero(T)
        else
            g[j] += inp[:b,i]*brcf[i]*bvkinv[i]
        end
    end

    f2c = A.f2cloop.f2c
    rcf = A.rcf
    vkinv = A.vkinv
    @inbounds for i in ifaces(A)
        j1,j2 = f2c[i]
        rc1,rc2 = rcf[i]
        v1,v2 = vkinv[i]
        g[j1] += inp[i]*rc1*v1
        g[j2] -= inp[i]*rc2*v2
    end

    δ = zero(Float64)
    @inbounds for i in ibfaces(A)
        j = bf2c[i][1]
        if typeof(bcond[bt[i]]) <: Neumman
            out[:b,i] = zero(T)
        else
            out[:b,i] = g[j] ⋅ brcf[i]
        end
        δ += inp[:b,i] ⋅ out[:b,i]
    end

    @inbounds for i in ifaces(A)
        j1,j2 = f2c[i]
        rc1,rc2 = rcf[i]
        out[i] = (g[j1] ⋅ rc1) - (g[j2] ⋅ rc2)
        δ += inp[i] ⋅ out[i]
    end

    return δ
end

function evalP(A::MimeticGrad,out::FaceVector{T},inp::FaceVector{T}) where {T}
    brcf = A.brcf
    bvkinv = A.bvkinv
    bcond = A.bcond
    bt = A.f2cloop.bft

    δ = zero(Float64)
    @inbounds for i in ibfaces(A)
        if typeof(bcond[bt[i]]) <: Neumman
            out[:b,i] = zero(T)
            inp[:b,i] = zero(T)
        else
            out[:b,i] = inp[:b,i]/((brcf[i]⋅brcf[i])*bvkinv[i])
        end
        δ += out[:b,i] ⋅  inp[:b,i]
    end

    rcf = A.rcf
    vkinv = A.vkinv
    @inbounds for i in ifaces(A)
        rc1,rc2 = rcf[i]
        v1,v2 = vkinv[i]
        out[i] = inp[i]/((rc1⋅rc1)*v1 + (rc2⋅rc2)*v2)
        δ += out[i] ⋅ inp[i]
    end

    return δ
end

function Base.A_mul_B!(out::FaceVector{T},A::MimeticGrad,inp::FaceVector{T}) where {T}
    bf2c = A.f2cloop.bf2c
    brcf = A.brcf
    g = A.g
    bvkinv = A.bvkinv
    fill!(g,zero(eltype(g)))

    @inbounds for i in ibfaces(A)
        j = bf2c[i][1]
        g[j] += inp[:b,i]*brcf[i]*bvkinv[i]
    end

    f2c = A.f2cloop.f2c
    rcf = A.rcf
    vkinv = A.vkinv
    @inbounds for i in ifaces(A)
        j1,j2 = f2c[i]
        rc1,rc2 = rcf[i]
        v1,v2 = vkinv[i]
        g[j1] += inp[i]*rc1*v1
        g[j2] -= inp[i]*rc2*v2
    end

    @inbounds for i in ibfaces(A)
        j = bf2c[i][1]
        out[:b,i] = g[j] ⋅ brcf[i]
    end

    @inbounds for i in ifaces(A)
        j1,j2 = f2c[i]
        rc1,rc2 = rcf[i]
        out[i] = (g[j1] ⋅ rc1) - (g[j2] ⋅ rc2)
    end

    return out
end

struct ImplicitMethodB{MeshType,BcondType,Ktype,T} <: AbstractMatrixLike
    m::MeshType
    bcond::BcondType
    k::Ktype
    dtrhoc::T
end

function evalT(A::ImplicitMethodB,out::AbstractArray{T},inp::AbstractArray{T}) where {T}

    fill!(out,zero(T))

    floops = A.m.f2cloops
    bf2c = floops.bf2c
    bcond = A.bcond
    bt = floops.bft
    bfc = floops.bfc
    bfn = floops.bfn
    k = A.k
    bccenter = floops.bccenter
    bcv = floops.bcv

    NBF = nbfaces(floops)
    @inbounds for i=1:NBF
        j = bf2c[i][1]
        bch = bcond[bt[i]]
        if typeof(bch) <: Neumman
            qf = 0.0
        else
            AfoLac = (bfn[i]⋅bfn[i])/(bfn[i]⋅(bfc[i] - bccenter[i]))
            qf = -k[j]*inp[j]*AfoLac
        end
        out[j] += qf*bcv[i]
    end

    f2c = floops.f2c
    fn = floops.fn
    ccenter = floops.ccenter
    cv = floops.cv
    AoL = floops.AoL

    NF = nfaces(floops)
    @inbounds for i=1:NF
        el = f2c[i]
        j1 = el[1]
        j2 = el[2]
        To = inp[j1]
        Ta = inp[j2]
        AfoLac = AoL[i]#(n⋅n)/(n⋅(c[2] - c[1]))
        qf = k[i]*(Ta-To)*AfoLac

        v = cv[i]
        out[j1] += qf*v[1]
        out[j2] -= qf*v[2]
    end

    δ = zero(T)
    dtrhoc = A.dtrhoc
    @inbounds for i in linearindices(out)
        out[i] = inp[i] - dtrhoc*out[i]
        δ += out[i]*inp[i]
    end

    return δ
end

function evalP(A::ImplicitMethodB,out::AbstractArray{T},inp::AbstractArray{T}) where {T}

    fill!(out,zero(T))

    floops = A.m.f2cloops
    bf2c = floops.bf2c
    bcond = A.bcond
    bt = floops.bft
    bfc = floops.bfc
    bfn = floops.bfn
    k = A.k
    bccenter = floops.bccenter
    bcv = floops.bcv

    NBF = nbfaces(floops)
    @inbounds for i=1:NBF
        j = bf2c[i][1]
        bch = bcond[bt[i]]
        if typeof(bch) <: Neumman
            qf = 0.0
        else
            AfoLac = (bfn[i]⋅bfn[i])/(bfn[i]⋅(bfc[i] - bccenter[i]))
            qf = -k[j]*AfoLac
        end
        out[j] += qf*bcv[i]
    end

    f2c = floops.f2c
    fn = floops.fn
    ccenter = floops.ccenter
    cv = floops.cv
    AoL = floops.AoL

    NF = nfaces(floops)
    @inbounds for i=1:NF
        el = f2c[i]
        j1 = el[1]
        j2 = el[2]
        AfoLac = AoL[i]
        #qf = k[i]*(Ta-To)*AfoLac

        v = cv[i]
        out[j1] -= k[j1]*AfoLac*v[1]
        out[j2] -= k[j2]*AfoLac*v[2]
    end

    δ = zero(T)
    dtrhoc = A.dtrhoc
    @inbounds for i in linearindices(out)
        out[i] = inp[i]/(1 - dtrhoc*out[i])
        δ += out[i]*inp[i]
    end

    return δ
end

function Base.A_mul_B!(out::AbstractArray{T},A::ImplicitMethodB,inp::AbstractArray{T}) where {T}
    floops = A.m.f2cloops
    bf2c = floops.bf2c
    bcond = A.bcond
    bt = floops.bft
    bfc = floops.bfc
    bfn = floops.bfn
    k = A.k
    bccenter = floops.bccenter
    bcv = floops.bcv

    fill!(out,zero(T))

    NBF = nbfaces(floops)
     @inbounds for i=1:NBF
         j = bf2c[i][1]
         qf = bpart(bcond[bt[i]],i,j,bfc,bfn,k,bccenter,inp)
         out[j] += qf*bcv[i]
     end

    f2c = floops.f2c
    fn = floops.fn
    ccenter = floops.ccenter
    cv = floops.cv

    NF = nfaces(floops)
    @inbounds for i=1:NF
        el = f2c[i]
        j1 = el[1]
        j2 = el[2]
        To = inp[j1]
        Ta = inp[j2]
        n = fn[i]
        c = ccenter[i]
        AfoLac = (n⋅n)/(n⋅(c[2] - c[1]))
        qf = k[i]*(Ta-To)*AfoLac

        v = cv[i]
        out[j1] += qf*v[1]
        out[j2] -= qf*v[2]
    end

    dtrhoc = A.dtrhoc
    @inbounds for i in linearindices(out)
        out[i] = inp[i] - dtrhoc*out[i]
    end

    return nothing
end

# Solve for a * I() + ∇ ⋅ ( b * ∇())
struct aIpDbG{MeshType,BType} <: AbstractMatrixLike # (aI() + Div(b*Grad()))()
    m::MeshType
    a::ConstVec{Float64}
    boundary_indices::BType
    bface_loop_value::Vector{Float64}
    face_loop_value::Vector{Tuple{Float64,Float64}}
    diag::Vector{Float64}
end

function aIpDbG(a,k,d,m)
    #a = ConstVec(d[:density]/d[:dt])
    #k = -0.5*d[:viscosity]
    boudary, bfl = set_boundary_things(k,d,m)
    fl = set_face_loop(k,d,m)
    diag = set_diag(a,k,d,m)
    return aIpDbG{typeof(m),typeof(boudary)}(m,a,boudary,bfl,fl,diag)
end

function set_boundary_things(k,d,m)
    floops = m.f2cloops
    bf2c = floops.bf2c
    bcond = uboundary_conditions(d)
    bt = floops.bft
    bfc = floops.bfc
    bfn = floops.bfn
    bccenter = floops.bccenter
    bcv = floops.bcv

    bfl = Vector{Float64}()
    boundary = ()
    NBF = nbfaces(floops)

    @inbounds for i=1:NBF
        j = bf2c[i][1]
        bch = bcond[bt[i]]
        if typeof(bch) <: Dirichlet
            AfoLac = (bfn[i]⋅bfn[i])/(bfn[i]⋅(bfc[i] - bccenter[i]))
            push!(bfl,-k*bcv[i]*AfoLac)
            boundary = (boundary..., j)
        end
    end
    return boundary, bfl
end

function set_face_loop(k,d,m)
    floops = m.f2cloops
    fn = floops.fn
    f2c = floops.f2c
    ccenter = floops.ccenter
    cv = floops.cv
    NF = nfaces(floops)
    fl = Vector{Tuple{Float64,Float64}}(NF)
    @inbounds for i=1:NF
        el = f2c[i]
        j1 = el[1]
        j2 = el[2]
        n = fn[i]
        c = ccenter[i]
        AfoLac = (n⋅n)/(n⋅(c[2] - c[1]))
        qf = k*AfoLac

        v = cv[i]
        fl[i] = (qf*v[1],qf*v[2])
    end
    return fl
end

function set_diag(a,k,d,m)
    out = zeros(length(m.cells))

    floops = m.f2cloops
    bf2c = floops.bf2c
    bcond = uboundary_conditions(d)
    bt = floops.bft
    bfc = floops.bfc
    bfn = floops.bfn
    bccenter = floops.bccenter
    bcv = floops.bcv

    NBF = nbfaces(floops)
    @inbounds for i=1:NBF
        j = bf2c[i][1]
        bch = bcond[bt[i]]
        if typeof(bch) <: Neumman
            qf = 0.0
        else
            AfoLac = (bfn[i]⋅bfn[i])/(bfn[i]⋅(bfc[i] - bccenter[i]))
            qf = -k*AfoLac
        end
        out[j] += qf*bcv[i]
    end

    f2c = floops.f2c
    fn = floops.fn
    ccenter = floops.ccenter
    cv = floops.cv

    NF = nfaces(floops)
    @inbounds for i=1:NF
        el = f2c[i]
        j1 = el[1]
        j2 = el[2]
        n = fn[i]
        c = ccenter[i]
        AfoLac = (n⋅n)/(n⋅(c[2] - c[1]))

        v = cv[i]
        out[j1] -= k*AfoLac*v[1]
        out[j2] -= k*AfoLac*v[2]
    end

    @inbounds for i in linearindices(out)
        out[i] = 1/(a[i] + out[i])
    end
    return out
end

function evalT(A::aIpDbG,out::AbstractArray{T},inp::AbstractArray{T}) where {T}

    fill!(out,zero(T))
    #bface_loop_value = -b[j] * AfoLac[i] * bcv[i]
    bface_loop_value = A.bface_loop_value
    @inbounds for (i,j) in enumerate(A.boundary_indices)
        out[j] += inp[j]*bface_loop_value[i]
    end

    floops = A.m.f2cloops
    f2c = floops.f2c
    NF = nfaces(floops)
    face_loop_value = A.face_loop_value
    @inbounds for i=1:NF
        el = f2c[i]
        j1 = el[1]
        j2 = el[2]
        To = inp[j1]
        Ta = inp[j2]
        ΔT = Ta - To

        f = face_loop_value[i] #face_loop_value = k[i]*AfoLac .* (cv[i])
        out[j1] += f[1]*ΔT
        out[j2] -= f[2]*ΔT
    end

    δ = zero(Float64)
    a = A.a
    @inbounds for i in linearindices(out)
        out[i] = a[i]*inp[i] + out[i]
        δ += out[i]⋅inp[i]
    end

    return δ
end

function evalP(A::aIpDbG,out::AbstractArray{T},inp::AbstractArray{T}) where {T}
    diag = A.diag
    δ = zero(Float64)
    @inbounds for i in linearindices(out)
        out[i] = inp[i]*diag[i]
        δ += out[i]⋅inp[i]
    end
    return δ
end

struct PoissonP{MeshType} <: AbstractMatrixLike # (Div(b*Grad()))()
    m::MeshType
    diag::Vector{Float64}
    k::ConstVec{Float64}
    g::Vector{Vec2D{Float64}}
end

function PoissonP(d,m)
    diag = set_diag_poisson(d,m)
    k = ConstVec(d[:dt]/d[:density])
    g = Vector{Vec2D{Float64}}(length(m.cells))
    return PoissonP{typeof(m)}(m,diag,k,g)
end

function set_diag_poisson(d,m)
    out = zeros(length(m.cells))
    g = Vec2DArray{Float64}(length(m.cells))
    floops = m.f2cloops
    bf2c = floops.bf2c
    bt = floops.bft
    bfc = floops.bfc
    bfn = floops.bfn
    k = d[:dt]/d[:density]
    bccenter = floops.bccenter
    bcv = floops.bcv

    f2c = floops.f2c
    fn = floops.fn
    ccenter = floops.ccenter
    cv = floops.cv

    NBF = FiniteVolumeMesh.nbfaces(floops)
    @inbounds for i=1:NBF
        j = bf2c[i][1]
        flux = 0.5 * bfn[i] * bcv[i]
        g[j] += flux
    end

    NF = FiniteVolumeMesh.nfaces(floops)
    @inbounds for i=1:NF
        j1,j2 = f2c[i]
        flux = 0.5 * fn[i]
        v = cv[i]
        g[j1] += flux*v[1]
        g[j2] += flux*v[2]
    end

    @inbounds for i in linearindices(out)
        out[i] = k * (g[i]⋅g[i])
    end

    NBF = FiniteVolumeMesh.nbfaces(floops)
    @inbounds for i=1:NBF
        j = bf2c[i][1]
        flux = 0.25 * k * (bfn[i]⋅bfn[i]) * bcv[i]
        out[j] -= flux
    end

    NF = FiniteVolumeMesh.nfaces(floops)
    @inbounds for i=1:NF
        j1,j2 = f2c[i]
        flux = 0.25 * k * (fn[i]⋅fn[i])
        v = cv[i]
        out[j1] -= flux*v[1]
        out[j2] -= flux*v[2]
    end

    @inbounds for i in linearindices(out)
        out[i] = 1/out[i]
    end

    return out
end

function pressure_gradient!(out::Vector, f, f2c,fn,cv,bf2c,bfn,bcv,floop::FiniteVolumeMesh.Face2CellLoop)

    fill!(out,zero(eltype(out)))
    NBF = FiniteVolumeMesh.nbfaces(floop)
    @inbounds for i=1:NBF
        j = bf2c[i][1]
        flux = bfn[i] * f[j] * bcv[i]
        out[j] += flux
    end

    NF = FiniteVolumeMesh.nfaces(floop)

    @inbounds for i=1:NF
        j1,j2 = f2c[i]
        flux = fn[i] * ((f[j1] + f[j2])*0.5)
        v = cv[i]
        out[j1] += flux*v[1]
        out[j2] -= flux*v[2]
    end
end

@inline pressure_gradient!(out, f, floop) = pressure_gradient!(out, f, floop.f2c,floop.fn,floop.cv,floop.bf2c,floop.bfn,floop.bcv,floop)

function pressure_div!(out::Vector, f, a,f2c,fn,cv,bf2c,bfn,bcv,floop::FiniteVolumeMesh.Face2CellLoop)

    fill!(out,zero(eltype(out)))
    NF = FiniteVolumeMesh.nfaces(floop)

    @inbounds for i=1:NF
        j1,j2 = f2c[i]
        flux = a * (fn[i] ⋅ ((f[j1] + f[j2])*0.5))
        v = cv[i]
        out[j1] += flux*v[1]
        out[j2] -= flux*v[2]
    end
end

@inline pressure_div!(out, f, a,floop) = pressure_div!(out, f, a,floop.f2c,floop.fn,floop.cv,floop.bf2c,floop.bfn,floop.bcv,floop)

function evalT(A::PoissonP,out::AbstractArray{T},inp::AbstractArray{T}) where {T}

    fill!(out,zero(T))
    g = A.g
    pressure_gradient!(g,inp,A.m.f2cloops)
    pressure_div!(out,g,A.k[1],A.m.f2cloops)

    δ = zero(Float64)

    @inbounds for i in linearindices(out)
        δ += out[i]⋅inp[i]
    end

    return δ

end

function evalP(A::PoissonP,out::AbstractArray{T},inp::AbstractArray{T}) where {T}

    diag = A.diag

    δ = zero(Float64)

    @inbounds for i in linearindices(out)
        out[i] = inp[i]*diag[i]
        δ += out[i]⋅inp[i]
    end

    return δ
end


text:
const IDENTITY_VA = gensym()

@testset "Basics" begin
    mock(IDENTITY_VA, identity) do id
        identity(10)
        @test called_once_with(id, 10)
    end

    mock(IDENTITY_VA, identity) do id
        identity(1, 2, 3)
        identity()
        @test called(id)
        @test called_with(id, 1, 2, 3)
        @test !called_once(id)
        @test ncalls(id) == 2
        @test has_call(id, Call(1, 2, 3))
        @test has_calls(id, Call(1, 2, 3), Call())
    end
end

@testset "Non-Mock mocks" begin
    @test mock(_id -> identity(2) == 4, identity => x -> 2x)
end

@testset "Specific methods" begin
    mock((+, Float64, Int) => Mock((a, b) -> 2a + b)) do plus
        @test 1 + 1 == 2
        @test 2.0 + 1 == 5.0
        @test called_once_with(plus, 2.0, 1)
    end
end

@testset "Multiple mocks" begin
    mock(*, /) do mul, div
        @test 1 * 1 != 1
        @test 1 / 1 != 1
        @test called_once_with(mul, 1, 1)
        @test called_once_with(div, 1, 1)
    end
end

@testset "Varargs" begin
    varargs(::Int, ::Int, ::String, ::String, ::String, ::Bool...) = true
    varargs(::Any) = false

    mock((varargs, Vararg{Int, 2}, Vararg{String, 3}, Vararg{Bool})) do va
        @test varargs(0, 0, "", "", "") !== true
        @test varargs(0, 0, "", "", "", false, false) !== true
        @test !varargs(0)
        @test ncalls(va) == 2
        @test has_calls(va, Call(0, 0, "", "", ""), Call(0, 0, "", "", "", false, false))
    end
end

@testset "Parametric types" begin
    params(::Vector{Int}) = 1
    params(::Vector{Bool}) = 2
    params(::Vector{<:AbstractString}) = 3
    params(::Vector{T}) where T <: Number = 4

    mock((params, Vector{Int})) do p
        @test params([1]) != 1
        @test params([true]) == 2
        @test called_once_with(p, [1])
    end

    mock((params, Vector{<:AbstractString})) do p
        @test params([""]) != 3
        @test params([strip("")]) != 3
        @test params([1.0]) == 4
        @test ncalls(p) == 2
        @test has_calls(p, Call([""]), Call([strip("")]))
    end

    mock((params, Vector{<:Number})) do p
        @test params([""]) == 3
        @test params([1.0]) != 4
        @test called_once_with(p, [1.0])
    end
end

@testset "mock does not overwrite methods" begin
    # https://github.com/fredrikekre/jlpkg/blob/3b1c2400932dbe13fa7c3cba92bde3842315976c/src/cli.jl#L151-L160
    o = JLOptions()
    if o.warn_overwrite == 0
        args = map(n -> n === :warn_overwrite ? 1 : getfield(o, n), fieldnames(JLOptions))
        unsafe_store!(cglobal(:jl_options, JLOptions), JLOptions(args...))
    end
    ctx = gensym()
    mock(identity, ctx, identity)
    out = @capture_err mock(identity, ctx, identity)
    @test isempty(out)
end

@testset "Reusing Context" begin
    f(x) = strip(uppercase(x))
    # If the method checks aren't working properly, this will throw.
    ctx = gensym()
    @test mock(_f -> f(" hi "), ctx, strip => identity) == " HI "
    @test mock(_f -> f(" hi "), ctx, uppercase => identity) == "hi"
end

@testset "Any context name is valid" begin
    # This will throw if the context names aren't independent of the user-supplied names.
    @test mock(_id -> true, :mock, identity)
end

@testset "Filters" begin
    @testset "Maximum/minimum depth" begin
        f(x) = identity(x)
        g(x) = f(x)
        h(x) = g(x)

        mock(IDENTITY_VA, identity; filters=[max_depth(3)]) do id
            @test f(1) != 1
            @test g(2) != 2
            @test h(3) == 3
            @test ncalls(id) == 2 && has_calls(id, Call(1), Call(2))
        end

        mock(IDENTITY_VA, identity; filters=[min_depth(3)]) do id
            @test f(1) == 1
            @test g(2) != 2
            @test h(3) != 3
            @test ncalls(id) == 2 && has_calls(id, Call(2), Call(3))
        end
    end

    @testset "Exclude/include" begin
        @eval module Bar
        a(x) = identity(x)
        b(x) = identity(x)
        end
        c(x) = identity(x)
        d(x) = identity(x)

        mock(IDENTITY_VA, identity; filters=[excluding(Bar, c)]) do id
            @test Bar.a(1) == 1
            @test Bar.b(2) == 2
            @test c(3) == 3
            @test d(4) != 4
        end

        mock(IDENTITY_VA, identity; filters=[including(Bar, c)]) do id
            @test Bar.a(1) != 1
            @test Bar.b(2) != 2
            @test c(3) != 3
            @test d(4) == 4
        end
    end
end

@testset "Keyword arguments" begin
    foo(; kwargs...) = get(kwargs, :foo, nothing)
    bar(; kwargs...) = foo(; kwargs...)
    baz(; kwargs...) = bar(; kwargs...)

    @testset "Keyword arguments are passed to mocked functions" begin
        mock(foo) do f
            @test bar(; foo=:bar) !== :bar
            @test called_once_with(f; foo=:bar)
        end
    end

    @testset "Keyword arguments are discarded when recursing" begin
        ctx = gensym()

        mock(ctx, bar; filters=[including()]) do _b
            @test_logs (:warn, "Discarding keyword arguments") baz(; foo=:baz)
            @test @suppress baz(; foo=:baz) === nothing
        end

        mock(ctx, foo) do f
            @test_logs (:warn, "Discarding keyword arguments") baz(; foo=:baz)
            result = @suppress baz(; foo=:baz)
            @test result !== nothing && result !== :baz
            @test ncalls(f) == 2 && has_calls(f, Call(), Call())
        end
    end

    @testset "Metadata records original functions and not keyword wrappers" begin
        @eval begin
            kw_f(; x) = kw_g(; x=x)
            kw_g(; x) = 2x
            kw_h = (; x) -> kw_i(; x=x)
            kw_i = (; x) -> 2x
        end
        kw_j(; x) = kw_k(; x=x)
        kw_k(; x) = 2x

        test(f; broken::Bool=false) = function(m::SimpleMock.Metadata)
            eq = SimpleMock.current_function(m) === f
            if broken
                @test_broken eq
            else
                @test eq
            end
            return true
        end

        # Regular functions.
        @test mock(_g -> kw_f(; x=1), kw_g; filters=[test(kw_f)]) != 2
        # Anonymous functions.
        @test mock(_i -> kw_h(; x=1), kw_i; filters=[test(kw_h)]) != 2
        # Closures.
        @test mock(_k -> kw_j(; x=1), kw_k; filters=[test(kw_j; broken=true)]) != 2
    end
end


text:
# p1.jl - convergence of fourth-order finite differences

# For various N, set up grid in [-pi,pi] and function u(x):
Nvec = 2.^(3:12);
clf(); axes([.1,.4,.8,.5]);
for N = Nvec
    h = 2*pi/N; x = -pi + (1:N)*h;
    u = @. exp(sin(x)^2);
    uprime = @. 2*sin(x)*cos(x)*u;

    # Construct sparse fourth-order differentiation matrix:
    e = ones(N);
    D = sparse(1:N,[2:N;1],2*e/3) - sparse(1:N,[3:N;1:2],e/12);
    D = (D-D')/h;

    # Plot max(abs(D*u-uprime)):
    error = norm(D*u-uprime,Inf);
    loglog(N,error,"k.",markersize=6);
end
grid(true); xlabel("N"); ylabel("error");
title("Convergence of fourth-order finite differences");
loglog(Nvec,1.0./Nvec.^4,"--");
text(105,5e-8,L"N^{-4}",fontsize=18);


text:
##################################################
# Random net and data draw from nets
##################################################

nets = [MCBN.random_net_weights(x) for x=[5,10,40]]
graphs = map(MCBN.net2graph, nets)
draws = map(x->MCBN.draw_data(x,10), nets)
#klds = Float64[]

global g1,g2,d1,d2,bns1,bns2,bns1c,bns2c
for i=1:100
    g1 = MCBN.random_net_weights(15)
    g2 = MCBN.random_net_weights(15)

    d1 = MCBN.draw_data(g1,50)
    d2 = MCBN.draw_data(g2,50)

    bns1 = MCBN.BayesNetSampler(g1,d1)
    bns1c = MCBN.BayesNetSampler(g1,d2)

    bns2 = MCBN.BayesNetSampler(g2,d2)
    bns2c = MCBN.BayesNetSampler(g2,d1)

    diff1 = MCBN.energy(bns1) - MCBN.energy(bns1c)
    diff2 = MCBN.energy(bns2) - MCBN.energy(bns2c)
    @test (diff1 + diff2) < 0
    # We do it this way because sometimes both networks will actually favor the
    # same dataset, but usually the 'wrong' favoritism is much smaller
    # magnitude than the correct one. So summing them up will make this
    # apparent (and prevent failures)

    #push!(klds, MCBN.kld(g1,g2))
    i%10==0 && print("\rSuccessfully completed round $i of energy checks")
end
println("")

D = 4
template = zeros(Int,D,D)
cpds = MCBN.BayesNetDAI(D)

function factor_klds(bnd1::MCBN.BayesNetDAI, bnd2::MCBN.BayesNetDAI)
    0 #FIXME
end

function prior(bns::MCBN.BayesNetSampler)
    p_cpds=0.2
    p_structural=0.2
    limparent=4
    d = size(bns.mat,1)
    tot = 0
    for i=1:d
        for j=1:d
            tot += bns.mat[i,j]
        end
        if tot > bns.limparent + 1
            return 1e10
        end
        tot = 0
    end

    adj_diff = sum(abs(bns.mat[bns.x,bns.x] - template))
    e_struct = p_structural * adj_diff

    # energy contribution from cpds
    cpd_diff = factor_klds(bns.bnd, cpds)
    e_cpd = p_cpds * cpd_diff
    return e_cpd + e_struct
end

gold = MCBN.random_net_weights(D)
data = MCBN.draw_data(gold,10)
samc2 = SAMC.MHRecord(MCBN.BayesNetSampler(D, data, prior), burn = 100, thin=20)
MCBN.check_bns(bns)
for i=1:1000
  SAMC.sample!(samc2, 10, verbose=0)
  MCBN.check_bns(bns)
end
println("Passed 1000 check_bns points every 10 samples")

# TODO: Test factor_klds


text:
#function logsumexp(x)
#    max_x = maximum(x)
#    exp_x = exp.(x .- max_x)
#    log(sum(exp_x)) .+ max_x
#end

function mbar_equation!(F, f_k, K, N_max, N_k, u_kln, idx, log_wi_jn, log_term, FlogN)
    for k = 1:K
        FlogN[k] = log(Float64(N_k[k])) + f_k[k]
    end
    for i = 1:K
        for j = 1:K
            for n = 1:N_k[j]
                ui_xjn = u_kln[j, i, n];
                for k = 1:K
                    uk_xjn = u_kln[j, k, n];
                    log_term[k] = FlogN[k] - (uk_xjn - ui_xjn);
                end
                log_wi_jn[j, n] = -logsumexp(log_term);
            end
        end
        F[i] = f_k[i] + logsumexp(log_wi_jn[idx])
    end
    F[:] = F[:] .- F[1]
end

"""
    mbar(u_kl; tol=1e-8, iterations_self=10) -> F

Estimates the free energy differences between different thermodynamnic states by using the Multistate Bennet Acceptance Ratio Method (MBAR).
Let K be # of thermodynamic states. `u_kl` is a K x K Array whose elements are reduced bias-factors or
potential energies of state k evaluated at state l. Solvs the MBAR equations by applying self-consistent iterations `iterations_self` times,
followed by the Newton-Raphson iterations. The iteration stops when the difference from the previous falls below the threshold `tol`.

Returns dimensionless free energies `f_k` whose elements corresponds to thermodynamics state k.

# Example
```julia-repl
see Juypter notebooks in MDToolbox.jl/examples/
```

# References
```
M. R. Shirts and J. D. Chodera, J. Chem. Phys. 129, 124105 (2008).
```
"""
function mbar(u_kl; tol=1e-8, iterations_self=10)
    # K: number of umbrella windows
    K, L = size(u_kl)

    # N_k: number of data in k-th umbrella window
    N_k = zeros(Int64, K)
    for k = 1:K
        N_k[k] = length(u_kl[k, 1])
    end
    N_max = maximum(N_k)

    # conversion from array of array (u_kl) to array (u_kln)
    u_kln = zeros(Float64, K, K, N_max)
    for k = 1:K
        for l = 1:K
            u_kln[k, l, 1:N_k[k]] .= u_kl[k, l]
        end
    end

    dummy_kn = zeros(Int64, K, N_max)
    for k = 1:K
      dummy_kn[k, 1:N_k[k]] .= 1
    end
    idx = dummy_kn .> 0;

    # solve the MBAR equation
    #log_wi_jn = zeros(Float64, K, N_max)
    #log_term = zeros(Float64, K)
    #FlogN = zeros(Float64, K)
    ## f!(F, x) = mbar_equation!(F, x, K, N_max, N_k, u_kln, idx, log_wi_jn, log_term, FlogN)
    ## x_init = zeros(Float64, K)
    ## sol = nlsolve(f!, x_init, autodiff=:forward, ftol=tol, iterations=iterations)
    ## #sol = nlsolve(f!, x_init, method = :anderson, ftol=tol, iterations=iterations)
    ## f_k = sol.zero
    ## f_k .- f_k[1]

    #F = zeros(Float64, K)
    #x_init = zeros(Float64, K)
    #for n = 1:1000
    #    mbar_equation!(F, x_init, K, N_max, N_k, u_kln, idx, log_wi_jn, log_term, FlogN)
    #    x_init = x_init - F
    #    x_init = x_init .- x_init[1]
    #end
    #f_k = x_init

    # solve the MBAR equation by self-consistent iterations
    f_k = zeros(Float64, K)
    f_k_new = zeros(Float64, K)

    log_wi_jn = zeros(Float64, (K, N_max))
    for j = 1:K
      log_wi_jn[j, 1:N_k[j]] .= 1
    end
    index = log_wi_jn .> 0.5

    for count_iteration = 1:iterations_self
        for i = 1:K
            log_wi_jn = mbar_log_wi_jn(N_k, f_k, u_kln, u_kln[:, i, :], K, N_max)
            f_k_new[i] = - logsumexp_1d(log_wi_jn[index])
        end

        f_k_new .= f_k_new .- f_k_new[1]
        check_convergence = maximum(abs.(f_k_new .- f_k)) ./ std(f_k_new)
        f_k .= f_k_new

        @printf "iteration = %d  delta = %e  tol = %e\n" count_iteration check_convergence tol
        @printf "free energies = "
        for k = 1:K
            @printf " %f" f_k[k]
        end
        @printf "\n"
        @printf "\n"
    end

    first_gamma = 0.1;
    gamma = 1.0;

    check_convergence = Inf64
    W_nk = zeros(Float64, sum(N_k), K)

    count_iteration = iterations_self
    while check_convergence > tol
        f_k_new .= f_k
        for i = 1:K
          log_wi_jn = mbar_log_wi_jn(N_k, f_k, u_kln, u_kln[:, i, :], K, N_max)
          W_nk[:, i] .= exp.(log_wi_jn[index] .+ f_k[i])
        end

        g = zeros(Float64, K-1, 1) # gradient
        H = zeros(Float64, K-1, K-1) # hessian
        for i = 1:(K-1)
            g[i] = N_k[i+1] - N_k[i+1] * sum(W_nk[:, i+1])
            H[i, i] = - sum(N_k[i+1] .* W_nk[:, i+1] .* (1.0 .- N_k[i+1] .* W_nk[:, i+1]))
            for j = 1:(i-1)
                H[i, j] = sum((N_k[i+1] .* W_nk[:, i+1]) .* (N_k[j+1] .* W_nk[:, j+1]))
                H[j, i] = H[i, j]
            end
        end

        Hinvg = pinv(H, atol=1.0e-10) * g
        for k = 1:(K-1)
            if count_iteration == iterations_self
                f_k_new[k+1] = f_k_new[k+1] - first_gamma * Hinvg[k]
            else
                f_k_new[k+1] = f_k_new[k+1] - gamma * Hinvg[k]
            end
        end

        check_convergence = maximum(abs.(f_k_new .- f_k)) ./ std(f_k_new)
        f_k .= f_k_new

        count_iteration = count_iteration + 1
        @printf "iteration =%d  delta = %e  tol = %e\n" count_iteration check_convergence tol
        @printf "free energies = "
        for k = 1:K
            @printf " %f" f_k[k]
        end
        @printf "\n"
        @printf "\n"
    end

    return f_k
end

"""
    mbar_weight(u_kl, f_k, u_k=nothing)

Calculate weights for computing expectations at the target thermodynamics state (ensemble)
by using the Multistate Bennet Acceptance Ratio Method (MBAR).
Let K be # of thermodynamic states. `u_kl` is a K x K Array whose elements are reduced bias-factors or
potential energies of state k evaluated at state l. `f_k` is dimension dimensionless free energies `f_k`
obtained by prior MBAR calculation. If `u_k` (reduced potential energies of the target state) is given,
weights for that state are computed. By default `u_k = nothing`, in thise case zero potential
energies are used for computing weights (useful for umbrella sampling data).

Returns weights `w_k` computed at the target thermodynamic state.

# Example
```julia-repl
see Juypter notebooks in MDToolbox.jl/examples/
```

# References
```
M. R. Shirts and J. D. Chodera, J. Chem. Phys. 129, 124105 (2008).
```
"""
function mbar_weight(u_kl, f_k, u_k=nothing)
    # K: number of umbrella windows
    K, L = size(u_kl)

    # N_k: number of data in k-th umbrella window
    N_k = zeros(Int64, K)
    for k = 1:K
        N_k[k] = length(u_kl[k, 1])
    end
    N_max = maximum(N_k)

    # conversion from array of array (u_kl) to array (u_kln)
    u_kln = zeros(Float64, K, K, N_max)
    for k = 1:K
        for l = 1:K
            u_kln[k, l, 1:N_k[k]] .= u_kl[k, l]
        end
    end

    # conversion from cell (u_k) to array (u_kn)
    u_kn = zeros(Float64, K, N_max)
    for k = 1:K
        if u_k === nothing
            u_kn[1, 1:N_k[k]] .= zero(Float64)
        else
            u_kn[k, 1:N_k[k]] .= u_k[k]
        end
    end

    log_w_kn = zeros(Float64, K, N_max)
    for k = 1:K
      log_w_kn[k, 1:N_k[k]] .= 1.0
    end
    idx = log_w_kn .> 0.5;

    log_w_kn = mbar_log_wi_jn(N_k, f_k, u_kln, u_kn, K, N_max)
    log_w_n  = log_w_kn[idx]

    s = logsumexp_1d(log_w_n)
    w_k = Vector{Vector{Float64}}(undef, K)
    for k = 1:K
      w_k[k] = exp.((log_w_kn[k, 1:N_k[k]] .- s))
    end

    return w_k
end

function ChainRulesCore.rrule(::typeof(mbar_weight), u_kl, f_k, u_k)
    w_k = mbar_weight(u_kl, f_k, u_k)
    function mbar_weight_pullback(dw_k)
        du_k = deepcopy(w_k)
        for k = 1:length(w_k)
            for n = 1:length(w_k[k])
                du_k[k][n] = dw_k[k][n] * (- w_k[k][n])
            end
        end
        return NoTangent(), ZeroTangent(), NoTangent(), du_k
    end
    return w_k, mbar_weight_pullback
end

# MATLAB-style coding
function mbar_log_wi_jn(N_k, f_k, u_kln, u_kn, K, N_max)
    log_wi_jn = zeros(Float64, (K, N_max))
    for k = 1:K
        x = repeat(log.(N_k), 1, N_k[k]) .+ repeat(f_k, 1, N_k[k]) .- (u_kln[k, :, 1:N_k[k]] .- repeat(u_kn[k:k, 1:N_k[k]], K, 1))
        log_wi_jn[k:k, 1:N_k[k]] .= - logsumexp_over_row(x)
    end
    return log_wi_jn
end

# FORTRAN-style coding
function mbar_log_wi_jn2(N_k, f_k, u_kln, u_kn, K, N_max)
    log_wi_jn = zeros(Float64, (K, N_max))
    FlogN = zeros(Float64, K)
    for k = 1:K
        FlogN[k] = log(N_k[k]) + f_k[k]
    end

    log_term = zeros(Float64, K)
    #Threads.@threads for k = 1:K
    for k = 1:K
        for n = 1:N_k[k]
            max_log_term = -Inf64

            u_k = u_kn[k, n]
            for l = 1:K
                u_l = u_kln[k, l, n]
                log_term[l] = FlogN[l] - (u_l - u_k)
                if log_term[l] > max_log_term
                    max_log_term = log_term[l]
                end
            end

            term_sum = 0.0
            for l = 1:K
                term_sum += exp(log_term[l] - max_log_term)
            end
            log_sum = log(term_sum) + max_log_term
            log_wi_jn[k, n] = - log_sum
        end
    end
    return log_wi_jn
end

function logsumexp_1d(x)
    max_x = maximum(x)
    exp_x = exp.(x .- max_x)
    s = log(sum(exp_x)) + max_x
    return s
end

function logsumexp_over_row(x)
    max_x = maximum(x, dims=1)
    exp_x = exp.(x .- max_x)
    s = log.(sum(exp_x, dims=1)) .+ max_x
    return s
end


text:
struct FrequencyBin
	f :: Float64
	df :: Float64
end

@with_kw mutable struct Trajectory{T}
	t :: Vector{T} = zeros(T, 0)
	f :: Vector{Float64} = zeros(Float64, 0)
	#
	pos :: Int = 0. # position in the genome
	state :: Int = 0.
	final_state :: Symbol = :missing # :lost, :fixed, :missing
	#
	# index where freqbin is reached
	timetofreq :: Dict{FrequencyBin, T}  = Dict{FT.FrequencyBin, eltype(t)}()
	ϕpos :: Vector{Union{Missing, Float64}} = ones(Missing, 0) # contribution of the position to fitness
	ϕgenotype :: Vector{Union{Missing, Float64}} = ones(Missing, 0) # average fitness of genotypes carying this mutation
end


text:
module TestWiringDiagramAlgorithms

using Test
using Catlab.Doctrines
using Catlab.WiringDiagrams

A, B, C, D = Ob(FreeBiproductCategory, :A, :B, :C, :D)
I = munit(FreeBiproductCategory.Ob)
f = Hom(:f, A, B)
g = Hom(:g, B, C)
h = Hom(:h, C, D)

# Diagonals and codiagonals
###########################

junction_diagram(args...) = singleton_diagram(Junction(args...))

# Add junctions for copies.
d = to_wiring_diagram(compose(f, mcopy(B)))
original = copy(d)
junctioned = compose(to_wiring_diagram(f), junction_diagram(:B,1,2))
@test add_junctions!(d) == junctioned
@test rem_junctions!(d) == original

d = to_wiring_diagram(compose(mcopy(A), otimes(f,f)))
original = copy(d)
junctioned = compose(junction_diagram(:A,1,2), to_wiring_diagram(otimes(f,f)))
@test add_junctions!(d) == junctioned
@test rem_junctions!(d) == original

# Add junctions for merges.
d = to_wiring_diagram(compose(mmerge(A), f))
original = copy(d)
junctioned = compose(junction_diagram(:A,2,1), to_wiring_diagram(f))
@test is_permuted_equal(add_junctions!(d), junctioned, [2,1])
@test rem_junctions!(d) == original

d = to_wiring_diagram(compose(otimes(f,f), mmerge(B)))
original = copy(d)
junctioned = compose(to_wiring_diagram(otimes(f,f)), junction_diagram(:B,2,1))
@test is_permuted_equal(add_junctions!(d), junctioned, [2,3,1])
@test rem_junctions!(d) == original

# Add junctions for deletions.
d = to_wiring_diagram(compose(f, delete(B)))
original = copy(d)
junctioned = compose(to_wiring_diagram(f), junction_diagram(:B,1,0))
@test add_junctions!(d) == junctioned
@test rem_junctions!(d) == original

# Add junctions for creations.
d = to_wiring_diagram(compose(create(A), f))
original = copy(d)
junctioned = compose(junction_diagram(:A,0,1), to_wiring_diagram(f))
@test is_permuted_equal(add_junctions!(d), junctioned, [2,1])
@test rem_junctions!(d) == original

# Add junctions for copies, merges, deletions, and creations, all at once.
d = to_wiring_diagram(compose(create(A),f,mcopy(B),mmerge(B),g,delete(C)))
original = copy(d)
junctioned = compose(
  junction_diagram(:A,0,1),
  to_wiring_diagram(f),
  junction_diagram(:B,1,2),
  junction_diagram(:B,2,1),
  to_wiring_diagram(g),
  junction_diagram(:C,1,0)
)
d = add_junctions!(d)
# XXX: An isomorphism test would be more convenient.
perm = [ findfirst([b] .== boxes(d)) for b in boxes(junctioned) ]
@test is_permuted_equal(d, junctioned, perm)
@test rem_junctions!(d) == original

# Normalize copies.
d = to_wiring_diagram(compose(mcopy(A), otimes(f,f)))
normalized = to_wiring_diagram(compose(f, mcopy(B)))
@test normalize_copy!(d) == normalized

d = to_wiring_diagram(compose(f, mcopy(B), otimes(g,g)))
normalize_copy!(d)
normalized = to_wiring_diagram(compose(f, g, mcopy(C)))
perm = sortperm(boxes(d); by=box->box.value)
@test is_permuted_equal(d, normalized, perm)

d = to_wiring_diagram(compose(mcopy(A), otimes(f,f), otimes(g,g)))
normalize_copy!(d)
perm = sortperm(boxes(d); by=box->box.value)
@test is_permuted_equal(d, normalized, perm)

# Normalize deletions.
d = to_wiring_diagram(f)
@test normalize_delete!(d) == to_wiring_diagram(f)

d = WiringDiagram(I,I)
add_box!(d, f)
@test normalize_delete!(d) == WiringDiagram(I,I)

d = WiringDiagram(A, B)
fv = add_box!(d, f)
gv = add_box!(d, g)
hv = add_box!(d, h)
add_wires!(d, [
  (input_id(d),1) => (fv,1),
  (fv,1) => (gv,1),
  (gv,1) => (hv,1),
  (fv,1) => (output_id(d),1),
])
@test normalize_delete!(d) == to_wiring_diagram(f)

# Normalize wiring diagrams representing morphisms in a cartesian category.
d = to_wiring_diagram(compose(
  mcopy(A),
  otimes(id(A),mcopy(A)),
  otimes(f,f,f),
  otimes(id(B), id(B), compose(g, delete(C)))
))
normalized = to_wiring_diagram(compose(f, mcopy(B)))
@test normalize_cartesian!(d) == normalized

# Layout
########

d = to_wiring_diagram(f)
@test crossing_minimization_by_sort(d, box_ids(d), sources=[input_id(d)]) == [1]

d = to_wiring_diagram(otimes(f,g))
@test crossing_minimization_by_sort(d, box_ids(d), sources=[input_id(d)]) == [1,2]

d = WiringDiagram(I,I)
fv1, fv2 = add_box!(d, f), add_box!(d, f)
gv1, gv2 = add_box!(d, g), add_box!(d, g)
add_wires!(d, [ Wire((fv1,1) => (gv1,1)), Wire((fv2,1),(gv2,1)) ])
@test crossing_minimization_by_sort(d, [gv1,gv2], sources=[fv1,fv2]) == [1,2]
@test crossing_minimization_by_sort(d, [fv1,fv2], targets=[gv1,gv2]) == [1,2]

end


text:
include("models/base_models.jl")
include("models/models.jl")


text:
# License for this file: MIT (expat)
# Copyright 2020-2021, DLR Institute of System Dynamics and Control
# Author: Martin Otter, DLR-SR
#
# Provide information about the structure of the ODE or DAE equations.

import DataFrames
import OrderedCollections
using  LinearAlgebra
using  StaticArrays
import TimerOutputs
import RecursiveFactorization

export StateCategory, XD, XALG, XLAMBDA, XMUE
export ResidualCategory, FD, FC_ALG, FC_LOW_HIGH, FC_MUE

export LinearEquations

export EquationInfo, StateElementInfo, update_equationInfo!, get_stateNames, get_x_table
export initialStateVector, updateEquationInfo!

export EquationInfoStatus, MANUAL, CODE_GENERATION, SOLVER_MODEL

"""
    isFixedLengthStartOrInit(startOrInit::Any, name::AbstractString)

Returns true, if `startOrInit` (so start value, init value, `nothing`) characterizes a value
with fixed length, so `length(v_startOrInit)` is fixed after compilation (= either a number or a StaticArray).
Note, if a start-value is not defined (`startOrInit = nothing`) a scalar default is used.

An error is triggered, if `startOrInit` is neither a `Number` nor an `AbstractVector` nor `nothing`.
In this case argument `name` is used in the error message (currently, EquationAndStateInfo is restricted
to numbers and vectors. In a subsequent version, two and more-dimensional arrays might be supported).

# Examples
```
isFixedLengthStartOrInit(1.0)                 # = true
isFixedLengthStartOrInit(SVector{3}(1,2,3))   # = true
isFixedLengthStartOrInit([1,2,3])             # = false (size can be changed after compilation)
isFixedLengthStartOrInit(nothing)             # = true (is assumed to be a scalar with value zero)
isFixedLengthStartOrInit(missing)             # error
```
"""
isFixedLengthStartOrInit(startOrInit, name::AbstractString) = begin
    if !(startOrInit isa Nothing || startOrInit isa Number || startOrInit isa AbstractVector)
        error("Start or init value $startOrInit of variable $name\nis neither a subtype of Number, subtype of AbstractVector, nor is nothing")
    end
    !(startOrInit isa AbstractVector && !(startOrInit isa StaticArray))
end

"""
    @enum StateCategory

Type of a state.

| value     | description                                                                       |
|:----------|:----------------------------------------------------------------------------------|
| `XD`      | Differential variable (derivative of variable appears in equations)               |
| `XALG`    | Algebraic variable that is included in x                                          |
| `XLAMBDA` | Algebraic variable that is included in der_x                                      |
| `XMUE`    | Algebraic variable used for stabilization and included in der_x (exact value = 0) |
"""
@enum StateCategory XD=1 XALG XLAMBDA XMUE

"""
    @enum ResidualCategory

Type of a residual (see also [`StateCategory`](@ref) for the type of a state).

| value         | description                                                              |
|:--------------|:-------------------------------------------------------------------------|
| `FD`          | Differential equation (depends on XD, XLAMBDA, XMUE variables)           |
| `FC_ALG`      | Constraint equation that is purely algebraic (depends on XALG variables) |
| `FC_LOW_HIGH` | Constraint equation on lowest and on highest derivative level            |
| `FC_MUE`      | Constraint equation that is associated with a stabilizer XMUE            |
"""
@enum ResidualCategory FD=1 FC_ALG FC_LOW_HIGH FC_MUE

const niter_max = 20  # Maximum number of fixed-point iterations to solve A*x = b

"""
    leq = LinearEquations{FloatType}(x_names::Vector{String}, x_vec_julia_names::AbstractVector,
                                     x_lengths::Vector{Int}, nx_fixedLength::Int, A_is_constant::Bool;
                                     useRecursiveFactorizationUptoSize = 0)

Define linear equation system "A*x=b" with `x::Vector{FloatType}`.

- `x` is constructed from a set of scalar or vector variables i with names `x_names[i]`
   and length `x_length[i]` (that is `length(x) = sum(x_lengths)`).
  x_names[1:nx_fixedLength] are elements with fixed lengths (dimensions do not change after compilation).
  x_names[nx_fixedLength+1:end] are vector-valued elements where the dimensions may changer after compilation.

- `x_vec_julia_names` are the Julia names of the vector-valued elements.

- If `A_is_constant = true` then `A` is a matrix that is constant after initialization

- If length(x) <= useRecursiveFactorizationUptoSize, then linear equation systems will be solved with
  `RecursiveFactorization.jl` instead of the default `lu!(..)` and `ldiv!(..)`.

For details how to use this constructor, see [`LinearEquationsIteration!`](@ref).
"""
mutable struct LinearEquations{FloatType <: Real}
    odeMode::Bool                     # Set from the calling function after LinearEquations was instantiated (default: true)
                                      # = true (standard mode): Compute "x" from equation "residuals = A*x - b"
                                      # = false (DAE solver): During events (including initialization)
                                      #   compute "x" as for odeMode=true. Outside of events:
                                      #   (1) "x" is set from the outside (= der(x_dae) provided by DAE solver)
                                      #   (2) Compute "residuals" from "residuals := A*x - b"
                                      #   (3) From the outside copy "residuals" into "residuals_dae" of the DAE solver.

    A_is_constant::Bool               # = true, if A-matrix is constant
    x_names::Vector{String}           # Names of the x-variables
    x_vec_julia_names::AbstractVector # Julia names of the vector-valued x-variables
                                      # (needed to generate code for
                                      #     if _m.storeResult
                                      #        x_vec_julia_names[1] = deepcopy(x_vec_julia_names[1])
                                      #        ...)
    x_lengths::Vector{Int}            # Lengths of the x-variables (sum(x_lengths) = length(x))
    nx_fixedLength::Int               # x_names[1:nx_fixedLength] are elements with fixed lengths (after compilation)
                                      # x_names[nx_fixedLength+1:end] are vector elements where the size may change after compilation
    x_vec::Vector{Vector{FloatType}}  # x_vec[1] is the value vector of x_names[nx_fixedLength+1]
                                      # x_vec[2] is the value vector of x_names[nx_fixedLength+2]
                                      # <..>

    A::Matrix{FloatType}
    b::Vector{FloatType}
    x::Vector{FloatType}              # Values of iteration variables
    pivots::Vector{Int}               # Pivot vector if recursiveFactorization = true
    residuals::Vector{FloatType}      # Values of the residuals FloatType vector; length(residuals) = length(x)

    # Iteration status of for-loop
    mode::Int       # Operational mode (see function LinearEquationsIteration!)
    niter::Int      # Current number of iterations in the fix point iteration scheme
    niter_max::Int  # Maximum number of iterations
    success::Bool   # = true, if solution of A*x = b is successfully computed
                    # = false, if solution is not computed; continue with fix point iteration

    # For warning message if niter > niter_max
    inconsistentPositive::Vector{String}
    inconsistentNegative::Vector{String}

    # Constructed during initialization
    useRecursiveFactorizationUptoSize::Int
    useRecursiveFactorization::Bool             # = true, if RecursiveFactorization.jl shall be used to solve the linear equation system

    luA::LU{FloatType,Array{FloatType,2}}       # lu-Decomposition of A

    function LinearEquations{FloatType}(x_names::Vector{String}, x_vec_julia_names::AbstractVector,
                                        x_lengths::Vector{Int}, nx_fixedLength::Int, A_is_constant::Bool;
                                        useRecursiveFactorizationUptoSize::Int = 0) where {FloatType <: Real}
        @assert(length(x_names) > 0)
        @assert(length(x_names) == length(x_lengths))
        nx = sum(x_lengths)
        @assert(nx > 0)
        useRecursiveFactorization = nx <= useRecursiveFactorizationUptoSize

        # Allocate value storage for vector elements
        nx_vec = length(x_names) - nx_fixedLength
        x_vec  = fill(FloatType[], nx_vec)
        j      = nx_fixedLength
        for i = 1:nx_vec
            j += 1
            x_vec[i] = zeros(FloatType, x_lengths[j])
        end

        new(true, A_is_constant, x_names, Any[], x_lengths, nx_fixedLength, x_vec,
            zeros(FloatType,nx,nx), zeros(FloatType,nx), zeros(FloatType,nx), fill(0,nx), zeros(FloatType,nx),
            -2, 0, niter_max, false, String[], String[],
            useRecursiveFactorizationUptoSize, useRecursiveFactorization)
    end
end
LinearEquations(args...) = LinearEquations{Float64}(args...)

"""
    b = copy_x_into_x_vec!(leq)

If vector valued elements of x-vector, copy vector valued elements of leq.x
into vector leq.x_vec. The function returns true
"""
function copy_x_into_x_vec!(leq::LinearEquations{FloatType})::Bool where {FloatType}
    i_vec = 0
    i_x   = leq.nx_fixedLength
    for i = leq.nx_fixedLength+1:length(leq.x_names)
        i_vec += 1
        x_vec = leq.x_vec[i_vec]
        for j = 1:leq.x_lengths[i]
            x_vec[j] = leq.x[i_x+j]
        end
        i_x += leq.x_lengths[i]
    end
    return true
end

"""
     iterating = LinearEquationsIteration!(leq, isInitial, [solve, isStoreResult,] time, timer)

This function solves a linear equation system in residual form "residual = A*x - b"
by iterating with a while loop over this system (arguments `solve, isStoreResult` are optional
and have only an effect if leq.odeMode=false, that is if the linear equation system is
solved in DAE mode):

```julia
function getDerivatives!(_der_x, _x, _m, _time)::Nothing
    _leq::Union{Nothing,LinearEquations{FloatType}} = nothing
    ...
    _leq      = _m.linearEquations[<nr>]   # leq::LinearEquations{FloatType}
    _leq.mode = -3  # initializes the iteration
    while LinearEquationsIteration!(_leq, _m.isInitial, _m.solve_leq, _m.storeResult, _m.time, _m.timer)
        x1 = _leq.x[1]
        x2 = _leq.x[2]
        x3 = _leq.x_vec[1]
        ...
        v_solved1 = f(x1, x2, ..., positive(x1,.., _leq))
        v_solved2 = f(x1, x2, ..., v_solved1)
        ...
        appendResidual!(_leq.residuals, < getResiduals(x1,x2,...) > )
    end
    ...
end
```

Note, A and b can be functions of linear event operators, such as `positive(c_i'*x - d_i)`.
In this case this system is solved by a fixed point iteration scheme,
that is `A*x = b` is solved until potentially present
`positive(c_i'*x - d_i)` calls are consistent to x
(`c_i'*x - d_i` and positive(..) must have the same sign).
Iteration is terminated after 20 iterations if no convergence is found.
In this case a warning message is triggered and simulation is continued
(it might be that simulation is still successful, even if `x` and
`positive(..)` are temporarily not consistent to each other).

The current values of A,x,b,residuals are stored in leq.
If A is fixed after initialization (leq.A_is_constant = true), then A is computed
only once at initialization, the LU decomposition of A is stored in leq
and used in subsequent calls to solve the equation system.

# Input arguments

- `leq::LinearEquations{FloatType}`: Instance of `LinearEquations`.
- `isInitial::Bool`: = true: Called during initialization.
- `solve::Bool`: = true: leq.x is computed by LinearEquationsIteration!.
                 = false: leq.x has been set by calling environment
                          (typically when called from DAE integrator).
                          Note, at events and at terminate, solve must be true).
- `isStoreResult::Bool`: = true: Called at a communication point (store results).
- `time`: Actual time instant (only used for warning/error messages).
- `timer::TimerOutputs.TimerOutput`: Timer used to measure the solution of the linear equation system.

# Output arguments

- iterating::Bool: = true : Continue iteration of while-loop.
                   = false: Immediately break while-loop after this function was terminated.

# Enhancing efficiency

Variable `leq.mode::Int` can be accessed (read-only) in the body of the while-loop to enhance efficiency.
For example, if residuals are computed in a function and the function evaluation is expensive.

```
leq.mode = -2: Residuals might be computed, but they are not used.
leq.mode = -1: Compute "residuals .= A*x - b"
leq.mode =  0: Compute "residuals .= A*0 - b"
leq.mode >  0: Compute "residuals .= A*e_i - b"   # e_i = unit vector with i = leq.mode
```

# Hidden argument `leq.mode::Int` on input

```julia
leq.mode = -3  # LinearEquationsIteration! is called the first time
               if leq.odeMode || solve
                  # ODE mode or DAE mode at an event (solve "x" from equation "residuals = A*x - b")
                  # Initialize fixed point iteration or continue fixed point iteration (if in DAE mode)
                  leq.niter     = 0      # Number of event iterations
                  leq.success   = false  # Event iteration was not yet successful
                  leq.x        .= 0
                  leq.mode      = 0      # Compute "residuals .= A*0 - b"
               else # solve = false
                  # DAE mode (but not at an event)
                  leq.niter     = 0
                  leq.success   = true
                  leq.x_is_zero = false
                  if isStoreResult
                     leq.mode = -2  # Residuals might be computed, but they are not used.
                  else
                     leq.mode = -1  # Compute "residuals .= A*x - b"
                  end
               end
               return true  # Continue while-loop

leq.mode = -2  # Terminate while-loop or initialize next event iteration
               if leq.success           # Set from leq.mode=-2 or leq.mode>=0 or from event operators, such as positive!(...)
                  return false          # Terminate while-loop
               elseif leq.niter > 20
                  <warning>             # Event iteration failed (variables are not consistent)
                  return false          # Terminate while-loop
               else
                  # Initialize next event iteration
                  leq.x        .= 0
                  leq.x_is_zero = true
                  leq.mode      = 0     # Compute "residuals .= A*0 - b"
                  return true           # Continue while-loop
               end

leq.mode = -1: @assert(!leq.odeMode && !solve)
               # DAE mode, but not at an event
               return false      # Terminate while-loop

leq.mode =  0: @assert(leq.odeMode || solve)
               # ODE mode or DAE mode at an event (solve "x" from equation "residuals = A*x - b")
               leq.b    = -leq.residuals
               leq.x[1] = 1.0
               leq.mode = 1      # Compute "residuals := A*e_1 - b"
               return true       # Continue while-loop

leq.mode >  0: @assert(leq.odeMode || solve)
               j          = leq.mode
               leq.A[:,j] = leq.residuals + leq.b
               leq.x[j]   = 0.0
               if j < nx
                   leq.x[j+1] = 1.0
                   leq.mode   = j+1      # Compute A[:,j+1] after the next iteration
                   return true           # Continue while-loop
               elseif j == nx
                   leq.x = solve(leq.A,leq.b) # Solve linear equation system A*x = b for x.
                   leq.mode    = -1      # Compute all variables IN the next iteration as function of x
                   leq.niter  += 1       # Increment number of iterations to solve A*x = b
                   leq.success = true    # Terminate for-loop at the beginning of the next iteration,
                                         # provided no positive(..) call changes its value.
                                         # (leq.success is set to false in positive(..), if the return value changes).
                   return true           # Continue while-loop
               end
```
"""
LinearEquationsIteration!(leq::LinearEquations, isInitial, time, timer) = LinearEquationsIteration!(leq, isInitial, true, false, time, timer)
function LinearEquationsIteration!(leq::LinearEquations{FloatType}, isInitial::Bool, solve::Bool,
                                   isStoreResult::Bool, time, timer)::Bool where {FloatType}
    mode = leq.mode
    nx   = length(leq.x)

    if mode == -3
        # LinearEquationsIteration! is called the first time in the current model evaluation
        leq.niter = 0      # Number of event iterations
        empty!(leq.inconsistentPositive)
        empty!(leq.inconsistentNegative)

        if leq.odeMode || solve
           # ODE mode or DAE mode at an event (solve "x" from equation "residuals = A*x - b")
           # Initialize fixed point iteration or continue fixed point iteration (if in DAE mode)
           leq.success  = false  # Event iteration was not yet successful
           leq.x       .= 0
           leq.mode     = 0      # Initialize leq and compute "residuals .= A*0 - b"
        else
           # DAE mode (but not at an event)
           leq.success   = true
           if isStoreResult
              leq.mode = -2   # Residuals might be computed, but they are not used.
           else
              leq.mode = -1   # Initialize leq and compute "residuals .= A*x - b"
           end
        end
        empty!(leq.residuals)
        return copy_x_into_x_vec!(leq)  # Continue while-loop

    elseif mode == -2
        # Terminate while-loop or initialize next event iteration
        if leq.success   # Set from leq.mode=-2 or leq.mode>=0 or from event operators, such as positive!(...)
            return false # Terminate while-loop
        elseif leq.niter > niter_max
            str = ""
            if length(leq.inconsistentPositive) > 0
                str = str * "positive(expr) is inconsistent for expr = $(leq.inconsistentPositive)."
            end
            if length(leq.inconsistentNegative) > 0
                if length(str) > 0
                    str = str * "\n"
                end
                str = str * "negative(expr) is inconsistent for expr = $(leq.inconsistentNegative)."
            end
            @warn "At time = $time, no consistent solution found for mixed linear equation system.\n" *
                  "Simulation is continued although some variables might not be correct at this time instant.\n$str"
            return false       # Terminate while-loop
        end
        leq.x        .= 0
        leq.mode      = 0      # Compute "residuals .= A*0 - b"
        empty!(leq.residuals)
        return copy_x_into_x_vec!(leq) # Continue while-loop

    elseif mode < -3 || mode > nx
        @goto ERROR
    end

    x = leq.x
    A = leq.A
    b = leq.b
    residuals = leq.residuals

    if length(residuals) != length(x)
        error("Function LinearEquationsIteration! wrongly used:\n",
              "length(leq.residuals) = ", length(leq.residuals), ", length(leq.x) = ", length(leq.x))
    end

    if mode == -1
        @assert(!leq.odeMode && !solve)
        return false   # Terminate while-loop (leq.residuals must be copied into DAE residuals)

    elseif !leq.A_is_constant || isInitial  # A is not constant or A is constant and isInitial = true
        if mode == 0
            # residuals = A*x - b -> b = -residuals)
            for i = 1:nx
                b[i] = -residuals[i]
            end
            leq.mode = 1
            x[1] = max(FloatType(1e-3), abs(b[1]))     #convert(FloatType, 1)
            empty!(leq.residuals)
            return copy_x_into_x_vec!(leq)
        end

        # residuals = A*e_j*x_j - b -> A[:,j] = (residuals + b)/x[j]
        j = mode
        for i = 1:nx
            A[i,j] = (residuals[i] + b[i])/x[j]
        end
        x[j] = 0

        if j < nx
            leq.mode += 1
            x[leq.mode] = max(FloatType(1e-3), abs(b[leq.mode]))    # convert(FloatType, 1)
            empty!(leq.residuals)
            return copy_x_into_x_vec!(leq)
        end

        # Solve linear equation system
        if nx == 1
            x[1] = b[1]/A[1,1]
            if !isfinite(x[1])
                error("Linear scalar equation system is singular resulting in: ", leq.x_names[1], " = ", x[1])
            end
        else
            x .= b
            if leq.useRecursiveFactorization
                Modia.TimerOutputs.@timeit timer "Modia LinearEquationsIteration! (solve A*x=b Rec.Fac.)" begin
                    leq.luA = RecursiveFactorization.lu!(A, leq.pivots)
                    ldiv!(leq.luA, x)
                end
            else
                Modia.TimerOutputs.@timeit timer "Modia LinearEquationsIteration! (solve A*x=b)" begin
                    leq.luA = lu!(A)
                    ldiv!(leq.luA, x)
                end
            end
        end

    elseif leq.A_is_constant && !isInitial  # isInitial=false, LU decomposition of A is available in leq.luA
        for i = 1:nx
            x[i] = -residuals[i]
        end

        # Solve linear equation system
        if nx == 1
            x[1] = x[1]/A[1,1]
            if !isfinite(x[1])
                error("Linear scalar equation system is singular resulting in: ", leq.x_names[1], " = ", x[1])
            end
        else
            ldiv!(leq.luA, x)
        end

    else
        @goto ERROR
    end

    # Linear equation system solved
    leq.mode    = -2   # Compute all variables IN the next iteration as function of x
    leq.niter  += 1    # Increment number of iterations to solve A*x = b
    leq.success = true # Terminate for-loop at the beginning of the next iteration,
                       # provided no positive(..) call changes its value.
                       # (leq.success is set to false in positive(..), if the return value changes).

    empty!(leq.residuals)
    return copy_x_into_x_vec!(leq)

    @label ERROR
    error("Should not occur (Bug in file Modia/src/EquationAndStateInfo.jl):\n",
          "   time           = $time\n",
          "   isInitial      = $isInitial,\n",
          "   solve          = $solve,\n",
          "   isStoreResult  = $isStoreREsult,\n",
          "   leq.odeMode    = $(leq.odeMode),\n",
          "   leq.mode       = $(leq.mode),\n",
          "   leq.x_names    = $(leq.x_names),\n",
          "   leq.x_lengths  = $(leq.x_lengths),\n")
end

"""
    @enum EquationInfoStatus

Status of an EquationInfo instance:

- `MANUAL`: Is defined manually. The following variables in x_info::Vector{StateElementInfo}
   are **not** defined:
   x_names_julia, der_x_names_julia, length, unit, startIndex.
   Also variables nx and x_infoByIndex are not defined.
   With function [`update_equationInfo!`](@ref), the variables length, unit, startIndex
   in x_info::Vector{StateElementInfo} are computed, as well as nx and x_infoByIndex.

- `CODE_GENERATION`: Is constructed during code generation with getSortedAndSolvedAST(..).
  The following variables in x_info::Vector{StateElementInfo}
  are **not** defined: startIndex.
  Also variables nx and x_infoByIndex are not defined.
  With function [`update_equationInfo!`](@ref), the variables startIndex
  in x_info::Vector{StateElementInfo} are computed, as well as nx and x_infoByIndex.

- `SOLVER_MODEL`: Is used during simulation in a SolverModel. With function
  [`update_equationInfo!`](@ref), missing variables are constructed depending
  on the information given with `MANUAL` or `CODE_GENERATION` and the actual
  modelValues instance (here unit and length information is available).
"""
@enum EquationInfoStatus MANUAL=1 CODE_GENERATION SOLVER_MODEL

"""
    xe_info = StateElementInfo(...)

Return an instance of the mutable struct `StateElementInfo` that defines the information
for one element of the state vector. There are three constructors:

- Default constructor (all variables under section Arguments are given;
  used to write/read the complete information).

- All variables under section Arguments are given with exception of startIndex
  (used for `EquationInfoStatus = CODE_GENERATION`).

- StateElementInfo(x_name, der_x_name, stateCategory; fixed=nothing, nominal=NaN, unbounded=false)
  (used for `EquationInfoStatus = MANUAL`):

# Arguments
- x_name: Name of x-element or "" if no name (if stateCatebory = XLAMBDA or XMUE)
- x_name_julia: Julia name of x-element in getDerivatives!/getResiduals! function
  or "" if not needed (since no code generation).
- der_x_name: Name of der_x-element or "" if either `der(x_name)` or if no name
  (if stateCategory = XALG).
- der_x_name_julia: Julia name of der_x-element in getDerivatives!/getResiduals! function
  or "" if not needed (since no code generation).
- stateCategory::StateCategory: Category of the state
- length: length of x-element (or -1 if not yet known)
- unit: unit of XD, XALG (x_name) or XLAMBDA, XMUE (der_x_name) or "" if not yet known)
- fixed: false (= guess value) or true (= not changed by initialization).
         Only relevant for ode=false, otherwise ignored.
- nominal: Nominal value (NaN if determined via start value)
- unbounded: false or true
- startIndex: start index of x-element with respect to x-vector
              or -1 if not yet known.
"""
mutable struct StateElementInfo
    x_name::String                # Modia name of x-element or "" if no name (for )
    x_name_julia                  # Julia name of x-element in getDerivatives! function
                                  # or "", if not needed (since no code generation).
    der_x_name::String            # Modia name of der_x-element or "" if either "der(x_name)" or if no name,
    der_x_name_julia              # Julia name of der_x-element in getDerivatives! function
                                  # or not needed (since no code generation)
    stateCategory::StateCategory  # category of the state
    unit::String                  # unit of x-element as string (or "" if not yet known)
    startOrInit::Any              # start or init value or nothing depending on fixed
                                  # (init   : if value and fixed=true
                                  #  start  : if value and fixed=false
                                  #  nothing: neither start nor init value; fixed=false)
    fixed::Bool                   #
    nominal::Float64              # nominal value (NaN if determined via start value)
    unbounded::Bool               # false or true
    fixedLength::Bool             # = true, if length of value is fixed at compile time (= scalar or StaticArray)
                                  # = false, if length of value is determined before initialization
    scalar::Bool                  # = true, if scalar
                                  # = false, if vector
    length::Int                   # length of x-element
    startIndex::Int               # start index of element with respect to x-vector (if fixedLength=true)
                                  # or with respect to x_vec-vector (if fixedLength=false)
end

# Constructor for code-generation
StateElementInfo(x_name, x_name_julia, der_x_name, der_x_name_julia,
                 stateCategory, unit, startOrInit, fixed, nominal, unbounded) = StateElementInfo(
                 x_name, x_name_julia, der_x_name, der_x_name_julia,
                 stateCategory, unit, startOrInit, fixed, nominal, unbounded,
                 isFixedLengthStartOrInit(startOrInit, x_name), !(startOrInit isa AbstractArray),
                 startOrInit isa Nothing ? 1 : length(startOrInit), -1)

# Constructor for reading StateElementInfo after code-generation (x_name_julia and der_x_name_julia are not included
#StateElementInfo(x_name, der_x_name,
#                 stateCategory, unit, startOrInit, fixed, nominal, unbounded) = StateElementInfo(
#                 x_name, :(), der_x_name, :(),
#                 stateCategory, unit, startOrInit, fixed, nominal, unbounded, -1)

# Constructor for manually generated equation info.
#StateElementInfo(x_name, der_x_name, stateCategory=XD; fixed=false, nominal=NaN, unbounded=false) =
#    StateElementInfo(x_name, :(), der_x_name, :(), stateCategory, "", nothing, fixed, nominal, unbounded, -1)

function Base.show(io::IO, xe_info::StateElementInfo)
    print(io, "Modia.StateElementInfo(")
    show( io, xe_info.x_name)
	#print(io, ",")
    #show( io, xe_info.x_name_julia)
	print(io, ",")
    show( io, xe_info.der_x_name)
	#print(io, ",")
    #show( io, xe_info.der_x_name_julia)
    print(io, ",Modia.", xe_info.stateCategory)
	print(io, ",")
    show( io, xe_info.unit)
    print(io, ",", xe_info.startOrInit)
    print(io, ",", xe_info.fixed)
    print(io, ",", xe_info.nominal)
    print(io, ",", xe_info.unbounded)
    print(io, ",", xe_info.fixedLength)
    print(io, ",", xe_info.scalar)
    print(io, ",", xe_info.length)
    print(io, ",", xe_info.startIndex)
    print(io, ")")
    return nothing
end

"""
    x_table = get_x_table(x_info::Vector{StateElementInfo})

Return the state element info as DataFrames.DataFrame table, for example
to print it as:

```
show(x_table, allrows=true, allcols=true, summary=false, eltypes=false)
```
"""
function get_x_table(x_info::Vector{StateElementInfo})
    x_table = DataFrames.DataFrame(name=String[], unit=String[], startOrInit=[], fixed=Bool[], nominal=Float64[], unbounded=Bool[])

    for xe_info in x_info
        push!(x_table, (xe_info.x_name, xe_info.unit, xe_info.startOrInit, xe_info.fixed, xe_info.nominal, xe_info.unbounded))
    end

    return x_table
end

"""
    nx = stateVectorLength(x_info::Vector{StateElementInfo})

Return the length of the state vector
"""
function stateVectorLength(x_info::Vector{StateElementInfo})::Int
    nx = 0
    for xe_info in x_info
        nx += xe_info.length
    end
    return nx
end

"""
    eqInfo = EquationInfo(;
                status               = MANUAL,
                ode                  = true,
                nz                   = 0,
                x_info               = StateElementInfo[],
                nx                   = -1,
                nxFixedLength        = nx,
                residualCategories   = ResidualCategory[],
                linearEquations      = Tuple{Vector{String},AbstractVector,Vector{Int},Int,Bool}[],
                vSolvedWithFixedTrue = String[],
                defaultParameterAndStartValues = nothing,
                ResultType = nothing,
                ResultTypeHasFloatType = false)

Return instance `eqInfo` that defines the information for the equation system.

# Arguments
- status::EquationInfoStatus: Defines the variables that have a value.
- ode: = true if ODE-interface (`getDerivatives!`),
       = false if DAE-interface (`getResiduals!`).
- nz: Number of zero crossing functions.
- x_info: Vector of StateElementInfo elements provding info for every x-element
- residualCategories: If ode=true, length(residualCategories) = 0.
             If ode=false: residualCategories[i] is the `ResidualCategory`](@ref) of x-element "i".
- `linearEquations::Vector{Tuple{Vector{String},AbstractVector,Vector{Int},Int,Bool}}`:
               linearEquations[i] defines a
               Modia.LinearEquations system, where the first tuple value
               is a vector of the names of the unknowns,
               the second tuple value is a vector of the Julia names of the vector-valued elements,
               the third tuple value is a vector with the lengths of the unknowns,
               the fourth tuple value is the number of residuals and the fifth tuple value
               defines whether the coefficient matrix A
               has only constant entries (=true) or not (=false).
- `vSolvedWithFixedTrue::Vector{String}`: Vector of variables that are computed
               from other variables and have `fixed=true`. During initialization
               it is checked whether the calculated values and the start values of
               these variables agree. If this is not the case, an error is triggered.
- `defaultParameterAndStartValues`: Dictionary of default paramter and default start values.
- `ResultType::AbstractModelValues`: ModelValues type that shall be used for the result generation
               (= ModelValues struct without parameters).
- `ResultTypeHasFloatType::Bool`: = false, if `ResultType` is not parameterized.
                                  = true, if `ResultType` has parameter `FloatType`.
"""
mutable struct EquationInfo
    status::EquationInfoStatus
    ode::Bool                                      # = true if ODE model interface, otherwise DAE model interface
    nz::Int                                        # Number of crossing functions
    x_info::Vector{StateElementInfo}
    residualCategories::Vector{ResidualCategory}   # If ode=true, length(residualCategories) = 0
                                                   # If ode=false, residualCategories[j] is the ResidualCategory of residual[j]
    linearEquations::Vector{Tuple{Vector{String},AbstractVector,Vector{Int},Int,Bool}}
    vSolvedWithFixedTrue::Vector{String}
    nx::Int                                        # = length(x) or -1 if not yet known
    nxFixedLength::Int                             # x_info[1:nxFixedLengt] are states with fixed length (does not change after compilation)
    x_infoByIndex::Vector{Int}                     # i = x_infoByIndex[j] -> x_info[i]
                                                   # or empty vector, if not yet known.
    x_dict::OrderedCollections.OrderedDict{String,Int}      # x_dict[x_name] returns the index of x_name with respect to x_info
    der_x_dict::OrderedCollections.OrderedDict{String,Int}  # der_x_dict[der_x_name]  returns the index of der_x_name with respect to x_info
    defaultParameterAndStartValues::Union{AbstractDict,Nothing}
    ResultType
    ResultTypeHasFloatType::Bool
end

EquationInfo(; status                = MANUAL,
               ode                   = true,
               nz                    = 0,
               x_info                = StateElementInfo[],
               residualCategories    = ResidualCategory[],
               linearEquations       = Tuple{Vector{String},AbstractVector,Vector{Int},Int,Bool}[],
               vSolvedWithFixedTrue  = String[],
               nxFixedLength         = -1,
               x_infoByIndex         = Int[],
               defaultParameterAndStartValues::Union{AbstractDict,Nothing} = nothing,
               ResultType = nothing,
               ResultTypeHasFloatType = false) = EquationInfo(status, ode, nz, x_info,
                                                    residualCategories, linearEquations,
                                                    vSolvedWithFixedTrue, stateVectorLength(x_info), -1, x_infoByIndex,
                                                    OrderedCollections.OrderedDict{String,Int}(),
                                                    OrderedCollections.OrderedDict{String,Int}(),
                                                    defaultParameterAndStartValues,
                                                    ResultType, ResultTypeHasFloatType)

"""
    initEquationInfo!(eqInfo::EquationInfo)

Set eqInfo.x_dict, eqInfo.der_x_dict, eqInfo.nx and eqInfo.x_info[:].startIndex
"""
function initEquationInfo!(eqInfo::EquationInfo)::Nothing
    x_dict     = eqInfo.x_dict
    der_x_dict = eqInfo.der_x_dict
    startIndex = 1
    for (i, xi_info) in enumerate(eqInfo.x_info)
        x_dict[xi_info.x_name] = i
        der_x_dict[xi_info.der_x_name] = i
        xi_info.startIndex = startIndex
        startIndex += xi_info.length
    end
    eqInfo.nx = startIndex - 1
    return nothing
end

"""
    x_init = initialStateVector(eqInfo::EquationInfo, FloatType)

Return initial state vector `Vector{FloatType}` with stripped units.
"""
function initialStateVector(eqInfo::EquationInfo, FloatType::Type)::Vector{FloatType}
    x_start = fill(FloatType(0), eqInfo.nx)
    startIndex = 1
    for xe_info in eqInfo.x_info
        if xe_info.scalar
            if !(xe_info.startOrInit isa Nothing)
                x_start[startIndex] = FloatType(ustrip(xe_info.startOrInit))
                startIndex += 1
            end
        else
            for i = 1:xe_info.length
                x_start[startIndex] = FloatType(ustrip(xe_info.startOrInit[i]))
                startIndex += 1
            end
        end
    end
    return x_start
end

"""
    x_start = updateEquationInfo!(eqInfo::EquationInfo, FloatType)

Set eqInfo.x_dict, eqInfo.der_x_dict, eqInfo.nx and eqInfo.x_info[:].startIndex
"""
function updateEquationInfo!(eqInfo::EquationInfo, FloatType::Type)::Vector{FloatType}
    nxFixedLength = eqInfo.nxFixedLength
    if nxFixedLength == 0
        startIndex = 1
    else
        xi_info = eqInfo.x_info[nxFixedLength]
        startIndex = xi_info.startIndex + xi_info.length
    end

    x_info = eqInfo.x_info
    for i = nxFixedLength+1:length(x_info)
        xi_info = x_info[i]
        xi_info.length     = length(xi_info.startOrInit)
        xi_info.startIndex = startIndex
        startIndex        += xi_info.length
    end
    eqInfo.nx = startIndex - 1

    return initialStateVector(eqInfo, FloatType)
end

get_x_startIndexAndLength(eqInfo::EquationInfo, name::String) = begin
    xe_info = eqInfo.x_info[ eqInfo.x_dict[name] ]
    (xe_info.startIndex, xe_info.length)
end

"""
    nvec = get_equationSizes(equationInfo)

Return the sizes of the linear equations as `nvec::Vector{Int}`.
(e.g. `nvec=[10,3]` means that there are two linear equation systems with
size 10 and size 3).
"""
get_equationSizes(eqInfo::EquationInfo) = Int[leq[3] for leq in eqInfo.linearEquations]

function Base.show(io::IO, eqInfo::EquationInfo; indent=4)
    indentation  = repeat(" ", indent)
    indentation2 = repeat(" ", 2*indent)
    indentation3 = repeat(" ", 3*indent)
    println(io, "Modia.EquationInfo(")
    println(io, indentation2, "status = Modia.", eqInfo.status, ",")
    println(io, indentation2, "ode = ", eqInfo.ode, ",")
    if eqInfo.nz > 0
        println(io, indentation2, "nz = ", eqInfo.nz, ",")
    end
    print(io, indentation2, "x_info = Modia.StateElementInfo[")
    for (i, xe_info) in enumerate(eqInfo.x_info)
        if i == 1
            print(io, "\n", indentation3)
        else
            print(io, ",\n", indentation3)
        end
        show(io, xe_info)
    end
    print(io,"]")

    if length(eqInfo.residualCategories) > 0
        print(io, ",\n", indentation2, "residualCategories = [")
        for (i, rcat) in enumerate(eqInfo.residualCategories)
            if i > 1
                print(io, ",")
            end
            show(io, rcat)
        end
        print(io,"]")
    end

    leqs = eqInfo.linearEquations
    if length(leqs) > 0
        println(io, ",\n", indentation2, "linearEquations = [")
        for (i, leq) in enumerate(leqs)
            print(io, indentation3, "(", leq[1], ",\n",
                      indentation3, " ", leq[2], ",\n",
                      indentation3, " ", leq[3], ", ", leq[4], ")")
            if i < length(leqs)
                print(io, ",\n")
            end
        end
        print(io, "]")
    end

    if length(eqInfo.vSolvedWithFixedTrue) > 0
        print(io, ",\n", indentation2, "vSolvedWithFixedTrue = ")
        show(io, eqInfo.vSolvedWithFixedTrue)
    end

    if eqInfo.nx > 0
        print(io, ",\n", indentation2, "nx = ", eqInfo.nx)
    end

    if length(eqInfo.x_infoByIndex) > 0
        print(io, ",\n", indentation2, "x_infoByIndex = ")
        show(io, eqInfo.x_infoByIndex)
    end

    if !isnothing(eqInfo.defaultParameterAndStartValues)
        print(io, ",\n", indentation2, "defaultParameterAndStartValues = ")
        show(io, eqInfo.defaultParameterAndStartValues, indent=12, finalLineBreak=false)
    end

    if !isnothing(eqInfo.ResultType)
        print(io, ",\n", indentation2, "ResultType = ")
        show(io, eqInfo.ResultType)

        if eqInfo.ResultTypeHasFloatType
            print(io, ",\n", indentation2, "ResultTypeHasFloatType = ")
            show(io, eqInfo.ResultTypeHasFloatType)
        end
    end

    println(io, "\n", indentation, ")")
    return nothing
end

"""
    names = get_stateNames(equationInfo::EquationInfo)

Return the names of the states defined in `equationInfo` as a Vector of strings.
"""
get_stateNames(eqInfo::EquationInfo) = String[xi_info.x_name for xi_info in eqInfo.x_info]

"""
    names = get_xNames(eqInfo::EquationInfo)

Return the names of all elements of the x-vector as a vector of strings.
"""
function get_xNames(eqInfo::EquationInfo)::Vector{String}
    xNames = Vector{String}(undef, eqInfo.nx)
    for xe_info in eqInfo.x_info
        if xe_info.length == 1
            xNames[xe_info.startIndex] = xe_info.x_name
        else
            for i = 1:xe_info.length
                xNames[xe_info.startIndex+i-1] = xe_info.x_name*"["*string(i)*"]"
            end
        end
    end
    return xNames
end

#=
"""
    update_equationInfo!(eqInfo::EquationInfo)

Update equationInfo as needed for [`EquationInfoStatus`](@ref)` = SOLVER_MODEL`.
"""
function update_equationInfo!(eqInfo::EquationInfo, modelValues::Modia.AbstractModelValues)::Nothing
    x_info = eqInfo.x_info

    if eqInfo.status == MANUAL
        # Determine length and unit
		count = 0
        for xi_info in x_info
            if xi_info.x_name != "" # XD or XALG
                (component,name) = Modia.get_modelValuesAndName(modelValues, xi_info.x_name)
            elseif xi_info.der_x_name != ""  # XLAMBDA or XMUE
                (component,name) = Modia.get_modelValuesAndName(modelValues, xi_info.der_x_name)
            else
				# Pure algebraic with dummy equation der(x) = -x, x(0) = 0
				xi_info.length = 1
				count = count+1
				continue
            end

	        xi_type        = fieldtype(typeof(component), name)
			xi_numberType  = Modia.numberType(xi_type)
            xi_info.unit   = replace(string(unit(xi_numberType)), " " => "*")
            xi_info.length = Base.length(xi_type)  # Base.length(::Type{<:Number})=1 is defined in Modia, since not defined in Base
        end
		if count == 1 && length(x_info) > 1 || count >  1
			error("Error in update_equationInfo!: x_info is wrong.\n",
				  "x_info = $x_info")
		end
    end

    # Set startIndex in x_info and compute nx
    startIndex = 1
    for xi_info in x_info
        xi_info.startIndex = startIndex
        startIndex += xi_info.length
    end
    nx = startIndex - 1
    eqInfo.nx = nx

    # Set x_infoByIndex::Vector{Int}
    eqInfo.x_infoByIndex = zeros(Int,nx)
    i1 = 0
    i2 = 0
    for (i,xi_info) in enumerate(x_info)
        i1 = xi_info.startIndex
        i2 = i1 + xi_info.length - 1
        for j = i1:i2
            eqInfo.x_infoByIndex[j] = i
        end
    end

    eqInfo.status = SOLVER_MODEL

    return nothing
end
=#


text:
using Visualize, GeometryTypes, ModernGL
using Visualize: GLRasterization, TextUniforms, get_texture_atlas, Sprite2
using Visualize.GLRasterization: show!, destroy!
using GPUArrays

resolution = (1024, 1024)
window = GLFWWindow(Area => resolution)
GLBackend.init(window[Window])

for event in Visualize.NativeWindowEvents
    add!(window, event)
end
show!(window)
window[Visualize.Open] = true

proj = orthographicprojection(SimpleRectangle(0, 0, resolution...), -10_000f0, 10_000f0)
uniforms = TextUniforms(
    proj,
    Vec4f0(1, 0, 0, 1),
    Vec4f0(1, 0, 1, 1)
)

atlas = get_texture_atlas();

text = Visualize.Text(
    Sprite2[],
    Char[],
    atlas,
    [0],
    Visualize.defaultfont(),
    45f0,
    0f0,
    Vec4f0(0f0, 0f0, 0f0, 1f0),
    Point2f0(0),
    1.5
);

print(text, 'c')
print(text, " nice test bro!")
text.color = Vec4f0(1f0, 0f0, 0f0, 1f0)
print(text, " Muashahahaha@")

uniform_buff = UniformBuffer(uniforms);
vbo = VertexArray(text.data, face_type = Face{1, OffsetInteger{1, GLint}});
Visualize.vertex_main(text.data[1], uniforms, ())
atlas_tex = GPUArrays.buffer(GLRasterization.atlas_texture(atlas));
draw_particles = GLRasterizer(
    vbo, (uniform_buff, atlas_tex),
    Visualize.vertex_main, Visualize.fragment_main;
    geometryshader = Visualize.geometry_main
)

glDisable(GL_DEPTH_TEST)
glClearColor(1, 1, 1, 0)
GLAbstraction.enabletransparency()

@async begin
    while isopen(window)
        GLFW.PollEvents()
        glViewport(0, 0, widths(window[Area])...)
        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
        draw_particles(vbo, (uniform_buff, atlas_tex))
        GLRasterization.swapbuffers!(window)
        sleep(0.001)
        yield()
    end
    destroy!(window)
end