text:
include("./artificial_data.jl")

using CellwiseEditingDifferentiation
using Gadfly

hyperpars = Dict{Symbol, Any}(:λ=>fill(10.0, nsites))

prior = Dict{Symbol, Any}(:v=>Function[(m::Float64, v::Float64, a::Float64, b::Float64)->
                          vpcprior(m, v, a, b, hyperpars[:λ][i]) for i in 1:nsites])

target = Dict{Symbol, Any}(:v=>Array(Function, nsites))
for i in 1:nsites
  target[:v][i] = function (v::Float64)
    a, b = beta_pars_from_mv(data[:m][i], v)
    exp(prior[:v][i](data[:m][i], v, a, b))
  end
end

support = Dict{Symbol, Any}(:v=>Any[0.0001:0.0001:0.9999*data[:m][i]*(1-data[:m][i]) for i in 1:nsites])

# vprior_xmin = fill(0., nsites)
# vprior_xmax = fill(0.03, nsites)
vprior_ymin = fill(0., nsites)
vprior_ymax = fill(40., nsites)

for i in 1:nsites
  vpdf = Float64[target[:v][i](x) for x in support[:v][i]]
  c = quadgk(target[:v][i], support[:v][i][1], support[:v][i][end])

  colors = Dict{Symbol, Color.RGB{Float64}}(:posterior=>color("red"), :prior=>color("blue"))

  layers = Layer[]

  push!(layers, layer(
    x=vec(readdlm(joinpath(OUTDIR, @sprintf("vchain_%s_site%02d.txt", string(simulationid), i)))),
    Stat.histogram(bincount=50, density=true),
    Geom.line,
    Theme(default_color=colors[:posterior])
  )[1])

  push!(layers, layer(
    x=collect(support[:v][i]),
    y=vpdf/c[1],
    Geom.line,
    Theme(default_color=colors[:prior])
  )[1])

  vplot = plot(
    layers,
    Guide.xlabel("v$i"),
    Guide.title("Histogram of v$i")
    # Guide.manual_color_key("Distribution", [string(k) for k in keys(colors)], [c for c in values(colors)]),
    # Coord.Cartesian(xmin=vprior_xmin[i], xmax=vprior_xmax[i], ymin=vprior_ymin[i], ymax=vprior_ymax[i])
  )

  draw(PDF(joinpath(OUTDIR, @sprintf("vhist_%s_site%02d.pdf", string(simulationid), i)), 4inch, 3inch), vplot)
end


text:
function conv2!(p,ip,At::StridedMatrix{T},Bt::StridedMatrix{T},C::StridedMatrix{T},
                A::StridedMatrix{T},B::StridedMatrix{T}) where T <: Complex
    sa,sb = size(A),size(B)

    At[:,:] = zero(T)
    Bt[:,:] = zero(T)
    At[1:sa[1],1:sa[2]] .= A
    Bt[1:sb[1],1:sb[2]] .= B

    A_mul_B!(C,ip,(p*At).*(p*Bt))
    return C
end
conv2!(A::StridedMatrix{T},B::StridedMatrix{T}) where {T<:Integer} =
       round.(Int,conv2(float(A), float(B)))
conv2!(u::StridedVector{T},v::StridedVector{T},A::StridedMatrix{T}) where {T<:Integer} =
       round.(Int, conv2(float(u),float(v),float(A)))

function density!(dp,p,tt,params)
    plan,iplan,At,Bt,C,vv,f = getindex.(params,[:plan,:iplan,:At,:Bt,:C,:vv,:f])
    nn = size(p)[1]
    kk = repmat(-round(Int, nn/2):round(Int, nn/2-1),1,nn)
    ll = copy(kk)'
    inds = (-round(Int,nn/2):round(Int,nn/2-1))+nn+1
    dp[:,:] .= (-im*2π*(f[1]*kk.+f[2]*ll).*p
                 .- 4π^2*(  conv2!(plan,iplan,At,Bt,C,p.*(kk.^2.+ll.^2),tt)[inds,inds]
                          .+conv2!(plan,iplan,At,Bt,C,p.*kk,tt.*kk)[inds,inds]
                          .+conv2!(plan,iplan,At,Bt,C,p.*ll,tt.*ll)[inds,inds])
                 .-4π^2*(  conv2!(plan,iplan,At,Bt,C,p,vv.*(kk.^2.+ll.^2))[inds,inds]
                         .+conv2!(plan,iplan,At,Bt,C,p.*kk,vv.*kk)[inds,inds]
                         .+conv2!(plan,iplan,At,Bt,C,p.*ll,vv.*ll)[inds,inds]))
    dp
end

function matching_indices(current_indices,center_indices)
    current_indices[[in(index, center_indices) for index in current_indices]]
end

function density_pjac!(pjac,p,tt,params)
    plan,iplan,At,Bt,C,vv,f = getindex.(params,[:plan,:iplan,:At,:Bt,:C,:vv,:f])
    nn = size(p)[1]
    nd2 = round(Int,nn/2)
    nd4 = round(Int,nn/4)
    kk = repmat((-nd2:nd2-1)*(1+0im),1,nn)
    ll = copy(kk)'
    inds = (-nd2:(nd2-1))+nn+1

    pjac[:,:] .= -im*2π*Diagonal((f[1]*kk.+f[2]*ll)[:])
    δ_i_j = zeros(eltype(p),nn,nn)
    for j in 1:nn^2
        δ_i_j *= zero(eltype(p))
        δ_i_j[j] = one(eltype(p))

        pjac[:,j] .+= (.-4π^2*(  conv2!(plan,iplan,At,Bt,C,δ_i_j.*(kk.^2 .+ll.^2),tt)[inds,inds]
                                .+conv2!(plan,iplan,At,Bt,C,δ_i_j.*kk,tt.*kk)[inds,inds]
                                .+conv2!(plan,iplan,At,Bt,C,δ_i_j.*ll,tt.*ll)[inds,inds])
                       .-4π^2*(  conv2!(plan,iplan,At,Bt,C,δ_i_j,vv.*(kk.^2 .+ll.^2))[inds,inds]
                               .+conv2!(plan,iplan,At,Bt,C,δ_i_j.*kk,vv.*kk)[inds,inds]
                               .+conv2!(plan,iplan,At,Bt,C,δ_i_j.*ll,vv.*ll)[inds,inds]))[:]
    end
    pjac
end

function density_pjac_approx!(pjac,p,tt,params,filt_matrix)
    plan,iplan,At,Bt,C,vv,f = getindex.(params,[:plan,:iplan,:At,:Bt,:C,:vv,:f])
    nn = size(p)[1]
    nd2 = round(Int,nn/2)
    nd4 = round(Int,nn/4)
    kk = repmat((-nd2:nd2-1)*(1+0im),1,nn)
    ll = copy(kk)'
    inds = (-nd2:(nd2-1))+nn+1

    pjac[:,:] .= -im*2π*Diagonal((f[1]*kk.+f[2]*ll)[:])

    δ = zeros(eltype(p),nn,nn)
    δ[1] = 1
    v1 = (.-4π^2*(  conv2!(plan,iplan,At,Bt,C,δ.*(kk.^2 .+ll.^2),tt)[inds,inds]
                                .+conv2!(plan,iplan,At,Bt,C,δ.*kk,tt.*kk)[inds,inds]
                                .+conv2!(plan,iplan,At,Bt,C,δ.*ll,tt.*ll)[inds,inds])
                       .-4π^2*(  conv2!(plan,iplan,At,Bt,C,δ,vv.*(kk.^2 .+ll.^2))[inds,inds]
                               .+conv2!(plan,iplan,At,Bt,C,δ.*kk,vv.*kk)[inds,inds]
                               .+conv2!(plan,iplan,At,Bt,C,δ.*ll,vv.*ll)[inds,inds]))[:]

    δ[1] = zero(eltype(p))
    δ[end] = 1
    v2 = (.-4π^2*(  conv2!(plan,iplan,At,Bt,C,δ.*(kk.^2 .+ll.^2),tt)[inds,inds]
                                .+conv2!(plan,iplan,At,Bt,C,δ.*kk,tt.*kk)[inds,inds]
                                .+conv2!(plan,iplan,At,Bt,C,δ.*ll,tt.*ll)[inds,inds])
                       .-4π^2*(  conv2!(plan,iplan,At,Bt,C,δ,vv.*(kk.^2 .+ll.^2))[inds,inds]
                               .+conv2!(plan,iplan,At,Bt,C,δ.*kk,vv.*kk)[inds,inds]
                               .+conv2!(plan,iplan,At,Bt,C,δ.*ll,vv.*ll)[inds,inds]))[:]
    v2[1] = v1[1]

    pjac[:,:] .+= (filt_matrix.*ToeplitzMatrices.Toeplitz(v1,v2))[:,:]
 end

function density_ttjac!(ttjac,p,tt,params)
    plan,iplan,At,Bt,C,vv,f = getindex.(params,[:plan,:iplan,:At,:Bt,:C,:vv,:f])
    nn = size(p)[1]
    kk = repmat(-round(Int, nn/2):round(Int,nn/2-1),1,nn)
    kk = convert(Array{eltype(p),2},kk)
    ll = copy(kk)'
    inds = -round(Int,nn/2):round(Int,nn/2-1)
    δ_i_j = zeros(eltype(p),nn,nn)
    for j in 1:nn^2
        # δ_i_j[:,:] .= zeros(eltype(p),nn,nn)
        δ_i_j[:,:] *= zero(eltype(p))
        δ_i_j[j] = one(eltype(p))
        ttjac[:,j] .= (.- 4π^2*(  conv2!(plan,iplan,At,Bt,C,p.*(kk.^2 .+ll.^2),δ_i_j)[inds+nn+1,inds+nn+1]
                                .+conv2!(plan,iplan,At,Bt,C,p.*kk,δ_i_j.*kk)[inds+nn+1,inds+nn+1]
                                .+conv2!(plan,iplan,At,Bt,C,p.*ll,δ_i_j.*ll)[inds+nn+1,inds+nn+1]))[:]
    end
    ttjac
end

function density_ttjac_approx!(ttjac,p,tt,params,filt_matrix)
    plan,iplan,At,Bt,C,vv,f = getindex.(params,[:plan,:iplan,:At,:Bt,:C,:vv,:f])
    nn = size(p)[1]
    nd2 = round(Int,nn/2)
    kk = repmat((-nd2:nd2-1)*(1+0im),1,nn)
    ll = copy(kk)'
    inds = (-nd2:nd2-1)+nn+1

    δ = zeros(eltype(p),nn,nn)
    δ[1] = one(eltype(p))
    v1 = (.- 4π^2*(  conv2!(plan,iplan,At,Bt,C,p.*(kk.^2 .+ll.^2),δ)[inds,inds]
                                .+conv2!(plan,iplan,At,Bt,C,p.*kk,δ.*kk)[inds,inds]
                                .+conv2!(plan,iplan,At,Bt,C,p.*ll,δ.*ll)[inds,inds]))[:]

    δ[1] = zero(eltype(p))
    δ[end] = one(eltype(p))
    v2 = (.- 4π^2*(  conv2!(plan,iplan,At,Bt,C,p.*(kk.^2 .+ll.^2),δ)[inds,inds]
                   .+conv2!(plan,iplan,At,Bt,C,p.*kk,δ.*kk)[inds,inds]
                   .+conv2!(plan,iplan,At,Bt,C,p.*ll,δ.*ll)[inds,inds]))[:]
    v2[1] = v1[1]

    ttjac[:,:] .= (filt_matrix.*ToeplitzMatrices.Toeplitz(v1,v2))[:,:]
end

function temperature!(dtt,tt,p,params)
    plan,iplan,At,Bt,C,vv,A,B,f,T0 = getindex.(params,[:plan,:iplan,:At,:Bt,:C,:vv,:A,:B,:f,:T0])
    nn = size(p)[1]
    nd2 = round(Int,nn/2)
    kk = repmat(-round(Int,nn/2):round(Int,nn/2-1),1,nn)
    ll = copy(kk)'
    delta_k_l_0 = zeros(eltype(p),nn,nn)
    delta_k_l_0[nd2+1,nd2+1] = one(eltype(p))
    inds = (-nd2:nd2-1)+nn+1

    dtt[:, :] .= (  A*(f[1]^2+f[2]^2)*p-2π*im*A*conv2!(plan,iplan,At,Bt,C,(kk*f[1]+ll*f[2]).*p,tt)[inds,inds]
                  - 4π*im*A*conv2!(plan,iplan,At,Bt,C,p,(f[1]*kk+f[2]*ll).*vv)[inds,inds]
                  - 4π^2*A*conv2!(plan,iplan,At,Bt,C,tt,conv2!(plan,iplan,At,Bt,C,p.*kk,vv.*kk)[inds,inds]
                                    +conv2!(plan,iplan,At,Bt,C,p.*ll,vv.*ll)[inds,inds])[inds,inds]
                  - 4π^2*A*conv2!(plan,iplan,At,Bt,C,p,conv2!(plan,iplan,At,Bt,C,vv.*kk,vv.*kk)[inds,inds]
                                   +conv2!(plan,iplan,At,Bt,C,vv.*ll,vv.*ll)[inds,inds])[inds,inds]
                  - 4π^2*B[1]*tt.*(kk.^2 .+ll.^2)
                  .- B[2]*(tt.-T0*delta_k_l_0)
    )
end

function temperature_pjac!(pjac,tt,p,params)
    plan,iplan,At,Bt,C,vv,A,B,f = getindex.(params,[:plan,:iplan,:At,:Bt,:C,:vv,:A,:B,:f])
    nn = size(p)[1]
    nd2 = round(Int,nn/2)
    kk = repmat(-round(Int,nn/2):round(Int,nn/2-1),1,nn)
    ll = copy(kk)'
    delta_k_l_0 = zeros(eltype(p),nn,nn)
    delta_k_l_0[nd2+1,nd2+1] = one(eltype(p))
    inds = (-nd2:nd2-1)+nn+1
    pjac[:,:] .= A*(f[1]^2+f[2]^2)*Diagonal(ones(nn^2))

    δ_i_j = zeros(eltype(p),nn,nn)
    for j in 1:nn^2
        δ_i_j[:,:] .= zeros(eltype(p),nn,nn)
        δ_i_j[j] = one(eltype(p))
        pjac[:,j] .+= ( -2π*im*A*conv2!(plan,iplan,At,Bt,C,(kk*f[1]+ll*f[2]).*δ_i_j,tt)[inds,inds]
                      - 4π*im*A*conv2!(plan,iplan,At,Bt,C,δ_i_j,(f[1]*kk+f[2]*ll).*vv)[inds,inds]
                      - 4π^2*A*conv2!(plan,iplan,At,Bt,C,tt,conv2!(plan,iplan,At,Bt,C,δ_i_j.*kk, vv.*kk)[inds,inds]
                                        +conv2!(plan,iplan,At,Bt,C,δ_i_j.*ll,vv.*ll)[inds,inds])[inds,inds]
                      - 4π^2*A*conv2!(plan,iplan,At,Bt,C,δ_i_j,conv2!(plan,iplan,At,Bt,C,vv.*kk,vv.*kk)[inds,inds]
                                       +conv2!(plan,iplan,At,Bt,C,vv.*ll,vv.*ll)[inds,inds])[inds,inds])[:]
    end
    pjac
end

function temperature_pjac_approx!(pjac,tt,p,params,filt_matrix)
    plan,iplan,At,Bt,C,vv,A,B,f = getindex.(params,[:plan,:iplan,:At,:Bt,:C,:vv,:A,:B,:f])
    nn = size(p)[1]
    nd2 = round(Int,nn/2)
    kk = repmat(-round(Int,nn/2):round(Int,nn/2-1),1,nn)
    ll = copy(kk)'
    delta_k_l_0 = zeros(eltype(p),nn,nn)
    delta_k_l_0[nd2+1,nd2+1] = one(eltype(p))
    inds = (-nd2:nd2-1)+nn+1
    pjac[:,:] .= A*(f[1]^2+f[2]^2)*Diagonal(ones(nn^2))

    δ = zeros(eltype(p),nn,nn)
    δ[1] = one(eltype(p))
    v1 = ( -2π*im*A*conv2!(plan,iplan,At,Bt,C,(kk*f[1].+ll*f[2]).*δ,tt)[inds,inds]
                      .- 4π*im*A*conv2!(plan,iplan,At,Bt,C,δ,(f[1]*kk+f[2]*ll).*vv)[inds,inds]
                      .- 4π^2*A*conv2!(plan,iplan,At,Bt,C,tt,conv2!(plan,iplan,At,Bt,C,δ.*kk, vv.*kk)[inds,inds]
                                        .+conv2!(plan,iplan,At,Bt,C,δ.*ll,vv.*ll)[inds,inds])[inds,inds]
                      .- 4π^2*A*conv2!(plan,iplan,At,Bt,C,δ,conv2!(plan,iplan,At,Bt,C,vv.*kk,vv.*kk)[inds,inds]
                                       .+conv2!(plan,iplan,At,Bt,C,vv.*ll,vv.*ll)[inds,inds])[inds,inds])[:]

    δ[1] *= zero(eltype(p))
    δ[end] = one(eltype(p))
    v2 = ( -2π*im*A*conv2!(plan,iplan,At,Bt,C,(kk*f[1]+ll*f[2]).*δ,tt)[inds,inds]
                      .- 4π*im*A*conv2!(plan,iplan,At,Bt,C,δ,(f[1]*kk+f[2]*ll).*vv)[inds,inds]
                      .- 4π^2*A*conv2!(plan,iplan,At,Bt,C,tt,conv2!(plan,iplan,At,Bt,C,δ.*kk, vv.*kk)[inds,inds]
                                        .+conv2!(plan,iplan,At,Bt,C,δ.*ll,vv.*ll)[inds,inds])[inds,inds]
                      .- 4π^2*A*conv2!(plan,iplan,At,Bt,C,δ,conv2!(plan,iplan,At,Bt,C,vv.*kk,vv.*kk)[inds,inds]
                                       .+conv2!(plan,iplan,At,Bt,C,vv.*ll,vv.*ll)[inds,inds])[inds,inds])[:]
    v2[1] = v1[1]

    pjac[:,:] .+= filt_matrix.*ToeplitzMatrices.Toeplitz(v1,v2)[:,:]
end

function temperature_ttjac!(ttjac,tt,p,params)
    plan,iplan,At,Bt,C,vv,A,B,f = getindex.(params,[:plan,:iplan,:At,:Bt,:C,:vv,:A,:B,:f])
    nn = size(p)[1]
    nd2 = round(Int,nn/2)
    kk = repmat(-round(Int,nn/2):round(Int,nn/2-1),1,nn)
    ll = copy(kk)'
    inds = (-nd2:nd2-1)+nn+1
    delta_k_l_0 = zeros(eltype(p),nn,nn)
    delta_k_l_0[nd2+1,nd2+1] = one(eltype(p))

    ttjac[:,:] .= Diagonal(-4π^2*B[1]*(kk.^2 .+ll.^2)[:]-B[2])
    δ_i_j = zeros(eltype(p),nn,nn)
    for j in 1:nn^2
        δ_i_j[:,:] .= zeros(eltype(p),nn,nn)
        δ_i_j[j] = one(eltype(p))

        ttjac[:,j] .+= ( .-2π*im*A*conv2!(plan,iplan,At,Bt,C,(kk*f[1]+ll*f[2]).*p,δ_i_j)[inds,inds]
                         .-4π^2*A*conv2!(plan,iplan,At,Bt,C,δ_i_j,conv2!(plan,iplan,At,Bt,C,p.*kk,vv.*kk)[inds,inds]
                                           .+conv2!(plan,iplan,At,Bt,C,p.*ll,vv.*ll)[inds,inds])[inds,inds])[:]
    end
    ttjac
end

function temperature_ttjac_approx!(ttjac,tt,p,params,filt_matrix)
    plan,iplan,At,Bt,C,vv,A,B,f = getindex.(params,[:plan,:iplan,:At,:Bt,:C,:vv,:A,:B,:f])
    nn = size(p)[1]
    nd2 = round(Int,nn/2)
    kk = repmat(-round(Int,nn/2):round(Int,nn/2-1),1,nn)
    ll = copy(kk)'
    inds = (-nd2:nd2-1)+nn+1
    delta_k_l_0 = zeros(eltype(p),nn,nn)
    delta_k_l_0[nd2+1,nd2+1] = one(eltype(p))

    ttjac[:,:] .= Diagonal(-4π^2*B[1]*(kk.^2 .+ll.^2)[:]-B[2])

    δ = zeros(eltype(p),nn,nn)
    δ[1] = one(eltype(p))
    v1 = ( .-2π*im*A*conv2!(plan,iplan,At,Bt,C,(kk*f[1].+ll*f[2]).*p,δ)[inds,inds]
           .-4π^2*A*conv2!(plan,iplan,At,Bt,C,δ,conv2!(plan,iplan,At,Bt,C,p.*kk,vv.*kk)[inds,inds]
                           .+conv2!(plan,iplan,At,Bt,C,p.*ll,vv.*ll)[inds,inds])[inds,inds])[:]

    δ[1] = zero(eltype(p))
    δ[end] = one(eltype(p))
    v2 = (.-2π*im*A*conv2!(plan,iplan,At,Bt,C,(kk*f[1].+ll*f[2]).*p,δ)[inds,inds]
           .-4π^2*A*conv2!(plan,iplan,At,Bt,C,δ,conv2!(plan,iplan,At,Bt,C,p.*kk,vv.*kk)[inds,inds]
                           .+conv2!(plan,iplan,At,Bt,C,p.*ll,vv.*ll)[inds,inds])[inds,inds])[:]
    v2[1] = v1[1]

    ttjac[:,:] .+= (filt_matrix.*ToeplitzMatrices.Toeplitz(v1,v2))[:,:]
end

function energy(ϕ, vv, A)
    potential_energy = zero(eltype(ϕ))
    nn = size(vv)[1]
    inds = -round(Int,nn/2):round(Int,nn/2-1)
    tmp = round(Int,nn/2)
    for k in inds
        for l in inds
            for m in inds
                for n in inds
                    (k+m == 0 && n+l == 0) && (potential_energy += ϕ[k+tmp+1,l+tmp+1]*vv[m+tmp+1,n+tmp+1])
                end
            end
        end
    end

    thermal_energy = ϕ[tmp+1, nn+1+tmp]/A

    potential_energy + thermal_energy
end

function real_current(ϕ,params)
    vv = params[:vv]
    nn = size(ϕ)[1]
    nd2 = round(Int,nn/2)
    kk = repmat(-round(Int,nn/2):round(Int,nn/2-1),1,nn)
    ll = copy(kk)'
    p = ϕ[1:nn,1:nn]
    tt = ϕ[1:nn,nn+1:2nn]

    ∇p = -im*2π*[kk.*p,ll.*p]
    ∇tt = -im*2π*[kk.*tt,ll.*tt]
    ∇vv = -im*2π*[kk.*vv,ll.*vv]

    P = fft(fftshift(p))
    T = fft(fftshift(tt))
    ∇P = [fft(fftshift(∇p[1])),fft(fftshift(∇p[2]))]
    ∇T = [fft(fftshift(∇tt[1])),fft(fftshift(∇tt[2]))]
    ∇V = [fft(fftshift(∇vv[1])),fft(fftshift(∇vv[2]))]

    J = -[P.*(∇V[1].-params[:f][1]).+T.*∇P[1],P.*(∇V[2].-params[:f][2]).+T.*∇P[2]]
end

ext: jl


text:
# To add a new cell, type '# %%'
# To add a new markdown cell, type '# %% [markdown]'
# %%
using HDF5
using Plots
using Plots.PlotMeasures
using GR
using QuadGK
using Printf
using LaTeXStrings
using DelimitedFiles
ENV["GKSwstype"] = "100"
# %%
H = 0.25 # height from galactic plane in kpc
rlim = 5 # limit of radius
kmsinkpcyr = 3.1536e7/3.085677581e16 # 1 km/sec in kpc/yr
# timestep = np.linspace(0.01,0.8,80)
time = 300 # snapshot
rootdirectory = "/home/oku/SimulationData/isogal"

models = [
    # ["Osaka2019_isogal", "Osaka2019 (Shimizu+19)"],
        #   , "geodome_model/geodome_original"\
    # ["geodome_model/ver_19.11.1", "Geodesic dome model & Cioffi+ 1988"],
    # ["geodome_model/OKU2020","Geodesic dome model & Athena fitting"],
    # ["centroid_model/ver07271_NoMomCeiling","Centroid model & Athena fitting (alpha = 0)"],
    # ["centroid_model/ver07272_nfb1","Centroid model & Athena fitting (nfb = 1)"],
    # ["centroid_model/ver07272_SFE001","Centroid model & Athena fitting (SFE = 0.01)"],
    # ["centroid_model/ver07311","Centroid model & Athena fitting (alpha = -1)"],
    # ["centroid_model/ver07311_fdens-2","Centroid model & Athena fitting (alpha = -2)"],
    # ["centroid_model/ver08041_alpha-1","Centroid model & Athena fitting (new, alpha = -1)"],
    # ["centroid_model/ver08041_alpha-2","Centroid model & Athena fitting (new, alpha = -2)"],
    # ["centroid_model/ver08041_NoThermal","Centroid model (alpha = -1)"],
    # ["centroid_model/ver08191","Centroid model (alpha = -0.24)"],
    # ["centroid_model/ver07272_CHEVALIER1974","Centroid model & Cioffi+ 1988"],
    # ["centroid_model/ver09011_n1f2a0","Centroid model"],
    # ["centroid_model/ver09011_n1f2a024","Centroid model (new, alpha = -0.24)"],
    # ["centroid_model/ver09011_n1f2a050","Centroid model (new, alpha = -0.50)"],
    # ["centroid_model/ver09011_n1f2a075","Centroid model (new, alpha = -0.75)"],
    # ["centroid_model/ver09011_n1f2a1","Centroid model (new, alpha = -1.00)"],
    # ["centroid_model/ver12151","Centroid"],
    # ["ss_model/ver01092", "Spherical superbubble"],
    # ["ss_model/ver01121", "Spherical superbubble"],
    # ["ss_model/ver01193", "Spherical superbubble"],
    # ["ss_model/ver03221", "Spherical superbubble"],
    # ["ss_model/ver12211", "Spherical superbubble"],
    # ["ss_model/ver03311", "with cool-off"],
    ["ss_model/ver03311/fiducial", "fiducial"],
    ["ss_model/ver03311/med", "high-reso"],
    # ["ss_model/ver03311_NoCoolOff", "without cool-off"],
    # ["NoFB", "No feedback"],
    ]
# %%
function getcenter(sph)
  x = y = z = mtot = 0
  for i in 1:length(sph.mass)
    x += sph.mass[i] * sph.pos[1, i]
    y += sph.mass[i] * sph.pos[2, i]
    z += sph.mass[i] * sph.pos[3, i]
    mtot += sph.mass[i]
  end
  x /= mtot
  y /= mtot
  z /= mtot
  return [x, y, z]
end

function main(galpos)
  Volume = zeros(2)
  Mtot = 0
  for i in 1:length(sph.mass)
    x = sph.pos[1,i] - galpos[1]
    y = sph.pos[2,i] - galpos[2]
    z = sph.pos[3,i] - galpos[3]
    # if (x^2 + y^2 < (rlim + sph.hsml[i])^2) && (z^2 < (H + sph.hsml[i])^2) && x > 0
    if (x^2 + y^2 < rlim^2) && (z^2 < H^2)
      index = 0
      if sph.temp[i] > 1e5
        index = 1 # hot
      else
        index = 2 # cold
      end
      Volume[index] += sph.mass[i]
      # Volume[index] += sph.mass[i]/sph.rho[i]
      Mtot += sph.mass[i]
      if sph.u[i]/sph.press[i] * (5/3 - 1) < 0
        @show sph.u[i]
        @assert 0
      end
    end
  end
  # Vtot = pi*rlim^2*2*H
  Vtot = Volume[1] + Volume[2]
  Volume[1] = Volume[1]/Vtot
  Volume[2] = Volume[2]/Vtot
  @show Volume[1], length(sph.mass)
  return Volume
end
struct Sph
    pos::Array{Float32, 2}
    mass::Array{Float32, 1}
    rho::Array{Float32, 1}
    u::Array{Float32, 1}
    press::Array{Float32, 1}
    temp::Array{Float32, 1}
    hsml::Array{Float32, 1}
end
# %%
Volume = [[[0.0 for j in 1:time] for k in [1, 2]] for i in 1:length(models)]
for j in 1:time
  for i in 1:length(models)
        path2snapshot = @sprintf("%s/%s/snapshot_%03d/snapshot_%03d.hdf5", rootdirectory, models[i][1], j, j)
        h5open(path2snapshot, "r") do fp
            pos = read(fp,"PartType0/Coordinates")
            mass = read(fp, "PartType0/Masses")
            rho = read(fp, "PartType0/Density")
            u = read(fp, "PartType0/InternalEnergy")
            press = read(fp, "PartType0/Pressure")
            temp = read(fp, "PartType0/Temperature")
            hsml = read(fp, "PartType0/SmoothingLength")
            global sph = Sph(pos, mass, rho, u, press, temp, hsml)
        end
        galpos = getcenter(sph)

        result = main(galpos)
        Volume[i][1][j] = result[1]
        Volume[i][2][j] = result[2]
  end
  @show j
end

gr()
font = Plots.font(16)
lfont = Plots.font(12)
Plots.plot(legend=:bottomleft,xlim=(0, time), ylim=(1e-8, 1e-1), yscale=:log10, xlabel="Time [Gyr]", ylabel=L"f_{\rm hot}", xtickfont=font, ytickfont=font, xlabelfont=font, ylabelfont=font, legendfont=lfont, guidefont=font, topmargin=3mm, bottommargin=3mm, rightmargin=3mm)
for i in 1:length(models)
Plots.plot!(Volume[i][1], label="hot phase", markersize=5, markercolor=:red)
Plots.plot!(Volume[i][2], label="cold phase", markersize=5, markercolor=:blue)
end
Plots.savefig("Energyloading.png")
# Plots.savefig("results/Energyloading.pdf")


text:
using PSCOPF

using .PSCOPFFixtures

using Test
using JuMP
using Dates
using DataStructures
using Printf

@testset verbose=true "test_bilevel_pilotables" begin

    TS = [DateTime("2015-01-01T11:00:00")]
    ech = DateTime("2015-01-01T07:00:00")
    next_ech = DateTime("2015-01-01T07:30:00")
    function create_instance(load_1, load_2,
                            limit::Float64=35.,
                            logs=nothing)
        network = PSCOPFFixtures.network_2buses(limit=limit)
        # Pilotables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                                0., 200.,
                                                0., 10.,
                                                Dates.Second(0), Dates.Second(0))
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_2", "prod_2_1", PSCOPF.Networks.PILOTABLE,
                                                0., 200.,
                                                0., 50.,
                                                Dates.Second(0), Dates.Second(0))

        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", TS[1], "S1", load_1)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_2", TS[1], "S1", load_2)

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                    SortedDict{String,PSCOPF.GeneratorState}(), #gen_initial_state
                                    uncertainties, nothing, logs)

        return context
    end

    #=
    TS: [11h, 11h15]
    S: [S1,S2]
                        bus 1                   bus 2
                        |                      |
    (limitable) prod_1_1|       "1_2"          |prod_2_1
    Pmin=0, Pmax=200    |                      |Pmin=0, Pmax=200
    Csta=0, Cprop=10    |                      |Csta=0, Cprop=50
                        |----------------------|
                        |         35           |
                        |                      |
                        |                      |
           load(bus_1)  |                      |load(bus_2)
      S1: 30            |                      | S1: 30

    Market :
        The market uses prod_1_1 to satisfy all the demand
        This causes no RSO constraints

    TSO : does not to take any action
        lol_min = p_global_loss_of_load = 0.
        prod_1_1 : p_imposition_min=0.  p_imposition_max=200 (unchanged bounds)
        prod_2_1 : p_imposition_min=0.  p_imposition_max=200 (unchanged bounds)

    =#
    @testset "no_risk_of_breaking_rso_constraint" begin

        context = create_instance(30., 30.,
                                35.)

        tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(USE_UNITS_PROP_COST_AS_TSO_BOUNDING_COST=false))
        firmness = PSCOPF.compute_firmness(tso,
                                            ech, next_ech,
                                            TS, context)
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

        #TSO RSO constraints are OK
        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

        #Market EOD constraints are OK
        @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

        #TSO does not need to bound pilotables
        #prod_1_1 is not bound
        @test value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"]) < 1e-09
        @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"])
        #prod_2_1 is not bound
        @test value(result.upper.pilotable_model.p_imposition_min["prod_2_1",TS[1],"S1"]) < 1e-09
        @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_2_1",TS[1],"S1"])

        #Market chooses the levels respecting the bounds for pilotable production
        @test 60. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"])
        @test value(result.lower.pilotable_model.p_injected["prod_2_1",TS[1],"S1"]) < 1e-09

        @test value(PSCOPF.get_upper_obj_expr(result)) < 1e-09
        @test value(PSCOPF.get_lower_obj_expr(result)) ≈ (60. * 10.)

    end

    #=
    TS: [11h, 11h15]
    S: [S1,S2]
                        bus 1                   bus 2
                        |                      |
    (limitable) prod_1_1|       "1_2"          |prod_2_1
    Pmin=0, Pmax=200    |                      |Pmin=0, Pmax=200
    Csta=0, Cprop=10    |                      |Csta=0, Cprop=50
                        |----------------------|
                        |         35           |
                        |                      |
                        |                      |
           load(bus_1)  |                      |load(bus_2)
      S1: 10            |                      | S1: 100

    Market :
        market would use prod_1_1 at 110MW cause it is cheaper
        This would cause a RSO constraint (flow of 100 > 35)

    TSO needs to anticipate to prevent that:
        option 1 : oblige prod_1_1 to produce at most 10+35=45 => P_1_1 in [0,45]
            => cost : 200-45=155
        option 2 : oblige prod_2_1 to produce at least 100-35=65 => P_2_1 in [65,200]
            => cost : 65

    Only obliging one of the options + EOD constraint is enough.
    The cheaper option is chosen : option 2
        lol_min = p_global_loss_of_load = 0.
        prod_1_1 : p_imposition_min=0.  p_imposition_max=200 (unchanged bounds)
        prod_2_1 : p_imposition_min=65.  p_imposition_max=200

    Market solution :
        lol = p_loss_of_load = 0.
        prod_1_1 : 45.
        prod_2_1 : 65.
    =#
    @testset "rso_constraint_requires_setting_pilotables_bounds" begin

        context = create_instance(10., 100.,
                                35.)

        tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(USE_UNITS_PROP_COST_AS_TSO_BOUNDING_COST=false))
        firmness = PSCOPF.compute_firmness(tso,
                                            ech, next_ech,
                                            TS, context)
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

        #TSO RSO constraints are OK
        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

        #Market EOD constraints are OK
        @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

        #TSO sets the bounds for pilotable production respecting units' pmin and pmax
        #prod_1_1 is not bound
        @test value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"]) < 1e-09
        @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"])
        #prod_2_1 bounding + EOD, guarantees RSO constraints
        @test 65. ≈ value(result.upper.pilotable_model.p_imposition_min["prod_2_1",TS[1],"S1"])
        @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_2_1",TS[1],"S1"])

        #Market chooses the levels respecting the bounds for pilotable production
        @test 45. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"])
        @test 65. ≈ value(result.lower.pilotable_model.p_injected["prod_2_1",TS[1],"S1"])

        @test value(PSCOPF.get_upper_obj_expr(result)) ≈ ( 65 )
        @test value(PSCOPF.get_lower_obj_expr(result)) ≈ ( 45*10 + 65*50. )

    end

    #=
    TS: [11h, 11h15]
    S: [S1,S2]
                        bus 1                   bus 2
                        |                      |
    (limitable) prod_1_1|       "1_2"          |prod_2_1
    Pmin=0, Pmax=200    |                      |Pmin=0, Pmax=200
    Csta=0, Cprop=10    |                      |Csta=0, Cprop=50
                        |----------------------|
                        |         500          |
                        |                      |
                        |                      |
           load(bus_1)  |                      |load(bus_2)
      S1: 400           |                      | S1: 20

    Market :
        Possible prod : 400,  Demand : 420
        market will use prod_1_1 and prod_2_1 (=> 400MW)
        market will need to cut 20MW
        This will not cause RSO constraints violation

    TSO does not need to intervene
        lol_min = p_global_loss_of_load = 0.
        prod_1_1 : p_imposition_min=0.  p_imposition_max=200 (unchanged bounds)
        prod_2_1 : p_imposition_min=0.  p_imposition_max=200 (unchanged bounds)

        The TSO will locate the LoL on either of bus 1 or 2
        p_loss_of_load[bus1] + p_loss_of_load[bus2] = 20

    Market solution :
        lol = p_loss_of_load = 20.
        prod_1_1 : 200.
        prod_2_1 : 200.

    =#
    @testset "eod_problem_requires_cutting_conso" begin

        context = create_instance(400., 20.,
                                500.)

        tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(USE_UNITS_PROP_COST_AS_TSO_BOUNDING_COST=false))
        firmness = PSCOPF.compute_firmness(tso,
                                            ech, next_ech,
                                            TS, context)
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK

        #TSO RSO constraints can be respected without capping or losing load OK
        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

        #Market needs to cut conso to assure EOD constraint
        @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test 20. ≈ value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"])

        # RSO locates cut conso : here any combination will do
        @test 20. ≈ ( value(result.upper.lol_model.p_loss_of_load["bus_1", TS[1],"S1"])
                    + value(result.upper.lol_model.p_loss_of_load["bus_2", TS[1],"S1"]) )

        #TSO does not need to bound pilotables
        #prod_1_1 is not bound
        @test value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"]) < 1e-09
        @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"])
        #prod_2_1 is not bound
        @test value(result.upper.pilotable_model.p_imposition_min["prod_2_1",TS[1],"S1"]) < 1e-09
        @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_2_1",TS[1],"S1"])

        #Market chooses the levels respecting the bounds for pilotable production
        @test 200. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"])
        @test 200. ≈ value(result.lower.pilotable_model.p_injected["prod_2_1",TS[1],"S1"])

        @test value(PSCOPF.get_upper_obj_expr(result)) ≈ (20 *tso.configs.TSO_LOL_PENALTY)
        @test value(PSCOPF.get_lower_obj_expr(result)) ≈ ( 200*10. + 200*50.
                                                        + 20 *tso.configs.MARKET_LOL_PENALTY)

    end

    #=
    TS: [11h, 11h15]
    S: [S1,S2]
                        bus 1                   bus 2
                        |                      |
    (limitable) prod_1_1|                      |prod_2_1
    Pmin=0, Pmax=200    |                      |Pmin=0, Pmax=200
    Csta=0, Cprop=10    |        "1_2"         |Csta=0, Cprop=50
                        |----------------------|
                        |         180          |
                        |                      |
                        |                      |
           load(bus_1)  |                      |load(bus_2)
      S1: 400           |                      | S1: 20

    Market :
        Possible prod : 400,  Demand : 420
        market will use prod_1_1 and prod_2_1 (=> 400MW)
        market will need to cut 20MW

    If load is cut from bus 1 => no problem
    If load is cut from bus 2 => RSO constraint violation

    TSO needs to prevent this by locating the cut conso on bus 2
        option 1 : limit prod2 to 0-20
        lol_min = p_global_loss_of_load = 0.
        prod_1_1 : p_imposition_min=0.  p_imposition_max=200 (unchanged bounds)
        prod_2_1 : p_imposition_min=0.  p_imposition_max=200 (unchanged bounds)

        The TSO will locate the LoL on bus 1 to avoid the constraint violation
        p_loss_of_load[bus1] = 20
        p_loss_of_load[bus2] = 0

    Market solution :
        lol = p_loss_of_load = 20.
        prod_1_1 : 200.
        prod_2_1 : 200.

    =#
    @testset "eod_and_rso_problems_require_cutting_conso_and_locating_it" begin

        context = create_instance(400., 20.,
                                180.)

        tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(USE_UNITS_PROP_COST_AS_TSO_BOUNDING_COST=false))
        firmness = PSCOPF.compute_firmness(tso,
                                            ech, next_ech,
                                            TS, context)
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK

        #TSO RSO constraints can be respected without capping or losing load OK
        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

        #Market needs to cut conso to assure EOD constraint
        @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test 20. ≈ value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"])

        # RSO locates cut conso : due to flow limit, conso must be cut on bus_1
        @test 20. ≈ value(result.upper.lol_model.p_loss_of_load["bus_1", TS[1],"S1"])
        @test value(result.upper.lol_model.p_loss_of_load["bus_2", TS[1],"S1"]) < 1e-09

        #TSO does not need to bound pilotables
        #prod_1_1 is not bound
        @test value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"]) < 1e-09
        @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"])
        #prod_2_1 is not bound
        @test value(result.upper.pilotable_model.p_imposition_min["prod_2_1",TS[1],"S1"]) < 1e-09
        @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_2_1",TS[1],"S1"])

        #Market chooses the levels respecting the bounds for pilotable production
        @test 200. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"])
        @test 200. ≈ value(result.lower.pilotable_model.p_injected["prod_2_1",TS[1],"S1"])

        @test value(PSCOPF.get_upper_obj_expr(result)) ≈ (20 * tso.configs.TSO_LOL_PENALTY) #LoL due to Market

    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
"""
    fit(estimator, X, y[, test...]; [verbosity = 1, is_row_major = false])

Fit the `estimator` with features data `X` and label `y` using the X-y pairs in `test` as
validation sets.

Return a dictionary with an entry for each validation set. Each entry of the dictionary is another
dictionary with an entry for each validation metric in the `estimator`. Each of these entries is an
array that holds the validation metric's value at each iteration.

# Arguments
* `estimator::LGBMEstimator`: the estimator to be fit.
* `X::Matrix{TX<:Real}`: the features data.
* `y::Vector{Ty<:Real}`: the labels.
* `test::Tuple{Matrix{TX},Vector{Ty}}...`: optionally contains one or more tuples of X-y pairs of
    the same types as `X` and `y` that should be used as validation sets.
* `verbosity::Integer`: keyword argument that controls LightGBM's verbosity. `< 0` for fatal logs
    only, `0` includes warning logs, `1` includes info logs, and `> 1` includes debug logs.
* `is_row_major::Bool`: keyword argument that indicates whether or not `X` is row-major. `true`
    indicates that it is row-major, `false` indicates that it is column-major (Julia's default).
* `weights::Vector{Tw<:Real}`: the training weights.
* `init_score::Vector{Ti<:Real}`: the init scores.
"""
function fit(estimator::LGBMEstimator, X::Matrix{TX},
    y::Vector{Ty}, test::Tuple{Matrix{TX},Vector{Ty}}...; verbosity::Integer = 1,
    is_row_major = false, weights::Vector{Tw} = Vector{Float32}(),
    init_score::Vector{Ti} = Vector{Float64}()) where {TX<:Real,Ty<:Real,Tw<:Real,Ti<:Real}

    start_time = now()

    log_debug(verbosity, "Started creating LGBM training dataset\n")
    ds_parameters = stringifyparams(estimator, DATASETPARAMS)
    train_ds = LGBM_DatasetCreateFromMat(X, ds_parameters, is_row_major)
    LGBM_DatasetSetField(train_ds, "label", y)
    if length(weights) > 0
        LGBM_DatasetSetField(train_ds, "weight", weights)
    end
    if length(init_score) > 0
        LGBM_DatasetSetField(train_ds, "init_score", init_score)
    end

    log_debug(verbosity, "Started creating LGBM booster\n")
    bst_parameters = stringifyparams(estimator, BOOSTERPARAMS) * " verbosity=$verbosity"
    estimator.booster = LGBM_BoosterCreate(train_ds, bst_parameters)

    n_tests = length(test)
    tests_names = Array{String}(undef,n_tests)

    if n_tests > 0
        log_debug(verbosity, "Started creating LGBM test datasets\n")
        @inbounds for (test_idx, test_entry) in enumerate(test)
            tests_names[test_idx] = "test_$(test_idx)"
            test_ds = LGBM_DatasetCreateFromMat(test_entry[1], ds_parameters, train_ds, is_row_major)
            LGBM_DatasetSetField(test_ds, "label", test_entry[2])
            LGBM_BoosterAddValidData(estimator.booster, test_ds)
        end
    end

    log_debug(verbosity, "Started training...\n")
    results = train(estimator, tests_names, verbosity, start_time)
    # estimator.model = readlines("$(tempdir)/model.txt")

    return results
end

function train(estimator::LGBMEstimator, tests_names::Vector{String}, verbosity::Integer,
               start_time::DateTime)
    results = Dict{String,Dict{String,Vector{Float64}}}()
    n_tests = length(tests_names)
    metrics = LGBM_BoosterGetEvalNames(estimator.booster)
    n_metrics = length(metrics)
    bigger_is_better = [ifelse(in(metric, MAXIMIZE_METRICS), 1., -1.) for metric in metrics]
    best_score = fill(-Inf, (n_metrics, n_tests))
    best_iter = fill(1, (n_metrics, n_tests))

    for iter in 1:estimator.num_iterations
        is_finished = LGBM_BoosterUpdateOneIter(estimator.booster)
        log_debug(verbosity, Dates.CompoundPeriod(now() - start_time),
                  " elapsed, finished iteration ", iter, "\n")
        if is_finished == 0
            is_finished = eval_metrics!(results, estimator, tests_names, iter, n_metrics,
                                        verbosity, bigger_is_better, best_score, best_iter, metrics)
        else
            shrinkresults!(results, iter - 1)
            log_info(verbosity, "Stopped training because there are no more leaves that meet the ",
                     "split requirements.")
        end
        is_finished == 1 && return results
    end
    return results
end

function eval_metrics!(results::Dict{String,Dict{String,Vector{Float64}}},
                       estimator::LGBMEstimator, tests_names::Vector{String}, iter::Integer,
                       n_metrics::Integer, verbosity::Integer,
                       bigger_is_better::Vector{Float64}, best_score::Matrix{Float64},
                       best_iter::Matrix{Int}, metrics::Vector{String})
    if (iter - 1) % estimator.metric_freq == 0
        if estimator.is_training_metric
            scores = LGBM_BoosterGetEval(estimator.booster, 0)
            store_scores!(results, estimator, iter, "training", scores, metrics)
            print_scores(estimator, iter, "training", n_metrics, scores, metrics, verbosity)
        end
    end

    if (iter - 1) % estimator.metric_freq == 0 || estimator.early_stopping_round > 0
        for (test_idx, test_name) in enumerate(tests_names)
            scores = LGBM_BoosterGetEval(estimator.booster, test_idx)

            # Check if progress should be stored and/or printed
            if (iter - 1) % estimator.metric_freq == 0
                store_scores!(results, estimator, iter, test_name, scores, metrics)
                print_scores(estimator, iter, test_name, n_metrics, scores, metrics, verbosity)
            end

            # Check if early stopping is called for
            @inbounds for metric_idx in eachindex(metrics)
                maximize_score = bigger_is_better[metric_idx] * scores[metric_idx]
                if maximize_score > best_score[metric_idx, test_idx]
                    best_score[metric_idx, test_idx] = maximize_score
                    best_iter[metric_idx, test_idx] = iter
                elseif iter - best_iter[metric_idx, test_idx] >= estimator.early_stopping_round
                    shrinkresults!(results, best_iter[metric_idx, test_idx])
                    log_info(verbosity, "Early stopping at iteration ", iter,
                             ", the best iteration round is ", best_iter[metric_idx, test_idx], "\n")
                    return 1
                end
            end
        end
    end

    return 0
end

function store_scores!(results::Dict{String,Dict{String,Vector{Float64}}},
                       estimator::LGBMEstimator, iter::Integer, evalname::String,
                       scores::Vector{Cdouble}, metrics::Vector{String})
    for (metric_idx, metric_name) in enumerate(metrics)
        if !haskey(results, evalname)
            num_evals = cld(estimator.num_iterations, estimator.metric_freq)
            results[evalname] = Dict{String,Vector{Float64}}()
            results[evalname][metric_name] = Array{Float64}(undef,num_evals)
        elseif !haskey(results[evalname], metric_name)
            num_evals = cld(estimator.num_iterations, estimator.metric_freq)
            results[evalname][metric_name] = Array{Float64}(undef,num_evals)
        end
        eval_idx = cld(iter, estimator.metric_freq)
        results[evalname][metric_name][eval_idx] = scores[metric_idx]
    end

    return nothing
end

function print_scores(estimator::LGBMEstimator, iter::Integer, name::String, n_metrics::Integer,
                      scores::Vector{Cdouble}, metrics::Vector{String}, verbosity::Integer)
    log_info(verbosity, "Iteration: ", iter, ", ", name, "'s ")
    for (metric_idx, metric_name) in enumerate(metrics)
        log_info(verbosity, metric_name, ": ", scores[metric_idx])
        metric_idx < n_metrics && log_info(verbosity, ", ")
    end
    log_info(verbosity, "\n")
end

function stringifyparams(estimator::LGBMEstimator, params::Vector{Symbol})
    paramstring = ""
    n_params = length(params)
    valid_names = fieldnames(typeof(estimator))
    for (param_idx, param_name) in enumerate(params)
        if in(param_name, valid_names)
            param_value = getfield(estimator, param_name)
            if !isempty(param_value)
                # Convert parameters that contain indices to C's zero-based indices.
                if in(param_name, INDEXPARAMS)
                    param_value -= 1
                end

                if typeof(param_value) <: Array
                    n_entries = length(param_value)
                    if n_entries >= 1
                        paramstring = string(paramstring, param_name, "=", param_value[1])
                        for entry_idx in 2:n_entries
                            paramstring = string(paramstring, ",", param_value[entry_idx])
                        end
                        paramstring = string(paramstring, " ")
                    end
                else
                    paramstring = string(paramstring, param_name, "=", param_value, " ")
                end
            end
        end
    end
    return paramstring[1:end - 1]
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    IdDict([itr])

`IdDict{K,V}()` constructs a hash table using [`objectid`](@ref) as hash and
`===` as equality with keys of type `K` and values of type `V`.

See [`Dict`](@ref) for further help. In the example below, The `Dict`
keys are all `isequal` and therefore get hashed the same, so they get overwritten.
The `IdDict` hashes by object-id, and thus preserves the 3 different keys.

# Examples
```julia-repl
julia> Dict(true => "yes", 1 => "no", 1.0 => "maybe")
Dict{Real, String} with 1 entry:
  1.0 => "maybe"

julia> IdDict(true => "yes", 1 => "no", 1.0 => "maybe")
IdDict{Any, String} with 3 entries:
  true => "yes"
  1.0  => "maybe"
  1    => "no"
```
"""
mutable struct IdDict{K,V} <: AbstractDict{K,V}
    ht::Vector{Any}
    count::Int
    ndel::Int
    IdDict{K,V}() where {K, V} = new{K,V}(Vector{Any}(undef, 32), 0, 0)

    function IdDict{K,V}(itr) where {K, V}
        d = IdDict{K,V}()
        for (k,v) in itr; d[k] = v; end
        d
    end

    function IdDict{K,V}(pairs::Pair...) where {K, V}
        d = IdDict{K,V}()
        sizehint!(d, length(pairs))
        for (k,v) in pairs; d[k] = v; end
        d
    end

    IdDict{K,V}(d::IdDict{K,V}) where {K, V} = new{K,V}(copy(d.ht), d.count, d.ndel)
end

IdDict() = IdDict{Any,Any}()
IdDict(kv::Tuple{}) = IdDict()

IdDict(ps::Pair{K,V}...)           where {K,V} = IdDict{K,V}(ps)
IdDict(ps::Pair{K}...)             where {K}   = IdDict{K,Any}(ps)
IdDict(ps::(Pair{K,V} where K)...) where {V}   = IdDict{Any,V}(ps)
IdDict(ps::Pair...)                            = IdDict{Any,Any}(ps)

function IdDict(kv)
    try
        dict_with_eltype((K, V) -> IdDict{K, V}, kv, eltype(kv))
    catch
        if !applicable(iterate, kv) || !all(x->isa(x,Union{Tuple,Pair}),kv)
            throw(ArgumentError(
                "IdDict(kv): kv needs to be an iterator of tuples or pairs"))
        else
            rethrow()
        end
    end
end

empty(d::IdDict, ::Type{K}, ::Type{V}) where {K, V} = IdDict{K,V}()

function rehash!(d::IdDict, newsz = length(d.ht))
    d.ht = ccall(:jl_idtable_rehash, Vector{Any}, (Any, Csize_t), d.ht, newsz)
    d
end

function sizehint!(d::IdDict, newsz)
    newsz = _tablesz(newsz*2)  # *2 for keys and values in same array
    oldsz = length(d.ht)
    # grow at least 25%
    if newsz < (oldsz*5)>>2
        return d
    end
    rehash!(d, newsz)
end

function setindex!(d::IdDict{K,V}, @nospecialize(val), @nospecialize(key)) where {K, V}
    !isa(key, K) && throw(ArgumentError("$(limitrepr(key)) is not a valid key for type $K"))
    if !(val isa V) # avoid a dynamic call
        val = convert(V, val)
    end
    if d.ndel >= ((3*length(d.ht))>>2)
        rehash!(d, max(length(d.ht)>>1, 32))
        d.ndel = 0
    end
    inserted = RefValue{Cint}(0)
    d.ht = ccall(:jl_eqtable_put, Array{Any,1}, (Any, Any, Any, Ptr{Cint}), d.ht, key, val, inserted)
    d.count += inserted[]
    return d
end

function get(d::IdDict{K,V}, @nospecialize(key), @nospecialize(default)) where {K, V}
    val = ccall(:jl_eqtable_get, Any, (Any, Any, Any), d.ht, key, default)
    val === default ? default : val::V
end

function getindex(d::IdDict{K,V}, @nospecialize(key)) where {K, V}
    val = ccall(:jl_eqtable_get, Any, (Any, Any, Any), d.ht, key, secret_table_token)
    val === secret_table_token && throw(KeyError(key))
    return val::V
end

function pop!(d::IdDict{K,V}, @nospecialize(key), @nospecialize(default)) where {K, V}
    found = RefValue{Cint}(0)
    val = ccall(:jl_eqtable_pop, Any, (Any, Any, Any, Ptr{Cint}), d.ht, key, default, found)
    if found[] === Cint(0)
        return default
    else
        d.count -= 1
        d.ndel += 1
        return val::V
    end
end

function pop!(d::IdDict{K,V}, @nospecialize(key)) where {K, V}
    val = pop!(d, key, secret_table_token)
    val === secret_table_token && throw(KeyError(key))
    return val::V
end

function delete!(d::IdDict{K}, @nospecialize(key)) where K
    pop!(d, key, secret_table_token)
    d
end

function empty!(d::IdDict)
    resize!(d.ht, 32)
    ccall(:memset, Ptr{Cvoid}, (Ptr{Cvoid}, Cint, Csize_t), d.ht, 0, sizeof(d.ht))
    d.ndel = 0
    d.count = 0
    return d
end

_oidd_nextind(a, i) = reinterpret(Int, ccall(:jl_eqtable_nextind, Csize_t, (Any, Csize_t), a, i))

function iterate(d::IdDict{K,V}, idx=0) where {K, V}
    idx = _oidd_nextind(d.ht, idx)
    idx == -1 && return nothing
    return (Pair{K, V}(d.ht[idx + 1]::K, d.ht[idx + 2]::V), idx + 2)
end

length(d::IdDict) = d.count

copy(d::IdDict) = typeof(d)(d)

function get!(d::IdDict{K,V}, @nospecialize(key), @nospecialize(default)) where {K, V}
    val = ccall(:jl_eqtable_get, Any, (Any, Any, Any), d.ht, key, secret_table_token)
    if val === secret_table_token
        val = isa(default, V) ? default : convert(V, default)
        setindex!(d, val, key)
        return val
    else
        return val::V
    end
end

function get(default::Callable, d::IdDict{K,V}, @nospecialize(key)) where {K, V}
    val = ccall(:jl_eqtable_get, Any, (Any, Any, Any), d.ht, key, secret_table_token)
    if val === secret_table_token
        return default()
    else
        return val::V
    end
end

function get!(default::Callable, d::IdDict{K,V}, @nospecialize(key)) where {K, V}
    val = ccall(:jl_eqtable_get, Any, (Any, Any, Any), d.ht, key, secret_table_token)
    if val === secret_table_token
        val = default()
        if !isa(val, V)
            val = convert(V, val)::V
        end
        setindex!(d, val, key)
        return val
    else
        return val::V
    end
end

in(@nospecialize(k), v::KeySet{<:Any,<:IdDict}) = get(v.dict, k, secret_table_token) !== secret_table_token

# For some AbstractDict types, it is safe to implement filter!
# by deleting keys during iteration.
filter!(f, d::IdDict) = filter_in_one_pass!(f, d)

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
# readme example

addprocs(2)
using Caesar, RoME, TransformUtils, Distributions

# load scene and ROV model (might experience UDP packet loss LCM buffer not set)
vc = startdefaultvisualization()
sc1 = loadmodel(:scene01); sc1(vc)
rovt = loadmodel(:rov); rovt(vc)

initCov = 0.001*eye(6); [initCov[i,i] = 0.00001 for i in 4:6];
odoCov = 0.0001*eye(6); [odoCov[i,i] = 0.00001 for i in 4:6];
rangecov, bearingcov = 3e-4, 2e-3

# start and add to a factor graph
fg = identitypose6fg(initCov=initCov)
tf = SE3([0.0;0.7;0.0], Euler(pi/4,0.0,0.0) )
addOdoFG!(fg, Pose3Pose3(MvNormal(veeEuler(tf), odoCov) ) )

visualizeallposes!(vc, fg, drawlandms=false)

addLinearArrayConstraint(fg, (4.0, 0.0), :x0, :l1, rangecov=rangecov,bearingcov=bearingcov)
visualizeDensityMesh!(vc, fg, :l1)
addLinearArrayConstraint(fg, (4.0, 0.0), :x1, :l1, rangecov=rangecov,bearingcov=bearingcov)

solveTree!(fg)

visualize(fg, vc, drawlandms=true, densitymeshes=[:l1;:x2])

#


text:
#rump's formula

using Unums
include("../methods/unum_optimizer.jl")

vtypes = [Float16, Float32, Float64, BigFloat]

function type_rump(T, x, y)
  A = T(333.75)
  B = T(11)
  C = T(121)
  D = T(2)
  E = T(5.5)

  A * y^6 + (x^2) * (B * x^2 * y^2 - y^6 - C * y^4 - D) + E * y^8 + x / (D * y)
end

for T in vtypes
  x = convert(T, 77617)
  y = convert(T, 33096)

  z = Float64(type_rump(T, x, y))

  println("rump's formula result in $(T) : $(z)")
end

function u_rump(T)
  x = T(77617)
  y = T(33096)

  type_rump(T, x, y)
end

r = optimize(u_rump, 0.1, 1, 1, verbose = true)

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
function validate_poly(p)
    if first(p) != last(p)
        error("Polygon should have first and last elements equal")
    end
end


text:
abstract type Shaft <: GeneratorComponent end

@def shaft_ports begin
    state_input = Vector{Int64}()
    inner_input = [τe_var, τm_var]
end

"""
Parameters of single mass shaft model. Typically represents the rotor mass.

# Conmutable structor
```julia
SingleMass(H, D)
```

# Arguments
* `H`   ::Float64 : Rotor inertia constant in MWs/MVA
* `D`   ::Float64 : Rotor natural damping in pu

"""
mutable struct SingleMass <: Shaft
    H::Float64
    D::Float64
    n_states::Int64
    states::Vector{Symbol}
    ports::Ports

    function SingleMass(H::Float64,
                        D::Float64)

        n_states = 2
        states = [:δ, :ω]
        @shaft_ports

        new(H,
            D,
            n_states,
            states,
            Ports(state_input, inner_input))
    end

end

"""
Parameters of 5 mass-spring shaft model.
It contains a High-Pressure (HP) steam turbine, Intermediate-Pressure (IP)
steam turbine, Low-Pressure (LP) steam turbine, the Rotor and
an Exciter (EX) mover.

# Conmutable structor
```julia
FiveMassShaft(H, H_hp, H_ip, H_lp, H_ex,
              D, D_hp, D_ip, D_lp, D_ex,
              D_12, D_23, D_34, D_45,
              K_hp, K_ip, K_lp, K_ex)
```

# Arguments
* `H`   ::Float64 : Rotor inertia constant in MWs/MVA
* `H_hp`::Float64 : High pressure turbine inertia constant in MWs/MVA
* `H_ip`::Float64 : Intermediate pressure turbine inertia constant in MWs/MVA
* `H_lp`::Float64 : Low pressure turbine inertia constant in MWs/MVA
* `H_ex`::Float64 : Exciter inertia constant in MWs/MVA
* `D`   ::Float64 : Rotor natural damping in pu
* `D_hp`::Float64 : High pressure turbine natural damping in pu
* `D_ip`::Float64 : Intermediate pressure turbine natural damping in pu
* `D_lp`::Float64 : Low pressure turbine natural damping in pu
* `D_ex`::Float64 : Exciter natural damping in pu
* `D_12`::Float64 : High-Intermediate pressure turbine damping
* `D_23`::Float64 : Intermediate-Low pressure turbine damping
* `D_34`::Float64 : Low pressure turbine-Rotor damping
* `D_45`::Float64 : Rotor-Exciter damping
* `K_hp`::Float64 : High pressure turbine angle coefficient
* `K_ip`::Float64 : Intermediate pressure turbine angle coefficient
* `K_lp`::Float64 : Low pressure turbine angle coefficient
* `K_ex`::Float64 : Exciter angle coefficient

"""
mutable struct FiveMassShaft <: Shaft
    H::Float64
    H_hp::Float64
    H_ip::Float64
    H_lp::Float64
    H_ex::Float64
    D::Float64
    D_hp::Float64
    D_ip::Float64
    D_lp::Float64
    D_ex::Float64
    D_12::Float64
    D_23::Float64
    D_34::Float64
    D_45::Float64
    K_hp::Float64
    K_ip::Float64
    K_lp::Float64
    K_ex::Float64
    n_states::Int64
    states::Vector{Symbol}
    ports::Ports

    function FiveMassShaft(H::Float64,
                           H_hp::Float64,
                           H_ip::Float64,
                           H_lp::Float64,
                           H_ex::Float64,
                           D::Float64,
                           D_hp::Float64,
                           D_ip::Float64,
                           D_lp::Float64,
                           D_ex::Float64,
                           D_12::Float64,
                           D_23::Float64,
                           D_34::Float64,
                           D_45::Float64,
                           K_hp::Float64,
                           K_ip::Float64,
                           K_lp::Float64,
                           K_ex::Float64)

        n_states = 10
        states = [:δ, :ω,
                  :δ_hp, :ω_hp,
                  :δ_ip, :ω_ip,
                  :δ_lp, :ω_lp,
                  :δ_ex, :ω_ex]
        @shaft_ports

        new(H,
            H_hp,
            H_ip,
            H_lp,
            H_ex,
            D,
            D_hp,
            D_ip,
            D_lp,
            D_ex,
            D_12,
            D_23,
            D_34,
            D_45,
            K_hp,
            K_ip,
            K_lp,
            K_ex,
            n_states,
            states,
            Ports(state_input, inner_input))
    end

end


text:
using Aqua
using LogicCircuits
using Test

@testset "Aqua tests" begin
    Aqua.test_all(LogicCircuits,
                    ambiguities = false,
                    unbound_args= false)
end


text:
mutable struct Queue
    c::Int64 # number of servers
    adist::Distribution # arrival distribution
    sdist::Distribution # service distribution (single server)
end
Queue() = Queue(2, Exponential(1/λ), Exponential(1/μ))

mutable struct Sim
    ncust::Int64 # total number of customers generated
    timelimit::Int64 # time limit for simulation
    seed::Int64 # rng seed
end
Sim() = Sim(1_000, 10_000, 8710)

# define customer behavior
@resumable function customer(env::Environment, server::Resource, id::Integer, t_a::Float64, d_s::Distribution, starttime::Vector{Float64}, departtime::Vector{Float64})
    @yield timeout(env, t_a) # customer arrives
    # println("Customer $id arrived: ", now(env))
    @yield request(server) # customer starts service
    # println("Customer $id entered service: ", now(env))
    starttime[id] = now(env)
    @yield timeout(env, rand(d_s)) # server is busy
    @yield release(server) # customer exits service
    # println("Customer $id exited service: ", now(env))
    departtime[id] = now(env)
end

function runsim(s::Sim, q::Queue)
    ncust = s.ncust
    seed = s.seed
    c = q.c
    adist = q.adist
    sdist = q.sdist
    srand(seed) # set rng seed
    # output of queue simulation
    # arrive order, arrive time, start service time, depart order, depart time
    arriveorder = collect(1:ncust)
    arrivetime = cumsum(rand(adist, ncust))
    starttime = zeros(ncust)
    departtime = zeros(ncust)
    # setup and run simulation
    queuesim = Simulation() # initialize simulation environment
    server = Resource(queuesim, c) # initialize servers
    for i = 1:ncust # initialize customers
        @process customer(queuesim, server, i, arrivetime[i], sdist, starttime, departtime)
    end
    run(queuesim) # run simulation
    depart = hcat(sortrows(hcat(departtime, arriveorder)), collect(1:ncust))
    departorder = convert(Vector{Int64}, sortrows(depart, by = x -> (x[2],x[1],x[3]))[:,3])
    DataFrame(aorder = arriveorder, atime = arrivetime, stime = starttime, dorder = departorder, dtime = departtime)
end

# Disorder Function
# Create noisy data for robust methods
# Input: queue sim output
# Output: true and noisy arrival and departure data
function disorder(simout::DataFrame, num_servers::Int64)
  # actual data from simulation output
  A = simout[:atime] # arrive time
  D_num = simout[:dorder] # depart_order
  D_sort = simout[:dtime] # depart_time
  D_raw = sort(D_sort) # sorted depart_time
  # constants and data structures
  max_custs = length(A)
  results = DataFrame()
  # Introduce partial observation error on cust depart order and time
  D_num_err = zeros(max_custs)
  D_num_err[1:max_custs] = D_num[1:max_custs] # order error
  D_err = zeros(max_custs) # misordered depart times
  D_err[1:max_custs] = D_sort[1:max_custs] # misordered depart times
  for i = 1:(max_custs - num_servers - 1)
    if rand() < 0.1
      swap = i + rand(1:num_servers)
      temp1 = D_num_err[i]
      temp2 = D_err[i]
      D_num_err[i] = D_num_err[swap]
      D_num_err[swap] = temp1
      D_err[i] = D_err[swap]
      D_err[swap] = temp2
    end
  end
  len = size(D_raw, 1)
  D_raw_meas = zeros(len) # measure error in depart times (original order)
  D_raw_meas[1] = D_raw[1] / 2 + rand() * (mean(D_raw[1:2]) - D_raw[1] / 2)
  D_raw_meas[len] = D_raw[len]
  for i = 2:(len - 1)
    a = mean(D_raw[(i-1):i])
    b = mean(D_raw[i:(i+1)])
    D_raw_meas[i] = a + rand() * (b - a)
  end
  D_meas = zeros(len)
  incr = 0
  for ind in simout[:dorder]
      incr += 1
      D_meas[incr] = D_raw_meas[ind]
  end
  DataFrame(ArriveOrder = 1:max_custs, Arrive = A[1:max_custs],
    Depart = D_sort[1:max_custs], DepartMeas = D_meas[1:max_custs],
    DepartErr = D_err[1:max_custs], DepartOrder = D_num[1:max_custs],
    DepartOrderErr = D_num_err[1:max_custs])
end


text:
# Note: we use the 100s place of the enum to represent the number of dimensions
# And the 10s place to represent which dimension is the observation (N), 0 for no observation dimension
@enum ImageFormat begin
    L=100;
    CL=200; LC; HW; WH; NL=210; LN=220;
    HWC=300; WHC; CHW; CWH; NCL=310; NHW; NWH; NLC; CLN=330; LCN; HWN; WHN;
    NHWC=410; NWHC; NCHW; NCWH; HWCN=440; WHCN; CHWN; CWHN;
end

expected_ndims(x::ImageFormat) = Int(x) ÷ 100 # get the 100's digit
obs_dim(x::ImageFormat) = (Int(x) % 100) ÷ 10 # get the 10's digit
const strip_obs = Dict(
    NL=>L, LN=>L, NCL=>CL, CLN=>CL, NLC=>LC, LCN=>LC,
    NHW=>HW, HWN=>HW, NWH=>WH, WHN=>WH,
    NHWC=>HWC, NWHC=>WHC, NCHW=>CHW, NCWH=>CWH, HWCN=>HWC, WHCN=>WHC, CHWN=>CHW, CWHN=>CWH
)

"""
    log_images(logger::TBLogger, name::AbstractString, imgArrays::AbstractArray, format::ImageFormat; step = nothing)

Log multiple images using `Array` of images and format
- imgArrays: `Array` of images, e.g. Array{Array{Float64, 3}, 1}. `Array` of images can be multidimensional.
- format: format which applies to each image in the `Array` of images. It can be one of
  {L, LN, NL,
   CL, LC, NCL, NLC, CLN, LCN,
   HW, WH, HWC, WHC, CHW, CWH,
   HWN, WHN, NHW, NWH,
   HWCN, WHCN, CHWN, CWHN, NHWC, NWHC, NCHW, NCWH}
"""
function log_images(logger::TBLogger, name::AbstractString, imgArrays::AbstractArray; step = nothing)
    @assert isa(first(imgArrays), AbstractArray{<:Colorant}) "Please specify format"
    for (n, imgArray) in enumerate(imgArrays)
        log_image(logger, name*"/$n", imgArray, step = step)
    end
end

function log_images(logger::TBLogger, name::AbstractString, imgArrays::AbstractArray, format::ImageFormat; step = nothing)
    for (n, imgArray) in enumerate(imgArrays)
        log_image(logger, name*"/$n", imgArray, format, step = step)
    end
end
"""
    log_image(logger::TBLogger, name::AbstractString, imgArray::AbstractArray, format::ImageFormat, step = nothing)

Log an image using image data and format
- imgArray: image data. A 1-D, 2-D or 3-D `Array` of pixel values. pixel values can be Real [0, 1] or Integer[0, 255]
- format: format of the image. It can be one of
  {L, LN, NL,
   CL, LC, NCL, NLC, CLN, LCN,
   HW, WH, HWC, WHC, CHW, CWH,
   HWN, WHN, NHW, NWH,
   HWCN, WHCN, CHWN, CWHN, NHWC, NWHC, NCHW, NCWH}
  - L: Length
  - C: Channel/Color
  - H: Height
  - W: Width
  - N: Observation
"""
function log_image(logger::TBLogger, name::AbstractString, img::AbstractArray{<:Colorant}; step = nothing)
    dimensions = ndims(img)
    #if it has 3 dimensions, eg MRI image, log each channel vth same name
    #so that one can slide through channels
    if dimensions == 3
        #3rd d is channel dim as observed in testimages
        channels = size(img, 3)
        for c in 1:channels-1
            log_image(logger, name, img[:, :, c], step = step)
        end
        img = img[:, :, channels]
    end
    summ = SummaryCollection(image_summary(name, img))
    write_event(logger.file, make_event(logger, summ, step=step))
end

function log_image(logger::TBLogger, name::AbstractString, imgArray::AbstractArray, format::ImageFormat; step=nothing)
    imgArray = channelview(imgArray)
    dims = ndims(imgArray)
    @assert dims == expected_ndims(format)
    obsdim = obs_dim(format)
    if iszero(obsdim)
        summ = SummaryCollection()
        push!(summ.value, image_summary(name, imgArray, format))
        write_event(logger.file, make_event(logger, summ, step=step))
    else
        format = strip_obs[format]
        index = collect("[:"* ",:"^(dims-1) *"]")
        index[2*obsdim] = 'g'
        index = join(index)
        global gimgArray = imgArray
        nth_img = Meta.parse("gimgArray$(index)")
        for n in 1:size(imgArray, obsdim)
            global g = n
            log_image(logger, name*"/$n", eval(nth_img), format, step = step)
        end
    end
end

function image_summary(name::AbstractString, img::AbstractArray{<:Colorant})
    #image is of type AbstractArray{<:Colorant}
    getcolorspaceid(::Type{<:Gray}) = 1
    getcolorspaceid(::Type{<:GrayA}) = 2
    getcolorspaceid(::Type{<:RGB}) = 3
    getcolorspaceid(::Type{<:RGBA}) = 4
    getcolorspaceid(::Type{<:BGRA}) = 6
    colorspace = getcolorspaceid(eltype(img))
    dimensions = ndims(img)
    if dimensions == 1
        #Grayscale/color array
        height = 1
        width = size(img, 1)
    elseif dimensions == 2
        #Grayscale/color matrix
        height, width = size(img)
    else
        throw("Too many dimensions")
    end
    #save image in a buffer
    io = IOBuffer()
    save(Stream(format"PNG", io), img)
    #read from buffer to obtain encoded string of the image
    eis = io.data
    imgsumm = Summary_Image(height = height, width = width, colorspace = colorspace, encoded_image_string = eis)
    Summary_Value(tag = name, image = imgsumm)
end

function image_summary(name::AbstractString, imgArray::AbstractArray, format::ImageFormat)
    #logger and step are only relevant when using explicit function `log_image` and format contains N
    #unpack RGB, RGBA value to channels using channelview
    imgArray = channelview(imgArray)
    #if data contains integer numbers, scale them to 0-255
    if isa(first(imgArray), Integer)
        imgArray = (imgArray./255)
    end
    #convert all values to `Float64` for uniformity
    imgArray = Float64.(imgArray)
    #scale all values to 0-1
    imgArray = (imgArray./(max(maximum(imgArray), 1)))
    #convert any format to CHW
    imgArray =
    format == L   ? reshape(imgArray, (1, 1, size(imgArray, 1))) :
    format == CL  ? reshape(imgArray, (size(imgArray, 1), 1, size(imgArray, 2))) :
    format == LC  ? reshape(transpose(imgArray), (size(imgArray, 2), 1, size(imgArray, 1))) :
    format == HW  ? reshape(imgArray, (1, size(imgArray, 1), size(imgArray, 2))) :
    format == WH  ? reshape(transpose(imgArray), (1, size(imgArray, 2), size(imgArray, 1))) :
    format == HWC ? permutedims(imgArray, (3, 1, 2)) :
    format == WHC ? permutedims(imgArray, (3, 2, 1)) :
    format == CHW ? imgArray :
    format == CWH ? permutedims(imgArray, (1, 3, 2)) :
    #== else ==#    throw("Invalid format")
    channels, height, width = size(imgArray)
    color =
    channels == 1 ? Gray :
    channels == 2 ? GrayA :
    channels == 3 ? RGB :
    channels == 4 ? RGBA :
    #== else ==#    throw("Too many channels")
    #if it is a single channel Array, convert it to HW
    if color == Gray
        imgArray = imgArray[1, :, :]
    end
    #convert Array to PNG pass it to image_summary
    img = colorview(color, imgArray)
    image_summary(name, img)
end


text:
struct blosc2_metalayer_unsafe
  name::Cstring #!< The metalayer identifier for Blosc client (e.g. Caterva).
  content::Ptr{UInt8} #!< The serialized (msgpack preferably) content of the metalayer.
  content_len::Int32 #!< The length in bytes of the content.
end

"""
    blosc2_storage

This struct is meant for holding storage parameters for a for a blosc2 container,
allowing to specify, for example, how to interpret the contents included in the schunk.
"""
struct blosc2_storage

    contiguous ::Bool #!< Whether the chunks are contiguous or sparse.
    urlpath ::Cstring #!< The path for persistent storage. If NULL, that means in-memory.
    cparams ::Ptr{blosc2_cparams} #!< The compression params when creating a schunk. If NULL, sensible defaults are used depending on the context.
    dparams ::Ptr{blosc2_dparams} #!< The decompression params when creating a schunk. If NULL, sensible defaults are used depending on the context.
    io::Ptr{blosc2_io} #!< Input/output backend.

    blosc2_storage(;
        contiguous = false,
        urlpath = C_NULL,
        cparams = C_NULL,
        dparams = C_NULL,
        io = C_NULL
    ) = new(contiguous, urlpath, cparams, dparams, io)
end

"""
    blosc2_schunk

This struct is the standard container for Blosc 2 compressed data.

This is essentially a container for Blosc 1 chunks of compressed data, and it allows to overcome the 32-bit limitation in Blosc 1.
Optionally, a blosc2_frame can be attached so as to store the compressed chunks contiguously.
"""
struct blosc2_schunk
    version ::UInt8
    compcode ::UInt8 #!< The default compressor. Each chunk can override this.
    compcode_meta ::UInt8 #!< The default compressor metadata. Each chunk can override this.
    clevel ::UInt8 #!< The compression level and other compress params.
    typesize ::Int32 #!< The type size.
    blocksize ::Int32 #!< The requested size of the compressed blocks (0; meaning automatic).
    chunksize ::Int32 #!< Size of each chunk. 0 if not a fixed chunksize.
    filters ::FiltersTuple #!< The (sequence of) filters.  8-bit per filter.
    filters_meta ::FiltersTuple #!< Metadata for filters. 8-bit per meta-slot.
    nchunks::Int32 #!< Number of chunks in super-chunk.
    nbytes ::Int64 #!< The data size + metadata size + header size (uncompressed).
    cbytes ::Int64 #!< The data size + metadata size + header size (compressed).
    data ::Ptr{Ptr{UInt8}} #! Pointer to chunk data pointers buffer.
    data_len ::Csize_t #!< Length of the chunk data pointers buffer.
    storage ::Ptr{blosc2_storage} #!< Pointer to storage info.
    frame ::Ptr{Cvoid} #!< Pointer to frame used as store for chunks.
    cctx ::Ptr{blosc2_context} #!< Context for compression
    dctx ::Ptr{blosc2_context} #!< Context for decompression.
    metalayers ::NTuple{BLOSC2_MAX_METALAYERS, blosc2_metalayer_unsafe} #!< The array of metalayers.
    nmetalayers ::Int16 #!< The number of metalayers in the super-chunk
    vlmetalayers ::NTuple{BLOSC2_MAX_VLMETALAYERS, blosc2_metalayer_unsafe} #<! The array of variable-length metalayers.
    nvlmetalayers ::Int16 #!< The number of variable-length metalayers.
    udbtune ::Ptr{Cvoid}
end

"""
    blosc2_schunk_new(storage::blosc2_storage)::Ptr{blosc2_storage}

Create a new super-chunk.

# Remark

In case that `storage.urlpath` is not NULL, the data is stored on-disk. If the data file(s) exist, they are overwritten.

# Parameters

-  `storage` – The storage properties.

"""
function blosc2_schunk_new(storage::blosc2_storage = blosc2_storage())
    rstorage = Ref(storage)
    @ccall lib.blosc2_schunk_new(rstorage::Ptr{blosc2_storage})::Ptr{blosc2_schunk}
end

"""
    blosc2_schunk_free(schunk::blosc2_schunk)

Release resources from a super-chunk.

# Remark
All the memory resources attached to the super-frame are freed. If the super-chunk is on-disk, the data continues there for a later re-opening.

# Returns
0 if success
"""
function blosc2_schunk_free(schunk::Ptr{blosc2_schunk})
    @ccall lib.blosc2_schunk_free(schunk::Ptr{blosc2_schunk})::Cint
end

"""
    blosc2_schunk_append_buffer(schunk::Ptr{blosc2_schunk}, src, nbytes)

Append a src data buffer to a super-chunk.

# Parameters

- `schunk` – The super-chunk where data will be appended.

- `src` – The buffer of data to compress.

- `nbytes` – The size of the src buffer.

# Returns

The number of chunks in super-chunk. If some problem is detected, this number will be negative.
"""
function blosc2_schunk_append_buffer(schunk::Ptr{blosc2_schunk}, src, nbytes)
    @ccall lib.blosc2_schunk_append_buffer(schunk::Ptr{blosc2_schunk}, src::Ptr{Cvoid}, nbytes::Int32)::Cint
end

"""
    blosc2_schunk_decompress_chunk(schunk::Ptr{blosc2_schunk}, nchunk, dest, nbytes)

Decompress and return the `nchunk` chunk of a super-chunk.

If the chunk is uncompressed successfully, it is put in the `dest`

# Arguments

- `schunk` – The super-chunk from where the chunk will be decompressed.

- `nchunk` – The chunk to be decompressed (0 indexed).

- `dest` – The buffer where the decompressed data will be put.

- `nbytes` – The size of the dest.

# Returns

The size of the decompressed chunk or 0 if it is non-initialized. If some problem is detected, a negative code is returned instead.
"""
function blosc2_schunk_decompress_chunk(schunk::Ptr{blosc2_schunk}, nchunk, dest, nbytes)
    @ccall lib.blosc2_schunk_decompress_chunk(schunk::Ptr{blosc2_schunk}, nchunk::Cint, dest::Ptr{Cvoid}, nbytes::Int32)::Cint
end

"""
    blosc2_schunk_get_chunk(schunk::Ptr{blosc2_schunk}, nchunk)::Vector{UInt8}

Return a compressed chunk that is part of a super-chunk in the chunk parameter.

# Arguments

- `schunk` – The super-chunk from where to extract a chunk.

- `nchunk` – The chunk to be extracted (0 indexed).

C signature `int blosc2_schunk_get_chunk(blosc2_schunk *schunk, int nchunk, uint8_t **chunk, bool *needs_free)`
"""
function blosc2_schunk_get_chunk(schunk::Ptr{blosc2_schunk}, nchunk)::Vector{UInt8}
   chunk_ptr = Ref{Ptr{UInt8}}()
   needs_free = Ref{Bool}()
   r = @ccall lib.blosc2_schunk_get_chunk(schunk::Ptr{blosc2_schunk}, nchunk::Cint, chunk_ptr::Ptr{Ptr{UInt8}}, needs_free::Ptr{Bool})::Cint
   r < 0 && error("blosc2_schunk_get_chunk error")
   r == 0 && return Vector{UInt8}(undef, 0)
   return unsafe_wrap(Array, chunk_ptr[], r, own = needs_free[])
end

"""
    blosc2_schunk_append_chunk(schunk::Ptr{blosc2_schunk}, chunk, copy)

Append an existing chunk o a super-chunk.

# Arguments

- `schunk` – The super-chunk where the chunk will be appended.

- `chunk` – The chunk to append. An internal copy is made, so chunk can be reused or freed if desired.

- `copy` – Whether the chunk should be copied internally or can be used as-is.

#Returns

The number of chunks in super-chunk. If some problem is detected, this number will be negative.
"""
function blosc2_schunk_append_chunk(schunk::Ptr{blosc2_schunk}, chunk, copy)
    @ccall lib.blosc2_schunk_append_chunk(schunk::Ptr{blosc2_schunk}, chunk::Ptr{UInt8}, copy::Bool)::Cint
end

"""
    blosc2_schunk_insert_chunk(schunk::Ptr{blosc2_schunk}, nchunk, chunk, copy)

Insert a chunk at a specific position in a super-chunk.

# Arguments

- `schunk` – The super-chunk where the chunk will be appended.

- `nchunk` – The position where the chunk will be inserted.(zero indexed)

- `chunk` – The chunk to insert. If an internal copy is made, the chunk can be reused or freed if desired.

- `copy` – Whether the chunk should be copied internally or can be used as-is.

# Returns

The number of chunks in super-chunk. If some problem is detected, this number will be negative.

"""
function blosc2_schunk_insert_chunk(schunk::Ptr{blosc2_schunk}, nchunk, chunk, copy)
    @ccall lib.blosc2_schunk_insert_chunk(schunk::Ptr{blosc2_schunk}, nchunk::Cint, chunk::Ptr{UInt8}, copy::Bool)::Cint
end

"""
    blosc2_schunk_update_chunk(schunk::ptr{blosc2_schunk}, nchunk, chunk, copy)

Update a chunk at a specific position in a super-chunk.

# Arguments

- `schunk` – The super-chunk where the chunk will be updated.

- `nchunk` – The position where the chunk will be updated.(zero indexed)

- `chunk` – The new chunk. If an internal copy is made, the chunk can be reused or freed if desired.

- `copy` – Whether the chunk should be copied internally or can be used as-is.

# Returns

The number of chunks in super-chunk. If some problem is detected, this number will be negative.
"""
function blosc2_schunk_update_chunk(schunk::Ptr{blosc2_schunk}, nchunk, chunk::Vector{UInt8}, copy)
    @ccall lib.blosc2_schunk_update_chunk(schunk::Ptr{blosc2_schunk}, nchunk::Cint, chunk::Ptr{UInt8}, copy::Bool)::Cint
end

"""
    blosc2_schunk_delete_chunk(schunk::ptr{blosc2_schunk}, nchunk)

Delete a chunk at a specific position in a super-chunk.

Parameters
- `schunk` – The super-chunk where the chunk will be deleted.

- `nchunk` – The position where the chunk will be deleted. (zero indexed)

# Returns

The number of chunks in super-chunk. If some problem is detected, this number will be negative.
"""
function blosc2_schunk_delete_chunk(schunk::Ptr{blosc2_schunk}, nchunk)
    @ccall lib.blosc2_schunk_delete_chunk(schunk::Ptr{blosc2_schunk}, nchunk::Cint)::Cint
end

"""
    blosc2_schunk_reorder_offsets(schunk::Ptr{blosc2_schunk}, offsets_order::Vector{Cint})

Reorder the chunk offsets of an existing super-chunk.

# Arguments

- `schunk` – The super-chunk whose chunk offsets are to be reordered.

- `offsets_order` – The new order of the chunk offsets (zero ordered)

# Returns

0 if suceeds. Else a negative code is returned.
"""
function blosc2_schunk_reorder_offsets(schunk::Ptr{blosc2_schunk}, offsets_order::Vector{Cint})
    @ccall lib.blosc2_schunk_reorder_offsets(schunk::Ptr{blosc2_schunk}, offsets_order::Ptr{Cint})::Cint
end

"""
    blosc2_schunk_get_cparams(schunk::Ptr{blosc2_schunk})

Return the cparams associated to a super-chunk.
"""
function blosc2_schunk_get_cparams(schunk::Ptr{blosc2_schunk})
   params_ptr = Ref{Ptr{blosc2_cparams}}()
   r = @ccall lib.blosc2_schunk_get_cparams(schunk::Ptr{blosc2_schunk}, params_ptr::Ptr{Ptr{blosc2_cparams}})::Cint
   r < 0 && error("blosc2_schunk_get_cparams error")
   return unsafe_wrap(Array, params_ptr[], 1, own = true)[1]
end

"""
    blosc2_schunk_get_dparams(schunk::Ptr{blosc2_schunk})

Return the dparams associated to a super-chunk.
"""
function blosc2_schunk_get_dparams(schunk::Ptr{blosc2_schunk})
   params_ptr = Ref{Ptr{blosc2_dparams}}()
   r = @ccall lib.blosc2_schunk_get_dparams(schunk::Ptr{blosc2_schunk}, params_ptr::Ptr{Ptr{blosc2_dparams}})::Cint
   r < 0 && error("blosc2_schunk_get_dparams error")
   return unsafe_wrap(Array, params_ptr[], 1, own = true)[1]
end


text:
@doc raw"""
`TableauEPRK`: Tableau of an Explicit Partitioned Runge-Kutta method
```math
\begin{align*}
V_{n,i} &= \hphantom{-} \dfrac{\partial H}{\partial p} (Q_{n,i}, P_{n,i}) , &
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , &
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} , \\
F_{n,i} &= - \dfrac{\partial H}{\partial q} (Q_{n,i}, P_{n,i}) , &
P_{n,i} &= p_{n} + h  \sum \limits_{i=1}^{s} \bar{a}_{ij} \, F_{n,j} , &
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{b}_{i} \, F_{n,i} ,
\end{align*}
```
usually satisfying the symplecticity conditions
```math
\begin{align*}
b_{i} \bar{a}_{ij} + b_{j} a_{ji} &= b_{i} b_{j} , &
\bar{b}_i &= b_i .
\end{align*}
```
"""
struct TableauEPRK{T} <: AbstractTableauPRK{T}
    @HeaderTableau

    q::CoefficientsRK{T}
    p::CoefficientsRK{T}

    function TableauEPRK{T}(name, o, q, p) where {T}
        @assert q.s==p.s
        # TODO check that both tableaus are lower triangular and that only one element
        #      a_q[i,i] or a_p[i,i] is non-zero for all i.
        new(name, o, q.s, q, p)
    end
end

function TableauEPRK(name::Symbol, order::Int, q::CoefficientsRK{T}, p::CoefficientsRK{T}) where {T}
    TableauEPRK{T}(name, order, q, p)
end

function TableauEPRK(name::Symbol, order::Int, q::CoefficientsRK{T}) where {T}
    TableauEPRK{T}(name, order, q, q)
end

# TODO function readAbstractTableauEPRKFromFile(dir::AbstractString, name::AbstractString)

"Parameters for right-hand side function of explicit partitioned Runge-Kutta methods."
struct ParametersEPRK{DT, TT, ET <: PODE{DT,TT}, D, S} <: Parameters{DT,TT}
    equ::ET
    tab::TableauEPRK{TT}
    Δt::TT
end

function ParametersEPRK(equ::ET, tab::TableauEPRK{TT}, Δt::TT) where {DT, TT, ET <: PODE{DT,TT}}
    ParametersEPRK{DT, TT, ET, equ.d, tab.s}(equ, tab, Δt)
end

"Explicit Runge-Kutta integrator cache."
struct IntegratorCacheEPRK{DT,D,S} <: PODEIntegratorCache{DT,D}
    Q::OffsetArray{Array{DT,1},1,Array{Array{DT,1},1}}
    P::OffsetArray{Array{DT,1},1,Array{Array{DT,1},1}}

    Y::Vector{Vector{DT}}
    Z::Vector{Vector{DT}}
    V::Vector{Vector{DT}}
    F::Vector{Vector{DT}}

    function IntegratorCacheEPRK{DT,D,S}() where {DT,D,S}
        Q = create_internal_stage_vector_with_zero(DT, D, S)
        P = create_internal_stage_vector_with_zero(DT, D, S)

        Y = create_internal_stage_vector(DT, D, S)
        Z = create_internal_stage_vector(DT, D, S)
        V = create_internal_stage_vector(DT, D, S)
        F = create_internal_stage_vector(DT, D, S)

        new(Q, P, Y, Z, V, F)
    end
end

"Explicit partitioned Runge-Kutta integrator."
struct IntegratorEPRK{DT, TT, PT <: ParametersEPRK{DT,TT}, D, S} <: IntegratorPRK{DT,TT}
    params::PT
    cache::IntegratorCacheEPRK{DT,D,S}

    function IntegratorEPRK(equation::PODE{DT,TT}, tableau::TableauEPRK{TT}, Δt::TT) where {DT,TT}
        D = equation.d
        S = tableau.s

        # create params
        params = ParametersEPRK(equation, tableau, Δt)

        # create cache
        cache = IntegratorCacheEPRK{DT,D,S}()

        # create integrators
        new{DT, TT, typeof(params), D, S}(params, cache)
    end
end

"Compute Q stages of explicit partitioned Runge-Kutta methods."
function computeStageQ!(int::IntegratorEPRK{DT,TT}, cache::IntegratorCacheEPRK{DT}, i::Int, jmax::Int, t::TT) where {DT,TT}
    fill!(cache.Y[i], zero(DT))
    for j in 1:jmax
        for k in eachdim(int)
            cache.Y[i][k] += tableau(int).q.a[i,j] * cache.V[j][k]
        end
    end
    for k in eachdim(int)
        cache.Q[i][k] = cache.Q[0][k] + timestep(int) * cache.Y[i][k]
    end
    equation(int).f(t, cache.Q[i], cache.P[jmax], cache.F[i])
end

"Compute P stages of explicit partitioned Runge-Kutta methods."
function computeStageP!(int::IntegratorEPRK{DT,TT}, cache::IntegratorCacheEPRK{DT}, i::Int, jmax::Int, t::TT) where {DT,TT}
    fill!(cache.Z[i], zero(DT))
    for j in 1:jmax
        for k in eachdim(int)
            cache.Z[i][k] += tableau(int).p.a[i,j] * cache.F[j][k]
        end
    end
    for k in eachdim(int)
        cache.P[i][k] = cache.P[0][k] + timestep(int) * cache.Z[i][k]
    end
    equation(int).v(t, cache.Q[jmax], cache.P[i], cache.V[i])
end

"Integrate partitioned ODE with explicit partitioned Runge-Kutta integrator."
function integrate_step!(int::IntegratorEPRK{DT,TT}, sol::AtomicSolutionPODE{DT,TT}) where {DT,TT}
    local tqᵢ::TT
    local tpᵢ::TT

    # store previous solution
    int.cache.Q[0] .= sol.q
    int.cache.P[0] .= sol.p

    # reset cache
    reset!(sol, timestep(int))

    # compute internal stages
    for i in eachstage(int)
        tqᵢ = sol.t̅ + timestep(int) * tableau(int).q.c[i]
        tpᵢ = sol.t̅ + timestep(int) * tableau(int).p.c[i]

        if tableau(int).q.a[i,i] ≠ zero(TT) && tableau(int).p.a[i,i] ≠ zero(TT)
            error("This is an implicit method!")
        elseif tableau(int).q.a[i,i] ≠ zero(TT)
            computeStageP!(int, int.cache, i, i-1, tpᵢ)
            computeStageQ!(int, int.cache, i, i, tqᵢ)
        elseif tableau(int).p.a[i,i] ≠ zero(TT)
            computeStageQ!(int, int.cache, i, i-1, tqᵢ)
            computeStageP!(int, int.cache, i, i, tpᵢ)
        else
            computeStageQ!(int, int.cache, i, i-1, tqᵢ)
            computeStageP!(int, int.cache, i, i-1, tpᵢ)
        end
    end

    # compute final update
    update_solution!(sol.q, sol.q̃, int.cache.V, tableau(int).q.b, tableau(int).q.b̂, timestep(int))
    update_solution!(sol.p, sol.p̃, int.cache.F, tableau(int).p.b, tableau(int).p.b̂, timestep(int))
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1PodAffinity <: SwaggerModel
    preferredDuringSchedulingIgnoredDuringExecution::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1WeightedPodAffinityTerm} } # spec name: preferredDuringSchedulingIgnoredDuringExecution
    requiredDuringSchedulingIgnoredDuringExecution::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1PodAffinityTerm} } # spec name: requiredDuringSchedulingIgnoredDuringExecution

    function IoK8sApiCoreV1PodAffinity(;preferredDuringSchedulingIgnoredDuringExecution=nothing, requiredDuringSchedulingIgnoredDuringExecution=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PodAffinity, Symbol("preferredDuringSchedulingIgnoredDuringExecution"), preferredDuringSchedulingIgnoredDuringExecution)
        setfield!(o, Symbol("preferredDuringSchedulingIgnoredDuringExecution"), preferredDuringSchedulingIgnoredDuringExecution)
        validate_property(IoK8sApiCoreV1PodAffinity, Symbol("requiredDuringSchedulingIgnoredDuringExecution"), requiredDuringSchedulingIgnoredDuringExecution)
        setfield!(o, Symbol("requiredDuringSchedulingIgnoredDuringExecution"), requiredDuringSchedulingIgnoredDuringExecution)
        o
    end
end # type IoK8sApiCoreV1PodAffinity

const _property_map_IoK8sApiCoreV1PodAffinity = Dict{Symbol,Symbol}(Symbol("preferredDuringSchedulingIgnoredDuringExecution")=>Symbol("preferredDuringSchedulingIgnoredDuringExecution"), Symbol("requiredDuringSchedulingIgnoredDuringExecution")=>Symbol("requiredDuringSchedulingIgnoredDuringExecution"))
const _property_types_IoK8sApiCoreV1PodAffinity = Dict{Symbol,String}(Symbol("preferredDuringSchedulingIgnoredDuringExecution")=>"Vector{IoK8sApiCoreV1WeightedPodAffinityTerm}", Symbol("requiredDuringSchedulingIgnoredDuringExecution")=>"Vector{IoK8sApiCoreV1PodAffinityTerm}")
Base.propertynames(::Type{ IoK8sApiCoreV1PodAffinity }) = collect(keys(_property_map_IoK8sApiCoreV1PodAffinity))
Swagger.property_type(::Type{ IoK8sApiCoreV1PodAffinity }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1PodAffinity[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PodAffinity }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PodAffinity[property_name]

function check_required(o::IoK8sApiCoreV1PodAffinity)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1PodAffinity }, name::Symbol, val)
end


text:
"""
Compute the Sweldens quadrature rule for integrals involving a refinable function.

The method is described in:
W. Sweldens and R. Piessens, "Quadrature formulae and asymptotic error expansions
for wavelet approximations of smooth functions", SIAM J. Numer. Anal. 31(4),
pp. 1240-1264, 1994.

We return a quadrature rule that with points spaced 1/M apart, symmetric with
respect to the origin, and including the origin itself.
"""
function refinable_quadrature(coefficients, M = length(coefficients)-1)
    @assert sum(coefficients) ≈ sqrt(2)
    L = length(coefficients)
    # Perform calculations in BigFloat for greater accuracy
    if isodd(L) || (iseven(L) && iseven(M))
        t = range(-BigFloat(1.),stop = BigFloat(1.), length = M+1)
    else
        # We have to shift the grid in order to include the center point
        t = range(-BigFloat(1.),stop = BigFloat(1.), length = 2M+1)[2:2:end-1]
    end
    A = chebyshev_vandermonde_matrix(t)
    m = chebyshev_moments(BigFloat.(coefficients), length(t)-1)
    w = A\m
    (L-1) * t/2, map(eltype(coefficients), w)[:]
end

function chebyshev_vandermonde_matrix(t::AbstractVector{T}) where {T}
    M = length(t)-1
    A = zeros(T, M+1, M+1)
    for i in 0:M
        A[i+1,:] = chebyshev_eval.(i, t)
    end
    A
end

"Evaluate the chebyshev polynomial of degree p in x."
chebyshev_eval(p, x) = real(cos(p*acos(x)))

"""
Calculate the modified moments of the scaling function defined by coefficients.
"""
function chebyshev_moments(coefficients::AbstractVector, M::Int)
    T = eltype(coefficients)
    m = zeros(T,M+1)
    m[1] = 1
    L = length(coefficients)

    w = zeros(T,M+1,M+1,L)
    for k in 1:L
        w[:,:,k] = wcoefficients(T, L, k, M)
    end

    for p in 1:M
        d = zero(T)
        for i in 0:p-1
            e = zero(T)
            for k in 1:L
                # % CAVE: our coefficients differ by a factor 1/sqrt(2) from Sweldens!
                e = e + 1/sqrt(T(2))*coefficients[k]*w[p+1,i+1,k]
            end
            d = d + e*m[i+1]
        end
        m[p+1] = d/(2^p-1)
    end

    m
end

# Calculate the `w' coefficients in the algorithm by Sweldens and Piessens.
# This is based on Appendix A in the above-mentioned reference.
function wcoefficients(::Type{T}, L, k, M::Int) where {T}
    lambda = T(2)*(k-1)/(T(L)-1)-1

    w = zeros(T,M+1,M+1)

    w[1,1] = 1
    w[2,1] = lambda
    w[2,2] = 1
    if M > 1
        w[3,1] = 2*lambda^2-3
        w[3,2] = 4*lambda
        w[3,3] = 1
    end

    # p+1 is the degree of the next polynomial
    for p in 2:M-1
        w[p+2,1] = w[p+1,2] + 2*lambda*w[p+1,1] - 4*w[p,1]
        w[p+2,2] = 2*w[p+1,1]+w[p+1,3]+2*lambda*w[p+1,2]-4*w[p,2]
        for i in 2:p-1
            w[p+2,i+1] = w[p+1,i]+w[p+1,i+2]+2*lambda*w[p+1,i+1]-4*w[p,i+1]
        end
        w[p+2,p+1] = w[p+1,p] + 2*lambda*w[p+1,p+1]
        w[p+2,p+2] = w[p+1,p+1]
    end

    w
end

bspline_refinable_coefficients(degree::Int, ::Type{T} = Float64) where {T} =
    [binomial(degree+1,k) for k in 0:degree+1] / 2^(degree+one(T)/2)

################################
# Definition of QBF quadrature
################################

export QuadQBF
"""
QuadQBF uses specialized quadrature routines that incorporate the basis function of the discretization
into their weight function. This means only evaluation of the Green's function are required.
The evaluations are equispaced and may be shared for neighbouring elements in the discretization matrix.
"""
struct QuadQBF{T} <: DomainIntegrals.IntervalRule{T}
    coef    ::  Array{T,1}  #   Coefficients of the two-scale relation
    a       ::  T
    b       ::  T
    x       ::  Array{T,1}  #   1d quadrature points
    w       ::  Array{T,1}  #   1d quadrature weights
end

QuadQBF(args...; options...) = QuadQBF{Float64}(args...; options...)
QuadQBF(coef::AbstractArray{T}, args...; options...) where {T} = QuadQBF{T}(coef, args...; options...)
QuadQBF{T}(n::Int = 1; oversamplingfactor = 2) where {T} =
    QuadQBF{T}(n, oversamplingfactor*(n+1))
QuadQBF{T}(coef::AbstractVector; oversamplingfactor = 2) where {T} =
    QuadQBF{T}(coef, oversamplingfactor*length(coef)-1)

QuadQBF{T}(n::Int, M::Int) where {T} = QuadQBF{T}(bspline_refinable_coefficients(n, T), M)

function QuadQBF{T}(coef::AbstractVector, M::Int) where {T}
    x, w = refinable_quadrature(coef, M)
    QuadQBF{T}(coef, -(length(coef)-one(T))/2, (length(coef)-one(T))/2, x, w)
end

leftpoint(q::QuadQBF) = q.a
rightpoint(q::QuadQBF) = q.b

points(q::QuadQBF) = q.x
weights(q::QuadQBF) = q.w

DomainIntegrals.domain(q::QuadQBF) = q.a..q.b

splineorder(q::QuadQBF) = length(q.coef)-2

# Map a point x from the interval [a,b] linearly to the interval [c,d]
mapx(x, a, b, c, d) = c + (x-a)/(b-a)*(d-c)

function qbf_quadrature(f, dict, idx, measure, domain, qbf_a, qbf_b, qbf_x, qbf_w)
    a, b = extrema(domain)
    T = typeof(f(a))
    z = zero(T)
    for i = 1:length(qbf_x)
        t = mapx(qbf_x[i], qbf_a, qbf_b, a, b)
        z += qbf_w[i] * f(t) * unsafe_weightfun(measure, t)
    end
    z * sqrt((b-a) / (qbf_b-qbf_a))
end

function qbf_quadrature2(f, dict1, idx1, measure1, domain1, dict2, idx2, measure2, domain2, qbf_a, qbf_b, qbf_x, qbf_w)
    a1, b1 = extrema(domain1)
    a2, b2 = extrema(domain2)
    T = typeof(f(a1,a2))
    z = zero(T)
    for i = 1:length(qbf_x)
        t1 = mapx(qbf_x[i], qbf_a, qbf_b, a1, b1)
        for j = 1:length(qbf_x)
            t2 = mapx(qbf_x[j], qbf_a, qbf_b, a2, b2)
            z += qbf_w[i] * qbf_w[j] * f(t2,t1) * unsafe_weightfun(measure1, t1) * unsafe_weightfun(measure2, t2)
        end
    end
    z * sqrt((b1-a1)*(b2-a2)) / (qbf_b-qbf_a)
end


text:
module Core

using DistributionsAD, Bijectors
using Libtask, ForwardDiff, Random
using Distributions, LinearAlgebra
using ..Utilities, Reexport
using Tracker: Tracker
using ..Turing: Turing
using DynamicPPL: Model, AbstractSampler, Sampler, SampleFromPrior
using LinearAlgebra: copytri!
using Bijectors: PDMatDistribution
import Bijectors: link, invlink
using AdvancedVI
using StatsFuns: logsumexp, softmax
@reexport using DynamicPPL
using Requires

import AdvancedPS
import ZygoteRules

include("container.jl")
include("ad.jl")

function __init__()
    @require ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" begin
        include("compat/reversediff.jl")
        export ReverseDiffAD, getrdcache, setrdcache, emptyrdcache
    end
end

export  @model,
        @varname,
        generate_observe,
        translate_tilde!,
        get_vars,
        get_data,
        get_default_values,
        ParticleContainer,
        Particle,
        Trace,
        fork,
        forkr,
        current_trace,
        getweights,
        getweight,
        effectiveSampleSize,
        sweep!,
        ResampleWithESSThreshold,
        ADBackend,
        setadbackend,
        setadsafe,
        ForwardDiffAD,
        TrackerAD,
        ZygoteAD,
        value,
        gradient_logp,
        CHUNKSIZE,
        ADBACKEND,
        setchunksize,
        verifygrad,
        @logprob_str,
        @prob_str

end # module