text:
# %% chaoticity map
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

using OrdinaryDiffEq
diffeq = (alg = Vern9(), abstol = 1e-9, reltol = 1e-9)

fig = figure(figsize = (figx/2, figy))
ax = gca()
hh = Systems.henonheiles()
ics = Systems.henonheiles_ics(0.13, 15)
cmap = matplotlib.cm.get_cmap("viridis")

scres = nothing
for ic in ics
    psos = poincaresos(hh, (1, 0.0), 2000; u0 = ic)
    λ = lyapunov(hh, 10000; u0 = ic, diffeq...)
    λmax = 0.06
    v = clamp(λ/λmax, 0, 1)
    global scres = ax.scatter(psos[:, 2], psos[:, 4];
        color = cmap(v), s = 5, edgecolors = "0.5", linewidths = 0.1,
        vmin = 0, vmax = 1,
    )
end

cb = colorbar(scres)
cb.set_ticks([0, 1])
cb.set_ticklabels(["0", "0.06"])
cb.set_label("chaoticity (\$\\lambda_1\$)"; labelpad=-50)
ax.set_xticks(-0.5:0.4:0.8)
ax.set_xlabel("\$y\$", labelpad = -20)
ax.set_yticks(-0.5:0.3:0.5)
ax.set_ylabel("\$p_y\$", labelpad = -10)
fig.tight_layout(; pad = 0.25)
fig.subplots_adjust(bottom = 0.16, top = 0.95, left = 0.17, wspace = 0.01)
wsave(plotsdir("3", "chaoticity"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
module UNet

# To do: better doc

using Flux

export chcat, uchain, unet

include("utils.jl")
include("uchain.jl")
include("unet.jl")

end

ext: jl


text:
__precompile__()

module RootTips

using NDFiles, Gtk, ProgressMeter

export main, batch_main

head = Gtk.GtkWindow("TrackRoots", 20, 20, false, true, visible=false, modal=true)

include(joinpath(Pkg.dir("RootTips"), "src", "utils.jl"))

include(joinpath(Pkg.dir("RootTips"), "src", "startPoints.jl"))

function main(ndfile::String = open_dialog("Pick an `.nd` file", head, ("*.nd",)))
    stages = nd2stages(ndfile)
    get_startpoints.(stages)
    run(`julia -e 'using TrackRoots; main(ARGS[1])' $ndfile`)
    # print(ndfile)
    # open(joinpath(tempdir(), "JULIA_NDFILE"), "w") do o
        # println(o, ndfile)
    # end
    # nothing
end

# batch

isnd(file::String) = last(splitext(file)) == ".nd"

function findall_nd(home::String)
    ndfiles = String[]
    for (root, dirs, files) in walkdir(home)
        filter!(isnd, files)
        append!(ndfiles, joinpath.(home, root, files))
    end
    return ndfiles
end

function batch_main(home::String = open_dialog("Pick a folder where all the `.nd` files are", action=GtkFileChooserAction.SELECT_FOLDER))
    ndfiles = findall_nd(home)
    n = length(ndfiles)
    # open(joinpath(tempdir(), "JULIA_NDFILES"), "w") do o
        info("Found $n `.nd` files")
        @showprogress 1 "Collecting start points…" for i in 1:n
            stages = nd2stages(ndfiles[i])
            get_startpoints.(stages)
            # println(o, ndfile)
        end
    # end
    run(`julia -e 'using TrackRoots; batch_main(ARGS[1])' $home`)
end

end # module


text:
import CuArrays: cu
import Yao.YaoArrayRegister: _measure, measure, measure!, measure_collapseto!, measure_remove!
import Yao.YaoBase: batch_normalize!
import Yao: expect

export cpu, cu, GPUReg

cu(reg::ArrayReg{B}) where B = ArrayReg{B}(CuArray(reg.state))
cpu(reg::ArrayReg{B}) where B = ArrayReg{B}(collect(reg.state))
const GPUReg{B, T, MT} = ArrayReg{B, T, MT} where MT<:CuArray

function batch_normalize!(s::CuSubArr, p::Real=2)
    p!=2 && throw(ArgumentError("p must be 2!"))
    s./=norm2(s, dims=1)
    s
end

@inline function tri2ij(l::Int)
    i = ceil(Int, sqrt(2*l+0.25)-0.5)
    j = l-i*(i-1)÷2
    i+1,j
end

############### MEASURE ##################
function measure(::ComputationalBasis, reg::GPUReg{1}, ::AllLocs; rng::AbstractRNG=Random.GLOBAL_RNG, nshots::Int=1)
    _measure(rng, reg |> probs |> Vector, nshots)
end

# TODO: optimize the batch dimension using parallel sampling
function measure(::ComputationalBasis, reg::GPUReg{B}, ::AllLocs; rng::AbstractRNG=Random.GLOBAL_RNG, nshots::Int=1) where B
    regm = reg |> rank3
    pl = dropdims(mapreduce(abs2, +, regm, dims=2), dims=2)
    _measure(rng, pl |> Matrix, nshots)
end

function measure!(::RemoveMeasured, ::ComputationalBasis, reg::GPUReg{B}, ::AllLocs; rng::AbstractRNG=Random.GLOBAL_RNG) where B
    regm = reg |> rank3
    nregm = similar(regm, 1<<nremain(reg), B)
    pl = dropdims(mapreduce(abs2, +, regm, dims=2), dims=2)
    pl_cpu = pl |> Matrix
    res_cpu = map(ib->_measure(rng, view(pl_cpu, :, ib), 1)[], 1:B)
    res = CuArray(res_cpu)
    CI = Base.CartesianIndices(nregm)
    @inline function kernel(nregm, regm, res, pl)
        state = (blockIdx().x-1) * blockDim().x + threadIdx().x
        if state <= length(nregm)
            @inbounds i,j = CI[state].I
            @inbounds r = Int(res[j])+1
            @inbounds nregm[i,j] = regm[r,i,j]/CUDAnative.sqrt(pl[r, j])
        end
        return
    end
    X, Y = cudiv(length(nregm))
    @cuda threads=X blocks=Y kernel(nregm, regm, res, pl)
    reg.state = reshape(nregm,1,:)
    B == 1 ? Array(res)[] : res
end

function measure!(::NoPostProcess, ::ComputationalBasis, reg::GPUReg{B, T}, ::AllLocs; rng::AbstractRNG=Random.GLOBAL_RNG) where {B, T}
    regm = reg |> rank3
    pl = dropdims(mapreduce(abs2, +, regm, dims=2), dims=2)
    pl_cpu = pl |> Matrix
    res_cpu = map(ib->_measure(rng, view(pl_cpu, :, ib), 1)[], 1:B)
    res = CuArray(res_cpu)
    CI = Base.CartesianIndices(regm)

    @inline function kernel(regm, res, pl)
        state = (blockIdx().x-1) * blockDim().x + threadIdx().x
        if state <= length(regm)
            @inbounds k,i,j = CI[state].I
            @inbounds rind = Int(res[j]) + 1
            @inbounds regm[k,i,j] = k==rind ? regm[k,i,j]/CUDAnative.sqrt(pl[k, j]) : T(0)
        end
        return
    end

    X, Y = cudiv(length(regm))
    @cuda threads=X blocks=Y kernel(regm, res, pl)
    B == 1 ? Array(res)[] : res
end

function measure!(rst::ResetTo, ::ComputationalBasis, reg::GPUReg{B, T}, ::AllLocs; rng::AbstractRNG=Random.GLOBAL_RNG) where {B, T}
    regm = reg |> rank3
    pl = dropdims(mapreduce(abs2, +, regm, dims=2), dims=2)
    pl_cpu = pl |> Matrix
    res_cpu = map(ib->_measure(rng, view(pl_cpu, :, ib), 1)[], 1:B)
    res = CuArray(res_cpu)
    CI = Base.CartesianIndices(regm)

    @inline function kernel(regm, res, pl, val)
        state = (blockIdx().x-1) * blockDim().x + threadIdx().x
        if state <= length(regm)
            @inbounds k,i,j = CI[state].I
            @inbounds rind = Int(res[j]) + 1
            @inbounds k==val+1 && (regm[k,i,j] = regm[rind,i,j]/CUDAnative.sqrt(pl[rind, j]))
            CuArrays.sync_threads()
            @inbounds k!=val+1 && (regm[k,i,j] = 0)
        end
        return
    end

    X, Y = cudiv(length(regm))
    @cuda threads=X blocks=Y kernel(regm, res, pl, rst.x)
    B == 1 ? Array(res)[] : res
end

import Yao.YaoArrayRegister: insert_qubits!, join
function YaoBase.batched_kron(A::Union{CuArray{T1, 3}, Adjoint{<:Any, <:CuArray{T1, 3}}}, B::Union{CuArray{T2, 3}, Adjoint{<:Any, <:CuArray{T2, 3}}}) where {T1 ,T2}
    res = CuArrays.zeros(promote_type(T1,T2), size(A,1)*size(B, 1), size(A,2)*size(B,2), size(A, 3))
    CI = Base.CartesianIndices(res)
    @inline function kernel(res, A, B)
        state = (blockIdx().x-1) * blockDim().x + threadIdx().x
        if state <= length(res)
            @inbounds i,j,b = CI[state].I
            i_A = (i-1) ÷ size(B,1) + 1
            j_A = (j-1) ÷ size(B,2) + 1
            i_B = (i-1) % size(B,1) + 1
            j_B = (j-1) % size(B,2) + 1
            (@inbounds res[state] = A[i_A, j_A, b]*B[i_B, j_B, b])
        end
        return
    end

    X, Y = cudiv(length(res))
    @cuda threads=X blocks=Y kernel(res, A, B)
    res
end

function join(reg1::GPUReg{B}, reg2::GPUReg{B}) where {B}
    s1 = reg1 |> rank3
    s2 = reg2 |> rank3
    state = YaoBase.batched_kron(s1, s2)
    ArrayReg{B}(copy(reshape(state, size(state, 1), :)))
end

export insert_qubits!
function insert_qubits!(reg::GPUReg{B}, loc::Int; nqubits::Int=1) where B
    na = nactive(reg)
    focus!(reg, 1:loc-1)
    reg2 = join(zero_state(nqubits; nbatch=B) |> cu, reg) |> relax! |> focus!((1:na+nqubits)...)
    reg.state = reg2.state
    reg
end

#=
for FUNC in [:measure!, :measure!]
    @eval function $FUNC(rng::AbstractRNG, op::AbstractBlock, reg::GPUReg, al::AllLocs; kwargs...) where B
        E, V = eigen!(mat(op) |> Matrix)
        ei = Eigen(E|>cu, V|>cu)
        $FUNC(rng::AbstractRNG, ei, reg, al; kwargs...)
    end
end
=#


text:
using CLIMA.VariableTemplates

import CLIMA.DGmethods: BalanceLaw, vars_aux, vars_state, vars_gradient,
                        vars_diffusive, flux_nondiffusive!, flux_diffusive!,
                        source!, wavespeed, boundary_state!,
                        gradvariables!,
                        diffusive!, init_aux!, init_state!,
                        init_ode_state, LocalGeometry

struct MMSModel{dim} <: BalanceLaw
end

vars_aux(::MMSModel,T) = @vars(x1::T, x2::T, x3::T)
vars_state(::MMSModel, T) = @vars(ρ::T, ρu::T, ρv::T, ρw::T, ρe::T)
vars_gradient(::MMSModel, T) = @vars(u::T, v::T, w::T)
vars_diffusive(::MMSModel, T) = @vars(τ11::T, τ22::T, τ33::T, τ12::T, τ13::T, τ23::T)

function flux_nondiffusive!(::MMSModel, flux::Grad, state::Vars,
                            auxstate::Vars, t::Real)
  # preflux
  T = eltype(flux)
  γ = T(γ_exact)
  ρinv = 1 / state.ρ
  u, v, w = ρinv * state.ρu, ρinv * state.ρv, ρinv * state.ρw
  P = (γ-1)*(state.ρe - ρinv * (state.ρu^2 + state.ρv^2 + state.ρw^2) / 2)

  # invisc terms
  flux.ρ  = SVector(state.ρu          , state.ρv          , state.ρw)
  flux.ρu = SVector(u * state.ρu  + P , v * state.ρu      , w * state.ρu)
  flux.ρv = SVector(u * state.ρv      , v * state.ρv + P  , w * state.ρv)
  flux.ρw = SVector(u * state.ρw      , v * state.ρw      , w * state.ρw + P)
  flux.ρe = SVector(u * (state.ρe + P), v * (state.ρe + P), w * (state.ρe + P))
end

function flux_diffusive!(::MMSModel, flux::Grad, state::Vars,
                         diffusive::Vars, auxstate::Vars, t::Real)
  ρinv = 1 / state.ρ
  u, v, w = ρinv * state.ρu, ρinv * state.ρv, ρinv * state.ρw

  # viscous terms
  flux.ρu -= SVector(diffusive.τ11, diffusive.τ12, diffusive.τ13)
  flux.ρv -= SVector(diffusive.τ12, diffusive.τ22, diffusive.τ23)
  flux.ρw -= SVector(diffusive.τ13, diffusive.τ23, diffusive.τ33)

  flux.ρe -= SVector(u * diffusive.τ11 + v * diffusive.τ12 + w * diffusive.τ13,
                     u * diffusive.τ12 + v * diffusive.τ22 + w * diffusive.τ23,
                     u * diffusive.τ13 + v * diffusive.τ23 + w * diffusive.τ33)
end

function gradvariables!(::MMSModel, transformstate::Vars, state::Vars, auxstate::Vars, t::Real)
  ρinv = 1 / state.ρ
  transformstate.u = ρinv * state.ρu
  transformstate.v = ρinv * state.ρv
  transformstate.w = ρinv * state.ρw
end

function diffusive!(::MMSModel, diffusive::Vars, ∇transform::Grad, state::Vars, auxstate::Vars, t::Real)
  T = eltype(diffusive)
  μ = T(μ_exact)

  dudx, dudy, dudz = ∇transform.u
  dvdx, dvdy, dvdz = ∇transform.v
  dwdx, dwdy, dwdz = ∇transform.w

  # strains
  ϵ11 = dudx
  ϵ22 = dvdy
  ϵ33 = dwdz
  ϵ12 = (dudy + dvdx) / 2
  ϵ13 = (dudz + dwdx) / 2
  ϵ23 = (dvdz + dwdy) / 2

  # deviatoric stresses
  diffusive.τ11 = 2μ * (ϵ11 - (ϵ11 + ϵ22 + ϵ33) / 3)
  diffusive.τ22 = 2μ * (ϵ22 - (ϵ11 + ϵ22 + ϵ33) / 3)
  diffusive.τ33 = 2μ * (ϵ33 - (ϵ11 + ϵ22 + ϵ33) / 3)
  diffusive.τ12 = 2μ * ϵ12
  diffusive.τ13 = 2μ * ϵ13
  diffusive.τ23 = 2μ * ϵ23
end

function source!(::MMSModel{dim}, source::Vars, state::Vars, aux::Vars, t::Real) where {dim}
  source.ρ  = Sρ_g(t, aux.x1, aux.x2, aux.x3, Val(dim))
  source.ρu = SU_g(t, aux.x1, aux.x2, aux.x3, Val(dim))
  source.ρv = SV_g(t, aux.x1, aux.x2, aux.x3, Val(dim))
  source.ρw = SW_g(t, aux.x1, aux.x2, aux.x3, Val(dim))
  source.ρe = SE_g(t, aux.x1, aux.x2, aux.x3, Val(dim))
end

function wavespeed(::MMSModel, nM, state::Vars, aux::Vars, t::Real)
  T = eltype(state)
  γ = T(γ_exact)
  ρinv = 1 / state.ρ
  u, v, w = ρinv * state.ρu, ρinv * state.ρv, ρinv * state.ρw
  P = (γ-1)*(state.ρe - ρinv * (state.ρu^2 + state.ρv^2 + state.ρw^2) / 2)
  return abs(nM[1] * u + nM[2] * v + nM[3] * w) + sqrt(ρinv * γ * P)
end

function boundary_state!(::Rusanov, bl::MMSModel, stateP::Vars, auxP::Vars, nM,
                         stateM::Vars, auxM::Vars, bctype, t, _...)
  init_state!(bl, stateP, auxP, (auxM.x1, auxM.x2, auxM.x3), t)
end

# FIXME: This is probably not right....
boundary_state!(::CentralGradPenalty, bl::MMSModel, _...) = nothing

function boundary_state!(::CentralNumericalFluxDiffusive, bl::MMSModel,
                         stateP::Vars, diffP::Vars, auxP::Vars, nM,
                         stateM::Vars, diffM::Vars, auxM::Vars, bctype, t, _...)
  init_state!(bl, stateP, auxP, (auxM.x1, auxM.x2, auxM.x3), t)
end

function init_aux!(::MMSModel, aux::Vars, g::LocalGeometry)
  x1,x2,x3 = g.coord
  aux.x1 = x1
  aux.x2 = x2
  aux.x3 = x3
end

function init_state!(bl::MMSModel{dim}, state::Vars, aux::Vars, (x1,x2,x3), t) where {dim}
  state.ρ = ρ_g(t, x1, x2, x3, Val(dim))
  state.ρu = U_g(t, x1, x2, x3, Val(dim))
  state.ρv = V_g(t, x1, x2, x3, Val(dim))
  state.ρw = W_g(t, x1, x2, x3, Val(dim))
  state.ρe = E_g(t, x1, x2, x3, Val(dim))
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
@time using HDF5
@time import InsarTimeseries, MapImages, Sario
repostr = expanduser("~/repos/")
@time include(repostr * "/InsarTimeseries.jl/scripts/plotting.jl")
@time include(repostr * "/InsarTimeseries.jl/scripts/find_abs_shift.jl")
@time include(repostr * "/InsarTimeseries.jl/scripts/point_analysis.jl")


text:
export invert, OptInterface, Opt_Optim, Opt_OptimPackNextGen

abstract type OptInterface end  # To accomodate multiple optimizers which are incompatible
struct Opt_Optim <: OptInterface end  #
struct Opt_OptimPackNextGen <: OptInterface end  #

"""
    invert(measured, rec0, forward; <keyword arguments>)

Tries to invert the `forward` model. `forward` is a function taking
an input with the shape of `rec0` and returns an object which has the
same shape as `measured`

Multiple keyword arguments can be specified for different loss functions,
regularizers and mappings.

# Arguments
- `loss=Poisson()`: the loss function being compatible to compare with `measured`.
- `regularizer=nothing`: A regularizer function, same form as `loss`.
- `λ=0.05`: A float indicating the total weighting of the regularizer with
    respect to the global loss function
- `mapping=Non_negative()`: Applies a mapping of the optimizer weight. Default is a
              parabola which achieves a non-negativity constraint.
- `iterations=10`: Specifies a number of iterations after the optimization.
    definitely should stop. Will be overwritten if `opt_options` is provided.
- `opt_package=Opt_Optim`: decides which backend for the optimizer is used.
- `opt=LBFGS()`: The chosen optimizer which must fit to `opt_package`.
- `opt_options=nothing`: Can be a options file required by Optim.jl. Will overwrite iterations.
- `debug_f=nothing`: A debug function which must take a single argument, the current reconstruction.
"""
function invert(measured, rec0, forward;
                iterations=10, λ=eltype(rec0)(0.05),
                regularizer=nothing,
                opt=LBFGS(linesearch=LineSearches.BackTracking()),
                opt_options=nothing,
                mapping=Non_negative(),
                loss=Poisson(),
                real_gradient=true,
                debug_f=nothing,
                opt_package=Opt_Optim)

    # if not special options are given, just restrict iterations
    if isa(opt_package, Opt_Optim) && opt_options == nothing
        opt_options = Optim.Options(iterations=iterations)
    end

    # Get the mapping functions to achieve constraints
    # like non negativity
    mf, m_invf = get_mapping(mapping)
    regularizer = get_regularizer(regularizer, eltype(rec0))

    debug_f_n(x) = let
        if isnothing(debug_f)
            identity(x)
        else
            debug_f(mf(x))
        end
    end

    storage_μ = deepcopy(measured)
    function total_loss(rec)
        # handle if there is a provided mapping function
        mf_rec = mf(rec)
        forward_v = forward(mf_rec)
        loss_v = sum(loss(forward_v, measured, storage_μ))
        loss_v += λ .* regularizer(mf_rec)
        return loss_v
    end
    # nice precompilation before calling Zygote etc.
    Base.invokelatest(total_loss, rec0)

    # see here https://github.com/FluxML/Zygote.jl/issues/342
    take_real_or_not(g) = real_gradient ? real.(g) : g
    take_real_or_not(g::AbstractArray{<:Real}) = g

    # this is the function which will be provided to Optimize
    # check Optim's documentation for the purpose of F and Get
    # but simply speaking F is the loss value and G it's gradient
    # depending whether one of them is nothing, we skip some computations
    # we need to call Base.invokelatest because the regularizer is a function
    # generated at runtime with eval.
    # This leads to the common "world age problem" in Julia
    # for more details on that check:
    # https://discourse.julialang.org/t/dynamically-create-a-function-initial-idea-with-eval-failed-due-to-world-age-issue/49139/17
    function fg!(F, G, rec)

        # Zygote calculates both derivative and loss, therefore do everything in one step
        if G != nothing
            # apply debug function
            debug_f_n(rec)

            y, back = Base.invokelatest(Zygote._pullback, total_loss, rec)
            # calculate gradient
            G .= take_real_or_not(Base.invokelatest(back, 1)[2])
            if F != nothing
                return y
            end
        end
        if F != nothing
            return Base.invokelatest(total_loss, rec)
        end
    end

    if isa(opt_package, Type{Opt_Optim})
        opt_options = Optim.Options(iterations=iterations)

        # do the optimization with LBGFS
        res = Optim.optimize(Optim.only_fg!(fg!), rec0, opt, opt_options)
        res_out = mf(Optim.minimizer(res))
    # supports a different interface as for example used in OptimPackNextGen for the function 'vmlmb!'
    elseif isa(opt_package, Type{Opt_OptimPackNextGen})
        res = copy(rec0)

        if isnothing(opt_options)
            opt((x,g) -> fg!(true, g, x), res; maxiter=iterations)
        else
            opt((x,g) -> fg!(true, g, x), res; maxiter=iterations, opt_options...)
        end
        res_out = mf(res)
    else
        error("Unknown optimizer interface $(typeof(opt_package))")
    end

    return res_out, res
end


text:
memphis_data_dir = joinpath(@__DIR__, "..", "data","memphis_house_location")
china_data_dir = joinpath(@__DIR__, "..", "data","covid19_china")


text:
## (c) 2013--2014 David A. van Leeuwen, (c) 2005--2012, Dan Ellis

## Freely adapted from Dan Ellis's rastamat package.
## We've kept routine names the same, but the interface has changed a bit.

## we haven't implemented rasta filtering, yet, in fact.  These routines are a minimum for
## encoding HTK-style mfccs

# powspec tested against octave with simple vectors

using DSP
using FFTW
using LinearAlgebra

function powspec(x::Vector{T}, sr::Real=8000.0; wintime=0.025, steptime=0.01, dither=true) where {T<:AbstractFloat}
    nwin = round(Integer, wintime * sr)
    nstep = round(Integer, steptime * sr)

    nfft = 2 .^ Integer((ceil(log2(nwin))))
    window = hamming(nwin)      # overrule default in specgram which is hanning in Octave
    noverlap = nwin - nstep

    y = spectrogram(x .* (1<<15), nwin, noverlap, nfft=nfft, fs=sr, window=window, onesided=true).power
    ## for compability with previous specgram method, remove the last frequency and scale
    y = y[1:end-1, :] ##  * sumabs2(window) * sr / 2
    y .+= dither * nwin / (sum(abs2, window) * sr / 2) ## OK with julia 0.5, 0.6 interpretation as broadcast!

    return y
end

# audspec tested against octave with simple vectors for all fbtypes
function audspec(x::Matrix{T}, sr::Real=16000.0; nfilts=ceil(Int, hz2bark(sr/2)), fbtype=:bark,
                 minfreq=0., maxfreq=sr/2, sumpower=true, bwidth=1.0) where {T<:AbstractFloat}
    nfreqs, nframes = size(x)
    nfft = 2(nfreqs-1)
    if fbtype == :bark
        wts = fft2barkmx(nfft, nfilts, sr=sr, width=bwidth, minfreq=minfreq, maxfreq=maxfreq)
    elseif fbtype == :mel
        wts = fft2melmx(nfft, nfilts, sr=sr, width=bwidth, minfreq=minfreq, maxfreq=maxfreq)
    elseif fbtype == :htkmel
        wts = fft2melmx(nfft, nfilts, sr=sr, width=bwidth, minfreq=minfreq, maxfreq=maxfreq,
                        htkmel=true, constamp=true)
    elseif fbtype == :fcmel
        wts = fft2melmx(nfft, nfilts, sr=sr, width=bwidth, minfreq=minfreq, maxfreq=maxfreq,
                        htkmel=true, constamp=false)
    else
        error("Unknown filterbank type ", fbtype)
    end
    wts = wts[:, 1:nfreqs]
    if sumpower
        return wts * x
    else
        return (wts * sqrt.(x)).^2
    end
end

function fft2barkmx(nfft::Int, nfilts::Int; sr=8000.0, width=1.0, minfreq=0., maxfreq=sr/2)
    hnfft = nfft >> 1
    minbark = hz2bark(minfreq)
    nyqbark = hz2bark(maxfreq) - minbark
    wts=zeros(nfilts, nfft)
    stepbark = nyqbark/(nfilts-1)
    binbarks = hz2bark.((0:hnfft) * sr / nfft)
    for i in 1:nfilts
        midbark = minbark + (i-1) * stepbark
        lof = (binbarks .- midbark) / width .- 0.5
        hif = (binbarks .- midbark) / width .+ 0.5
        logwts = min.(0, hif, -2.5lof)
        wts[i, 1:1+hnfft] = 10.0.^logwts
    end
    return wts
end

## Hynek's formula
hz2bark(f) = 6asinh(f / 600)
bark2hz(bark) = 600 * sinh(bark / 6)

function fft2melmx(nfft::Int, nfilts::Int; sr=8000.0, width=1.0, minfreq=0.0, maxfreq=sr/2, htkmel=false, constamp=false)
    wts = zeros(nfilts, nfft)
    # Center freqs of each DFT bin
    fftfreqs = collect(0:nfft-1) / nfft * sr;
    # 'Center freqs' of mel bands - uniformly spaced between limits
    minmel = hz2mel(minfreq, htkmel);
    maxmel = hz2mel(maxfreq, htkmel);
    binfreqs = mel2hz(minmel .+ collect(0:(nfilts+1)) / (nfilts + 1) * (maxmel - minmel), htkmel);
##    binbin = iround(binfrqs/sr*(nfft-1));

    for i in 1:nfilts
        fs = binfreqs[i .+ (0:2)]
        # scale by width
        fs = fs[2] .+ (fs .- fs[2])width
        # lower and upper slopes for all bins
        loslope = (fftfreqs .- fs[1]) / (fs[2] - fs[1])
        hislope = (fs[3] .- fftfreqs) / (fs[3] - fs[2])
        # then intersect them with each other and zero
        wts[i,:] = max.(0, min.(loslope, hislope))
    end

    if !constamp
        ## unclear what this does...
        ## Slaney-style mel is scaled to be approx constant E per channel
        wts = broadcast(*, 2 ./ ((binfreqs[3:nfilts+2]) - binfreqs[1:nfilts]), wts)
    end
    # Make sure 2nd half of DFT is zero
    wts[:, (nfft>>1)+1:nfft] .= 0.
    return wts
end

function hz2mel(f::Vector{T}, htk=false) where {T<:AbstractFloat}
    if htk
        return 2595 .* log10.(1 .+ f / 700)
    else
        f0 = 0.0
        fsp = 200/3
        brkfrq = 1000.0
        brkpt = (brkfrq - f0) / fsp
        logstep = exp(log(6.4) / 27)
        linpts = f .< brkfrq
        z = zeros(size(f))      # prevent InexactError() by making these Float64
        z[findall(linpts)] = f[findall(linpts)]/brkfrq ./ log(logstep)
        z[findall(.!linpts)] = brkpt .+ log.(f[findall(.!linpts)] / brkfrq) ./ log(logstep)
    end
    return z
end
hz2mel(f::AbstractFloat, htk=false)  = hz2mel([f], htk)[1]

function mel2hz(z::Vector{T}, htk=false) where {T<:AbstractFloat}
    if htk
        f = 700 .* (10 .^ (z ./ 2595) .- 1)
    else
        f0 = 0.0
        fsp = 200/3
        brkfrq = 1000.0
        brkpt = (brkfrq - f0) / fsp
        logstep = exp(log(6.4) / 27)
        linpts = z .< brkpt
        f = similar(z)
        f[linpts] = f0 .+ fsp * z[linpts]
        f[.!linpts] = brkfrq .* exp.(log.(logstep) * (z[.!linpts] .- brkpt))
    end
    return f
end

function postaud(x::Matrix{T}, fmax::Real, fbtype=:bark, broaden=false) where {T<:AbstractFloat}
    (nbands,nframes) = size(x)
    nfpts = nbands + 2broaden
    if fbtype == :bark
        bandcfhz = bark2hz.(range(0, hz2bark(fmax), length=nfpts))
    elseif fbtype == :mel
        bandcfhz = mel2hz(range(0, hz2mel(fmax), length=nfpts))
    elseif fbtype == :htkmel || fbtype == :fcmel
        bandcfhz = mel2hz(range(0, hz2mel(fmax,1), length=nfpts),1);
    else
        error("Unknown filterbank type")
    end
    # Remove extremal bands (the ones that will be duplicated)
    bandcfhz = bandcfhz[1+broaden:nfpts-broaden];
    # Hynek's magic equal-loudness-curve formula
    fsq = bandcfhz.^2
    ftmp = fsq + 1.6e5
    eql = ((fsq ./ ftmp).^2) .* ((fsq + 1.44e6) ./ (fsq + 9.61e6))
    # weight the critical bands
    z = broadcast(*, eql, x)
    # cube root compress
    z .^= 0.33
    # replicate first and last band (because they are unreliable as calculated)
    if broaden
        z = z[[1, 1:nbands, nbands], :];
    else
        z = z[[2, 2:(nbands-1), nbands-1], :]
    end
    return z
end

function dolpc(x::Array{T}, modelorder::Int=8) where {T<:AbstractFloat}
    nbands, nframes = size(x)
    r = real(ifft(vcat(x, x[collect(nbands-1:-1:2), :]), 1)[1:nbands, :])
    # Find LPC coeffs by durbin
    y, e = levinson(r, modelorder)
    # Normalize each poly by gain
    y ./= e
end

function lpc2cep(a::Array{T}, ncep::Int=0) where {T<:AbstractFloat}
    nlpc, nc = size(a)
    order = nlpc - 1
    if ncep == 0
        ncep = nlpc
    end
    c = zeros(ncep, nc)
    # Code copied from HSigP.c: LPC2Cepstrum
    # First cep is log(Error) from Durbin
    c[1,:] = -log.(a[1, :])
    # Renormalize lpc A coeffs
    a ./= a[1, :]
    for n in 2:ncep
        sum = zero(T)
        for m in 2:n
            sum += (n - m) * a[m, :] .* c[n - m + 1, :]
        end
        c[n, :] = -a[n, :] - sum / (n - 1)
    end
    return c
end

function spec2cep(spec::Array{T}, ncep::Int=13, dcttype::Int=2) where {T<:AbstractFloat}
    # no discrete cosine transform option
    dcttype == -1 && return log(spec)

    (nr, nc) = size(spec)
    dctm = zeros(typeof(spec[1]), ncep, nr)
    if 1 < dcttype < 4          # type 2,3
        for i in 1:ncep
            dctm[i, :] = cos.((i - 1) * collect(1:2:2nr-1)π / (2nr)) * √(2/nr)
        end
        if dcttype == 2
            dctm[1, :] /= √2
        end
    elseif dcttype == 4           # type 4
        for i in 1:ncep
            dctm[i, :] = 2cos.((i-1) * collect(1:nr)π/(nr+1))
            dctm[i, 1] += 1
            dctm[i, nr] += (-1)^(i-1)
        end
        dctm /= 2(nr + 1)
    else                        # type 1
        for i in 1:ncep
            dctm[i, :] = cos.((i-1) * collect(0:nr-1)π / (nr - 1)) / (nr - 1)
        end
        dctm[:, [1, nr]] /= 2
    end
    return dctm * log.(spec)
end

function lifter(x::Array{T}, lift::Real=0.6, invs=false) where {T<:AbstractFloat}
    (ncep, nf) = size(x)
    if lift == 0
        return x
    end
    if lift > 0
        if lift > 10
            error("Too high lift number")
        end
        liftw = [1; collect(1:ncep-1).^lift]
    else
        # Hack to support HTK liftering
        if !isa(lift, Integer)
            error("Negative lift must be interger")
        end
        lift = -lift            # strictly speaking unnecessary...
        liftw = vcat(1, (1 .+ lift/2 * sin.(collect(1:ncep-1)π / lift)))
    end
    if invs
        liftw = 1 ./ liftw
    end
    y = broadcast(*, x, liftw)
    return y
end

## Freely after octave's implementation, by Paul Kienzle <pkienzle@users.sf.net>
## Permission granted to usee this in a MIT license on 20 dec 2013 by the author Paul Kienzle:
## "You are welcome to move my octave code from GPL to MIT like core Julia."
## untested
## only returns a, v
function levinson(acf::Vector{T}, p::Int) where {T<:Real}
    if length(acf) < 1
        error("empty autocorrelation function")
    end
    if p < 0
        error("negative model order")
    end
    if p < 100
        ## direct solution [O(p^3), but no loops so slightly faster for small p]
        ## Kay & Marple Eqn (2.39)
        R = toeplitz(acf[1:p], conj(acf[1:p]))
        a = R \ -acf[2:p+1]
        pushfirst!(a, 1)
        v = real(a.*conj(acf[1:p+1]))
    else
        ## durbin-levinson [O(p^2), so significantly faster for large p]
        ## Kay & Marple Eqns (2.42-2.46)
        ref = zeros(p)
        g = -acf[2] / acf[1]
        a = [g]
        v = real((1 - abs2(g)) * acf[1])
        ref[1] = g
        for t=2:p
            g = -(acf[t+1] + a ⋅ acf[t:-1:2]) / v
            a = [a + g * conj(a[t-1:-1:1]), g]
            v *= 1 - abs2(g)
            ref[t] = g
        end
        pushfirst!(a, 1)
    end
    return (a, v)
end

function levinson(acf::Array{T}, p::Int) where {T<:Real}
    (nr,nc) = size(acf)
    a = zeros(p + 1, nc)
    v = zeros(p + 1, nc)
    for i in 1:nc
        a[:,i], v[:,i] = levinson(acf[:,i], p)
    end
    return (a, v)
end

## Freely after octave's implementation, ver 3.2.4, by jwe && jh
## skipped sparse implementation
function toeplitz(c::Vector{T}, r::Vector{T}=c) where {T<:Real}
    nc = length(r)
    nr = length(c)
    res = zeros(typeof(c[1]), nr, nc)
    if nc == 0 || nr == 0
        return res
    end
    if r[1] != c[1]
        ## warn
    end
    if typeof(c) <: Complex
        conj!(c)
        c[1] = conj(c[1])       # bug in julia?
    end
    ## if issparse(c) && ispsparse(r)
    data = [r[end:-1:2], c]
    for (i, start) in zip(1:nc, nc:-1:1)
        res[:,i] = data[start:start+nr-1]
    end
    return res
end


text:
using JuliaForBioinformatics
using Test

@testset "JuliaForBioinformatics.jl" begin

    @testset "Type tree" begin

        @test show_type_tree(Number) === nothing

        @test_throws ArgumentError show_type_tree(Any)
    end
end

max_stars_repo_name: diegozea/JuliaForBioinformatics

max_issues_repo_name: diegozea/JuliaForBioinformatics

max_forks_repo_name: diegozea/JuliaForBioinformatics


text:
"""
	EBayesTarget

Abstract type that describe Empirical Bayes estimands (which we want to estimate or conduct inference for).
"""
abstract type EBayesTarget end

# maybe have a distinction between 1D EBayes Target and 2D.

# functions so that EBayesTargets can play nicely with vectorization
broadcastable(target::EBayesTarget) = Ref(target)

"""
	extrema(target::EBayesTarget)

Returns a tuple `(a,b)` such that ` a ≦ target ≦ b` is always true. By default this
is ``(-\\infty,+\\infty)``, but e.g., for a `PriorTailProbability` it returns `(0,1)`.
"""
function extrema(target::EBayesTarget) #fallback
	(-Inf, Inf)
end

function pretty_label(target::EBayesTarget)
	""
end

function (targets::AbstractVector{<:EBayesTarget})(prior)
	[target(prior) for target in targets]
end

"""
	PosteriorEBayesTarget <: EBayesTarget

Abstract type for Empirical Bayes estimands that take the form ``E_G[h(\\mu) \\mid Z_i = z]`` for some function ``h``.
"""
abstract type PosteriorEBayesTarget <: EBayesTarget end

"""
	LinearEBayesTarget <: EBayesTarget

Abstract type for Empirical Bayes estimands that are linear functionals of the prior ``G``,
i.e., they take the form ``L(G)`` for some function linear functional ``L``.
"""
abstract type LinearEBayesTarget <: EBayesTarget end

"""
	cf(::LinearEBayesTarget, t)

The characteristic function of ``L(\\cdot)``, a `LinearEBayesTarget`, which we define as follows:

For ``L(\\cdot)`` which may be written as ``L(G) = \\int \\psi(\\mu)dG\\mu``
(for a measurable function ``\\psi``) this returns the Fourier transform of ``\\psi``
evaluated at t, i.e., ``\\psi^*(t) = \\int \\exp(it x)\\psi(x)dx``.

Note that ``\\psi^*(t)`` is such that for distributions ``G`` with density ``g``
(and ``g^*`` the Fourier Transform of ``g``) the following holds:

```math
L(G) = \\frac{1}{2\\pi}\\int g^*(\\mu)\\psi^*(\\mu) d\\mu
```
"""
function cf(::LinearEBayesTarget, t) end

location(target::LinearEBayesTarget) = target.Z

response(target::EBayesTarget) = _response(location(target))

function _response(Z)
	isa(Z, EBayesSample) ? response(Z) : Z
end

# general way of computing these for Mixtures.
function (target::MinimaxCalibratedEBayes.LinearEBayesTarget)(prior::MixtureModel)
    prior_probs = probs(prior)
    componentwise_target =[target(prior_j) for prior_j in components(prior)]
    dot(prior_probs, componentwise_target)
end

#--- MarginalDensityTarget

"""
	MarginalDensityTarget(Z::StandardNormalSample) <: LinearEBayesTarget

## Example call
```julia
MarginalDensityTaget(StandardNormalSample(2.0))
```
## Description
Describes the marginal density evaluated at ``Z=z``  (e.g. ``Z=2`` in the example above)
of a sample drawn from the hierarchical model
```math
\\mu \\sim G, Z \\sim \\mathcal{N}(0,1)
```
In other words, letting ``\\phi`` the Standard Normal pdf
```math
L(G) = \\phi \\star dG(z)
```
Note that `2.0` has to be wrapped inside `StandardNormalSample(2.0)` since this target
depends not only on `G` and the location, but also on the likelihood. Additional
likelihoods will be added in the future.
"""
struct MarginalDensityTarget{NS} <: LinearEBayesTarget
    Z::NS
end

pretty_label(target::MarginalDensityTarget) = L"f(x)"

function cf(target::MarginalDensityTarget{<:StandardNormalSample}, t)
    error_dbn = Normal(response(location(target))) #TODO...
    cf(error_dbn, t)
end

function riesz_representer(target::MarginalDensityTarget{<:StandardNormalSample}, t)
    error_dbn = Normal(response(location(target))) #TODO...
    pdf(error_dbn, t)
end

function extrema(target::MarginalDensityTarget{<:StandardNormalSample})
	(0, 1/sqrt(2π))
end

"""
	PriorDensityTarget(z::Float64) <: LinearEBayesTarget

## Example call
```julia
PriorDensityTarget(2.0)
```
## Description
This is the evaluation functional of the density of ``G`` at `z`, i.e.,
``L(G) = G'(z) = g(z)`` or in Julia code `L(G) = pdf(G, z)`.
"""
struct PriorDensityTarget <: LinearEBayesTarget
    x::Float64
end

pretty_label(target::PriorDensityTarget) = L"g(x)"

location(target::PriorDensityTarget) = target.x

function cf(target::PriorDensityTarget, t)
    exp(im*location(target)*t)
end

function (target::PriorDensityTarget)(prior)
    pdf(prior, location(target))
end

# hard coding this since almost everything else mixture related is handled as above.
function (target::PriorDensityTarget)(prior::MixtureModel)
	pdf(prior, location(target))
end

function extrema(target::PriorDensityTarget)
	(0, Inf)
end

#--- PriorTailProbability
"""
	PriorTailProbability(cutoff::Float64) <: LinearEBayesTarget

## Example call
```julia
PriorTailProbability(2.0)
```
## Description
This is the evaluation functional of the tail probability of ``G`` at `cutoff`, i.e.,
``L(G) = 1-G(\\text{cutoff})`` or in Julia code `L(G) = ccdf(G, z)`.
"""
Base.@kwdef struct PriorTailProbability <: LinearEBayesTarget
    cutoff::Float64 = 0.0
end

function extrema(target::PriorTailProbability)
	(0, 1)
end

function riesz_representer(target::PriorTailProbability, t)
    one(Float64)*(t >= target.cutoff)
end

function (target::PriorTailProbability)(prior)
	ccdf(prior, target.cutoff)
end

#function integration_domain(t::LinearInferenceTarget)
#end

"""
	LFSRNumerator(Z::StandardNormalSample) <: LinearEBayesTarget

## Example call
```julia
LFSRNumerator(StandardNormalSample(2.0))
```
## Description
Describes the linear functional
```math
L(G) =  \\int \\mathbf 1(\\mu \\geq 0) \\phi(z-\\mu) dG(\\mu)
```
This is used as an intermediate step in representing the local false sign rate, c.f. [`PosteriorTarget`](@ref)
and [`LFSR`](@ref).
"""
struct LFSRNumerator{NS <: StandardNormalSample} <: LinearEBayesTarget
    Z::NS
end

function cf(target::LFSRNumerator{<:StandardNormalSample}, t)
    x = response(target)
    exp(im*t*x- t^2/2)*(1+im*erfi((t-im*x)/sqrt(2)))/2
end

function riesz_representer(target::LFSRNumerator{<:StandardNormalSample}, t)
    x = response(target)
    pdf(Normal(), x - t)*(t>=0)
end

#--------- PosteriorMeanNumerator ---------------------------------

"""
	PosteriorMeanNumerator(Z::StandardNormalSample) <: LinearEBayesTarget

## Example call
```julia
PosteriorMeanNumerator(StandardNormalSample(2.0))
```
## Description
Describes the linear functional
```math
L(G) =  \\int \\mu \\phi(z-\\mu) dG(\\mu)
```
This is used as an intermediate step in representing the posterior mean, c.f. [`PosteriorTarget`](@ref)
and [`PosteriorMean`](@ref)
"""
struct PosteriorMeanNumerator{NS <: StandardNormalSample} <: LinearEBayesTarget
    Z::NS
end

function cf(target::PosteriorMeanNumerator{<:StandardNormalSample}, t)
    x = response(location(target))
    cf(Normal(x), t)*(x + im*t)
end

function riesz_representer(target::PosteriorMeanNumerator{<:StandardNormalSample}, t)
    x = response(location(target))
    pdf(Normal(), x - t)*t
end

#----------------------- Posterior Targets------------------------------------------
"""
	PosteriorTarget(num_target::LinearEBayesTarget) <: EBayesTarget

Type for Empirical Bayes estimands that take the form:

```math
E_G[h(\\mu) \\mid Z_i = z] = \\frac{ \\int h(\\mu) \\phi(z-\\mu) dG(\\mu)}{\\int \\phi(z-\\mu) dG(\\mu}
```
`num_target` is a `LinearEBayesTarget` that represents the numerator of the above estimands, i.e., of
``L(G)``.
"""
struct PosteriorTarget{NT <: LinearEBayesTarget} <: PosteriorEBayesTarget
    num_target::NT
end

"""
	PosteriorMean(Z)

## Example call
```julia
PosteriorMean(StandardNormalSample(2.0))
```
## Description
Shortcut for `PosteriorTarget(PosteriorMeanNumerator(Z))`, c.f. [`PosteriorMeanNumerator`](@ref)
"""
PosteriorMean(Z) = PosteriorTarget(PosteriorMeanNumerator(Z))

"""
	LFSR(Z)

## Example call
```julia
LFSR(StandardNormalSample(2.0))
```
## Description
Shortcut for `PosteriorTarget(LFSRNumerator(Z))`, c.f. [`LFSRNumerator`](@ref)
"""
LFSR(Z) = PosteriorTarget(LFSRNumerator(Z))

location(target::PosteriorTarget) = location(target.num_target)

function (post_target::PosteriorTarget)(prior)
    num = post_target.num_target(prior)
    denom = MarginalDensityTarget(location(post_target.num_target))(prior)
    num/denom
end

function MarginalDensityTarget(post_target::PosteriorTarget)
	MarginalDensityTarget(location(post_target.num_target))
end

function extrema(target::PosteriorTarget{LF} where LF<:LFSRNumerator)
	(0.0,1.0)
end

function pretty_label(target::PosteriorTarget{LF} where LF<:LFSRNumerator)
	L"\theta(x) := \Pr[\mu \geq 0 \mid X=x]"
end

function pretty_label(target::PosteriorTarget{PM} where PM<:PosteriorMeanNumerator)
	L"\theta(x) := E[\mu | X=x]"
end

#------------------ Calibration---------------------------------

#---------Should probably just turn this all into one function-----

Base.@kwdef struct CalibratedTarget{T <: PosteriorTarget, S<:Real} <: LinearEBayesTarget
    posterior_target::T
    θ̄::S #Pilot
	denom::S #denominator pilot
end

function MarginalDensityTarget(calib_target::CalibratedTarget)
	MarginalDensityTarget(calib_target.posterior_target)
end

# special semantics here.
function (calib_target::CalibratedTarget)(prior)
    numerator_target = calib_target.posterior_target.num_target
    denominator_target = MarginalDensityTarget(calib_target)
	θ̄ = calib_target.θ̄
    numerator_target(prior) - θ̄*denominator_target(prior)
end

# hardcoding
function (calib_target::CalibratedTarget)(prior::MixtureModel)
    numerator_target = calib_target.posterior_target.num_target
    denominator_target = MarginalDensityTarget(calib_target)
	θ̄ = calib_target.θ̄
    numerator_target(prior) - θ̄*denominator_target(prior)
end

#function (target::CalibratedTarget)(prior)
#	target.post.num_target(prior) -
#end

#function riesz_representer(target::CalibratedTarget, t)
#    x = target.num.x
#    θ̄ = target.θ̄
#    riesz_representer(target.num, t) - θ̄*riesz_representer(MarginalDensityTarget(x),t)
#end
#------


text:
"""
    DRO

A structure representing the DRO model.

Fields:
- `d::Constraint`:  type of the constraint
- `q::Vector`: probability distribution (vector)
- `c::Vector` vector of weights
- `ε::Real`:
- `cmin::Real`: minimum of the vector c
- `cmax::Real`: maximum of the vector c
- `Imax::Vector`: set of indexes, which fulfill c[i] == cmax
"""
struct DRO <: Model
    d::Constraint
    q::Vector
    c::Vector
    ε::Real
    cmin::Real
    cmax::Real
    Imax::Vector

    function DRO(d::Constraint, q::Vector{<:Real}, c::Vector{<:Real}, ε::Real)
        cmin = minimum(c)
        cmax = maximum(c)

        @assert length(q) == length(c)  "The length of the vector `q` and `c` must be the same."
        @assert all(q .>= 0)            "All components of the vector `q` must be nonnegative."
        @assert sum(q) ≈ 1 atol = 1e-10 "The `sum(q) = $(sum(q))` of the vector q must be equal to 1."
        @assert cmin != cmax            "The vector `c` must be non-constant."
        @assert ε > 0                   "The parameter `ε = $(ε)` must be greater than 0."

        return new(d, q, c, ε, cmin, cmax, findall(c .== cmax))
    end
end

"""
    Simplex

A structure representing the DRO model.

Fields:

"""
struct Simplex <: Model
    q::Vector
    lb::Vector
    ub::Vector

    function Simplex(q::Vector{<:Real}, lb::Vector{<:Real}, ub::Vector{<:Real})

        @assert sum(lb) <= 1   "The feasible set is empty: sum(lb) = $(sum(lb)) > 1"
        @assert sum(ub) >= 1   "The feasible set is empty: sum(ub) = $(sum(ub)) < 1"
        @assert all(ub .>= lb) "Upper bounds `ub` mus be must greater or equl to lower bounds"

        return new(q, lb, ub)
    end
end

max_stars_repo_name: VaclavMacha/Projections

max_issues_repo_name: VaclavMacha/Projections

max_forks_repo_name: VaclavMacha/Projections


text:
using Makie, MakieThemes, Colors
AbstractPlotting.set_theme!(ggthemr(:earth))

using SampledSignals, DSP, LinearAlgebra, PortAudio
using Observables

using FileIO: load, save, loadstreaming, savestreaming
import LibSndFile

function getTestSampleBuf()
  wav_dir = "../data/reaper-files"
  files = glob("*.wav", wav_dir)
  return load(files[1])
end

"Return a scene which visualizes the wave form of an n channel SampleBuf"
function visualize(samples::SampleBuf{T, N})::Scene where {T<:AbstractFloat, N<:Any}
  track_scenes = []
  for track in eachcol(samples)
    push!(track_scenes, lines(track; colormap=:bmw))
  end

  return vbox(track_scenes...)
end

function visualize(periodogram::DSP.Periodograms.Periodogram)::Scene
  return lines(freq(periodogram), power(periodogram))
end

"A live visualizer which updates the returned scene if the input buffer changes."
function visualize(live_sample_buffer::Observable{SampleBuf{T, N}})::Scene where {T<:AbstractFloat, N<:Any}
  n_tracks = size(live_sample_buffer.val)[2]
  track_scenes = []
  for iₜ in 1:n_tracks
    scene_for_track = lines(lift(buffer -> buffer[:, iₜ], live_sample_buffer); colormap=:bmw)
    push!(track_scenes, scene_for_track)

    axis = scene_for_track[Axis]
    axis[:names][:axisnames] = ("Sample #", "Amplitude")
  end

  return hbox(track_scenes...)
end

function visualize(spectrogram::DSP.Periodograms.Spectrogram)
  time_window_index = Node(1)
  scene = lines(freq(spectrogram), lift(tᵢ -> power(spectrogram)[:, tᵢ], time_window_index))
  axis = scene[Axis]
  axis[:names][:axisnames] = ("Frequency", "Power")
  xlims!(scene, (0, 5000))

  parent_scene = Scene(resolution = (1920, 1080))
  hbox(scene; parent = parent_scene)
  @async begin
    spectrogram_time = time(spectrogram)
      while !isopen(parent_scene) # wait for screen to be open
          sleep(0.01)
      end
      while isopen(parent_scene)
        print("-")
        for i in 1:length(spectrogram_time)
          time_window_index[] = i
          sleep(step(spectrogram_time))
        end
      end
  end
  return parent_scene
end


text:
module TestView

using ProfileVega
import VegaLite
using Test
using Profile

function profile_test(n)
    for i = 1:n
        A = randn(100,100,20)
        m = maximum(A)
        Am = mapslices(sum, A; dims=2)
        B = A[:,:,5]
        Bsort = mapslices(sort, B; dims=1)
        b = rand(100)
        C = B.*b
    end
end

ProfileVega.set_options(width = 800, height = 400)

@testset "ProfileVega View" begin

    p1 = @profview profile_test(2)

    @test p1 isa VegaLite.VLSpec
    @test p1.width == 800
    @test p1.height == 400

    ProfileVega.set_default_size(100, 200)
    p2 = @profview profile_test(2)

    @test p2.width == 100
    @test p2.height == 200

    p3 = @profview profile_test(2) width = 500
    @test p3.width == 500
    @test p3.height == 200

    p4 = @profview profile_test(2) height = 600
    @test p4.width == 100
    @test p4.height == 600

    p5 = @profview profile_test(2) height = 600 width = 400
    @test p5.width == 400
    @test p5.height == 600

    Profile.clear()
    @profile profile_test(10)
    p6 = ProfileVega.view()
    @test p6.width == 100
    @test p6.height == 200

    p7 = ProfileVega.view(width = 400, height = 500)
    @test p7.width == 400
    @test p7.height == 500
end

end # module

max_stars_repo_name: davidanthoff/ProfileVega

max_issues_repo_name: davidanthoff/ProfileVega

max_forks_repo_name: davidanthoff/ProfileVega


text:
abstract type AbstractContact{T} <: AbstractObject{T} end

"""
    mutable struct Contact{T, G, MT} <: AbstractContact{T}

Contact of a [`SolidStateDetector`](@ref).

For the simulation of the [`ElectricPotential`](@ref), all contacts are fixed to a constant potential value.

## Parametric types:
* `T`: Precision type.
* `G`: Type of `geometry`.
* `MT`: Type of `material`.

## Fields
* `potential::T`: Potential (in V) to which the contact will be fixed during the calculation of the [`ElectricPotential`](@ref).
* `material::MT`: Material of the contact.
* `id::Int`: Unique id that will unambiguously identify the contact.
* `name::String`: Custom name for the contact, relevant for plotting.
* `geometry::G`: Geometry of the contact, see [Constructive Solid Geometry (CSG)](@ref).

## Definition in Configuration File

A `Contact` is defined in the configuration file through an entry in the `contacts` array of a detector.
It needs `id`, `potential` and `geometry` and can optionally be given a `name` and `material`.

An example definition of contacts looks like this:
```yaml
contacts:
  - name: "n+ contact"
    id: 1
    potential: 5000V
    material: HPGe # optional
    geometry: # ....
  - name: "p+ contact"
    id: 2
    potential: 0
    material: HPGe #optional
    geometry: # ....
```
"""
struct Contact{T,G,MT} <: AbstractContact{T}
    potential::T
    material::MT
    id::Int
    name::String
    geometry::G
end

function Contact{T}(dict::Union{Dict{String,Any}, Dict{Any, Any}}, input_units::NamedTuple, outer_transformations)::Contact{T} where {T <: SSDFloat}
    id::Int = haskey(dict, "id") ? dict["id"] : -1
    material = haskey(dict, "material") ? material_properties[materials[dict["material"]]] : material_properties[materials["HPGe"]]
    name = haskey(dict,"name") ? dict["name"] : ""
    inner_transformations = parse_CSG_transformation(T, dict, input_units)
    transformations = combine_transformations(inner_transformations, outer_transformations)
    geometry = Geometry(T, dict["geometry"], input_units, transformations)
    return Contact( T(dict["potential"]), material, id, name, geometry )
end

function println(io::IO, d::Contact)
    println("\t________"*"Contact $(d.id)"*"________\n")
    println("\t---General Properties---")
    println("\t-Potential: \t\t $(d.potential) V")
    println("\t-Contact Material: \t $(d.material.name)")
    println()
end
print(io::IO, d::Contact{T}) where {T} = print(io, "Contact $(d.id) - $(d.potential) V - $(d.material.name)")
show(io::IO, d::Contact) = print(io, d)
show(io::IO,::MIME"text/plain", d::Contact) = show(io, d)


text:
# # System Independent Types
# @assert Cint == Int32
# @assert Cdouble == Float64

# Must match with AlgorithmParameters.h in HGS-CVRP: https://github.com/vidalt/HGS-CVRP
Base.@kwdef mutable struct AlgorithmParameters
    nbGranular :: Int32 = 20
    mu :: Int32 = 25
    lambda :: Int32 = 40
    nbElite :: Int32 = 4
    nbClose :: Int32 = 5
    targetFeasible :: Float64 = 0.2
    seed :: Int32 = 0
    nbIter :: Int32 = 20000
    timeLimit :: Float64 = 0.0
    useSwapStar :: Int32 = 1 # 1 = true, 0 = false
end

mutable struct C_SolutionRoute
    length :: Cint
    path :: Ptr{Cint}
end

mutable struct C_Solution
    cost :: Cdouble
    time :: Cdouble
    n_routes :: Cint
    routes :: Ptr{C_SolutionRoute}
end

mutable struct RoutingSolution
    cost :: Float64
    time :: Float64
    routes :: Vector{Vector{Int}}

    function RoutingSolution(sol_ptr::Ptr{C_Solution})
        sol = unsafe_load(sol_ptr)
        routes = Vector{Vector{Int}}(undef, sol.n_routes)
        for i in eachindex(routes)
            r = unsafe_load(sol.routes, i)
            path = unsafe_wrap(Array, r.path, r.length)
            routes[i] = unsafe_wrap(Array, r.path, r.length) .+ 1 # customer numbering offset
        end
        return new(sol.cost, sol.time, routes)
    end
end

function convert_destroy(c_sol_ptr::Ptr{C_Solution})
    # copy C structs to Julia structs
    j_sol = RoutingSolution(c_sol_ptr)

    # don't forget to delete pointers, important to avoid memory leaks
    ccall((:delete_solution, LIBHGSCVRP), Cvoid, (Ptr{C_Solution},), c_sol_ptr)

    return j_sol
end

"""
    function c_api_solve_cvrp(
        n::Integer,
        x::Vector,
        y::Vector,
        service_time::Vector,
        demand::Vector,
        vehicle_capacity::Float64,
        duration_limit::Float64,
        isRoundingInteger::Bool,
        isDurationConstraint::Bool,
        n_vehicles::Integer,
        parameters::AlgorithmParameters,
        verbose::Bool
    )

`x`, `y`, `service_time`, `demand` all require at least `n` elements.
Anything after `n` elements are ignored.

The first element should be the depot.
The next `n-1` elements should be customers.

"""
function c_api_solve_cvrp(
    n::Integer,
    x::Vector,
    y::Vector,
    service_time::Vector,
    demand::Vector,
    vehicle_capacity::Real,
    duration_limit::Real,
    isRoundingInteger::Bool,
    isDurationConstraint::Bool,
    n_vehicles::Integer,
    parameters::AlgorithmParameters,
    verbose::Bool
)

    @assert service_time[1] == 0.0
    @assert demand[1] == 0.0

    c_solution_ptr = ccall(
        (:solve_cvrp, LIBHGSCVRP),
        Ptr{C_Solution},
        (
            Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
            Cdouble, Cdouble, Cchar, Cchar,
            Cint, Ptr{AlgorithmParameters}, Cchar
        ),
        n, Cdouble.(x), Cdouble.(y), Cdouble.(service_time), Cdouble.(demand),
        vehicle_capacity, duration_limit, isRoundingInteger, isDurationConstraint,
        n_vehicles, Ref(parameters), verbose
    )

    return convert_destroy(c_solution_ptr)
end

function c_api_solve_cvrp_dist_mtx(
    n::Integer,
    x::Vector,
    y::Vector,
    dist_mtx::Matrix, # row-first matrix as in C
    service_time::Vector,
    demand::Vector,
    vehicle_capacity::Real,
    duration_limit::Real,
    isDurationConstraint::Bool,
    n_vehicles::Integer,
    parameters::AlgorithmParameters,
    verbose::Bool
)

    @assert service_time[1] == 0.0
    @assert demand[1] == 0.0

    if length(x) == length(y) == n
        x_ptr = Cdouble.(x)
        y_ptr = Cdouble.(y)
    else
        x_ptr = y_ptr = C_NULL
    end

    c_solution_ptr = ccall(
        (:solve_cvrp_dist_mtx, LIBHGSCVRP),
        Ptr{C_Solution},
        (
            Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
            Cdouble, Cdouble, Cchar,
            Cint, Ptr{AlgorithmParameters}, Cchar
        ),
        n, x_ptr, y_ptr, Cdouble.(dist_mtx), Cdouble.(service_time), Cdouble.(demand),
        vehicle_capacity, duration_limit, isDurationConstraint,
        n_vehicles, Ref(parameters), verbose
    )

    return convert_destroy(c_solution_ptr)
end


text:
using Test

@testset "AlgebraicNets" begin
  include("AlgebraicNets.jl")
end

@testset "MathFormulas.jl" begin
  include("MathFormulas.jl")
end


text:
module mstraight_duct_H8_1
using FinEtools
using Compat.Test
function test()
    t0  =  time()

    rho = 1.21*phun("kg/m^3");# mass density
    c  = 343.0*phun("m/s");# sound speed
    bulk =  c^2*rho;
    omega =  54.5901*phun("rev/s")
    vn0 =  -1.0*phun("m/s")
    Lx = 10.0*phun("m");# length of the box, millimeters
    Ly = 1.0*phun("m"); # length of the box, millimeters
    n = 20;#number of elements along the length

    # println("""
    #
    # Straight duct with anechoic termination.
    # Example from Boundary element acoustics: Fundamentals and computer codes, TW Wu, page 44.
    # Both real and imaginary components of the pressure should have amplitude of
    # rho*c = $(rho*c).
    #
    # Hexahedral mesh.
    # """)

    fens,fes  =  H8block(Lx,Ly,Ly,n,2,2); # Mesh
    bfes  =  meshboundary(fes)
    L0 = selectelem(fens,bfes, facing = true, direction = [-1.0 0.0 0.0])
    L10 = selectelem(fens,bfes, facing = true, direction = [+1.0 0.0 0.0])
    nLx = selectnode(fens,box = [0.0 Lx  0.0 0.0 0.0 0.0], inflate = Lx/1.0e5)

    geom  =  NodalField(fens.xyz)
    P  =  NodalField(zeros(Complex128,size(fens.xyz,1),1))

    numberdofs!(P)

    material = MatAcoustFluid(bulk,rho)
    femm  =  FEMMAcoust(IntegData(fes, GaussRule(3, 2)), material)

    S  =  acousticstiffness(femm, geom, P);
    C  =  acousticmass(femm, geom, P);

    E10femm  =  FEMMAcoustSurf(IntegData(subset(bfes,L10),GaussRule(2, 2)), material)
    D  =  acousticABC(E10femm, geom, P);

    E0femm  =  FEMMBase(IntegData(subset(bfes,L0), GaussRule(2,  2)))
    fi  =  ForceIntensity(-1.0im*omega*rho*vn0);
    F  =  distribloads(E0femm, geom, P, fi, 2);

    p = (-omega^2*S +omega*1.0im*D + C)\F
    scattersysvec!(P, p[:])

    # println("Pressure amplitude bounds: ")
    # println("  real $(minimum(real(P.values)))/$(maximum(real(P.values)))")
    # println("  imag $(minimum(imag(P.values)))/$(maximum(imag(P.values)))")
    #
    # println("Total time elapsed  =  ",time() - t0,"s")

    File  =   "straight_duct.vtk"
    scalars = real(P.values);
    vtkexportmesh(File, connasarray(fes), geom.values, FinEtools.MeshExportModule.H8;
    scalars = [("Pressure", scalars)])
    try rm(File) catch end
    # @async run(`"paraview.exe" $File`)

    ref=rho*c
    @test abs(minimum(real(P.values))-(-ref))/ref < 0.02
    @test abs(minimum(imag(P.values))-(-ref))/ref < 0.02
    @test abs(maximum(real(P.values))-(+ref))/ref < 0.02
    @test abs(maximum(imag(P.values))-(+ref))/ref < 0.02

    # plotly()
    # ix = sortperm(geom.values[nLx,1])
    # plot(geom.values[nLx,1][ix], real(P.values)[nLx][ix], color = :blue, label = "real")
    # plot!(geom.values[nLx,1][ix], imag(P.values)[nLx][ix], color = :red, label  =  "imag")
    # plot!(title = "Straight duct with anechoic termination",
    # xlabel = "x", ylabel = "Pressure")
    # gui()

    true

end
end
using .mstraight_duct_H8_1
mstraight_duct_H8_1.test()


text:
export CGPEvolution

import Cambrian.populate, Cambrian.evaluate

mutable struct CGPEvolution{T} <: Cambrian.AbstractEvolution
    config::NamedTuple
    logger::CambrianLogger
    population::Array{T}
    fitness::Function
    gen::Int
end

populate(e::CGPEvolution) = Cambrian.oneplus_populate(e)
evaluate(e::CGPEvolution) = Cambrian.fitness_evaluate(e, e.fitness)

function CGPEvolution(cfg::NamedTuple, fitness::Function;
                      logfile=string("logs/", cfg.id, ".csv"))
    logger = CambrianLogger(logfile)
    population = Cambrian.initialize(CGPInd, cfg)
    CGPEvolution(cfg, logger, population, fitness, 0)
end


text:
# ============================================================================ #
function plot_with_error(x::AbstractArray, y::Vector, yerr::Vector,
    col::AbstractString, ax=nothing; linewidth=4.0, args...)
    plot_with_error(x, y, yerr, parse(Colorant, col), ax, args...)
end
# ---------------------------------------------------------------------------- #
function plot_with_error(x::AbstractArray, y::Vector, yerr::Vector,
    col::ColorTypes.RGB, ax=nothing; linewidth=4.0, args...)
    if ax == nothing
        ax = default_axes()
    end
    col_array = [col.r, col.g, col.b]
    fcol, ecol = shading_color(col_array)
    ax.fill_between(x, y.-yerr, y.+yerr, facecolor=fcol, edgecolor=ecol)
    ax.plot(x, y, "-", color=col_array, linewidth=linewidth, args...)
end
# ---------------------------------------------------------------------------- #
function shading_color(col::Vector{T}) where {T<:Number}
    ferr = 8.0
    fedge = 0.25
    orig = convert(HSV, RGB(col...))
    hsv = HSV(orig.h, orig.s/ferr, 1.0 - (abs(1.0 - orig.s)^ferr/ferr))
    col_err = hsv2rgb(hsv)
    col_edge = (1.0 - fedge) * col_err + fedge * col
    return col_err, col_edge
end
# ---------------------------------------------------------------------------- #
function hsv2rgb(x::ColorTypes.HSV{T}) where {T<:Number}
    # xv = HSV(min([x.h, x.s, x.v], T(1.0))...)
    rgb = convert(RGB, x)
    return T[rgb.r, rgb.g, rgb.b]
end
# ============================================================================ #
plot_hist(x::AbstractVector) = _plot_hist(hist(x)...)
plot_hist(x::AbstractVector, n::Integer) = _plot_hist(hist(x, n)...)
plot_hist(x::AbstractVector, e::AbstractVector) = _plot_hist(hist(x, e)...)
# ============================================================================ #
edges2centers(edges::AbstractVector) = edges[1:end-1] .+ (step(edges)/2.0)
# ============================================================================ #
function _plot_hist(edges::AbstractVector, counts::AbstractVector)
    return qbar(edges2centers(edges), counts, step(edges))
end
# ============================================================================ #
function vline(x::Number; args...)
    plot([x, x], ylim(); args...)
end
# --------------------------------------------------------------------------- #
function vline(x::AbstractArray; args...)
    for item in x
        vline(item; args...)
    end
end
# ============================================================================ #
function hline(y::Number; args...)
    plot(xlim(), [y, y]; args...)
end
# --------------------------------------------------------------------------- #
function hline(y::AbstractArray; args...)
    for item in y
        hline(item; args...)
    end
end
# ============================================================================ #
function qplot(x, y)
    h = default_figure()
    p = plot(x, y, color="blue", linewidth=3.0)

    return h
end
# ============================================================================ #
function qplot(y)
    x = 1:length(y)
    h = qplot(x, y)
    return h
end
# ============================================================================ #
function qbar(x, y, width, col="black")
    h = default_figure()

    #width = absolute bar width
    # add error with yerr=...
    p = bar(x, y, width, color=col)

    return h
end
# ============================================================================ #
function qbar(y)
    x = 1:length(y)
    h = qbar(x-0.5, y, 1)
    return h
end
# ============================================================================ #
function closeall()
    plt.close()
end
# ============================================================================ #
default_figure() = default_figure(figure())
# ---------------------------------------------------------------------------- #
function default_figure(h::PyPlot.Figure)

    if !plt.fignum_exists(h.number)
        h = figure()
    end

    h.set_facecolor("white")
    ax = h.get_axes()
    if !isempty(ax)
        for x in ax
            delaxes(x)
        end
    end

    ax = h.add_axes(default_axes())

    return h
end
# ============================================================================ #
function default_axes(ax=nothing)
    if ax == nothing
        ax = PyPlot.axes()
    end
    #set ticks to face out
    ax.tick_params(direction="out", length=8.0, width=4.0)

    #turn off top and right axes
    ax.spines["right"].set_visible(false)
    ax.spines["top"].set_visible(false)

    #remove top and right tick marks
    tmp = ax.get_xaxis()
    tmp.tick_bottom()

    tmp = ax.get_yaxis()
    tmp.tick_left()

    ax.spines["left"].set_linewidth(4.0)
    ax.spines["bottom"].set_linewidth(4.0)

    return ax
end
# ============================================================================ #

# ============================================================================ #
# NOTES
# ============================================================================ #
#  #update figure
#  h.canvas.draw()
#
#  #close all figures
#  plt.close()


text:
# Convenience function for wavelength conversion
@inline aa_to_invum(wave::Real) = 10000 / wave

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

# Optical coefficients
const ccm89_ca = [1.0, 0.17699, -0.50447, -0.02427, 0.72085, 0.01979, -0.7753, 0.32999, 0.0]
const ccm89_cb = [0.0, 1.41338, 2.28305, 1.07233, -5.38434, -0.62251, 5.3026, -2.09002, 0.0]

const od94_ca = [1.0, 0.104, -0.609, 0.701, 1.137, -1.718, -0.827, 1.647, -0.505]
const od94_cb = [0.0, 1.952, 2.908, -3.989, -7.985, 11.102, 5.491, -10.805, 3.347]

"""
    CCM89(;Rv=3.1)

Clayton, Cardelli and Mathis (1989) dust law.

Returns E(B-V) in magnitudes at the given wavelength relative to the extinction
at 5494.5 Å. The default support is [1000, 33333]. Outside of that range this will return 0. `Rv` is the selective extinction and is valid over [2, 6]. A typical value for the Milky Way is 3.1.

# References
[Clayton,Cardelli and Mathis (1989)](https://ui.adsabs.harvard.edu/abs/1989ApJ...345..245C)
"""
@with_kw struct CCM89 <: ExtinctionLaw
    Rv::Float64 = 3.1
end

function (law::CCM89)(wave::T) where T
    checkbounds(law, wave) || return zero(float(T))
    x = aa_to_invum(wave)
    return ccm89_invum(x, law.Rv, ccm89_ca, ccm89_cb)
end

bounds(::Type{CCM89}) = (1000, 33333)

"""
    OD94(;Rv=3.1)

O'Donnell (1994) dust law.

This is identical to the Clayton, Cardelli and Mathis (1989) dust law, except
for different coefficients used in the optical (3030.3 Å to 9090.9 Å).

# References
[O'Donnell (1994)](https://ui.adsabs.harvard.edu/abs/1994ApJ...422..158O)

# See Also
[`CCM89`](@ref)
"""
@with_kw struct OD94 <: ExtinctionLaw
    Rv::Float64 = 3.1
end

function (law::OD94)(wave::T) where T
    checkbounds(law, wave) || return zero(float(T))
    x = aa_to_invum(wave)
    return ccm89_invum(x, law.Rv, od94_ca, od94_cb)
end

bounds(::Type{OD94}) = (1000, 33333)

"""
    DustExtinction.ccm89_invum(x, Rv, c_a, c_b)

The algorithm used for the [`CCM89`](@ref) extinction law, given inverse microns, Rv, and a set of coefficients for use in the optical (only difference between ccm89 and od94). For more information, seek the original paper.
"""
function ccm89_invum(x::Real, Rv::Real, c_a::Vector{<:Real}, c_b::Vector{<:Real})
    if x < 0.3
        error("out of bounds of CCM89, support is over $(bounds(CCM89)) angstrom")
    elseif x < 1.1  # Near IR
        y = x^1.61
        a = 0.574y
        b = -0.527y
    elseif x < 3.3  # Optical
        y = x - 1.82
        yn = 1.0
        a = @evalpoly y c_a[1] c_a[2] c_a[3] c_a[4] c_a[5] c_a[6] c_a[7] c_a[8] c_a[9]
        b = @evalpoly y c_b[1] c_b[2] c_b[3] c_b[4] c_b[5] c_b[6] c_b[7] c_b[8] c_b[9]
    elseif x < 8.0  # NUV
        a =  1.752 - 0.316x - (0.104 / ((x - 4.67)^2 + 0.341))
        b = -3.090 + 1.825x + (1.206 / ((x - 4.62)^2 + 0.263))
        if x > 5.9 # Far NUV
            y = x - 5.9
            a += @evalpoly y 0.0 0.0 -0.04473 -0.009779
            b += @evalpoly y 0.0 0.0 0.213 0.1207
        end
    elseif x ≤ 10.0 # FUV
        y = x - 8.0
        a = @evalpoly y -1.073 -0.628 0.137 -0.07
        b = @evalpoly y 13.67 4.257 -0.42 0.374
    else
        error("out of bounds of CCM89, support is over $(bounds(CCM89)) angstrom")
    end
    return a + b / Rv
end

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

"""
    CAL00(;Rv=4.05)

Calzetti et al. (2000) Dust Law.

Returns E(B-V) in magnitudes at the given wavelength. `λ` is the wavelength in Å
 and has support over [1200, 22000]. Outside of that range this will return 0.

Calzetti et al. (2000) developed a recipe for dereddening the spectra of
galaxies where massive stars dominate the radiation output. They found the best
 fit value for such galaxies was 4.05±0.80.

# References
[Calzetti et al. (2000)](https://ui.adsabs.harvard.edu/abs/2000ApJ...533..682C)
"""
@with_kw struct CAL00 <: ExtinctionLaw
    Rv::Float64 = 4.05
end
function (law::CAL00)(wave::T) where T
    checkbounds(law, wave) || return zero(float(T))
    x = aa_to_invum(wave)
    return cal00_invum(x, law.Rv)
end

bounds(::Type{CAL00}) = (1200, 22000)

"""
    DustExtinction.cal00_invum(x, Rv)

The algorithm used for the [`CAL00`](@ref) extinction law, given inverse microns and Rv. For more information, seek the original paper.
"""
function cal00_invum(x::Real, Rv::Real)
    if x > 1 / 0.12
        error("out of bounds of CAL00, support is over $(bounds(CAL00)) angstrom")
    elseif x > 1 / 0.63
        k = @evalpoly x -2.156 1.509 -0.198 0.011
    elseif x > 1 / 2.2
        k = @evalpoly x -1.857 1.040
    else
        error("out of bounds of CAL00, support is over $(bounds(CAL00)) angstrom")
    end
    return 1.0 + 2.659 * k / Rv
end

"""
    VCG04(;Rv=3.1)

Valencic, Clayton, & Gordon (2004) dust law.

This model applies to the UV spectral region all the way to 912 Å.
This model was not derived for the optical or NIR.

# References
[Valencic, Clayton, & Gordon (2004)](https://ui.adsabs.harvard.edu/abs/2004ApJ...616..912V/)
"""
@with_kw struct VCG04 <: ExtinctionLaw
    Rv::Float64 = 3.1
end

function (law::VCG04)(wave::T) where T
    checkbounds(law, wave) || return zero(float(T))
    x = aa_to_invum(wave)
    return vcg04_invum(x, law.Rv)
end

bounds(::Type{VCG04}) = (1250, 3030.3)

"""
    DustExtinction.vcg04_invum(x, Rv)

The algorithm used for the [`VCG04`](@ref) extinction law, given inverse microns and Rv. For more information, seek the original paper.
"""
function vcg04_invum(x::Real, Rv::Real)
    if 3.3 ≤ x ≤ 8.0  # NUV
        a = 1.808 - 0.215 * x - 0.134 / ((x - 4.558)^2 + 0.566)
        b = -2.350 + 1.403 * x + 1.103 / ((x - 4.587)^2 + 0.263)
    else
        error("out of bounds of VCG04, support is over $(bounds(VCG04)) angstrom")
    end
    if 5.9 ≤ x ≤ 8.0  # far-NUV
        y = x - 5.9
        a += @evalpoly y 0.0 0.0 -0.0077 -0.0030
        b += @evalpoly y 0.0 0.0 0.2060 0.0550
    end

    return a + b / Rv
end

"""
    GCC09(;Rv=3.1)

Gordon, Cartledge, & Clayton (2009) dust law.

This model applies to the UV spectral region all the way to 909.09 Å.
This model was not derived for the optical or NIR.

# References
[Gordon, Cartledge, & Clayton (2009)](https://ui.adsabs.harvard.edu/abs/2009ApJ...705.1320G/)
"""
@with_kw struct GCC09 <: ExtinctionLaw
    Rv::Float64 = 3.1
end

function (law::GCC09)(wave::T) where T
    checkbounds(law, wave) || return zero(float(T))
    x = aa_to_invum(wave)
    return gcc09_invum(x, law.Rv)
end

bounds(::Type{GCC09}) = (909.09, 3030.3)

"""
    DustExtinction.gcc09_invum(x, Rv)

The algorithm used for the [`GCC09`](@ref) extinction law, given inverse microns and Rv. For more information, seek the original paper.
"""
function gcc09_invum(x::Real, Rv::Real)
    if 3.3 ≤ x ≤ 11.0  # NUV
        a = 1.894 - 0.373 * x - 0.0101 / ((x - 4.57)^2 + 0.0384)
        b = -3.490 + 2.057 * x + 0.706 / ((x - 4.59)^2 + 0.169)
    else # out of bounds
        error("out of bounds of GCC09, support is over $(bounds(GCC09)) angstrom")
    end
    if 5.9 ≤ x ≤ 11.0  # far-NUV
        y = x - 5.9
        a += @evalpoly y 0.0 0.0 -0.110 -0.0100
        b += @evalpoly y 0.0 0.0 0.531 0.0544
    end

    return a + b / Rv
end


text:
"""
Submodule providing helper functions for generating plots.
"""
module Plotter
using ..CarrierCapture: potential, conf_coord

export plot_pots, plot_pot!, plot_ccs, plot_cc!

using Plots, LaTeXStrings

"""
    plot_pot!(pot; lplt_wf = false, plt = nothing, color = Nothing, label = "", scale_factor = 2e-2)

Plots `potential` with its' wave functions if `lplt_wf` is `true`.
"""
function plot_pot!(pot::potential; lplt_wf = false, plt = nothing, color = Nothing, label = "", scale_factor = 2e-2)
    if plt == nothing
        plt = plot()
    end
    label = if label == "" pot.name else label end
    color = if color == nothing "black" else color end

    # plot wave functions
    if lplt_wf
        ϵ = pot.ϵ; χ = pot.χ
        for i = 1:length(ϵ)
            plot!(plt, pot.Q, χ[i, :] * scale_factor .+ ϵ[i], lw = 0,
                  fillrange = [χ[i, :] * 0 .+ ϵ[i], χ[i, :] * scale_factor .+ ϵ[i]],
                  color = color, alpha = 0.5, label = "")
        end
    end

    # plot function
    plot!(plt, pot.Q, pot.E, lw = 4, color = color, label = label)

    # plot data
    if size(pot.QE_data)[1] > 1
        scatter!(plt, pot.QE_data.Q, pot.QE_data.E, color = color, label = "")
    end
    xaxis!(L"\ Q (amu^{1\/2} Å) \ (^{}$$"); yaxis!(L"\ Energy  (eV) \ (^{}$$")

    return plt
end

"""
    plot_cc!(cc; plt = nothing, color = nothing, label = "")

Plots capture coefficient `cc.capt_coeff` as a function of `1000/T`.
"""
function plot_cc!(cc::conf_coord; plt = nothing, color = nothing, label = "")
    if plt == nothing
        plt = plot()
    end
    label = if label == "" cc.name else label end
    # color = if color==Nothing "black" else color end
    plot!(plt, 1000 ./ cc.temperature, cc.capt_coeff, lw = 4, label = label)
    xaxis!(L"\ 1000\/T (K^{-1}) \ (^{}$$"); yaxis!(L"C (cm^{3}\/s) \ (^{}$$", :log10)
    return plt
end

"""
Depreciated:
Function for plotting `potential` in cli
"""
function plot_pots(pots::Dict{String,potential}, plot_cfg; output_fig = "potential.pdf")
    plt = plot(legend = :bottomleft)

    for (name, pot) in pots
        plot_pot!(pot, lplt_wf = true, plt = plt)
    end

    Emin = get(plot_cfg, "Emin", minimum(hcat([pot.E for (name, pot) in pots]...)))
    Emax = get(plot_cfg, "Emax", maximum(hcat([pot.E for (name, pot) in pots]...)))
    ylims!(Emin, Emax)

    Qmin = get(plot_cfg, "Qmin", minimum(hcat([pot.Q for (name, pot) in pots]...)))
    Qmax = get(plot_cfg, "Qmax", maximum(hcat([pot.Q for (name, pot) in pots]...)))
    xlims!(Qmin, Qmax)

    savefig(plt, output_fig)
end

"""
Depreciated:
Function for plotting `conf_coord` in cli
"""
function plot_ccs(ccs::Array{conf_coord,1}, plot_cfg::Dict; output_fig = "captcoeff.pdf")
    plt = plot(legend = :bottomleft)
    for cc in ccs
        plot_cc!(cc; plt = plt)
    end

    Cmin = get(plot_cfg, "Cmin", minimum(hcat([cc.capt_coeff for cc in ccs]...)) / 2.)
    Cmax = get(plot_cfg, "Cmax", maximum(hcat([cc.capt_coeff for cc in ccs]...)) * 2.)
    ylims!(Cmin, Cmax)

    invTmin = get(plot_cfg, "invTmin", 1000 / maximum(hcat([cc.temperature for cc in ccs]...)))
    invTmax = get(plot_cfg, "invTmax", 1000 / minimum(hcat([cc.temperature for cc in ccs]...)))
    xlims!(invTmin, invTmax)

    savefig(plt, output_fig)
end

end


text:
using AtomEnergyLevels, Printf

const eV = 27.211384

function slater_ts(Z)
    IPs = []
    for (l, subshell) in enumerate(atom[Z])
        for nᵣ in eachindex(subshell)
            config_ts = deepcopy(atom[Z])
            config_ts[l][nᵣ] -= 0.5
            results = lda(Z, conf = config_ts).orbitals
            push!(IPs, (results[(nᵣ = nᵣ - 1, l = l - 1)].ϵᵢ, l - 1, nᵣ - 1))
        end
    end
    @printf("\n== Electron binding energies of the %s atom, eV ==\n\n", keys(atom)[Z])
    for level in sort(IPs, by = x -> first(x))
        ϵᵢ, l, nᵣ = level
        n = nᵣ + l + 1
        @printf("\t%i%s\t%12.2f\n", n, shells[l + 1], ϵᵢ*eV)
    end
end

slater_ts(36)
#=
== Electron binding energies of the Kr atom, eV ==

        1S         -14097.00
        2S          -1838.84
        2P          -1676.35
        3S           -266.07
        3P           -205.35
        3D            -96.10
        4S            -27.95
        4P            -14.40
=#


text:
using StatsBase: sample

export pair_ring, pair_square
export rotor, merged_rotor, rotorset
export variational_circuit
export rand_single_gate, rand_gate, rand_circuit

################## Entangler ###################
"""
    pair_ring(n::Int) -> Vector

Pair ring.
"""
pair_ring(n::Int) = [i=>mod(i, n)+1 for i=1:n]

"""
    pair_square(m::Int, n::Int) -> Vector

Pair square.
"""
function pair_square(m::Int, n::Int; periodic=false)
    res = Vector{Pair{Int, Int}}(undef, (m-!periodic)*n+m*(n-!periodic))
    li = LinearIndices((m, n))
    k = 1
    for i = 1:2:m, j=1:n
        if periodic || i<m
            res[k] = li[i, j] => li[i%m+1, j]
            k+=1
        end
    end
    for i = 2:2:m, j=1:n
        if periodic || i<m
            res[k] = li[i, j] => li[i%m+1, j]
            k+=1
        end
    end
    for i = 1:m, j=1:2:n
        if periodic || j<n
            res[k] = li[i, j] => li[i, j%n+1]
            k+=1
        end
    end
    for i = 1:m, j=2:2:n
        if periodic || j<n
            res[k] = li[i, j] => li[i, j%n+1]
            k+=1
        end
    end
    res
end

###################### rotor and rotorset #####################
"""
    merged_rotor(noleading::Bool=false, notrailing::Bool=false) -> ChainBlock{1, ComplexF64}

Single qubit arbitrary rotation unit, set parameters notrailing, noleading true to remove trailing and leading Z gates.

!!! note

    Here, `merged` means `Rz(η)⋅Rx(θ)⋅Rz(ξ)` are multiplied first, this kind of operation if now allowed in differentiable
    circuit with back-propagation (`:BP`) mode (just because we are lazy to implement it!).
    But is a welcoming component in quantum differentiation.
"""
merged_rotor(noleading::Bool=false, notrailing::Bool=false) = noleading ? (notrailing ? Rx(0) : chain(Rx(0), Rz(0))) : (notrailing ? chain(Rz(0), Rx(0)) : chain(Rz(0), Rx(0), Rz(0)))

"""
    rotor(nbit::Int, ibit::Int, noleading::Bool=false, notrailing::Bool=false) -> ChainBlock{nbit, ComplexF64}

Arbitrary rotation unit (put in `nbit` space), set parameters notrailing, noleading true to remove trailing and leading Z gates.
"""
function rotor(nbit::Int, ibit::Int, noleading::Bool=false, notrailing::Bool=false)
    rt = chain(nbit, [put(nbit, ibit=>Rz(0.0)), put(nbit, ibit=>Rx(0.0)), put(nbit, ibit=>Rz(0.0))])
    noleading && popfirst!(rt)
    notrailing && pop!(rt)
    rt
end

rotorset(::Val{:Merged}, nbit::Int, noleading::Bool=false, notrailing::Bool=false) = chain(nbit, [put(nbit, j=>merged_rotor(noleading, notrailing)) for j=1:nbit])
rotorset(::Val{:Split}, nbit::Int, noleading::Bool=false, notrailing::Bool=false) = chain(nbit, [rotor(nbit, j, noleading, notrailing) for j=1:nbit])
rotorset(mode::Symbol, nbit::Int, noleading::Bool=false, notrailing::Bool=false) = rotorset(Val(mode), nbit, noleading, notrailing)

"""
    variational_circuit(nbit[, nlayer][, pairs]; mode=:Split, do_cache=false, entangler=cnot)

A kind of widely used differentiable quantum circuit, angles in the circuit is randomely initialized.

    * pairs: list of `Pair`s for entanglers in a layer, default to `pair_ring` structure,
    * mode: :Split or :Merged,
    * do_cache: cache the entangler matrix,
    * entangler: a constructor returns a two qubit gate, `f(n,i,j) -> gate`.
        The default value is `cnot(n,i,j)`.

ref:
    1. Kandala, A., Mezzacapo, A., Temme, K., Takita, M., Chow, J. M., & Gambetta, J. M. (2017).
       Hardware-efficient Quantum Optimizer for Small Molecules and Quantum Magnets. Nature Publishing Group, 549(7671), 242–246.
       https://doi.org/10.1038/nature23879.
"""
function variational_circuit(nbit, nlayer, pairs; mode=:Split, do_cache=false, entangler=cnot)
    circuit = chain(nbit)

    ent = chain(nbit, entangler(nbit, i, j) for (i, j) in pairs)
    if do_cache
        ent = ent |> cache
    end
    has_param = nparameters(ent) != 0
    for i = 1:(nlayer + 1)
        i!=1 && push!(circuit, has_param ? deepcopy(ent) : ent)
        push!(circuit, rotorset(mode, nbit, i==1, i==nlayer+1))
    end
    circuit
end

variational_circuit(n::Int; kwargs...) = variational_circuit(n, 3, pair_ring(n); kwargs...)

variational_circuit(nbit::Int, nlayer::Int; kwargs...) = variational_circuit(nbit, nlayer, pair_ring(nbit), kwargs...)

############### Completely random circuits (for testing and demo) ################
randlocs(nbit::Int, mbit::Int) = sample(1:nbit, mbit, replace=false)
const SINGLE_GATES = [X, Y, Z, H, Rx, Ry, Rz, shift, phase]

rand_single_gate(ngate::Int) = [rand_single_gates() for i=1:ngate]
function rand_single_gate()
    gate = rand(SINGLE_GATES)
    gate isa AbstractBlock ? gate : gate(rand()*2π)
end

"""
    rand_gate(nbit::Int, mbit::Int, [ngate::Int]) -> AbstractBlock

random nbit gate.
"""
rand_gate(nbit::Int, mbit::Int) = rand_gate(nbit, Val(mbit))
rand_gate(nbit::Int, mbit::Int, ngate::Int) = [rand_gate(nbit, mbit) for i=1:ngate]
rand_gate(nbit::Int, ::Val{1}) = put(nbit, rand(1:nbit)=>rand_single_gate())
function rand_gate(nbit::Int, ::Val{M}) where M
    locs = randlocs(nbit, M)
    control(nbit, locs[1:M-1], last(locs)=>rand_single_gate())
end

function rand_circuit(nbit::Int; p1 = (nbit==1 ? 1.0 : 0.66), ngate=5*nbit)
    c = chain(nbit)
    for i=1:ngate
        if rand() < p1
            push!(c, rand_gate(nbit, 1))
        else
            push!(c, rand_gate(nbit, 2))
        end
    end
    c
end


text:
function run!{T,NI}(method::Union{ODEString, StaticString}, E, dE, x0::Path{T,NI})
   # read all the parameters
   @unpack tol, maxtol, maxnit, precon_scheme, path_traverse, fixed_ends, verbose = method
   @unpack precon, precon_prep! = precon_scheme
   @unpack direction = path_traverse
   # initialise variables
   x = x0.x

   nit = 0
   numdE, numE = 0, 0
   log = PathLog()
   # and just start looping
   file =[]
   if verbose >= 4
       dt = Dates.format(now(), "d-m-yyyy_HH:MM")
       file = open("log_$(dt).txt", "w")
   end
   xout, log, alpha = odesolve(solver(method),
               (X, P, nit) -> forces(P, typeof(x0), X, dE, precon_scheme,
                                       direction(NI, nit), fixed_ends),
               vec(x), log; file = file,
               g = (X, P) -> redistribute(X, typeof(x0), P, precon_scheme),
               tol = tol, maxtol = maxtol, maxnit=maxnit,
               P = precon,
               precon_prep! = (P, X) -> precon_prep!(P, convert(typeof(x0), X)),
               method = "$(typeof(method))",
               verbose = verbose )

   x_return = verbose < 4 ? convert(typeof(x0), xout[end]) : [convert(typeof(x0), xout_n) for xout_n in xout]
   return x_return, log, alpha
end

function run!{T,NI}(method::AccelString, E, dE, ddE, x0::Path{T,NI})
   # read all the parameters
   @unpack tol, maxtol, maxnit, precon_scheme, path_traverse, fixed_ends, verbose = method
   @unpack precon, precon_prep! = precon_scheme
   @unpack direction = path_traverse
   # initialise variables
   x = x0.x

   nit = 0
   numdE, numE = 0, 0
   log = PathLog()
   # and just start looping
   file =[]
   if verbose >= 4
       dt = Dates.format(now(), "d-m-yyyy_HH:MM")
       file = open("log_$(dt).txt", "w")
   end

   xout, log, alpha = odesolve(solver(method),
               (X, P, nit) -> forces(P, typeof(x0), X, dE, precon_scheme,
                                       direction(NI, nit), fixed_ends),
               (X, P) -> jacobian(P, typeof(x0), X, dE, ddE),
               vec(x), log; file = file,
               g = (X, P) -> redistribute(X, typeof(x0), P, precon_scheme),
               tol = tol, maxtol = maxtol, maxnit = maxnit,
               P = precon,
               precon_prep! = (P, X) -> precon_prep!(P, convert(typeof(x0), X)),
               method = "$(typeof(method))",
               verbose = verbose )

   x_return = verbose < 4 ? convert(typeof(x0), xout[end]) : [convert(typeof(x0), xout_n) for xout_n in xout]
   return x_return, log, alpha
end

# forcing term for string method
function forces{T,NI}(precon, path_type::Type{Path{T,NI}}, X::Vector{Float64}, dE,
                  precon_scheme, direction, fixed_ends::Bool)

   x = convert(path_type, X)
   dxds = deepcopy(x)

   # preconditioner
   Np = size(precon, 1)
   P(i) = precon[mod(i-1,Np)+1, 1]
   P(i, j) = precon[mod(i-1,Np)+1, mod(j-1,Np)+1]

   # find tangents
   ds = [dist(precon_scheme, P, x, i) for i=1:length(x)-1]
   param = [0; [sum(ds[1:i]) for i in 1:length(ds)]]
   param /= param[end]; param[end] = 1.
   parametrise!(dxds, x, ds, parametrisation = param)
   dxds ./= point_norm(precon_scheme, P, dxds)
   dxds[1] = zeros(dxds[1]); dxds[end] = zeros(dxds[1])

   # potential gradient
   dE0_temp = []
   if !fixed_ends
      dE0_temp = [dE(x[i]) for i in direction]
      cost = length(param)
   else
      dE0_temp = [[zeros(x[1])]; [dE(x[i]) for i in direction[2:end-1]];
                  [zeros(x[1])]]
      cost = length(param) - 2
   end
   dE0 = [dE0_temp[i] for i in direction]

   # projecting out tangent term of potential gradient
   dE0⟂ = proj_grad(precon_scheme, P, dE0, dxds)

   # collecting force term
   F = forcing(precon_scheme, precon, dE0⟂)

   # residual error
   res = maxres(precon_scheme, P, dE0⟂)

   return F, res, cost, (X, Y) -> dot_P(precon_scheme, convert(path_type, X), P, convert(path_type, Y))
end

function jacobian{T,NI}(precon, path_type::Type{Path{T,NI}}, X::Vector{Float64},
   dE, ddE)

   x = convert(path_type, X)

   # preconditioner
   Np = size(precon, 1)
   P(i) = precon[mod(i-1,Np)+1, 1]
   P(i, j) = precon[mod(i-1,Np)+1, mod(j-1,Np)+1]

   hessian = ddE.precon; hessian_prep! = ddE.precon_prep!
   hessian = hessian_prep!(hessian, x)
   H(i) = hessian[mod(i-1,Np)+1, 1]
   H(i, j) = hessian[mod(i-1,Np)+1, mod(j-1,Np)+1]

   N = length(x); M = length(x[1])
   O = zeros(M, M); J = fill(O,(N, N))
   [J[n,n-1] = ∂Fⁿ⁻(x, n, dE, P) for n=2:N-1]
   if Np==1 && P(1)==I
      [J[n,n] = δFⁿ(x, n, H, P, H) for n=1:N]
   else
      [J[n,n] = δFⁿ(x, n, H, P, n -> I) for n=1:N]
   end
   [J[n,n+1] = ∂Fⁿ⁺(x, n, dE, P) for n=2:N-1]

   return ref(J)
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 81888440-4d42-11eb-1bf5-0f0108717ed6
using Markdown

# ╔═╡ ee1dacd0-52e6-11eb-3486-0bb6ee49c951
using StatsPlots

# ╔═╡ e4abc46e-52e6-11eb-19e0-19a94e1257c6
using Random, Distributions

# ╔═╡ 5cce7220-846c-11eb-0a99-7bed14a9f937
using DataFrames

# ╔═╡ 8aa500ee-52e5-11eb-034e-15d346c9433e
begin
	co2_mean = 1420 # g/kg residue
	co2_std = 240
	ch4_mean = 5.5
	ch4_std = 5.7
	n2o_mean = 0.09
	n2o_std = 0.04;

	ch4_co2_eq_lower = 25; ch4_co2_eq_upper = 28; #unitless
	n2o_co2_eq_lower = 265; n2o_co2_eq_upper = 298;

	rcr_lower = 1.4; rcr_upper = 1.8; # kg residue / kg rice
	billion_tonnes_to_million_tonnes = 10^3;
	india_annual_emission_billion_tonnes = 2.46;
	india_annual_emission_million_tonnes = india_annual_emission_billion_tonnes * billion_tonnes_to_million_tonnes;

	Random.seed!(123)
	n_samples = 100;

	co2_amt = rand(truncated(Normal(co2_mean, co2_std), 0, co2_mean * 3), n_samples); # g/kg residue
	ch4_amt = rand(truncated(Normal(ch4_mean, ch4_std), 0, ch4_mean * 3), n_samples);
	ch4_co2eq_factor = rand(Uniform(ch4_co2_eq_lower, ch4_co2_eq_upper), n_samples);
	ch4_co2eq = ch4_amt .* ch4_co2eq_factor;
	n2o_amt = rand(truncated(Normal(n2o_mean, n2o_std), 0, n2o_mean * 3), n_samples);
	n2o_co2eq_factor = rand(Uniform(n2o_co2_eq_lower, n2o_co2_eq_upper), n_samples);
	n2o_co2eq = n2o_amt .* n2o_co2eq_factor;

	residue_per_kg_rice = rand(Uniform(rcr_lower, rcr_upper), n_samples);

	co2eq = co2_amt + ch4_co2eq + n2o_co2eq # co2eq / kg residue
	co2eq_per_kg_rice = co2eq .* residue_per_kg_rice # co2eq / kg rice

	residue_amt = 14; # units of million tonnes
	residue_amt_3_lower = 15; residue_amt_3_upper = 22.5;
	residue_amt_3 = rand(Uniform(residue_amt_3_lower, residue_amt_3_upper), n_samples);
	residue_amt_2 = 23;

	total_co2 = co2_amt .* residue_amt / 1e3 # g/kg residue * million tonnes residue * 1kg/1000g = million tonnes
	total_co2eq = co2eq .* residue_amt / 1e3
	total_co2_2 = co2_amt .* residue_amt_2 / 1e3
	total_co2eq_2 = co2eq .* residue_amt_2 / 1e3
	total_co2_3 = co2_amt .* residue_amt_3 / 1e3
	total_co2eq_3 = co2eq .* residue_amt_3 / 1e3

	co2_df = DataFrame(Amount = vcat(total_co2, total_co2_2, total_co2_3),
	Case = vcat(repeat([1], n_samples),
		repeat([2], n_samples), repeat([3], n_samples)),
	Quantity = repeat(["CO2 only"], n_samples * 3));

	total_df = DataFrame(Amount = vcat(total_co2eq, total_co2eq_2, total_co2eq_3),
	Case = vcat(repeat([1], n_samples),
		repeat([2], n_samples), repeat([3], n_samples)),
	Quantity = repeat(["Total"], n_samples * 3));
end;

# ╔═╡ e986fbd0-4d42-11eb-3342-ff3bb8eff1ad
md"# What is the global warming impact of crop stubble burning in Northern India?

By: Srikiran Chandrasekaran, Dileep Kishore

To clear the fields of the residue from the previous rice crop and make way for the subsequent wheat crop, some farmers in Northern India burn what's left in the field. Called crop residue burning or stubble burning, this happens in the months of October and November [1]. Manually clearing the residue, or using machines built for the purpose, may cost more [10].

We were interested in understanding the global warming impact of stubble burning. It contributes to global warming in at least two ways - through the release of greenhouse gases (GHGs) and the obvious one - the heat from burning the stubble. We were interested in the former.

## How is greenhouse gas global warming potential measured?

GHGs like carbon dioxide ($CO_2$) and methane ($CH_4$) have different propensities to cause global warming. Hence global warming potential is measured in CO2eq or, $CO_2$ equivalents. That is, any contribution from say, methane, is understood as the contribution coming from an equivalent amount of $CO_2$.

For example, suppose a process releases 10g of $CO_2$ and 2g of $CH_4$. And suppose that the global warming propensity of $CH_4$ is 3 times that of $CO_2$. Then this process contributes $10 + 2 \times 3 = 16$g CO2eq.

## Which GHGs, and how much?

Before understanding the overall global warming effects, its useful to know which GHGs are released when burning crop stubble. Here we looked at carbon dioxide, methane and nitrous oxide ($N_2O$). While the latter two are addressed in the analysis by the FAO [3], $CO_2$ was not.

To understand how much of each of these gases are produced per kg of crop residue, we looked at the meta analysis by Meinrat Andreae [2] and found the following

Gas | Mean emission (g/kg residue) | Standard deviation emission (g/kg residue) | Number of samples on which esitmate is based
---|---|---|---
$CO_2$ | $(co2_mean) | $(co2_std) | 25
$CH_4$ | $(ch4_mean) | $(ch4_std) | 17
$N_2O$ | $(n2o_mean) | $(n2o_std) | 5

And we also have the CO2eq data for the two gases from the following sources

Gas | CO2eq (g $CO_2$ / g gas) | Source
---|---|---
$CH_4$ | $(ch4_co2_eq_upper) | [4]
$N_2O$ | $(n2o_co2_eq_lower) | \"
$CH_4$ | $(ch4_co2_eq_lower) | [5]
$N_2O$ | $(n2o_co2_eq_upper) | \"

### Global warming impact per kg of residue

The numbers above are sufficient to understand the global warming impact of burning 1kg of crop residue. We assumed a normal distribution for the gas amounts (since we have mean and std) and a uniform distribution for the CO2eq factors (since we have two different values).

"

# ╔═╡ d6b15000-7f83-11eb-1fca-478cddbecf54
begin
	density([log10.(co2_amt), log10.(ch4_amt), log10.(n2o_amt)], label=["Carbon Dioxide" "Methane" "Nitrous Oxide"], width=2, ylabel="Density", xlabel="Log10(GHG emission) (g/kg residue)", title="Distribution of GHG emission/kg residue burnt", fill=true, legend=:topleft, alpha=0.5)
end

# ╔═╡ f4406d20-586b-11eb-0d35-51da6f73592e
md"Since the X axis is on the log scale, we observe that $CO_2$ emission is $\approx$2 orders of magnitude greater than $CH_4$ emission, which is an order greater than $N_2O$ emission. While these are simply the amounts, the impact on global warming is better reflected by the chart below."

# ╔═╡ 55efed92-7f84-11eb-0c2e-cd38094e461d
begin
	# gr(size=(2500,3000))
	density([log10.(co2_amt), log10.(ch4_co2eq), log10.(n2o_co2eq), log10.(co2eq)], width=2, ylabel="Density", xlabel="Log10(CO2eq) (g/kg residue)", title="Distribution of GHG emission gCO2eq/kg residue burnt", legend=:topleft, fill=true, alpha=0.5, label=["Carbon Dioxide" "Methane" "Nitrous Oxide" "Total"])
end

# ╔═╡ 925ff0a0-586e-11eb-24c5-3be428cb076f
md"

So it does appear that the impact of methane and nitrous oxide is higher than surmised purely based on their amounts. However it is only a fraction of the impact of $CO_2$ alone, based on the chosen data.

### Putting these numbers in perspective

To put these numbers in perspective, we looked at the $CO_2$ output of growing rice.

- According to a publication [7] and the subsequent analysis by Our World in Data [8], growing one kg of *rice* creates approximately 4kg $CO_2$. Note that this estimate does not talk about the residue.
- In comparison, CO2 from burning residue is : $(round(mean(co2eq) / 1000, digits=2)) $$\pm$$ $(round(std(co2eq) / 1000, digits=2)) (mean $$\pm$$ SD) kg CO2eq for one kg of *residue* (from the plot above).
- Using a conversion factor of $(rcr_lower) - $(rcr_upper) kg residue per kg rice crop [1], we have to $(round(mean(co2eq_per_kg_rice) / 1000, digits=2)) $$\pm$$ $(round(std(co2eq_per_kg_rice) / 1000, digits=2)) kg CO2eq per kg of rice grown (mean $$\pm$$ SD), from stubble burning.

It appears that growing rice with stubble burning creates $\approx$50% more CO2 ($(round(mean(co2eq_per_kg_rice) / 1000, digits=2)) kg compared to 4kg $$CO_2$$),than growing rice without stubble burning!

## Overall impact on global warming

To understand the overall impact of stubble burning on global warming, we found estimates for the amount of agricultural residue that was burned in 2016.

Case | Amount of residue (million tonnes)/ year | Source
---|---|---
1 | $(residue_amt) | [2]
2 | $(residue_amt_2) | [3]
3 | $(residue_amt_3_lower) - $(residue_amt_3_upper) | [6]

Computing the impact of these residue amounts is simply multiplying one or more of these residue amounts with the emission per kg from the previous section. This gives

"

# ╔═╡ 745cd7d0-846f-11eb-1cab-e1ec4e7d7fad
begin
	@df co2_df violin(:Case, :Amount, group = :Quantity, side=:left,
		label="CO2 only", xticks=[1, 2, 3], xlabel="Cases", ylabel="CO2eq (million tonnes)", title="Distribution of gCO2eq from crop residue burning")
	@df total_df violin!(:Case, :Amount, group = :Quantity, side=:right, label="Total")
end

# ╔═╡ c2cb8850-86c5-11eb-2c30-8168e9df6a5a
md"We trust the estimate of of case number 1 most, since the authors of that publication follow a first principled approach to estimate the amount of residue, and use many data sources to avoid underestimations. For example, they fill in the gaps from satellite data with surveys.

### Putting these numbers in perspective

According to Our World in Data [9], the total CO2 emission of India is $india_annual_emission_billion_tonnes billion tonnes in 2017. The CO2eq from case number 1 as a fraction of this quantity evaluates to $(round(mean(total_co2eq) / (india_annual_emission_million_tonnes) * 100, digits=2)) $$\pm$$
$(round(std(total_co2eq) / (india_annual_emission_million_tonnes) * 100, digits=2)) % (mean $$\pm$$ SD) or almost 1%! This is a lower estimate than case 2, where we are looking at $(round(mean(total_co2eq_2) / (india_annual_emission_million_tonnes) * 100, digits=2)) % $$\pm$$ $(round(std(total_co2eq_2) / (india_annual_emission_million_tonnes) * 100, digits=2)) % (mean $$\pm$$ SD).

## Summary of findings

1. The bulk of GHG emissions from stubble burning is from the direct contribution of $CO_2$, with small contributions from nitrous oxide and methane.
1. Compared to growing rice without stubble burning, growing rice with burning of crop residue releases $\approx$50% more CO2eq of GHG.
1. Accounting for how much rice is grown with stubble burning, the GHG contribution of stubble burning was $\approx$ 1% of the total GHG emissions of India in 2017.

> Note, we are not experts in this field. So results to be interpreted with caution.

"

# ╔═╡ 701820f0-52e5-11eb-1cff-195defc1e22a
md"

## Links

[1]: <https://doi.org/10.1016/j.aeaoa.2020.100091>

[2]: <https://doi.org/10.5194/acp-19-8523-2019>

[3]: <http://www.fao.org/faostat/en/#data/GB/metadata>

[4]: <https://cdiac.ess-dive.lbl.gov/pns/current_ghg.html>

[5]: <https://climatechangeconnection.org/emissions/co2-equivalents/>

[6]: <https://doi.org/10.1038/s41598-019-52799-x>

[7]: <https://doi.org/10.1126/science.aaq0216>

[8]: <https://ourworldindata.org/environmental-impacts-of-food>

[9]: <https://ourworldindata.org/co2/country/india?country=~IND>

[10]: <https://indianexpress.com/article/india/stubble-burning-punjab-farmers-amarinder-singh-ngt-air-pollution-4897240/>"

# ╔═╡ 8685b4b0-52e5-11eb-1414-bfb2439d473c

# ╔═╡ 77ed9760-52e5-11eb-3534-d512d858bd9d

# ╔═╡ Cell order:
# ╠═81888440-4d42-11eb-1bf5-0f0108717ed6
# ╠═ee1dacd0-52e6-11eb-3486-0bb6ee49c951
# ╠═e4abc46e-52e6-11eb-19e0-19a94e1257c6
# ╠═5cce7220-846c-11eb-0a99-7bed14a9f937
# ╟─e986fbd0-4d42-11eb-3342-ff3bb8eff1ad
# ╟─8aa500ee-52e5-11eb-034e-15d346c9433e
# ╟─d6b15000-7f83-11eb-1fca-478cddbecf54
# ╟─f4406d20-586b-11eb-0d35-51da6f73592e
# ╟─55efed92-7f84-11eb-0c2e-cd38094e461d
# ╟─925ff0a0-586e-11eb-24c5-3be428cb076f
# ╟─745cd7d0-846f-11eb-1cab-e1ec4e7d7fad
# ╟─c2cb8850-86c5-11eb-2c30-8168e9df6a5a
# ╟─701820f0-52e5-11eb-1cff-195defc1e22a
# ╟─8685b4b0-52e5-11eb-1414-bfb2439d473c
# ╟─77ed9760-52e5-11eb-3534-d512d858bd9d

max_stars_repo_name: quantumbrake/ghg-dashboard

max_issues_repo_name: quantumbrake/ghg-dashboard

max_forks_repo_name: quantumbrake/ghg-dashboard


text:
function iterate!(job::BasicMCJob, ::Type{RAM}, ::Type{Univariate})
  job.sstate.count += 1

  if job.tuner.verbose
    job.sstate.tune.proposed += 1
  end

  job.sstate.randnsample = randn()
  job.sstate.pstate.value = job.pstate.value+job.sstate.S*job.sstate.randnsample

  job.parameter.logtarget!(job.sstate.pstate)

  job.sstate.ratio = job.sstate.pstate.logtarget-job.pstate.logtarget

  if (job.sstate.ratio > 0 || (job.sstate.ratio > log(rand())))
    job.pstate.value = job.sstate.pstate.value

    job.pstate.logtarget = job.sstate.pstate.logtarget
    if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
      job.pstate.loglikelihood = job.sstate.pstate.loglikelihood
    end
    if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
      job.pstate.logprior = job.sstate.pstate.logprior
    end

    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = true
      end
    end

    if job.tuner.verbose
      job.sstate.tune.accepted += 1
    end
  else
    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = false
      end
    end
  end

  if job.tuner.verbose
    if job.sstate.tune.totproposed <= job.range.burnin && mod(job.sstate.tune.proposed, job.tuner.period) == 0
      rate!(job.sstate.tune)
      println(
        "Burnin iteration ",
        job.fmt_iter(job.sstate.tune.totproposed),
        " of ",
        job.range.burnin,
        ": ",
        job.fmt_perc(100*job.sstate.tune.rate),
        " % acceptance rate"
      )
      reset_burnin!(job.sstate.tune)
    end
  end

  job.sstate.η = min(1, job.sstate.count^(-job.sampler.γ))
  job.sstate.SST = job.sstate.η*(min(1, exp(job.sstate.ratio))-job.sampler.targetrate)
  job.sstate.SST = abs2(job.sstate.S)*(1+job.sstate.SST)
  job.sstate.S = chol(job.sstate.SST)
end

function iterate!(job::BasicMCJob, ::Type{RAM}, ::Type{Multivariate})
  job.sstate.count += 1

  if job.tuner.verbose
    job.sstate.tune.proposed += 1
  end

  job.sstate.randnsample[:] = randn(job.pstate.size)
  job.sstate.pstate.value[:] = job.pstate.value+job.sstate.S*job.sstate.randnsample

  job.parameter.logtarget!(job.sstate.pstate)

  job.sstate.ratio = job.sstate.pstate.logtarget-job.pstate.logtarget

  if (job.sstate.ratio > 0 || (job.sstate.ratio > log(rand())))
    job.pstate.value = copy(job.sstate.pstate.value)

    job.pstate.logtarget = job.sstate.pstate.logtarget
    if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
      job.pstate.loglikelihood = job.sstate.pstate.loglikelihood
    end
    if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
      job.pstate.logprior = job.sstate.pstate.logprior
    end

    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = true
      end
    end

    if job.tuner.verbose
      job.sstate.tune.accepted += 1
    end
  else
    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = false
      end
    end
  end

  if job.tuner.verbose
    if job.sstate.tune.totproposed <= job.range.burnin && mod(job.sstate.tune.proposed, job.tuner.period) == 0
      rate!(job.sstate.tune)
      println(
        "Burnin iteration ",
        job.fmt_iter(job.sstate.tune.totproposed),
        " of ",
        job.range.burnin,
        ": ",
        job.fmt_perc(100*job.sstate.tune.rate),
        " % acceptance rate"
      )
      reset_burnin!(job.sstate.tune)
    end
  end

  job.sstate.η = min(1, job.pstate.size*job.sstate.count^(-job.sampler.γ))
  job.sstate.SST[:, :] = (
    job.sstate.randnsample*job.sstate.randnsample'/dot(job.sstate.randnsample, job.sstate.randnsample)*
    job.sstate.η*(min(1, exp(job.sstate.ratio))-job.sampler.targetrate)
  )
  job.sstate.SST[:, :] = job.sstate.S*(eye(job.pstate.size)+job.sstate.SST)*job.sstate.S'
  job.sstate.S[:, :] = ctranspose(chol(Hermitian(job.sstate.SST)))
end


text:
"""
    get_powers(x, y, powers)

Obtain a DataFrame of the powers and cross products of x and y

# Arguments

- x = 1d array of floats
- y = 1d array of floats of same length as `x`
- powers = 1d array of integers indicating the powers to include. Defaults to
  [1, 2, 3, 4] to get everything up to the quadric term

# Value

DataFrame of the powers and cross products of `x` and `y`, as specified by the
powers argument.

"""
function get_powers(x::Array{Float64,1}, y::Array{Float64,1},
                    powers::Array{Int64,1}=[1, 2, 3, 4])

    # Check that x and y are the same length
    @assert length(x) == length(y)

    # Initialze DataFrame to return
    out = DataFrame()

    # Get powers of x and y
    for i in 1:length(powers)
        out[!, Symbol(string("X", powers[i]))] = x.^powers[i]
            out[!, Symbol(string("Y", powers[i]))] = y.^powers[i]
    end

    # Get cross products
    for i in 1:(length(powers)-1)
        for j in 1:(length(powers)-i)
            out[!, Symbol(string("X", powers[i], "Y", powers[j]))] =
                out[:, Symbol(string("X", powers[i]))] .*
                out[:, Symbol(string("Y", powers[j]))]
        end
    end

    return out
end