text:
# NOTE: this is a quick and dirty implementation of the Stillinger-Weber model;
# in the future it will be good to do this more elegantly / more efficiently
# there is huge overhead in this code that can be
#   significantly improved in terms of performance, but for now we just
#   want a correct and readable code

# [Stillinger/Weber, PRB 1985]
# ---------------------------------------------------
# v2 = ϵ f2(r_ij / σ); f2(r) = A (B r^{-p} - r^{-q}) exp( (r-a)^{-1} )
# v3 = ϵ f3(ri/σ, rj/σ, rk/σ); f3 = h(rij, rij, Θjik) + ... + ...
# h(rij, rik, Θjik) = λ exp[ γ (rij-a)^{-1} + γ (rik-a)^{-1} ] * (cosΘjik+1/3)^2
#       >>>
# V2 = 0.5 * ϵ * A (B r^{-p} - r^{-q}) * exp( (r-a)^{-1} )
# V3 = √ϵ * λ exp[ γ (r-a)^{-1} ]
#
# Parameters from QUIP database:
# -------------------------------
# <per_pair_data atnum_i="14" atnum_j="14" AA="7.049556277" BB="0.6022245584"
         # p="4" q="0" a="1.80" sigma="2.0951" eps="2.1675" />
# <per_triplet_data atnum_c="14" atnum_j="14" atnum_k="14" lambda="21.0"
         # gamma="1.20" eps="2.1675" />

using ForwardDiff
export StillingerWeber

# using JuLIP.Potentials: evaluate_dd, @D, @DD
# import JuLIP.Potentials: evaluate, evaluate_d

"""
`bondangle(S1, S2) -> (dot(S1, S2) + 1.0/3.0)^2`

* not this assumes that `S1, S2` are normalised
* see `bondangle_d` for the derivative
"""
bondangle(S1, S2) = (dot(S1, S2) + 1.0/3.0)^2

"""
`b := bondangle(S1, S2)` then

`bondangle_d(S1, S2, r1, r2) -> b, db1, db2`

where `dbi` is the derivative of `b` w.r.t. `Ri` where `Si= Ri/ri`.
"""
function bondangle_d(S1, S2, r1, r2)
   d = dot(S1, S2)
   b1 = (1.0/r1) * S2 - (d/r1) * S1
   b2 = (1.0/r2) * S1 - (d/r2) * S2
   return (d+1./3.)^2, 2.0*(d+1./3.)*b1, 2.0*(d+1./3.)*b2
end

function _ad_bondangle_(R)
   R1, R2 = R[1:3], R[4:6]
   r1, r2 = norm(R1), norm(R2)
   return bondangle(R1/r1, R2/r2)
end

# TODO: need a faster implementation of bondangle_dd
function bondangle_dd(R1, R2)
   R = [R1; R2]
   hh = ForwardDiff.hessian(_ad_bondangle_, R)
   h = zeros(JMatF, 2,2)
   h[1,1] = JMatF(hh[1:3,1:3])
   h[1,2] = JMatF(hh[1:3, 4:6])
   h[2,1] = JMatF(hh[4:6,1:3])
   h[2,2] = JMatF(hh[4:6, 4:6])
   return h
end

@pot struct StillingerWeber{P1,P2} <: SitePotential
   V2::P1
   V3::P2
end

"""
Stillinger-Weber potential with parameters for Si.

Functional form and default parameters match the original SW potential
from [Stillinger/Weber, PRB 1985].

The `StillingerWeber` type can also by "abused" to generate arbitrary
bond-angle potentials of the form
   ∑_{i,j} V2(rij) + ∑_{i,j,k} V3(rij) V3(rik) (cos Θijk + 1/3)^2

Constructor admits the following key-word parameters:
`ϵ=2.1675, σ = 2.0951, A=7.049556277, B=0.6022245584,
                  p = 4, a = 1.8, λ=21.0, γ=1.20`

which enter the potential as follows:
```
V2(r) = 0.5 * ϵ * A * (B * (r/σ)^(-p) - 1.0) * exp(1.0 / (r/σ - a))
V3(r) = sqrt(ϵ * λ) * exp(γ / (r/σ - a))
```
"""
StillingerWeber

cutoff(calc::StillingerWeber) = max(cutoff(calc.V2), cutoff(calc.V3))

# TODO: brittle StillingerWeber, make λ = 42.0

function StillingerWeber(; brittle = false,
               ϵ=2.1675, σ = 2.0951, A=7.049556277, B=0.6022245584,
               p = 4, a = 1.8, λ = brittle ? 42.0 : 21.0, γ=1.20 )
   cutoff = a*σ-1e-2
   V2 = @analytic(r -> (ϵ*A) * (B*(r/σ)^(-p) - 1.0) * exp(1.0/(r/σ - a))) *
         HS(cutoff)
   V3 = @analytic(r -> sqrt(ϵ * λ) * exp( γ / (r/σ - a) )) * HS(cutoff)
   return StillingerWeber(V2, V3)
end

function evaluate(calc::StillingerWeber, r, R)
   # two-body contributions
   if length(r) == 0
      return 0.0
   end
   Es = 0.5 * sum( calc.V2(s) for s in r )
   # three-body contributions
   S = [ R1/r1 for (R1,r1) in zip(R, r) ]
   V3 = [ calc.V3(r1) for r1 in r ]
   for i1 = 1:(length(r)-1), i2 = (i1+1):length(r)
      Es += V3[i1] * V3[i2] * bondangle(S[i1], S[i2])
   end
   return Es
end

# function energy(calc::StillingerWeber, at::ASEAtoms)
#    nlist = neighbourlist(at, cutoff(calc))
#    # 2-body contribution
#    E = sum(calc.V2, nlist.r)
#    # 3-body contribution
#    V3 = [calc.V3(r)  for r in nlist.r]
#    n = 0
#    for idx = 1:length(at)
#       n += 1
#       a = n    # TODO: this should be built into the neighbourlist datastructure
#       while n < length(nlist) && nlist.i[n+1] == idx
#          n += 1
#       end
#       b = n
#       for i1 = a:(b-1), i2 = (i1+1):(b)
#          E += V3[i1] * V3[i2] * bondangle(nlist.R[i1]/nlist.r[i1], nlist.R[i2]/nlist.r[i2])
#       end
#    end
#
#    return E
# end

function evaluate_d(calc::StillingerWeber, r, R)
   # two-body terms
   if length(r) == 0
      return JVecF[]
   end
   dEs = [ 0.5 * grad(calc.V2, ri, Ri) for (ri, Ri) in zip(r, R) ]
   # three-body terms
   S = [ R1/r1 for (R1,r1) in zip(R, r) ]
   V3 = [calc.V3(s) for s in r]
   gV3 = [ grad(calc.V3, r1, R1) for (r1, R1) in zip(r, R) ]
   for i1 = 1:(length(r)-1), i2 = (i1+1):length(r)
      a, b1, b2 = bondangle_d(S[i1], S[i2], r[i1], r[i2])
      dEs[i1] += (V3[i1] * V3[i2]) * b1 + (V3[i2] * a) * gV3[i1]
      dEs[i2] += (V3[i1] * V3[i2]) * b2 + (V3[i1] * a) * gV3[i2]
   end
   return dEs
end

# function forces(calc::StillingerWeber, at::ASEAtoms)
#    nlist = neighbourlist(at, cutoff(calc))
#
#    # pair potential contribution to forces
#    dE = zerovecs(length(at))
#    for n = 1:length(nlist)
#       dE[nlist.i[n]] += 2 * grad(calc.V2, nlist.r[n], nlist.R[n])
#    end
#
#    # 3-body contribution
#    V3 = [calc.V3(r)  for r in nlist.r]
#    dV3 = [(@D calc.V3(r))/r  for r in nlist.r]
#    n = 0
#    for idx = 1:length(at)
#       n += 1
#       a = n
#       while n < length(nlist) && nlist.i[n+1] == idx
#          n += 1
#       end
#       b = n
#       for i1 = a:(b-1), i2 = (i1+1):(b)
#          α, b1, b2 = bondangle_d(nlist.R[i1]/nlist.r[i1],
#                            nlist.R[i2]/nlist.r[i2], nlist.r[i1], nlist.r[i2])
#          f1 = (V3[i1] * V3[i2]) * b1 + ((V3[i2] * α) * dV3[i1]) * nlist.R[i1]
#          f2 = (V3[i1] * V3[i2]) * b2 + ((V3[i1] * α) * dV3[i2]) * nlist.R[i2]
#          dE[nlist.j[i1]] -= f1
#          dE[nlist.i[i1]] += f1
#          dE[nlist.j[i2]] -= f2
#          dE[nlist.i[i2]] += f2
#       end
#    end
#    return dE
# end

function _ad_dV(V::StillingerWeber, R_dofs)
   R = vecs( reshape(R_dofs, 3, length(R_dofs) ÷ 3) )
   r = norm.(R)
   dV = evaluate_d(V, r, R)
   return mat(dV)[:]
end

function _ad_ddV(V::StillingerWeber, r, R)
   ddV = ForwardDiff.jacobian( Rdofs -> _ad_dV(V, Rdofs), mat(R)[:] )
   # convert into a block-format
   n = length(r)
   hV = zeros(JMatF, n, n)
   for i = 1:n, j = 1:n
      hV[i, j] = ddV[ ((i-1)*3)+(1:3), ((j-1)*3)+(1:3) ]
   end
   return hV
end

hess(V::StillingerWeber, r, R) = _ad_ddV(V, r, R)

# function hess(V::StillingerWeber, r, R)
#    n = length(r)
#    hV = zeros(JMatF, n, n)
#
#    # two-body contributions
#    for (i, (r_i, R_i)) in enumerate(zip(r, R))
#       hV[i,i] += hess(V.V2, r_i, R_i)
#    end
#
#    # three-body terms
#    S = [ R1/r1 for (R1,r1) in zip(R, r) ]
#    V3 = [ V.V3(r1) for r1 in r ]
#    dV3 = [ grad(V.V3, r1, R1) for (r1, R1) in zip(r, R) ]
#    hV3 = [ hess(V.V3, r1, R1) for (r1, R1) in zip(r, R) ]
#
#    for i1 = 1:(length(r)-1), i2 = (i1+1):length(r)
#       # Es += V3[i1] * V3[i2] * bondangle(S[i1], S[i2])
#       # precompute quantities
#       ψ, Dψ_i1, Dψ_i2 = bondangle_d(S[i1], S[i2], r[i1], r[i2])
#       Hψ = bondangle_dd(R[i1], R[i2])  # <<<< this should be SLOW (AD)
#       # assemble local hessian contributions
#       hV[i1,i1] +=
#          hV3[i1] * V3[i2] * ψ       +   dV3[i1] * V3[i2] * Dψ_i1' +
#          Dψ_i1 * V3[i2] * dV3[i1]'  +   V3[i1] * V3[i2] * Hψ[1,1]
#       hV[i2,i2] +=
#          V3[i2] * hV3[i2] * ψ       +   V3[i1] * dV3[i2] * Dψ_i2' +
#          Dψ_i2 * V3[i1] * dV3[i2]'  +   V3[i1] * V3[i2] * Hψ[2,2]
#       hV[i1,i2] +=
#          dV3[i1] * dV3[i2]' * ψ     +   V3[i1] * Dψ_i1 * dV3[i2]' +
#          dV3[i1] * V3[i2] * Dψ_i2'  +   V3[i1] * V3[i2] * Hψ[1,2]
#       hV[i2, i1] +=
#          dV3[i2] * dV3[i1]' * ψ     +   V3[i1] * Dψ_i2 * dV3[i1]' +
#          dV3[i2] * V3[i1] * Dψ_i1'  +   V3[i1] * V3[i2] * Hψ[2,1]
#    end
#    return hV
# end

function precon(V::StillingerWeber, r, R)
   n = length(r)
   pV = zeros(JMatF, n, n)

   # two-body contributions
   for (i, (r1, R1)) in enumerate(zip(r, R))
      pV[i,i] += precon(V.V2, r1, R1)
   end

   # three-body terms
   S = [ R1/r1 for (R1,r1) in zip(R, r) ]
   V3 = V.V3.(r)
   for i1 = 1:(n-1), i2 = (i1+1):n
      Θ = dot(S[i1], S[i2])
      dΘ1 = (1.0/r[i1]) * S[i2] - (Θ/r[i1]) * S[i1]
      dΘ2 = (1.0/r[i2]) * S[i1] - (Θ/r[i2]) * S[i2]
      # ψ = (Θ + 1/3)^2, ψ' = (Θ + 1/3), ψ'' = 2.0
      a = abs((V3[i1] * V3[i2] * 2.0))
      pV[i1,i2] += a * dΘ1 * dΘ2'
      pV[i1,i1] += a * dΘ1 * dΘ1'
      pV[i2, i2] += a * dΘ2 * dΘ2'
      pV[i2, i1] += a * dΘ2 * dΘ1'
   end

   return pV
end


text:
using DifferentialEquations
using DiffEqSensitivity
using LinearAlgebra
using Zygote
using Statistics
using Distributions
using Plots
using Random

# set random seed for reproducible results
Random.seed!(1)

## ----- setup true solution ------------------------------------

# define function for RHS of lorenz system
function lorenz!(du, u, p, t)
    x,y,z = u
    σ,ρ,β = p

    du[1] = σ*(y - x)
    du[2] = x*(ρ - z) - y
    du[3] = x*y - β*z
end

θ = [10.0, 28.0, (8/3)]  # true parameter values
dt = 0.01  # time step for simulation (this may be unnecessary)
tend = 10.0  # max time for total simulation
ts_sim = (0.0, tend)
u0True = [1.0, 1.0, 1.0]  # true initial condition

# integrate to obtain the true solution
prob_true = ODEProblem(lorenz!, u0True, ts_sim, θ, saveat=dt)
sol_true = solve(prob_true)

# visualize the data
# px = plot(sol_true, vars=(0,1), color=:black, alpha=0.75, linewidth=2, label="true system")
# xlabel!("t")
# ylabel!("x(t)")
# py = plot(sol_true, vars=(0,2), color=:black, alpha=0.75, linewidth=2, label="true system")
# xlabel!("t")
# ylabel!("y(t)")

# pz = plot(sol_true, vars=(0,3), color=:black, alpha=0.75, linewidth=2, label="true system")
# xlabel!("t")
# ylabel!("z(t)")

# train_plot = plot(px, py, pz, layout=(3,1))

## ----- setup observation data ------------------------------------

# set up observation uncertainties
σ_m = 0.15  # i.e. observation std
# set up observation error covariance matrix
R = σ_m^2 .* I(3)

# set up details for measurement points
dt_m = 0.2
tm_end = 2.0  # final observation time
ts_m = dt_m:dt_m:tm_end

# define observation operator
function h(u)
    w = u
    return w
end

# define observation operator jacobian. We should probably generate this symbolically
function Dh(u)
    n = size(u, 1)
    D = I(n)
    return D
end

# populate observation vector
w = hcat([h(sol_true(tᵢ))+rand(MvNormal(zeros(3), R)) for tᵢ ∈ ts_m]...)

# update the graphs with observation stuff
# plot!(px, ts_m, w[1, :], seriestype=:scatter, color="light green", label="observation" )
# plot!(py, ts_m, w[2, :], seriestype=:scatter, color="light green", label="observation" )
# plot!(pz, ts_m, w[3, :], seriestype=:scatter, color="light green", label="observation" )

# train_plot = plot(px, py, pz, layout=(3,1))

## ----- perform the assimilation ------------------------------------

# set up data assimilation
u0b = [2.0; 3.0; 4.0]  # guess at initial condition
σ_b = 0.1  # model error std
B = σ_b^2 .* I(3)  # initial model error covariance
Q = 0.0 * I(3)  # process noise (assumed zero)

# integrate to obtain the true solution
prob_no_analysis = ODEProblem(lorenz!, u0b, ts_sim, θ, saveat=dt)
sol_no_analysis = solve(prob_no_analysis)

# set up arrays for solution outputs... we should probably pre-allocate these
ua = zeros(3, size(0:dt:tend, 1))
ua[:,1] = u0b

# function for getting indices from time value
function timeindex(t, Δt)
    return Int(round(t/Δt)) + 1
end

# helper function to update the analysis result
function update!(sol, t_now, t_next, dt)
    idx_now = timeindex(t_now, dt)
    idx_next = timeindex(t_next, dt)

    ua[:, idx_now:idx_next] = hcat(sol(t_now:dt:t_next)...)
end
Zygote.@nograd update!  # tell zygote to ignore update! when computing derivatives

function model_forward!(uₖ)
    prob = ODEProblem(lorenz!, uₖ, (t_now, t_next), θ)
    sol = solve(prob, Tsit5(), reltol=1e-6,abstol=1e-6)

    update!(sol, t_now, t_next, dt)

    return sol[:, end]
end

# ------------- perform the analysis ---------------------------------------------
# times = hcat(collect(0:dt_m:tm_end)..., tend)
times = 0:dt_m:tend

t_now = times[1]
t_next = times[2]
#model_forward!(u0b)

res = Zygote.withjacobian(model_forward!, u0b)
res.val
res.grad[1]

for k ∈ 1:(length(times)-1)
    t_now = times[k]
    t_next = times[k+1]

    idx_now = timeindex(t_now, dt)
    u_now = ua[:, idx_now]

    # u_next, DM = model_forward!(u_now)
    u_next, DM = Zygote.withjacobian(model_forward!, u_now)
    DM = DM[1]

    # update error covariance matrix
    B = DM * B * DM' + Q

    if t_next ∈ ts_m
        println("doing the update!")
        n = size(ua, 1)
        idx_next = timeindex(t_next, dt)
        idx_h = timeindex(t_next, dt_m) - 1 # -1 since this doesn't include 0

        # compute the Kalman gain matrix
        DH = Dh(ua[:,idx_next])
        println("\t Dₕ=$(DH)")
        denom = DH*B*DH' + R
        println("\t B:= $(B)")
        println("\t denom:= $(denom)")
        K = B*DH'*inv(denom)

        println("\tK= $(K)")
        # perform the analysis

        println("\tua = $(ua[:, idx_next])")
        ua[:,idx_next] .= ua[:, idx_next] + K*(w[:, idx_h] - h(ua[:, idx_next]))

        println("\tua = $(ua[:, idx_next])")
        B = (I(n) - K*DH)*B
    end
end

var_names = ["x(t)", "y(t)", "z(t)"]
plots = []
for i ∈ 1:3
    p = plot(sol_true, vars=(0,i), color=:black, alpha=0.75, linewidth=2, label="true system")
    plot!(p, ts_m, w[i, :], seriestype=:scatter, color="light green", label="observation" )
    plot!(p, sol_no_analysis, vars=(0, i), linestyle=:dashdot, color="purple", label="forecast")
    plot!(p, 0:dt:tend, ua[i, :], color="royal blue", alpha=1, linewidth = 1, label="analysis")
    xlabel!(p, "t")
    ylabel!(p, var_names[i])
    push!(plots, p)
end

final_plot = plot(plots..., layout=(3,1))

savefig("EKF_using_zygote.svg")
savefig("EKF_using_zygote.png")
display(final_plot)

ext: jl


text:
module LaTeXDatax
using  Unitful, UnitfulLatexify, Latexify

export @datax

"""
```julia
@datax
```

Print the arguments to a data file to be read by pgfkeys. Best used with the
[`datax` LaTeX package](https://ctan.org/pkg/datax). Variables can be supplied
either by name or as assignments, and meta-arguments are supplied using the
`:=` operator.

The meta-arguments include all the keyword arguments from `Latexify.jl` and
`UnitfulLatexify.jl` (for instance `unitformat := :siunitx`, `fmt := "%.2e"`),
as well as a few extra:

# Meta-arguments
* `filename`: The path to a file that will be written.
* `io`: an `IO` object to write to instead of a file. Overrides the `filename`
  argument. If neither this nor `filename` is given, `stdout` is used.
* `permissions`: Defaults to `"w"`. Can be given as `"a"` to append to a file
  instead of overwriting. Only meaningful with a `filename` argument.

# Examples
```julia
julia> a = 2;
julia> b = 3.2u"m";
julia> @datax a b c=3*a d=27 unitformat:=:siunitx
\\pgfkeyssetvalue{/datax/a}{\\num{2}}
\\pgfkeyssetvalue{/datax/b}{\\SI{3.2}{\\meter}}
\\pgfkeyssetvalue{/datax/c}{\\num{6}}
\\pgfkeyssetvalue{/datax/d}{\\num{27}}

```
"""
macro datax(args...)
    esc(datax_helper(args...))
end

function datax_helper(args...)
    metaargs = Expr[]
    names = Symbol[]
    values = Any[]
    for a in args
        if a isa Symbol
            push!(names, a)
            push!(values, a)
            continue
        end
        if a.head == :(=)
            push!(names, a.args[1])
            push!(values, a.args[2])
            continue

        end
        if a.head == :(:=)
            push!(metaargs, Expr(:kw,a.args[1],a.args[2]))
            continue
        end
        error("I don't know what to do with argument $a")
    end

    names = Expr(:tuple, QuoteNode.(names)...)
    values = Expr(:tuple, values...)
    quote
        LaTeXDatax.datax($names,$values;$(metaargs...))
        nothing
    end
end

function datax(names, values; kwargs...)
    if haskey(kwargs, :io)
        datax(kwargs[:io], names, values; kwargs...)
        return nothing
    end
    if haskey(kwargs, :filename)
        datax(kwargs[:filename], names, values; kwargs...)
        return nothing
    end
    datax(stdout, names, values; kwargs...)
    return nothing
end

datax(io::IO, names, values; kwargs...) = printkeyval.(Ref(io),names,values; kwargs...)

function datax(filename::String, names, values; permissions="w", kwargs...)
    if haskey(kwargs,:overwrite)
        if kwargs[:overwrite]
            permissions = "w"
        else
            permissions = "a"
        end
    end
    open(filename, permissions) do io
        permissions == "w" && println(io, "% Autogenerated by LaTeXDatax.jl, will be overwritten")
        datax(io, names, values; kwargs...)
    end
end

function printkeyval(io::IO, name, value; kwargs...)
    print(io, "\\pgfkeyssetvalue{/datax/", name, "}{")
    printdata(io, value; kwargs...)
    print(io, "}\n")
    return nothing
end

printdata(io::IO, v::String; kwargs...) = print(io, v)
printdata(io::IO, v::Number; kwargs...) = print(io, latexify(v*u"one"; kwargs...))
printdata(io::IO, v; kwargs...) = print(io, latexify(v; kwargs...))

end # module


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using Distributions, StatsPlots, Roots

tstat(μ, X) = (mean(X) - μ)/√(var(X) / length(X))
pvalue(μ, X) = 2ccdf(TDist(length(X) - 1), abs(tstat(μ, X)))

function ci(X, α)
    X̄ = mean(X)
    c = quantile(TDist(length(X) - 1), 1 - α/2)
    d = √(var(X) / length(X))
    X̄ - c*d, X̄ + c*d
end

isininterval(x, int) = first(int) ≤ x ≤ last(int)
isinci(μ, X, α) = isininterval(μ, ci(X, α))
prob_trueisinci(dist, n, α; L=10^6) = mean(isinci(mean(dist), rand(dist, n), α) for _ in 1:L)

# %%
ci_roots(X, α) = find_zeros(μ -> log(pvalue(μ, X)) - log(α), -1e2, 1e2)

# %%
@show X = rand(10)
@show ci(X, 0.05)
@show ci_roots(X, 0.05);

# %%
isinci_roots(μ, X, α) = isininterval(μ, ci_roots(X, α))
prob_trueisinci_roots(dist, n, α; L=10^4) = mean(isinci_roots(mean(dist), rand(dist, n), α) for _ in 1:L)

# %%
@time prob_trueisinci_roots(Normal(), 20, 0.05)

# %%

max_stars_repo_path: 0013/probability of true value being in CI (2).jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using Test, OrdinaryDiffEq

using SparseArrays, LinearSolve
using LinearAlgebra, Random
N = 30
AA = sprand(MersenneTwister(12),  N, N, 0.5)
mm = sprand(MersenneTwister(123), N, N, 0.5)
A = DiffEqArrayOperator(AA)
M = DiffEqArrayOperator(mm'mm)
u0 = ones(N)
prob = ODEProblem(ODEFunction(A; mass_matrix=M), u0, (0.0, 1.0))

for alg in [Rosenbrock23(), Rosenbrock23(linsolve=KLUFactorization())]
  sol = solve(prob, alg)
  @test sol.destats.njacs == 0
  @test sol.destats.nw == 1
end

## OOP

foop(u, p, t) = jac(u, p, t) * u
jac(u, p, t) = spdiagm(0=>p)
paramjac(u, p, t) = SparseArrays.spdiagm(0=>u)

n = 2
p = collect(1.0:n)
u0 = ones(n)
tspan = [0.0, 1]
odef = ODEFunction(foop; jac=jac, jac_prototype=jac(u0, p, 0.0), paramjac=paramjac)

function g_helper(p; alg=Rosenbrock23(linsolve=LUFactorization()))
    prob = ODEProblem(odef, u0, tspan, p)
    soln = Array(solve(prob, alg; u0=prob.u0, p=prob.p, abstol=1e-4, reltol=1e-4))[:, end]
    return soln
end
function g(p; kwargs...)
    soln = g_helper(p; kwargs...)
    return sum(soln)
end

@test isapprox(exp.(p), g_helper(p; alg=Rosenbrock23(linsolve=KLUFactorization())); atol=1e-3, rtol=1e-3)
@test isapprox(exp.(p), g_helper(p; alg=Rosenbrock23(linsolve=UMFPACKFactorization())); atol=1e-1, rtol=1e-1)
@test isapprox(exp.(p), g_helper(p; alg=Rosenbrock23(linsolve=KrylovJL_GMRES())); atol=1e-1, rtol=1e-1)

@test isapprox(exp.(p), g_helper(p; alg=Rodas4(linsolve=KLUFactorization())); atol=1e-3, rtol=1e-3)
@test isapprox(exp.(p), g_helper(p; alg=Rodas4(linsolve=UMFPACKFactorization())); atol=1e-1, rtol=1e-1)
@test isapprox(exp.(p), g_helper(p; alg=Rodas4(linsolve=KrylovJL_GMRES())); atol=1e-1, rtol=1e-1)

@test isapprox(exp.(p), g_helper(p; alg=Rodas5(linsolve=KLUFactorization())); atol=1e-3, rtol=1e-3)
@test isapprox(exp.(p), g_helper(p; alg=Rodas5(linsolve=UMFPACKFactorization())); atol=1e-1, rtol=1e-1)
@test isapprox(exp.(p), g_helper(p; alg=Rodas5(linsolve=KrylovJL_GMRES())); atol=1e-1, rtol=1e-1)

@test isapprox(exp.(p), g_helper(p; alg=ImplicitEuler(linsolve=KLUFactorization())); atol=1e-1, rtol=1e-1)
@test isapprox(exp.(p), g_helper(p; alg=ImplicitEuler(linsolve=UMFPACKFactorization())); atol=1e-1, rtol=1e-1)
@test isapprox(exp.(p), g_helper(p; alg=ImplicitEuler(linsolve=KrylovJL_GMRES())); atol=1e-1, rtol=1e-1)

@test isapprox(exp.(p), g_helper(p; alg=TRBDF2(linsolve=KLUFactorization())); atol=1e-1, rtol=1e-1)
@test isapprox(exp.(p), g_helper(p; alg=TRBDF2(linsolve=UMFPACKFactorization())); atol=1e-1, rtol=1e-1)
@test isapprox(exp.(p), g_helper(p; alg=TRBDF2(linsolve=KrylovJL_GMRES())); atol=1e-1, rtol=1e-1)

@test isapprox(exp.(p), g_helper(p; alg=KenCarp47(linsolve=KLUFactorization())); atol=1e-1, rtol=1e-1)
@test isapprox(exp.(p), g_helper(p; alg=KenCarp47(linsolve=UMFPACKFactorization())); atol=1e-1, rtol=1e-1)
@test isapprox(exp.(p), g_helper(p; alg=KenCarp47(linsolve=KrylovJL_GMRES())); atol=1e-1, rtol=1e-1)


text:
# TextDiagnosticPrinter
function TextDiagnosticPrinter(opts::DiagnosticOptions)
    @check_ptrs opts
    status = Ref{CXInit_Error}(CXInit_NoError)
    consumer = clang_TextDiagnosticPrinter_create(opts, status)
    @assert status[] == CXInit_NoError
    return TextDiagnosticPrinter(consumer)
end


text:
using Statistics
using JLD2
using Printf
using GLMakie
using Oceananigans
using Oceananigans.Units
using Oceananigans.Fields: fill_halo_regions!, ReducedField

using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom
using CUDA: @allowscalar

#####
##### Grid
#####

latitude = (-84.375, 84.375)
Δφ = latitude[2] - latitude[1]

# 2.8125 degree resolution
Nx = 128
Ny = 60
Nz = 18

using DataDeps

path = "https://github.com/CliMA/OceananigansArtifacts.jl/raw/main/lat_lon_bathymetry_and_fluxes/"

dh = DataDep("near_global_lat_lon",
    "Forcing data for global latitude longitude simulation", path * "bathymetry_lat_lon_128x60_FP32.bin"
)

DataDeps.register(dh)

datadep"near_global_lat_lon"

bathymetry_data = Array{Float32}(undef, Nx*Ny)
bathymetry_path = @datadep_str "near_global_lat_lon/bathymetry_lat_lon_128x60_FP32.bin"
read!(bathymetry_path, bathymetry_data)

bathymetry_data = bswap.(bathymetry_data) |> Array{Float64}
bathymetry = reshape(bathymetry_data, Nx, Ny)

Nmonths = 12
bytes = sizeof(Float32) * Nx * Ny

H = 3600.0
#bathymetry = - H .* (bathymetry .< -10)

# A spherical domain
underlying_grid = LatitudeLongitudeGrid(size = (Nx, Ny, Nz),
                                        longitude = (-180, 180),
                                        latitude = latitude,
                                        halo = (3, 3, 3),
                                        z = (-H, 0))

grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bathymetry))

function geographic2cartesian(λ, φ; r=1)
    λ = cat(λ, λ[1:1] .+ 360, dims=1)

    Nλ = length(λ)
    Nφ = length(φ)

    λ = repeat(reshape(λ, Nλ, 1), 1, Nφ)
    φ = repeat(reshape(φ, 1, Nφ), Nλ, 1)

    λ_azimuthal = λ .+ 180  # Convert to λ ∈ [0°, 360°]
    φ_azimuthal = 90 .- φ   # Convert to φ ∈ [0°, 180°] (0° at north pole)

    x = @. r * cosd(λ_azimuthal) * sind(φ_azimuthal)
    y = @. r * sind(λ_azimuthal) * sind(φ_azimuthal)
    z = @. r * cosd(φ_azimuthal)

    return x, y, z
end

λu, ϕu, ru = nodes((Face, Center, Center), grid)
λv, ϕv, rv = nodes((Center, Face, Center), grid)
λc, ϕc, rc = nodes((Center, Center, Center), grid)

xu, yu, zu = geographic2cartesian(λu, ϕu)
xv, yv, zv = geographic2cartesian(λv, ϕv)
xc, yc, zc = geographic2cartesian(λc, ϕc)

output_prefix = "annual_cycle_global_lat_lon_128_60_18_temp"

surface_file = jldopen(output_prefix * "_surface.jld2")
bottom_file = jldopen(output_prefix * "_bottom.jld2")

iterations = parse.(Int, keys(surface_file["timeseries/t"]))

iter = Node(0)

# Continents
land_ccc = bathymetry .> - underlying_grid.Δzᵃᵃᶜ / 2

function mask_land_cfc(v)
    land_cfc = cat(land_ccc, land_ccc[:, end:end], dims=2)
    v[land_cfc] .= NaN
    return v
end

function mask_land(data)
    data[land_ccc] .= NaN
    return data
end

ηi(iter) = mask_land(surface_file["timeseries/η/" * string(iter)][:, :, 1])
ui(iter) = mask_land(surface_file["timeseries/u/" * string(iter)][:, :, 1])
vi(iter) = mask_land_cfc(surface_file["timeseries/v/" * string(iter)][:, :, 1])
Ti(iter) = mask_land(surface_file["timeseries/T/" * string(iter)][:, :, 1])
ti(iter) = prettytime(surface_file["timeseries/t/" * string(iter)])

ubi(iter) = mask_land(bottom_file["timeseries/u/" * string(iter)][:, :, 1])
vbi(iter) = mask_land_cfc(bottom_file["timeseries/v/" * string(iter)][:, :, 1])

uri = ReducedField(Face, Center, Nothing, CPU(), grid, dims=3)
vri = ReducedField(Center, Face, Nothing, CPU(), grid, dims=3)
spi = ReducedField(Center, Center, Nothing, CPU(), grid, dims=3)

bathymetry = cat(bathymetry, bathymetry[1:1, :], dims=1)
land = bathymetry .> -10
water = @. !(land)

continents = 1.0 .* land
continents[water] .= NaN

sp = @lift begin
    uri .= ui($iter)
    vri .= vi($iter)
    fill_halo_regions!(uri, CPU())
    fill_halo_regions!(vri, CPU())
    spi .= sqrt(uri^2 + vri^2)
    fill_halo_regions!(spi, CPU())
    sp = spi[1:Nx+1, 1:Ny, 1]
    # sp[land] .= NaN
    return sp
end

η = @lift ηi($iter)
u = @lift ui($iter)
v = @lift vi($iter)
T = @lift Ti($iter)

T_sphere = @lift begin
    T = cat(Ti($iter), Ti($iter)[1:1, :], dims=1)
    T[land] .= NaN
    return T
end

v_sphere = @lift begin
    v = cat(vi($iter)[:, 1:Ny], vi($iter)[1:1, 1:Ny], dims=1)
    v[land] .= NaN
    return v
end

u_sphere = @lift begin
    u = cat(ui($iter), ui($iter)[1:1, :], dims=1)
    u[land] .= NaN
    return u
end

ub = @lift ubi($iter)
vb = @lift vbi($iter)

max_η = 2
min_η = - max_η
max_u = 0.1
min_u = - max_u
max_T = 32
min_T = 0

dλbg = 0.01
dφbg = 0.01
λbg = range(-180, stop=180 - dλbg, step=dλbg)
φbg = range(-89.5, stop=89.5 - dφbg, step=dφbg)
xbg, ybg, zbg = geographic2cartesian(λbg, φbg)

sphere_continents!(ax_T) = surface!(ax_T, xbg, ybg, zbg, color=fill("#080", 1, 1))

#####
##### Sphere
#####

sphere_fig = Figure(resolution = (150, 90))

graylim = 1.0
α = 0.99

#surface!(ax, xc, yc, zc, color=sp, colormap=:blues, colorrange=(0, max_u),
#         show_axis=false, shading=false, ssao=true)

# Temperature
ax_T = sphere_fig[1, 1:3] = LScene(sphere_fig)

surface!(ax_T, xc, yc, zc, color=T_sphere, colormap=:thermal, colorrange=(0, max_T),
         show_axis=false, shading=false, ssao=true)

sphere_continents!(ax_T)

rotate_cam!(ax_T.scene, (π/8, π/6, 0))

# Meridional velocity
ax_v = sphere_fig[1, 4:6] = LScene(sphere_fig)

surface!(ax_v, xc, yc, zc, color=u_sphere, colormap=:balance, colorrange=(min_u, max_u),
         show_axis=false, shading=false, ssao=true)

sphere_continents!(ax_v)

rotate_cam!(ax_v.scene, (π/8, π/6, 0))

plot_title = @lift "t = " * ti($iter)
supertitle = sphere_fig[0, 3:4] = Label(sphere_fig, plot_title, textsize=40)

#v_title = sphere_fig[0, 5] = Label(sphere_fig, "Ocean meridional velocity", textsize=20)
#T_title = sphere_fig[0, 2] = Label(sphere_fig, "Ocean temperature", textsize=20)

display(sphere_fig)

save_interval = 5days
full_rotation_savepoints = round(Int, 10year / save_interval)
dθ = 2π / full_rotation_savepoints
dϕ = π/4 / full_rotation_savepoints

function dϕi(i)
    savepoint = i/iterations[end] * length(iterations)
    mod_savepoint = mod1(savepoint, full_rotation_savepoints)
    halfway_savepoint = 1/2 * length(iterations)
    return mod_i < halfway_savepoint ? dϕ : - dϕ
end

GLMakie.record(sphere_fig, output_prefix * "_sphere.mp4", iterations, framerate=8) do i
    @info "Plotting iteration $i of $(iterations[end])..."
    iter[] = i
    rotate_cam!(ax_T.scene, cameracontrols(ax_T.scene), (dϕi(i), dθ, 0))
    rotate_cam!(ax_v.scene, cameracontrols(ax_v.scene), (dϕi(i), dθ, 0))
end

#=
#####
##### Meridional velocity
#####

plane_fig = Figure(resolution = (1200, 900))

# kwargs = NamedTuple()
kwargs = (; interpolate=true)

function continents!(ax)
    Nλ = length(λc)
    Nφ = length(ϕc)
    λ = repeat(reshape(λc, Nλ, 1), 1, Nφ)
    φ = repeat(reshape(ϕc, 1, Nφ), Nλ, 1)
    surface!(ax, λ, φ, zeros(Nx, Ny) .- 0.02, color=fill("#6e7f80", 1, 1))
    return nothing
end

ax = Axis(plane_fig[1, 1], title="Free surface displacement (m)")
continents!(ax)
hm = heatmap!(ax, λc, ϕc, η; colorrange=(min_η, max_η), colormap=:balance, kwargs...)
cb = Colorbar(plane_fig[1, 2], hm)

ax = Axis(plane_fig[2, 1], title="Sea surface temperature (ᵒC)")
continents!(ax)
hm = heatmap!(ax, λc, ϕc, T; colorrange=(min_T, max_T), colormap=:thermal, kwargs...)
cb = Colorbar(plane_fig[2, 2], hm)

ax = Axis(plane_fig[1, 3], title="East-west surface velocity (m s⁻¹)")
continents!(ax)
hm = heatmap!(ax, λc, ϕc, u; colorrange=(min_u, max_u), colormap=:balance, kwargs...)
cb = Colorbar(plane_fig[1, 4], hm)

ax = Axis(plane_fig[2, 3], title="North-south surface velocity (m s⁻¹)")
continents!(ax)
hm = heatmap!(ax, λc, ϕv, v; colorrange=(min_u, max_u), colormap=:balance, kwargs...)
cb = Colorbar(plane_fig[2, 4], hm)

ax = Axis(plane_fig[3, 1], title="East-west bottom velocity (m s⁻¹)")
continents!(ax)
hm = heatmap!(ax, λc, ϕc, ub; colorrange=(min_u, max_u), colormap=:balance, kwargs...)
cb = Colorbar(plane_fig[3, 2], hm)

ax = Axis(plane_fig[3, 3], title="North-south bottom velocity (m s⁻¹)")
continents!(ax)
hm = heatmap!(ax, λc, ϕv, vb; colorrange=(min_u, max_u), colormap=:balance, kwargs...)
cb = Colorbar(plane_fig[3, 4], hm)

title_str = @lift "Earth day = " * ti($iter)
ax_t = plane_fig[0, :] = Label(plane_fig, title_str)

GLMakie.record(plane_fig, output_prefix * "_plane.mp4", iterations, framerate=8) do i
    @info "Plotting iteration $i of $(iterations[end])..."
    iter[] = i
end

display(plane_fig)
=#

close(surface_file)
close(bottom_file)


text:
module ThriftUtilsTests

using Thrift
using Test
import Thrift: meta
using Base.Threads

struct _enum_TestEnum
    BOOLEAN::Int32
    INT32::Int32
    INT64::Int32
    INT96::Int32
    FLOAT::Int32
    DOUBLE::Int32
    BYTE_ARRAY::Int32
    FIXED_LEN_BYTE_ARRAY::Int32
end
const TestEnum = _enum_TestEnum(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4), Int32(5), Int32(6), Int32(7))

function test_enum()
    @testset "enum" begin
        @test enumstr(TestEnum, TestEnum.BOOLEAN) == "BOOLEAN"
        @test_throws ErrorException enumstr(TestEnum, Int32(11))
    end
end

function test_container_check()
    @testset "iscontainer/isplain" begin
        for T in (Bool, UInt8, Float64, Int16, Int32, Int64, Vector{UInt8}, String)
            @test Thrift.isplain(T)
            @test !Thrift.iscontainer(T)
        end
        for T in (Dict, Set, Any)
            @test !Thrift.isplain(T)
            @test Thrift.iscontainer(T)
        end
    end
end

mutable struct TestMetaAllTypes <: Thrift.TMsg
    meta::ThriftMeta
    values::Dict{Symbol,Any}

    function TestMetaAllTypes(; kwargs...)
        obj = new(__meta__TestMetaAllTypes, Dict{Symbol,Any}())
        values = obj.values
        symdict = obj.meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtype
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        Thrift.setdefaultproperties!(obj)
        obj
    end
end # mutable struct TestMetaAllTypes

const __meta__TestMetaAllTypes = meta(TestMetaAllTypes,
    Symbol[:bool_val,:byte_val,:i16_val,:i32_val,:i64_val,:double_val,:string_val],
    Type[Bool, UInt8, Int16, Int32, Int64, Float64, String],
    Symbol[],
    Int[],
    Dict{Symbol,Any}()
)

function Base.getproperty(obj::TestMetaAllTypes, name::Symbol)
    if name === :bool_val
        return (obj.values[name])::Bool
    elseif name === :byte_val
        return (obj.values[name])::UInt8
    elseif name === :i16_val
        return (obj.values[name])::Int16
    elseif name === :i32_val
        return (obj.values[name])::Int32
    elseif name === :i64_val
        return (obj.values[name])::Int64
    elseif name === :double_val
        return (obj.values[name])::Float64
    elseif name === :string_val
        return (obj.values[name])::String
    else
        getfield(obj, name)
    end
end

meta(::Type{TestMetaAllTypes}) = __meta__TestMetaAllTypes

mutable struct AllTypesDefault <: Thrift.TMsg
    meta::ThriftMeta
    values::Dict{Symbol,Any}

    function AllTypesDefault(; kwargs...)
        obj = new(__meta__AllTypesDefault, Dict{Symbol,Any}())
        values = obj.values
        symdict = obj.meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtype
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        Thrift.setdefaultproperties!(obj)
        obj
    end
end # mutable struct AllTypesDefault

const __meta__AllTypesDefault = meta(AllTypesDefault,
  Symbol[:bool_val,:byte_val,:i16_val,:i32_val,:i64_val,:double_val,:string_val,:map_val,:list_val,:set_val],
  Type[Bool,UInt8,Int16,Int32,Int64,Float64,String,Dict{Int32,Int16},Vector{Int16},Set{UInt8}],
  Symbol[:bool_val,:byte_val,:i16_val,:i32_val,:i64_val,:double_val,:string_val,:map_val,:list_val,:set_val],
  Int[],
  Dict{Symbol,Any}(:bool_val => true, :byte_val => UInt8(1), :i16_val => Int16(10), :i32_val => Int32(20), :i64_val => Int64(30), :double_val => Float64(10.1), :string_val => "hello world", :map_val => Dict(Int32(1) => Int16(10), Int32(2) => Int16(20)), :list_val => Int16[1, 2, 3], :set_val => union!(Set{UInt8}(), UInt8[3, 4, 5]))
)

function Base.getproperty(obj::AllTypesDefault, name::Symbol)
    if name === :bool_val
        return (obj.values[name])::Bool
    elseif name === :byte_val
        return (obj.values[name])::UInt8
    elseif name === :i16_val
        return (obj.values[name])::Int16
    elseif name === :i32_val
        return (obj.values[name])::Int32
    elseif name === :i64_val
        return (obj.values[name])::Int64
    elseif name === :double_val
        return (obj.values[name])::Float64
    elseif name === :string_val
        return (obj.values[name])::String
    elseif name === :map_val
        return (obj.values[name])::Dict{Int32,Int16}
    elseif name === :list_val
        return (obj.values[name])::Vector{Int16}
    elseif name === :set_val
        return (obj.values[name])::Set{UInt8}
    else
        getfield(obj, name)
    end
end

meta(::Type{AllTypesDefault}) = __meta__AllTypesDefault

function test_meta()
    @testset "test metadata" begin
        @test !Thrift.isplain(TestMetaAllTypes)
        @test Thrift.iscontainer(TestMetaAllTypes)

        types = TestMetaAllTypes()
        @test !isfilled(types)
        @test !hasproperty(types, :bool_val)

        types.bool_val = true
        @test hasproperty(types, :bool_val)
        @test !isfilled(types)

        types.byte_val = 1
        types.i16_val = 1
        types.i32_val = 1
        types.i64_val = 1
        types.double_val = 1.1
        types.string_val = "1"

        @test isfilled(types)
        @test isinitialized(types)

        types2 = TestMetaAllTypes()
        @test !isfilled(types2)
        copy!(types2, types)
        @test isfilled(types2)

        types3 = AllTypesDefault()
        @test isfilled(types3)
    end

    nothing
end

function test_parallel_readwrite()
    types = TestMetaAllTypes()
    types.bool_val = true
    types.byte_val = 1
    types.i16_val = 1
    types.i32_val = 1
    types.i64_val = 1
    types.double_val = 1.1
    types.string_val = "1"

    iob = IOBuffer()
    transport = TFileTransport(iob)
    Thrift.write(Thrift.TCompactProtocol(transport), types)
    transport = TFileTransport(IOBuffer(take!(iob)))
    types_read = Thrift.read(Thrift.TCompactProtocol(transport), TestMetaAllTypes)

    for name in propertynames(types)
        @test getproperty(types_read, name) == getproperty(types, name)
    end
end

function test_zigzag()
    testcases = [
        0 => (nbytes=1, encval=0),
        -1 => (nbytes=1, encval=1),
        1 => (nbytes=1, encval=2),
        -2 => (nbytes=1, encval=3),
        2147483647 => (nbytes=5, encval=4294967294),
        -2147483648 => (nbytes=5, encval=4294967295)
    ]

    for T in (Int64, Int32)
        iob = PipeBuffer()
        for kv in testcases
            @test Thrift._write_zigzag(iob, T(kv[1])) == kv[2].nbytes
        end
        for kv in testcases
            read_val = Thrift._read_zigzag(iob, T)
            @test read_val === T(kv[1])
        end
    end
end

@testset "utility functions" begin
    test_enum()
    test_container_check()
    test_meta()
    test_zigzag()
end

@testset "parallel read write" begin
    Threads.@threads for idx in 1:10
        test_parallel_readwrite()
    end

    @sync begin
        for idx in 1:10
            @async test_parallel_readwrite()
        end
    end
end

end


text:
"""
    build_Problem!(OPFpbs, scenario::String)

Return a polynomial problem in complex variables with polynomial constraints
Step 0 : verify that there is at most one infered generator variable for each node
Step 1 : create variables in `OPFpbs[scenario].mp` according to formulations in `OPFpbs[scenario].mp` and add them in optimization problem `pb_opt`
Step 2 : create power balance constraints for each node in optimization problem `pb_opt` according to elements associated to the node
Step 3 : create objective in optimization problem `pb_opt` iterating on all nodal and link elements

# Arguments
- `OPFpbs` : dictionary of scenario (string) associated to structure Scenario
- `scenario` : String which must be a key of OPFpbs

# Output
- `pb_opt` : a polynomial optimization problem

# Example
```
julia > OPFpbs = build_OPFproblems(MatpowerInput, "instances\\matpower\\WB2.m")
julia > OPF = build_Problem!(OPFpbs, "BaseCase")
julia > print(OPF)
▶ variables: VOLT_1 VOLT_2
▶ objective:
(192.3) * conj(VOLT_1) * VOLT_1 + (-96.15384615384615 - 480.7692307692308im) * VOLT_1 * conj(VOLT_2) + (-96.15384615384615 + 480.7692307692308im) * conj(VOLT_1) * VOLT_2
▶ constraints:
LOAD_BUS_2: (-96.15384615384615 - 480.7692307692308im) * conj(VOLT_1) * VOLT_2 + (96.15384615384615 + 480.7692307692308im) * conj(VOLT_2) * VOLT_2 = -350.0f0 + 350.0f0im
UNIT_BUS_1: 0.0f0 - 400.0f0im < (96.15384615384615 + 480.7692307692308im) * conj(VOLT_1) * VOLT_1 + (-96.15384615384615 - 480.7692307692308im) * VOLT_1 * conj(VOLT_2) < 600.0f0 + 400.0f0im
   VOLTM_1: 0.9025 < conj(VOLT_1) * VOLT_1 < 1.1025
   VOLTM_2: 0.9025 < conj(VOLT_2) * VOLT_2 < 1.056784

```
"""

function build_Problem!(OPFpbs#=::OPFProblems=#, scenario::String)

  ds = OPFpbs[scenario].ds
  gs = OPFpbs[scenario].gs
  mp = OPFpbs[scenario].mp

  # 1. Creation des variables
  for bus in get_buses(OPFpbs, scenario)
    bus_elems = mp.node_formulations[bus]
    bus_vars = mp.node_vars[bus]
    for (elemid, elem) in ds.bus[bus]
      create_vars!(elem, bus, elemid, bus_elems[elemid], bus_vars, scenario)
    end
  end

  for link in get_links(OPFpbs, scenario)
    link_elems = mp.link_formulations[link]
    link_vars = mp.link_vars[link]
    for (elemid, elem) in ds.link[link]
      create_vars!(elem, link, elemid, link_elems[elemid], link_vars, scenario)
    end
  end

  # 2. Création des bilans
  pb_opt = Problem()
  for (node, vars) in mp.node_vars
    for (elem, var) in vars
      add_variable!(pb_opt, var)
    end
  end

  for (link, vars) in mp.link_vars
    for (elem, var) in vars
      add_variable!(pb_opt, var)
    end
  end

  for (busid, bus_elems) in ds.bus
    bus_elems_formulations = mp.node_formulations[busid]
    bus_elems_var = mp.node_vars[busid]

    cstrnames, Snode, lb, ub = get_Snodal(busid, bus_elems, bus_elems_formulations, bus_elems_var)
    S_balance = Snode + get_Slinks_in(busid, ds, gs, mp) + get_Slinks_out(busid, ds, gs, mp)

    add_constraint!(pb_opt, cstrname_nodal_balance(cstrnames, scenario, busid), lb << S_balance << ub)
  end

  ## 3. Création des contraintes
  for bus in get_buses(OPFpbs, scenario)
    bus_elems_formulations = mp.node_formulations[bus]
    bus_vars = mp.node_vars[bus]
    for (elemid, elem) in ds.bus[bus]
      elem_formulation = bus_elems_formulations[elemid]
      constraints = constraint(elem, bus, elemid, elem_formulation, bus_vars, scenario, OPFpbs)
      for (cstrname, cstr) in constraints
    	   add_constraint!(pb_opt, get_cstrname(scenario, bus, elemid, cstrname), cstr)
      end
    end
  end

  for link in get_links(OPFpbs, scenario)
    link_elems_formulations = mp.link_formulations[link]
    link_vars = mp.link_vars[link]
    for (elemid, elem) in ds.link[link]
      elem_formulation = link_elems_formulations[elemid]
      constraints = constraint(elem, link, elemid, elem_formulation, link_vars, scenario, OPFpbs)
      for (cstrname, cstr) in constraints
        add_constraint!(pb_opt, get_cstrname(scenario, link, elemid, cstrname), cstr)
      end
    end
  end

    ## 4. Création de l'objectif
    objective = Polynomial()
    for bus in get_buses(OPFpbs, scenario)
      bus_elems = ds.bus[bus]
      bus_elems_formulations = mp.node_formulations[bus]
      bus_elems_var = mp.node_vars[bus]

      _, Snode, lb, ub = get_Snodal(bus, bus_elems, bus_elems_formulations, bus_elems_var)
      add!(Snode, get_Slinks_in(bus, ds, gs, mp) + get_Slinks_out(bus, ds, gs, mp))

      for (elemid, elem) in bus_elems
        gencost::Polynomial = cost(elem, bus, elemid, bus_elems_formulations[elemid], bus_elems_var, Snode, lb, ub)
        add!(objective, gencost)
      end
    end

    set_objective!(pb_opt, objective)
  return pb_opt
end

"""
    get_Snodal(bus::String, bus_elems::Dict{String, Any}, bus_elems_formulations::Dict{String, Symbol}, bus_elems_var::Dict{String, Variable})

Return the total injected power at `bus` coming from nodal elements with bounds. The injected power for each specific element is computed using `Snodal`.

# Arguments
- `bus::String` :
- `bus_elems::Dict{String,Any}`:
- `bus_elems_formulations` :
- `bus_elems_var` :

# Output

#Examples
```
Instance Matpower WB2 with two nodes
julia > get_Snodal("BUS_1", )
```
"""
function get_Snodal(bus::String, bus_elems::Dict{String, Any}, bus_elems_formulations::Dict{String, Symbol}, bus_elems_var::Dict{String, Variable})
  cstrnames, Sres, lbres, ubres = [Set{String}(), Polynomial(), 0, 0]
  for (elemid, element) in bus_elems
    elem_formulation = bus_elems_formulations[elemid]
    cstrname_, Sres_, lbres_, ubres_ = Snodal(element, bus, elemid, elem_formulation, bus_elems_var)
    Sres, lbres, ubres = [Sres, lbres, ubres] + [Sres_, lbres_, ubres_]
    push!(cstrnames, cstrname_)
  end

  # Remove all constraints names equal to ""
  dense_cstrnames = Set{String}()
  for cstr in cstrnames
    if cstr != ""
      push!(dense_cstrnames, cstr)
    end
  end
  return [dense_cstrnames, Sres, lbres, ubres]
end

function get_Slinks_in(bus, ds::DataSource, gs::GridStructure, mp::MathematicalProgramming)
  Sres = Polynomial()
  for link in get_links_in(bus, gs)
    link_elems_formulation = mp.link_formulations[link]
    link_vars = mp.link_vars[link]
    for (elemid, elem) in ds.link[link]
      elem_formulation = link_elems_formulation[elemid]
      add!(Sres, Sdest(elem, link, elemid, elem_formulation, link_vars))
    end
  end
  return Sres
end

function get_Slinks_out(bus, ds::DataSource, gs::GridStructure, mp::MathematicalProgramming)
  Sres = Polynomial()
  for link in get_links_out(bus, gs)
    link_elems_formulations = mp.link_formulations[link]
    link_vars = mp.link_vars[link]
    for (elemid, elem) in ds.link[link]
      elem_formulation = link_elems_formulations[elemid]
      add!(Sres, Sorig(elem, link, elemid, elem_formulation, link_vars))
    end
  end
  return Sres
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
# UCEvaluation: UC Evaluation by scenarios
function m_UCEval(v, pexp, rup, rdw, RESData)

    NetLoad = collect(sum(Pd,dims=1)[t] - sum(mean(RESData[t,:])) for t in T)

    mUCEval = Model(() -> Gurobi.Optimizer(GUROBI_ENV))
        set_optimizer_attribute(mUCEval, "OutputFlag", 0)

        # Variables
        @variables(mUCEval, begin
            0. ≤ p[g in G,T] ≤ GenData.Pmax[g]
            0. ≤ r_p[G,T]
            0. ≤ r_m[G,T]
            0. ≤ l_s[t in T] ≤ sum(Pd,dims=1)[t] # load shedding
            0. ≤ w_s[t in T]  # wind spillage
        end)

        # Constraints
        # 2nd stage
        # Power balance
        @constraint(mUCEval, PowerBalance_Second[t in T],
                            NetLoad[t] - l_s[t]== sum(p[:,t]) - w_s[t])

            # Power limits
            @constraints(mUCEval, begin
                PMax2nd[g in G, t in T], p[g,t] ≤ v[g,t]*GenData.Pmax[g]
                PMin2nd[g in G, t in T], p[g,t] ≥ v[g,t]*GenData.Pmin[g]
                RplusLim[g in G, t in T], r_p[g,t] ≤ rup[g,t]
                RminusLim[g in G, t in T], r_m[g,t] ≤ rdw[g,t]
                Ramps_2ndUp[g in G, t in T], p[g,t] - p[g,T[t-1]] ≤ GenData.Rup[g]
                Ramps_2ndDw[g in G, t in T], -GenData.Rdw[g] ≤ p[g,t] - p[g,T[t-1]]
                Spillage_Limits[t in T], w_s[t] ≤ mean(RESData[t,:])

                end)

            # Generated Power
            @constraint(mUCEval, PGen_2nd[g in G, t in T],
                    p[g,t] - r_p[g,t] + r_m[g,t] == pexp[g,t])

        ### Objective
        @objective(mUCEval, Min, sum(C_ls*l_s[t] + C_ws*w_s[t] + sum(GenData.Cplus[g]*r_p[g,t] - GenData.Cminus[g]*r_m[g,t] for g in G) for t in T) )

    optimize!(mUCEval)
    return mUCEval
end

max_stars_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_issues_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_forks_repo_name: gonzalez-alvaro/MultidimensionalGAPM


text:
module OpenStreetMapPBF
using Logging

include("proto/OSMPBF_pb.jl")
include("model.jl")
include("scan.jl")
include("load.jl")

export scan_pbf, scan_nodes, scan_ways, scan_relations, load_pbf, Node, Way, Relation, way, node, relation
end


text:
# # Distribution of calibration error estimates
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/distribution.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/distribution.ipynb)

# ## Introduction
#
# This example is taken from the publication
# "Calibration tests in multi-class classification: A unifying framework" by Widmann,
# Lindsten, and Zachariah.
#
# We estimate calibration errors of the model
# ```math
# \begin{aligned}
#    g(X) &\sim \textrm{Dir}(\alpha),\\
#    Z &\sim \textrm{Ber}(\pi),\\
#    Y \,|\, g(X) = \gamma, Z = 1 &\sim \textrm{Categorical}(\beta),\\
#    Y \,|\, g(X) = \gamma, Z = 0 &\sim \textrm{Categorical}(\gamma),
# \end{aligned}
# ```
# where $\alpha \in \mathbb{R}_{>0}^m$ determines the distribution of
# predictions $g(X)$, $\pi > 0$ determines the degree of miscalibration, and
# $\beta$ defines a fixed categorical distribution.
#
# Here we consider only the choices $\alpha = (0.1, \ldots, 0.1)$, mimicking a
# distribution after training that is pushed towards the edges of the
# probability simplex, and $\beta = (1, 0, \ldots, 0)$.
#
# In our experiments we sample 250 predictions from the Dirichlet distribution
# $\textrm{Dir}(\alpha)$, and then we generate corresponding labels according to
# the model stated above, for different choices of $\pi$ and number of classes $m$.
#
# We evaluate the standard estimators of expected calibration error (ECE) based on a
# uniform binning scheme and a data-dependent binning scheme, and the biased estimator of
# the squared kernel calibration error (SKCE), the quadratic unbiased estimator of
# the SKCE, and the linear unbiased estimator of the SKCE for a specific choice of
# matrix-valued kernels.
#
# The sampling procedure and the evaluation are repeated 100 times, to obtain a sample
# of 100 estimates for each considered setting of $\pi$ and $m$.
#
# For our choice of $\alpha$ and $\beta$, the analytical ECE with respect to the
# total variation distance $\|.\|_{\mathrm{TV}}$ is
# ```math
# \mathbb{E}[\|.\|_{\mathrm{TV}},g] = \frac{\pi(m-1)}{m}.
# ```

# ## Setup

using CalibrationErrors
using Distances
using Distributions
using StatsBase

using LinearAlgebra
using Random

using Plots
gr(fmt = :png, dpi = 600)

# ## Estimates
#
#

function estimates(rng::AbstractRNG, estimator, π::Real, m::Int)
    ## check arguments
    m > 0 || throw(ArgumentError("number of classes must be positive"))
    zero(π) <= π <= one(π) ||
        throw(ArgumentError("probability π must be between 0 and 1"))

    ## cache array for predictions, modified predictions, and labels
    predictions = Matrix{Float64}(undef, m, 250)
    labels = Vector{Int}(undef, 250)
    data = (predictions, labels)

    ## define sampler of predictions
    sampler_predictions = sampler(Dirichlet(m, 0.1))

    ## initialize estimates
    estimates = Vector{Float64}(undef, 100)

    ## for each run
    @inbounds for i in eachindex(estimates)
        ## sample predictions
        rand!(rng, sampler_predictions, predictions)

        ## sample labels
        @inbounds for j in eachindex(labels)
            if rand(rng) < π
                labels[j] = 1
            else
                labels[j] = sample(rng, Weights(view(predictions, :, j), 1))
            end
        end

        ## evaluate estimator
        estimates[i] = estimator isa CalibrationErrors.CalibrationErrorEstimator ?
            calibrationerror(estimator, data) : calibrationerror(estimator(data), data)
    end

    estimates
end

# We use a helper function to run the experiment for all desired parameter settings.

struct EstimatesSet{E}
    estimator::E
    m::Vector{Int}
    π::Vector{Float64}
    estimates::Vector{Vector{Float64}}
end

function estimates(rng::AbstractRNG, estimator)
    ## create arrays
    mvec = Vector{Int}(undef, 0)
    πvec = Vector{Float64}(undef, 0)
    estimatesvec = Vector{Vector{Float64}}(undef, 0)

    ## for all combinations of m and π
    for m in (2, 10, 100), π in (0.0, 0.5, 1.0)
        ## compute estimates
        push!(estimatesvec, estimates(Random.GLOBAL_RNG, estimator, π, m))

        ## save m and π
        push!(mvec, m)
        push!(πvec, π)
    end

    EstimatesSet(estimator, mvec, πvec, estimatesvec)
end

estimates(estimator) = estimates(Random.GLOBAL_RNG, estimator)

# As mentioned above, we can calculate the analytic expected calibration error. For the squared
# kernel calibration error, we take the mean of the estimates of the unbiased quadratic
# estimator as approximation of the true value.

# We provide simple histogram plots of our results. The mean value of the
# estimates is indicated by a solid black vertical line and the analytic
# calibration error is visualized as a dashed red line.

@recipe function f(set::EstimatesSet)
    ## default settings
    layout := (3, 3)
    legend := false
    xlabel := "calibration estimate"
    ylabel := "# runs"
    size --> (1080, 960)
    seriestype := :histogram

    ## add subplots
    for i in 1:length(set.estimates)
        ## retrieve data
        m = set.m[i]
        π = set.π[i]
        estimates = set.estimates[i]

        ## plot histogram of estimates
        @series begin
            subplot := i
            x := estimates
        end

        ## indicate analytic calibration error for ECE
        if set.estimator isa ECE
            @series begin
                seriestype := :vline
                subplot := i
                linewidth := 2
                color := :red
                x := [π * (m - 1) / m]
            end
        end

        ## indicate mean of estimates
        @series begin
            title := "$m classes, pi = $π"
            seriestype := :vline
            subplot := i
            linewidth := 2
            color := :black
            x := [mean(estimates)]
        end
    end
end

# ## Kernel choice

# We use an identity matrix scaled with an exponential kernel as matrix-valued
# kernel. The bandwidth of the exponential kernel is set according to the
# median heuristic.

function kernel(data)
    ## compute median TV distance
    predictions = data[1]
    γ = inv(median(pairwise(TotalVariation(), predictions, dims=2)))

    ## create kernel
    UniformScalingKernel(ExponentialKernel(γ, TotalVariation()))
end

# ## Expected calibration error
#
# ### Uniform binning
#
# We start by analyzing the expected calibration error (ECE).
# For our estimation we use 10 bins of uniform width in each dimension.

Random.seed!(1234)
data = estimates(ECE(UniformBinning(10), TotalVariation()))
plot(data)

# ### Non-uniform binning
#
# We repeat our experiments with a different data-dependent binning scheme. This
# time the bins will be computed dynamically by splitting the predictions at the
# median of the classes with the highest variance, as long as the number of bins
# does not exceed a given threshold and the number of samples per bin is above
# a certain lower bound. In our experiments we do not impose any restriction on
# the number of bins but only stop splitting if the number of samples is less
# than 10.

Random.seed!(1234)
data = estimates(ECE(MedianVarianceBinning(10), TotalVariation()))
plot(data)

# ## Biased estimator of the squared kernel calibration error
#

Random.seed!(1234)
data = estimates(x -> BiasedSKCE(kernel(x)))
plot(data)

# ## Unbiased quadratic estimator of the squared kernel calibration error
#

Random.seed!(1234)
data = estimates(x -> QuadraticUnbiasedSKCE(kernel(x)))
plot(data)

# ## Unbiased linear estimator of the squared kernel calibration error
#

Random.seed!(1234)
data = estimates(x -> LinearUnbiasedSKCE(kernel(x)))
plot(data)


text:
function sesam3C_protect(;name, k_N=60.0)
    # extends sesam3C by assuming only a portion of residues accessible to enzymatic
    # degradation. That portion is in Langmuir-equilibrium with a protected pool
    @parameters t
    D = Differential(t)

    @parameters ϵ ϵ_tvr κ_E a_E  m  τ
    @parameters k_L  k_R  k_mN_L k_mN_R k_N=k_N
    @parameters α_R0
    @parameters ρ_NBtvr = 0.0 # used in N component, only relevant for variants implementing microbial nutrient resorption
    @parameters K_eqR=1e-3 Q_max=10_000.0

    @variables (begin
        B(t),  L(t),   R(t),  cumresp(t), Ra(t), # accessible portion of R
        dB(t), dL(t), dR(t), r_tot(t),
        syn_Enz(t), tvr_Enz(t), r_M(t), tvr_B(t), tvr_B0(t),
        E_L(t), E_R(t),
        dec_LPot(t), dec_L(t), dec_RPot(t),
        dec_R(t), u_C(t),
        C_synBCt(t), C_synBC(t), r_G(t),
        r_tvr(t),
        r_B(t), r_GEnz(t), r_O(t),
        # need to be defined by component across all elements
        α_L(t), α_R(t),
        # need to be specified by coupled system:
        i_L(t), syn_B(t),
        k_Rtot(t)
    end)

    eqs = [
        D(B) ~ dB, dB ~ syn_B - tvr_B,
        D(L) ~ dL, dL ~ -dec_L + i_L,
        D(R) ~ dR, dR ~ -dec_R + ϵ_tvr*tvr_B + (1-κ_E)*syn_Enz,
        syn_Enz ~ a_E*B, tvr_Enz ~ syn_Enz,
        r_M ~ m*B,
        tvr_B ~ τ*B,
        tvr_B0 ~ tvr_B, # used in N component, differs from tvr_B in microbila nutrient resorption variants
        dec_LPot ~ k_L * L,
        dec_L ~ dec_LPot*(α_L * syn_Enz)/(k_mN_L + α_L*syn_Enz),
        # the next two lines are different from sesam3
        # note the new observable Ra and paramters Q_max, K_eqR
        dec_RPot ~ k_R * Ra,
        Ra ~ compute_accessible_langmuir(R, Q_max, K_eqR),
        k_Rtot ~ k_R * Ra/R,
        dec_R ~ dec_RPot*(α_R * syn_Enz)/(k_mN_R + α_R*syn_Enz),
        u_C ~ dec_L + dec_R + κ_E*syn_Enz,
        C_synBCt ~ u_C - syn_Enz/ϵ - r_M,
        #C_synBC ~ IfElse.ifelse(C_synBCt > 0.0, ϵ*C_synBCt, C_synBCt),
        C_synBC ~ C_synBCt - ((C_synBCt > 0.0) * (1-ϵ)*C_synBCt),
        #syn_B ~ C_synBC # define in all-elements model
        r_O ~ u_C - (syn_Enz/ϵ +  syn_B + r_G + r_M),
        #D(cumresp) ~ r_B + r_tvr,
        r_tvr ~ (1-ϵ_tvr)*tvr_B,
        r_B ~ r_GEnz + r_G + r_M + r_O,
        r_GEnz ~ (1-ϵ)/ϵ * syn_Enz,
        #r_G ~ IfElse.ifelse(syn_B.val > 0.0, (1-ϵ)/ϵ * syn_B, 0.0),
        r_G ~ (syn_B.val > 0.0) * (1-ϵ)/ϵ * syn_B,
        r_tot ~ r_B + r_tvr,
        E_L ~ (α_L * syn_Enz)/k_N,
        E_R ~ (α_R * syn_Enz)/k_N,
        ]
    ODESystem(eqs; name)
end

function compute_accessible_langmuir(R, Q_max, K_eqR)
    # positive solution of a quadratic equation
    # see inst/pluto/R_langmuir.jl
    # provide Q_max and K_eqR as Ref-values to not interfere with symbolics when registering
    phalv = (Q_max - R + 1/K_eqR)/2
	q = -R/K_eqR
	-phalv + sqrt(phalv * phalv - q)
end
@register_symbolic compute_accessible_langmuir(t, Q_max, K_eqR)

sesam3_protect(args...;kwargs...) = sesam3CN(
    args...;
    sN = sN=sesam3N(name=:sN, sC = sesam3C_protect(name=:sC)),
    kwargs...
)


text:
export alias_elimination, flatten

function flatten(sys::ODESystem)
    if isempty(sys.systems)
        return sys
    else
        return ODESystem(equations(sys),
                         independent_variable(sys),
                         states(sys),
                         parameters(sys),
                         observed=observed(sys))
    end
end

using SymbolicUtils: Rewriters

function fixpoint_sub(x, dict)
    y = substitute(x, dict)
    while !isequal(x, y)
        y = x
        x = substitute(y, dict)
    end

    return x
end

function substitute_aliases(diffeqs, dict)
    lhss(diffeqs) .~ fixpoint_sub.(rhss(diffeqs), (dict,))
end

# Note that we reduce parameters, too
# i.e. `2param = 3` will be reduced away
isvar(s::Sym) = true
isvar(s::Term) = isvar(operation(s))
isvar(s::Any) = false

function get_α_x(αx)
    if isvar(αx)
        return 1, αx
    elseif istree(αx) && operation(αx) === (*)
        args = arguments(αx)
        nums = []
        syms = []
        for arg in args
            isvar(arg) ? push!(syms, arg) : push!(nums, arg)
        end

        if length(syms) == 1
            return prod(nums), syms[1]
        end
    else
        return nothing
    end
end

function alias_elimination(sys::ODESystem)
    eqs = vcat(equations(sys), observed(sys))
    neweqs = Equation[]; sizehint!(neweqs, length(eqs))
    subs = Pair[]
    diff_vars = filter(!isnothing, map(eqs) do eq
            if isdiffeq(eq)
                arguments(eq.lhs)[1]
            else
                nothing
            end
        end) |> Set

    deps = Set()
    for (i, eq) in enumerate(eqs)
        # only substitute when the variable is algebraic
        if isdiffeq(eq)
            push!(neweqs, eq)
            continue
        end

        maybe_alias = isalias = false
        res_left = get_α_x(eq.lhs)
        if !isnothing(res_left) && !(res_left[2] in diff_vars)
            # `α x = rhs` => `x = rhs / α`
            α, x = res_left
            sub = x => _isone(α) ? eq.rhs : eq.rhs / α
            maybe_alias = true
        else
            res_right = get_α_x(eq.rhs)
            if !isnothing(res_right) && !(res_right[2] in diff_vars)
                # `lhs = β y` => `y = lhs / β`
                β, y = res_right
                sub =  y => _isone(β) ? eq.lhs : β * eq.lhs
                maybe_alias = true
            end
        end

        if maybe_alias
            l, r = sub
            # alias equations shouldn't introduce cycles
            if !(l in deps) && isempty(intersect(deps, vars(r)))
                push!(deps, l)
                push!(subs, sub)
                isalias = true
            end
        end

        if !isalias
            neweq = _iszero(eq.lhs) ? eq : 0 ~ eq.rhs - eq.lhs
            push!(neweqs, neweq)
        end
    end

    eqs′ = substitute_aliases(neweqs, Dict(subs))

    alias_vars = first.(subs)
    sys_states = states(sys)
    alias_eqs = alias_vars .~ last.(subs)
    #alias_eqs = topsort_equations(alias_eqs, sys_states)

    newstates = setdiff(sys_states, alias_vars)
    ODESystem(eqs′, sys.iv, newstates, parameters(sys), observed=alias_eqs)
end

"""
$(SIGNATURES)

Use Kahn's algorithm to topologically sort observed equations.

Example:
```julia
julia> @variables t x(t) y(t) z(t) k(t)
(t, x(t), y(t), z(t), k(t))

julia> eqs = [
           x ~ y + z
           z ~ 2
           y ~ 2z + k
       ];

julia> ModelingToolkit.topsort_equations(eqs, [x, y, z, k])
3-element Vector{Equation}:
 Equation(z(t), 2)
 Equation(y(t), k(t) + 2z(t))
 Equation(x(t), y(t) + z(t))
```
"""
function topsort_equations(eqs, states; check=true)
    graph, assigns = observed2graph(eqs, states)
    neqs = length(eqs)
    degrees = zeros(Int, neqs)

    for 𝑠eq in 1:length(eqs); var = assigns[𝑠eq]
        for 𝑑eq in 𝑑neighbors(graph, var)
            # 𝑠eq => 𝑑eq
            degrees[𝑑eq] += 1
        end
    end

    q = Queue{Int}(neqs)
    for (i, d) in enumerate(degrees)
        d == 0 && enqueue!(q, i)
    end

    idx = 0
    ordered_eqs = similar(eqs, 0); sizehint!(ordered_eqs, neqs)
    while !isempty(q)
        𝑠eq = dequeue!(q)
        idx+=1
        push!(ordered_eqs, eqs[𝑠eq])
        var = assigns[𝑠eq]
        for 𝑑eq in 𝑑neighbors(graph, var)
            degree = degrees[𝑑eq] = degrees[𝑑eq] - 1
            degree == 0 && enqueue!(q, 𝑑eq)
        end
    end

    (check && idx != neqs) && throw(ArgumentError("The equations have at least one cycle."))

    return ordered_eqs
end

function observed2graph(eqs, states)
    graph = BipartiteGraph(length(eqs), length(states))
    v2j = Dict(states .=> 1:length(states))

    # `assigns: eq -> var`, `eq` defines `var`
    assigns = similar(eqs, Int)

    for (i, eq) in enumerate(eqs)
        lhs_j = get(v2j, eq.lhs, nothing)
        lhs_j === nothing && throw(ArgumentError("The lhs $(eq.lhs) of $eq, doesn't appear in states."))
        assigns[i] = lhs_j
        vs = vars(eq.rhs)
        for v in vs
            j = get(v2j, v, nothing)
            j !== nothing && add_edge!(graph, i, j)
        end
    end

    return graph, assigns
end


text:
@doc_str """
    transitiveclosure!(g, selflooped=false)

Compute the transitive closure of a directed graph, using the Floyd-Warshall
algorithm. If `selflooped` is true, add self loops to the graph.

### Performance
Time complexity is \\mathcal{O}(|V|^3).

### Implementation Notes
This version of the function modifies the original graph.
"""
function transitiveclosure! end
@traitfn function transitiveclosure!(g::::IsDirected, selflooped=false)
    for k in vertices(g)
        for i in vertices(g)
            i == k && continue
            for j in vertices(g)
                j == k && continue
                if (has_edge(g, i, k) && has_edge(g, k, j))
                    if (i != j || selflooped)
                        add_edge!(g, i, j)
                    end
                end
            end
        end
    end
    return g
end

"""
    transitiveclosure(g, selflooped=false)

Compute the transitive closure of a directed graph, using the Floyd-Warshall
algorithm. Return a graph representing the transitive closure. If `selflooped`
is `true`, add self loops to the graph.

### Performance
Time complexity is \\mathcal{O}(|V|^3).
"""
function transitiveclosure(g::DiGraph, selflooped = false)
    copyg = copy(g)
    return transitiveclosure!(copyg, selflooped)
end


text:
#######################################
# UdpTest module name
#######################################
module UdpTest
    using Sockets
    import Base.show

    export server, client, show_info

    function __init__()
        inizia()
    end

    function inizia()
        println("Init")
        global sock = UDPSocket()
        global port=5005
        for a in ARGS
            b=split(a)
            if occursin("-p",a) || occursin("--port",a)
                try
                    port=split(a)[2]
                catch
                    println("Using defalut port $port")
                end
            end

        end
    end

    function show_info(message, remote, prefix)
        try
            prefix=prefix
            remote_ip=remote
            remote_port="-"#remote[2]
            message=message
            info = "$prefix $remote_ip : $remote_port ==> $message"
            println(info)
            return 0
        catch e
            println(e)
            return 1
        end
    end
"""
function server(port::Int)
    start to listen on port and reply any Received message
"""
    function server(port)
        try
            bind(sock , IPv4(0), port)
        @async while true
                    remote, message = recvfrom(sock)
                    println("messaggio ricevuto $remote")
                    show_info(message, remote, "Get a message from")
                    send(sock,remote,port,message)
                end
        catch e
            println(e)
            return 1
        end
    end

"""
    function client(host::String, port::Int)
    require a message to send to the HOST on the PORT.
        Eventually show the received replay.
"""
    function client(host, port)
        try
            bind(sock, ip"0.0.0.0", port)
            host_ip = getaddrinfo(host)
            while true
                print("Send Message: ")
                line = readline()
                send(sock,host_ip,port,line)
                remote, line = recvfrom(sock)
                show_info(line, host_ip, "Received reply from")
            end
        catch e
            println(e)
            return 1
        end
    end

end # module


text:
struct GLMResult{S, T}
    x::AbstractVector{S}
    y::AbstractVector{T}
    l::AbstractVector{T}
    u::AbstractVector{T}
end

function linear(x::AbstractVector{T1}, y::AbstractVector{T2}; n_points = 100) where {T1, T2}
    try
        lin_model = GLM.lm(@formula(Y ~ X), (X=x, Y=y))
        x_min, x_max = extrema(x)
        x_new = range(x_min, x_max, length = n_points)
        y_new, lower, upper = GLM.predict(lin_model,
         [ones(T1, n_points) x_new], interval=:confidence)
        # the GLM predictions always return matrices
       return GLMResult(x_new, vec(y_new), vec(lower), vec(upper))
    catch e
        error("Linear fit not possible for the given data")
    end
end

convert_arguments(P::PlotFunc, g::GLMResult) = PlotSpec{ShadedLine}(g.x, g.y, g.l, g.u)

struct Smooth{S, T}
    x::Vector{S}
    y::Vector{T}
end

convert_arguments(P::PlotFunc, l::Smooth) = PlotSpec{Lines}(Point2f0.(l.x,l.y))

function smooth(x, y; length = 100, kwargs...)
    model = loess(x, y; kwargs...)
    min, max = extrema(x)
    us = collect(range(min, stop = max, length = length))
    vs = Loess.predict(model, us)
    return Smooth(us, vs)
end

smooth(; kwargs...) = (args...) -> smooth(args...; kwargs...)

max_stars_repo_name: piever/StatMakie

max_issues_repo_name: piever/StatMakie

max_forks_repo_name: piever/StatMakie


text:
using Documenter, PiCamera

makedocs(
    modules = [PiCamera],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "yakir12",
    sitename = "PiCamera.jl",
    pages = Any["index.md"]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(
    repo = "github.com/yakir12/PiCamera.jl.git",
)


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

struct DistributionDensity{
    D<:ContinuousDistribution,
    B<:AbstractVarBounds
} <: DistLikeDensity
    dist::D
    bounds::B
end

DistributionDensity(d::Distribution) = DistributionDensity(d, dist_param_bounds(d))

Base.convert(::Type{AbstractDensity}, d::ContinuousDistribution) = DistributionDensity(d)
Base.convert(::Type{DistLikeDensity}, d::ContinuousDistribution) = DistributionDensity(d)

Base.convert(::Type{Distribution}, d::DistributionDensity) = d.dist
Base.convert(::Type{ContinuousDistribution}, d::DistributionDensity) = d.dist

Base.parent(density::DistributionDensity) = density.dist

Base.:(==)(a::DistributionDensity, b::DistributionDensity) = a.dist == b.dist && a.bounds == b.bounds

function eval_logval_unchecked(density::DistributionDensity{<:Distribution{Univariate,Continuous}}, v::Real)
    d = density.dist
    logd = logpdf(d, v)
    R = typeof(logd)
    # ToDo: Move these workarounds somewhere else? Still necessary at all?
    if isnan(logd)
        if isinf(v)
            # Weibull yields NaN logpdf at infinity (Distributions.jl issue #1197), possibly others too,
            # so force to -Inf (there should never be any probability mass at infinity):
            convert(R, -Inf)
        elseif v ≈ minimum(d)
            # Weibull yields NaN logpdf at 0 (Distributions.jl issue #1197), possibly others too,
            # so move an epsilon away from minimum:
            convert(R, logpdf(d, minimum(d) + eps(typeof(v))))
        elseif v ≈ maximum(d)
            # Likewise at maxiumum:
            convert(R, logpdf(d, maximum(d) - eps(typeof(v))))
        else
            logd
        end
    else
        logd
    end
end

eval_logval_unchecked(density::DistributionDensity, v::Any) = Distributions.logpdf(density.dist, v)

ValueShapes.varshape(density::DistributionDensity) = varshape(density.dist)

ValueShapes.unshaped(density::DistributionDensity) = DistributionDensity(unshaped(density.dist))

(shape::AbstractValueShape)(density::DistributionDensity) = DistributionDensity(shape(density.dist))

# For user convenience, don't use within BAT:
@inline Random.rand(rng::AbstractRNG, density::DistributionDensity) = rand(rng, density.dist)
@inline Random.rand(rng::AbstractRNG, density::DistributionDensity, dims::Dims) = rand(rng, density.dist, dims)
@inline Random.rand(rng::AbstractRNG, density::DistributionDensity, dims::Integer...) = rand(rng, density.dist, dims...)

Distributions.sampler(density::DistributionDensity) = Distributions.sampler(density.dist)
bat_sampler(density::DistributionDensity) = bat_sampler(density.dist)

Statistics.cov(density::DistributionDensity{<:MultivariateDistribution}) = cov(density.dist)

var_bounds(density::DistributionDensity) = density.bounds

dist_param_bounds(d::Distribution{Univariate,Continuous}) =
    HyperRectBounds([minimum(d)], [maximum(d)])

dist_param_bounds(d::Distribution{Multivariate,Continuous}) =
    HyperRectBounds(fill(_default_PT(-Inf), length(d)), fill(_default_PT(+Inf), length(d)))

dist_param_bounds(d::ReshapedDist) = dist_param_bounds(unshaped(d))

dist_param_bounds(d::StandardMvUniform) =
    HyperRectBounds(fill(_default_PT(Float32(0)), length(d)), fill(_default_PT(Float32(1)), length(d)))

dist_param_bounds(d::Product{Continuous}) =
    HyperRectBounds(minimum.(d.v), maximum.(d.v))

dist_param_bounds(d::ConstValueDist) = HyperRectBounds(Int32[], Int32[])

dist_param_bounds(d::NamedTupleDist) = vcat(map(x -> dist_param_bounds(x), values(d))...)
dist_param_bounds(d::ValueShapes.UnshapedNTD) = dist_param_bounds(d.shaped)

dist_param_bounds(d::HierarchicalDistribution) =
    HyperRectBounds(fill(_default_PT(-Inf), length(d)), fill(_default_PT(+Inf), length(d)))

const StandardUniformDensity = Union{
    DistributionDensity{<:StandardUvUniform},
    DistributionDensity{<:StandardMvUniform}
}

const StandardNormalDensity= Union{
    DistributionDensity{<:StandardUvNormal},
    DistributionDensity{<:StandardMvNormal}
}


text:
# Type
"""
	(type) Arc{T<:AnyComplex} in the complex plane

Each `Arc` type is parameterized according to the common type of its complex input arguments.
"""
struct Arc{T<:AnyComplex} <: AbstractCurve
	circle::Circle{T}
	start::Float64  # specified as positive fraction of 1 ccw rotation from positive real
	delta::Float64
end

# Untyped constructors
"""
	Arc(C,start,delta)

Consruct an arc that is the part of the Circle `C` starting at parameter value `start` and ending at `start+delta`. The values are expressed in terms of fractions of a complete circle. The `start` value should be in [0,1), and `delta` should be in [-1,1].

	Arc(a,b,c)

Construct the arc starting at point `a`, passing through `b`, and ending at `c`. If the three points are collinear, a `Segment` is returned instead.
"""
function Arc(C::Circle{T},start::Real,delta::Real) where T<:AnyComplex
	if delta < 0
		Arc{T}(reverse(C),float(-start),float(-delta))
	else
		Arc{T}(C,float(start),float(delta))
	end
end

# Construct from 3 points
function Arc(a::Number,m::Number,b::Number)
	a,m,b = promote(complex(float(a)),m,b)
	C = Circle(a,m,b)
	if isa(C,Line)  # collinear
		Segment(a,b)
	else
		α,β = a-C.center,b-C.center
		ti = mod(angle(α)/(2π),1)
		delta = mod(angle(β/α)/(2π),1) # force into (0,1)
		# which of the two circle pieces do we use?
		if !C.ccw
			ti = 1-ti
			delta = 1-delta
		end
		Arc(C,ti,delta)
	end
end

# Complex type converters
for ctype in [:Spherical,:Polar,:Complex]
	@eval begin
		function $ctype(A::Arc{T}) where T<:AnyComplex
			Arc($ctype(A.circle),A.start,A.delta)
		end
	end
end

# Required methods
function point(A::Arc,t::Real)
	s = scaleto(A.start,A.start+A.delta,t)
	point(A.circle,s)
end
arclength(A::Arc) = arclength(A.circle)*A.delta
(C::Arc)(t::Real) = point(C,t)

"""
	arg(A::Arc,z)

Find the parameter argument `t` such that `A(t)==z` is true.

This gives undefined results if `z` is not actually on the arc.
"""
function arg(A::Arc,z::Number)
	tc = arg(A.circle,z)
	t = mod(tc-A.start+DEFAULT[:tol],1)-DEFAULT[:tol]
	t/A.delta
end

unittangent(A::Arc,t::Real) = unittangent(A.circle,A.start+t*A.delta)
function tangent(A::Arc{T},t::Real) where T <: AnyComplex
	T( tangent(A.circle,A.start+t*A.delta)*A.delta )
end

+(A::Arc,z::Number) = Arc(A.circle+z,A.start,A.delta)
function -(A::Arc)
	ti = mod(A.start+0.5,1)
	Arc(-A.circle,ti,A.delta)
end

function *(A::Arc,z::Number)
	phi = angle(z)/(2*pi)
	ti = A.circle.ccw ? mod(A.start+phi,1) : mod(A.start-phi,1)
	Arc(z*A.circle,ti,A.delta)
end

"""
	inv(A)

Invert the arc `A` through the origin. In general the inverse is an `Arc`, though the result is a `Segment` if the arc's circle passes through the origin.
"""
function inv(A::Arc)
	w = 1 ./ point(A,[0,0.5,1])
	Arc(w...)
end

# Other methods
isfinite(::Arc) = true
conj(A::Arc) = Arc(conj(A(0)),conj(A(0.5)),conj(A(1)))
reverse(A::Arc) = Arc(A(1),A(0.5),A(0))

"""
	isapprox(A1::Arc,A2::Arc; tol=<default>)
	A1 ≈ A2

Determine if `A1` and `A2` represent the same arc, irrespective of the type or values of its parameters. Identity is determined by agreement within `tol`, which is interpreted as the weaker of absolute and relative differences.
"""
function isapprox(A1::Arc,A2::Arc;tol=DEFAULT[:tol])
	if isapprox(A1.circle,A2.circle,tol=tol)
		z1 = point(A1,[0,1])
		z2 = point(A2,[0,1])
		appx(u,v) = isapprox(u,v,atol=tol,rtol=tol)
		return appx(z1,z2) || appx(z1,reverse(z2))
	else
		return false
	end
end

"""
	dist(z,A::Arc)

Compute the distance from number `z` to the arc `A`.
"""
function dist(z::Number,A::Arc)
	C = A.circle
	# result depends on whether the radius to z intersects the arc
	ti,del = A.start,A.delta
	if mod(arg(C,z)-ti,1) ≤ del
		return dist(z,C)
	else
		return min( abs(z-point(A,0)), abs(z-point(A,1)) )
	end
end

"""
	closest(z,A::Arc)

Find the point on arc `A` that lies closest to `z`.
"""
function closest(z::Number,A::Arc)
	ζ = z - A.circle.center
	d = A.delta/2
	# rotate arc to position symmetric about positive Re axis
	ϕ = angle( ζ/exp(2im*pi*(A.start+d)) ) / (2π)
	if ϕ > d
		point(A,1)
	elseif ϕ < -d
		point(A,0)
	else
		A.circle.center + A.circle.radius*sign(ζ)
	end
end

function show(io::IO,A::Arc{T}) where {T}
	print(IOContext(io,:compact=>true),"Arc(",A(0.0),"...",A(1.0),")")
end
function show(io::IO,::MIME"text/plain",A::Arc{T}) where {T}
	print(io,"Arc{$T} in the complex plane:\n   fraction ",A.delta," of (",A.circle,") starting at ",A.start)
end


text:
using SimpleTasks
using Base.Test

include("SimpleTasksTests.jl")

max_forks_repo_name: seung-lab/SimpleTasks


text:
function test_fba(model, obj_ider, cost_ider = nothing; summary = false)
    fbaout = isnothing(cost_ider) ? ChLP.fba(model, obj_ider) : ChLP.fba(model, obj_ider, cost_ider)
    ChU.tagprintln_inmw("FBA SOLUTION",
        "\nobj_ider:         ", obj_ider,
        "\nfba obj_val:      ", ChU.av(model, fbaout, obj_ider),
        "\nmax exp obj_val:  ", maximum(FolsomData.val("D")),
        "\ncost_ider:        ", isnothing(cost_ider) ? "Not included" : cost_ider,
        "\nfba cost_val:     ", isnothing(cost_ider) ? "Not included" : ChU.av(model, fbaout, cost_ider),
    )
    summary && ChU.summary(model, fbaout)
    nothing
end

function test_fba(exp::Int, model, obj_ider, cost_ider = nothing; summary = false)
    fbaout = isnothing(cost_ider) ? ChLP.fba(model, obj_ider) : ChLP.fba(model, obj_ider, cost_ider)
    ChU.tagprintln_inmw("FBA SOLUTION",
        "\nobj_ider:         ", obj_ider,
        "\nfba obj_val:      ", ChU.av(model, fbaout, obj_ider),
        "\nexp obj_val:      ", FolsomData.val("D", exp),
        "\ncost_ider:        ", isnothing(cost_ider) ? "Not included" : cost_ider,
        "\nfba cost_val:     ", isnothing(cost_ider) ? "Not included" : ChU.av(model, fbaout, cost_ider),
    )
    summary && ChU.summary(model, fbaout)
    nothing
end


text:
"""
    std_res_differential!(gc)
compute the gradient of residual vector ∇resβ (standardized residual) with respect to beta. For Normal it will be X.
"""
function std_res_differential!(gc::Union{GLMCopulaVCObs{T, D, Link}, GLMCopulaARObs{T, D, Link}, GLMCopulaCSObs{T, D, Link}}) where {T <: BlasReal, D<:Normal{T}, Link}
    gc.∇resβ .= gc.X
    gc.∇resβ .*= -one(T)
    nothing
end

"""
    std_res_differential!(gc)
compute the gradient of residual vector ∇resβ (standardized residual) with respect to beta, for Negative Binomial.
"""
function std_res_differential!(gc::Union{NBCopulaVCObs{T, D, Link}, NBCopulaARObs{T, D, Link}, NBCopulaCSObs{T, D, Link}}
    ) where {T<: BlasReal, D<:NegativeBinomial{T}, Link}
    @turbo for i in 1:gc.p
         for j in 1:gc.n
            gc.∇resβ[j, i] = -inv(sqrt(gc.varμ[j])) * gc.dμ[j] * gc.X[j, i] - (0.5 * inv(gc.varμ[j])) * gc.res[j] * (gc.μ[j] * inv(gc.d.r) + (1 + inv(gc.d.r) * gc.μ[j])) *  gc.dμ[j] * gc.X[j, i]
        end
    end
    nothing
end

"""
    std_res_differential!(gc)
compute the gradient of residual vector ∇resβ (standardized residual) with respect to beta, for Poisson.
"""
function std_res_differential!(gc::Poisson_Bernoulli_VCObs{T, VD, VL}) where {T<: BlasReal, VD, VL}
    p1 = Integer(gc.p/2)
    @inbounds for i in 1:p1
            # first is poisson
            gc.∇resβ[1, i] = gc.X[1, i]
            gc.∇resβ[1, i] *= -(inv(sqrt(gc.varμ[1])) + (0.5 * inv(gc.varμ[1])) * gc.res[1]) * gc.dμ[1]
        end
    @inbounds for i in (p1 + 1):gc.p
            # second is bernoulli
            gc.∇resβ[2, i] = -sqrt(gc.varμ[2]) * gc.X[2, i] - (0.5 * gc.res[2] * (1 - (2 * gc.μ[2])) * gc.X[2, i])
        end
end

"""
    std_res_differential!(gc)
compute the gradient of residual vector ∇resβ (standardized residual) with respect to beta, for Poisson.
"""
function std_res_differential!(gc::Union{GLMCopulaVCObs{T, D, Link}, GLMCopulaARObs{T, D, Link}, GLMCopulaCSObs{T, D, Link}}) where {T<: BlasReal, D<:Poisson{T}, Link}
    @turbo for i in 1:gc.p
        for j in 1:gc.n
            gc.∇resβ[j, i] = gc.X[j, i]
            gc.∇resβ[j, i] *= -(inv(sqrt(gc.varμ[j])) + (0.5 * inv(gc.varμ[j])) * gc.res[j]) * gc.dμ[j]
        end
    end
    nothing
end

"""
    std_res_differential!(gc)
compute the gradient of residual vector ∇resβ (standardized residual) with respect to beta, for Bernoulli.
"""
function std_res_differential!(gc::Union{GLMCopulaVCObs{T, D, Link}, GLMCopulaARObs{T, D, Link}, GLMCopulaCSObs{T, D, Link}}) where {T<: BlasReal, D<:Bernoulli{T}, Link}
    @turbo for i in 1:gc.p
         for j in 1:gc.n
            gc.∇resβ[j, i] = -sqrt(gc.varμ[j]) * gc.X[j, i] - (0.5 * gc.res[j] * (1 - (2 * gc.μ[j])) * gc.X[j, i])
        end
    end
    nothing
end

"""
    glm_gradient(gc, β)
Calculates the gradient with respect to beta for our the glm portion for one obs. Keeps the residuals standardized.
"""
function glm_gradient(gc::Union{GLMCopulaVCObs, NBCopulaVCObs, GLMCopulaARObs, NBCopulaARObs, NBCopulaCSObs, Poisson_Bernoulli_VCObs, GLMCopulaCSObs})
    gc.storage_n .= gc.w1 .* gc.res
    mul!(gc.storage_p1, transpose(gc.X), gc.storage_n)
    gc.storage_p1
end

"""
    glm_hessian(gc)

Compute the part of the hessian relevant to the glm density with respect to beta for a single obs
"""
function glm_hessian(gc::Union{GLMCopulaVCObs, NBCopulaVCObs, GLMCopulaARObs, NBCopulaARObs, NBCopulaCSObs, Poisson_Bernoulli_VCObs, GLMCopulaCSObs})
    mul!(gc.storage_np, Diagonal(gc.w2), gc.X)
    BLAS.gemm!('T', 'N', -1.0, gc.X, gc.storage_np, 0.0, gc.storage_pp)
end

"""
    update_∇θ!(gcm)

Update θ gradient for Newton's Algorithm, given β.
"""
function update_∇θ!(
    gcm::Union{GLMCopulaVCModel{T, D, Link}, NBCopulaVCModel{T, D, Link}}) where {T <: BlasReal, D, Link}
    @inbounds for i in eachindex(gcm.data)
        standardize_res!(gcm.data[i])            # standardize the residuals GLM variance(μ)
        QuasiCopula.update_quadform!(gcm.data[i]) # with standardized residuals
        gcm.QF[i, :] .= gcm.data[i].q
    end
    BLAS.gemv!('N', T(1), gcm.QF,  gcm.θ, T(0), gcm.storage_n)
    gcm.storage_n .= inv.(1 .+ gcm.storage_n)
    BLAS.gemv!('N', T(1), gcm.TR,  gcm.θ, T(0), gcm.storage_n2)
    gcm.storage_n2 .= inv.(1 .+ gcm.storage_n2)
    BLAS.gemv!('T', T(1), gcm.QF, gcm.storage_n, T(0), gcm.∇θ)
    BLAS.gemv!('T', -T(1), gcm.TR, gcm.storage_n2, T(1), gcm.∇θ)
end

"""
    update_Hθ!(gcm)

Update θ Hessian for Newton's Algorithm, given β.
"""
function update_Hθ!(
    gcm::Union{GLMCopulaVCModel{T, D, Link}, NBCopulaVCModel{T, D, Link}}) where {T <: BlasReal, D, Link}
    fill!(gcm.Hθ, 0.0)
    @inbounds for j in 1:gcm.m
        @simd for i in 1:gcm.n
            gcm.hess1[j, i] = gcm.QF[i, j] * gcm.storage_n[i]
            gcm.hess2[j, i] = gcm.TR[i, j] * gcm.storage_n2[i]
        end
    end
    BLAS.gemm!('N', 'T', -T(1), gcm.hess1, gcm.hess1, T(0), gcm.Hθ)
    BLAS.gemm!('N', 'T', T(1), gcm.hess2, gcm.hess2, T(1), gcm.Hθ)
end


text:
module Client

using HTTP, JSON3, Base64
using ..Model, ..Auth

const SERVER = Ref{String}("localhost")
const AUTH_TOKEN = Ref{String}()

function createUser(username, password)
    body = (; username, password=base64encode(password))
    resp = HTTP.post(string(SERVER[], "/user"), [], JSON3.write(body); cookies=true)
    if HTTP.hasheader(resp, Auth.JWT_TOKEN_COOKIE_NAME)
        AUTH_TOKEN[] = HTTP.header(resp, Auth.JWT_TOKEN_COOKIE_NAME)
    end
    return JSON3.read(resp.body, User)
end

function loginUser(username, password)
    body = (; username, password=base64encode(password))
    resp = HTTP.post(string(SERVER[], "/user/login"), [], JSON3.write(body); cookies=true)
    if HTTP.hasheader(resp, Auth.JWT_TOKEN_COOKIE_NAME)
        AUTH_TOKEN[] = HTTP.header(resp, Auth.JWT_TOKEN_COOKIE_NAME)
    end
    return JSON3.read(resp.body, User)
end

function createAlbum(name, artist, year, songs)
    body = (; name, artist, year, songs)
    resp = HTTP.post(string(SERVER[], "/album"), [Auth.JWT_TOKEN_COOKIE_NAME => AUTH_TOKEN[]], JSON3.write(body); cookies=true)
    return JSON3.read(resp.body, Album)
end

function getAlbum(id)
    resp = HTTP.get(string(SERVER[], "/album/$id"), [Auth.JWT_TOKEN_COOKIE_NAME => AUTH_TOKEN[]]; cookies=true)
    return JSON3.read(resp.body, Album)
end

function updateAlbum(album)
    resp = HTTP.put(string(SERVER[], "/album/$(album.id)"), [Auth.JWT_TOKEN_COOKIE_NAME => AUTH_TOKEN[]], JSON3.write(album); cookies=true)
    return JSON3.read(resp.body, Album)
end

function deleteAlbum(id)
    resp = HTTP.delete(string(SERVER[], "/album/$id"), [Auth.JWT_TOKEN_COOKIE_NAME => AUTH_TOKEN[]]; cookies=true)
    return
end

function pickAlbumToListen()
    resp = HTTP.get(string(SERVER[], "/"), [Auth.JWT_TOKEN_COOKIE_NAME => AUTH_TOKEN[]]; cookies=true)
    return JSON3.read(resp.body, Album)
end

end # module

max_stars_repo_name: vanhp/MicroServiceJulia

max_issues_repo_name: vanhp/MicroServiceJulia

max_forks_repo_name: vanhp/MicroServiceJulia


text:
__precompile__()
module CodecLz4
using Compat
using TranscodingStreams: TranscodingStream, Memory, Error
using TranscodingStreams
export LZ4Compressor, LZ4CompressorStream,
    LZ4Decompressor, LZ4DecompressorStream,
    BlockSizeID, default_size, max64KB, max256KB, max1MB, max4MB,
    BlockMode, block_linked, block_independent,
    FrameType, normal_frame, skippable_frame

depsjl = joinpath(@__DIR__, "..", "deps", "deps.jl")
if isfile(depsjl)
    include(depsjl)
else
    error("CodecLz4 not properly installed. Please run Pkg.build(\"CodecLz4\") and restart julia")
end

struct LZ4Exception <: Exception
    src::AbstractString
    msg::AbstractString
end

function Base.showerror(io::IO, ex::LZ4Exception, bt; backtrace=false)
    printstyled(io, "$(ex.src): $(ex.msg)", color=Base.error_color())
end

include("lz4frame.jl")
include("stream_compression.jl")

end


text:
"Holds the coefficients of a projected Gauss-Legendre Runge-Kutta method."
struct CoefficientsPGLRK{T} <: AbstractCoefficients{T}
    @HeaderCoefficientsRK
    @CoefficientsRK

    P::Matrix{T}
    Q::Matrix{T}
    X::Matrix{T}
    W::Matrix{T}
    A::Matrix{T}

    function CoefficientsPGLRK{T}(name,o,s,a,b,c,P,X,W) where {T}
        @assert T <: Real
        @assert isa(name, Symbol)
        @assert isa(o, Integer)
        @assert isa(s, Integer)
        @assert s ≥ 2 "Number of stages must be ≥ 2"
        @assert s==size(a,1)==size(a,2)==length(b)==length(c)
        @assert s==size(P,1)==size(P,2)
        @assert s==size(X,1)==size(X,2)
        @assert s==size(W,1)==size(W,2)

        Q = inv(P)
        A = zero(a)
        B = zero(a)

        simd_mult!(B, W, Q)
        simd_mult!(A, P, B)

        new(name,o,s,a,b,c,P,Q,X,W,A)
    end
end

function CoefficientsPGLRK(name::Symbol, order::Int, a::Matrix{T}, b::Vector{T}, c::Vector{T}, P::Matrix{T}, X::Matrix{T}, W::Matrix{T}) where {T}
    CoefficientsPGLRK{T}(name, order, length(c), a, b, c, P, X, W)
end

Base.hash(tab::CoefficientsPGLRK, h::UInt) = hash(tab.o, hash(tab.s, hash(tab.a, hash(tab.b, hash(tab.c, hash(tab.P, hash(tab.Q, hash(tab.X, hash(tab.W, hash(tab.A, h))))))))))

Base.:(==)(tab1::CoefficientsPGLRK, tab2::CoefficientsPGLRK) = (tab1.o == tab2.o
                                                             && tab1.s == tab2.s
                                                             && tab1.a == tab2.a
                                                             && tab1.b == tab2.b
                                                             && tab1.c == tab2.c
                                                             && tab1.P == tab2.P
                                                             && tab1.Q == tab2.Q
                                                             && tab1.X == tab2.X
                                                             && tab1.W == tab2.W
                                                             && tab1.A == tab2.A)

Base.isequal(tab1::CoefficientsPGLRK{T1}, tab2::CoefficientsPGLRK{T2}) where {T1,T2} = (tab1 == tab2 && T1 == T2 && typeof(tab1) == typeof(tab2))

"Print Runge-Kutta coefficients."
function Base.show(io::IO, tab::CoefficientsPGLRK)
    print(io, "Projected Gauss-Legendre Runge-Kutta Coefficients ", tab.name, " with ", tab.s, " stages and order ", tab.o)
    print(io, "  a = ", tab.a)
    print(io, "  b = ", tab.b)
    print(io, "  c = ", tab.c)
    print(io, "  P = ", tab.P)
    print(io, "  Q = ", tab.Q)
    print(io, "  X = ", tab.X)
    print(io, "  W = ", tab.W)
    print(io, "  A = ", tab.A)
end

function getTableauPGLRK(coeff::CoefficientsPGLRK{T}, λ::T, a::Matrix{T}) where {T}
    a .= coeff.a .+ λ .* coeff.A
end

"Parameters for right-hand side function of variational partitioned Runge-Kutta methods."
mutable struct ParametersPGLRK{DT,TT,D,S,ET} <: Parameters{DT,TT}
    equ::ET
    tab::CoefficientsPGLRK{TT}
    Δt::TT

    A_q::Array{TT, 3}
    A_p::Array{TT, 3}
    a_q::Matrix{TT}
    a_p::Matrix{TT}

    t̅::TT

    q̅::Vector{DT}
    p̅::Vector{DT}

    function ParametersPGLRK{DT,TT,D,S,ET}(equ, tab, Δt) where {DT,TT,D,S,ET}
        # create coefficient matrices
        A_q = zeros(TT, S, S, D)
        A_p = zeros(TT, S, S, D)
        a_q = zero(tab.a)
        a_p = zero(tab.a)

        new(equ, tab, Δt, A_q, A_p, a_q, a_p, zero(TT), zeros(DT,D), zeros(DT,D))
    end
end

mutable struct IntegratorCachePGLRK{DT,D,S} <: IODEIntegratorCache{DT,D}
    θ::Vector{DT}
    θ̅::Vector{DT}
    λ::Vector{DT}
    λ̅::Vector{DT}

    q̃::Vector{DT}
    p̃::Vector{DT}
    ṽ::Vector{DT}
    f̃::Vector{DT}
    θ̃::Vector{DT}
    λ̃::Vector{DT}
    s̃::Vector{DT}

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

    function IntegratorCachePGLRK{DT,D,S}() where {DT,D,S}
        # create solution vectors
        θ = zeros(DT,D)
        θ̅ = zeros(DT,D)
        λ = zeros(DT,D)
        λ̅ = zeros(DT,D)

        # create temporary vectors
        q̃ = zeros(DT,D)
        p̃ = zeros(DT,D)
        ṽ = zeros(DT,D)
        f̃ = zeros(DT,D)
        θ̃ = zeros(DT,D)
        λ̃ = zeros(DT,D)
        s̃ = zeros(DT,D)

        # create internal stage vectors
        Q = create_internal_stage_vector(DT, D, S)
        P = create_internal_stage_vector(DT, D, S)
        V = create_internal_stage_vector(DT, D, S)
        F = create_internal_stage_vector(DT, D, S)
        Y = create_internal_stage_vector(DT, D, S)
        Z = create_internal_stage_vector(DT, D, S)

        new(θ, θ̅, λ, λ̅,
            q̃, p̃, ṽ, f̃, θ̃, λ̃, s̃,
            Q, P, V, F, Y, Z)
    end
end

"Variational partitioned Runge-Kutta integrator."
struct IntegratorPGLRK{DT, TT, PT <: ParametersPGLRK{DT,TT},
                               ST <: NonlinearSolver{DT},
                               IT <: InitialGuessPODE{DT,TT}, N, D, S} <: IntegratorPRK{DT,TT}
    params::PT
    solver::ST
    iguess::IT
    cache::IntegratorCachePGLRK{DT,D,S}
end

function IntegratorPGLRK(equation::IODE{DT,TT,ΑT,FT,GT,VT,N}, tableau::CoefficientsPGLRK{TT}, Δt::TT) where {DT,TT,ΑT,FT,GT,VT,N}
    D = equation.d
    M = equation.n
    S = tableau.s

    # create params
    params = ParametersPGLRK{DT,TT,D,S,typeof(equation)}(equation, tableau, Δt)

    # create solver
    solver = create_nonlinear_solver(DT, D*(S+3), params)

    # create initial guess
    iguess = InitialGuessPODE(get_config(:ig_interpolation), equation, Δt)

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

    # create integrator
    IntegratorPGLRK{DT, TT, typeof(params), typeof(solver), typeof(iguess), N, D, S}(params, solver, iguess, cache)
end

@inline nstages(integrator::IntegratorPGLRK{DT,TT,PT,ST,IT,N,D,S}) where {DT,TT,PT,ST,IT,N,D,S} = S

function update_params!(params::ParametersPGLRK, sol::AtomicSolutionPODE)
    # set time for nonlinear solver and copy previous solution
    params.t̅  = sol.t
    params.q̅ .= sol.q
    params.p̅ .= sol.p
end

function compute_stages!(x::Vector{ST}, cache::IntegratorCachePGLRK{ST}, params::ParametersPGLRK) where {ST}
    compute_stages!(x, cache.Q, cache.V,
                       cache.P, cache.F,
                       cache.Y, cache.Z,
                       cache.q̃, cache.p̃,
                       cache.θ̃, cache.λ̃,
                       cache.ṽ, cache.f̃,
                       params)
end

function compute_stages!(x::Vector{ST}, Q::Vector{Vector{ST}}, V::Vector{Vector{ST}},
                                        P::Vector{Vector{ST}}, F::Vector{Vector{ST}},
                                        Y::Vector{Vector{ST}}, Z::Vector{Vector{ST}},
                                        q::Vector{ST}, p::Vector{ST},
                                        θ::Vector{ST}, λ::Vector{ST},
                                        y::Vector{ST}, z::Vector{ST},
                                        params::ParametersPGLRK{DT,TT,D,S}) where {ST,DT,TT,D,S}

    local tᵢ::TT
    local t₀::TT = params.t̅
    local t₁::TT = params.t̅ + params.Δt

    # copy x to V
    for i in 1:S
        for k in 1:D
            V[i][k] = x[D*(i-1)+k]
        end
    end

    # copy x to λ, q and p
    for k in 1:D
        q[k] = x[D*(S+0)+k]
        p[k] = x[D*(S+1)+k]
        λ[k] = x[D*(S+2)+k]
    end

    # compute tableaus
    for k in 1:D
        getTableauPGLRK(params.tab, λ[k], params.a_q)
        # get_symplectic_conjugate_coefficients(params.a_q, params.tab.b, params.a_p)
        for j=1:S
            for i=1:S
                params.A_q[i,j,k] = params.a_q[i,j]
                params.A_p[i,j,k] = params.a_q[i,j]
            end
        end
    end

    # compute Y and Q
    for i in 1:S
        for k in 1:D
            Y[i][k] = 0
            for j in 1:S
                Y[i][k] += params.A_q[i,j,k] * V[j][k]
            end
            Q[i][k] = params.q̅[k] + params.Δt * Y[i][k]
        end
    end

    # compute P=ϑ(Q,V) and F=f(Q,V)
    for i in 1:S
        tᵢ = params.t̅ + params.Δt * params.tab.c[i]
        params.equ.ϑ(tᵢ, Q[i], V[i], P[i])
        params.equ.f(tᵢ, Q[i], V[i], F[i])
    end

    # compute Z
    for i in 1:S
        for k in 1:D
            Z[i][k] = 0
            for j in 1:S
                Z[i][k] += params.A_p[i,j,k] * F[j][k]
            end
        end
    end

    # compute y and z
    y .= 0
    z .= 0
    for k in 1:D
        for j in 1:S
            y[k] += params.tab.b[j] * V[j][k]
            z[k] += params.tab.b[j] * F[j][k]
        end
    end

    # compute θ=ϑ(q,λ)
    params.equ.ϑ(t₁, q, λ, θ)
end

"Compute stages of variational partitioned Runge-Kutta methods."
@generated function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersPGLRK{DT,TT,D,S}) where {ST,DT,TT,D,S}

    cache = IntegratorCachePGLRK{ST,D,S}()

    quote
        compute_stages!(x, $cache, params)

        # compute b = [P-p-AF]
        for i in 1:S
            for k in 1:D
                b[D*(i-1)+k] = $cache.P[i][k] - params.p̅[k] - params.Δt * $cache.Z[i][k]
            end
        end

        # compute b = [q̅-q-bV, p̅-p-bF, p̅-α(q̅)]
        for k in 1:D
            b[D*(S+0)+k] = $cache.q̃[k] - params.q̅[k] - params.Δt * $cache.ṽ[k]
            b[D*(S+1)+k] = $cache.p̃[k] - params.p̅[k] - params.Δt * $cache.f̃[k]
            b[D*(S+2)+k] = $cache.p̃[k] - $cache.θ̃[k]
        end
    end
end

function initialize!(int::IntegratorPGLRK, sol::AtomicSolutionPODE)
    sol.t̅ = sol.t - timestep(int)

    equation(int).v(sol.t, sol.q, sol.p, sol.v)
    equation(int).f(sol.t, sol.q, sol.p, sol.f)

    initialize!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f,
                            sol.t̅, sol.q̅, sol.p̅, sol.v̅, sol.f̅)
end

function initial_guess!(int::IntegratorPGLRK{DT,TT}, sol::AtomicSolutionPODE{DT,TT}) where {DT,TT}
    for i in eachstage(int)
        evaluate!(int.iguess, sol.q, sol.p, sol.v, sol.f,
                              sol.q̅, sol.p̅, sol.v̅, sol.f̅,
                              int.cache.q̃, int.cache.ṽ,
                              tableau(int).c[i])

        for k in eachdim(int)
            int.solver.x[ndims(int)*(i-1)+k] = int.cache.ṽ[k]
        end
    end

    evaluate!(int.iguess, sol.q, sol.p, sol.v, sol.f,
                          sol.q̅, sol.p̅, sol.v̅, sol.f̅,
                          int.cache.q̃, int.cache.p̃,
                          one(TT), one(TT))

    for k in eachdim(int)
        int.solver.x[ndims(int)*(nstages(int)+0)+k] = int.cache.q̃[k]
    end
    for k in eachdim(int)
        int.solver.x[ndims(int)*(nstages(int)+1)+k] = int.cache.p̃[k]
    end
    for k in eachdim(int)
        int.solver.x[ndims(int)*(nstages(int)+2)+k] = 0
    end
end

"Integrate PODE with variational partitioned Runge-Kutta integrator."
function integrate_step!(int::IntegratorPGLRK{DT,TT}, sol::AtomicSolutionPODE{DT,TT}) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int.params, sol)

    # compute initial guess
    initial_guess!(int, sol)

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

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector fields at internal stages
    compute_stages!(int.solver.x, int.cache, int.params)

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

    # copy solution to initial guess
    update!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f)
end