text:
# comparing the iSARA kernel to canonical kernels
using SARA
using SARA: iSARA, random_sampling
using CovarianceFunctions
const Kernel = CovarianceFunctions
using GaussianDistributions

################# load data
include("inner_load_data.jl")
using SARA: estimate_noise

################# setup kernels for comparison ##############
# θ are stripe-dependent rescaling parameters
# LSA-bias is encoded in last parameter tuple
# l is length scale for all kernels
function get_kernels(l, θ)
    k1 = Kernel.Lengthscale(Kernel.EQ(), l)
    k2 = Kernel.Lengthscale(Kernel.MaternP(2), l)
    k3 = iSARA(l)

    # only LSA-inspired rescaling
    k4 = iSARA(l, θ[end])  # last parameter tuple

    # with all rescaling parameters
    k5 = iSARA(l, θ)
    kernels = [k1, k2, k3, k4, k5]
    return kernels
end

############### generate sampling paths with different kernels
function sampling_path(policy, G::Gaussian, x::AbstractVector, y::AbstractVector, σ::Real, nsample::Int)
    C = G
    rmse, r2 = zeros(nsample), zeros(nsample)
    indices = zeros(Int, 0)
    sstot = sum(abs2, y .- mean(y)) # total sum of squares
    for i in 1:nsample # random sampling
        xi = policy(C, x)
        append!(indices, findfirst(==(xi), x))
        sort!(indices)
        C = G | @views (x[indices], y[indices], σ^2)
        μ = mean(C, x)
        δ = y - μ
        rmse[i] = sqrt(mean(δ.^2)) # RMSD (root mean squared deviation)
        r2[i] = 1 - (sum(abs2, δ) / sstot) # coefficient of determination
    end
    return rmse, r2
end

# on a single stripe using random sampling
function kernel_benchmark!(rmse, r2, kernels, x, y, σ)
    nsample, nexp = size(rmse)[1:2]
    for (j, k) in enumerate(kernels)
        G = Gaussian(k)
        Base.Threads.@threads for i in 1:nexp
            d, r = sampling_path(random_sampling, G, x, y, σ, nsample)
            rmse[:, i, j] = d
            r2[:, i, j] = r
        end
    end
    return rmse
end

function kernel_benchmark_driver(position, optical, rescaling_parameters,
                                 σ::Real, y_normalization::Real)
    nstripes = length(position)
    nkernel = 5
    nsample = 64
    lengthscales = collect(range(.05, .25, step = .05))
    nlength = length(lengthscales)
    nexp = 128 # number of independent experiments per stripe
    rmse = zeros(nsample, nexp, nkernel, nlength, nstripes)
    r2 = zeros(nsample, nexp, nkernel, nlength, nstripes)
    for i in 1:nstripes
        println("stripe $i of $nstripes")
        θ = rescaling_parameters[i]
        for (j, l) in enumerate(lengthscales)
            kernels = get_kernels(l, θ)

            rmse_i = @view rmse[:, :, :, j, i]
            r2_i = @view r2[:, :, :, j, i]

            x = position[i]
            y = optical[i][1] # for simplicity, benchmark with first and most significant optical coefficient
            y .-= mean(y) # normalized and centered
            y /= y_normalization
            kernel_benchmark!(rmse_i, r2_i, kernels, x, y, σ)
        end
    end

    if true
        f = h5open("inner_kernel_benchmark.h5", "w")
        f["rmse"] = rmse
        f["r2"] = r2
        f["lengthscales"] = lengthscales
        f["nexp"] = nexp
        f["nsample"] = nsample
        f["nlength"] = nlength
        f["nstripes"] = nstripes
        f["nkernel"] = nkernel
        f["kernel names"] = ["EQ", "Matern_2", "iSARA", "iSARA+LSA", "iSARA+RGB"]
        close(f)
    end
    return rmse, r2
end

# path = "/Users/sebastianament/Documents/SEA/XRD Analysis/SARA/Bi2O3_19F44_01/"
path = ""
file = "Bi2O3_19F44_01_inner_loop_data.json"
f = JSON.parsefile(path * file)
position, optical, rescaling_parameters, T_max, log10_τ = load_data(f)

# estimating standard deviation of noise of normalized coefficient
x = position[1]
y = optical[1][1]
y .-= mean(y)
y_normalization = std(y)
# y_normalization = maximum(abs, y)
y /= y_normalization
σ = estimate_noise(x, y)
rmse, r2 = kernel_benchmark_driver(position, optical, rescaling_parameters, σ, y_normalization)

# mrmse = reshape(mean(rmse, dims = (2, 5)), 64, 5)
# mr2 = reshape(mean(r2, dims = (2, 5)), 64, 5)

# c = [x[1] for x in optical]


text:
# https://adventofcode.com/2021/day/2
using AdventOfCode

input = readlines("2021/data/day_2.txt")

function part_1(input)
    h,d=0,0
    for command in input
        X = parse(Int, match(r"[a-z]* (?<num>\d+)", command)[:num])
        if occursin(r"forward.*",command) h+=X end
        if occursin(r"down.*",command) d+=X end
        if occursin(r"up.*",command) d-=X end
    end
    h*d
end

function part_2(input)
    h,d,aim=0,0,0
    for command in input
        X = parse(Int, match(r"[a-z]* (?<num>\d+)", command)[:num])
        if occursin(r"forward.*",command) h+=X;d+=X*aim end
        if occursin(r"down.*",command) aim+=X end
        if occursin(r"up.*",command) aim-=X end
    end
    h*d
end

ext: jl


text:
function *(x::Sonum16,y::Sonum16)
    if signbit(x)
        if signbit(y)
            @inbounds return TableMul16S[Int(-x)+1,Int(-y)+1]
        else
            @inbounds return -TableMul16S[Int(-x)+1,Int(y)+1]
        end
    else
        if signbit(y)
            @inbounds return -TableMul16S[Int(x)+1,Int(-y)+1]
        else
            @inbounds return TableMul16S[Int(x)+1,Int(y)+1]
        end
    end
end

function +(x::Sonum16,y::Sonum16)
    if signbit(x)
        if signbit(y)   # -a-b = -(a+b)
            @inbounds return -TableAdd16S[Int(-x)+1,Int(-y)+1]
        else            # -a+b = b-a

            # anti-symmetric: check for size
            a = Int(y)+1
            b = Int(-x)+1

            if a > b
                @inbounds return -TableSub16S[b,a]
            else
                @inbounds return TableSub16S[a,b]
            end

        end
    else
        if signbit(y)   # a-b
            #return TableSub16[Int(x)+1,Int(-y)+1]

            # anti-symmetric: check for size
            a = Int(x)+1
            b = Int(-y)+1

            if a > b
                @inbounds return -TableSub16S[b,a]
            else
                @inbounds return TableSub16S[a,b]
            end
        else
            @inbounds return TableAdd16S[Int(x)+1,Int(y)+1]
        end
    end
end

function -(x::Sonum16,y::Sonum16)
    if signbit(x)
        if signbit(y)   # -a--b = b-a
            # return TableSub16[Int(-y)+1,Int(-x)+1]

            # anti-symmetric: check for size
            a = Int(-y)+1
            b = Int(-x)+1

            if a > b
                @inbounds return -TableSub16S[b,a]
            else
                @inbounds return TableSub16S[a,b]
            end

        else            # -a-b = -(a+b)
            @inbounds return -TableAdd16S[Int(-x)+1,Int(y)+1]
        end
    else
        if signbit(y)   # a--b = a+b
            @inbounds return TableAdd16S[Int(x)+1,Int(-y)+1]
        else            # a-b
            # return TableSub16[Int(x)+1,Int(y)+1]

            # anti-symmetric: check for size
            a = Int(x)+1
            b = Int(y)+1

            if a > b
                @inbounds return -TableSub16S[b,a]
            else
                @inbounds return TableSub16S[a,b]
            end
        end
    end
end

function /(x::Sonum16,y::Sonum16)
    if signbit(x)
        if signbit(y)
            @inbounds return TableDiv16[Int(-x)+1,Int(-y)+1]
        else
            @inbounds return -TableDiv16[Int(-x)+1,Int(y)+1]
        end
    else
        if signbit(y)
            @inbounds return -TableDiv16[Int(x)+1,Int(-y)+1]
        else
            @inbounds return TableDiv16[Int(x)+1,Int(y)+1]
        end
    end
end

# This increases the error as two roundings are applied
# function /(x::Sonum16,y::Sonum16)
#     return x*inv(y)
# end

function sqrt(x::Sonum16)
    if signbit(x)
        return notareal(Sonum16)
    else
        @inbounds return ListSqrt16[Int(x)+1]
    end
end

function inv(x::Sonum16)
    if signbit(x)
        @inbounds return -ListInv16[Int(-x)+1]
    else
        @inbounds return ListInv16[Int(x)+1]
    end
end

# Table-less functions via Float64 conversion - slow but memory efficient
# function *(x::Sonum16,y::Sonum16)
#     return Sonum16(Float64(x)*Float64(y))
# end
#
# function +(x::Sonum16,y::Sonum16)
#     return Sonum16(Float64(x)+Float64(y))
# end
#
# function -(x::Sonum16,y::Sonum16)
#     return Sonum16(Float64(x)-Float64(y))
# end
#
# function /(x::Sonum16,y::Sonum16)
#     return Sonum16(Float64(x)/Float64(y))
# end


text:
# """
#     MCSampling

# Sampling method for the Monte Carlo integration.
# """
abstract type MCSampling{T} end
Base.eltype(::MCSampling{T}) where T = eltype(T)

"""
    UniformSampling(a, b)

Uniform sampling for the Monte Carlo integration, in the hypercube `[a, b]^2`.
"""
struct UniformSampling{A} <: MCSampling{A}
    a::A
    b::A
end
@functor UniformSampling

function (mc_sample::UniformSampling{T})(x_mc, kwargs...) where T
    Tel = eltype(T)
    rand!(x_mc)
    m = (mc_sample.b + mc_sample.a) ./ convert(Tel,2)
    x_mc .= (x_mc .- convert(Tel,0.5)) .* (mc_sample.b - mc_sample.a) .+ m
end

"""
    NormalSampling(σ)
    NormalSampling(σ, shifted)

Normal sampling method for the Monte Carlo integration.

# Arguments
* `σ`: the standard devation of the sampling
* `shifted` : if true, the integration is shifted by `x`. Defaults to false.
"""
struct NormalSampling{T} <: MCSampling{T}
    σ::T
    shifted::Bool # if true, we shift integration by x when invoking mc_sample::MCSampling(x)
end
@functor NormalSampling

NormalSampling(σ) = NormalSampling(σ,false)

function (mc_sample::NormalSampling)(x_mc)
    randn!(x_mc)
    x_mc .*=  mc_sample.σ
end

function (mc_sample::NormalSampling)(x_mc, x)
    mc_sample(x_mc)
    mc_sample.shifted ? x_mc .+= x : nothing
end

struct NoSampling <: MCSampling{Nothing} end

(mc_sample::NoSampling)(x...) = nothing

function _integrate(::MCS) where {MCS <: MCSampling}
    if MCS <: NoSampling
        return false
    else
        return true
    end
end


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

# Trimmed-down code from Pkg2 with only the essential parts needed by loadmeta.jl

module Pkg2

# DIR
const DIR_NAME = ".julia"
_pkgroot() = abspath(get(ENV,"JULIA_PKGDIR",joinpath(homedir(),DIR_NAME)))
isversioned(p::AbstractString) = ((x,y) = (VERSION.major, VERSION.minor); basename(p) == "v$x.$y")

function dir()
    b = _pkgroot()
    x, y = VERSION.major, VERSION.minor
    d = joinpath(b,"v$x.$y")
    if isdir(d) || !isdir(b) || !isdir(joinpath(b, "METADATA"))
        return d
    end
    return b
end
dir(pkg::AbstractString...) = normpath(dir(),pkg...)

include("types.jl")
include("reqs.jl")

end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
function Experiment(
    ::Val{:JuliaRL},
    ::Val{:IQN},
    ::Val{:CartPole},
    ::Nothing;
    save_dir = nothing,
    seed = 123,
)
    if isnothing(save_dir)
        t = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
        save_dir = joinpath(pwd(), "checkpoints", "JuliaRL_IQN_CartPole_$(t)")
    end

    lg = TBLogger(joinpath(save_dir, "tb_log"), min_level = Logging.Info)
    rng = StableRNG(seed)

    env = CartPoleEnv(; T = Float32, rng = rng)
    ns, na = length(state(env)), length(action_space(env))

    init = glorot_uniform(rng)
    Nₑₘ = 16
    n_hidden = 64
    κ = 1.0f0

    nn_creator() =
        ImplicitQuantileNet(
            ψ = Dense(ns, n_hidden, relu; init = init),
            ϕ = Dense(Nₑₘ, n_hidden, relu; init = init),
            header = Dense(n_hidden, na; init = init),
        ) |> cpu

    agent = Agent(
        policy = QBasedPolicy(
            learner = IQNLearner(
                approximator = NeuralNetworkApproximator(
                    model = nn_creator(),
                    optimizer = ADAM(0.001),
                ),
                target_approximator = NeuralNetworkApproximator(model = nn_creator()),
                κ = κ,
                N = 8,
                N′ = 8,
                Nₑₘ = Nₑₘ,
                K = 32,
                γ = 0.99f0,
                stack_size = nothing,
                batch_size = 32,
                update_horizon = 1,
                min_replay_history = 100,
                update_freq = 1,
                target_update_freq = 100,
                default_priority = 1.0f2,
                rng = rng,
                device_rng = rng,
            ),
            explorer = EpsilonGreedyExplorer(
                kind = :exp,
                ϵ_stable = 0.01,
                decay_steps = 500,
                rng = rng,
            ),
        ),
        trajectory = CircularArrayPSARTTrajectory(
            capacity = 1000,
            state = Vector{Float32} => (ns,),
        ),
    )

    stop_condition = StopAfterStep(10_000)

    total_reward_per_episode = TotalRewardPerEpisode()
    time_per_step = TimePerStep()
    hook = ComposedHook(
        total_reward_per_episode,
        time_per_step,
        DoEveryNStep() do t, agent, env
            if agent.policy.learner.update_step % agent.policy.learner.update_freq == 0
                with_logger(lg) do
                    @info "training" loss = agent.policy.learner.loss
                end
            end
        end,
        DoEveryNEpisode() do t, agent, env
            with_logger(lg) do
                @info "training" reward = total_reward_per_episode.rewards[end] log_step_increment =
                    0
            end
        end,
    )

    description = """
    This experiment uses the `IQNLearner` method with a `ImplicitQuantileNet`.
    The testing environment is CartPoleEnv.
    """

    Experiment(agent, env, stop_condition, hook, description)
end


text:
"""
    buildcircuit(M::Union{MPS,MPO}, gates::Vector{<:Tuple};
                 noise = nothing)

Generates a vector of (gate) `ITensor`, from a vector of `Tuple`
associated with a list of quantum gates.
If noise is nontrivial, the corresponding Kraus operators are
added to each gate as a tensor with an extra (Kraus) index.
"""
function buildcircuit(
  M::Union{MPS,MPO,ITensor},
  circuit::Union{Tuple,Vector{<:Any}};
  noise::Union{Nothing, Tuple, NamedTuple} = nothing
)
  circuit_tensors = ITensor[]
  if circuit isa Tuple
    circuit = [circuit]
  end
  if !isnothing(noise)
    circuit = insertnoise(circuit, noise)
  end
  for g in circuit
    push!(circuit_tensors, gate(M, g))
  end
  return circuit_tensors
end

function buildcircuit(M::Union{MPS,MPO,ITensor}, circuit::Vector{Vector{<:Any}}; kwargs...)
  return buildcircuit(M, vcat(circuit...); kwargs...)
end

"""
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
                            TENSOR NETWORK SIMULATOR
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
"""

"""
    runcircuit(M::Union{MPS,MPO}, gate_tensors::Vector{<:ITensor};
               kwargs...)

Apply the circuit to a state (wavefunction/densitymatrix) from a list of tensors.
"""
function runcircuit(
  M::Union{MPS,MPO},
  circuit_tensors::Vector{<:ITensor};
  apply_dag=nothing,
  cutoff=1e-15,
  maxdim=10_000,
  svd_alg="divide_and_conquer",
  move_sites_back::Bool=true,
  kwargs...,
)

  # Check if gate_tensors contains Kraus operators
  inds_sizes = [length(inds(g)) for g in circuit_tensors]
  noiseflag = any(x -> x % 2 == 1, inds_sizes)

  if apply_dag == false && noiseflag == true
    error("noise simulation requires apply_dag=true")
  end

  # Default mode (apply_dag = nothing)
  if isnothing(apply_dag)
    # Noisy evolution: MPS/MPO -> MPO
    if noiseflag
      # If M is an MPS, |ψ⟩ -> ρ = |ψ⟩⟨ψ| (MPS -> MPO)
      ρ = (typeof(M) == MPS ? MPO(M) : M)
      # ρ -> ε(ρ) (MPO -> MPO, conjugate evolution)
      return apply(
        circuit_tensors,
        ρ;
        apply_dag=true,
        cutoff=cutoff,
        maxdim=maxdim,
        svd_alg=svd_alg,
        move_sites_back=move_sites_back,
      )

      # Pure state evolution
    else
      # |ψ⟩ -> U |ψ⟩ (MPS -> MPS)
      #  ρ  -> U ρ U† (MPO -> MPO, conjugate evolution)
      if M isa MPS
        return apply(
          circuit_tensors,
          M;
          cutoff=cutoff,
          maxdim=maxdim,
          svd_alg=svd_alg,
          move_sites_back=move_sites_back,
        )
      else
        return apply(
          circuit_tensors,
          M;
          apply_dag=true,
          cutoff=cutoff,
          maxdim=maxdim,
          svd_alg=svd_alg,
          move_sites_back=move_sites_back,
        )
      end
    end
    # Custom mode (apply_dag = true / false)
  else
    if M isa MPO
      # apply_dag = true:  ρ -> U ρ U† (MPO -> MPO, conjugate evolution)
      # apply_dag = false: ρ -> U ρ (MPO -> MPO)
      return apply(
        circuit_tensors,
        M;
        apply_dag=apply_dag,
        cutoff=cutoff,
        maxdim=maxdim,
        svd_alg=svd_alg,
        move_sites_back=move_sites_back,
      )
    elseif M isa MPS
      # apply_dag = true:  ψ -> U ψ -> ρ = (U ψ) (ψ† U†) (MPS -> MPO, conjugate)
      # apply_dag = false: ψ -> U ψ (MPS -> MPS)
      Mc = apply(
        circuit_tensors,
        M;
        cutoff=cutoff,
        maxdim=maxdim,
        svd_alg=svd_alg,
        move_sites_back=move_sites_back,
      )
      if apply_dag
        Mc = MPO(Mc)
      end
      return Mc
    else
      error("Input state must be an MPS or an MPO")
    end
  end
end

"""
     runcircuit(M::Union{MPS, MPO}, gates::Union{Tuple,Vector{<:Any}};
                noise = nothing, kwargs...)

    runcircuit(M::Union{MPS, MPO}, gates::Vector{<:Vector{<:Any}};
               observer! = nothing,
               move_sites_back_before_measurements::Bool = false,
               noise = nothing,
               kwargs...)

Apply a quantum circuit to a state (wavefunction or density matrix)

If an MPS `|ψ⟩` is input, there are three possible modes:

1. By default (`noise = nothing` and `apply_dag = nothing`), the evolution `U|ψ⟩` is performed.
2. If `noise` is set to something nontrivial, the mixed evolution `ε(|ψ⟩⟨ψ|)` is performed.
   Example: `noise = ("amplitude_damping", (γ = 0.1,))` (amplitude damping channel with decay rate `γ = 0.1`)
3. If `noise = nothing` and `apply_dag = true`, the evolution `U|ψ⟩⟨ψ|U†` is performed.

If an MPO `ρ` is input, there are three possible modes:

1. By default (`noise = nothing` and `apply_dag = nothing`), the evolution `U ρ U†` is performed.
2. If `noise` is set to something nontrivial, the evolution `ε(ρ)` is performed.
3. If `noise = nothing` and `apply_dag = false`, the evolution `Uρ` is performed.

If an `Observer` is provided as input, and `circuit` is made out of a sequence
of layers of gates, performs a measurement of the observables contained in
Observer, after the application of each layer.
"""
runcircuit(M::Union{MPS, MPO,ITensor}, circuit::Union{Tuple,Vector{<:Any}};
           noise = nothing, kwargs...) =
  runcircuit(M, buildcircuit(M, circuit; noise = noise); kwargs...)

function runcircuit(
  M::Union{MPS,MPO},
  circuit::Vector{<:Vector{<:Any}};
  (observer!)=nothing,
  move_sites_back_before_measurements::Bool=false,
  noise = nothing,
  outputlevel = 1,
  outputpath = nothing,
  savestate = false,
  print_metrics = [],
  kwargs...,
)

  # is the observer is not provided, apply the vectorized circuit
  isnothing(observer!) && return runcircuit(M, vcat(circuit...); noise=noise, kwargs...)

  # issue warning if there are custom functions and the sites are not being moved back
  if !isnothing(observer!) && move_sites_back_before_measurements == false
    println("--------------")
    println(" WARNING")
    println(
      "\nA custom function is being measured during the gate evolution,\nbut the MPS sites at depth D are not being restored to their \nlocation at depth D-1. If any custom measurement is dependent \non the specific site ordering of the MPS, it should take the \nre-ordering into account. The MPS can be restored to the \ncanonical index ordering by setting:\n\n `move_sites_back_before_measurements = true`\n",
    )
    println("--------------\n")
  end
  M0 = copy(M)
  # record the initial configuration of the indices
  s = siteinds(M)
  for l in 1:length(circuit)
    layer = circuit[l]
    t = @elapsed begin
      M = runcircuit(M, layer; noise = noise,
                     move_sites_back=move_sites_back_before_measurements,
                     kwargs...)
      if !isnothing(observer!)
        update!(observer!, M; sites = s)
      end
    end
    if outputlevel ≥ 1
      @printf("%-4d  ", l)
      @printf("χ = %-4d  ", maxlinkdim(M))
      #TODO add the truncation error here
      printobserver(observer!, print_metrics)
      @printf("elapsed = %-4.3fs", t)
      println()
    end
    if !isnothing(outputpath)
      observerpath = outputpath * "_observer.jld2"
      save(observerpath, observer!)
      if savestate
        statepath = outputpath * "_state.h5"
        h5rewrite(statepath) do fout
          write(fout, "state", M)
        end
      end
    end
  end
  if move_sites_back_before_measurements == false
    new_s = siteinds(M)
    ns = 1:length(M)
    ñs = [findsite(M0, i) for i in new_s]
    M = movesites(M, ns .=> ñs; kwargs...)
  end
  return M
end

"""
    runcircuit(N::Int, gates::Vector{<:Tuple};
               process = false,
               noise = nothing,
               cutoff = 1e-15,
               maxdim = 10000,
               svd_alg = "divide_and_conquer")

Run the circuit corresponding to a list of quantum gates on a system of `N` qubits.
The starting state is generated automatically based on the flags `process`, `noise`, and `apply_dag`.

1. By default (`noise = nothing`, `apply_dag = nothing`, and `process = false`),
   the evolution `U|ψ⟩` is performed where the starting state is set to `|ψ⟩ = |000...⟩`.
   The MPS `U|000...⟩` is returned.
2. If `noise` is set to something nontrivial, the mixed evolution `ε(|ψ⟩⟨ψ|)` is performed,
   where the starting state is set to `|ψ⟩ = |000...⟩`.
   The MPO `ε(|000...⟩⟨000...|)` is returned.
3. If `process = true` and `noise = nothing`, the evolution `U 1̂` is performed,
   where the starting state `1̂ = (1⊗1⊗1⊗…⊗1)`. The MPO approximation for the unitary
   represented by the set of gates is returned.
4. If `process = true` and `noise` is set to something nontrivial, the function returns the Choi matrix
   `Λ = ε⊗1̂(|ξ⟩⟨ξ|)`, where `|ξ⟩= ⨂ⱼ |00⟩ⱼ+|11⟩ⱼ`, approximated by a MPO with 4 site indices,
   two for the input and two for the output Hilbert space of the quantum channel.
"""
function runcircuit(sites::Vector{<:Index},
                    circuit::Union{Tuple,Vector{<:Any},Vector{Vector{<:Any}}};
                    process::Bool = false, noise = nothing,
                    full_representation::Bool = false,
                    kwargs...)

  # Unitary operator for the circuit
  if process && isnothing(noise)
    U₀ = full_representation ? prod(productoperator(sites)) : productoperator(sites)
    return runcircuit(U₀, circuit;
                      noise = nothing,
                      apply_dag = false,
                      kwargs...)
  end
  # Choi matrix
  if process && !isnothing(noise)
    return choimatrix(sites, circuit;
                      noise = noise,
                      full_representation = full_representation,
                      kwargs...)
  end

  M₀ = full_representation ? prod(productstate(sites)) : productstate(sites)
  return runcircuit(M₀, circuit; noise = noise, kwargs...)
end

runcircuit(N::Int, circuit::Any; kwargs...) =
  runcircuit(siteinds("Qubit", N), circuit; kwargs...)

runcircuit(circuit::Any; kwargs...) =
  runcircuit(nqubits(circuit), circuit; kwargs...)

"""
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
                               EXACT SIMULATOR
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
"""

"""
    runcircuit(M::ITensor,gate_tensors::Vector{ <: ITensor}; kwargs...)

Apply the circuit to a ITensor from a list of tensors.
"""

function runcircuit(T::ITensor, circuit_tensors::Vector{<:ITensor}; apply_dag = nothing, kwargs...)
  # Check if gate_tensors contains Kraus operators
  inds_sizes = [length(inds(g)) for g in circuit_tensors]
  noiseflag = any(x -> x % 2 == 1 , inds_sizes)

  if apply_dag == false && noiseflag == true
    error("noise simulation requires apply_dag=true")
  end
  # Default mode (apply_dag = nothing)
  if isnothing(apply_dag)
    # Noisy evolution: MPS/MPO -> MPO
    if noiseflag
      T = is_operator(T) ? T : prime(T, "Site") * dag(T)
      # ρ -> ε(ρ) (MPO -> MPO, conjugate evolution)
      return apply(circuit_tensors, T; apply_dag = true)
    # Pure state evolution
    else
      # |ψ⟩ -> U |ψ⟩ (MPS -> MPS)
      #  ρ  -> U ρ U† (MPO -> MPO, conjugate evolution)
      if !is_operator(T)
        return apply(circuit_tensors, T)
      else
        return apply(circuit_tensors, T; apply_dag = true)
      end
    end
  # Custom mode (apply_dag = true / false)
  else
    Tc = apply(circuit_tensors, T; apply_dag = apply_dag)
    !is_operator(T) && apply_dag && return prime(Tc, "Site") * dag(Tc)
    return Tc
  end
end

"""
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
                                  CHOI MATRIX
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
"""

"""

Compute the Choi matrix `Λ  = ε ⊗ 1̂(|ξ⟩⟨ξ|)`, where `|ξ⟩= ⨂ⱼ |00⟩ⱼ+|11⟩ⱼ`,
and `ε`` is a quantum channel built out of a set of quantum gates and
a local noise model. Returns a MPO with `N` tensor having 4 sites indices.
"""
choimatrix(circuit::Vector{<:Any}; kwargs...) =
  choimatrix(nqubits(circuit), circuit; kwargs...)

choimatrix(sites::Union{Int, Vector{<:Index}}, args...; kwargs...) =
  choimatrix(unitary_mpo_to_choi_mpo(productoperator(sites)), args...; kwargs...)

function choimatrix(Λ0::MPO, circuit::Vector{<:Any};
                    noise = nothing, cutoff = 1e-15, maxdim = 10000,
                    svd_alg = "divide_and_conquer", full_representation = false)
  N = length(Λ0)
  if isnothing(noise)
    error("choi matrix requires noise")
  end
  # if circuit has layer structure, transform to vectorized circuit
  circuit = (circuit isa Vector{<:Vector} ? vcat(circuit...) : circuit)

  # TODO: simplify by building the circuit directly from the Choi MPO.
  s = [firstind(Λ0[j], tags = "Output") for j in 1:length(Λ0)]
  compiler = productoperator(s)
  prime!(compiler,-1,tags = "Qubit")
  circuit_tensors = buildcircuit(compiler, circuit; noise = noise)

  # contract to compute the Choi matrix
  Λ₀ = full_representation ? prod(Λ0) : Λ0
  return runcircuit(Λ₀, circuit_tensors; apply_dag = true, cutoff = cutoff,
                    maxdim = maxdim, svd_alg = svd_alg)
end


text:
"""
    SingleLevel <: ERA5Variable

Abstract supertype for Single-Level variables.  Contains the following fields:
- `varID` : The variable ID, that is also the identifier in the NetCDF files
- `lname` : The variable long-name, which is used to specify retrievals from CDS
- `vname` : The full-name of the variable
- `units` : The units of the variable
"""
abstract type SingleLevel <: ERA5Variable end

"""
    SingleVariable <: SingleLevel

Subtype for Single-Level variables that can be directly retrieved from the CDS
"""
struct SingleVariable{ST<:AbstractString} <: SingleLevel
    varID :: ST
    lname :: ST
    vname :: ST
    units :: ST
end

"""
    SingleCustom <: SingleLevel

Subtype for custom user-defined Single-Level variables
"""
struct SingleCustom{ST<:AbstractString} <: SingleLevel
    varID :: ST
    lname :: ST
    vname :: ST
    units :: ST
end

"""
    SingleVariable(
        varID :: AbstractString,
        ST = String,
    ) -> evar :: SingleLevel

Retrieve the basic properties of the Single-Level variable defined by `varID` and put them in the `evar` SingleLevel type structure.

Arguments
=========

- `varID` : variable ID (in string format) used in the NetCDF file
"""
function SingleVariable(
    varID :: AbstractString,
    ST = String,
)

    isSingle(varID)

    @info "$(modulelog()) - Retrieving information for the SingleVariable defined by the ID \"$varID\""
    vlist,flist = listSingles();      ind = findall(varID.==vlist)[1]
    vtype = replace(flist[ind],".txt"=>"")
    fname = joinpath(DEPOT_PATH[1],"files","ERA5Reanalysis",flist[ind])
    vlist = listera5variables(fname); ind = findall(varID.==vlist)[1]

    IDinfo = readdlm(fname,',',comments=true,comment_char='#')[ind,:]
    varID,lname,vname,units = IDinfo[[1,2,3,4]]

    if vtype == "singlevariable"
          return SingleVariable{ST}(varID,lname,vname,units)
    else; return SingleCustom{ST}(varID,lname,vname,units)
    end

end

"""
    SingleVariable(
        ST = String;
        varID :: AbstractString,
        lname :: AbstractString = "",
        vname :: AbstractString,
        units :: AbstractString,
        inCDS :: Bool = true
    ) -> evar :: SingleLevel

Create a custom Single-Level variable that is not in the default list exported by ERA5Reanalysis.jl.  These variables are either available in the CDS store (whereby they can be both downloaded analyzed), or not (in which case means that they were separately calculated from other variables and analyzed).

Keyword Arguments
=================

- `varID` : variable ID (in string format) used in the NetCDF file
- `lname` : long-name for variable (used in specifying variable for CDS downloads)
- `vname` : user-defined variable name
- `units` : user-defined units of the variable
- `inCDS` : Boolean that indicates if this variable is available on the CDS store.  True if available.
"""
function SingleVariable(
    ST = String;
    varID :: AbstractString,
    lname :: AbstractString = "",
    vname :: AbstractString,
    units :: AbstractString,
    inCDS :: Bool = true
)

    if isSingle(varID,throw=false)
        error("$(modulelog()) - The SingleVariable \"$(varID)\" has already been defined,please use another identifier.")
    else
        @info "$(modulelog()) - Adding the SingleVariable \"$(varID)\" to the list."
    end

    if inCDS

        open(joinpath(
            DEPOT_PATH[1],"files","ERA5Reanalysis","singlevariable.txt"
        ),"a") do io
            write(io,"$varID,$lname,$vname,$units\n")
        end
        return SingleVariable{ST}(varID,lname,vname,units)

    else

        open(joinpath(DEPOT_PATH[1],"files","ERA5Reanalysis","singlecustom.txt"),"a") do io
            write(io,"$varID,$lname,$vname,$units\n")
        end
        return SingleCustom{ST}(varID,lname,vname,units)

    end

end

"""
    listSingles() -> varlist :: Array{ST}, fidlist :: Array{ST}

List all Single-Level Variables and the files the data are stored in.

Output
======

- `varlist` : List of all the Single-Level Variable IDs
- `fidlist` : List of the files that the Single-Level Variable information is stored in
"""
function listSingles()

    flist   = ["singlevariable.txt","singlecustom.txt"]
    varlist = []
    fidlist = []

    for fname in flist
        copyera5variables(fname)
        vvec = listera5variables(joinpath(DEPOT_PATH[1],"files","ERA5Reanalysis",fname))
        varlist = vcat(varlist,vvec); fvec = fill(fname,length(vvec))
        fidlist = vcat(fidlist,fvec)
    end

    return varlist,fidlist

end

"""
    isSingle(
        varID :: AbstractString;
        throw :: Bool = true
    ) -> tf :: Bool

Extracts information of the Single-Level Variable with the ID `varID`.  If no Single-Level Variable with this ID exists, an error is thrown.

Arguments
=========

- `RegID` : The keyword ID that will be used to identify the Single-Level Variable.
        If the ID is not valid (i.e. not being used), then an error will be thrown.
- `throw` : If `true`, then throws an error if `RegID` is not a valid Single-Level Variable identifier instead of returning the Boolean `tf`

Returns
=======

- `tf` : True / False
"""
function isSingle(
    varID :: AbstractString;
    throw :: Bool = true
)

    @info "$(modulelog()) - Checking to see if the SingleVariable ID \"$varID\" is in use"
    vlist,_ = listSingles()

    if sum(vlist.==varID) == 0
        if throw
            error("$(modulelog()) - \"$(varID)\" is not a valid SingleVariable identifier, use the function SingleVariable() or SingleCustom() to add this ERA5Variable to the list.")
        else
            @warn "$(modulelog()) - \"$(varID)\" is not a valid SingleVariable identifier, use the function SingleVariable() or SingleCustom() to add this ERA5Variable to the list."
            return false
        end
    else
        @info "$(modulelog()) - The SingleVariable ID \"$varID\" is already in use"
        return true
    end

end

"""
    rmSingle( varID :: AbstractString ) -> nothing

Remove the Single-Level Variable with the ID `varID` from the lists.

Arguments
=========

- `RegID` : The keyword ID that will be used to identify the Single-Level Variable that is to be removed
"""
function rmSingle(varID::AbstractString)

    if isSingle(varID,throw=false)
        disable_logging(Logging.Warn)
        rmERA5Variable(SingleVariable(varID))
        disable_logging(Logging.Debug)
        @info "$(modulelog()) - Successfully removed the Single-Level variable defined by \"$(varID)\""
    else
        @warn "$(modulelog()) - No Single-Level variable defined by \"$(varID)\" exists, please make sure you specified the correct variable ID"
    end

    return nothing

end

"""
    resetSingles( allfiles :: Bool ) -> nothing

Reset the list of Single-Level variables.

Arguments
=========

- `allfiles` : If false, only get rid of all the SingleCustom variables, but if true, then the SingleVariable list will be reset back to the default for ERA5Reanalysis.jl
"""
function resetSingles(;allfiles=false)

    if allfiles
        @info "$(modulelog()) - Resetting both the master and custom lists of ERA5 SingleVariables back to the default"
        flist = ["singlevariable.txt","singlecustom.txt"]
    else
        @info "$(modulelog()) - Resetting the custom lists of ERA5 variables back to the default"
        flist = ["singlecustom.txt"]
    end

    for fname in flist
        copyera5variables(fname,overwrite=true)
    end

    return nothing

end

function tableSingles()

    jfol = joinpath(DEPOT_PATH[1],"files","ERA5Reanalysis"); mkpath(jfol);
    fvar = ["SingleVariable","SingleCustom"]
    fmat = []

    for fname in fvar
        fid  = joinpath(jfol,"$(lowercase(fname)).txt")
        try
            vmat = readdlm(fid,',',comments=true,comment_char='#')
            nvar = size(vmat,1); ff = fill(fname,nvar)
            vmat = cat(ff,vmat[:,[1,3,4,2]],dims=2)
            fmat = cat(fmat,vmat,dims=1)
        catch
        end
    end

    head = ["Variable Type","ID","Name","Units","ERA5 Long-Name"];

    pretty_table(
        fmat,head,
        alignment=[:c,:c,:l,:c,:l],
        crop = :none, tf = tf_compact
    );

end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/functions/nonlinear/affine.jl
# Defines buffered structures to store quadratic functions:
# AffineFunctionIneq, AffineFunctionEq, as well as the
# lower_interval_bound, interval_bound, and eliminate_fixed_variables!
# functions associated with each structure.
#############################################################################

###
### Structure definitions
###

"""
$(TYPEDEF)

Current only used for bound tightening. Stores a representation
of an affine inequality.
"""
mutable struct AffineFunctionIneq <: AbstractEAGOConstraint
    terms::Vector{Tuple{Float64,Int}}
    constant::Float64
    len::Int
end

"""
$(TYPEDEF)

Current only used for bound tightening. Stores a representation
of an affine equality.
"""
mutable struct AffineFunctionEq <: AbstractEAGOConstraint
    terms::Vector{Tuple{Float64,Int}}
    constant::Float64
    len::Int
end

###
### Constructor definitions
###
AffineFunctionIneq() = AffineFunctionIneq(Tuple{Float64,Int}[], 0.0, 0)
function AffineFunctionIneq(func::SAF, set::LT)
    terms = map(x -> (x.coefficient, x.variable_index.value), func.terms)
    return AffineFunctionIneq(terms, func.constant - set.upper, length(func.terms))
end

function AffineFunctionIneq(func::SAF, set::GT)
    terms = map(x -> (-x.coefficient, x.variable_index.value), func.terms)
    return AffineFunctionIneq(terms, set.lower - func.constant, length(func.terms))
end

AffineFunctionEq() = AffineFunctionEq(Tuple{Float64,Int}[], 0.0, 0)
function AffineFunctionEq(func::SAF, set::ET)
    terms = map(x -> (x.coefficient, x.variable_index.value), func.terms)
    return AffineFunctionEq(terms, func.constant - set.value, length(func.terms))
end

###
### Parsing definitions
###

function eliminate_fixed_variables!(f::T, v::Vector{VariableInfo}) where T <: Union{AffineFunctionIneq,
                                                                                AffineFunctionEq}
    deleted_count = 0
    index = 1
    while i + deleted_count <= f.len
        coeff, indx = @inbounds f.terms[i]
        variable_info = @inbounds v[indx]
        if variable_info.is_fixed
            f.constant += coeff*variable_info.lower_bound
            deleteat!(f.terms, i)
            deleted_count += 1
        else
            i += 1
        end
    end
    f.len -= deleted_count
    return nothing
end


text:
include("simulator.jl")

# ## Bayesian wrapper
using Gen, Genify, Distributions

trunc_normal(μ, σ, lb, ub) =
    Genify.WrappedDistribution(truncated(Normal(μ, σ), lb, ub))

# Set up fixed model parameters
function create_params(; β=0.5)
    Ns = [50; 50; 50]
    migration_rates = [0.95 0.04 0.01; 0.10 0.85 0.05; 0.01 0.04 0.95]
    β_det = β * ones(3)
    β_und = β_det ./ 10
    death_rate = 0.0
    Is = [10; 0; 0]
    return @dict(Ns, migration_rates, β_det, β_und, death_rate, Is)
end

# Define observation model
@gen function observe(model::ABM, noise::Float64=5.0)
    obs = zeros(length(nodes(model)) * 2)
    for city in nodes(model)
        infected = count(a.status == :I for a in get_node_agents(city, model))
        recovered = count(a.status == :R for a in get_node_agents(city, model))
        obs[[city*2-1, city*2]] = [infected, recovered]
        {:infected => city} ~ trunc_normal(infected, noise, 0, Inf)()
        {:recovered => city} ~ trunc_normal(recovered, noise, 0, Inf)()
    end
    for (i, agent_id) in enumerate(model.scheduler(model))
        agent = model.agents[agent_id]
        probs = [c == agent.pos ? 0.99 : 0.01/(model.C-1) for c in 1:model.C]
        {:location => i} ~ Gen.categorical(probs)
    end
    return obs
end

# Genify the Agents.step! method
abm_step! = genify(Agents.step!, AgentBasedModel, Any, Int, Bool; recurse=true)

# Define Unfold step
@gen function bayesian_sir_step!(t::Int, _::Nothing,
                                model::AgentBasedModel, noise::Float64)
    {:agents} ~ abm_step!(model, agent_step!, 1, true)
    {:obs} ~ observe(model, noise)
    return nothing
end

# Wrap SIR step model in Gen model with parameter uncertainty
@gen (static) function bayesian_sir(T::Int, noise::Float64)
    β ~ uniform_continuous(0, 2.0)
    params = create_params(β=β)
    model = model_initiation(; params...)
    {:step} ~ Unfold(bayesian_sir_step!)(T, nothing, model, noise)
    return model
end

load_generated_functions()

# Test run the model
params = create_params()
model = model_initiation(; params...)
step!(model, agent_step!, 1, true)
abm_step!(model, agent_step!, 1, true)


text:
using Dates, Printf
using Plots
using LongwaveModePropagator
using LMPTools

# Ground

tx = TRANSMITTER[:NAA]
rx = Receiver("Boulder", 40.01, -105.244, 0.0, VerticalDipole())

sigmas = [get_sigma(y, x) for y in LAT, x in LON]
heatmap(LON, LAT, sigmas, clims=(0, 0.11))

epsilons = [get_epsilon(y, x) for y in LAT, x in LON]
heatmap(LON, LAT, epsilons, clims=(0, 20))

pts = [forward(tx.longitude, tx.latitude, az, d) for d in distances]
heatmap(LON, LAT, sigmas, clims=(0, 0.11), xlims=(-110, -65), ylims=(35, 48))
scatter!(getfield.(pts, :lon), getfield.(pts, :lat))

# Magnetic field

# Compare to LWPC
tx = Transmitter{VerticalDipole}("fdtdnaa", 44.646, -67.281, VerticalDipole(), Frequency(24e3), 100e3)
rx = Receiver("Boulder", 40.01, -105.244, 0.0, VerticalDipole())
grounds, distances = groundsegments(tx.latitude, tx.longitude, rx.latitude, rx.longitude;
    require_permittivity_change=false)
bfields = igrf(tx, rx, 2020, distances)

species = Species(LMP.QE, LMP.ME, z->waitprofile(z, 82, 0.55), electroncollisionfrequency)
ground = GROUND[10]
wvg = SegmentedWaveguide([HomogeneousWaveguide(bfields[i], species, ground, distances[i]) for i in 1:length(distances)])
sampler = GroundSampler(0:5e3:2000e3, Fields.Ez)
e, a, p = propagate(wvg, tx, sampler)

ld, la, lp = readlog(LMPTools.project_path("data", "bfield.log"))

plot(sampler.distance/1000, a, label="LMP",
     ylabel="Amplitude (dB)", xlabel="Range (km)", linewidth=1.5)
plot!(ld, la, label="LWPC", linewidth=1.5)

# Zenith angle

dt = DateTime(2021, 2, 1)
lats = 15:89
lons = -160:-60
szas = [zenithangle(la, lo, dt) for la in lats, lo in lons]
heatmap(lons, lats, szas,
        color=cgrad(:starrynight, [50, 70, 80, 85, 90, 95, 100, 110, 130]/180, rev=true), clims=(0, 180))

function gpmap(file, lats::AbstractRange, lons::AbstractRange, data)
    size(data) == (length(lats), length(lons)) || throw(ArgumentError("data not compatible with lats, lons"))
    issorted(lons) && issorted(lats) || throw(ArgumentError("lats and lons must be sorted"))

    δlat = step(lats)/2
    δlon = step(lons)/2

    open(file, "w") do f
        for j in eachindex(lons)
            for i in eachindex(lats)

                str = @sprintf("%f, %f, %f\n", lons[j]-δlon, lats[i]-δlat, szas[i,j])
                write(f, str)
                str = @sprintf("%f, %f, %f\n", lons[j]-δlon, lats[i]+δlat, szas[i,j])
                write(f, str)
                str = @sprintf("%f, %f, %f\n", lons[j]+δlon, lats[i]+δlat, szas[i,j])
                write(f, str)
                str = @sprintf("%f, %f, %f\n", lons[j]+δlon, lats[i]-δlat, szas[i,j])
                write(f, str)

                write(f, "\n")
            end
        end
    end
end

dt = DateTime(2021, 2, 1)
lats = 35:60
lons = -110:-75
szas = [zenithangle(la, lo, dt) for la in lats, lo in lons]
gpmap("szas.csv", lats, lons, szas)

function gppm3d(file, lats, lons, data)
    lon=vec([lo for la in lats, lo in lons])
    lat=vec([la for la in lats, lo in lons])
    sza=vec(szas)

    open(file, "w") do f
        for j in eachindex(lons)
            for i in eachindex(lats)
                str = @sprintf("%f, %f, %f\n", lons[j], lats[i], szas[i,j])
                write(f, str)
            end
            write(f, "\n")
        end
    end
end

dt = DateTime(2021, 2, 1)
lats = 40:65
lons = -145:-60
szas = [zenithangle(la, lo, dt) for la in lats, lo in lons]
gppm3d("szas_3col.csv", lats, lons, szas)

# Ionospheres

dt = DateTime(2021, 2, 1)
lats = 15:89
lons = -160:-60
x = [ferguson(la, zenithangle(la, lo, dt), dt) for la in lats, lo in lons]
heatmap(lons, lats, getindex.(x,1), color=:amp, clims=(68, 90))

dt = DateTime(2020, 3, 1, 2)
szas = [zenithangle(la, lo, dt) for la in lats, lo in lons]

heatmap(lons, lats, szas,
        color=cgrad(:starrynight, [50, 70, 80, 85, 90, 95, 100, 110, 130]/180, rev=true), clims=(0, 180))

hprimes, betas = flatlinearterminator(szas)
heatmap(lons, lats, hprimes, color=:amp, clims=(68, 90))
heatmap(lons, lats, betas, color=:tempo, clims=(0.2, 0.6))

ionos = smoothterminator.(szas)
heatmap(lons, lats, getindex.(ionos,1), color=:amp, clims=(68, 90))
heatmap(lons, lats, getindex.(ionos,2), color=:tempo, clims=(0.2, 0.6))


text:
using HTTP2.Frame

@test decode(IOBuffer(encode(ContinuationFrame(false, 0x51, b"test")))) ==
    ContinuationFrame(false, 0x51, b"test")
@test decode(IOBuffer(encode(DataFrame(0x51, false, b"test")))) ==
    DataFrame(0x51, false, b"test")
@test decode(IOBuffer(encode(GoawayFrame(0x4f, 0x4, b"test")))) ==
    GoawayFrame(0x4f, 0x4, b"test")
@test decode(IOBuffer(encode(HeadersFrame(false, false, true, 0x51,
                                          false, 0x50,
                                          0x1, b"test")))) ==
                                              HeadersFrame(false, false, true, 0x51,
                                                           false, 0x50,
                                                           0x1, b"test")
@test decode(IOBuffer(encode(PingFrame(false, b"testtest")))) ==
    PingFrame(false, b"testtest")
@test decode(IOBuffer(encode(PriorityFrame(0x51, false, 0x50, 0x2)))) ==
    PriorityFrame(0x51, false, 0x50, 0x2)
@test decode(IOBuffer(encode(PushPromiseFrame(false, 0x51, 0x54, b"test")))) ==
    PushPromiseFrame(false, 0x51, 0x54, b"test")
@test decode(IOBuffer(encode(RstStreamFrame(0x51, 0x4)))) ==
    RstStreamFrame(0x51, 0x4)
@test decode(IOBuffer(encode(SettingsFrame(true, nothing)))) == SettingsFrame(true, nothing)
@test decode(IOBuffer(encode(WindowUpdateFrame(0x51, 0x2)))) ==
    WindowUpdateFrame(0x51, 0x2)


text:
# by lazarusA # HIDE
#using GLMakie, Random # HIDE
using CairoMakie, Random
CairoMakie.activate!() #HIDE
let
    Random.seed!(123)
    n = 15
    x, y, color = rand(n), rand(n), rand(n)
    cmaps = [:cool, :viridis, :inferno, :thermal]

    function FigGridHeatSharedCbarH()
        fig = Figure(resolution = (600, 400))
        axes = []
        c = 1
        for i in 1:1, j in 1:2
            ax = Axis(fig[i, j],aspect=1,xgridvisible = false,ygridvisible = false)
            pnts = heatmap!(rand(10,10), colormap= :thermal, colorrange=(0, 1))
            ax.xticks = [1,10]
            ax.yticks = [1,10]
            cbar = Colorbar(fig, pnts, vertical = false, flipaxis = false,
            width = Relative(3/4), height = 15, tickwidth = 2, ticklabelsize = 14)
            cbar.ticks = [0,0.5,1]
            fig[2, 1:2] = cbar
            c+=1
            push!(axes, ax)
        end
        hideydecorations!(axes[2], ticks = false)
        fig
    end
    fig = FigGridHeatSharedCbarH()
    fig
    save(joinpath(@__DIR__, "output", "heatmap1SharedCbar.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using LinearAlgebra
using MAT
using AdFem
using PyPlot
using Random
using SparseArrays
Random.seed!(118)

# geometry setup in domain [0,1]^2
solid_left = 0.4
solid_right = 0.6
solid_top = 0.5
solid_bottom = 0.505

# chip_left = 0.48
# chip_right = 0.52
# chip_top = 0.5
# chip_bottom = 0.505

k_mold = 0.014531
k_chipdie = 2.60475
k_air = 0.64357
nu = 0.47893 # equal to 1/Re
power_source = 0.8246295 #82.46295 = 1.0e6 divide by air rho cp   #0.0619 = 1.0e6 divide by chip die rho cp
buoyance_coef = 299102.83

u_std = 0.001
p_std = 0.000001225
T_infty = 300

function k_exact(x, y)
    1 + 1 / (1 + x^2)
end

m = 50
n = 50
h = 1/n
NT = 7    # number of iterations for Newton's method, 8 is good for m=400

# compute solid indices and chip indices
solid_fem_idx = Array{Int64, 1}([])
solid_fvm_idx = Array{Int64, 1}([])
# chip_fem_idx = Array{Int64, 1}([])
# chip_fvm_idx = Array{Int64, 1}([])
# chip_fem_top_idx = Array{Int64, 1}([])

for i = 1:(m+1)
    for j = 1:(n+1)
        if (i-1)*h >= solid_left-1e-9 && (i-1)*h <= solid_right+1e-9 && (j-1)*h >= solid_top-1e-9 && (j-1)*h <= solid_bottom+1e-9
            # print(i, j)
            global solid_fem_idx = [solid_fem_idx; (j-1)*(m+1)+i]
            # if (i-1)*h >= chip_left-1e-9 && (i-1)*h <= chip_right+1e-9 && (j-1)*h >= chip_top-1e-9 && (j-1)*h <= chip_bottom+1e-9
            #     global chip_fem_idx = [chip_fem_idx; (j-1)*(m+1)+i]
            # end
            # if (i-1)*h >= chip_left-1e-9 && (i-1)*h <= chip_right+1e-9 && (j-1)*h >= chip_top-1e-9 && (j-1)*h <= chip_top+1e-9
            #     global chip_fem_top_idx = [chip_fem_top_idx; (j-1)*(m+1)+i]
            # end
        end
    end
end

for i = 1:m
    for j = 1:n
        if (i-1)*h + h/2 >= solid_left-1e-9 && (i-1)*h + h/2 <= solid_right+1e-9 &&
            (j-1)*h + h/2 >= solid_top-1e-9 && (j-1)*h + h/2 <= solid_bottom+1e-9
            global solid_fvm_idx = [solid_fvm_idx; (j-1)*m+i]
            # if (i-1)*h + h/2 >= chip_left-1e-9 && (i-1)*h + h/2 <= chip_right+1e-9 && (j-1)*h + h/2 >= chip_top-1e-9 && (j-1)*h + h/2<= chip_bottom+1e-9
            #     global chip_fvm_idx = [chip_fvm_idx; (j-1)*m+i]
            # end
        end
    end
end

# initialize space varying k and heat source
xy = fem_nodes(m, n, h)
solid_x, solid_y = xy[solid_fem_idx, 1], xy[solid_fem_idx, 2]
k_solid = @. k_exact(solid_x, solid_y)
# k_solid=stack(k_solid)

k_fem = k_air * constant(ones((m+1)*(n+1)))
k_fem = scatter_update(k_fem, solid_fem_idx, k_solid)
# k_fem = scatter_update(k_fem, chip_fem_idx, k_chip)
kgauss = fem_to_gauss_points(k_fem, m, n, h)

heat_source_fem = zeros((m+1)*(n+1))
heat_source_fem[solid_fem_idx] .= power_source
# heat_source_fem[chip_fem_top_idx] .= 82.46295

heat_source_gauss = fem_to_gauss_points(heat_source_fem, m, n, h)

# chip_gauss_idx = [ 4 .* chip_fvm_idx; 4 .* chip_fvm_idx .- 1; 4 .* chip_fvm_idx .- 2; 4 .* chip_fvm_idx .- 3]
# heat_source_gauss = zeros(4*m*n)
# heat_source_gauss[chip_gauss_idx] .= power_source

# F1 = compute_fem_source_term1(eval_f_on_gauss_pts(ffunc_, m, n, h), m, n, h)
# F2 = compute_fem_source_term1(eval_f_on_gauss_pts(gfunc_, m, n, h), m, n, h)
# H = h^2*eval_f_on_fvm_pts(hfunc_, m, n, h)
B = constant(compute_interaction_matrix(m, n, h))

# compute F
Laplace = nu * constant(compute_fem_laplace_matrix1(m, n, h))
heat_source = constant(compute_fem_source_term1(heat_source_gauss, m, n, h))

LaplaceK = constant(compute_fem_laplace_matrix1(kgauss, m, n, h))
# xy = fem_nodes(m, n, h)
# x, y = xy[:, 1], xy[:, 2]
# k = @. k_nn(x, y); k=stack(k)
# kgauss = fem_to_gauss_points(k, m, n, h)
# LaplaceK = compute_fem_laplace_matrix1(kgauss, m, n, h)

bd = bcnode("all", m, n, h)

# only apply Dirichlet to velocity; set left bottom two points to zero to fix rank deficient problem for pressure

bd = [bd; bd .+ (m+1)*(n+1);
     2*(m+1)*(n+1)+1; 2*(m+1)*(n+1)+m;
    #  (2*(m+1)*(n+1)+m*n )+1:(2*(m+1)*(n+1)+m*n )+m+1]
     bd .+ (2*(m+1)*(n+1)+m*n )]

# add solid region into boundary condition for u, v, p
bd = [bd; solid_fem_idx; solid_fem_idx .+ (m+1)*(n+1); solid_fvm_idx .+ 2(m+1)*(n+1)]

function compute_residual(S)
    u, v, p, T = S[1:(m+1)*(n+1)],
        S[(m+1)*(n+1)+1:2(m+1)*(n+1)],
        S[2(m+1)*(n+1)+1:2(m+1)*(n+1)+m*n],
        S[2(m+1)*(n+1)+m*n+1:end]
    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    interaction = compute_interaction_term(p, m, n, h) # julia kernel needed
    f1 = compute_fem_source_term1(ugauss.*ux, m, n, h)
    f2 = compute_fem_source_term1(vgauss.*uy, m, n, h)
    f3 = -interaction[1:(m+1)*(n+1)]
    f4 = Laplace*u
    # f5 = -F1
    F = f1 + f2 + f3 + f4 #+ f5

    g1 = compute_fem_source_term1(ugauss.*vx, m, n, h)
    g2 = compute_fem_source_term1(vgauss.*vy, m, n, h)
    g3 = -interaction[(m+1)*(n+1)+1:end]
    g4 = Laplace*v
    # g5 = -F2
    T_gauss = fem_to_gauss_points(T, m, n, h)
    buoyance_term = - buoyance_coef * compute_fem_source_term1(T_gauss, m, n, h)

    G = g1 + g2 + g3 + g4 + buoyance_term #+ g5

    H0 = -B * [u;v] # + H

    T0 = LaplaceK * T + compute_fem_advection_matrix1(ugauss,vgauss, m, n, h) * T - heat_source
    R = [F;G;H0;T0]
    return R
end

function compute_jacobian(S)
    u, v, p, T = S[1:(m+1)*(n+1)],
        S[(m+1)*(n+1)+1:2(m+1)*(n+1)],
        S[2(m+1)*(n+1)+1:2(m+1)*(n+1)+m*n],
        S[2(m+1)*(n+1)+m*n+1:end]

    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    M1 = constant(compute_fem_mass_matrix1(ux, m, n, h))
    M2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h)) # a julia kernel needed
    M3 = Laplace
    Fu = M1 + M2 + M3

    Fv = constant(compute_fem_mass_matrix1(uy, m, n, h))

    N1 = constant(compute_fem_mass_matrix1(vy, m, n, h))
    N2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h))
    N3 = Laplace
    Gv = N1 + N2 + N3

    Gu = constant(compute_fem_mass_matrix1(vx, m, n, h))

    M = LaplaceK + constant(compute_fem_advection_matrix1(ugauss,vgauss, m, n, h))

    gradT = eval_grad_on_gauss_pts1(T, m, n, h)
    Tx, Ty = gradT[:,1], gradT[:,2]
    DU_TX = constant(compute_fem_mass_matrix1(Tx, m, n, h))       # (m+1)*(n+1), (m+1)*(n+1)
    DV_TY = constant(compute_fem_mass_matrix1(Ty, m, n, h))       # (m+1)*(n+1), (m+1)*(n+1)

    T_mat = -buoyance_coef * constant(compute_fem_mass_matrix1(m, n, h))
    T_mat = [SparseTensor(spzeros((m+1)*(n+1), (m+1)*(n+1))); T_mat]

    J0 = [Fu Fv
          Gu Gv]

    J1 = [J0 -B' T_mat
        -B spdiag(zeros(size(B,1))) SparseTensor(spzeros(m*n, (m+1)*(n+1)))]

    N = 2*(m+1)*(n+1) + m*n
    J = [J1
        [DU_TX DV_TY SparseTensor(spzeros((m+1)*(n+1), m*n)) M]]
end

function solve_steady_cavityflow_one_step(S)
    residual = compute_residual(S)
    J = compute_jacobian(S)

    J, _ = fem_impose_Dirichlet_boundary_condition1(J, bd, m, n, h)
    residual = scatter_update(residual, bd, zeros(length(bd)))    # residual[bd] .= 0.0 in Tensorflow syntax

    d = J\residual
    residual_norm = norm(residual)
    op = tf.print("residual norm", residual_norm)
    d = bind(d, op)
    S_new = S - d
    return S_new
end

function condition(i, S_arr)
    i <= NT + 1
end

function body(i, S_arr)
    S = read(S_arr, i-1)
    op = tf.print("i=",i)
    i = bind(i, op)
    S_new = solve_steady_cavityflow_one_step(S)
    S_arr = write(S_arr, i, S_new)
    return i+1, S_arr
end

S_arr = TensorArray(NT+1)
S_arr = write(S_arr, 1, zeros(m*n+3*(m+1)*(n+1)))

i = constant(2, dtype=Int32)

_, S = while_loop(condition, body, [i, S_arr])
S = set_shape(stack(S), (NT+1, 2*(m+1)*(n+1)+m*n+(m+1)*(n+1)))

sess = Session(); init(sess)
output = run(sess, S)

matwrite("xzchipv2_figures1/data_m$(m)_n$n.mat",
    Dict(
        "V"=>output[end, :]
    ))

u0, v0, p0, t0 = zeros((m+1)*(n+1)), zeros((m+1)*(n+1)), zeros(m*n), zeros((m+1)*(n+1))

figure(figsize=(10,10))
subplot(221)
title("u velocity")
visualize_scalar_on_fem_points(output[NT+1, 1:(m+1)*(n+1)] .* u_std, m, n, h);gca().invert_yaxis()
subplot(222)
title("v velocity")
visualize_scalar_on_fem_points(output[NT+1, (m+1)*(n+1)+1:2*(m+1)*(n+1)] .* u_std, m, n, h);gca().invert_yaxis()
subplot(223)
visualize_scalar_on_fvm_points(output[NT+1, 2*(m+1)*(n+1)+1:2*(m+1)*(n+1)+m*n] .* p_std, m, n, h);gca().invert_yaxis()
title("pressure")
subplot(224)
title("temperature")
visualize_scalar_on_fem_points(output[NT+1, 2*(m+1)*(n+1)+m*n+1:end].* T_infty .+ T_infty, m, n, h);gca().invert_yaxis()
tight_layout()
savefig("xzchipv2_figures1/forward_soln.png")

print("Solution range:",
    "\n [u velocity] \t min:", minimum(output[NT+1, 1:(m+1)*(n+1)] .* u_std), ",\t max:", maximum(output[NT+1, 1:(m+1)*(n+1)] .* u_std),
    "\n [v velocity] \t min:", minimum(output[NT+1, (m+1)*(n+1)+1:2*(m+1)*(n+1)] .* u_std), ",\t max:", maximum(output[NT+1, (m+1)*(n+1)+1:2*(m+1)*(n+1)] .* u_std),
    "\n [pressure]   \t min:", minimum(output[NT+1, 2*(m+1)*(n+1)+1:2*(m+1)*(n+1)+m*n] .* p_std), ",\t max:", maximum(output[NT+1, 2*(m+1)*(n+1)+1:2*(m+1)*(n+1)+m*n] .* p_std),
    "\n [temperature]\t min:", minimum(output[NT+1, 2*(m+1)*(n+1)+m*n+1:end].* T_infty .+ T_infty), ",\t\t\t max:", maximum(output[NT+1, 2*(m+1)*(n+1)+m*n+1:end].* T_infty .+ T_infty))

# separate plots
# figure();visualize_scalar_on_fem_points(output[NT+1, 1:(m+1)*(n+1)] .* u_std, m, n, h);gca().invert_yaxis()
# figure();visualize_scalar_on_fem_points(output[NT+1, (m+1)*(n+1)+1:2*(m+1)*(n+1)] .* u_std, m, n, h);gca().invert_yaxis()
# figure();visualize_scalar_on_fvm_points(output[NT+1, 2*(m+1)*(n+1)+1:2*(m+1)*(n+1)+m*n] .* p_std, m, n, h);gca().invert_yaxis()
# figure();visualize_scalar_on_fem_points(output[NT+1, 2*(m+1)*(n+1)+m*n+1:end].* T_infty .+ T_infty, m, n, h);gca().invert_yaxis()


text:
using IntervalWavelets
using Test

@testset "Filters" begin
    filter_coefficients = rand(3)
    f = Filter(filter_coefficients, -1:1)

    @test support(f) == -1:1
    @test support_boundaries(f) == (-1, 1)
    @test length(f) == 3

    @test f[-1] == filter_coefficients[1]
    @test f[0] == filter_coefficients[2]
    @test f[1] == filter_coefficients[3]

    @test f[2] == 0.0
    @test f[-2] == 0.0
end

@testset "Interior filters" begin
    p = rand(2:8)

    @testset "Construct interior filters $phase" for phase in [:min, :symmlet]
        h = interior_filter(p, phase)
        @test vanishing_moments(h) == p
        @test length(h) == 2*p
    end

    @testset "Daubechies filters are normalized in l2 $phase" for phase in [:min, :symmlet]
        h = interior_filter(p, phase)
        @test sum(coefficients(h).^2) ≈ 1 atol = sqrt(eps())
    end

    @testset "Failure constructing interior filters" begin
        @test_throws DomainError interior_filter(0, :symmlet)
        @test_throws DomainError interior_filter(9, :symmlet)

        @test_throws MethodError interior_filter(p, :foo)
    end
end

@testset "Boundary filters" begin
    @testset "Construct boundary filter" for s in [IntervalWavelets.LEFT, IntervalWavelets.RIGHT]
        p = 2

        b = boundary_filters(p, s)
        @test side(b) == s
        @test vanishing_moments(b) == p
        @test length(filters(b)) == p
    end

    @testset "Failure constructing boundary filter" begin
        @test_throws DomainError boundary_filters(0, IntervalWavelets.LEFT)
        @test_throws DomainError boundary_filters(9, IntervalWavelets.LEFT)
    end
end


text:
module Trajectory

include("minimumjerk.jl")

end # module


text:
function accuracy(y::AbstractVector{T}, ŷ::AbstractVector{S}) where {T,S}
    @assert length(y) == length(ŷ)
    return sum(y .== ŷ) / length(y)
end

function accuracy(y::AbstractVector{T}, ŷ::AbstractMatrix{S}) where {T,S}
    max_idxs = dropdims(mapslices(x -> findmax(x)[2], ŷ; dims=1), dims=1)
    return accuracy(y, max_idxs)
end

# true labels on rows, predicted labels on columns
function confusion_matrix(y::AbstractVector{T}, ŷ::AbstractVector{T}; norm=false) where T
    labels = sort(collect(Set(y)))
    index = Dict(l => i for (i, l) in enumerate(labels))
    num_classes = length(labels)
    mat = zeros(Int, num_classes, num_classes)
    for (t, p) in zip(y, ŷ)
        mat[index[t], index[p]] += 1
    end
    if norm
        mat = mat ./ sum(mat, dims=2)
    end
    return mat
end

function recall(y::AbstractVector{T}, ŷ::AbstractVector{T}) where T
    tn, fn, fp, tp = confusion_matrix(y, ŷ)
    return tp / (tp + fn)
end

function precision(y::AbstractVector{T}, ŷ::AbstractVector{T}) where T
    tn, fn, fp, tp = confusion_matrix(y, ŷ)
    return tp / (tp + fp)
end


text:
using DataFrames

"""
    estimate_volume(
        model,
        p_fun,
        df::SubDataFrame,
        bounds,
        args...;
        n_sim = 10_000,
        parm_names,
        kwargs...
    )

Estimate volume of region with an eillipsoid and hit or miss bias adjustment.

# Arguments
- `model`: a model function that returns predictions given a vector of parameters
- `p_fun`: a function that that evaluates the qualitative data pattern
- `df::SubDataFrame`: a SubDataFrame containing MCMC samples in columns identities by `parm_names`
    and column containing data patterns named `pattern`
- `bounds`: a vector of tuples representing (lowerbound, upperbound) for each dimension in
- `args...`: arguments passed to `model` and `p_fun`

# Keywords

- `n_sim=10_000`: the number of samples for hit or miss bias adjustment
- `parm_names`: a vector of symbols representing parameter columns in `df`
- `kwargs...`: additional keyword arguments passed to `model` or `p_fun`
"""
function estimate_volume(
    model,
    p_fun,
    df::SubDataFrame,
    bounds,
    args...;
    n_sim = 10_000,
    parm_names,
    kwargs...
    )

    points = Array(df[:,parm_names])
    target = df.pattern[1]

    return estimate_volume(
        model,
        p_fun,
        points,
        target,
        bounds,
        args...;
        parm_names,
        n_sim,
        kwargs...
    )
end

max_stars_repo_name: itsdfish/ParameterSpacePartitions

max_forks_repo_name: itsdfish/ParameterSpacePartitions


text:
#=
This file is auto-generated. Do not edit.
=#

#! format: off

"""
    mutable struct HydroEnergyReservoir <: HydroGen
        name::String
        available::Bool
        bus::Bus
        active_power::Float64
        reactive_power::Float64
        rating::Float64
        prime_mover::PrimeMovers
        active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        reactive_power_limits::Union{Nothing, Min_Max}
        ramp_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
        time_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
        base_power::Float64
        storage_capacity::Float64
        inflow::Float64
        initial_storage::Float64
        operation_cost::OperationalCost
        storage_target::Float64
        conversion_factor::Float64
        time_at_status::Float64
        services::Vector{Service}
        dynamic_injector::Union{Nothing, DynamicInjection}
        ext::Dict{String, Any}
        time_series_container::InfrastructureSystems.TimeSeriesContainer
        internal::InfrastructureSystemsInternal
    end

# Arguments
- `name::String`
- `available::Bool`
- `bus::Bus`
- `active_power::Float64`
- `reactive_power::Float64`, validation range: `reactive_power_limits`, action if invalid: `warn`
- `rating::Float64`: Thermal limited MVA Power Output of the unit. <= Capacity, validation range: `(0, nothing)`, action if invalid: `error`
- `prime_mover::PrimeMovers`: Prime mover technology according to EIA 923
- `active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`
- `reactive_power_limits::Union{Nothing, Min_Max}`, action if invalid: `warn`
- `ramp_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}`: ramp up and ramp down limits in MW (in component base per unit) per minute, validation range: `(0, nothing)`, action if invalid: `error`
- `time_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}`: Minimum up and Minimum down time limits in hours, validation range: `(0, nothing)`, action if invalid: `error`
- `base_power::Float64`: Base power of the unit in MVA, validation range: `(0, nothing)`, action if invalid: `warn`
- `storage_capacity::Float64`: Maximum storage capacity in the reservoir (units can be p.u-hr or m^3)., validation range: `(0, nothing)`, action if invalid: `error`
- `inflow::Float64`: Baseline inflow into the reservoir (units can be p.u. or m^3/hr), validation range: `(0, nothing)`, action if invalid: `error`
- `initial_storage::Float64`: Initial storage capacity in the reservoir (units can be p.u-hr or m^3)., validation range: `(0, nothing)`, action if invalid: `error`
- `operation_cost::OperationalCost`: Operation Cost of Generation [`OperationalCost`](@ref)
- `storage_target::Float64`: Storage target at the end of simulation as ratio of storage capacity.
- `conversion_factor::Float64`: Conversion factor from flow/volume to energy: m^3 -> p.u-hr.
- `time_at_status::Float64`
- `services::Vector{Service}`: Services that this device contributes to
- `dynamic_injector::Union{Nothing, DynamicInjection}`: corresponding dynamic injection device
- `ext::Dict{String, Any}`
- `time_series_container::InfrastructureSystems.TimeSeriesContainer`: internal time_series storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct HydroEnergyReservoir <: HydroGen
    name::String
    available::Bool
    bus::Bus
    active_power::Float64
    reactive_power::Float64
    "Thermal limited MVA Power Output of the unit. <= Capacity"
    rating::Float64
    "Prime mover technology according to EIA 923"
    prime_mover::PrimeMovers
    active_power_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    reactive_power_limits::Union{Nothing, Min_Max}
    "ramp up and ramp down limits in MW (in component base per unit) per minute"
    ramp_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    "Minimum up and Minimum down time limits in hours"
    time_limits::Union{Nothing, NamedTuple{(:up, :down), Tuple{Float64, Float64}}}
    "Base power of the unit in MVA"
    base_power::Float64
    "Maximum storage capacity in the reservoir (units can be p.u-hr or m^3)."
    storage_capacity::Float64
    "Baseline inflow into the reservoir (units can be p.u. or m^3/hr)"
    inflow::Float64
    "Initial storage capacity in the reservoir (units can be p.u-hr or m^3)."
    initial_storage::Float64
    "Operation Cost of Generation [`OperationalCost`](@ref)"
    operation_cost::OperationalCost
    "Storage target at the end of simulation as ratio of storage capacity."
    storage_target::Float64
    "Conversion factor from flow/volume to energy: m^3 -> p.u-hr."
    conversion_factor::Float64
    time_at_status::Float64
    "Services that this device contributes to"
    services::Vector{Service}
    "corresponding dynamic injection device"
    dynamic_injector::Union{Nothing, DynamicInjection}
    ext::Dict{String, Any}
    "internal time_series storage"
    time_series_container::InfrastructureSystems.TimeSeriesContainer
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function HydroEnergyReservoir(name, available, bus, active_power, reactive_power, rating, prime_mover, active_power_limits, reactive_power_limits, ramp_limits, time_limits, base_power, storage_capacity, inflow, initial_storage, operation_cost=TwoPartCost(0.0, 0.0), storage_target=1.0, conversion_factor=1.0, time_at_status=INFINITE_TIME, services=Device[], dynamic_injector=nothing, ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), )
    HydroEnergyReservoir(name, available, bus, active_power, reactive_power, rating, prime_mover, active_power_limits, reactive_power_limits, ramp_limits, time_limits, base_power, storage_capacity, inflow, initial_storage, operation_cost, storage_target, conversion_factor, time_at_status, services, dynamic_injector, ext, time_series_container, InfrastructureSystemsInternal(), )
end

function HydroEnergyReservoir(; name, available, bus, active_power, reactive_power, rating, prime_mover, active_power_limits, reactive_power_limits, ramp_limits, time_limits, base_power, storage_capacity, inflow, initial_storage, operation_cost=TwoPartCost(0.0, 0.0), storage_target=1.0, conversion_factor=1.0, time_at_status=INFINITE_TIME, services=Device[], dynamic_injector=nothing, ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), internal=InfrastructureSystemsInternal(), )
    HydroEnergyReservoir(name, available, bus, active_power, reactive_power, rating, prime_mover, active_power_limits, reactive_power_limits, ramp_limits, time_limits, base_power, storage_capacity, inflow, initial_storage, operation_cost, storage_target, conversion_factor, time_at_status, services, dynamic_injector, ext, time_series_container, internal, )
end

# Constructor for demo purposes; non-functional.
function HydroEnergyReservoir(::Nothing)
    HydroEnergyReservoir(;
        name="init",
        available=false,
        bus=Bus(nothing),
        active_power=0.0,
        reactive_power=0.0,
        rating=0.0,
        prime_mover=PrimeMovers.HY,
        active_power_limits=(min=0.0, max=0.0),
        reactive_power_limits=nothing,
        ramp_limits=nothing,
        time_limits=nothing,
        base_power=0.0,
        storage_capacity=0.0,
        inflow=0.0,
        initial_storage=0.0,
        operation_cost=TwoPartCost(nothing),
        storage_target=0.0,
        conversion_factor=0.0,
        time_at_status=INFINITE_TIME,
        services=Device[],
        dynamic_injector=nothing,
        ext=Dict{String, Any}(),
        time_series_container=InfrastructureSystems.TimeSeriesContainer(),
    )
end

"""Get [`HydroEnergyReservoir`](@ref) `name`."""
get_name(value::HydroEnergyReservoir) = value.name
"""Get [`HydroEnergyReservoir`](@ref) `available`."""
get_available(value::HydroEnergyReservoir) = value.available
"""Get [`HydroEnergyReservoir`](@ref) `bus`."""
get_bus(value::HydroEnergyReservoir) = value.bus
"""Get [`HydroEnergyReservoir`](@ref) `active_power`."""
get_active_power(value::HydroEnergyReservoir) = get_value(value, value.active_power)
"""Get [`HydroEnergyReservoir`](@ref) `reactive_power`."""
get_reactive_power(value::HydroEnergyReservoir) = get_value(value, value.reactive_power)
"""Get [`HydroEnergyReservoir`](@ref) `rating`."""
get_rating(value::HydroEnergyReservoir) = get_value(value, value.rating)
"""Get [`HydroEnergyReservoir`](@ref) `prime_mover`."""
get_prime_mover(value::HydroEnergyReservoir) = value.prime_mover
"""Get [`HydroEnergyReservoir`](@ref) `active_power_limits`."""
get_active_power_limits(value::HydroEnergyReservoir) = get_value(value, value.active_power_limits)
"""Get [`HydroEnergyReservoir`](@ref) `reactive_power_limits`."""
get_reactive_power_limits(value::HydroEnergyReservoir) = get_value(value, value.reactive_power_limits)
"""Get [`HydroEnergyReservoir`](@ref) `ramp_limits`."""
get_ramp_limits(value::HydroEnergyReservoir) = get_value(value, value.ramp_limits)
"""Get [`HydroEnergyReservoir`](@ref) `time_limits`."""
get_time_limits(value::HydroEnergyReservoir) = value.time_limits
"""Get [`HydroEnergyReservoir`](@ref) `base_power`."""
get_base_power(value::HydroEnergyReservoir) = value.base_power
"""Get [`HydroEnergyReservoir`](@ref) `storage_capacity`."""
get_storage_capacity(value::HydroEnergyReservoir) = get_value(value, value.storage_capacity)
"""Get [`HydroEnergyReservoir`](@ref) `inflow`."""
get_inflow(value::HydroEnergyReservoir) = get_value(value, value.inflow)
"""Get [`HydroEnergyReservoir`](@ref) `initial_storage`."""
get_initial_storage(value::HydroEnergyReservoir) = get_value(value, value.initial_storage)
"""Get [`HydroEnergyReservoir`](@ref) `operation_cost`."""
get_operation_cost(value::HydroEnergyReservoir) = value.operation_cost
"""Get [`HydroEnergyReservoir`](@ref) `storage_target`."""
get_storage_target(value::HydroEnergyReservoir) = value.storage_target
"""Get [`HydroEnergyReservoir`](@ref) `conversion_factor`."""
get_conversion_factor(value::HydroEnergyReservoir) = value.conversion_factor
"""Get [`HydroEnergyReservoir`](@ref) `time_at_status`."""
get_time_at_status(value::HydroEnergyReservoir) = value.time_at_status
"""Get [`HydroEnergyReservoir`](@ref) `services`."""
get_services(value::HydroEnergyReservoir) = value.services
"""Get [`HydroEnergyReservoir`](@ref) `dynamic_injector`."""
get_dynamic_injector(value::HydroEnergyReservoir) = value.dynamic_injector
"""Get [`HydroEnergyReservoir`](@ref) `ext`."""
get_ext(value::HydroEnergyReservoir) = value.ext
"""Get [`HydroEnergyReservoir`](@ref) `time_series_container`."""
get_time_series_container(value::HydroEnergyReservoir) = value.time_series_container
"""Get [`HydroEnergyReservoir`](@ref) `internal`."""
get_internal(value::HydroEnergyReservoir) = value.internal

"""Set [`HydroEnergyReservoir`](@ref) `available`."""
set_available!(value::HydroEnergyReservoir, val) = value.available = val
"""Set [`HydroEnergyReservoir`](@ref) `bus`."""
set_bus!(value::HydroEnergyReservoir, val) = value.bus = val
"""Set [`HydroEnergyReservoir`](@ref) `active_power`."""
set_active_power!(value::HydroEnergyReservoir, val) = value.active_power = set_value(value, val)
"""Set [`HydroEnergyReservoir`](@ref) `reactive_power`."""
set_reactive_power!(value::HydroEnergyReservoir, val) = value.reactive_power = set_value(value, val)
"""Set [`HydroEnergyReservoir`](@ref) `rating`."""
set_rating!(value::HydroEnergyReservoir, val) = value.rating = set_value(value, val)
"""Set [`HydroEnergyReservoir`](@ref) `prime_mover`."""
set_prime_mover!(value::HydroEnergyReservoir, val) = value.prime_mover = val
"""Set [`HydroEnergyReservoir`](@ref) `active_power_limits`."""
set_active_power_limits!(value::HydroEnergyReservoir, val) = value.active_power_limits = set_value(value, val)
"""Set [`HydroEnergyReservoir`](@ref) `reactive_power_limits`."""
set_reactive_power_limits!(value::HydroEnergyReservoir, val) = value.reactive_power_limits = set_value(value, val)
"""Set [`HydroEnergyReservoir`](@ref) `ramp_limits`."""
set_ramp_limits!(value::HydroEnergyReservoir, val) = value.ramp_limits = set_value(value, val)
"""Set [`HydroEnergyReservoir`](@ref) `time_limits`."""
set_time_limits!(value::HydroEnergyReservoir, val) = value.time_limits = val
"""Set [`HydroEnergyReservoir`](@ref) `base_power`."""
set_base_power!(value::HydroEnergyReservoir, val) = value.base_power = val
"""Set [`HydroEnergyReservoir`](@ref) `storage_capacity`."""
set_storage_capacity!(value::HydroEnergyReservoir, val) = value.storage_capacity = set_value(value, val)
"""Set [`HydroEnergyReservoir`](@ref) `inflow`."""
set_inflow!(value::HydroEnergyReservoir, val) = value.inflow = set_value(value, val)
"""Set [`HydroEnergyReservoir`](@ref) `initial_storage`."""
set_initial_storage!(value::HydroEnergyReservoir, val) = value.initial_storage = set_value(value, val)
"""Set [`HydroEnergyReservoir`](@ref) `operation_cost`."""
set_operation_cost!(value::HydroEnergyReservoir, val) = value.operation_cost = val
"""Set [`HydroEnergyReservoir`](@ref) `storage_target`."""
set_storage_target!(value::HydroEnergyReservoir, val) = value.storage_target = val
"""Set [`HydroEnergyReservoir`](@ref) `conversion_factor`."""
set_conversion_factor!(value::HydroEnergyReservoir, val) = value.conversion_factor = val
"""Set [`HydroEnergyReservoir`](@ref) `time_at_status`."""
set_time_at_status!(value::HydroEnergyReservoir, val) = value.time_at_status = val
"""Set [`HydroEnergyReservoir`](@ref) `services`."""
set_services!(value::HydroEnergyReservoir, val) = value.services = val
"""Set [`HydroEnergyReservoir`](@ref) `ext`."""
set_ext!(value::HydroEnergyReservoir, val) = value.ext = val
"""Set [`HydroEnergyReservoir`](@ref) `time_series_container`."""
set_time_series_container!(value::HydroEnergyReservoir, val) = value.time_series_container = val


text:
@test length(Dojo.module_dir()) > 0
@test Dojo.scn(1.0) == " 1.0e+0"
@test Dojo.scn(0.0) == " 0.0e+0"
@test Dojo.scn(Inf) == " Inf"
@test Dojo.scn(123.1,
    digits=0) == " 1e+2"


text:
# tests for `get_optimal` method
# testing only the case f_2p
# using LikelihoodProfiler

@testset "testing get_optimal() for :LN_PRAXIS and f_2p, identifiable" begin
    res0 = get_optimal(
        [4., 1.],
        f_2p;
        silent = true,
        loss_tol = 1e-3,
        scan_tol = 0.,
        local_alg = :LN_PRAXIS
    )
    @test isapprox(res0.loss, 5., atol=1e-3)
    @test res0.ret == :FTOL_REACHED

    res1 = get_optimal(
        [4., 1.],
        f_2p;
        silent = true,
        loss_tol = 0.,
        scan_tol = 1e-3,
        local_alg = :LN_PRAXIS
    )
    @test isapprox(res1.params[1], 3., atol=1e-3)
    @test isapprox(res1.params[2], 4., atol=1e-3)
    @test res1.ret == :XTOL_REACHED

    res2 = get_optimal(
        [4., 1.],
        f_2p;
        silent = true,
        loss_tol = 0.,
        scan_tol = 0.,
        max_iter = 10,
        local_alg = :LN_PRAXIS
    )
    @test res2.ret == :MAXEVAL_REACHED

    res3 = get_optimal(
        [4., 1.],
        f_2p;
        silent = true,
        loss_tol = 1e-3,
        scan_tol = 0.,
        local_alg = :LN_PRAXIS,
        scale = [:log, :log]
    )
    @test isapprox(res3.loss, 5., atol=1e-3)
    @test res3.ret == :FTOL_REACHED
end

@testset "testing get_optimal() for :LN_PRAXIS and f_2p, nonidentifiable" begin
    res0 = get_optimal(
        [4., 1., 1.],
        f_2p;
        silent = true,
        loss_tol = 1e-3,
        scan_tol = 0.,
        local_alg = :LN_PRAXIS
    )
    @test isapprox(res0.loss, 5., atol=1e-3)
    @test res0.ret == :FTOL_REACHED

    res1 = get_optimal(
        [4., 1., 1.],
        f_2p;
        silent = true,
        loss_tol = 0.,
        scan_tol = 1e-3,
        local_alg = :LN_PRAXIS
    )
    @test isapprox(res1.params[1], 3., atol=1e-3)
    @test isapprox(res1.params[2], 4., atol=1e-3)
    @test res1.ret == :XTOL_REACHED

    res2 = get_optimal(
        [4., 1., 1.],
        f_2p;
        silent = true,
        loss_tol = 0.,
        scan_tol = 0.,
        max_iter = 10,
        local_alg = :LN_PRAXIS
    )
    @test res2.ret == :MAXEVAL_REACHED

    res3 = get_optimal(
        [4., 1., 1.],
        f_2p;
        silent = true,
        loss_tol = 1e-3,
        scan_tol = 0.,
        local_alg = :LN_PRAXIS,
        scale = [:log, :log, :log]
    )
    @test isapprox(res3.loss, 5., atol=1e-3)
    @test res3.ret == :FTOL_REACHED
end

f_2p(x) = 5.0 + (x[1]-3.0)^2 + (x[2]-4.0)^2
grad_2p(x) = [2*(x[1]-3.), 2*(x[2]-4.), 0.]

@testset "testing get_optimal() for :LD_MMA and f_2p, nonidentifiable" begin
    res0 = get_optimal(
        [4., 1., 1.],
        f_2p;
        silent = true,
        loss_tol = 1e-3,
        scan_tol = 0.,
        local_alg = :LD_MMA,
        loss_grad = grad_2p
    )
    @test isapprox(res0.loss, 5., atol=1e-3)
    @test res0.ret == :FTOL_REACHED

    res1 = get_optimal(
        [4., 1., 1.],
        f_2p;
        silent = true,
        loss_tol = 0.,
        scan_tol = 1e-3,
        local_alg = :LD_MMA,
        loss_grad = grad_2p
    )
    @test isapprox(res1.params[1], 3., atol=1e-3)
    @test isapprox(res1.params[2], 4., atol=1e-3)
    @test res1.ret == :XTOL_REACHED

    res2 = get_optimal(
        [4., 1., 1.],
        f_2p;
        silent = true,
        loss_tol = 0.,
        scan_tol = 0.,
        max_iter = 10,
        local_alg = :LD_MMA,
        loss_grad = grad_2p
    )
    @test res2.ret == :MAXEVAL_REACHED

    res3 = get_optimal(
        [4., 1., 1.],
        f_2p;
        silent = true,
        loss_tol = 1e-3,
        scan_tol = 0.,
        local_alg = :LD_MMA,
        scale = [:log, :log, :log],
        loss_grad = grad_2p
    )
    @test isapprox(res3.loss, 5., atol=1e-3)
    @test res3.ret == :FTOL_REACHED
end

@testset "testing get_optimal() for :LD_MMA and f_2p, nonidentifiable, finite" begin
    res0 = get_optimal(
        [4., 1., 1.],
        f_2p;
        silent = true,
        loss_tol = 1e-3,
        scan_tol = 0.,
        local_alg = :LD_MMA,
        loss_grad = :FINITE # :AUTODIFF
    )
    @test isapprox(res0.loss, 5., atol=1e-3)
    @test res0.ret == :FTOL_REACHED

    res1 = get_optimal(
        [4., 1., 1.],
        f_2p;
        silent = true,
        loss_tol = 0.,
        scan_tol = 1e-3,
        local_alg = :LD_MMA,
        loss_grad = :FINITE
    )
    @test isapprox(res1.params[1], 3., atol=1e-3)
    @test isapprox(res1.params[2], 4., atol=1e-3)
    @test res1.ret == :XTOL_REACHED

    res2 = get_optimal(
        [4., 1., 1.],
        f_2p;
        silent = true,
        loss_tol = 0.,
        scan_tol = 0.,
        max_iter = 10,
        local_alg = :LD_MMA,
        loss_grad = :FINITE
    )
    @test res2.ret == :MAXEVAL_REACHED

    res3 = get_optimal(
        [4., 1., 1.],
        f_2p;
        silent = true,
        loss_tol = 1e-3,
        scan_tol = 0.,
        local_alg = :LD_MMA,
        scale = [:log, :log, :log],
        loss_grad = :FINITE
    )
    @test isapprox(res3.loss, 5., atol=1e-3)
    @test res3.ret == :FTOL_REACHED
end

@testset "testing get_optimal() for :LD_MMA and f_2p, nonidentifiable, autodiff" begin
    res0 = get_optimal(
        [4., 1., 1.],
        f_2p;
        silent = true,
        loss_tol = 1e-3,
        scan_tol = 0.,
        local_alg = :LD_MMA,
        loss_grad = :AUTODIFF
    )
    @test isapprox(res0.loss, 5., atol=1e-3)
    @test res0.ret == :FTOL_REACHED

    res1 = get_optimal(
        [4., 1., 1.],
        f_2p;
        silent = true,
        loss_tol = 0.,
        scan_tol = 1e-3,
        local_alg = :LD_MMA,
        loss_grad = :AUTODIFF
    )
    @test isapprox(res1.params[1], 3., atol=1e-3)
    @test isapprox(res1.params[2], 4., atol=1e-3)
    @test res1.ret == :XTOL_REACHED

    res2 = get_optimal(
        [4., 1., 1.],
        f_2p;
        silent = true,
        loss_tol = 0.,
        scan_tol = 0.,
        max_iter = 10,
        local_alg = :LD_MMA,
        loss_grad = :AUTODIFF
    )
    @test res2.ret == :MAXEVAL_REACHED

    res3 = get_optimal(
        [4., 1., 1.],
        f_2p;
        silent = true,
        loss_tol = 1e-3,
        scan_tol = 0.,
        local_alg = :LD_MMA,
        scale = [:log, :log, :log],
        loss_grad = :AUTODIFF
    )
    @test isapprox(res3.loss, 5., atol=1e-3)
    @test res3.ret == :FTOL_REACHED
end


text:
function log(a::Real, x::Real, ε::Real)
    z = x
    t = 1
    y = 0
    while z > a || z < 1/a || t > ε
        if z > a
            z /= a
            y += t
        elseif z < 1/a
            z *= a
            y -= t
        else
            t /= 2
            z *= z
        end
    end
    return y
end


text:
# See also https://github.com/JuliaLang/julia/pull/35778

struct Singleton{T}
    value::T
end

Base.ndims(::Singleton) = 0
Base.ndims(::Type{<:Singleton}) = 0
Base.size(::Singleton) = ()
Base.axes(::Singleton) = ()
Base.length(::Singleton) = 1
Base.eltype(::Type{Singleton{T}}) where {T} = T
Base.getindex(s::Singleton) = s.value
Base.getindex(s::Singleton, ::CartesianIndex{0}) = s.value
Base.iterate(s::Singleton) = (s.value, nothing)
Base.iterate(::Singleton, ::Nothing) = nothing

Base.IteratorSize(::Type{<:Singleton}) = Base.HasLength()
Base.IteratorEltype(::Type{Singleton{T}}) where {T} = Base.HasEltype()


text:
#= Load these if running these tests in terminal
import Pkg; Pkg.activate("../")  # CytofDensityEstimation
include("runtests.jl")
=#

ENV["GKSwstype"] = "nul"  # For StatsPlots

@testset "Test fit" begin
  K = 3
  # use_big_data = true
  use_big_data = false

  if use_big_data
    data = CDE.Model.Data([fill(-Inf, 30000); randn(70000)],
                          [fill(-Inf, 20000); randn(60000)])
  else
    # data = CDE.Model.Data([-Inf, 1.0, 2.0], [-1.0, -Inf, 3.0, 4.0])
    data = CDE.Model.Data([fill(-Inf, 300); randn(700) .+ 1],
                          [fill(-Inf, 300); randn(700) .+ 1])
    # NOTE: When the two distributions are the same, E[beta|data,param]
    # approaches 50% because the parameters gammaC, gammaT, etaC, etaT can be
    # qualitatively similar. So, p needs to be informatively small apriori?
  end

  prior = CDE.Model.Prior(K, p=Beta(1, 99))
  println(prior)
  state = CDE.Model.State(data, prior)
  tuners = CDE.Model.Tuners(K)

  chain, laststate, summarystats = CDE.fit(state, data, prior, tuners,
                                           nsamps=[500], nburn=0)

  CDE.Model.printsummary(chain, summarystats)

  println("Make plots ...")
  imgdir = "img/test"; mkpath(imgdir)
  @time CDE.Model.plotpostsummary(chain, summarystats, data.yC, data.yT,
                                  imgdir, bw_postpred=0.2,
                                  ygrid=CDE.Model.default_ygrid())
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
#!/usr/bin/env julia

using SphericalFourierBesselDecompositions
SFB = SphericalFourierBesselDecompositions

using Test
using LinearAlgebra

@testset "Modes" begin
    @testset "AnlmModes" begin
        rmin = 100.0
        rmax = 10000.0
        kmax = 0.01
        #rmin = 500.0
        #rmax = 3000.0
        #kmax = 0.2

        @time modes = SFB.AnlmModes(kmax, rmin, rmax, cache=false)
        @show modes.kmax modes.nmax modes.lmax modes.nside
        @test 2*modes.nside >= modes.lmax
        @test maximum(modes.lmax_n) == modes.lmax
        @test maximum(modes.nmax_l) == modes.nmax
        @show SFB.getnlmsize(modes)
        @time for n=1:modes.nmax, l=0:modes.lmax_n[n], m=0:l
            idx = SFB.getidx(modes, n, l, m)
            n2, l2, m2 = SFB.getnlm(modes, idx)
            @test n == n2
            @test l == l2
            @test m == m2
        end
    end

    @testset "ClnnModes" begin
        rmin = 500.0
        rmax = 2000.0
        kmax = 0.02
        @time anlmmodes = SFB.AnlmModes(kmax, rmin, rmax)
        #@time anlmmodes = SFB.AnlmModes(2, 0, rmin, rmax)
        @time clnnmodes = SFB.ClnnModes(anlmmodes, Δnmax=4)
        @show anlmmodes.nmax anlmmodes.lmax anlmmodes.nside
        @show clnnmodes.Δnmax clnnmodes.Δnmax_l
        @show SFB.getlnnsize(clnnmodes)
        @test maximum(clnnmodes.Δnmax_l) == clnnmodes.Δnmax
        idx2 = 1
        @time for l=0:clnnmodes.amodes.lmax, Δn=0:clnnmodes.Δnmax_l[l+1], n̄=1:clnnmodes.amodes.nmax_l[l+1]-Δn
            idx = SFB.getidx(clnnmodes, l, n̄+Δn, n̄)
            idx3 = SFB.getidx(clnnmodes, l, n̄, n̄+Δn)
            @test idx == idx3
            l2, n2, n2′ = SFB.getlnn(clnnmodes, idx)
            Δn2 = n2′ - n2
            n̄2 = n2
            #@show idx, l,Δn,n̄, l2,Δn2,n̄2
            @test l2 <= clnnmodes.amodes.lmax
            @test Δn2 <= clnnmodes.Δnmax_l[l+1]
            @test n̄2 <= clnnmodes.amodes.nmax_l[l+1] - Δn
            @test l == l2
            @test Δn == Δn2
            @test n̄ == n̄2
            @test idx == idx2
            idx2 += 1
        end
    end

    @testset "ClnnBinnedModes" begin
        rmin = 1370.0
        rmax = 3540.0
        kmax = 0.04
        amodes = SFB.AnlmModes(kmax, rmin, rmax)
        cmodes = SFB.ClnnModes(amodes, Δnmax=0)
        @time w̃, v = SFB.bandpower_binning_weights(cmodes, Δℓ=18)
        @time bcmodes = SFB.ClnnBinnedModes(w̃, v, cmodes)
        @show unique(sort(bcmodes.LKK[1,:]))
        incomplete_bins = SFB.get_incomplete_bins(w̃)
        @show incomplete_bins
        @show bcmodes.LKK[2,incomplete_bins]
        @show maximum(kmax .- bcmodes.LKK[2,incomplete_bins])
        @test all(bcmodes.LKK[2,incomplete_bins] .>= kmax - 0.002884)

        rmin = 500.0
        rmax = 2000.0
        nmax = 10
        lmax = 11
        amodes = SFB.AnlmModes(nmax, lmax, rmin, rmax)
        cmodes = SFB.ClnnModes(amodes, Δnmax=0)

        @time w̃, v = SFB.bandpower_binning_weights(cmodes, Δℓ=1)
        @time bcmodes = SFB.ClnnBinnedModes(w̃, v, cmodes)
        @show bcmodes.LKK[1,:]
        @show unique(bcmodes.LKK[1,:])
        @test all(unique(bcmodes.LKK[1,:]) .≈ 0:lmax)

        @time w̃, v = SFB.bandpower_binning_weights(cmodes, Δℓ=2)
        @time bcmodes = SFB.ClnnBinnedModes(w̃, v, cmodes)
        @test all(unique(bcmodes.LKK[1,:]) .≈ 0.5:2:lmax)

        @time w̃, v = SFB.bandpower_binning_weights(cmodes, Δℓ=3)
        @time bcmodes = SFB.ClnnBinnedModes(w̃, v, cmodes)
        @show bcmodes.LKK[1,:]
        @show unique(bcmodes.LKK[1,:])
        @show unique(SFB.getlkk(bcmodes)[1,:])
        @test all(unique(bcmodes.LKK[1,:]) .≈ 1:3:lmax)

        rmin = 500.0
        rmax = 2000.0
        nmax = 10
        lmax = 10
        amodes = SFB.AnlmModes(nmax, lmax, rmin, rmax)
        cmodes = SFB.ClnnModes(amodes, Δnmax=0)

        w̃, v = SFB.bandpower_binning_weights(cmodes, Δℓ=1)
        bcmodes = SFB.ClnnBinnedModes(w̃, v, cmodes)
        @test all(unique(bcmodes.LKK[1,:]) .≈ 0:lmax)
        @test w̃ == I
        @test v == I

        @show bcmodes.LKK[1,:]
        @show bcmodes.LKK[2,:]
        @show bcmodes.LKK[3,:]
        @test bcmodes.LKK[2,:] == bcmodes.LKK[3,:]
        @show unique(bcmodes.LKK[1,:])
        @show unique(bcmodes.LKK[2,:])
        @show unique(bcmodes.LKK[3,:])

        i = 91
        @show "====",i
        l, k, k′ = SFB.getlkk(bcmodes, i)
        idx = SFB.getidxapprox(bcmodes, l, k, k′)
        L, K, K′ = SFB.getlkk(bcmodes, idx)
        @show (l,k,k′),i (l,K,K′),idx
        @test i == idx

        i = 12
        @show "====",i
        l, k, k′ = SFB.getlkk(bcmodes, i)
        idx = SFB.getidxapprox(bcmodes, l, k, k′)
        L, K, K′ = SFB.getlkk(bcmodes, idx)
        @show (l,k,k′),i (l,K,K′),idx
        @test i == idx

        for i=1:SFB.getlnnsize(bcmodes)
            l, k, k′ = SFB.getlkk(bcmodes, i)
            idx = SFB.getidxapprox(bcmodes, l, k, k′)
            @test i == idx
        end
    end
end

# vim: set sw=4 et sts=4 :


text:
"""
    GammaTrickHomotopy(start, target [, γ])
    GammaTrickHomotopy(start, target [, seed])

Construct the homotopy (1-t)target + tγ(t)start, with t from 1 to 0.
The gamma trick is a special type of homotopy interpolation.
One chooses an ``γ`` uniformly from the (complex) unit circle.
This yields a more generic interpolation.

Possible to create it optionally with a seed.
"""
struct GammaTrickHomotopy{T<:Number, Start<:AbstractPolySystem{T}, Target<:AbstractPolySystem{T}} <: AbstractHomotopy{T}
    start::Start
    target::Target
    γ::Complex128

    function GammaTrickHomotopy{T, Start, Target}(start::AbstractPolySystem{T}, target::AbstractPolySystem{T}, γ::Complex128) where {T<:Number, Start<:AbstractPolySystem{T}, Target<:AbstractPolySystem{T}}
        if (homogenized(start) != homogenized(target))
            return error("start and target have to be both either homogenized or not")
        end
        if (nvariables(start) != nvariables(target))
            return error("The number of unknowns of start and target system doesn't match." *
                "The start system of the homotopy has $(nvariables(start)) unknowns " *
                "and the target system has $(nvariables(target)) unknowns.")
        end

        if (length(target) != length(start))
            return error("The number of equations of start and target system doesn't match." *
                "The start system of the homotopy has $(length(start)) equations " *
                "and the target system has $(length(target)) equations.")
        end

        new(start, target, γ)
    end
end

function GammaTrickHomotopy(start::Start,target::Target, γ::Complex128) where {T<:Number, Start<:AbstractPolySystem{T}, Target<:AbstractPolySystem{T}}
    GammaTrickHomotopy{T, Start, Target}(start,target,γ)
end
# promote if they don't have the same coefficients
function GammaTrickHomotopy(start::AbstractPolySystem,target::AbstractPolySystem, γ::Complex128)
    T = promote_type(coefftype(start), coefftype(target))
    Start = promote_type(T, typeof(start))
    Target = promote_type(T, typeof(target))
    GammaTrickHomotopy{T, Start, Target}(convert(Start, start), convert(Target, target), γ)
end
GammaTrickHomotopy(s, t, γ::Complex128) = GammaTrickHomotopy(PolySystem(s), PolySystem(t), γ)

GammaTrickHomotopy(start, target) = GammaTrickHomotopy(start, target, exp(im * (rand() * 2π - π)))
function GammaTrickHomotopy(start, target, seed::Int)
    srand(seed)
    θ = rand() * 2π - π
    GammaTrickHomotopy(start, target, exp(im * θ))
end

function Base.promote_type(::Type{GammaTrickHomotopy{T, Start, Target}}, ::Type{S}) where {S,T, Start, Target}
    GammaTrickHomotopy{promote_type(T,S), promote_type(Start, S), promote_type(Target, S)}
end
function Base.promote_type(::Type{GammaTrickHomotopy{T, Start1, Target1}}, ::GammaTrickHomotopy{S, Start2, Target2}) where {S,T, Start1, Start2, Target1, Target2}
    GammaTrickHomotopy{promote_type(T,S), promote_type(Start1, Start2), promote_type(Target1, Target2)}
end

function Base.convert(::Type{GammaTrickHomotopy{T, Start, Target}}, H::GammaTrickHomotopy) where {T, Start, Target}
    GammaTrickHomotopy{T, Start, Target}(convert(Start, H.start), convert(Target, H.target), H.γ)
end

function Base.show(io::IO, H::GammaTrickHomotopy)
    print(io, typeof(H), ":\n")
    print(io, "* γ: ", H.γ, "\n")
    println(io, "* start:")
    println(io, H.start)
    println(io, "* target:")
    print(io, H.target)
end

function evaluate(H::GammaTrickHomotopy, x, t)
    (1 - t) * evaluate(H.target, x) + t * H.γ * evaluate(H.start, x)
end
(H::GammaTrickHomotopy)(x,t) = evaluate(H,x,t)

function substitute(H::GammaTrickHomotopy, pair::Pair{Symbol,<:Number})
    GammaTrickHomotopy(substitute(H.start, pair), substitute(H.target, pair), H.γ)
end

function removepoly(H::GammaTrickHomotopy, i::Int)
    GammaTrickHomotopy(removepoly(H.start,i), removepoly(H.target,i), H.γ)
end

"""
    differentiate(H::GammaTrickHomotopy)

Differentiates `H` w.r.t. the place and returns an evaluation function
`(x,t) -> t * H.γ * J_s(x) + (1 - t) * J_t(x)`
where `J_s` is the jacobian of the start system and
`J_t` is the jacobian of the target system.
"""
function differentiate(H::GammaTrickHomotopy)
    J_s = differentiate(H.start)
    J_t = differentiate(H.target)

    (x, t) -> t .* H.γ .* J_s(x) .+ (1 - t) .* J_t(x)
end

"""
    dt(H::GammaTrickHomotopy)

Differentiates `H` w.r.t. the the time and returns an evaluation function
`(x,t) -> T`
"""
dt(H::GammaTrickHomotopy) = (x, ::Number) -> H.γ * evaluate(H.start, x) .- evaluate(H.target, x)

"""
    homogenize(H::GammaTrickHomotopy)

Make the start and target system of `H` homogenous.
"""
function homogenize(H::GammaTrickHomotopy)
    GammaTrickHomotopy(homogenize(H.start), homogenize(H.target), H.γ)
end

"""
    homogenized(H::GammaTrickHomotopy)

Checks whether `H` was homogenized.
"""
homogenized(H::GammaTrickHomotopy) = homogenized(H.start)

ishomogenous(H::GammaTrickHomotopy) = ishomogenous(H.start) && ishomogenous(H.target)

"""
    nvariables(H::GammaTrickHomotopy)

The number variables of `H`.
"""
function nvariables(H::GammaTrickHomotopy)
    # Constructor guarantees that start and target have the same number of variables
    nvariables(H.start)
end
function length(H::GammaTrickHomotopy)
    # Constructor guarantees that start and target have the same length
    length(H.start)
end
"""
    degrees(H::GammaTrickHomotopy)

The (total) degrees of the polynomials of the system `P`.
"""
degrees(H::GammaTrickHomotopy) = max.(degrees(H.start), degrees(H.target))
"""
    startsystem(H::GammaTrickHomotopy)

The startsystem of `H`.
"""
startsystem(H::GammaTrickHomotopy) = H.start
"""
    targetsystem(H::GammaTrickHomotopy)

The targetsystem of `H`.
"""
targetsystem(H::GammaTrickHomotopy) = H.target

"""
    weylnorm(H, t)

Computes the weylnorm of the homotopy to the given time `t`.

## Explanation
For ``H = (1-t)F+tγG`` we have
```math
\begin{align*}
<H,H> &= <(1-t)F+tγG,(1-t)F+tγG> \\
      &= <(1-t)F,(1-t)F+tγG> + <tγG,(1-t)F+tγG> \\
      &= <(1-t)F,(1-t)F> + <(1-t)F,tγG> + <tγG,(1-t)F> + <tγG,tγG> \\
      &= <(1-t)F,(1-t)F> + 2*real(<(1-t)F,tγG>) + <tγG,tγG> \\
      &= (1-t)^2<F,F> + (1-t)t*2*real(γ)*real(<F,G>) + abs(γ)*t^2<G,G> \\
\end{align*}
```
"""
function weylnorm(H::GammaTrickHomotopy{T}, t::Float64) where {T<:Complex}
    F = H.target
    G = H.start
    sqrt((1 - t)^2 * weyldot(F, F) + (1 - t) * t * 2 * real(H.γ) * real(weyldot(F, G)) + t^2 * abs2(H.γ) * weyldot(G, G))
end