text:
n1 = 1. # bimolecular 1
n2 = 2. # bimolecular 2
n3 = 3. # trimolecular
function cooperativity_index(n1,n2,n3)
    Σ_bi = n1 + n2
    Σ_tri = n3
    Δn = Σ_tri - Σ_bi
    # cooperativity index formulav (relative change)
    Coop_indx = Δn / Σ_bi
end

cooperativity_index(n1,n2,n3)

n1 = [1. ± 0.1, 1. ± 0.1]  # bimolecular 1
n2 = [2. ± 0.1, 2. ± 0.1]# bimolecular 2
n3 = [3. ± 0.1, 4. ± 0.1]

cooperativity_index.(n1,n2,n3)

function solve_with_sds(rates_sds)
    prob = ODEProblem(f,u₀,tspan,10 .^ rates_sds)
    # prob = remake(prob, p=10 .^ rates)
    sol = solve(prob,Vern7(),saveat=0.01, abs_tol=1e-10,rel_tol=1e-10)
    n_tot = reduce(vcat,sum(sol,dims=1))
    return n_tot
end
function solve_model(rates_sds)
    prob = ODEProblem(f,u₀,tspan,10 .^ rates_sds)
    # prob = remake(prob, p=10 .^ rates)
    sol = solve(prob,Vern7(),saveat=0.01, abs_tol=1e-10,rel_tol=1e-10)
    return sol
end

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

densities_tcr = [25.0, 38.0, 1.0] # TCR, pMHC, Ac
densities_tc1 = [15.0, 38.0, 10.0, 1.0] # TCR, pMHC, CD4, Ac
densities_tc2 = [8.0,  38.0, 17.0, 1.0]
densities_tc3 = [3.0,  38.0, 22.0, 1.0]

# rates_T = [-3.37638358011457 ± 0.0483562400944489
#             -0.302013187183543 ± 0.0607189115063286]
# rates_TC1 = [-1.72689171836867 ± 0.469108114325548
#             -0.646695212785203 ± 0.685834531823618]
# rates_TC2 = [-0.644293816490510 ± 0.408735939559103
#             0.617796010276307 ± 0.467542299161307]
# rates_TC3 = [0.891017383250648 ± 0.559684851912230
#             2.16934160248305 ± 0.568940078989975]

rates_T = [-3.37638358011457
            -0.302013187183543]
rates_TC1 = [-1.72689171836867
            -0.646695212785203]
rates_TC2 = [-0.644293816490510
            0.617796010276307]
rates_TC3 = [0.891017383250648
            2.16934160248305]

sd_data = CSV.read("C:/Users/stravaglino3/Downloads/julia_code/n_stds_for_julia.csv")
nT_sds = sd_data[:,1]
nTC1_sds = Array(sd_data[:,2:3])
nTC2_sds = Array(sd_data[:,4:5])
nTC3_sds = Array(sd_data[:,6:7])

tspan = (0.0,16.0)
t = collect(0.01:0.01:16)

gr()
# bimolecular TCR-pMHC bonds
u₀ = [0.0]
f = (du,u,p,t) -> bimolecular!(du,u,p,t,densities_tcr,[])
n_tcr = Array(solve_model(rates_T))'[2:end]
n_tcr_sd = n_tcr .± nT_sds

# TCR-pMHC-CD4 (1)
u₀ = [0.0, 0.0]
f = (du,u,p,t) -> trimolecular!(du,u,p,t,densities_tc1,10 .^ rates_tcr)
n_tc1 = Array(solve_model(rates_TC1))'[2:end,:]
n_tc1_sd = sum(n_tc1 .± nTC1_sds, dims=2)

# TCR-pMHC-CD4 (2)
f = (du,u,p,t) -> trimolecular!(du,u,p,t,densities_tc2,10 .^ rates_tcr)
n_tc2 = Array(solve_model(rates_TC2))'[2:end,:]
n_tc2_sd = sum(n_tc2 .± nTC2_sds, dims=2)
plot!(t,n_tc2_sd, errorbars=:ribbon, alpha=.1,c=:red)

# TCR-pMHC-CD4 (3)
f = (du,u,p,t) -> trimolecular!(du,u,p,t,densities_tc3,10 .^ rates_tcr)
n_tc3 = Array(solve_model(rates_TC3))'[2:end,:]
n_tc3_sd = sum(n_tc3 .± nTC3_sds, dims=2)
plot!(t,n_tc3_sd, errorbars=:ribbon, alpha=0.1,ribboncolor=:cyan)

plot(t, n_tcr_sd, errorbars=:ribbon, alpha=0.1,msc=:black)
plot!(t, n_tcr, lw=3, lc=:black)
plot!(t,n_tc1_sd, errorbars=:ribbon, alpha=0.1, msc=:blue)
plot!(t, sum(n_tc1,dims=2), lw=3, lc=:blue)
plot!(t,n_tc2_sd, errorbars=:ribbon, alpha=0.1, msc=:red)
plot!(t, sum(n_tc2,dims=2), lw=3, lc=:red)
plot!(t,n_tc3_sd, errorbars=:ribbon, alpha=0.1, msc=:cyan)
plot!(t, sum(n_tc3,dims=2), lw=3, lc=:cyan)

scatter!(n_se_tcr[:,1],n_se_tcr[:,2],yerror=n_se_tcr[:,3],markercolor=:black,markersize=5)
scatter!(n_se_1[:,1],n_se_1[:,2],yerror=n_se_1[:,3],markercolor=:blue,markersize=5)
scatter!(n_se_2[:,1],n_se_2[:,2],yerror=n_se_2[:,3],markercolor=:red,markersize=5)
h = scatter!(n_se_3[:,1],n_se_3[:,2],yerror=n_se_3[:,3],markercolor=:cyan,markersize=5,legend=false)

xlabel!("Contact Time tc")
ylabel!("<n>")
# savefig("all_n_with_error_bands")

#-------------------------------------------------------------------------------
# FIT CD4 only data

cd4_1 = CSV.read("cd4_data_from_muaz.csv")
t_cd41, Pa_cd41, n_se_cd41 = clean_data(cd4_1)
n_cd41 = Pa2n(Pa_cd41)

f = (du,u,p,t) -> bimolecular!(du,u,p,t,dens,cons) # enclose constants

# u0 = [0.0;0.0]
u₀ = [0.0]          # Initial condition
tspan = (0.0,8.0)  # Simulation interval
pinit = [1e-4,0.5]      # equation parameter. p = [k₁, k₋₁]
densities_cd41 = [1700., 850., 1.] #CD4, pMHC, A
cons = []

rates_cd41,hes_cd41,l_cd41,boot_cd41 = optimization(bimolecular!,n_cd41,densities_cd41',[],model_opt2,t_cd41)

# Setup the ODE problem, then solve
plotly()
t = collect(0.01:0.01:16)
tspan = (0.0,16.0)
f = (du,u,p,t) -> bimolecular!(du,u,p,t,densities_cd41,cons)
prob = ODEProblem(f,u₀,tspan,10 .^ rates_cd41)
sol_cd41 = Array(solve(prob,Vern7(),saveat=t,abstol=1e-10,reltol=1e-10))'
plot(t,sol_cd41,lc=:orange,lw=3)
scatter!(n_se_cd41[:,1],n_se_cd41[:,2],yerror=n_se_cd41[:,3],markercolor=:orange)

densities_tcr = [25.0, 38.0, 1.0] # TCR, pMHC, Ac
densities_tc1 = [15.0, 38.0, 10.0, 1.0] # TCR, pMHC, CD4, Ac
densities_tc2 = [8.0,  38.0, 17.0, 1.0]
densities_tc3 = [3.0,  38.0, 22.0, 1.0]
densities_cd41 = [1700., 850., 1.] #CD4, pMHC, A
cooperativity_index(n1,n2,n3)
#---------------------
#normalize everything to tc2

n_tcr_norm = n_tcr./ 25 .* 8
n_tc1_norm = sum(n_tc1,dims=2) ./15 .* 8 ./10 .*17
n_tc2_norm = sum(n_tc2,dims=2)
n_tc3_norm = sum(n_tc3,dims=2) ./3 .* 8 ./22 .* 17
n_cd41_norm = sol_cd41 ./1700 .*17 ./850 .*38
n_tcr_norm = n_tcr ./25. .*15.
n_tc1_norm = sum(n_tc1,dims=2)
n_tc2_norm = sum(n_tc2,dims=2) ./8. .*15. ./17. .*10.
n_tc3_norm = sum(n_tc3,dims=2) ./3. .*15. ./22. .*10.
n_cd41_norm = sol_cd41 ./1700. .*10. ./850. .*38

plotly()
t = collect(0.01:0.01:16)
plot(t,n_tcr_norm,lw=3,lc=:black)
plot!(t,n_tc1_norm,lw=3,lc=:blue)
plot!(t,n_tc2_norm,lw=3,lc=:red)
plot!(t,n_tc3_norm,lw=3,lc=:cyan)
plot!(t,n_cd41_norm,lw=3,lc=:orange,yscale=:linear,legend=false)
yaxis!("log<n>")
xlabel!("Contact Time tc")
savefig("n_log_scaled")

coop_indx_tc1 = cooperativity_index.(n_tcr_norm,n_cd41_norm,n_tc1_norm)
coop_indx_tc2 = cooperativity_index.(n_tcr_norm,n_cd41_norm,n_tc2_norm)
coop_indx_tc3 = cooperativity_index.(n_tcr_norm,n_cd41_norm,n_tc3_norm)
plot(t,coop_indx_tc1 .* 100,lw=3,lc=:blue)
plot!(t,coop_indx_tc2 .* 100,lw=3,lc=:red)
plot!(t,coop_indx_tc3 .* 100,lw=3,lc=:cyan,legend=false)
xlabel!("Contact Time tc")
ylabel!("Cooperative Index (%)")
# savefig("Cooperative Index.png")


text:
using ControlSystems, RobustAndOptimalControl
using Plots

"""
This is a simple SISO example which was used for debugging the implementation,
as this exact example was use in the lecture notes of the "Principles of Optimal
Control" cours of the MIT OpenCourseWare [1], however, this example is in discrete time

[1] https://ocw.mit.edu/courses/aeronautics-and-astronautics/16-323-principles-of-optimal-control-spring-2008/lecture-notes/lec15.pdf

The example can be set to visualize and save plots using the two variables
  makeplots - true/false (true if plots are to be generated, false for testing)
"""
makeplots = true

# Define the process
ts = 0.005
Gd = bilinearc2d(ss(tf([200], [0.025,1.0025,10.1,1])),ts)

# Sensitivity weight function
M, wB, A = 1.5, 10, 1e-4
WS = tf([1/M, wB],[1, wB*A])

# Output sensitivity weight function
WU = ss(0.1)

# Complementary sensitivity weight function
WT = []

# Create continuous time approximation of the process
Gc = bilineard2c(ss(Gd))

# Form augmented P dynamics in state-space
Pc = hinfpartition(Gc, WS, WU, WT)

# Check that the assumptions are satisfied
flag = hinfassumptions(Pc)

# Synthesize the H-infinity optimal controller
flag, Cc, γ = hinfsynthesize(Pc)

# Extract the transfer functions defining some signals of interest, but do so
# using discrete equivalent of the continuous time objects Pc, Cc and Gc
Pcl, S, CS, T = hinfsignals(
  bilinearc2d(Pc, ts),
  bilinearc2d(Gc, ts),
  bilinearc2d(Cc, ts)
)

# Visualize results
if makeplots
  specificationplot([S, CS, T], [WS, WU, WT], γ)
  specificationplot(Pcl, γ; s_labels=["\$\\sigma(P_{cl}(j\\omega))\$"], w_labels=["\$\\gamma\$"])
end

ext: jl


text:
R_CNFL = Tariff("CNFL",
                "Residential",
                [(0:200, 63.13), (200:300, 96.87), (300:1000, 100.14)],
                [(0:100, 0.0)],
                0.06,
                21.03,
                3.51,
                0.10,
                [0.992832136 0.999770704 1.00427654 0.983489891 1.009385917 1.012387984 0.979508228 0.994153679 0.97686098 1.015651323 1.015252586 1.016430032],
				1893.90);

CI_CNFL = Tariff("CNFL",
                "Commercial - Industrial",
                [(0:3000, 106.68), (3001:10000, 64.22)],
                [(0:8,  80426.32), (8:1000,  10053.29 )],
                0.06,
                21.03,
                3.51,
                0.10,
                [0.997759211, 1.029462815, 0.970880339, 1.004265761, 1.070739906, 1.039263159, 0.993586335, 1.002523373, 0.983852239, 0.950187922, 0.990476922, 0.967002018],
				1893.90);

TMT_CNFL = Tariff("CNFL",
                "Medium Voltage",
                [("punta", 54.42), ("valle", 27.21), ("noche", 19.59)],
                [("punta",  9542.79), ("valle",  6790.00), ("noche",  4310.42)],
                0.06,
                21.03,
                3.51,
                0.10,
                [1.024524936, 0.949492666, 1.025369474, 1.061316875, 1.023024878, 0.996315038, 0.992322444, 0.986442682, 0.947434961, 1.016134986, 1.006473428, 0.971147632],
				1893.90);

R_ICE = Tariff("ICE",
                "Residential",
                [(0:40.99,  3189.60), (41:200.99,  79.74), (201:1000, 143.71)],
                [(0:100, 0.0)],
                0.06,
                21.03,
                4.50,
                0.01,
                [0.992832136 0.999770704 1.00427654 0.983489891 1.009385917 1.012387984 0.979508228 0.994153679 0.97686098 1.015651323 1.015252586 1.016430032],
				3189.60);

CI_ICE = Tariff("ICE",
                "Commercial - Industrial",
                [(0:3000, 120.10), (3001:10000,  71.86)],
                [(0:1000,  11878.17)],
                0.06,
                21.03,
                4.50,
                0.01,
                [0.997759211, 1.029462815, 0.970880339, 1.004265761, 1.070739906, 1.039263159, 0.993586335, 1.002523373, 0.983852239, 0.950187922, 0.990476922, 0.967002018],
				3189.60);

TMT_ICE = Tariff("ICE",
                "Medium Voltage",
                [("punta",  68.90), ("valle",  25.60), ("noche",  15.76)],
                [("punta",  11181.25), ("valle",  7806.89), ("noche",  5000.50)],
                0.06,
                21.03,
                4.50,
                0.01,
                [1.024524936, 0.949492666, 1.025369474, 1.061316875, 1.023024878, 0.996315038, 0.992322444, 0.986442682, 0.947434961, 1.016134986, 1.006473428, 0.971147632],
				3189.60);

SD = Tariff("ICE",
                "Bulk System",
                [("punta",  53.98), ("valle",  44.22), ("noche",  37.54)],
                [("punta",  2863.19), ("valle",  2863.19), ("noche",  2863.19)],
                0.02,
                0.0,
                0.0,
                0.0,
                [0.0],
				3189.60);

max_stars_repo_name: arnavgautam/PVCR

max_issues_repo_name: arnavgautam/PVCR

max_forks_repo_name: arnavgautam/PVCR


text:
export @scope, @autodestroy

macro autodestroy(obj)
    esc(quote
        o = $obj
        push!(autoDestroyPool, o)
        o
    end)
end

macro scope(args...)
    local autoDestroyPool = []
    esc(quote
        pool = $autoDestroyPool
        (autoDestroyPool->$(args[end]))(pool)
        for i = length(pool) : -1 : 1
            pool[i].destroy()
        end
        empty!(pool)
    end)
end


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

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

howell1 = CSV.read(rel_path("..", "data", "Howell1.csv"), delim=';')
df = convert(DataFrame, howell1);

df2 = filter(row -> row[:age] >= 18, df);

mean_weight = mean(df2[:weight])
df2 = hcat(df2, df2[:weight] .- mean_weight)
rename!(df2, :x1 => :weight_c); # Rename our col :x1 => :weight_c
first(df2, 5)

weightsmodel = "
data {
 int < lower = 1 > N; // Sample size
 vector[N] height; // Predictor
 vector[N] weight; // Outcome
}

parameters {
 real alpha; // Intercept
 real beta; // Slope (regression coefficients)
 real < lower = 0 > sigma; // Error SD
}

model {
 height ~ normal(alpha + weight * beta , sigma);
}
";

stanmodel = Stanmodel(name="weights", monitors = ["alpha", "beta", "sigma"],model=weightsmodel,
  output_format=:mcmcchains);

heightsdata = Dict("N" => length(df2[:height]), "height" => df2[:height], "weight" => df2[:weight_c]);

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

chn.value[1:5,:,1]

nvals = [10, 50, 150, 352];

p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 4)
for i in 1:length(nvals)
  N = nvals[i]
  heightsdataN = [
    Dict("N" => N, "height" => df2[1:N, :height], "weight" => df2[1:N, :weight_c])
  ]
  rc, chnN, cnames = stan(stanmodel, heightsdataN, ProjDir, diagnostics=false,
    summary=false, CmdStanDir=CMDSTAN_HOME)

  rws, vars, chns = size(chnN)
  xi = -15.0:0.1:15.0
  alpha_vals = convert(Vector{Float64}, reshape(chnN.value[:, 1, :], (rws*chns)))
  beta_vals = convert(Vector{Float64}, reshape(chnN.value[:, 2, :], (rws*chns)))

  p[i] = scatter(df2[1:N, :weight_c], df2[1:N, :height], leg=false, xlab="weight_c")
  for j in 1:N
    yi = alpha_vals[j] .+ beta_vals[j]*xi
    plot!(p[i], xi, yi, title="N = $N")
  end
end
plot(p..., layout=(2, 2))

rws, vars, chns = size(chn)
mu_at_50 = link(50:10:50, chn, [1, 2], mean_weight);
density(mu_at_50)

mu = link(25:10:75, chn, [1, 2], mean_weight);

q = Vector{Plots.Plot{Plots.GRBackend}}(undef, size(mu, 1))
for i in 1:size(mu, 1)
  q[i] = density(mu[i], ylim=(0.0, 1.5),
    leg=false, title="mu_bar = $(round(mean(mu[i]), digits=1))")
end
plot(q..., layout=(2, 3), ticks=(3))

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


text:
# ---
# title: 1519. Number of Nodes in the Sub-Tree With the Same Label
# id: problem1519
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search, Breadth-first Search
# link: <https://leetcode.com/problems/number-of-nodes-in-the-sub-tree-with-the-same-label/description/>
# hidden: true
# ---
#
# Given a tree (i.e. a connected, undirected graph that has no cycles)
# consisting of `n` nodes numbered from `0` to `n - 1` and exactly `n - 1`
# `edges`. The **root** of the tree is the node `0`, and each node of the tree
# has **a label** which is a lower-case character given in the string `labels`
# (i.e. The node with the number `i` has the label `labels[i]`).
#
# The `edges` array is given on the form `edges[i] = [ai, bi]`, which means
# there is an edge between nodes `ai` and `bi` in the tree.
#
# Return _an array of size`n`_ where `ans[i]` is the number of nodes in the
# subtree of the `ith` node which have the same label as node `i`.
#
# A subtree of a tree `T` is the tree consisting of a node in `T` and all of its
# descendant nodes.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/07/01/q3e1.jpg)
#
#
#
#     Input: n = 7, edges = [[0,1],[0,2],[1,4],[1,5],[2,3],[2,6]], labels = "abaedcd"
#     Output: [2,1,1,1,1,1,1]
#     Explanation: Node 0 has label 'a' and its sub-tree has node 2 with label 'a' as well, thus the answer is 2. Notice that any node is part of its sub-tree.
#     Node 1 has a label 'b'. The sub-tree of node 1 contains nodes 1,4 and 5, as nodes 4 and 5 have different labels than node 1, the answer is just 1 (the node itself).
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/07/01/q3e2.jpg)
#
#
#
#     Input: n = 4, edges = [[0,1],[1,2],[0,3]], labels = "bbbb"
#     Output: [4,2,1,1]
#     Explanation: The sub-tree of node 2 contains only node 2, so the answer is 1.
#     The sub-tree of node 3 contains only node 3, so the answer is 1.
#     The sub-tree of node 1 contains nodes 1 and 2, both have label 'b', thus the answer is 2.
#     The sub-tree of node 0 contains nodes 0, 1, 2 and 3, all with label 'b', thus the answer is 4.
#
#
# **Example 3:**
#
# ![](https://assets.leetcode.com/uploads/2020/07/01/q3e3.jpg)
#
#
#
#     Input: n = 5, edges = [[0,1],[0,2],[1,3],[0,4]], labels = "aabab"
#     Output: [3,2,1,1,1]
#
#
# **Example 4:**
#
#
#
#     Input: n = 6, edges = [[0,1],[0,2],[1,3],[3,4],[4,5]], labels = "cbabaa"
#     Output: [1,2,1,1,2,1]
#
#
# **Example 5:**
#
#
#
#     Input: n = 7, edges = [[0,1],[1,2],[2,3],[3,4],[4,5],[5,6]], labels = "aaabaaa"
#     Output: [6,5,4,1,3,2,1]
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 10^5`
#   * `edges.length == n - 1`
#   * `edges[i].length == 2`
#   * `0 <= ai, bi < n`
#   * `ai != bi`
#   * `labels.length == n`
#   * `labels` is consisting of only of lower-case English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
export WSum

struct WSum{T}
    c::T
end
@functor WSum
(l::WSum)(x::AbstractVector) = sum(x.*l.c)
(l::WSum)(x::AbstractArray) = sum_autobatch(x.*l.c)

WSum(in::Integer) = WSum(Complex{STD_REAL_PREC}, in)
function WSum(T::Type, in::Integer;
              initb = glorot_uniform)
  return WSum(initb(T, in))
end

mutable struct WSumCache{Ta,Tb,Tc}
    σᵢₙ::Ta
    out::Tc
    δℒ::Tb
    valid::Bool
end

cache(l::WSum, arr_T, in_T, in_sz)  =
    WSumCache(similar(l.c, Complex{real(eltype(l.c))}),
              zero(Complex{real(eltype(l.c))}),
              similar(l.c, Complex{real(eltype(l.c))}, 1, length(l.c)),
              false)

function layer_out_type_size(l::WSum, in_T ,in_sz)
    out_T     = promote_type(in_T, eltype(l.c))
    return out_T, (1,)
end

function (l::WSum)(c::WSumCache, x)
    σ = copyto!(c.σᵢₙ, x)

    # TODO the dot product allocates
    c.out = sum(x.*l.c)
    return c.out
end

function backprop(∇, l::WSum, c::WSumCache, δℒ::Number)
    # compute the derivative
    ∇.c .= c.σᵢₙ

    # Backpropagate
    c.δℒ .= δℒ .* transpose(l.c)
    return c.δℒ
end


text:
@testset "Linear Discrete System" begin
    # Create some linear data
    A = [0.9 -0.2; 0.0 0.2]
    y = [[10.; -10.]]
    for i in 1:10
        push!(y, A*y[end])
    end
    X = hcat(y...)
    prob = DiscreteDataDrivenProblem(X, t = 1:11)

    for alg in [DMDPINV(), DMDSVD(), TOTALDMD()]
        # Returns a named
        estimator = solve(prob, alg , operator_only = true)
        @test Matrix(estimator.K) ≈ A
        @test eigvals(estimator.K) ≈ eigvals(A)
        @test estimator.C ≈ diagm(ones(2))
        @test isempty(estimator.B)

        res = solve(prob, alg , operator_only = false)
        m = metrics(res)
        @test all(m[:L₂] ./ size(X, 2) .< 3e-1)
        @test Matrix(result(res)) ≈ Matrix(estimator.K)
    end
end

@testset "Linear Continuous System" begin
    A = [-0.9 0.1; 0.0 -0.2]
    f(u, p, t) = A*u
    u0 = [10.0; -20.0]
    prob = ODEProblem(f, u0, (0.0, 10.0))
    sol = solve(prob, Tsit5(), saveat = 0.001)

    prob = DataDrivenProblem(sol)

    for alg in [DMDPINV(), DMDSVD(), TOTALDMD()]
        estimator = solve(prob, alg , operator_only = true)
        @test isapprox(Matrix(estimator.K), A, atol = 1e-2)
        @test isapprox(eigvals(estimator.K), eigvals(A), atol = 1e-2)
        @test estimator.C ≈ diagm(ones(2))
        @test isempty(estimator.B)
        res = solve(prob, alg , operator_only = false)
        m = metrics(res)
        @test all(m[:L₂] ./ size(X, 2) .< 3e-1)
        @test Matrix(result(res)) ≈ Matrix(estimator.K)
    end
end

@testset "Low Rank Continuous System" begin
    K̃ = -0.5*I + [0 0 -0.2; 0.1 0 -0.1; 0. -0.2 0]
    F = qr(randn(20, 3))
    Q = F.Q[:, 1:3]
    dudt(u, p, t) = K̃*u
    prob = ODEProblem(dudt, [10.0; 0.3; -5.0], (0.0, 10.0))
    sol_ = solve(prob, Tsit5(), saveat = 0.01)

    # True Rank is 3
    X = Q*sol_[:,:] + 1e-3*randn(20, 1001)
    DX = Q*sol_(sol_.t, Val{1})[:,:] + 1e-3*randn(20, 1001)
    ddprob = ContinuousDataDrivenProblem(X, sol_.t, DX = DX)

    for alg in [TOTALDMD(3, DMDPINV()); TOTALDMD(0.01, DMDSVD(3))]
        res = solve(ddprob, alg, digits = 2)
        K = Matrix(result(res))
        m = metrics(res)
        @test Q'*K*Q ≈ K̃ atol = 1e-1
        @test Q*K̃*Q' ≈ K atol = 1e-1
        @test all(m[:L₂] ./ size(X, 2) .< 1e-2)
    end
end

@testset "Big System" begin
    # Creates a big system which would resulting in a segfault otherwise
    X = rand([0, 1], 128, 936);
    T = collect(LinRange(0, 4.367058580858928, 936));
    problem = DiscreteDataDrivenProblem(X, T);
    @test_nowarn res2 = solve(problem, DMDSVD(), operator_only = true)
end


text:
module Microgrids

export simulate,
       NonDispatchables,
       Project, DieselGenerator, Battery, Photovoltaic, WindPower, Microgrid,
       TotalCosts, OperVarsTraj, OperVarsAggr,
       operation, aggregation, dispatch, production,
       annual_costs, economics

include("components.jl")
include("dispatch.jl")
include("production.jl")
include("operation.jl")
include("economics.jl")

function simulate(mg::Microgrid)
    # Run the microgrid operation
    opervarstraj = operation(mg)

    # Aggregate the operation variables
    opervarsaggr = aggregation(mg, opervarstraj)

    # Eval the microgrid costs
    costs = economics(mg, opervarsaggr)

    return (opervarstraj = opervarstraj, opervarsaggr = opervarsaggr, costs = costs)
end

end


text:
using GLMakie
using Oceananigans
using Oceananigans.Units

arch = CPU()
grid = RectilinearGrid(arch, size=128, z=(-5, 5), topology=(Flat, Flat, Bounded))

single_closure = VerticalScalarDiffusivity(VerticallyImplicitTimeDiscretization(), κ=1)

tuple_closure = (VerticalScalarDiffusivity(VerticallyImplicitTimeDiscretization(), κ=1/2),
                 VerticalScalarDiffusivity(VerticallyImplicitTimeDiscretization(), κ=1/2))

closure = tuple_closure

model = HydrostaticFreeSurfaceModel(; grid, closure, tracers=:c, buoyancy=nothing, coriolis=nothing)
set!(model, c = (x, y, z) -> exp(-z^2))
simulation = Simulation(model, Δt=1e-4, stop_iteration=100)

c_snapshots = []
c = model.tracers.c
getc(sim) = push!(c_snapshots, Array(interior(c, 1, 1, :)))
simulation.callbacks[:c] = Callback(getc, IterationInterval(10))

run!(simulation)

fig = Figure()
ax = Axis(fig[1, 1])
Nt = length(c_snapshots)
slider = Slider(fig[2, 1], range=1:Nt, startvalue=1)
n = slider.value
z = znodes(c)
c = @lift c_snapshots[$n]
lines!(ax, c, z)
display(fig)


text:
# ---
# title: 1109. Corporate Flight Bookings
# id: problem1109
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Math
# link: <https://leetcode.com/problems/corporate-flight-bookings/description/>
# hidden: true
# ---
#
# There are `n` flights, and they are labeled from `1` to `n`.
#
# We have a list of flight bookings.  The `i`-th booking `bookings[i] = [i, j,
# k]` means that we booked `k` seats from flights labeled `i` to `j` inclusive.
#
# Return an array `answer` of length `n`, representing the number of seats
# booked on each flight in order of their label.
#
#
#
# **Example 1:**
#
#
#
#     Input: bookings = [[1,2,10],[2,3,20],[2,5,25]], n = 5
#     Output: [10,55,45,25,25]
#
#
#
#
# **Constraints:**
#
#   * `1 <= bookings.length <= 20000`
#   * `1 <= bookings[i][0] <= bookings[i][1] <= n <= 20000`
#   * `1 <= bookings[i][2] <= 10000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using DifferentialEquations
using ModelingToolkit
using DataDrivenDiffEq
using LinearAlgebra, DiffEqSensitivity, Optim
using DiffEqFlux, Flux
using Plots
gr()

# function triNN!(du,u,p,t,dens,cons)
#     # unpack rates and constants
#     nᵣ,nₓ,n₃ = u
#     k₁,k₋₁,k₂,k₋₂ = cons
#     mᵣ,mₗ,mₓ,A = dens
#     z = L(u,p)
#     # model
#     du[1] = dnᵣ = A*k₁*mᵣ*mₗ - k₋₁*nᵣ + z[1]
#     du[2] = dnₓ = A*k₂*mₓ*mₗ - k₋₂*nₓ + z[2]
#     du[3] = dn₃ = z[1] + z[2]
#
# end
function bimolecular!(du,u,p,t,dens,cons)
    # unpack rates and constants
    nᵣ = u[1]
    k₁,k₋₁  = p
    mᵣ,mₗ,A = dens
    # model
    du[1] = dnᵣ = A*k₁*mᵣ*mₗ - k₋₁*nᵣ

end

dens = [100.,10.,1.]
cons = Float64[]
f = (du,u,p,t) -> bimolecular!(du,u,p,t,dens,cons)
# Define the experimental parameter
tspan = (0.0f0,10.0f0)
u0 = Float32[0.]
p_ = Float32[1e-4, 0.5]
prob = ODEProblem(f, u0,tspan, p_)
solution = solve(prob, Tsit5(), saveat = 0.1)

scatter(solution, alpha = 0.25)
plot!(solution, alpha = 0.5)

# Ideal data
X = Array(solution)
# Add noise to the data
println("Generate noisy data")
Xₙ = X + Float32(1e-2)*randn(eltype(X), size(X))

# function dudt_(u, p,t)
#     x = u
#     z = L(u,p)
#     [p_[1]*x + z[1],
#     -p_[4]*y + z[2]]
# end
function biNN!(du,u,p,t,dens,cons)
    # unpack rates and constants
    nᵣ = u[1]
    k₁,k₋₁  = p
    mᵣ,mₗ,A = dens
    z = L(u,p)
    # model
    du[1] = dnᵣ = z[1]

end

L = FastChain(FastDense(1, 32, tanh),FastDense(32, 32, tanh), FastDense(32, 1))
#L = FastChain(FastDense(3, 20, tanh),FastDense(20, 20, tanh), FastDense(20, 2))
p = initial_params(L)

#
fNN = (du,u,p,t) -> biNN!(du,u, p,t,dens,cons)
prob_nn = ODEProblem(fNN,u0, tspan, p)
sol_nn = solve(prob_nn, Tsit5(), saveat =  solution.t)

plot(solution)
plot!(sol_nn)
# summ = reduce(vcat,sum(sol_nn,dims=1))
# h = plot!(sol_nn.t,summ,linecolor=:black)

function predict(θ)

    tmp_prob = remake(prob_nn,u0=u0,p=θ)
    tmp_sol =  solve(tmp_prob, VCABM(), saveat = solution.t,
                  abstol=1e-5, reltol=1e-5)
                  # backsolve=true)
    # Σ_sol = sum(Array(tmp_sol),dims=1) # Note: this returns a row vector!
    Array(tmp_sol)
end

function loss(θ)
    pred = predict(θ)
    sum(abs2, Xₙ .- pred), pred
end

# Test
@time loss(p)

const losses = []

callback(θ,l,pred) = begin
    push!(losses, l)
    @show l
    if length(losses)%50==0
        println("Current loss after $(length(losses)) iterations: $(losses[end])")
    end
    h = plot(solution.t, pred')
    scatter!(h, solution.t, Xₙ')
    display(h)
    false
end

p = initial_params(L)
# First train with ADAM for better convergence
@time res1 = DiffEqFlux.sciml_train(loss, p, ADAM(0.01), cb=callback, maxiters = 40)

# Train with BFGS
res2 = DiffEqFlux.sciml_train(loss, res1.minimizer, BFGS(),
                                cb=callback, maxiters = 100)

# @time res1 = DiffEqFlux.sciml_train(loss, res2.minimizer, ADAM(0.001), cb=callback, maxiters = 20)

println("Final training loss after $(length(losses)) iterations: $(losses[end])")

# Plot the losses
imgloss= plot(losses, yaxis = :log, xaxis = :log, xlabel = "Iterations", ylabel = "Loss",fmt = :svg)
savefig(imgloss,"training_loss")
# Plot the data and the approximation
NNsolution = predict(res2.minimizer)
# Trained on noisy data vs real solution
plot(solution.t, NNsolution',linewidth=3,xlabel = "t (s)", ylabel = "n",fmt = :svg)
scatter!(solution.t, Xₙ')
savefig("fit")

L̂ = L(Xₙ,res2.minimizer)
## Sparse Identification
#
# # Create a Basis
# @variables u[1:2]
# # Lots of polynomials
# polys = Operation[1]
#
# for i ∈ 1:2
#     push!(polys, u[1]^i)
#     push!(polys, u[2]^i)
#     for j ∈ i:2
#         if i != j
#             push!(polys, (u[1]^i)*(u[2]^j))
#             push!(polys, u[2]^i*u[1]^i)
#         end
#     end
# end
#
# # And some other stuff
# h = [cos.(u)...; sin.(u)...; polys...]
# h = [ polys...]
# basis = Basis(h, u)
#
# # Create an optimizer for the SINDy problem
# opt = SR3()
# # Create the thresholds which should be used in the search process
# λ = exp10.(-7:0.1:3)
# # Target function to choose the results from; x = L0 of coefficients and L2-Error of the model
# g(x) = x[1] < 1 ? Inf : norm(x, 2)
# # Test on original data and without further knowledge
# # println("SINDy on full ideal, unavailable data")
# # Ψ = SINDy(Xₙ[:, :], DX[:, :], basis, λ, opt, g = g, maxiter = 10000) # Fail
# # println(Ψ)
# # print_equations(Ψ)
# #
# # # Test on ideal derivative data for unknown function ( not available )
# # println("SINDy on partial ideal, unavailable data")
# # Ψ = SINDy(Xₙ[:, 1:end], L̄[:, 1:end], basis, λ,opt, g = g, maxiter = 10000) # Succeed
# # println(Ψ)
# # print_equations(Ψ)
#
# # Test on uode derivative data
# println("SINDy on learned, partial, available data")
# Ψ = SINDy(Xₙ[:, 2:end], L̂[:, 2:end], basis, λ,  opt, g = g, maxiter = 100000, normalize = true, denoise = true) # Succeed
# println(Ψ)
# print_equations(Ψ)
#
# # Extract the parameter
# p̂ = parameters(Ψ)
# println("First parameter guess : $(p̂)")

@variables u[1:1]
h = Operation[u; u.^2; 1]
basis = Basis(h, u)

# Create the thresholds which should be used in the search process
λ = exp10.(-7:0.1:3)
# Target function to choose the results from; x = L0 of coefficients and L2-Error of the model
g(x) = x[1] < 1 ? Inf : norm(x, 2)

opt = SR3(3e-1, 1.0)
# Ψ = SInDy(X[:, 1:1000], DX[:, 1:1000], basis, maxiter = 10000, opt = opt, normalize = true)
Ψ = SINDy(Xₙ[:, 2:end], L̂[:, 2:end], basis, λ,  opt, g = g, maxiter = 10000, normalize = true) # Succeed
# Ψ = SINDy(Xₙ[:, 2:end], L̂[:, 2:end], basis, maxiter = 10000, opt = opt, normalize = true,denoise = true)
println(Ψ)
print_equations(Ψ)

sys = ODESystem(Ψ)
p = parameters(Ψ)

dudt = ODEFunction(sys)

estimator = ODEProblem(dudt, u0, tspan, p)
estimation = solve(estimator, Tsit5(), saveat = solution.t)

plot!(estimation,linecolor=:red,linewidth=3)
savefig("sindy_raw_estim")

using DiffEqParamEstim
using Optim
# cost_function = build_loss_objective(estimator,t,Xₙ,Tsit5(),maxiters=10000)

cost_function = build_loss_objective(estimator,Tsit5(),L2Loss(solution.t,Xₙ),
                                     maxiters=10000)#,verbose=false)
result = optimize(cost_function, p, BFGS())
tmp_estimator = remake(estimator, p=result.minimizer)
better_estimation = solve(tmp_estimator, Tsit5(), saveat = solution.t)

plot!(better_estimation,linecolor=:black,linestyle=:dash,linewidth=4)
savefig("sindy_better_estim")
#
# # The parameters are a bit off, but the equations are recovered
# # Start another SINDy run to get closer to the ground truth
# # Create function
# unknown_sys = ODESystem(Ψ)
# unknown_eq = ODEFunction(unknown_sys)
#
# # Just the equations
# b = Basis((u, p, t)->unknown_eq(u, [1.; 1.], t), u)
#
# # Retune for better parameters -> we could also use DiffEqFlux or other parameter estimation tools here.
# Ψf = SINDy(Xₙ[:, 2:end], L̂[:, 2:end], b, STRRidge(0.01), maxiter = 100, convergence_error = 1e-18) # Succeed
# println(Ψf)
# p̂ = parameters(Ψf)
# println("Second parameter guess : $(p̂)")
#
# # Create function
# recovered_sys = ODESystem(Ψf)
# recovered_eq = ODEFunction(recovered_sys)
# estimator = ODEProblem(dudt, u0, tspan, p)
# estimation = solve(estimator, Tsit5(), saveat = solution.t)
#
# plot!(estimation)
#
#
#

max_stars_repo_path: neuralODE modeling/bimolecular_nODE.jl

max_issues_repo_path: neuralODE modeling/bimolecular_nODE.jl

max_forks_repo_path: neuralODE modeling/bimolecular_nODE.jl


text:
struct DataCite <: DataRepo
end

base_url(::DataCite) = "https://api.datacite.org/works/"

description(repo::DataCite, mainpage) = miss_null(mainpage["attributes"]["description"])
data_fullname(::DataCite, mainpage) = mainpage["attributes"]["title"]
author(::DataCite, mainpage) = join.([[names[2] for names in value] for value in  mainpage["attributes"]["author"]], " ")
website(repo::DataCite, mainpage_url, mainpage) =  replace(mainpage_url, base_url(repo) => "https://doi.org/")
license(::DataCite, mainpage) = miss_null(mainpage["attributes"]["license"])
published_date(::DataCite, mainpage) = mainpage["attributes"]["published"]
dataset_cite(::DataCite, mainpage) = citation_text(mainpage["id"])

function paper_cite(::DataCite, mainpage)
    paper_cite = missing
    for related in mainpage["attributes"]["related-identifiers"]
        if related["relation-type-id"] == "IsSupplementTo"
            paper_cite = citation_text(related["related-identifier"])
            break
        end
    end
    paper_cite
end

function mainpage_url(repo::DataCite, dataname)
    if match_doi(dataname) != nothing
        identifier = match_doi(dataname)
        url = base_url(repo) * identifier
        JSON.parse(text_only(getpage(url).root))["data"], url
    else
        throw(GeneratorError(repo))
    end
end


text:
export estimate_error
export Errmeasure
export StandardSPMFErrmeasure
export ResidualErrmeasure
export DefaultErrmeasure
export EigvalReferenceErrmeasure
export ErrmeasureType

"""
    abstract type Errmeasure; end

Concrete subtypes of `Errmeasure` represent specific ways of
measuring the error of an eigenpair. NEP-solvers take such
an object as input. As a NEP-solver user,
you use the type as follows
```julia
julia> quasinewton(nep,errmeasure=ResidualErrmeasure(nep))
```
User-specified ways of measuring error can be given by
creating a new subtype of `Errmeasure` and using it as a
`errmeasure` keyword. You need to specify the way to
measure the error in the method `estimate_error`.

Note that in practice a `Function` can essentially
be used instead of a `Errmeasure`-object, which is a simple
way to have user-defined error measures.
See [`ErrmeasureType`](@ref).

See also: [`ErrmeasureType`](@ref), [`DefaultErrmeasure`](@ref), [`ResidualErrmeasure`](@ref), [`StandardSPMFErrmeasure`](@ref), [`estimate_error`](@ref), [`EigvalReferenceErrmeasure`](@ref).

"""
abstract type Errmeasure; end

"""
    ErrmeasureType = Union{Type{<:Errmeasure}, Function}

The `ErrmeasureType` represents (essentially) what you can insert in the
`errmeasure` keyword argument for most NEP-solvers. It can be
a function or an  [`Errmeasure`](@ref) object.
If it is a `Function` this function will be used to obtain
error estimate.

# Example
This shows how to compute a reference solution and
then use this as a reference solution. The error
in the second run will be effectively the
eigenvector error (appropriately normalized).
```julia
julia> using LinearAlgebra
julia> nep=nep_gallery("qdep0");
julia> (λref,vref)=quasinewton(nep,λ=-1,v=ones(size(nep,1)));
julia> myerrmeasure=(λ,v) -> norm(vref/vref[1]-v/v[1]);
julia> (λ,v)=quasinewton(nep,errmeasure=myerrmeasure,λ=-1.0 ,logger=1,tol=5e-13,v=ones(size(nep,1)));
Precomputing linsolver
iter 1 err:46.40296482739195 λ=-1.0 + 0.0im
iter 2 err:2.1592671533657883 λ=-0.7063330111559607 + 0.0im
iter 3 err:0.17079231439255405 λ=-0.8919579082730457 + 0.0im
iter 4 err:0.1633846991066227 λ=-1.0097584042560848 + 0.0im
iter 5 err:0.003434042059262583 λ=-1.0023823873044 + 0.0im
iter 6 err:0.0003182517281689052 λ=-1.0024660870524031 + 0.0im
iter 7 err:2.0105257231740345e-5 λ=-1.0024677891861997 + 0.0im
iter 8 err:1.618661190265619e-6 λ=-1.0024669496893164 + 0.0im
iter 9 err:1.233489068442819e-7 λ=-1.0024669918249076 + 0.0im
iter 10 err:9.44707811957546e-9 λ=-1.0024669883439166 + 0.0im
iter 11 err:7.867601351698812e-10 λ=-1.0024669886059947 + 0.0im
iter 12 err:0.0 λ=-1.002466988585764 + 0.0im

```
The eigenvalue error can be measured with the
[`EigvalReferenceErrmeasure`](@ref).

See also: [`Errmeasure`](@ref)
"""
ErrmeasureType = Union{Errmeasure, Function}

"""
    struct DefaultErrmeasure <: Errmeasure
    function DefaultErrmeasure(nep::NEP)

When you specify this `Errmeasure`, NEP-PACK tries to determine
a suitable `Errmeasure` based on the type of the `NEP`.
Note that this behavior may change in future versions.

See also: [`Errmeasure`](@ref)

"""
struct DefaultErrmeasure{X<:Errmeasure}<:Errmeasure;
    errm::X
    function DefaultErrmeasure(nep::NEP)
        if (nep isa AbstractSPMF)
            errm=StandardSPMFErrmeasure(nep);
        else
            errm=ResidualErrmeasure(nep);
        end
        return new{typeof(errm)}(errm);
    end
end

"""
    struct ResidualErrmeasure <: Errmeasure
    function ResidualErrmeasure(nep::NEP)

This `Errmeasure` species that the residual norm should be
used to measure the error.

See also: [`Errmeasure`](@ref)

"""
struct ResidualErrmeasure <: Errmeasure;
    nep::NEP
end

"""
    function estimate_error(E::ErrmeasureType,λ,v)

Returns the error estimate for the eigenpair `(λ,v)`. The way to measure
the error is specified in `E`, which can be an `Errmeasure` or a `Function`.

See also: [`Errmeasure`](@ref), [`ErrmeasureType`](@ref)

"""
function estimate_error(errm::ResidualErrmeasure, λ,v)
   return norm(compute_Mlincomb(errm.nep,λ,v))/norm(v);
end

function estimate_error(errm::DefaultErrmeasure, λ,v)
   return estimate_error(errm.errm,λ,v);
end

function estimate_error(errm::Function, λ,v)
   return errm(λ,v);
end

"""
    struct StandardSPMFErrmeasure <: Errmeasure
    function StandardSPMFErrmeasure(nep::AbstractSPMF)

This `Errmeasure` provides a way to compute the backward error.
The backward error estimate are only given for NEPs which are
subtypes of  `AbstractSPMF`. We use the Frobenius norm as the
matrix norm, since it is much cheaper to compute than the spectral
norm.

# Example
```julia
julia> nep=nep_gallery("qdep0");
julia> (λ,v)=quasinewton(nep,λ=-1,v=ones(size(nep,1)),errmeasure=StandardSPMFErrmeasure(nep),tol=1e-10,logger=1);
Precomputing linsolver
iter 1 err:0.022010375110869937 λ=-1.0 + 0.0im
iter 2 err:0.002515422247048546 λ=-0.7063330111559607 + 0.0im
iter 3 err:0.000892354247568813 λ=-0.8919579082730457 + 0.0im
iter 4 err:5.445678793151584e-5 λ=-1.0097584042560848 + 0.0im
iter 5 err:6.649967517409105e-7 λ=-1.0023823873044 + 0.0im
iter 6 err:1.0557281809769784e-8 λ=-1.0024660870524031 + 0.0im
iter 7 err:6.420125566431444e-9 λ=-1.0024677891861997 + 0.0im
iter 8 err:3.181093707909799e-10 λ=-1.0024669496893164 + 0.0im
iter 9 err:2.6368050026394416e-11 λ=-1.0024669918249076 + 0.0im

```
See also: [`Errmeasure`](@ref)

"""
struct StandardSPMFErrmeasure{X<:Real} <: Errmeasure
    nep::AbstractSPMF
    coeffs::Vector{X};
    function StandardSPMFErrmeasure(nep::AbstractSPMF)
        Av=get_Av(nep);
        # Note: norm(A) is a the frobenius norm in Julia
        coeffs=map(i->norm(Av[i]),1:size(Av,1));
        return new{eltype(coeffs)}(nep,coeffs);
    end
end

function estimate_error(errm::StandardSPMFErrmeasure, λ,v)
    Av=get_Av(errm.nep);
    fv=get_fv(errm.nep);
    denom=mapreduce(i->errm.coeffs[i]*abs(fv[i](λ)), +, 1:size(Av,1));
    return norm(compute_Mlincomb(errm.nep,λ,v))/(norm(v)*denom);
end

"""
    struct EigvalReferenceErrmeasure{X<:Number} <: Errmeasure
    function EigvalReferenceErrmeasure(nep,λref)

Use the difference between a precomputed λ-value (reference solution)
and the eigenvalue estimate
as the error measure.

# Example
```julia
julia> using LinearAlgebra
julia> nep=nep_gallery("qdep0");
julia> (λref,vref)=quasinewton(nep,λ=-1,v=ones(size(nep,1)));
julia> (λ,v)=quasinewton(nep,errmeasure=EigvalReferenceErrmeasure(nep,λref),λ=-1.0 ,logger=1,tol=5e-13,v=ones(size(nep,1)));
Precomputing linsolver
iter 1 err:0.002466988585763996 λ=-1.0 + 0.0im
iter 2 err:0.2961339774298033 λ=-0.7063330111559607 + 0.0im
iter 3 err:0.11050908031271833 λ=-0.8919579082730457 + 0.0im
iter 4 err:0.007291415670320767 λ=-1.0097584042560848 + 0.0im
iter 5 err:8.460128136400513e-5 λ=-1.0023823873044 + 0.0im
iter 6 err:9.015333608530796e-7 λ=-1.0024660870524031 + 0.0im
iter 7 err:8.006004357241636e-7 λ=-1.0024677891861997 + 0.0im
iter 8 err:3.889644761834177e-8 λ=-1.0024669496893164 + 0.0im
iter 9 err:3.2391436199930013e-9 λ=-1.0024669918249076 + 0.0im
iter 10 err:2.418474309706653e-10 λ=-1.0024669883439166 + 0.0im
iter 11 err:2.0230705999324528e-11 λ=-1.0024669886059947 + 0.0im
iter 12 err:0.0 λ=-1.002466988585764 + 0.0im
```

See also: [`Errmeasure`](@ref)

"""
struct EigvalReferenceErrmeasure{X<:Number} <: Errmeasure
    nep::NEP
    λref::X
    function EigvalReferenceErrmeasure(nep,λref)
        return new{typeof(λref)}(nep,λref);
    end
end

function estimate_error(errm::EigvalReferenceErrmeasure, λ,v)
    return abs(errm.λref-λ)
end


text:
using DashBootstrapComponents, DashHtmlComponents

PLOTLY_LOGO = "https://images.plot.ly/logo/new-branding/plotly-logomark.png";

search_bar = dbc_row(
    [
        dbc_col(dbc_input(type="search", placeholder="Search")),
        dbc_col(
            dbc_button("Search", color="primary", className="ml-2", n_clicks=0),
            width="auto",
        ),
    ],
    no_gutters=true,
    className="ml-auto flex-nowrap mt-3 mt-md-0",
    align="center",
);

navbar = dbc_navbar(
    [
        html_a(
            # Use row and col to control vertical alignment of logo / brand
            dbc_row(
                [
                    dbc_col(html_img(src=PLOTLY_LOGO, height="30px")),
                    dbc_col(dbc_navbarbrand("Navbar", className="ml-2")),
                ],
                align="center",
                no_gutters=true,
            ),
            href="https://plotly.com",
        ),
        dbc_navbartoggler(id="navbar-toggler", n_clicks=0),
        dbc_collapse(search_bar, id="navbar-collapse", navbar=true, is_open=false),
    ],
    color="dark",
    dark=true,
);

# add callback for toggling the collapse on small screens
callback!(
    app,
    Output("navbar-collapse", "is_open"),
    Input("navbar-toggler", "n_clicks"),
    State("navbar-collapse", "is_open"),
) do n, is_open
    return n > 0 ? is_open == 0 : is_open
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
import ..UncertainValues:
    AbstractUncertainValue,
    FittedDistribution,
    TheoreticalFittedUncertainScalar,
	UncertainScalarTheoreticalFit
import Distributions: Truncated
import StatsBase: quantile, std, mean
import Distributions: support

########################################################################
# Resampling without constraints
########################################################################
"""
	resample(uv::AbstractUncertainValue)

Sample the uncertain value once, drawing values from the entire support of the probability
distribution furnishing it.
"""
resample(uv::AbstractUncertainValue) = rand(uv.distribution)

"""
	resample(uv::AbstractUncertainValue, n::Int)

Sample the uncertain value `n` times, drawing values from the entire support of the
probability distribution furnishing it.
"""
resample(uv::AbstractUncertainValue, n::Int) =
    rand(uv.distribution, n)

"""
    resample(fd::FittedDistribution)

Resample a `FittedDistribution` instance once.
"""
resample(fd::FittedDistribution) = rand(fd.distribution)

"""
    resample(fd::FittedDistribution, n::Int)

Resample a `FittedDistribution` instance `n` times.
"""
resample(fd::FittedDistribution, n::Int) = rand(fd.distribution, n)

"""
    resample(ue::UncertainScalarTheoreticalFit)

Resample an `UncertainScalarTheoreticalFit` instance once.
"""
resample(ue::UncertainScalarTheoreticalFit) =
    rand(ue.distribution.distribution)

"""
    resample(ue::UncertainScalarTheoreticalFit, n::Int)

Resample an `UncertainScalarTheoreticalFit` instance `n` times.
"""
resample(ue::UncertainScalarTheoreticalFit, n::Int) =
    rand(ue.distribution.distribution, n)

"""
    resample(uv::TheoreticalFittedUncertainScalar, n::Int)

Resample an `TheoreticalFittedUncertainScalar` instance once.
"""
resample(uv::TheoreticalFittedUncertainScalar) =
    rand(uv.distribution.distribution)

"""
    resample(uv::TheoreticalFittedUncertainScalar, n::Int)

Resample an `TheoreticalFittedUncertainScalar` instance `n` times.
"""
resample(uv::TheoreticalFittedUncertainScalar, n::Int) =
    rand(uv.distribution.distribution, n)


text:
#################################################################
#
#    1st order derivation testing (loops)
#
#################################################################

### vector setting
    ex = quote
        a = zeros(2)
        for i in 1:2
            a[i] = x
        end
        sum(a)
    end
    compare(ex, 1.)

    b = [1,2,4,5]
    ex = quote
        a=zeros(1+4)
        for i in 1:4
            t = 4+3+2
            a[i] = b[i]+t-x
        end
        sum(a)
    end
    compare(ex, 1.)

    ex = quote
        a=zeros(1+3)
        for i in 1:4
            t = 4+3+2
            a[i] = b[i]*x+t
        end
        sum(a)
    end
    compare(ex, 1.)

### vector accumulation
    ex = quote
        a = zeros(2)
        for i in 1:2
            a[i] += x
        end
        sum(a)
    end
    compare(ex, 1.)

    ex = quote
        a=zeros(1+3)
        for i in 1:4
            t = 4+3+2
            a[i] += b[i]*x+t
        end
        sum(a)
    end
    compare(ex, 1.)

### same var setting (vector)
    ex = quote
        a = zeros(2)
        for i in 1:2
            a[1] = x
        end
        sum(a)
    end
    compare(ex, 2.)

    m.@deriv_rule %(x,y)      x     0
    m.@deriv_rule %(x,y)      y     0
    ex = quote
        a = zeros(2)
        for i in 1:4
            a[1 + i % 2] = x
        end
        sum(a)
    end
    compare(ex, 1.)

### same var accumulator (vector)
    ex = quote
        a = zeros(2)
        for i in 1:2
            a[1] += x
        end
        sum(a)
    end
    compare(ex, 2.)

### same var setting (scalar)
    ex = quote
        a = 0.
        for i in 1:4
            a = 4*x
        end
        a
    end
    compare(ex, 3.)

    ex = quote
        a = 0.
        for i in 1:4
            a = x
        end
        a
    end
    compare(ex, 3.)

    ex = quote
        a = 0.
        for i in 1:length(x)
            a = x[i]
        end
        a
    end
    compare(ex, [3., 2.])
    compare(ex, [1.])
    compare(ex, ones(10))

### same var accumulator (scalar)
    ex = quote
        a = 0.
        for i in 1:4
            a += x
        end
        a
    end
    compare(ex, 2.)

    ex = quote
        a = 0.
        for i in 1:4
            a += 2+x
        end
        a
    end
    compare(ex, 2.)

    ex = quote
        a = 0.
        for i in 1:3
            a += x^i
        end
        a
    end
    compare(ex, 2.)

    ex = quote
        a = 0.
        for i in 1:length(x)
            a += x[i]
        end
        a
    end
    compare(ex, [3., 2.])
    compare(ex, [1.])
    compare(ex, ones(10))

    ex = quote
        a = 0.
        for i in 1:length(x)
            a += x[1]^i
        end
        a+1
    end
    compare(ex, [3., 2.])
    compare(ex, [1.])
    compare(ex, ones(10))

### nested loops
    ex = quote
        a=0
        for i in 1:10
            for j in 1:10
                a += x
            end
        end
        a
    end
    compare(ex, 0.1)

    ex = quote
        a=0
        for i in 1:10
            for j in 1:10
                a += log(x) * sin(j)
            end
        end
        a
    end
    compare(ex, 0.1)

    ex = quote
        a=0
        for i in 1:10
            for j in 1:10
                a += (j < 4) * log(x) * sin(j)
            end
        end
        a
    end
    compare(ex, 0.1)


text:
using DashBootstrapComponents

alerts = html_div([
    dbc_alert(
        [
            "This is a primary alert with an ",
            html_a("example link", href = "#", className = "alert-link"),
        ],
        color = "primary",
    ),
    dbc_alert(
        [
            "This is a danger alert with an ",
            html_a("example link", href = "#", className = "alert-link"),
        ],
        color = "danger",
    ),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
export
    cost

############################################################################################
#                              COST METHODS                                                #
############################################################################################

"$(TYPEDSIGNATURES) Evaluate the cost at a knot point, and automatically handle terminal
knot point, multiplying by dt as necessary."
function stage_cost(cost::CostFunction, z::KnotPoint)
    if is_terminal(z)
        stage_cost(cost, state(z))
    else
        stage_cost(cost, state(z), control(z))*z.dt
    end
end

# Default to no integration
cost(obj, dyn_con::DynamicsConstraint{Q}, Z) where Q<:QuadratureRule = cost(obj, Z)

"Evaluate the cost for a trajectory (non-allocating)"
@inline function cost!(obj::Objective, Z::Traj)
    map!(stage_cost, obj.J, obj.cost, Z)
end

# In-place cost-expansion
function cost_expansion!(E::AbstractExpansion, cost::CostFunction, z::KnotPoint)
    cost_gradient!(E, cost, z)
    cost_hessian!(E, cost, z)
    E.x0 .= E.x  # copy cost-only gradient
    return nothing
end

function cost_expansion!(E::Vector{<:AbstractExpansion}, obj::Objective, Z::Traj)
    for k in eachindex(Z)
        z = Z[k]
        cost_expansion!(E[k], obj.cost[k], z)
    end
end

"```
Qx,Qu = cost_gradient(cost::CostFunction, z::KnotPoint)
```
Get Qx, Qu pieces of gradient of cost function, multiplied by dt"
function cost_gradient!(E, cost::CostFunction, z::KnotPoint)
    gradient!(E, cost, state(z), control(z))
    if is_terminal(z)
        dt_x = 1.0
        dt_u = 0.0
    else
        dt_x = z.dt
        dt_u = z.dt
    end
	E.x .*= dt_x
	E.u .*= dt_u
    return nothing
end

"```
Qxx,Quu,Qux = cost_hessian(cost::CostFunction, z::KnotPoint)
```
Get Qxx, Quu, Qux pieces of Hessian of cost function, multiplied by dt"
function cost_hessian!(E, cost::CostFunction, z::KnotPoint)
    hessian!(E, cost, state(z), control(z))
    if is_terminal(z)
        dt_x = 1.0
        dt_u = 0.0
    else
        dt_x = z.dt
        dt_u = z.dt
    end
	E.xx .*= dt_x
	E.uu .*= dt_u
	E.ux .*= dt_u
    return nothing
end

"""```
cost_gradient!(E::CostExpansion, obj::Objective, Z::Traj)
cost_gradient!(E::CostExpansion, obj::Objective, dyn_con::DynamicsConstraint{Q}, Z::Traj)
```
Calculate the cost gradient for an entire trajectory. If a dynamics constraint is given,
    use the appropriate integration rule, if defined.
"""
function cost_gradient!(E, obj::Objective, Z::Traj)
    N = length(Z)
    for k in eachindex(Z)
        cost_gradient!(E, obj[k], Z[k])
    end
end

# Default to no integration
cost_gradient!(E, obj, dyn_con::DynamicsConstraint{Q}, Z) where Q = cost_gradient!(E, obj, Z)

"```
cost_hessian!(E::CostExpansion, obj::Objective, Z::Traj)
```
 Calculate the cost Hessian for an entire trajectory"
function cost_hessian!(E, obj::Objective, Z::Traj)
    N = length(Z)
    for k in eachindex(Z)
        cost_hessian(E[k], obj[k], Z[k])
    end
end

"$(TYPEDSIGNATURES) Expand cost for entire trajectory"
function cost_expansion!(E, obj::Objective, Z::Traj)
    cost_gradient!(E, obj, Z)
    cost_hessian!(E, obj, Z)
end

"""
Compute the error expansion along the entire trajectory
	Assumes the cost expansion is already complete
"""
@inline error_expansion!(E::Vector{<:SizedCostExpansion}, model::AbstractModel, Z::Traj, G) = nothing
function error_expansion!(E::Vector{<:SizedCostExpansion}, model::RigidBody, Z::Traj, G)
    for k in eachindex(E)
        error_expansion!(E[k], model, Z[k], G[k])
    end
end

"Compute the error expansion at a single KnotPoint (only for state vectors not in Vector space)"
function error_expansion!(E::SizedCostExpansion{<:Any,N}, model::RigidBody, z, G) where N
    if N < 15
        G = SMatrix(G)
        E.u_  .= E.u
        E.uu_ .= E.uu
        E.ux_ .= SMatrix(E.ux)*G
        E.x_  .= G'*SVector(E.x)

        ∇²differential!(E.xx_, model, state(z), E.x0)
        E.xx_ .+= G'E.xx*G
    else
        E.u_ .= E.u
        E.uu_ .= E.uu
        mul!(E.ux_, E.ux, G)
        mul!(E.x_, Transpose(G), E.x)

        ∇²differential!(E.xx_, model, state(z), E.x0)
        mul!(E.tmp, E.xx, G)
        mul!(E.xx_, Transpose(G), E.tmp, 1.0, 1.0)
    end
end

"Copy the error expansion to another expansion"
function error_expansion!(E::AbstractExpansion, Q::SizedCostExpansion)
	E.x  .= Q.x
	E.u  .= Q.u
	E.xx .= Q.xx
	E.uu .= Q.uu
	E.ux .= Q.ux
end

"Get the error expansion"
@inline function error_expansion(E::SizedCostExpansion, model::RigidBody)
	return StaticExpansion(E.x_, E.xx_, E.u_, E.uu_, E.ux_)
end

"Get the error expansion (same as cost expansion)"
@inline function error_expansion(E::SizedCostExpansion, model::AbstractModel)
	return StaticExpansion(E.x, E.xx, E.u, E.uu, E.ux)
end

"Get cost expansion"
@inline function cost_expansion(E::SizedCostExpansion{<:Any,N,N}) where N
	return StaticExpansion(E.x, E.xx, E.u, E.uu, E.ux)
end


text:
# Make sure we're running from the examples directory.
cd(joinpath(@__DIR__(), ".."))

module Temp

using Eratosthenes # Bring in the top-level sim environment.
using Plots # For plotting (Pkg.add("Plots"); Pkg.add("PyPlot");)
import HDF5 # For loading the log file (Pkg.add("HDF5");)
using Printf

# Create the hierarchy of objects used by the sim.
scenario = setup("scenarios/reaction-wheel.yaml")

# Run it with text progress (e.g., from VS Code).
simulate(scenario) do k, n
    if (k-1) % 100 == 0
        @printf("Progress: % 5.1f%% done.\n", 100*(k-1)/(n-1))
    end
    return true
end

# Open the logs and plot some things.
if !isempty(scenario.sim.log_file)

    # Read some data from the logs. This do block takes care of
    # error-handling, so we can just read whatever we want and return it.
    t, q_BI, ω_BI_B, rw1, rw2, rw3 =
        HDF5.h5open(scenario.sim.log_file, "r") do logs
            (read(logs, "/smallsat1/body/state/time"),
             read(logs, "/smallsat1/body/state/data/q_BI"),
             read(logs, "/smallsat1/body/state/data/ω_BI_B"),
             read(logs, "/smallsat1/rw1/state/data/ω_RI_G"),
             read(logs, "/smallsat1/rw2/state/data/ω_RI_G"),
             read(logs, "/smallsat1/rw3/state/data/ω_RI_G"),
             )
        end

    # Scalars are logged as 1-by-n (don't know why). Convert to just n for
    # the sake of plotting.
    t = dropdims(t, dims=1)

    # Choose a friendly plotting package. PyPlot behaves nicely.
    pyplot()
    # plotlyjs()

    # Define each plot that we'll need.
    display(plot(t, transpose(q_BI),
                 label  = ["q1" "q2" "q3" "q4"],
                 xlabel = "Time (s)",
                 ylabel = "Attitude Quaternion"))
    display(plot(t, 180/π * transpose(ω_BI_B),
                 label  = ["ω1" "ω2" "ω3"],
                 xlabel = "Time (s)",
                 ylabel = "Rotation Rate (deg/s)"))

    display(plot(t, 30/π*[rw1[1,:] rw2[1,:] rw3[1,:]],
                 label  = ["rw1" "rw2" "rw3"],
                 xlabel = "Time (s)",
                 ylabel = "Wheel Rate (RPM)"))

end # plotting

end # Temp


text:
using Gadfly

using BlackBoxOptimizationBenchmarking
const BBOB = BlackBoxOptimizationBenchmarking

##

function plot_fun(f)
    r = 5

    c, n, nline = Geom.contour(levels=exp10.(range(-6, stop=10, length=40))), 800, 1000
    if any(f .== [BBOB.F16,BBOB.F17,BBOB.F18,BBOB.F19])
        c, n, nline = Geom.contour(levels=exp10.(range(-6, stop=8, length=20))), 500, 1000
    end

    p1 = plot(
        layer(x=[f.x_opt[1]],y=[f.x_opt[2]],Geom.point,Theme(default_color=colorant"red")),
        layer(
            z=(x,y)->f([x,y])-f.f_opt, x=range(-r, stop=r, length=n), y=range(-r, stop=r, length=n),
            c,
        ),
        Coord.cartesian(xmin=-r,xmax=r,ymin=-r,ymax=r),
        Guide.title(string(f))
    )

    # line plot
    x = collect(range(-r, stop=r, length=nline))

    p2 = plot(
        layer(x=x,y=[ f([xi,f.x_opt[2],f.x_opt[3]]) for xi in x]-f.f_opt, Geom.line),
        layer(x=x,y=[ f([f.x_opt[1],xi,f.x_opt[3]]) for xi in x]-f.f_opt, Geom.line,Theme(default_color=colorant"orange")),
        layer(x=x,y=[ f([f.x_opt[1],f.x_opt[2],xi]) for xi in x]-f.f_opt, Geom.line,Theme(default_color=colorant"red")),
        Coord.cartesian(xmin=-r,xmax=r),
        Guide.title(string(f)),
        Scale.y_log10
    )
    p1,p2
end

outdir = joinpath(Pkg.dir(),"BBOB","data","plots","functions")

for i in 20:length(enumerate(BBOB.BBOBFunction))

    f = getfield(BBOB,Symbol("F$i"))
    info(f)
    p1,p2 = plot_fun(f)

    draw(PDF(joinpath(outdir,"F$i.pdf"),16cm,16cm),p1)
    draw(SVG(joinpath(outdir,"F$i.svg"),16cm,16cm),p1)
    draw(PDF(joinpath(outdir,"1D","F$i.pdf"),16cm,12cm),p2)
    draw(SVG(joinpath(outdir,"1D","F$i.svg"),16cm,12cm),p2)
end

##


text:
using JuMP

m = Model()

type Recipe
	name::AbstractString
	time::Float64
	outs::Dict{AbstractString, Float64}
	ins::Dict{AbstractString, Float64}
	Recipe(name, time, speed) = new(name, time/speed, Dict{AbstractString, Float64}(), Dict{AbstractString, Float64}())
end

function production(r::Recipe, nfacs)
	@printf "%s: %0.1f\n" r.name nfacs
	@printf "\tProduces\n"
	for (k,v) in r.outs
		@printf "\t\t%0.2f %s per min / %0.2f per sec\n" 60v * nfacs k v*nfacs
	end
	@printf "\tConsumes\n"
	for (k,v) in r.ins
		@printf "\t\t%0.2f %s per min / %0.2f per sec\n" 60v * nfacs k v*nfacs
	end
end
@enum Rnum BPACK BATT WIRE GREEN BLUE RED SPEED1

recips = Dict{Rnum, Recipe}()

include("science_recipe.jl")

macro IN(rnum, iname, iqty)
	:(recips[$rnum].ins[$iname] = $iqty / recips[$rnum].time)
end

macro OUT(rnum, iname, iqty)
	:(recips[$rnum].outs[$iname] = $iqty / recips[$rnum].time)
end

recips[YPACK] = Recipe("high-tech-science-pack", 14, 0.5)
@IN YPACK "battery" 1
@IN YPACK "processing-unit" 3
@IN YPACK "speed-module" 1
@IN YPACK "copper-cable" 30
@OUT YPACK "high-tech-science-pack" 2.8

recips[BATT] = Recipe("battery", 5, 1.25)
@IN BATT "sulfuric-acid" 20
@IN BATT "iron-plate" 1
@IN BATT "copper-plate" 1
@OUT BATT "battery" 1

recips[BLUE] = Recipe("processing-unit", 10, 1.25)
@IN BLUE "electronic-circuit" 20
@IN BLUE "advanced-circuit" 2
@IN BLUE "sulfuric-acid" 5
@OUT BLUE "processing-unit" 1

recips[RED] = Recipe("advanced-circuit", 6, 1.25)
@IN RED "electronic-circuit" 2
@IN RED "plastic-bar" 2
@IN RED "copper-cable" 4
@OUT RED "advanced-circuit" 1

recips[GREEN] = Recipe("electronic-circuit", 0.5, 1.25)
@IN GREEN "iron-plate" 14
@IN GREEN "copper-cable" 3
@OUT GREEN "electronic-circuit" 1

recips[SPEED1] = Recipe("speed-module", 15, 1.25)
@IN SPEED1 "advanced-circuit" 5
@IN SPEED1 "electronic-circuit" 5
@OUT SPEED1  "speed-module" 1

recips[WIRE] = Recipe("copper-cable", 0.5, 1.25)
@IN WIRE "copper-plate" 1
@OUT WIRE "copper-cable" 2

@variable(m, YPACKs >= 0) #, Int)
@variable(m, BATTs >= 0) #, Int)
@variable(m, BLUEs >= 0) #, Int)
@variable(m, REDs >= 0) #, Int)
@variable(m, GREENs >= 0) #, Int)
@variable(m, SPEED1s >= 0) #, Int)
@variable(m, WIREs >= 0) #, Int)

macro FIN(Fn, F, R) :($Fn * recips[$F].ins[$R]) end
macro FOUT(Fn, F, R) :($Fn * recips[$F].outs[$R] )end

@constraint(m, YPACKs == 10)

@constraint(m, @FIN(YPACKs, YPACK, "copper-cable") + @FIN(GREENs, GREEN, "copper-cable") +  + @FIN(REDs, RED, "copper-cable") == @FOUT(WIREs, WIRE, "copper-cable"))
@constraint(m, @FIN(YPACKs, YPACK, "battery") == @FOUT(BATTs, BATT, "battery"))
@constraint(m, @FIN(YPACKs, YPACK, "processing-unit") == @FOUT(BLUEs, BLUE, "processing-unit"))
@constraint(m, @FIN(YPACKs, YPACK, "speed-module") == @FOUT(SPEED1s, SPEED1, "speed-module"))

@constraint(m, @FIN(BLUEs, BLUE, "electronic-circuit") + @FIN(REDs, RED, "electronic-circuit") + @FIN(SPEED1s, SPEED1, "electronic-circuit") + @FIN(SPEED1s, SPEED1, "electronic-circuit") == @FOUT(GREENs, GREEN, "electronic-circuit"))
@constraint(m, @FIN(BLUEs, BLUE, "advanced-circuit") + @FIN(SPEED1s, SPEED1, "advanced-circuit") == @FOUT(REDs, RED, "advanced-circuit"))

@objective(m, Min, sum([YPACKs, WIREs]))
solve(m)

production(recips[YPACK], getvalue(YPACKs))
production(recips[BATT], getvalue(BATTs))
production(recips[BLUE], getvalue(BLUEs))
production(recips[RED], getvalue(REDs))
production(recips[GREEN], getvalue(GREENs))
production(recips[WIRE], getvalue(WIREs))
production(recips[SPEED1], getvalue(SPEED1s))

max_stars_repo_name: lawless-m/factorio-utils

max_issues_repo_name: lawless-m/factorio-utils

max_forks_repo_name: lawless-m/factorio-utils


text:
module CDSAPI

using HTTP
using JSON
using Base64

export
    py2ju

"""
    retrieve(name, params, filename)

Retrieves data for `name` from the Climate Data Store
with the specified `params` and stores it in the current
directory as `filename`.
"""
function retrieve(name, params, filename)
    creds = Dict()
    open(joinpath(homedir(),".cdsapirc")) do f
        for line in readlines(f)
            key, val = strip.(split(line,':', limit=2))
            creds[key] = val
        end
    end

    apikey = string("Basic ", base64encode(creds["key"]))
    response = HTTP.request(
        "POST",
        creds["url"] * "/resources/$name",
        ["Authorization" => apikey],
        body=JSON.json(params),
        verbose=1)

    resp_dict = JSON.parse(String(response.body))
    data = Dict("state" => "queued")
    while data["state"] != "completed"
        data = HTTP.request("GET", creds["url"] * "/tasks/" * string(resp_dict["request_id"]),  ["Authorization" => apikey])
        data = JSON.parse(String(data.body))
        println("request queue status ", data["state"])
    end

    HTTP.download(data["location"], filename)
    return data
end

"""
    py2ju(dictstr)

Takes a Python dictionary as string and converts it into Julia's `Dict`

# Examples
```julia-repl
julia> str = \"""{
               'format': 'zip',
               'variable': 'surface_air_temperature',
               'product_type': 'climatology',
               'month': '08',
               'origin': 'era_interim',
           }\""";

julia> py2ju(str)
Dict{String,Any} with 5 entries:
  "format"       => "zip"
  "month"        => "08"
  "product_type" => "climatology"
  "variable"     => "surface_air_temperature"
  "origin"       => "era_interim"

```
"""
function py2ju(dictstr)
    dictstr_cpy = replace(dictstr, "'" => "\"")
    lastcomma_pos = findlast(",", dictstr_cpy).start

    # if there's no pair after the last comma
    if findnext(":", dictstr_cpy, lastcomma_pos) == nothing
        # remove the comma
        dictstr_cpy = dictstr_cpy[firstindex(dictstr_cpy):(lastcomma_pos - 1)] * dictstr_cpy[(lastcomma_pos + 1):lastindex(dictstr_cpy)]
    end

    # removes trailing comma from a list
    rx = r",[ \n\r\t]*\]"
    dictstr_cpy = replace(dictstr_cpy, rx => "]")

    return JSON.parse(dictstr_cpy)
end

end # module


text:
for (root, dirs, files,) in walkdir(dirname(@__FILE__); onerror=(x)->())
    for file in files
        if lowercase(strip(file)) == "manifest.toml"
            path = joinpath(root, file,)
            @info("Removing $(path)")
            try
                rm(path; force = true, recursive = true,)
            catch ex
                @warn("Ignoring error", exception=ex,)
            end
        end
    end
end

max_stars_repo_name: bcbi/BCBI-base

max_issues_repo_name: bcbi/BCBI-base

max_forks_repo_name: bcbi/BCBI-base


text:
using Test
using ClimateMachine.MoistThermodynamics
using NCDatasets
using Random
using RootSolvers
MT = MoistThermodynamics
using LinearAlgebra

using CLIMAParameters
using CLIMAParameters.Planet

struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

# Tolerances for tested quantities:
const atol_temperature = 1e-1
const atol_pressure = MSLP(param_set) * 2e-2
const atol_energy = 1e-1
const rtol_temperature = 1e-1
const rtol_pressure = 1e-1
const rtol_energy = 1e-1

float_types = [Float32, Float64]

include("profiles.jl")
include("data_tests.jl")

@testset "moist thermodynamics - isentropic processes" begin
    for FT in [Float64]
        _R_d = FT(R_d(param_set))
        _molmass_ratio = FT(molmass_ratio(param_set))
        _cp_d = FT(cp_d(param_set))
        _cp_v = FT(cp_v(param_set))
        _cp_l = FT(cp_l(param_set))
        _cv_d = FT(cv_d(param_set))
        _cv_v = FT(cv_v(param_set))
        _cv_l = FT(cv_l(param_set))
        _cv_i = FT(cv_i(param_set))
        _T_0 = FT(T_0(param_set))
        _e_int_v0 = FT(e_int_v0(param_set))
        _e_int_i0 = FT(e_int_i0(param_set))
        _LH_v0 = FT(LH_v0(param_set))
        _LH_s0 = FT(LH_s0(param_set))
        _cp_i = FT(cp_i(param_set))
        _LH_f0 = FT(LH_f0(param_set))
        _press_triple = FT(press_triple(param_set))
        _R_v = FT(R_v(param_set))
        _T_triple = FT(T_triple(param_set))
        _T_freeze = FT(T_freeze(param_set))
        _T_min = FT(T_min(param_set))
        _MSLP = FT(MSLP(param_set))
        _T_max = FT(T_max(param_set))
        _kappa_d = FT(kappa_d(param_set))

        # for FT in float_types
        z, e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice =
            tested_profiles(param_set, 50, FT)
        Φ = FT(1)
        Random.seed!(15)
        perturbation = FT(0.1) * rand(length(T))

        # TODO: Use reasonable values for ambient temperature/pressure
        T∞, p∞ = T .* perturbation, p .* perturbation
        @test air_temperature.(
            Ref(param_set),
            p,
            θ_liq_ice,
            Ref(DryAdiabaticProcess()),
        ) ≈ (p ./ _MSLP) .^ (_R_d / _cp_d) .* θ_liq_ice
        @test air_pressure_given_θ.(
            Ref(param_set),
            θ_liq_ice,
            Φ,
            Ref(DryAdiabaticProcess()),
        ) ≈ _MSLP .* (1 .- Φ ./ (θ_liq_ice .* _cp_d)) .^ (_cp_d / _R_d)
        @test air_pressure.(
            Ref(param_set),
            T,
            T∞,
            p∞,
            Ref(DryAdiabaticProcess()),
        ) ≈ p∞ .* (T ./ T∞) .^ (FT(1) / _kappa_d)
    end
end

@testset "moist thermodynamics - correctness" begin
    FT = Float64
    _R_d = FT(R_d(param_set))
    _molmass_ratio = FT(molmass_ratio(param_set))
    _cp_d = FT(cp_d(param_set))
    _cp_v = FT(cp_v(param_set))
    _cp_l = FT(cp_l(param_set))
    _cv_d = FT(cv_d(param_set))
    _cv_v = FT(cv_v(param_set))
    _cv_l = FT(cv_l(param_set))
    _cv_i = FT(cv_i(param_set))
    _T_0 = FT(T_0(param_set))
    _e_int_v0 = FT(e_int_v0(param_set))
    _e_int_i0 = FT(e_int_i0(param_set))
    _LH_v0 = FT(LH_v0(param_set))
    _LH_s0 = FT(LH_s0(param_set))
    _cp_i = FT(cp_i(param_set))
    _LH_f0 = FT(LH_f0(param_set))
    _press_triple = FT(press_triple(param_set))
    _R_v = FT(R_v(param_set))
    _T_triple = FT(T_triple(param_set))
    _T_freeze = FT(T_freeze(param_set))
    _T_min = FT(T_min(param_set))
    _MSLP = FT(MSLP(param_set))
    _T_max = FT(T_max(param_set))
    _kappa_d = FT(kappa_d(param_set))
    _T_icenuc = FT(T_icenuc(param_set))

    # ideal gas law
    @test air_pressure(param_set, FT(1), FT(1), PhasePartition(FT(1))) === _R_v
    @test air_pressure(
        param_set,
        FT(1),
        FT(2),
        PhasePartition(FT(1), FT(0.5), FT(0)),
    ) === _R_v
    @test air_pressure(param_set, FT(1), FT(1)) === _R_d
    @test air_pressure(param_set, FT(1), FT(2)) === 2 * _R_d
    @test air_density(param_set, FT(1), FT(1)) === 1 / _R_d
    @test air_density(param_set, FT(1), FT(2)) === 2 / _R_d

    # gas constants and heat capacities
    @test gas_constant_air(param_set, PhasePartition(FT(0))) === _R_d
    @test gas_constant_air(param_set, PhasePartition(FT(1))) === _R_v
    @test gas_constant_air(param_set, PhasePartition(FT(0.5), FT(0.5))) ≈
          _R_d / 2
    @test gas_constant_air(param_set, FT) == _R_d

    @test cp_m(param_set, PhasePartition(FT(0))) === _cp_d
    @test cp_m(param_set, PhasePartition(FT(1))) === _cp_v
    @test cp_m(param_set, PhasePartition(FT(1), FT(1))) === _cp_l
    @test cp_m(param_set, PhasePartition(FT(1), FT(0), FT(1))) === _cp_i
    @test cp_m(param_set, FT) == _cp_d

    @test cv_m(param_set, PhasePartition(FT(0))) === _cp_d - _R_d
    @test cv_m(param_set, PhasePartition(FT(1))) === _cp_v - _R_v
    @test cv_m(param_set, PhasePartition(FT(1), FT(1))) === _cv_l
    @test cv_m(param_set, PhasePartition(FT(1), FT(0), FT(1))) === _cv_i
    @test cv_m(param_set, FT) == _cv_d

    # speed of sound
    @test soundspeed_air(param_set, _T_0 + 20, PhasePartition(FT(0))) ==
          sqrt(_cp_d / _cv_d * _R_d * (_T_0 + 20))
    @test soundspeed_air(param_set, _T_0 + 100, PhasePartition(FT(1))) ==
          sqrt(_cp_v / _cv_v * _R_v * (_T_0 + 100))

    # specific latent heats
    @test latent_heat_vapor(param_set, _T_0) ≈ _LH_v0
    @test latent_heat_fusion(param_set, _T_0) ≈ _LH_f0
    @test latent_heat_sublim(param_set, _T_0) ≈ _LH_s0

    # saturation vapor pressure and specific humidity
    p = FT(1.e5)
    q_tot = FT(0.23)
    ρ = FT(1.0)
    ρ_v_triple = _press_triple / _R_v / _T_triple
    @test saturation_vapor_pressure(param_set, _T_triple, Liquid()) ≈
          _press_triple
    @test saturation_vapor_pressure(param_set, _T_triple, Ice()) ≈ _press_triple

    @test q_vap_saturation(param_set, _T_triple, ρ, PhasePartition(FT(0))) ==
          ρ_v_triple / ρ
    @test q_vap_saturation(
        param_set,
        _T_triple,
        ρ,
        PhasePartition(q_tot, q_tot),
    ) == ρ_v_triple / ρ

    @test q_vap_saturation_generic(param_set, _T_triple, ρ; phase = Liquid()) ==
          ρ_v_triple / ρ
    @test q_vap_saturation_generic(param_set, _T_triple, ρ; phase = Ice()) ==
          ρ_v_triple / ρ
    @test q_vap_saturation_generic(
        param_set,
        _T_triple - 20,
        ρ;
        phase = Liquid(),
    ) >= q_vap_saturation_generic(param_set, _T_triple - 20, ρ; phase = Ice())

    @test saturation_excess(param_set, _T_triple, ρ, PhasePartition(q_tot)) ==
          q_tot - ρ_v_triple / ρ
    @test saturation_excess(
        param_set,
        _T_triple,
        ρ,
        PhasePartition(q_tot / 1000),
    ) == 0.0

    # energy functions and inverse (temperature)
    T = FT(300)
    e_kin = FT(11)
    e_pot = FT(13)
    @test air_temperature(param_set, _cv_d * (T - _T_0)) === FT(T)
    @test air_temperature(
        param_set,
        _cv_d * (T - _T_0),
        PhasePartition(FT(0)),
    ) === FT(T)

    @test air_temperature(
        param_set,
        cv_m(param_set, PhasePartition(FT(0))) * (T - _T_0),
        PhasePartition(FT(0)),
    ) === FT(T)
    @test air_temperature(
        param_set,
        cv_m(param_set, PhasePartition(FT(q_tot))) * (T - _T_0) +
        q_tot * _e_int_v0,
        PhasePartition(q_tot),
    ) ≈ FT(T)

    @test total_energy(param_set, FT(e_kin), FT(e_pot), _T_0) ===
          FT(e_kin) + FT(e_pot)
    @test total_energy(param_set, FT(e_kin), FT(e_pot), FT(T)) ≈
          FT(e_kin) + FT(e_pot) + _cv_d * (T - _T_0)
    @test total_energy(param_set, FT(0), FT(0), _T_0, PhasePartition(q_tot)) ≈
          q_tot * _e_int_v0

    # phase partitioning in equilibrium
    q_liq = FT(0.1)
    T = FT(_T_icenuc - 10)
    ρ = FT(1.0)
    q_tot = FT(0.21)
    @test liquid_fraction(param_set, T) === FT(0)
    @test liquid_fraction(param_set, T, PhasePartition(q_tot, q_liq, q_liq)) ===
          FT(0.5)
    q = PhasePartition_equil(param_set, T, ρ, q_tot)
    @test q.liq ≈ FT(0)
    @test 0 < q.ice <= q_tot

    T = FT(_T_freeze + 10)
    ρ = FT(0.1)
    q_tot = FT(0.60)
    @test liquid_fraction(param_set, T) === FT(1)
    @test liquid_fraction(
        param_set,
        T,
        PhasePartition(q_tot, q_liq, q_liq / 2),
    ) === FT(2 / 3)
    q = PhasePartition_equil(param_set, T, ρ, q_tot)
    @test 0 < q.liq <= q_tot
    @test q.ice ≈ 0

    # saturation adjustment in equilibrium (i.e., given the thermodynamic
    # variables E_int, p, q_tot, compute the temperature and partitioning of the phases
    tol_T = 1e-1
    q_tot = FT(0)
    ρ = FT(1)
    @test MT.saturation_adjustment_SecantMethod(
        param_set,
        internal_energy_sat(param_set, 300.0, ρ, q_tot),
        ρ,
        q_tot,
        10,
        ResidualTolerance(1e-2),
    ) ≈ 300.0
    @test abs(
        MT.saturation_adjustment(
            param_set,
            internal_energy_sat(param_set, 300.0, ρ, q_tot),
            ρ,
            q_tot,
            10,
            ResidualTolerance(1e-2),
        ) - 300.0,
    ) < tol_T

    q_tot = FT(0.21)
    ρ = FT(0.1)
    @test MT.saturation_adjustment_SecantMethod(
        param_set,
        internal_energy_sat(param_set, 200.0, ρ, q_tot),
        ρ,
        q_tot,
        10,
        ResidualTolerance(1e-2),
    ) ≈ 200.0
    @test abs(
        MT.saturation_adjustment(
            param_set,
            internal_energy_sat(param_set, 200.0, ρ, q_tot),
            ρ,
            q_tot,
            10,
            ResidualTolerance(1e-2),
        ) - 200.0,
    ) < tol_T
    q = PhasePartition_equil(param_set, T, ρ, q_tot)
    @test q.tot - q.liq - q.ice ≈
          vapor_specific_humidity(q) ≈
          q_vap_saturation(param_set, T, ρ)

    ρ = FT(1)
    ρu = FT[1, 2, 3]
    ρe = FT(1100)
    e_pot = FT(93)
    @test internal_energy(ρ, ρe, ρu, e_pot) ≈ 1000.0

    # potential temperatures
    T = FT(300)
    @test liquid_ice_pottemp_given_pressure(param_set, T, _MSLP) === T
    @test liquid_ice_pottemp_given_pressure(param_set, T, _MSLP / 10) ≈
          T * 10^(_R_d / _cp_d)
    @test liquid_ice_pottemp_given_pressure(
        param_set,
        T,
        _MSLP / 10,
        PhasePartition(FT(1)),
    ) ≈ T * 10^(_R_v / _cp_v)

    # dry potential temperatures. FIXME: add correctness tests
    T = FT(300)
    p = FT(1.e5)
    q_tot = FT(0.23)
    @test dry_pottemp_given_pressure(param_set, T, p, PhasePartition(q_tot)) isa
          typeof(p)
    @test air_temperature_from_liquid_ice_pottemp_given_pressure(
        param_set,
        dry_pottemp_given_pressure(param_set, T, p, PhasePartition(q_tot)),
        p,
        PhasePartition(q_tot),
    ) ≈ T

    # Exner function. FIXME: add correctness tests
    p = FT(1.e5)
    q_tot = FT(0.23)
    @test exner_given_pressure(param_set, p, PhasePartition(q_tot)) isa
          typeof(p)
end

@testset "moist thermodynamics - default behavior accuracy" begin
    # Input arguments should be accurate within machine precision
    # Temperature is approximated via saturation adjustment, and should be within a physical tolerance

    for FT in float_types
        rtol = FT(1e-2)
        z, e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice =
            tested_profiles(param_set, 50, FT)

        # PhaseEquil
        ts_exact = PhaseEquil.(Ref(param_set), e_int, ρ, q_tot, 100, FT(1e-3))
        ts = PhaseEquil.(Ref(param_set), e_int, ρ, q_tot)
        # Should be machine accurate (because ts contains `e_int`,`ρ`,`q_tot`):
        @test all(
            getproperty.(PhasePartition.(ts), :tot) .≈
            getproperty.(PhasePartition.(ts_exact), :tot),
        )
        @test all(internal_energy.(ts) .≈ internal_energy.(ts_exact))
        @test all(air_density.(ts) .≈ air_density.(ts_exact))
        # Approximate (temperature must be computed via saturation adjustment):
        @test all(isapprox.(
            air_pressure.(ts),
            air_pressure.(ts_exact),
            rtol = rtol,
        ))
        @test all(isapprox.(
            air_temperature.(ts),
            air_temperature.(ts_exact),
            rtol = rtol,
        ))

        # PhaseEquil
        ts_exact =
            PhaseEquil.(
                Ref(param_set),
                e_int,
                ρ,
                q_tot,
                100,
                FT(1e-3),
                MT.saturation_adjustment_SecantMethod,
            )
        ts =
            PhaseEquil.(
                Ref(param_set),
                e_int,
                ρ,
                q_tot,
                35,
                FT(1e-1),
                MT.saturation_adjustment_SecantMethod,
            ) # Needs to be in sync with default
        # Should be machine accurate (because ts contains `e_int`,`ρ`,`q_tot`):
        @test all(
            getproperty.(PhasePartition.(ts), :tot) .≈
            getproperty.(PhasePartition.(ts_exact), :tot),
        )
        @test all(internal_energy.(ts) .≈ internal_energy.(ts_exact))
        @test all(air_density.(ts) .≈ air_density.(ts_exact))
        # Approximate (temperature must be computed via saturation adjustment):
        @test all(isapprox.(
            air_pressure.(ts),
            air_pressure.(ts_exact),
            rtol = rtol,
        ))
        @test all(isapprox.(
            air_temperature.(ts),
            air_temperature.(ts_exact),
            rtol = rtol,
        ))

        # LiquidIcePotTempSHumEquil
        ts_exact =
            LiquidIcePotTempSHumEquil.(
                Ref(param_set),
                θ_liq_ice,
                ρ,
                q_tot,
                40,
                FT(1e-3),
            )
        ts = LiquidIcePotTempSHumEquil.(Ref(param_set), θ_liq_ice, ρ, q_tot)
        # Should be machine accurate:
        @test all(air_density.(ts) .≈ air_density.(ts_exact))
        @test all(
            getproperty.(PhasePartition.(ts), :tot) .≈
            getproperty.(PhasePartition.(ts_exact), :tot),
        )
        # Approximate (temperature must be computed via saturation adjustment):
        @test all(isapprox.(
            internal_energy.(ts),
            internal_energy.(ts_exact),
            rtol = 10 * rtol,
        ))
        @test all(isapprox.(
            liquid_ice_pottemp.(ts),
            liquid_ice_pottemp.(ts_exact),
            rtol = rtol,
        ))
        @test all(isapprox.(
            air_temperature.(ts),
            air_temperature.(ts_exact),
            rtol = rtol,
        ))

        # LiquidIcePotTempSHumEquil_given_pressure
        ts_exact =
            LiquidIcePotTempSHumEquil_given_pressure.(
                Ref(param_set),
                θ_liq_ice,
                p,
                q_tot,
                40,
                FT(1e-3),
            )
        ts =
            LiquidIcePotTempSHumEquil_given_pressure.(
                Ref(param_set),
                θ_liq_ice,
                p,
                q_tot,
            )
        # Should be machine accurate:
        @test all(
            getproperty.(PhasePartition.(ts), :tot) .≈
            getproperty.(PhasePartition.(ts_exact), :tot),
        )
        # Approximate (temperature must be computed via saturation adjustment):
        @test all(isapprox.(
            air_density.(ts),
            air_density.(ts_exact),
            rtol = rtol,
        ))
        @test all(isapprox.(
            internal_energy.(ts),
            internal_energy.(ts_exact),
            rtol = rtol,
        ))
        @test all(isapprox.(
            liquid_ice_pottemp.(ts),
            liquid_ice_pottemp.(ts_exact),
            rtol = rtol,
        ))
        @test all(isapprox.(
            air_temperature.(ts),
            air_temperature.(ts_exact),
            rtol = rtol,
        ))

        # LiquidIcePotTempSHumNonEquil
        ts_exact =
            LiquidIcePotTempSHumNonEquil.(
                Ref(param_set),
                θ_liq_ice,
                ρ,
                q_pt,
                40,
                FT(1e-3),
            )
        ts = LiquidIcePotTempSHumNonEquil.(Ref(param_set), θ_liq_ice, ρ, q_pt)
        # Should be machine accurate:
        @test all(
            getproperty.(PhasePartition.(ts), :tot) .≈
            getproperty.(PhasePartition.(ts_exact), :tot),
        )
        @test all(
            getproperty.(PhasePartition.(ts), :liq) .≈
            getproperty.(PhasePartition.(ts_exact), :liq),
        )
        @test all(
            getproperty.(PhasePartition.(ts), :ice) .≈
            getproperty.(PhasePartition.(ts_exact), :ice),
        )
        @test all(air_density.(ts) .≈ air_density.(ts_exact))
        # Approximate (temperature must be computed via non-linear solve):
        @test all(isapprox.(
            internal_energy.(ts),
            internal_energy.(ts_exact),
            rtol = rtol,
        ))
        @test all(isapprox.(
            liquid_ice_pottemp.(ts),
            liquid_ice_pottemp.(ts_exact),
            rtol = rtol,
        ))
        @test all(isapprox.(
            air_temperature.(ts),
            air_temperature.(ts_exact),
            rtol = rtol,
        ))

    end

end

@testset "moist thermodynamics - constructor consistency" begin

    # Make sure `ThermodynamicState` arguments are returned unchanged

    for FT in float_types
        rtol = FT(1e-2)

        _MSLP = FT(MSLP(param_set))

        z, e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice =
            tested_profiles(param_set, 50, FT)

        # PhaseDry
        ts = PhaseDry.(Ref(param_set), e_int, ρ)
        @test all(internal_energy.(ts) .≈ e_int)
        @test all(air_density.(ts) .≈ ρ)

        ts = PhaseDry_given_pT.(Ref(param_set), p, T)
        @test all(internal_energy.(ts) .≈ internal_energy.(Ref(param_set), T))
        @test all(air_density.(ts) .≈ ρ)

        # PhaseEquil
        ts =
            PhaseEquil.(
                Ref(param_set),
                e_int,
                ρ,
                q_tot,
                40,
                FT(1e-1),
                Ref(MT.saturation_adjustment_SecantMethod),
            )
        @test all(internal_energy.(ts) .≈ e_int)
        @test all(getproperty.(PhasePartition.(ts), :tot) .≈ q_tot)
        @test all(air_density.(ts) .≈ ρ)

        ts = PhaseEquil.(Ref(param_set), e_int, ρ, q_tot)
        @test all(internal_energy.(ts) .≈ e_int)
        @test all(getproperty.(PhasePartition.(ts), :tot) .≈ q_tot)
        @test all(air_density.(ts) .≈ ρ)

        # PhaseNonEquil
        ts = PhaseNonEquil.(Ref(param_set), e_int, ρ, q_pt)
        @test all(internal_energy.(ts) .≈ e_int)
        @test all(
            getproperty.(PhasePartition.(ts), :tot) .≈ getproperty.(q_pt, :tot),
        )
        @test all(
            getproperty.(PhasePartition.(ts), :liq) .≈ getproperty.(q_pt, :liq),
        )
        @test all(
            getproperty.(PhasePartition.(ts), :ice) .≈ getproperty.(q_pt, :ice),
        )
        @test all(air_density.(ts) .≈ ρ)

        # air_temperature_from_liquid_ice_pottemp_given_pressure-liquid_ice_pottemp inverse
        θ_liq_ice_ =
            liquid_ice_pottemp_given_pressure.(Ref(param_set), T, p, q_pt)
        @test all(
            air_temperature_from_liquid_ice_pottemp_given_pressure.(
                Ref(param_set),
                θ_liq_ice_,
                p,
                q_pt,
            ) .≈ T,
        )

        # liquid_ice_pottemp-air_temperature_from_liquid_ice_pottemp_given_pressure inverse
        T =
            air_temperature_from_liquid_ice_pottemp_given_pressure.(
                Ref(param_set),
                θ_liq_ice,
                p,
                q_pt,
            )
        @test all(
            liquid_ice_pottemp_given_pressure.(Ref(param_set), T, p, q_pt) .≈
            θ_liq_ice,
        )

        # Accurate but expensive `LiquidIcePotTempSHumNonEquil` constructor (Non-linear temperature from θ_liq_ice)
        T_non_linear =
            air_temperature_from_liquid_ice_pottemp_non_linear.(
                Ref(param_set),
                θ_liq_ice,
                ρ,
                20,
                ResidualTolerance(FT(5e-5)),
                q_pt,
            )
        T_expansion =
            air_temperature_from_liquid_ice_pottemp.(
                Ref(param_set),
                θ_liq_ice,
                ρ,
                q_pt,
            )
        @test all(isapprox.(T_non_linear, T_expansion, rtol = rtol))
        e_int_ = internal_energy.(Ref(param_set), T_non_linear, q_pt)
        ts = PhaseNonEquil.(Ref(param_set), e_int_, ρ, q_pt)
        @test all(T_non_linear .≈ air_temperature.(ts))
        @test all(isapprox(θ_liq_ice, liquid_ice_pottemp.(ts), rtol = rtol))

        # LiquidIcePotTempSHumEquil
        ts =
            LiquidIcePotTempSHumEquil.(
                Ref(param_set),
                θ_liq_ice,
                ρ,
                q_tot,
                40,
                FT(1e-3),
            )
        @test all(isapprox.(liquid_ice_pottemp.(ts), θ_liq_ice, atol = 1e-1))
        @test all(isapprox.(air_density.(ts), ρ, rtol = rtol))
        @test all(getproperty.(PhasePartition.(ts), :tot) .≈ q_tot)

        # The LiquidIcePotTempSHumEquil_given_pressure constructor
        # passes the consistency test within sufficient physical precision,
        # however, it fails to satisfy the consistency test within machine
        # precision for the input pressure.

        # LiquidIcePotTempSHumEquil_given_pressure
        ts =
            LiquidIcePotTempSHumEquil_given_pressure.(
                Ref(param_set),
                θ_liq_ice,
                p,
                q_tot,
                35,
                FT(1e-3),
            )
        @test all(isapprox.(liquid_ice_pottemp.(ts), θ_liq_ice, atol = 1e-1))
        @test all(
            getproperty.(PhasePartition.(ts), :tot) .≈ getproperty.(q_pt, :tot),
        )
        @test all(isapprox.(air_pressure.(ts), p, atol = _MSLP * 2e-2))

        # LiquidIcePotTempSHumNonEquil_given_pressure
        ts =
            LiquidIcePotTempSHumNonEquil_given_pressure.(
                Ref(param_set),
                θ_liq_ice,
                p,
                q_pt,
            )
        @test all(liquid_ice_pottemp.(ts) .≈ θ_liq_ice)
        @test all(air_pressure.(ts) .≈ p)
        @test all(
            getproperty.(PhasePartition.(ts), :tot) .≈ getproperty.(q_pt, :tot),
        )
        @test all(
            getproperty.(PhasePartition.(ts), :liq) .≈ getproperty.(q_pt, :liq),
        )
        @test all(
            getproperty.(PhasePartition.(ts), :ice) .≈ getproperty.(q_pt, :ice),
        )

        # LiquidIcePotTempSHumNonEquil
        ts =
            LiquidIcePotTempSHumNonEquil.(
                Ref(param_set),
                θ_liq_ice,
                ρ,
                q_pt,
                5,
                FT(1e-3),
            )
        @test all(isapprox.(θ_liq_ice, liquid_ice_pottemp.(ts), rtol = rtol))
        @test all(air_density.(ts) .≈ ρ)
        @test all(
            getproperty.(PhasePartition.(ts), :tot) .≈ getproperty.(q_pt, :tot),
        )
        @test all(
            getproperty.(PhasePartition.(ts), :liq) .≈ getproperty.(q_pt, :liq),
        )
        @test all(
            getproperty.(PhasePartition.(ts), :ice) .≈ getproperty.(q_pt, :ice),
        )
    end

end

@testset "moist thermodynamics - type-stability" begin

    # NOTE: `Float32` saturation adjustment tends to have more difficulty
    # with converging to the same tolerances as `Float64`, so they're relaxed here.
    FT = Float32
    z, e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice =
        tested_profiles(param_set, 50, FT)

    ρu = FT[1.0, 2.0, 3.0]
    e_pot = FT(100.0)
    @test typeof.(internal_energy.(ρ, ρ .* e_int, Ref(ρu), Ref(e_pot))) ==
          typeof.(e_int)

    ts_dry = PhaseDry.(Ref(param_set), e_int, ρ)
    ts_dry_pT = PhaseDry_given_pT.(Ref(param_set), p, T)
    ts_eq = PhaseEquil.(Ref(param_set), e_int, ρ, q_tot, 15, FT(1e-1))
    ts_T =
        TemperatureSHumEquil.(
            Ref(param_set),
            air_temperature.(ts_dry),
            air_pressure.(ts_dry),
            q_tot,
        )
    ts_neq = PhaseNonEquil.(Ref(param_set), e_int, ρ, q_pt)
    ts_θ_liq_ice_eq =
        LiquidIcePotTempSHumEquil.(
            Ref(param_set),
            θ_liq_ice,
            ρ,
            q_tot,
            40,
            FT(1e-3),
        )
    ts_θ_liq_ice_eq_p =
        LiquidIcePotTempSHumEquil_given_pressure.(
            Ref(param_set),
            θ_liq_ice,
            p,
            q_tot,
            40,
            FT(1e-3),
        )
    ts_θ_liq_ice_neq =
        LiquidIcePotTempSHumNonEquil.(Ref(param_set), θ_liq_ice, ρ, q_pt)
    ts_θ_liq_ice_neq_p =
        LiquidIcePotTempSHumNonEquil_given_pressure.(
            Ref(param_set),
            θ_liq_ice,
            p,
            q_pt,
        )

    for ts in (
        ts_dry,
        ts_dry_pT,
        ts_eq,
        ts_T,
        ts_neq,
        ts_θ_liq_ice_eq,
        ts_θ_liq_ice_eq_p,
        ts_θ_liq_ice_neq,
        ts_θ_liq_ice_neq_p,
    )
        @test typeof.(soundspeed_air.(ts)) == typeof.(e_int)
        @test typeof.(gas_constant_air.(ts)) == typeof.(e_int)
        @test typeof.(vapor_specific_humidity.(ts)) == typeof.(e_int)
        @test typeof.(relative_humidity.(ts)) == typeof.(e_int)
        @test typeof.(air_pressure.(ts)) == typeof.(e_int)
        @test typeof.(air_density.(ts)) == typeof.(e_int)
        @test typeof.(total_specific_humidity.(ts)) == typeof.(e_int)
        @test typeof.(cp_m.(ts)) == typeof.(e_int)
        @test typeof.(cv_m.(ts)) == typeof.(e_int)
        @test typeof.(air_temperature.(ts)) == typeof.(e_int)
        @test typeof.(internal_energy_sat.(ts)) == typeof.(e_int)
        @test typeof.(internal_energy.(ts)) == typeof.(e_int)
        @test typeof.(latent_heat_vapor.(ts)) == typeof.(e_int)
        @test typeof.(latent_heat_sublim.(ts)) == typeof.(e_int)
        @test typeof.(latent_heat_fusion.(ts)) == typeof.(e_int)
        @test typeof.(q_vap_saturation.(ts)) == typeof.(e_int)
        @test typeof.(saturation_excess.(ts)) == typeof.(e_int)
        @test typeof.(liquid_fraction.(ts)) == typeof.(e_int)
        @test typeof.(liquid_ice_pottemp.(ts)) == typeof.(e_int)
        @test typeof.(dry_pottemp.(ts)) == typeof.(e_int)
        @test typeof.(exner.(ts)) == typeof.(e_int)
        @test typeof.(liquid_ice_pottemp_sat.(ts)) == typeof.(e_int)
        @test typeof.(specific_volume.(ts)) == typeof.(e_int)
        @test typeof.(virtual_pottemp.(ts)) == typeof.(e_int)
        @test eltype.(gas_constants.(ts)) == typeof.(e_int)
        @test typeof.(getproperty.(PhasePartition.(ts), :tot)) == typeof.(e_int)
    end

end

@testset "moist thermodynamics - dry limit" begin

    FT = Float64
    z, e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice =
        tested_profiles(param_set, 50, FT)

    # PhasePartition test is noisy, so do this only once:
    ts_dry = PhaseDry(param_set, first(e_int), first(ρ))
    ts_eq = PhaseEquil(param_set, first(e_int), first(ρ), typeof(first(ρ))(0))
    @test PhasePartition(ts_eq).tot ≈ PhasePartition(ts_dry).tot
    @test PhasePartition(ts_eq).liq ≈ PhasePartition(ts_dry).liq
    @test PhasePartition(ts_eq).ice ≈ PhasePartition(ts_dry).ice

    ts_dry = PhaseDry.(Ref(param_set), e_int, ρ)
    ts_eq = PhaseEquil.(Ref(param_set), e_int, ρ, q_tot .* 0)

    @test all(gas_constant_air.(ts_eq) .≈ gas_constant_air.(ts_dry))
    @test all(relative_humidity.(ts_eq) .≈ relative_humidity.(ts_dry))
    @test all(air_pressure.(ts_eq) .≈ air_pressure.(ts_dry))
    @test all(air_density.(ts_eq) .≈ air_density.(ts_dry))
    @test all(specific_volume.(ts_eq) .≈ specific_volume.(ts_dry))
    @test all(
        total_specific_humidity.(ts_eq) .≈ total_specific_humidity.(ts_dry),
    )
    @test all(cp_m.(ts_eq) .≈ cp_m.(ts_dry))
    @test all(cv_m.(ts_eq) .≈ cv_m.(ts_dry))
    @test all(air_temperature.(ts_eq) .≈ air_temperature.(ts_dry))
    @test all(internal_energy.(ts_eq) .≈ internal_energy.(ts_dry))
    @test all(internal_energy_sat.(ts_eq) .≈ internal_energy_sat.(ts_dry))
    @test all(soundspeed_air.(ts_eq) .≈ soundspeed_air.(ts_dry))
    @test all(latent_heat_vapor.(ts_eq) .≈ latent_heat_vapor.(ts_dry))
    @test all(latent_heat_sublim.(ts_eq) .≈ latent_heat_sublim.(ts_dry))
    @test all(latent_heat_fusion.(ts_eq) .≈ latent_heat_fusion.(ts_dry))
    @test all(q_vap_saturation.(ts_eq) .≈ q_vap_saturation.(ts_dry))
    @test all(saturation_excess.(ts_eq) .≈ saturation_excess.(ts_dry))
    @test all(liquid_fraction.(ts_eq) .≈ liquid_fraction.(ts_dry))
    @test all(liquid_ice_pottemp.(ts_eq) .≈ liquid_ice_pottemp.(ts_dry))
    @test all(dry_pottemp.(ts_eq) .≈ dry_pottemp.(ts_dry))
    @test all(virtual_pottemp.(ts_eq) .≈ virtual_pottemp.(ts_dry))
    @test all(liquid_ice_pottemp_sat.(ts_eq) .≈ liquid_ice_pottemp_sat.(ts_dry))
    @test all(exner.(ts_eq) .≈ exner.(ts_dry))

    @test all(
        saturation_vapor_pressure.(ts_eq, Ref(Ice())) .≈
        saturation_vapor_pressure.(ts_dry, Ref(Ice())),
    )
    @test all(
        saturation_vapor_pressure.(ts_eq, Ref(Liquid())) .≈
        saturation_vapor_pressure.(ts_dry, Ref(Liquid())),
    )
    @test all(first.(gas_constants.(ts_eq)) ≈ first.(gas_constants.(ts_dry)))
    @test all(last.(gas_constants.(ts_eq)) ≈ last.(gas_constants.(ts_dry)))

end


text:
struct TypeGenerationException <: Exception
	func::Symbol
	msg::AbstractString
	caught::Union{Exception, Nothing}
	function TypeGenerationException(func::Symbol, msg::AbstractString, caught::Union{Exception, Nothing}=nothing)
		new(func, msg, caught)
	end
end

# non-abstract types that the generator can build special rules for generating values
# is the same as DataGenerators.GENERATOR_SUPPORTED_CHOOSE_TYPES, but redefined here to avoid circular module dependencies
const DATAGENERATORS_SUPPORTED_TYPES = Type[Bool; Int8; Int16; Int32; Int64; UInt8; UInt16; UInt32; UInt64; Float16; Float32; Float64; String;]
# non-abstract types that the type translator can handle directly, either via Data Generator or itself
const DIRECTLY_SUPPORTED_TYPES = Type[DATAGENERATORS_SUPPORTED_TYPES; Tuple; DataType; Union;]

include("type_utilities.jl")
include("type_parse.jl")
include("type_transform.jl")
include("type_build.jl")
include("type_generator_utilities.jl")

function type_rules(t::Type, supporteddts::Vector{DataType}=Vector{DataType}(), rulenameprefix="")
	ast = parse_type(t, supporteddts)
	transform_type_ast(ast)
	transform_ast(ast) # this standard transform (which analysis reachability) isn't really needed for a type, but included for consistency with other translators
	build_type_rules(ast, rulenameprefix)
end

function type_generator(io::IO, genname::Symbol, t::Type, supporteddts::Vector{DataType}=Vector{DataType}())
	rules = type_rules(t, supporteddts)
	description = "an instance of type $(t)"
	output_generator(io, genname, description, rules)
end

type_generator(genname::Symbol, t::Type, supporteddts::Vector{DataType}=Vector{DataType}()) = include_generator(genname, type_generator, t, supporteddts)


text:
"""
    eom_rossler_lorenz() -> Function

Equations of motions for a 6D unidirectionally coupled Rössler-Lorenz system.
The 3D Rössler attractor is coupled to  a 3D Lorenz attractor such that the
second component  (``x_2``) of the Rössler system unidirectionally forces
the second component (``y_2``) of the Lorenz system. This example is from [1].

The implementation here allows for tuning the parameters of each subsystem by
introducing the constants `a₁`, `a₂`, `a₃`, `b₁`, `b₂`, `b₃`. The default
values of these constants are as defined in Krakovská et al. (2018).

The dynamics is generated by the vector field:

```math
\\begin{aligned}
\\dot x_1 &= a_1(x_2 + x_3) \\
\\dot x_2 &= a_2(x_1 + 0.2x_2) \\
\\dot x_3 &= a_2(0.2 + x_3(x_1 - a_3)) \\
\\dot y_1 &= b_1(y_2 - y_1) \\
\\dot y_2 &= y_1(b_2 - y_3) - y_2 +c(x_2)^2 \\
\\dot y_3 &= y_1y_2 - b_3y_3
\\end{aligned}
```

with the coupling constant ``c \\geq 0``.

# References
[Krakovská, Anna, et al. "Comparison of six methods for the detection of causality in a bivariate time series." Physical Review E 97.4 (2018):042207](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.97.042207)
"""
function eom_rossler_lorenz(u, p, t)
    a₁, a₂, a₃, b₁, b₂, b₃, c = (p...,)
    x₁, x₂, x₃, y₁, y₂, y₃ = (u...,)

    dx₁ = a₁*(x₂ + x₃)
    dx₂ = a₂*(x₁ + 0.2*x₂)
    dx₃ = a₂*(0.2 + x₃*(x₁ - a₃))
    dy₁ = b₁*(y₂ - y₁)
    dy₂ = y₁*(b₂ - y₃) - y₂ + c*x₂^2
    dy₃ = y₁*y₂ - b₃*y₃
    return SVector{6}(dx₁, dx₂, dx₃, dy₁, dy₂, dy₃)
end

function rossler_lorenz(u₀, a₁, a₂, a₃, b₁, b₂, b₃, c)
    p = [a₁, a₂, a₃, b₁, b₂, b₃, c]
    ContinuousDynamicalSystem(eom_rossler_lorenz, u₀, p)
end

"""
    rossler_lorenz(;u₀ = rand(6), a₁ = -6, a₂ = 6, a₃ = 2.0,
        b₁ = 10, b₂ = 28, b₃ = 8/3, c = 2) -> ContinuousDynamicalSystem

Initialise a Rössler-Lorenz system consisting of two independent 3D subsystems:
one Rössler attractor and one Lorenz attractor. They are coupled such that the
second component (`x₂`) of the Rössler system unidirectionally forces the
second component (`y₂`) of the Lorenz system.

The implementation here allows for tuning the parameters of each subsystem by
introducing the constants `a₁`, `a₂`, `a₃`, `b₁`, `b₂`, `b₃`. The default
values are as in Krakovská et al. (2018).

The dynamics is generated by the vector field:

```math
\\begin{aligned}
\\dot x_1 &= a_1(x_2 + x_3) \\
\\dot x_2 &= a_2(x_1 + 0.2x_2) \\
\\dot x_3 &= a_2(0.2 + x_3(x_1 - a_3)) \\
\\dot y_1 &= b_1(y_2 - y_1) \\
\\dot y_2 &= y_1(b_2 - y_3) - y_2 +c(x_2)^2 \\
\\dot y_3 &= y_1 y_2 - b_3y_3
\\end{aligned}
```

with the coupling constant ``c \\geq 0``.

# References
[Krakovská, Anna, et al. "Comparison of six methods for the detection of causality in a bivariate time series." Physical Review E 97.4 (2018):042207](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.97.042207)
"""
rossler_lorenz(;u₀ = rand(6), a₁ = -6, a₂ = 6, a₃ = 2.0, b₁ = 10, b₂ = 28, b₃ = 8/3, c = 2) =
    rossler_lorenz(u₀, a₁, a₂, a₃, b₁, b₂, b₃, c)


text:
include("base.jl")
include("circular_array_buffer.jl")
include("reservoir_array_buffer.jl")
include("device.jl")


text:
using Plots, JLD, Statistics, LaTeXStrings,StatsBase
number_list=[i for i in 0:9]
NList=[n for n in 100:100:10000]
STDList=[]
final_res=[]
for N in NList
    res=rand(number_list, N)
    push!(final_res, res)
    frequency_list=[]
    for i in number_list
        push!(frequency_list, length(findall(x->x==i, res)))
    end
    push!(STDList,std(frequency_list))
end
function Linear_fit(X, Y)
    A = [hcat(X) reshape(ones(length(X)), length(X), 1)]
    b = reshape(hcat(Y), length(Y), 1)
    line = (A \ b)
    return line
end
Line=Linear_fit(log.(NList), log.(STDList))
X = log.(NList)
Y = X .* Line[1,1] .+ Line[2,1]
plot(log.(NList),log.(STDList), xlabel="log(N)", ylabel="log(sigma)",
    title=L"Sigma\ versus\ Number\ of\ Iterations(100<N<10000, with\ steps\ of\ 10)",label=nothing,dpi=400,titlefontsize=11)
plot!(X,Y, color=:black,label = L"Y = %$(round(Line[1],digits=3))\ X + %$(round(Line[2],digits=3))")
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter6\\6.1\\Figs\\Sigma_N1.png")

plot(log.(NList),log.(STDList), xlabel="log(N)", ylabel="log(sigma)",
    title=L"Sigma\ versus\ Number\ of\ Iterations(100<N<10000, with\ steps\ of\ 100)", label=nothing,dpi=400,titlefontsize=11)
plot!(X,Y, color=:black,label = L"Y = %$(round(Line[1],digits=3))\ X + %$(round(Line[2],digits=3))")
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter6\\6.1\\Figs\\Sigma_N2.png")