text:
# This file was generated, do not modify it. # hide
subdf_setosa = gdf[1]
describe(subdf_setosa, :min, :mean, :max)


text:
abstract type Preconditioner end

mutable struct HMCOperatorXMLElement <: OperatorXMLElement
    el::XMLOrNothing
    param_provider::MyXMLElement
    grads::Array{<:MyXMLElement}
    already_made::Vector{Bool}
    weight::Real
    n_steps::Int
    step_size::Float64
    draw_variance::Float64
    auto_optimize::Bool
    check_grad::Int
    grad_tolerance::Float64
    geodesic::Bool
    transform::String
    mask::Vector{<:Real}
    preconditioning::Preconditioner

    function HMCOperatorXMLElement(param_provider::MyXMLElement,
                        grads::Array{<:MyXMLElement};
                        geodesic::Bool = false,
                        already_made::Vector{Bool} = fill(false, length(grads)),
                        transform::String = "")
        return new(nothing, param_provider, grads, already_made,
                    1.0, 5, 0.05, 1.0, true, 0, 1e-3, geodesic, transform,
                    Float64[],
                    NoPreconditioning())
    end

end

function make_xml(hmcxml::HMCOperatorXMLElement)
    name = hmcxml.geodesic ? bn.GEODESIC_HMC : bn.HMC
    el = new_element(name)
    set_attributes(el, [bn.WEIGHT => string(hmcxml.weight),
                        bn.NSTEPS => string(hmcxml.n_steps),
                        bn.STEP_SIZE => string(hmcxml.step_size),
                        bn.DRAW_VARIANCE => string(hmcxml.draw_variance),
                        bn.AUTO_OPTIMIZE => string(hmcxml.auto_optimize),
                        bn.GRAD_CHECK_COUNT => string(hmcxml.check_grad),
                        bn.GRAD_TOLERANCE => string(hmcxml.grad_tolerance)])
    set_attributes(el, preconditioningAttrs(hmcxml.preconditioning))

    jg_el = new_child(el, bn.JOINT_GRADIENT)

    for i = 1:length(hmcxml.grads)
        grad = hmcxml.grads[i]
        make_xml(grad)
        if hmcxml.already_made[i]
            add_ref_el(jg_el, grad.el)
        else
            add_child(jg_el, grad.el)
        end
    end

    add_ref_el(el, get_hmc_parameter(hmcxml.param_provider))
    if hmcxml.transform != ""
        t_el = new_child(el, bn.TRANSFORM)
        set_attribute(t_el, bn.TYPE, hmcxml.transform)
    end

    if length(hmcxml.mask) > 0
        m_el = new_child(el, bn.MASK)
        add_child(m_el, make_parameter(value=hmcxml.mask))
    end

    hmcxml.el = el
    return el
end

function get_parameter(op::HMCOperatorXMLElement)
    return get_hmc_parameter(op.param_provider)
end

################################################################################
## Factor model loadings gradients
################################################################################

## loadings prior
mutable struct LoadingsGradientXMLElement <: MyXMLElement
    el::XMLOrNothing
    ifxml::IntegratedFactorsXMLElement

    function LoadingsGradientXMLElement(ifxml::IntegratedFactorsXMLElement)

        return new(nothing, ifxml)
    end

end

function make_xml(gxml::LoadingsGradientXMLElement)
    ifxml = gxml.ifxml
    el = make_loadings_gradient(ifxml.loadings, ifxml.loadings_prior)

    gxml.el = el
    return el
end

## integrated likelihood
mutable struct FactorLoadingsGradientXMLElement <: MyXMLElement
    el::XMLOrNothing
    ifxml::IntegratedFactorsXMLElement
    tdlxml::TraitLikelihoodXMLElement
    name::String

    function FactorLoadingsGradientXMLElement(
            ifxml::IntegratedFactorsXMLElement,
            tdlxml::TraitLikelihoodXMLElement)

        return new(nothing, ifxml, tdlxml, bn.INTEGRATED_FACTOR_GRADIENT)
    end

end

function make_xml(flgxml::FactorLoadingsGradientXMLElement)
    make_xml(flgxml.ifxml)
    make_xml(flgxml.tdlxml)

    el = new_element(flgxml.name)
    add_ref_el(el, flgxml.ifxml.el)
    add_ref_el(el, flgxml.tdlxml.el)

    flgxml.el = el
    return el
end

mutable struct ScaledMatrixGradient <: MyXMLElement
    el::XMLOrNothing
    original_gradient::MyXMLElement
    component::String

    function ScaledMatrixGradient(
            original_gradient::MyXMLElement,
            component::String)
        return new(nothing, original_gradient, component)
    end
end

function make_xml(smg::ScaledMatrixGradient)
    el = new_element(bn.SCALED_MATRIX_GRADIENT)
    set_attribute(el, bn.COMPONENT, smg.component)
    sub_el = make_xml(smg.original_gradient)
    add_child(el, sub_el)
    smg.el = el
    return el
end

function NormalizedLoadingsGradientXMLElement(
    ifxml::IntegratedFactorsXMLElement,
    tdlxml::TraitLikelihoodXMLElement)
    grad = FactorLoadingsGradientXMLElement(ifxml, tdlxml)
    return ScaledMatrixGradient(grad, bn.MATRIX)
end

function ScaleLoadingsGradientXMLElement(
    ifxml::IntegratedFactorsXMLElement,
    tdlxml::TraitLikelihoodXMLElement)
    grad = FactorLoadingsGradientXMLElement(ifxml, tdlxml)
    return ScaledMatrixGradient(grad, bn.SCALE)
end

## sampled likelihood
mutable struct SampledLoadingsGradient <: MyXMLElement
    el::XMLOrNothing
    factor_model::LatentFactorModelXMLElement

    function SampledLoadingsGradient(
            factor_model::LatentFactorModelXMLElement)
        return new(nothing, factor_model)
    end
end

function make_xml(slg::SampledLoadingsGradient)
    el = new_element(bn.SAMPLED_LOADINGS_GRADIENT)
    fac_el = make_xml(slg.factor_model)
    add_ref_el(el, fac_el)
    slg.el = el
    return el
end

################################################################################
## preconditioning
################################################################################

mutable struct PreconditioningSchedule
    delay::Int
    update_frequency::Int
    max_update::Int
    memory::Int
end

function preconditioningAttrs(ps::PreconditioningSchedule)
    return Dict{String, String}(
        bn.PRECONDITIONING_UPDATE_FREQUENCY => string(ps.update_frequency),
        bn.PRECONDITIONING_MAX_UPDATE => string(ps.max_update),
        bn.PRECONDITIONING_DELAY => string(ps.delay),
        bn.PRECONDITIONING_MEMORY => string(ps.memory)
    )
end

function PreconditioningSchedule()
    return PreconditioningSchedule(100, 100, 0, 100)
end

struct NoPreconditioning <: Preconditioner
end

function preconditioningAttrs(::NoPreconditioning)
    return Dict{String, String}()
end

mutable struct AdaptiveDiagonalPreconditioning <: Preconditioner
    schedule::PreconditioningSchedule
end

function AdaptiveDiagonalPreconditioning()
    return AdaptiveDiagonalPreconditioning(PreconditioningSchedule())
end

function preconditioningAttrs(adp::AdaptiveDiagonalPreconditioning)
    attrs = preconditioningAttrs(adp.schedule)
    attrs[bn.PRECONDITIONING] = bn.ADAPTIVE_DIAGONAL
    return attrs
end

ext: jl


text:
"""
    abstract type AbstracTree end

Supertype for tree-like objects.
"""
abstract type AbstractTree end

"""
    children(t::AbstractTree)

Iterable collection of the node's children.
"""
function children(t::AbstractTree)
    abstractmethod(t)
end

"""
    parent(t::AbstractTree)

The node's parent. If `t` is a root, then `parent(t)==t`.
"""
function Base.parent(t::AbstractTree)
    abstractmethod(t)
end

"""
    isroot(t)

Return `true` if `t` is its own parent.
"""
function isroot(t)
    parent(t) == t
end

"""
    isleaf(t)

Return `true` if `t` has no children.
"""
function isleaf(t::AbstractTree)
    isempty(children(t))
end

"""
    depth(tree::AbstractTree,acc=0)

Recursive function to compute the depth of `node` in a a tree-like structure.

Overload this function if your structure has a more efficient way to compute
`depth` (e.g. if it stores it in a field).
"""
function depth(tree,acc=0)
    if isroot(tree)
        return acc
    else
        depth(parent(tree),acc+1)
    end
end

"""
    filter_tree(f,tree,isterminal=true)

Return a vector containing all the nodes of `tree` such that
`filter(node)==true`.  The argument `isterminal` can be used to control whether
to continue the search on `children` of nodes for which `f(node)==true`.
"""
function filter_tree(f,tree,isterminal=true)
    nodes = Vector{typeof(tree)}()
    filter_tree!(f,nodes,tree,isterminal)
end

"""
    filter_tree!(filter,nodes,tree,[isterminal=true])

Like [`filter_tree`](@ref), but appends results to `nodes`.
"""
function filter_tree!(f,nodes,tree,isterminal=true)
    if f(tree)
        push!(nodes,tree)
        # terminate the search along this path if terminal=true
        isterminal || map(x->filter_tree!(f,nodes,x,isterminal),children(tree))
    else
        # continue on on children
        map(x->filter_tree!(f,nodes,x,isterminal),children(tree))
    end
    return nodes
end

"""
    partition_by_depth(tree)

Given a `tree`, return a `partition` vector whose `i`-th entry stores all the nodes in
`tree` with `depth=i-1`. Empty nodes are not added to the partition.
"""
function partition_by_depth(tree)
    T         = eltype(tree)
    partition = Vector{Vector{T}}()
    depth = 0
    _partition_by_depth!(partition,tree,depth)
end

function _partition_by_depth!(partition,tree,depth)
    T = eltype(tree)
    if length(partition) < depth+1
         push!(partition,T[])
    end
    length(tree) > 0 && push!(partition[depth+1],tree)
    for chd in children(tree)
        _partition_by_depth!(partition,chd,depth+1)
    end
    return partition
end

"""
    partition_by_height(tree)

Given a `tree`, return a `partition` vector whose `i`-th entry stores all the nodes in
`tree` with `height=i-1`. The `height` of the tree is thus `lenth(partition)`,
with `partition(end)==tree`.
"""
function partition_by_height(tree)
    # TODO: how to do this more or less efficiently? One idea is to start at the
    # leaves, push them, then push their parents and recurse, making sure call
    # `unique!` as you go up in order to avoid duplicate parents.
end

# interface to AbstractTrees. No children is determined by an empty tuple for
# AbstractTrees.
AbstractTrees.children(t::AbstractTree) = isleaf(t) ? () : t.children
AbstractTrees.nodetype(t::AbstractTree) = typeof(t)

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
# # Register Basics
# ## Table of Contents
# * Construction and Storage
# * Basics Arithmatics
# * Fidelity and DensityMatrix
# * Batched Registers

using Yao
using LinearAlgebra

# ## Construction and Storage
# `AbstractRegister{B, T}` is abstract type that registers will subtype from. B is the batch size, T is the data type.
# Normally, we use a matrix as the `state` (with columns the batch and environment dimension) of a register, which is called `DefaultRegister{B, T}`.

# To initialize a quantum register, all you need is
# * `register(vec)`,
# * `zero_state(nbit)`,
# * `rand_state(nbit)`, both real and imaginary parts are random normal distributions,
# * `product_state(nbit, val=0)`, where val is an `Integer` as bitstring, e.g. `0b10011` or `19`,
# * `uniform_state(nbit)`, evenly distributed state, i.e. H|0>.

# e.g.

ψ1 = zero_state(5)
@show ψ1
@show nqubits(ψ1)
@show nactive(ψ1)   # number of activated qubits
@show nremain(ψ1)   # number of remaining qubits

ψ2 = ψ1 |> focus!(3,2,4)   # set activated qubits
@show ψ2
@show nqubits(ψ2)
@show nactive(ψ2)
@show nremain(ψ2)

@assert relax!(ψ2, (3,2,4)) == ψ1

# The total number of qubits here is 5, they are all acitve by default. `active` qubits are also called `system` qubits that are visible to operations, `remaining` qubits are the `environment`. `nremain == nqubits-nactive` always holds.
#
# **focus! & relax!**
# `focus!(reg, (3,2,4))` is equivalent to `reg |> focus!(3,2,4)`, which changes focused bits to `(3,2,4)`. Here from ψ1 -> ψ2, qubit line numbers change as
# `(active)(remaining): (1,2,3,4,5)() -> (3,2,4)(1,5)`
#
# `focus!` uses *relative positions*, which means it sees only active qubits and does not memorize original qubits positions. We take this convension to support **modulized design**. For example, if we want to insert a `QFT` blocks into some parent module, both the `QFT` and its parent do not need to know `original position`, which provides flexibility.
#
# `relax!` is the inverse process of `focus!`, `relax!(reg, (3,2,4))` will cancel the above operation. Here we have a second parameter since a register does not memorize original positions. This annoying feature can be circumvented using `focus!(reg, (3,2,4)) do ... end`, which will automatically restore your focus operation, see an example [here](@ref focusdo).
#
# Please also notice **APIs for changing lines order**
# * reorder!(reg, order), change lines order
# * reg |> invorder!, inverse lines order
#
# and
# * reg |> oneto(n), return a register `view`, with first `n` bits focused.

#---------------
# **Extending Registers**
# We can extend registers by either joining two registers or adding bits.

@assert product_state(3, 0b110) ⊗ product_state(3, 0b001) == product_state(6, 0b110001)

#--------------
reg = product_state(5, 0b11100)
@assert addbit!(copy(reg), 2) == product_state(7, 0b0011100) == zero_state(2) ⊗ reg

# **Storage**
# Let's dive into the storage of a register, there are three types `representation`s
# * `reg |> state`, matrix format, size = `(2^nactive, 2^nremain * nbatch)`
# * `reg |> rank3`, rank 3 tensor format, size = `(2^nactive, 2^nremain, nbatch)`
# * `reg |> hypercubic`, hypercubic format, size = `(2, 2, 2, ..., nbatch)`

# Here, we add a dimension `nbatch` to support parallism among registers.
# They are all different views of same memory. Please also check `statevec` and `relaxedvec` format, which prefer vectors whenever possible.

@show ψ1 |> state |> size
@show ψ1 |> rank3 |> size
@show ψ1 |> hypercubic |> size
@show ψ1 |> statevec |> size
@show ψ1 |> relaxedvec |> size;

# ### [Example](@id focusdo)
# multiply `|0>` by a random unitary operator on qubits `(3, 1, 5)` (relax the register afterwards).

using Yao.Intrinsics: rand_unitary

reg = zero_state(5)
focus!(reg, [3,1,5]) do r
    r.state = rand_unitary(8) * r.state
    r
end
@show reg.state;

# ## Basic Arithmatics
# `+, -, *, /, ⊗, '` are implemented.
#
# The adjoint of a register is also called `bra`, it can be used in calculating state overlap

ψ1 = rand_state(5)
ψ2 = rand_state(5)

# arithmatics
@show ψ1
@show ψ2
@show ψ3 = (0.3ψ1 + 2ψ2)/2 ⊗ ψ1
@assert ψ3 ≈ 0.15ψ1 ⊗ ψ1 + ψ2 ⊗ ψ1

# normalize ψ3
@assert ψ1 |> isnormalized && ψ2 |> isnormalized
@assert ψ3 |> isnormalized == false
@show ψ3 |> normalize! |> isnormalized

@show ψ3' * ψ3;

# ## Measure
# * `measure(reg; nshot=1)`, measure without collapsing state,
# * `measure!(reg)`, measure and collapse,
# * `measure_remove!(reg)`, measure focused bits and remove them,
# * `measure_reset!(reg, val=0)`, measure focused bits and reset them to some value,
# * `reset!(reg)`, collapse to specific value directly.
# * `select(reg, x)`, select subspace projected on specific basis, i.e. $|\phi\rangle = |x\rangle\langle x|\psi\rangle$.

#----------
# **measure**

@show product_state(5, 0b11001) |> measure  # please notice binary number `0b11001` is equivalent to `25`!
reg = rand_state(7)
@show measure(reg; nshot=5);          # measure multiple times

# **measure!**
reg = rand_state(7)
@show [measure!(reg) for i=1:5];  # measure! will collapse state

# **measure_reset!**
reg = rand_state(7)
@show [measure_reset!(reg, val=i*10) for i=1:5];   # measure_reset! will reset the measured bit to target state (default is `0`)

# **measure_remove!**
reg = rand_state(7)
@show measure_remove!(reg)
@show reg;

reg = rand_state(7)
@show measure_remove!(reg |> focus!(2,3))
@show reg;

# **select**
#
# select will allow you to get the disired measurement result, and collapse to that state.
# It is equivalent to calculating $|\phi\rangle = |x\rangle\langle x|\psi\rangle$.
reg = rand_state(9) |> focus!(1, 2, 3, 4)
@show ψ = select(reg, 0b1110)
@show ψ |> relax!;

## Fidelity and Density Matrix
ψ1 = rand_state(6)
ψ2 = rand_state(6)
@show fidelity(ψ1, ψ2)
@show tracedist(ψ1, ψ2)
@show ψ1 |> ρ
@show tracedist(ψ1 |> ρ, ψ2|> ρ);  # calculate trace distance using density matrix
@assert ψ1 |> probs ≈ dropdims(ψ1 |> ρ |> probs, dims=2)

# ## Batched Registers
#
# Most operations support batched register, which means running multiple registers in parallel.
ψ = rand_state(6, 3)
@show ψ
@show nbatch(ψ)
@show viewbatch(ψ, 2)  # this is a view of register at 2nd column of the batch dimension
@show repeat(ψ, 3);    # repeat registers in batch dimension

# **broadcasting along batch dimension**
@. ψ * 5 - 4 * ψ ≈ ψ
#------------------
X2 = put(5, 2=>X)       # X operator on 2nd bit, with total number of bit 5.
direct = copy(ψ) |> X2  # applying X2 directly
map(reg->reg |> X2, ψ)  # applying X2 using broadcasting, here X2 operator is applied inplace!
ψ .≈ direct


text:
struct Slot
  id::Symbol
end

Base.show(io::IO, s::Slot) = print(io, "@", s.id)

phislot(b, i) = Slot(Symbol(:phi_, b, :_, i))

# TODO: handle undef arguments properly.
function slots!(ir::IR)
  slots = Dict()
  for b in blocks(ir)
    # Block arguments
    if b.id != 1
      for (i, var) in enumerate(BasicBlock(b).args)
        slots[var] = phislot(b.id, i)
      end
      empty!(BasicBlock(b).args)
      empty!(BasicBlock(b).argtypes)
    end
    # Branches
    for br in BasicBlock(b).branches
      isreturn(br) && continue
      for (i, val) in enumerate(br.args)
        push!(b, :($(phislot(br.block, i)) = $val))
      end
      empty!(br.args)
    end
  end
  return varmap(x -> get(slots, x, x), ir)
end

using Core.Compiler: CodeInfo, SlotNumber

function slots!(ci::CodeInfo)
  ss = Dict{Slot,SlotNumber}()
  for i = 1:length(ci.code)
    ci.code[i] = MacroTools.prewalk(ci.code[i]) do x
      x isa Slot || return x
      haskey(ss, x) && return ss[x]
      push!(ci.slotnames, x.id)
      push!(ci.slotflags, 0x00)
      ss[x] = SlotNumber(length(ci.slotnames))
    end
  end
  return ci
end

function varargs!(meta, ir::IR, n = 0)
  isva = meta.method.isva
  argTs = ir.blocks[1].argtypes
  Ts = widenconst.(argTs[n+1:end])
  typed = !all(T -> T==Any, Ts)
  allTs = !isva ?
    Any[argTs[1:n]..., Tuple{Ts...}] :
    Any[argTs[1:n]..., typed ? Tuple{Ts[1:end-1]...,Ts[end].parameters...} : Any]
  args = copy(ir.blocks[1].args)
  deletearg!(ir, length(allTs)+1:length(args))
  empty!(argTs); append!(argTs, allTs)
  argmap = Dict{Variable,Any}()
  argis = 1:(length(Ts))
  newargs = map(reverse(argis)) do i
    argmap[args[i+n]] = pushfirst!(ir, Statement(nothing, type = Ts[i]))
  end |> reverse
  ir = varmap(x -> get(argmap, x, x), ir)
  for i in argis
    if isva && i == length(argis)
      i = length(Ts)
      xs, T = arguments(ir)[end], argTs[end]
      for _ = 1:i-1
        T = typed ? Tuple{T.parameters[2:end]...} : Any
        st = Statement(xcall(Base, :tail, xs), type = T)
        xs = insertafter!(ir, xs, st)
      end
      ir[newargs[i]] = xcall(:identity, xs)
    else
      ir[newargs[i]] = xcall(Base, :getfield, arguments(ir)[end], i)
    end
  end
  return ir
end

# TODO this is hacky and leaves `ir.defs` incorrect
function splicearg!(meta, ir::IR, name)
  args = arguments(ir)
  push!(ir.defs, (1, -1))
  arg = var(length(ir.defs))
  pushfirst!(args, arg)
  pushfirst!(ir.blocks[1].argtypes, Any)
  pushfirst!(meta.code.slotnames, name)
  return arg
end

function update!(ci::CodeInfo, ir::Core.Compiler.IRCode)
  Core.Compiler.replace_code_newstyle!(ci, ir, length(ir.argtypes)-1)
  ci.inferred = false
  ci.ssavaluetypes = length(ci.code)
  slots!(ci)
  fill!(ci.slotflags, 0)
  return ci
end

update!(meta, ir::IR) = update!(meta, Core.Compiler.IRCode(slots!(ir)))


text:
#This script will set up and solve the equations, in fourier space, for the gravitational wave
#perturbation h as a function of time. We solve (Dtt+k^2)h(t,k)=T(t,k).

#TODO: The integrator does not evolve in time, always spits dh = 0, fix it.

using FFTW
using LinearAlgebra
using Interpolations
using OrdinaryDiffEq
#Transform of only real functions. Gives back a matrix of 2 component vector for
# the momenta k and the matrix for the mode coefficients fk.
abstract type Parameters end

struct param{T<:Interpolations.GriddedInterpolation, TP<:Real} <: Parameters# where {T<:Interpolations.GriddedInterpolation, TP<:Real}
    f  :: T
    x  :: Array{TP,1}
    y  :: Array{TP,1}
    kx :: Array{TP,2}
    ky :: Array{TP,2}
    dt :: TP
end

function get_evol_variable(u::Array{T,2}, u_t::Array{T,2}) where {T<:Complex}
    N      = length(u)
    u_evol = im.*zeros(N)
    u_evol = reshape(u, N)
    append!(u_evol, reshape(u_t, N))

    u_evol
end

function get_tensors(u_evol::Array{T,1}, Nkx::N, Nky::N) where {T<:Complex, N<:Int}
    Ntot     = length(u_evol)
    idxsplit = Int(Ntot/2)
    h        = reshape(u_evol[1:idxsplit],Nkx,Nky)
    h_t      = reshape(u_evol[idxsplit+1:end],Nkx,Nky)

    h, h_t
end

function Fourier_Transform_2D(f::Array{T,2}) where {T<:Real}
    Nx   = length(f[:,1])
    Ny   = length(f[1,:])
    plan = plan_rfft(f)

    fk   = 1/(Nx*Ny).*(plan * f)
end

function initial_conditions(ff::VEVTimeSeries, kx::Array{T,2}, ky::Array{T,2}) where {T<:Real}
    Tk  = Fourier_Transform_2D(ff[1,:,:])
    Nkx = length(kx[:,1])
    Nky = length(ky[1,:])

    tol = 10^-10
    indeces = findall(abs.(Tk).<tol)
    for i in indeces
        Tk[i] = 0.0
    end

    h_t    = im.*zeros(Nkx,Nky)
    h      = Tk./(kx.^2 + ky.^2)
    h[1,1] = Tk[1,1]

    get_evol_variable(h,h_t)
end

function rhs(h_evol::Array{T,1}, param::Parameters, t::TP) where {T<:Complex, TP<:Real}
    f_inter = param.f(t,param.x,param.y)
    kx      = param.kx
    ky      = param.ky
    Nkx     = length(kx[:,1])
    Nky     = length(ky[1,:])
    Tk      = Fourier_Transform_2D(f_inter)
    tol     = 10^-10
    indeces = findall(abs.(Tk).<tol)
    for i in indeces
        Tk[i] = 0.0
    end
    h, h_t    = get_tensors(h_evol,Nkx,Nky)
    dh        = h_t
    dh_t      = Tk - (kx.^2 + ky.^2) .*h
    dh_t[1,1] = 0.0

    get_evol_variable(dh,dh_t)
end

#Computes the component associated to field in the bounbdary stress tensor.
#E.g. field=energy gives us h00 as solution.
function solve_GW(dirname::String, field::Symbol; dt::T, alg, option::String="static") where {T<:Real}
    ff       = VEVTimeSeries(dirname, field)
    tt, x, y = get_coords(ff,:,:,:)
    tspan    = (tt[1],tt[end])
    dx       = x[2]-x[1]
    dy       = y[2]-y[1]
    kxx      = 2π.*rfftfreq(length(x), 1/dx)
    kyy      = 2π.*fftfreq(length(y), 1/dy)
    Nkx      = length(kxx)
    Nky      = length(kyy)
    kx       = zeros(Nkx,Nky)
    ky       = zeros(Nkx,Nky)
    for j in eachindex(kyy)
        for i in eachindex(kxx)
            kx[i,j] = kxx[i]
            ky[i,j] = kyy[j]
        end
    end
    h0_evol  = initial_conditions(ff,kx,ky)
    ff_inter = interpolate((tt,x,y,), ff[:,:,:], Gridded(Linear()))
    prm      = param{typeof(ff_inter),typeof(x[1])}(ff_inter,x,y,kx,ky,dt)

    problem    = ODEProblem(rhs, h0_evol, tspan, prm)
    integrator = init(problem, alg, save_everystep=false, dt=dt, adaptive=false)

    h   = im.*zeros(Int(floor((tspan[2]-tspan[1])/dt+1)),Nkx,Nky)
    h_t = similar(h)
    Tk  = similar(h)
    n = 0
    for (u, t) in tuples(integrator)
        n += 1
        h[n,:,:], h_t[n,:,:] = get_tensors(u,Nkx,Nky)
        Tk[n,:,:]            = Fourier_Transform_2D(ff_inter(t,x,y))
        umax  = maximum(abs.(u[2:Nkx*Nky]))
        dumax = maximum(abs.(u[Nkx*Nky+2:end]))

        println("t = $t")
        println("\n")
        println("max |h| = $umax")
        println("max |h_t| = $dumax")
        println("------------------------------------------------------------")
        println("\n")
    end
    h, h_t, Tk, kxx, kyy
end


text:
"""
generates Z (full)  with M (natural mortality) and F (catch mortality)
"""
function genZ(M,F)
    (tmax,jmax)=size(M)
    Z=zeros(tmax,jmax)
    Z=M+F
    return(Z)
end

"""
generates G with M (natural mortality) and F (catch mortality), Z (full mortality)
"""
function genG(M,F)
    (tmax,jmax)=size(M)
    Z=genZ(M,F)
    G=zeros(tmax,jmax)
    for t in 1:tmax
        for j in 1:jmax
            G[t,j]=F[t,j]/Z[t,j]*(1-exp(-Z[t,j]))
        end
    end
    return(G)
end

#FILTERING

"""
filtering: returns a vector of deviations (sigma^2^f) for a cohort
sigma^2_0 - variance  for initial age
sigma^2_e -  variance  for the evidance vector (catches)
sigma^2_x -  variance  for the state vector (population)
l - the length of the cohort
"""
function sigma2f(sigma20,sigma2x,sigma2e,l)
    s2f=zeros(l)
    s2f[1]=sigma20*sigma2e/(sigma20+sigma2e)
    if l >1
        for i in 1:l-1
        s2f[i+1]=(s2f[i]+sigma2x)*sigma2e/(s2f[i]+sigma2x+sigma2e)
        end
    end
    return(s2f)
end

"""
filtering: returns a vector of means (mu^f) for a cohort
alpha - coefficients for the state x_{t+1}=x_t +alpha_t
beta - coefficients for the observation y_{t}=x_t +beta_t
 e - the evidance vector
mu_0 - mean for the initial age
sigma^2_0 - variance (deviation is the square root of its variance) for initial age
mu_e - mean for the evidence
sigma^2_e - variance  for the evidance  (catch)
mu_x - mean for the state
sigma^2_e - variance for the state (population)
l - the length of the cohort
"""
function muf(alpha, beta, e,mu0,sigma20,mux,sigma2x,mue,sigma2e)
    l=length(e)
    muf=zeros(l)
    muf[1]=mu0
    #muf[1]=x0+beta[1]+mue
    s2f= sigma2f(sigma20,sigma2e,sigma2x,l)
    if l>1
        for i in 1:l-1
            muf[i+1]=(muf[i]+alpha[i]+mux)+(e[i+1]-(muf[i]+alpha[i]+mux+beta[i+1]+mue))*
        (s2f[i]+sigma2x)/(s2f[i]+sigma2x+sigma2e)
        end
    end
        return(muf)
    end

"""
filtering: returns filtered cohorts represented as vecs: csvlogts1muf - means,csvlogts1s2f - sigmas^2 (deviations),
"""
function genmus2f(ts1,c1,yb,M,F,mu0,s20,mux,s2x,muy,s2y)
    tmax,jmax=size(ts1)

    csvlogts1= csmtx2csvs(csmtx(yb,tmax,jmax,map(x->log(x),ts1)))
    csvlogc1= csmtx2csvs(csmtx(yb,tmax,jmax,map(x->log(x),c1)))
    csvZ= csmtx2csvs(csmtx(yb,tmax,jmax,genZ(M,F)))
    csvlogG= csmtx2csvs(csmtx(yb,tmax,jmax,map(x->log(x),genG(M,F))))

    csvlogts1muf= deepcopy(csvlogts1)
    csvlogts1s2f= deepcopy(csvlogts1)

    for i in 1:tmax+jmax-1
        x=csvlogts1.vecs[i].vec
        e=csvlogc1.vecs[i].vec
        alpha=-csvZ.vecs[i].vec
        beta=csvlogG.vecs[i].vec

        vec=muf(alpha, beta, e,x[1],s20,mux,s2x,muy,s2y)
        vecs2=sigma2f(s20,s2x,s2y,length(e))

        csvlogts1muf.vecs[i].vec=vec
          csvlogts1s2f.vecs[i].vec=vecs2

       # if i == n
        #    println(vec)
       #     println(csvlogts1f.vecs[i].vec)
      #  println(x)

       # end

    end
    return(csvlogts1muf,csvlogts1s2f)

end

#PREDICTION

"""
prediction: returns the vector of deviation predictions (sigma^2^p_{t+1},..., sigma^2^p_{t+k+1}) for a cohort
sigma^2_x - variance for the state vector (population)
l - the length of the cohort
"""
function sigma2p(sigma2ft,sigma2x,k)
    sigma2pk=zeros(k+1)
    sigma2pk[1]=sigma2ft+sigma2x
    for t in 1:k
        sigma2pk[t+1]=sigma2pk[t]+sigma2x
    end
    return(sigma2pk)
end

"""
prediction: returns the vector of mean predictions (mu^p_{t+1}, ..., mu^p_{t+k+1}) for a cohort
alphak=(alpha^t, ..., alpha_{t+k}) - (virtual) coefficients for x_{t+k+1}=x_{t+k} +alpha_{t+k} -
muft - mu^f_t
mu_x - mean for the state
k >= 1 - the prediction step
"""
function mup(alphatk,mux,muft,k)
    mupk=zeros(k+1)
    mupk[1]=muft+alphatk[1]+mux
            for t in 1:k
            mupk[t+1]=mupk[t]+alphatk[t+1]+mux
        end
        return(mupk)
    end

"""
makes k predictions for all cohorts (a^e+1, .., a^e+k)  with yr>=yr2. cuts extra elements with age greater than jmax.
returns filtered cohorts represented as vecs: csvlogts1mup - means,csvlogts1s2p - sigmas^2 (deviations),
mup and s2p (means and sigmas^2 (deviations))
"""
function genmus2p(ts1,c1,yb,M,F,mu0,s20,mux,s2x,muy,s2y,k)
    tmax,jmax=size(ts1)

    csvlogts1= csmtx2csvs(csmtx(yb,tmax,jmax,map(x->log(x),ts1)))
    csvlogc1= csmtx2csvs(csmtx(yb,tmax,jmax,map(x->log(x),c1)))
    csvZ= csmtx2csvs(csmtx(yb,tmax,jmax,genZ(M,F)))
    csvlogG= csmtx2csvs(csmtx(yb,tmax,jmax,map(x->log(x),genG(M,F))))

    csvlogts1muf= deepcopy(csvlogts1)
    csvlogts1s2f= deepcopy(csvlogts1)
    csvlogts1mufp= deepcopy(csvlogts1)
     csvlogts1mufm= deepcopy(csvlogts1)

    csvlogts1mup= deepcopy(csvlogts1)
    csvlogts1s2p= deepcopy(csvlogts1)

    yr1=yb
    yr0=yr1-jmax+1
    yr2=yr1+tmax-jmax
    yr3=yr1+tmax-1

    for i in 1:tmax+jmax-1
        x=csvlogts1.vecs[i].vec
        e=csvlogc1.vecs[i].vec
        alpha=-csvZ.vecs[i].vec
        beta=csvlogG.vecs[i].vec
        yr=csvlogts1.vecs[i].yb

        vec=muf(alpha, beta, e,x[1],s20,mux,s2x,muy,s2y)
        vecs2=sigma2f(s20,s2x,s2y,length(e))

        csvlogts1muf.vecs[i].vec=vec
        csvlogts1s2f.vecs[i].vec=vecs2

       # if i == n
        #    println(vec)
       #     println(csvlogts1f.vecs[i].vec)
      #  println(x)
        csvlogts1mup.vecs[i].vec=vec

        if k > 0
        muzs=zeros(k)

        s2zs=zeros(k)

       # end
     if yr> yr2
        begin
            muzs[1],s2zs[1]=last(vec)+last(alpha)+mux,last(vecs2)+s2x

            for l in 2:k

                      muzs[l],s2zs[l]=
                        if length(vec)+l<=jmax
                        muzs[l-1]+last(alpha)+mux,s2zs[l-1]+s2x
                        else
                           muzs[l],s2zs[l]= muzs[jmax-length(vec)],s2zs[jmax-length(vec)]
                        end
            end
                       #if length(vec)+k+1 <= jmax
                        #    for l in length(vec)+k+1:jmax
                         # muzs[l],s2zs[l]=muzs[length(vec)+k],s2zs[length(vec)+k]
                       # end
                       # end

                    #println(i)
                    #println(yr)
         #println(muzs)
         #  println(s2zs)

      append!(vec,muzs)
      append!(vecs2,s2zs)
                    if length(vec)>jmax
                        vec, vecs2= vec[1:jmax], vecs2[1:jmax]
                    end
        end
        end
            csvlogts1mup.vecs[i].vec=vec
            csvlogts1s2p.vecs[i].vec=vecs2

    end
    end

    return(csvlogts1mup,csvlogts1s2p)

end

#SMOOTHING

"""
backward: returns a vector of deviations (sigma^2^b) for a cohort
sigma^2_0 - variance  for initial age
sigma^2_e -  variance  for the evidance vector (catches)
sigma^2_x -  variance  for the state vector (population)
l - the length of the cohort
"""
function sigma2b(sigma20,sigma2x,sigma2e,l)
    s2b=zeros(l)
    s2f=sigma2f(sigma20,sigma2x,sigma2e,l)
    s2b[l]=sigma2x+sigma2e
    if l >1
        i=l
        while i>1
        s2b[i-1]=(sigma2x*sigma2e +s2b[i]*(sigma2x+sigma2e))/(sigma2e+s2b[i])
            i=i-1
        end
    end
    return(s2b)
end

"""
smoothing: makes several cohorts represented as vecs: csvlogts1mus - means,csvlogts1s2s - sigmas^2 (deviations)
"""
function genmus2s(ts1,c1,yb,M,F,mu0,s20,mux,s2x,muy,s2y)
    tmax,jmax=size(ts1)

    csvlogts1= csmtx2csvs(csmtx(yb,tmax,jmax,map(x->log(x),ts1)))
    csvlogc1= csmtx2csvs(csmtx(yb,tmax,jmax,map(x->log(x),c1)))
    csvZ= csmtx2csvs(csmtx(yb,tmax,jmax,genZ(M,F)))
    csvlogG= csmtx2csvs(csmtx(yb,tmax,jmax,map(x->log(x),genG(M,F))))

    csvlogts1muf= deepcopy(csvlogts1)
    csvlogts1s2f= deepcopy(csvlogts1)

    csvlogts1mub= deepcopy(csvlogts1)
    csvlogts1s2b= deepcopy(csvlogts1)

    csvlogts1mus= deepcopy(csvlogts1)
    csvlogts1s2s= deepcopy(csvlogts1)

    for i in 1:tmax+jmax-1
        x=csvlogts1.vecs[i].vec
        e=csvlogc1.vecs[i].vec
        alpha=-csvZ.vecs[i].vec
        beta=csvlogG.vecs[i].vec

        li=length(e)

        #forward

        vecmuf=muf(alpha, beta, e,x[1],s20,mux,s2x,muy,s2y)
        vecs2f=sigma2f(s20,s2x,s2y,li)

        csvlogts1muf.vecs[i].vec=vecmuf
          csvlogts1s2f.vecs[i].vec=vecs2f

        #backward

        vecmub=zeros(li)
       vecs2b=sigma2b(s20,s2x,s2y,li)

        vecmub[li]=vecmuf[li]

       if li >1
            vecmub[li]=e[li]-alpha[li-1]-mux-beta[li]-muy
        i1=li

        while i1>2
        vecmub[i1-1]=((vecmub[i1]-alpha[i1-2]-mux)*s2y+(e[i1-1]-alpha[i1-2]-beta[i1-1]-mux)*vecs2b[i1])/
                (s2y+vecs2b[i1])
                i1=i1-1

        end
    end

        csvlogts1mub.vecs[i].vec=vecmub
          csvlogts1s2b.vecs[i].vec=vecs2b

     #smoothing

         vecmus=zeros(li)
        vecs2s= zeros(li)

        vecmus[li]=vecmuf[li]
        vecs2s[li]=vecs2f[li]

        if li >1

        for i2 in 1:li-1

            vecmus[i2]=(vecmuf[i2]*vecs2b[i2+1]+vecmub[i2+1]*vecs2f[i2])/(vecs2b[i2+1]+vecs2f[i2])
            vecs2s[i2]=(vecs2b[i2+1]*vecs2f[i2])/(vecs2b[i2+1]+vecs2f[i2])
                #println("i2= ",i2)
                #println("vecs2s=", vecs2s[i2],"  vecs2f=", vecs2f[i2]," vecs2b=", vecs2b[i2])
        end
            end
         csvlogts1mus.vecs[i].vec=vecmus
          csvlogts1s2s.vecs[i].vec=vecs2s

    end

    return(csvlogts1mus,csvlogts1s2s)

end

"""
at +p%
"""
function genQuantile(lmu,lsigma2,p)
    l=length(lmu)
    lq=zeros(l)
    for i in 1:l
        lq[i]=if lsigma2[i]>0.
           quantile(Distributions.Normal(lmu[i],sqrt(lsigma2[i])),p)
            else
            lmu[i]
        end

        end
    return(lq)
end

"""
at -p%
"""
function genQuantile1(lmu,lsigma2,p)
    l=length(lmu)
    lq=zeros(l)
    for i in 1:l
        lq[i]=if lsigma2[i]>0.
           2*lmu[i]-quantile(Distributions.Normal(lmu[i],sqrt(lsigma2[i])),p)
            else
            lmu[i]
        end

        end
    return(lq)
end

max_stars_repo_name: picestcode/Fishmetica

max_issues_repo_name: picestcode/Fishmetica

max_forks_repo_name: picestcode/Fishmetica


text:
module TextUnidecode

const table = Vector{Vector{String}}(undef, 256)

"""
Convert non-ascii characters to "good enough" ascii.
```jldoctest
julia> unidecode("南无阿弥陀佛")
"Nan Wu A Mi Tuo Fo"

julia> unidecode("あみだにょらい")
amidaniyorai
```
"""
function unidecode(str::String)::String
    new_string = Vector{String}()
    @inbounds for c in str
        code_point = codepoint(c)
        if code_point < 0x80
            push!(new_string, string(c))
            continue
        elseif code_point > 0xffff
            continue
        else
            # Chop off last two hex digits
            section = code_point >> 8
            # Last two hex digits
            pos = code_point % 256;
            cache = get_cache(section)
            push!(new_string, cache[pos + 1])
        end
    end
    rstrip(join(new_string))
end

function get_cache(section::UInt32)::Vector{String}
    if isassigned(table, section+1)
        table[section + 1]
    else
        path = joinpath(@__DIR__, "..", "resources", "X$(string(section, base = 16, pad = 3))" )
        table[section + 1] = readlines(path)
    end
end

export unidecode
end


text:
@deprecate sexp(::Type{IntSxp}, x) sexp(RClass{:integer}, x)
@deprecate sexp(::Type{RealSxp}, x) sexp(RClass{:numeric}, x)
@deprecate sexp(::Type{StrSxp}, x) sexp(RClass{:character}, x)
@deprecate sexp(::Type{RawSxp}, x) sexp(RClass{:raw}, x)
@deprecate sexp(::Type{VecSxp}, x) sexp(RClass{:list}, x)


text:
struct QuerySynthesisFunctionOracle <: Oracle
    f::Function
end

function ask_oracle(oracle::QuerySynthesisFunctionOracle, query_objects::Array{T, 2})::Vector{Symbol} where T <: Real
    return oracle.f(query_objects)
end


text:
using SamplingHS3

gm2_mix = []
for zi=1:10;
    found_es = false
    vd_pre   = NaN * ones(2,2)
    while(!found_es)
        vd = qr(rand(2,2))[1]
        vd_pre = vd*diagm(rand(2))*vd'
        found_es = issymmetric(vd_pre)
    end
    gpi = SamplingHS3.GaussParam( 8*rand(2) , vd_pre )
    gm2_mix = [ gm2_mix ;  gpi ]
end

gm2 = SamplingHS3.GaussMixture(gm2_mix,0.5+rand(10))
xx  = SamplingHS3.draw_gaussmixture( gm2 , 100000 )

using Plots
#SamplingHS3.gr_hist_2d( xx , nbins=40 )

bb_G = [eye(2);-eye(2)]
bb_h = [10.;10;10;10]

X_init = SamplingHS3.sample_unif_hr( bb_G , bb_h , zeros(2) , 1000 , 2000 )

# now sample dat with nested sampling
f_cost = (x) -> -log(SamplingHS3.pdf_gaussmixture_nonnormalized(gm2,x))

x_unif = SamplingHS3.unifsamp_box([-10. 10;-10 10],2000)
x_unif_cost = [  f_cost(x_unif[:,zi]) for zi in (1:2000) ]

x_unif_cost = [  SamplingHS3.mvnpdf([0.;0],1.0*eye(2),x_unif[:,zi:zi]) for zi in (1:2000) ]

using Plots
gr()
gui()
Plots.scatter( x_unif[1,:] , x_unif[2,:] ; zcolor=x_unif_cost[:],c=:lightrainbow,ms=4.0 )

# sample..
#ns_result = SamplingHS3.nested_sampling( X_init , f_cost , bb_G , bb_h , SamplingHS3.NSConfig(8,0.8,SamplingHS3.StopCriterion(max_time=3.),SamplingHS3.MCMCConfig()) )
ns_result = SamplingHS3.nested_sampling( X_init[:,:] , f_cost , bb_G , bb_h , SamplingHS3.NSConfig(8,0.95,SamplingHS3.StopCriterion(max_time=10.0),SamplingHS3.MCMCConfig()) )

x_live_samples_1 = SamplingHS3.reconstruct_live_samples(ns_result.ns_iter_data)
x_live  = x_live_samples_1[1]
xf_live = x_live_samples_1[2]
# plot live samples..
# P_subplots = plot(layout=(4,4) , size=(2000,2000))
for zi in 1:length(x_live)
    display( Plots.scatter(x_live[zi][1,:] , x_live[zi][2,:]  , ma = 0.75 , zcolor=(1:size(x_live[zi],2)) , m=2, msw=0 , c=:lightrainbow) )
    #scatter!( x_live[zi][1,:] , x_live[zi][2,:] , subplot = zi , ma = 0.75 , zcolor=(1:size(x_live[zi],2)) , m=2, msw=0 , c=:lightrainbow , lab="grad" )
end
gui()
gui(plot(P_subplots , size=(4000,4000)))

using Plots
gr()
for n in 1:100
  display(histogram(randn(10000)))
  sleep(0.05)
end

ns_x_sequence = reduce( (x,y) -> [x y] , map( x -> x.X_sorted_out , ns_result.ns_iter_data ) )
scatter( ns_x_sequence[1,:] , ns_x_sequence[2,:] , ma = 0.75 , zcolor=(1:size(ns_x_sequence,2)) , m=2, msw=0 , c=:lightrainbow , lab="grad" )

ssc_data = SamplingHS3.sim_shell_compression(ns_result.ns_iter_data)
ei_result = SamplingHS3.eval_integral(400.,ns_result.ns_iter_data)

rs_result = SamplingHS3.resample_ns_data( ns_result.ns_iter_data )

scatter!(rs_result[1][1,:],rs_result[1][2,:])

scatter(ns_result.ns_iter_data[16].X_sorted_out[1,:],ns_result.ns_iter_data[16].X_sorted_out[2,:])

scatter!( ns_result.X_sequence[1,1:500] , ns_result.X_sequence[2,1:500] )
scatter( ns_result.X_sequence[1,1:500] , ns_result.X_sequence[2,1:500] )
scatter!( ns_result.X_sequence[1,1000:1500] , ns_result.X_sequence[2,1000:1500] )
scatter( ns_result.X_sequence[1,1300:1500] , ns_result.X_sequence[2,1300:1500] )
plot!()

# create plot of some deocrrelation trajectories:
heatmap()

## 3 d test:

gm3_mix = []
for zi=1:10;
    found_es = false
    vd_pre   = NaN * ones(3,3)
    while(!found_es)
        vd = qr(rand(3,3))[1]
        vd_pre = vd*diagm(rand(3))*vd'
        found_es = issymmetric(vd_pre)
    end
    gpi = SamplingHS3.GaussParam( 8*rand(3) , vd_pre )
    gm3_mix = [ gm3_mix ;  gpi ]
end

gm3 = SamplingHS3.GaussMixture(gm3_mix,0.5+rand(10))
xx3  = SamplingHS3.draw_gaussmixture( gm3 , 100000 )
bb3_G = [eye(3);-eye(3)]
bb3_h = [10.;10;10;10;10;10]

f_cost_3d = (x) -> -log(SamplingHS3.pdf_gaussmixture_nonnormalized(gm3,x))
scatter(xx3[1,1:2000] , xx3[2 ,1:2000] , xx3[3,1:2000] , ma = 0.1 , m=2, msw=0 , c=:lightrainbow , lab="grad")

# sample..
X_init = SamplingHS3.sample_unif_hr( bb3_G , bb3_h , zeros(3) , 400 , 2000 )
ns_result = SamplingHS3.nested_sampling( X_init , f_cost_3d , bb3_G , bb3_h , SamplingHS3.NSConfig(8,0.8,SamplingHS3.StopCriterion(max_time=20.)) )

ns_x_sequence = reduce( (x,y) -> [x y] , map( x -> x.X_sorted_out , ns_result.ns_iter_data ) )
scatter( ns_x_sequence[1,:] , ns_x_sequence[2,:] , ns_x_sequence[3,:] , ma = 0.75 , zcolor=(1:size(ns_x_sequence,2)) , m=2, msw=0 , c=:lightrainbow , lab="grad" )
scatter( ns_x_sequence[1,end-3000:end] , ns_x_sequence[2,end-3000:end] , ns_x_sequence[3,end-3000:end] , ma = 0.75 , zcolor=(1:2001) , m=2, msw=0 , c=:lightrainbow , lab="grad" )

ssc_data = SamplingHS3.sim_shell_compression(ns_result.ns_iter_data)
ei_result = SamplingHS3.eval_integral(400.,ns_result.ns_iter_data)
rs_result = SamplingHS3.resample_ns_data( ns_result.ns_iter_data )

scatter!(rs_result[1][1,:],rs_result[1][2,:],rs_result[1][3,:])

#binned_a = SamplingHS3.bin2d( ns_x_sequence[[1,2],:] , ones(size(ns_x_sequence,2)) , 40,40 , x -> sum(x) )
binned_a = SamplingHS3.bin2d( ns_x_sequence[[1,2],1:3000] , ones(3000) , 40,40 , x -> sum(x) )
heatmap(binned_a)
binned_b = SamplingHS3.bin2d( ns_x_sequence[[1,3],1:3000] , ones(3000) , 40,40 , x -> sum(x) )
heatmap(binned_b)


text:
"""
Row-by-row sparsification using a quantile threshold
"""
function sparsify_signal!(U_guess; quantile_threshold=0.95)
    sparsify(x, eps) = abs(x) > eps ? x : 0;
    for i in 1:size(U_guess, 1)
        tmp = abs.(U_guess[i, :])
        epsilon = quantile(tmp[tmp.>0.0], quantile_threshold)
        for i2 in 1:size(U_guess,2)
            U_guess[i,i2] = sparsify(U_guess[i,i2], epsilon)
        end
    end
end

"""
Row-by-row sparsification using a quantile threshold
"""
function sparsify_signal(U_guess; quantile_threshold=0.95)
    for i in 1:size(U_guess, 1)
        tmp = abs.(U_guess[i, :])
        epsilon = quantile(tmp[tmp.>0.0], quantile_threshold)
        ind = abs.(U_guess[i, :]) .< epsilon
        U_guess[i, ind] .= 0.0
    end
    return U_guess
end

"""Keeps only largest 'n' terms"""
function keep_n_terms(X, n)
    ind = sortperm(vec(abs.(X)), rev=true)
    new_X = zeros(size(X))
    new_X[ind[1:n]] = X[ind[1:n]]
    return new_X
end

"""
Row-by-row sparsification using a hard threshold
"""
function threshold_signal!(U_guess; hard_threshold=1e-5)
    sparsify(x, eps) = abs(x) > eps ? x : 0;
    for i in 1:size(U_guess, 1)
        tmp = abs.(U_guess[i, :])
        for i2 in 1:size(U_guess,2)
            U_guess[i,i2] = sparsify(U_guess[i,i2], hard_threshold)
        end
    end
end

# """
# Simple implementation of sequential least squares thresholding for a linear
#     regression model
# """
# function sparse_regression(X::Matrix, y::Matrix;
#                            num_iter=2,
#                            quantile_threshold=0.1,
#                            hard_threshold=nothing,
#                            num_terms=nothing)
#     if size(y,1) > 1
#         A = zeros(size(y,1), size(X,1))
#         # For the 2d predictor case, just loop over rows
#         for i in 1:size(y,1)
#             if num_terms !== nothing
#                 # Constrain the number of terms row by row
#                 if length(num_terms)>1
#                     n = num_terms[i]
#                 else
#                     n = num_terms
#                 end
#             else
#                 n = nothing
#             end
#             A[i, :] = sparse_regression(X, y[i:i,:],
#                             num_iter=num_iter,
#                             quantile_threshold=quantile_threshold,
#                             hard_threshold=hard_threshold,
#                             num_terms=n)
#         end
#         return A
#     end
#     # If this is only one row, initialize with L2
#     A = y/X
#     # sparsify_signal!(A, quantile_threshold=quantile_threshold)
#     # Redo, removing the rows in X that were zeroed out
#     for i in 1:num_iter
#         if hard_threshold == nothing && num_terms == nothing
#             sparsify_signal!(A, quantile_threshold=quantile_threshold)
#         elseif num_terms == nothing
#             threshold_signal!(A, hard_threshold=hard_threshold)
#         else
#             A = keep_n_terms(A, num_terms)
#         end
#         which_rows = vec(abs.(A) .> 0)
#         A_subset = y/(X[which_rows, :])
#         A[which_rows] = A_subset
#         A[.!which_rows] .= 0.0 # Retain the previous zeros
#     end
#
#     return A
# end

"""
Removes outliers, by default 3 times the median
"""
function remove_outliers(dat::Vector, metric=median, max_val=3)
    m = 3*median(abs.(dat))
    return dat[abs.(dat) .< m]
end

"""
Replaces outliers (by default 3 times the median) with 'replace_val'

function replace_outliers(dat::Vector,
            replace_val=NaN, metric=median, max_val=3)
"""
function replace_outliers(dat::Vector,
            replace_val=NaN, metric=median, max_val=3)
    new_dat = copy(dat)
    m = 3*median(abs.(dat))
    for i in 1:length(dat)
        if dat[i] > m
            new_dat[i] = replace_val
        end
    end
    return new_dat
end

# export sparse_regression, sparsify_signal, threshold_signal!,
#         keep_n_terms
export sparsify_signal, threshold_signal!,
        keep_n_terms, remove_outliers, replace_outliers


text:
# modified by N Laohakunakorn
# University of Edinburgh 2020

include("Kinetics.jl");
include("Control.jl");

# ----------------------------------------------------------------------------------- #
# Copyright (c) 2017 Varnerlab
# School of Chemical Engineering Purdue University
# W. Lafayette IN 46907 USA

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #
#function MassBalances(dxdt_vector,x,t,data_dictionary)
function MassBalances(dxdt_vector,x,data_dictionary,t)
# ---------------------------------------------------------------------- #
# MassBalances.jl was generated using the Kwatee code generation system.
# Username: nicholas
# Type: CFPS-JULIA
# Version: 1.0
# Generation timestamp: 03-24-2017 13:01:36
#
# Arguments:
# t  - current time
# x  - state vector
# dxdt_vector - right hand side vector
# data_dictionary  - Data dictionary instance (holds model parameters)
# ---------------------------------------------------------------------- #

# Set RIBOSOME_START_CAT to RIBOSOME_START_CAT initial condition + RIBOSOME initial condition - RIBOSOME concentration
initial_condition = data_dictionary["INITIAL_CONDITION_ARRAY"]
x[76] = min(x[76], initial_condition[76]+initial_condition[77])
x[77] = min(x[77], initial_condition[76]+initial_condition[77])

# Correct negative x's = throws errors in control even if small -
idx = x.<0;
x[idx] = zeros(size(x[idx]));

# Call the kinetics function -
rate_vector = Kinetics(t,x,data_dictionary);

# Call the control function -
rate_vector = Control(t,x,rate_vector,data_dictionary);

# Encode the balance equations as a matrix vector product -
S = data_dictionary["STOICHIOMETRIC_MATRIX"];
tmp_vector = S*rate_vector;
number_of_states = length(tmp_vector);

for state_index in collect(1:number_of_states)
    dxdt_vector[state_index] = tmp_vector[state_index];
end

dxdt_vector[[138;143]] .= 0; # Set M_h_c and M_o2_c to steady-state

end

max_stars_repo_name: Laohakunakorn-Group/nllab-CFPSmetabolic

max_issues_repo_name: Laohakunakorn-Group/nllab-CFPSmetabolic

max_forks_repo_name: Laohakunakorn-Group/nllab-CFPSmetabolic


text:
function _count_by_state(
  events::EventObservations{T, M},
  state::DiseaseState,
  time::Float64) where {
    T <: DiseaseStateSequence,
    M <: ILM}
  n_ids = 0
  if state == State_I && State_R ∈ T
  # E/I at or before time and I/R after time or never
    for i = 1:individuals(events)
      n_ids += events[state][i] <= time && ((events[State_R][i] > time) || (events[State_R][i] === NaN))
    end
  else
    for i = 1:individuals(events)
      n_ids += events[state][i] <= time
    end
  end
  @debug "$n_ids individual(s) observered in state $state at t = $time"
  return n_ids
end

function _obs_curve(
  events::EventObservations{T, M},
  state::DiseaseState,
  min::Float64,
  max::Float64) where {
    T <: DiseaseStateSequence,
    M <: ILM}
  if min >= max
    @error "Minimum time must be less than maximum time"
  end
  local times
  if state == State_I && State_R ∈ T
    times = events[[State_I; State_R]][:]
  else
    times = events[state]
  end
  times = times[Ref(min) .< times .< Ref(max)]
  sort!(times)
  insert!(times, 1, min)
  push!(times, max)
  counts = _count_by_state.(Ref(events), Ref(state), times)
  return times, counts
end

@recipe function f(
  events::EventObservations{T, M},
  state::DiseaseState,
  min::Float64,
  max::Float64) where {
    T <: DiseaseStateSequence,
    M <: ILM}
  xguide --> "Time"
  yguide --> "N"
  linewidth --> 2.0
  linecolor --> _state_color(state)
  label --> ""
  seriestype --> :steppost
  _obs_curve(events, state, min, max)
end

@recipe function f(
  events::EventObservations{T, M},
  state::DiseaseState) where {
    T <: DiseaseStateSequence,
    M <: ILM}
  events, state, 0.0, maximum(events)
end

@recipe function f(
  events::EventObservations{T, M},
  min::Float64,
  max::Float64) where {
    T <: DiseaseStateSequence,
    M <: ILM}
  @series begin
    linecolor --> _state_color(State_I)
    label --> "I"
    events, State_I, min, max
  end
  if State_R ∈ T
    @series begin
      linecolor --> _state_color(State_R)
      label --> "R"
      events, State_R, min, max
    end
  end
end

@recipe function f(
  events::EventObservations{T, M}) where {
    T <: DiseaseStateSequence,
    M <: ILM}
  events, minimum(events), maximum(events)
end


text:
# Interval
# ========
#
# Base interval types and utilities.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# Note, just to be clear: this shadows IntervalTrees.Interval
"A genomic interval specifies interval with some associated metadata"
struct Interval{T} <: IntervalTrees.AbstractInterval{Int64}
    seqname::String
    first::Int64
    last::Int64
    strand::Strand
    metadata::T
end

function Interval(seqname::AbstractString, first::Integer, last::Integer,
                  strand::Union{Strand,Char}=STRAND_BOTH, metadata=nothing)
    return Interval{typeof(metadata)}(seqname, first, last, strand, metadata)
end

function Interval(seqname::AbstractString, range::UnitRange{T},
                  strand::Union{Strand,Char}=STRAND_BOTH, metadata=nothing) where T<:Integer
    return Interval{typeof(metadata)}(seqname, first(range), last(range), strand, metadata)
end

function BioCore.seqname(i::Interval)
    return i.seqname
end

function BioCore.metadata(i::Interval)
    return i.metadata
end

function strand(i::Interval)
    return i.strand
end

"""
    leftposition(i::Interval)

Return the leftmost position of `i`.
"""
function BioCore.leftposition(i::Interval)
    return i.first
end

"""
    rightposition(i::Interval)

Return the rightmost position of `i`.
"""
function BioCore.rightposition(i::Interval)
    return i.last
end

IntervalTrees.first(i::Interval) = i.first
IntervalTrees.last(i::Interval) = i.last

function Base.isless(a::Interval{T}, b::Interval{T},
                     seqname_isless::Function=isless) where T
    if a.seqname != b.seqname
        return seqname_isless(a.seqname, b.seqname)::Bool
    elseif a.first != b.first
        return a.first < b.first
    elseif a.last != b.last
        return a.last < b.last
    elseif a.strand != b.strand
        return a.strand < b.strand
    else
        return false
    end
end

"""
Check if two intervals are well ordered.

Intervals are considered well ordered if a.seqname <= b.seqnamend and
a.first <= b.first.
"""
function isordered(a::Interval{T}, b::Interval{T},
                   seqname_isless::Function=isless) where T
    if a.seqname != b.seqname
        return seqname_isless(a.seqname, b.seqname)::Bool
    elseif a.first != b.first
        return a.first < b.first
    else
        return true
    end
end

"""
Return true if interval `a` entirely precedes `b`.
"""
function precedes(a::Interval{T}, b::Interval{T},
                  seqname_isless::Function=isless) where T
    return (a.last < b.first && a.seqname == b.seqname) ||
        seqname_isless(a.seqname, b.seqname)::Bool
end

function Base.:(==)(a::Interval{T}, b::Interval{T}) where T
    return a.seqname  == b.seqname &&
           a.first    == b.first &&
           a.last     == b.last &&
           a.strand   == b.strand &&
           a.metadata == b.metadata
end

"Return true if interval `a` overlaps interval `b`, with no consideration to strand"
function BioCore.isoverlapping(a::Interval{S}, b::Interval{T}) where {S, T}
    return a.first <= b.last && b.first <= a.last && a.seqname == b.seqname
end

function Base.show(io::IO, i::Interval)
    if get(io, :compact, false)
        print(io, i.seqname, ":", i.first, "-", i.last, "  ", i.strand,
            "  ", i.metadata === nothing ? "nothing" : i.metadata)
    else
        println(io, summary(i), ':')
        println(io, "  sequence name: ", i.seqname)
        println(io, "  leftmost position: ", i.first)
        println(io, "  rightmost position: ", i.last)
        println(io, "  strand: ", i.strand)
          print(io, "  metadata: ",
            i.metadata === nothing ? "nothing" : i.metadata)
    end
end

function metadatatype(::Type{T}) where T
    return _metadatatype(eltype(T))
end

function metadatatype(x::Any)
    return metadatatype(typeof(x))
end

function _metadatatype(::Type{Interval{T}}) where T
    return T
end

function _metadatatype(::Type{T}) where T
    return T
end


text:
tomlpath = joinpath(@__DIR__, "hbv_config.toml")
config = Wflow.Config(tomlpath)

model = Wflow.initialize_hbv_model(config)
@unpack network = model

model = Wflow.update(model)

# test if the first timestep was written to the CSV file
flush(model.writer.csv_io)  # ensure the buffer is written fully to disk
@testset "CSV output" begin
    row = csv_first_row(model.writer.csv_path)

    @test row.time == DateTime("2000-01-01T00:00:00")
    @test row.Q ≈ 517.9761666607884f0
    @test row.temp_bycoord ≈ 2.965437173843384f0
    @test row.temp_byindex ≈ 1.1716821193695068f0
    @test row.Q_1 ≈ 501.6252231440609f0
    @test row.perc_33 ≈ 2.308000087738037f0
    @test row.perc_34 ≈ 1.8980000019073486f0
    @test row.perc_35 ≈ 2.7100000381469727f0
    @test row.perc_36 ≈ 3.818000078201294f0
    @test row.perc_37 ≈ 2.1440000534057617f0
end

@testset "first timestep" begin
    hbv = model.vertical
    @test hbv.tt[4377] ≈ 0.0
    @test model.clock.iteration == 2
    @test hbv.soilmoisture[4377] ≈ 134.35299682617188f0
    @test hbv.runoff[4377] ≈ 7.406898120121746f0
    @test hbv.soilevap[4377] == 0.0
    @test hbv.snow[4377] == 0.0
end

# run the second timestep
model = Wflow.update(model)

@testset "second timestep" begin
    hbv = model.vertical
    @test hbv.soilmoisture[4377] ≈ 134.35299682617188f0
    @test hbv.runoff[4377] ≈ 4.3533463f0
    @test hbv.soilevap[4377] == 0.0
    @test hbv.snow[4377] == 0.0
end

@testset "overland domain" begin
    q = model.lateral.land.q_av
    land = model.lateral.land
    @test sum(q) ≈ 3278.9423039643552f0
    @test q[10354] ≈ 0.2348646912841589f0
    @test land.volume[10354] ≈ 2057.7314802432425f0
    @test land.inwater[10354] ≈ 0.027351853491789667f0
    @test q[network.land.order[end]] ≈ 0.296794455575309f0
end

@testset "river flow" begin
    q = model.lateral.river.q_av
    @test sum(q) ≈ 51943.24072462076f0
    @test q[651] ≈ 5.693405508134953f0
    @test q[1056] ≈ 8.839306661545436f0
    @test q[network.river.order[end]] ≈ 363.5647479049738f0
end

benchmark = @benchmark Wflow.update(model)
trialmin = BenchmarkTools.minimum(benchmark)

println("HBV Model update")
print_benchmark(trialmin)

Wflow.close_files(model, delete_output = false)


text:
using Distributions

const AllContinuousDistributions = [
    Arcsine, Beta, BetaPrime, Cauchy, Chi, Chisq, Erlang, Exponential,
    FDist, Frechet, Gamma, GeneralizedExtremeValue, GeneralizedPareto,
    Gumbel, InverseGamma, InverseGaussian, Laplace, Levy, LogNormal,
    Logistic, Normal, NormalInverseGaussian, Pareto, Rayleigh,
    SymTriangularDist, TDist, TriangularDist, Uniform, VonMises, Weibull
]

const ContinuousDistributionsThatCanFit = filter(AllContinuousDistributions) do D
    try
        fit_mle(D, rand(10))
        true
    catch _err
        false
    end
end

function fitalldistr{T<:Real}(x::AbstractArray{T}, sortby = :BIC)
    n = length(x)
    results = Any[]
    for D in ContinuousDistributionsThatCanFit
        try
            d = fit_mle(D, x)
            loglik = loglikelihood(d, x)
            k = length(params(d))
            if sortby == :BIC
                criteria = -2 * loglik + k * log(n)
            elseif sortby == :AIC
                criteria = -2 * loglik + 2 * k
            elseif sortby == :AICc
                criteria = -2 * loglik + 2 * k + ((2*k*(k+1))/(n-k-1))
            elseif sortby == :NLogL
                criteria = -loglik
            end
            push!(results, (d, criteria))
        catch _err
        end
    end
    sort(results, by = t -> t[2])
end

function fitall{T<:Real}(x::AbstractArray{T}, sortby = :BIC)
    fitalldistr(x, sortby)[1][1]
end


text:
@testset "Contructor Tests" begin

@testset "DictTen Constructor Tests" begin

    n = 5
    valid_edges = [([1,2,3],1.0),([1,2,2],-1.0),([3,3,3],1.0)]
    unordered_edge = [([2,3,1],1.0)]
    negative_index = [([-1,2,3],1.0)]

    @testset "Hyperedge List Constructor" begin

        @test_throws ErrorException SSSTensor(unordered_edge) #unsorted indices
        @test_throws ErrorException SSSTensor(unordered_edge,n)
        @test_throws ErrorException SSSTensor(unordered_edge,1) # too small cubical dim
        @test_throws ErrorException SSSTensor(valid_edges,1)
        @test_throws ErrorException SSSTensor(negative_index,1)
        @test_throws ErrorException SSSTensor(negative_index,n) #neg index
        @test_throws ErrorException SSSTensor(negative_index)

    end

    @testset "Dense Tensor Constructor" begin

        non_sym_tensor = rand(2,2,2,2)
        sym_tensor = ones(2,2,2,2)
        @test_throws ErrorException SSSTensor(non_sym_tensor)
        @test_throws ErrorException SSSTensor(sym_tensor,1) # too small cubical dim
        @test_throws ErrorException SSSTensor(non_sym_tensor,1)

    end

    @testset "Dictionary Constructor" begin

        D_valid_edges = Dict(valid_edges)
        D_unordered_edge = Dict(unordered_edge)
        D_negative_index = Dict(negative_index)

        @test_throws ErrorException SSSTensor(D_unordered_edge) #unsorted indices
        @test_throws ErrorException SSSTensor(D_unordered_edge,n)
        @test_throws ErrorException SSSTensor(D_unordered_edge,1) # too small cubical dim
        @test_throws ErrorException SSSTensor(D_valid_edges,1)
        @test_throws ErrorException SSSTensor(D_negative_index,1)
        @test_throws ErrorException SSSTensor(D_negative_index,n) #neg index
        @test_throws ErrorException SSSTensor(D_negative_index)

    end

    @test_throws ErrorException SSSTensor()
end

@testset "COOTEN" begin

    valid_indices = [1 2 3; 2 2 3; 1 2 4]
    unsorted_indices = [1 2 3; 2 2 3; 4 2 1]
    zero_indexed_indices = [1 2 3; 2 2 3; 0 0 0]

    valid_values = [1.0,2.0,3.0]

    @testset "COOTEN Constructor" begin

        A = COOTen(valid_indices,rand(3),10)
        @test A.cubical_dimension == 10

         A = COOTen(valid_indices,rand(3))
        @test A.cubical_dimension == maximum(valid_indices)

        #all params
        @test_throws ErrorException COOTen(valid_indices, rand(2),10) #mis-matched dimensions
        @test_throws ErrorException COOTen(unsorted_indices, valid_values,10) #unsorted rows
        @test_throws ErrorException COOTen(zero_indexed_indices, valid_values,10) #zero indexed

        @test_throws ErrorException COOTen(zero_indexed_indices,valid_values,2)#invalid cubicaldim

        #without cubical dimension
        @test_throws ErrorException COOTen(valid_indices, rand(2)) #mis-matched dimensions
        @test_throws ErrorException COOTen(unsorted_indices, valid_values) #unsorted rows
        @test_throws ErrorException COOTen(zero_indexed_indices, valid_values) #zero indexed

        #without weights
        @test_throws ErrorException COOTen(unsorted_indices) #unsorted rows
        @test_throws ErrorException COOTen(zero_indexed_indices) #zero indexed

    end

    @test_throws ErrorException COOTen()

    @testset "COOTEN iterator" begin
        valid_indices = [1 2 3; 2 2 3; 1 2 4]
        valid_values = [1.0,2.0,3.0]
        A = COOTen(valid_indices,valid_values)
        @test length(A) == 3

    end

end

end

#=
"""-----------------------------------------------------------------------------
                     Compute Tensor eigenvectors in MATLAB
-----------------------------------------------------------------------------"""
function MATLAB_eigs(A::SSST.SSSTensor)

    data = Array{Array{Any,1}}(undef,length(A.edges))
    i = 1
    for (indices, val) in A.edges
        data[i] = [indices,SSST.multiplicity_factor(indices),val]
        i += 1
    end

    xs = mxarray(data)

    n = A.cubical_dimension

    mat"""
        disp(xs)
    """

    #=

    mat"""

    %build homogenous polynomial

    mpol('x',$n)
    mpol('f',1)

    f = $xs{1}{2}*$xs{1}{3}*prod(x($xs{1}{1}));

    for i =2:length(xs)
        f = f + $xs{i}{2}*$xs{i}{3}*prod(x($xs{i}{1}));
    end

    %call AReig solver

    [$lmd, $eigvec,  $info] = AReigSTensors(f, [], x, 3, 2)

    """
    println(lmd)
    println(eigvec)

    =#

end
=#


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    @assert precompile(Tuple{typeof(findfont),String})   # time: 0.12886831
    @assert precompile(Tuple{typeof(try_load),String})   # time: 0.033520337
    @assert precompile(Tuple{typeof(renderface),FTFont,Char,Int64})   # time: 0.019107351
end


text:
struct LinearPolicyWithFluxBasis{TP, TB} <: AbstractPolicy where {TP, TB}
    π::TP
    ϕ::TB
end

function params(π::LinearPolicyWithFluxBasis)
    return Flux.params(π.ϕ, params(π.π)...)
end

function paramsvec(π::LinearPolicyWithFluxBasis)
    xs = Zygote.Buffer([])
    m = π.ϕ
    fmap(m) do x
        x isa AbstractArray && push!(xs, x)
        return x
    end
    ps = params(π.π)
    for p in ps
        push!(xs, p)
    end
    return vcat(vec.(copy(xs))...)
end

function paramsfromvec!(π::LinearPolicyWithFluxBasis, θ)
    i = 0
    m = π.ϕ
    fmap(m) do x
        x isa AbstractArray || return x
        x .= reshape(@view θ[i.+(1:length(x))], size(x))
        i += length(x)
        return x
    end
    length(xs) == i || @warn "Expected $(i) params, got $(length(xs))"
    paramsfromvec!(π.π, @view θ[i+1:end])
    return π
end

function rrule(::typeof(paramsfromvec!), π::LinearPolicyWithFluxBasis, θ)
    i = 0
    m = π.ϕ
    m̄ = fmap(m) do x
        x isa AbstractArray || return x
        x .= reshape(@view θ[i.+(1:length(x))], size(x))
        i += length(x)
        return x
    end
    length(xs) == i || @warn "Expected $(i) params, got $(length(xs))"
    paramsfromvec!(π.π, @view θ[i+1:end])

    function paramsfromvec_linearwithflux!(ȳ)
        xs = Zygote.Buffer([])
        dm = ȳ.ϕ
        fmap(dm) do x
            x isa AbstractArray && push!(xs, x)
            return x
        end
        dpi = ȳ.π
        for d in dpi
            push!(xs, p)
        end
        dθ = vcat(vec.(copy(xs))...)
        return NoTangent(), ȳ, dθ
    end
    return π, paramsfromvec_linearwithflux!
end

# function _restructure(m, xs)
#     i = 0
#     m̄ = fmap(m) do x
#         x isa AbstractArray || return x
#         x = reshape(xs[i.+(1:length(x))], size(x))
#         i += length(x)
#         return x
#     end
#     length(xs) == i || @warn "Expected $(i) params, got $(length(xs))"
#     return m̄
# end

# @adjoint function _restructure(m, xs)
#     m̄, numel = _restructure(m, xs), length(xs)
#     function _restructure_pullback(dm)
#         xs′ = destructure(dm)[1]
#         numel == length(xs′) || @warn "Expected $(numel) params, got $(length(xs′))"
#         return (nothing, xs′)
#     end
#     return m̄, _restructure_pullback
# end

# """
#     destructure(m)
# Flatten a model's parameters into a single weight vector.
#     julia> m = Chain(Dense(10, 5, σ), Dense(5, 2), softmax)
#     Chain(Dense(10, 5, σ), Dense(5, 2), softmax)
#     julia> θ, re = destructure(m);
#     julia> θ
#     67-element Vector{Float32}:
#     -0.1407104
#     ...
# The second return value `re` allows you to reconstruct the original network after making
# modifications to the weight vector (for example, with a hypernetwork).
#     julia> re(θ .* 2)
#     Chain(Dense(10, 5, σ), Dense(5, 2), softmax)
# """
# function destructure(m)
#     xs = Zygote.Buffer([])
#     fmap(m) do x
#         x isa AbstractArray && push!(xs, x)
#         return x
#     end
#     return vcat(vec.(copy(xs))...), p -> _restructure(m, p)
# end

function (π::LinearPolicyWithFluxBasis)(s)
    feats = π.ϕ(s)
    return π.π(feats)
end

function logpdf(π::LinearPolicyWithFluxBasis, s, a)
    feats = π.ϕ(s)
    return logpdf(π.π, feats, a)
end

function grad_logpdf(π::LinearPolicyWithFluxBasis, s, a)
    ps = params(π)
    logp, gp = Zygote.pullback(()->logpdf(π, s, a), ps)
    tmp = gp(1)
    ψ = [tmp[p] for p in ps]
    return logp, ψ
end

function sample_with_trace(π::LinearPolicyWithFluxBasis, s)
    ps = params(π)
    action = []
    function f!(res, π)
        d = π(s)

        Zygote.ignore() do
            a = rand(d)
            push!(res, a)
        end
        return logpdf(d, res[1])
    end
    logp, back = Zygote.pullback(()->f!(action, π), ps)
    tmp = back(1)
    ψ = [tmp[p] for p in ps]
    return action[1], logp, ψ
end


text:
abstract type AbstractMDP <: AbstractDecisionMakingProblem end

struct MDP{TS,TA,TP,TR,TG,TD,TM,TF} <: AbstractMDP where {TS, TA, TP, TR, TG, TD, TM, TF}
    S::TS
    A::TA
    p::TP
    r::TR
    γ::TG
    d0::TD
    meta::TM
    render::TF
end

struct POMDP{TS,TA,TX,TP,TO,TR,TG,TD,TM,TF} <: AbstractMDP where {TS, TA, TX, TP, TO, TR, TG, TD, TM, TF}
    S::TS
    A::TA
    X::TX
    p::TP
    obs::TO
    r::TR
    γ::TG
    d0::TD
    meta::TM
    render::TF
end

struct SequentialProblem{TS,TX,TA,TP,TD,TM,TF} <: AbstractMDP where {TS, TX, TA, TP, TR, TG, TD, TM, TF}
    S::TS
    X::TX
    A::TA
    p::TP
    d0::TD
    meta::TM
    render::TF
end

export AbstractMDP, MDP, POMDP
export create_simple_chain, create_minimum_time_chain_task
export create_finitetime_cartpole, create_cartpole_balancetask

function sample(prob::SequentialProblem, s, a)
    # s′, x, r, γ = prob.p(s,a)
    # return s′, x, r, γ
    ret = prob.p(s,a)
    return ret
end

function sample(prob::MDP, s, a)
    s′ = prob.p(s,a)
    r = prob.r(s,a,s′)
    γ = prob.γ(s,a,s′)
    return s′, r, γ
end

function sample(prob::POMDP, s, a)
    s′ = prob.p(s,a)
    x = prob.obs(s′)
    r = prob.r(s,a,s′)
    γ = prob.γ(s,a,s′)
    return s′, x, r, γ
end

function sample_trajectory!(prob::MDP, τ, π)
    s = prob.d0()
    done = false
    while !done
        d = π(s)
        a = rand(d)
        logp = logpdf(d,a)
        s′, r, γ = sample(prob,s,a)
        push!(τ, s, a, logp, r)
        if γ == 0
            done = true
        end
        s = s′
    end
end

function sample_trajectory!(prob::POMDP, τ, π)
    s = prob.d0()
    x = prob.obs(s)
    done = false
    while !done
        d = π(x)
        a = rand(d)
        logp = logpdf(d,a)
        s′, x′, r, γ = sample(prob,s,a)
        push!(τ, x, a, logp, r)
        if γ == 0
            done = true
        end
        s = s′
        x = x′
    end
end

function sample_trajectory!(prob::SequentialProblem, τ, π)
    s, x = prob.d0()
    done = false
    while !done
        d = π(x)
        a = rand(d)
        logp = logpdf(d,a)
        s′, x′, r, γ = sample(prob,s,a)
        push!(τ, x, a, logp, r)
        if γ == 0
            done = true
        end
        s = s′
        x = x′
    end
end

function sample_objective(prob::MDP, π)
    s = prob.d0()
    G = 0.0
    γt = 1.0
    done = false
    while !done
        d = π(s)
        a = rand(d)
        logp = logpdf(d,a)
        s′, r, γ = sample(prob,s,a)
        G += γt * r
        γt *= γ
        if γ == 0
            done = true
        end
        s = s′
    end
    return G
end

function sample_objective(prob::POMDP, π)
    s = prob.d0()
    x = prob.obs(s)
    G = 0.0
    γt = 1.0
    done = false
    while !done
        d = π(x)
        a = rand(d)
        logp = logpdf(d,a)
        s′, x′, r, γ = sample(prob,s,a)
        G += γt * r
        γt *= γ
        if γ == 0
            done = true
        end
        s = s′
        x = x′
    end
    return G
end

function sample_objective(prob::SequentialProblem, π)
    s, x = prob.d0()
    G = 0.0
    γt = 1.0
    done = false
    while !done
        d = π(x)
        a = rand(d)
        logp = logpdf(d,a)
        s′, x′, r, γ = sample(prob,s,a)
        G += γt * r
        γt *= γ
        if γ == 0
            done = true
        end
        s = s′
        x = x′
    end
    return G
end

function meta_information(prob::AbstractMDP)
    return prob.meta
end

function render(prob::AbstractMDP, state, cleardata=false)
    return prob.render(state, cleardata)
end

include("chains.jl")
include("cartpole.jl")


text:
module Flatten

export to_tuple, to_vector, from_tuple, from_vector, wrap

function field_expressions(T, expr::Union{Expr, Symbol})
    expressions = Expr[]
    for (i, field) in enumerate(fieldnames(T))
        field_expr = Expr(:., expr, Expr(:quote, field))
        append!(expressions, field_expressions(fieldtype(T, i), field_expr))
    end
    expressions
end

field_expressions{T2 <: AbstractArray}(T::Type{T2}, expr::Union{Expr, Symbol}) = error("Cannot flatten variable-length objects like arrays. Replace any arrays with tuples if possible.")

function field_expressions{T2 <: Tuple}(T::Type{T2}, expr::Union{Expr, Symbol})
    expressions = Expr[]
    for i in 1:length(T.types)
        field_expr = Expr(:ref, expr, i)
        append!(expressions, field_expressions(fieldtype(T, i), field_expr))
    end
    expressions
end

function field_expressions(T::Type{Any}, expr::Union{Expr, Symbol})
    [expr]
end

function field_expressions{T2 <: Number}(T::Type{T2}, expr::Union{Expr, Symbol})
    [expr]
end

function to_tuple_internal(T)
    expr = Expr(:tuple)
    append!(expr.args, field_expressions(T, :(T)))
    expr
end

@generated function to_tuple(T)
    to_tuple_internal(T)
end

function all_field_types(T)
    field_types = DataType[]
    for i in 1:length(fieldnames(T))
        append!(field_types, all_field_types(fieldtype(T, i)))
    end
    field_types
end

function all_field_types{T2 <: Number}(T::Type{T2})
    [T]
end

function all_field_types{T2 <: Tuple}(T::Type{T2})
    field_types = DataType[]
    for i in 1:length(T.types)
        append!(field_types, all_field_types(T.types[i]))
    end
    field_types
end

function to_vector_internal(T)
    field_types = all_field_types(T)
    num_elements = length(field_types)
    element_type = reduce(promote_type, field_types)
    expr = quote
        v = Array{$(element_type)}($(num_elements))
    end
    for (i, field_expr) in enumerate(field_expressions(T, :(T)))
        push!(expr.args, Expr(:(=), Expr(:ref, :v, i), field_expr))
    end
    push!(expr.args, :v)
    return expr
end

@generated function to_vector(T)
    to_vector_internal(T)
end

type Counter
    value::Int

    Counter() = new(1)
end

function construct(T, counter)
    expr = Expr(:call, T)
    for subtype in T.types
        push!(expr.args, construct(subtype, counter))
    end
    expr
end

function construct{T <: Tuple}(::Type{T}, counter)
    expr = Expr(:tuple)
    for subtype in T.types
        push!(expr.args, construct(subtype, counter))
    end
    expr
end

function construct_element(counter)
    expr = Expr(:ref, :data, counter.value)
    counter.value += 1
    expr
end

function construct(T::TypeVar, counter)
    construct_element(counter)
end

function construct{T <: Number}(::Type{T}, counter)
    construct_element(counter)
end

function construct(::Type{Any}, counter)
    construct_element(counter)
end

function from_tuple_internal(T, data)
    construct(T, Counter())
end

@generated function from_tuple{T}(::Type{T}, data)
    from_tuple_internal(T, data)
end

from_vector(T, data) = from_tuple(T, data)

function wrap(func, InputType)
	return x -> to_vector(func(from_vector(InputType, x)))
end

end # module


text:
#Written By Johnathan Bizzano

using ..Functions
using ...HyperSphere

export AbstractMathFun, SingleVarMathmaticalFunction,
MFCollection, MFCall, verify_function, emit, FunctionEmitter, FunctionCall, functionheader, trainer, simplify, is_constant

abstract type AbstractMathFun{T, N} <: AbstractObject end
abstract type AbstractMathmaticalEquation{T, N} <: AbstractMathFun{T, N} end
abstract type SingleVarMathmaticalFunction{T} <: AbstractMathFun{T, 1} end

const MFCollection = Array{AbstractMathFun, 1}

functionheader(x::AbstractMathFun, var_count::Integer) = "f(x_" * join(1:var_count, ", x_") * ") = "

function (f::AbstractMathFun{T, N})(args...; Data_Type=T) where {T, N}
    f(MFCall{T, N}(tuple(args...,), Data_Type))
end

function (f::AbstractMathFun{T, N})(args::NTuple{N, T}; Data_Type=T) where {T, N}
    f(MFCall{T, N}(args, Data_Type))
end

verify_function(f::AbstractMathFun) = ()
simplify(f::AbstractMathFun) = f
simplify(f::AbstractMathFun, emit_constants::Bool) = f
is_constant(f::AbstractMathFun, emit_constants::Bool) = f isa ConstFunction && (emit_constants || f.force_emit)

struct MFCall{T, N}
    args::NTuple{N, T}
    Data_Type::Type

    MFCall{T, N}(args::NTuple{N, T}, Data_Type::Type) where {T, N} = new{T, N}(args, Data_Type)

    function (call::MFCall)(arg)
        call.Data_Type(arg)
    end
end

struct FunctionEmitter{T}
    ConstantPool::Vector{T}
    EmitConstants::Bool

    FunctionEmitter{T}(EmitConstants::Bool) where T = new{T}(zeros(T, 0), EmitConstants)

    function (emitter::FunctionEmitter{T})(constant::T, force_emit::Bool) where T
        if emitter.EmitConstants || force_emit
            return 0
        else
            push!(emitter.ConstantPool, constant)
            return length(emitter.ConstantPool)
        end
    end
end

struct FunctionCall{T, N} <: AbstractMathFun{T, N}
    method::Function
    ConstantPool::Array{T, 1}
    expr::Expr

    FunctionCall{T, N}(constantPool::Array{T, 1}, expr::Expr) where {T, N} = new{T, N}(eval(expr), constantPool)

    FunctionCall{T, N}(method::Function, constantPool::Array{T, 1}, expr::Expr) where {T, N} = new{T, N}(method, constantPool, expr)

    (f::FunctionCall)(args...) = f(collect(args), f.ConstantPool)

    (f::FunctionCall)(args::AbstractArray) = f(args, f.ConstantPool)

    function (f::FunctionCall{T, N})(args::AbstractArray, constants::AbstractArray) where {T, N}
        if length(constants) != 0
            if N != 0
                return f.method(args, constants)
            else
                return f.method(constants)
            end
        else
            if N != 0
                return f.method(args)
            else
                return f.method()
            end
        end
    end

    function Base.string(f::FunctionCall{T, N}) where {T, N}
        str = "Function{$T, $N}:" * string(f.method)
        if length(f.ConstantPool) != 0
            if N != 0
                str *= "(ARGS, CONSTANTS). With Constant Pool:" * string(f.ConstantPool)
            else
                str *= "(CONSTANTS). With Constant Pool:" * string(f.ConstantPool)
            end
        else
            if N != 0
                str *= "(ARGS)"
            else
                str *= "()"
            end
        end
        str
    end
end

function emit(f::AbstractMathFun{T, N}, Name::Symbol; EmitConstants::Bool = true) where {T, N}
    emitter = FunctionEmitter{T}(EmitConstants)
    funArgs = Expr(:call, Name)
    N != 0 && push!(funArgs.args, :ARGS)

    expr = Expr(:function, funArgs)
    EmitConstants || push!(funArgs.args, :CONSTANTS)
    push!(expr.args, Expr(:block, emit(simplify(f), emitter)))
    return FunctionCall{T, N}(eval(expr), emitter.ConstantPool, expr)
end

emit(f::AbstractMathFun{T, N}, emitter::FunctionEmitter{T}) where {T, N} = ()
Functions.string(f::AbstractMathFun{T, N}, emit_constants::Bool) where {T, N} = string(f, FunctionEmitter{T}(emit_constants))


text:
@testset "linear_regression.jl internal function" begin
    # Ŷ = 0.1 + 0.3 x₁ + 0.6 X₂ + 0.9 x₃
    x₁ = collect(1:10)
    x₂ = collect(11:20)
    x₃ = collect(21:30)
    ŷ = 0.1 .+ 0.3x₁ .+ 0.5x₂ .+ 0.7x₃
    df = DataFrame(Ŷ=ŷ, X₁=x₁, X₂=x₂, X₃=x₃)

    lrm = LinearRegressionModel(df, :Ŷ, [:X₁, :X₂, :X₃], [0.1, 0.3, 0.5, 0.7])

    @test ML101.g(lrm, 1) == 0.1 + 0.3*1 + 0.5*11 + 0.7*21
    @test ML101.loss(lrm) == 0
end

@testset "linear_regression.jl multipal regression" begin
    # Ŷ = 0.1 + 0.3 x₁ + 0.6 X₂ + 0.9 x₃
    x₁ = rand(10)
    x₂ = rand(10)
    x₃ = rand(10)
    ŷ = 0.1 .+ 0.3x₁ .+ 0.5x₂ .+ 0.7x₃
    df = DataFrame(Ŷ=ŷ, X₁=x₁, X₂=x₂, X₃=x₃)

    lrm = LinearRegressionModel(df, :Ŷ, [:X₁, :X₂, :X₃])
    fit!(lrm)

    @test all(isapprox.(lrm.argv, [0.1, 0.3, 0.5, 0.7], atol=1e-2))
    @test isapprox(
        predict(lrm, [0.1, 0.2, 0.3]),
        0.1*1 + 0.3*0.1 + 0.5*0.2 + 0.7*0.3,
        atol=1e-2
    )
end

@testset "linear_regression.jl linear regression" begin
    # Ŷ = 0.1 + 0.3 x₁
    x₁ = rand(10)
    ŷ = 0.1 .+ 0.3x₁
    df = DataFrame(Ŷ=ŷ, X₁=x₁)

    lrm = LinearRegressionModel(df, :Ŷ, :X₁)
    fit!(lrm)

    @test all(isapprox.(lrm.argv, [0.1, 0.3], atol=1e-2))
    @test isapprox(predict(lrm, 0.1), 0.1*1 + 0.3*0.1, atol=1e-2)
end


text:
using RasterDataSources, Test, Dates
using RasterDataSources: rasterurl, rastername, rasterpath

@testset "CHELSA Future BioClim CMIP5" begin
    @test rastername(CHELSA{Future{BioClim,CMIP5,CCSM4,RCP26}}, 5; date=Date(2050)) == "CHELSA_bio_mon_CCSM4_rcp26_r1i1p1_g025.nc_5_2041-2060_V1.2.tif"

    bioclim_path = joinpath(ENV["RASTERDATASOURCES_PATH"], "CHELSA", "Future", "BioClim", "RCP26", "CCSM4")
    @test rasterpath(CHELSA{Future{BioClim,CMIP5,CCSM4,RCP26}}) == bioclim_path

    raster_path = joinpath(bioclim_path, "CHELSA_bio_mon_CCSM4_rcp26_r1i1p1_g025.nc_5_2041-2060_V1.2.tif")

    @test getraster(CHELSA{Future{BioClim,CMIP5,CCSM4,RCP26}}, 5; date=Date(2050)) == raster_path
    @test getraster(CHELSA{Future{BioClim,CMIP5,CCSM4,RCP26}}, (5,); date=Date(2050)) == (bio5=raster_path,)
    @test getraster(CHELSA{Future{BioClim,CMIP5,CCSM4,RCP26}}, [5]; date=Date(2050)) == (bio5=raster_path,)
    @test isfile(raster_path)
end

@testset "CHELSA Future BioClim CMIP6" begin
    bioclim_name = "CHELSA_bio5_2041-2070_mri-esm2-0_ssp126_V.2.1.tif"
    @test rastername(CHELSA{Future{BioClim,CMIP6,MRIESM2,SSP126}}, 5; date=Date(2050)) == bioclim_name

    bioclim_path = joinpath(ENV["RASTERDATASOURCES_PATH"], "CHELSA", "Future", "BioClim", "SSP126", "MRIESM2")
    @test rasterpath(CHELSA{Future{BioClim,CMIP6,MRIESM2,SSP126}}) == bioclim_path

    raster_path = joinpath(bioclim_path, bioclim_name)
    raster_path2 = joinpath(bioclim_path, "CHELSA_bio5_2071-2100_mri-esm2-0_ssp126_V.2.1.tif")
    @test getraster(CHELSA{Future{BioClim,CMIP6,MRIESM2,SSP126}}, 5; date=Date(2050)) == raster_path
    @test getraster(CHELSA{Future{BioClim,CMIP6,MRIESM2,SSP126}}, (5,); date=Date(2050)) == (bio5=raster_path,)
    @test getraster(CHELSA{Future{BioClim,CMIP6,MRIESM2,SSP126}}, [5]; date=Date(2050)) == (bio5=raster_path,)
    @test getraster(CHELSA{Future{BioClim,CMIP6,MRIESM2,SSP126}}, (5,); date=[Date(2050)]) ==
        [(bio5=raster_path,)]
    @test isfile(raster_path)
end

@testset "CHELSA Future Climate CMIP5" begin
    tmax_name = "CHELSA_tasmax_mon_CCSM4_rcp60_r1i1p1_g025.nc_7_2041-2060_V1.2.tif"
    @test rastername(CHELSA{Future{Climate,CMIP5,CCSM4,RCP60}}, :tmax; date=Date(2050), month=7) == tmax_name

    climate_path = joinpath(ENV["RASTERDATASOURCES_PATH"], "CHELSA", "Future", "Climate", "RCP60", "CCSM4")
    @test rasterpath(CHELSA{Future{Climate,CMIP5,CCSM4,RCP60}}) == climate_path

    raster_path = joinpath(climate_path, tmax_name)
    @test rasterpath(CHELSA{Future{Climate,CMIP5,CCSM4,RCP60}}, :tmax; date=Date(2050), month=7) == raster_path
    @test rasterurl(CHELSA{Future{Climate,CMIP5,CCSM4,RCP45}}, :prec; date=Date(2050), month=6) |> string ==
        "https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/cmip5/2041-2060/prec/CHELSA_pr_mon_CCSM4_rcp45_r1i1p1_g025.nc_6_2041-2060.tif"
    @test rasterurl(CHELSA{Future{Climate,CMIP5,CCSM4,RCP45}}, :tmin; date=Date(2050), month=1) |> string ==
        "https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V1/cmip5/2041-2060/tmin/CHELSA_tasmin_mon_CCSM4_rcp45_r1i1p1_g025.nc_1_2041-2060_V1.2.tif"
    @test getraster(CHELSA{Future{Climate,CMIP5,CCSM4,RCP60}}, :tmax; date=Date(2050), month=7) == raster_path
    @test getraster(CHELSA{Future{Climate,CMIP5,CCSM4,RCP60}}, [:tmax]; date=Date(2050), month=7) == (tmax=raster_path,)
    @test getraster(CHELSA{Future{Climate,CMIP5,CCSM4,RCP60}}, (:tmax,); date=Date(2050), month=7:7) == [(tmax=raster_path,)]
    @test getraster(CHELSA{Future{Climate,CMIP5,CCSM4,RCP60}}, :tmax; date=[Date(2050)], month=7:7) == [[raster_path]]
    @test isfile(raster_path)
end

@testset "CHELSA Future Climate CMIP6" begin
    date_name =  "CHELSA_gfdl-esm4_r1i1p1f1_w5e5_ssp585_tas_01_2011_2040_norm.tif"
    date_name2 = "CHELSA_gfdl-esm4_r1i1p1f1_w5e5_ssp585_tas_01_2041_2070_norm.tif"
    date_name3 = "CHELSA_gfdl-esm4_r1i1p1f1_w5e5_ssp585_tas_01_2071_2100_norm.tif"
    month_name =  "CHELSA_gfdl-esm4_r1i1p1f1_w5e5_ssp585_tas_01_2011_2040_norm.tif"
    month_name2 = "CHELSA_gfdl-esm4_r1i1p1f1_w5e5_ssp585_tas_02_2011_2040_norm.tif"
    @test rastername(CHELSA{Future{Climate,CMIP6,GFDLESM4,SSP585}}, :temp; date=Date(2030), month=1) == date_name

    climate_path = joinpath(ENV["RASTERDATASOURCES_PATH"], "CHELSA", "Future", "Climate", "SSP585", "GFDLESM4")
    @test rasterpath(CHELSA{Future{Climate,CMIP6,GFDLESM4,SSP585}}) == climate_path
    date_path = joinpath(climate_path, date_name)
    date_path2 = joinpath(climate_path, date_name2)
    date_path3 = joinpath(climate_path, date_name3)
    month_path = joinpath(climate_path, month_name)
    month_path2 = joinpath(climate_path, month_name2)
    @test rasterpath(CHELSA{Future{Climate,CMIP6,GFDLESM4,SSP585}}, :temp; date=Date(2030), month=1) == date_path
    @test rasterpath(CHELSA{Future{Climate,CMIP6,GFDLESM4,SSP585}}, :temp; date=Date(2030), month=1) == date_path
    date_url = "https://os.zhdk.cloud.switch.ch/envicloud/chelsa/chelsa_V2/GLOBAL/climatologies/2011-2040/GFDL-ESM4/ssp585/tas/" * date_name
    @test rasterurl(CHELSA{Future{Climate,CMIP6,GFDLESM4,SSP585}}, :temp; date=Date(2030), month=1) |> string == date_url
    @test getraster(CHELSA{Future{Climate,CMIP6,GFDLESM4,SSP585}}, :temp; date=Date(2030), month=1) == date_path
    @test getraster(CHELSA{Future{Climate,CMIP6,GFDLESM4,SSP585}}, (:temp,); date=Date(2030), month=1) == (temp=date_path,)
    @test getraster(CHELSA{Future{Climate,CMIP6,GFDLESM4,SSP585}}, :temp; date=Date(2030), month=1) == date_path
    # Month is the inner vector
    @test getraster(CHELSA{Future{Climate,CMIP6,GFDLESM4,SSP585}}, [:temp]; date=[Date(2030)], month=1:2) ==
        [[(temp=month_path,), (temp=month_path2,)]]
    @test getraster(CHELSA{Future{Climate,CMIP6,GFDLESM4,SSP585}}, [:temp]; date=(Date(2030), Date(2090)), month=1:1) ==
        [[(temp=date_path,)], [(temp=date_path2,)], [(temp=date_path3,)]]

    @test isfile(date_path)
    @test isfile(date_path2)
    @test isfile(date_path3)
    @test isfile(month_path)
    @test isfile(month_path2)
end

max_forks_repo_name: cesaraustralia/RasterDataSources


text:
# A dummy wrapper around an out-of-place FFT plan to make it appear in-place
# This is needed for some generic FFT implementations, which do not have in-place plans
struct DummyInplace{opFFT}
    fft::opFFT
end
LinearAlgebra.mul!(Y, p::DummyInplace, X) = (Y .= mul!(similar(X), p.fft, X))
LinearAlgebra.ldiv!(Y, p::DummyInplace, X) = (Y .= ldiv!(similar(X), p.fft, X))

import Base: *, \, length
*(p::DummyInplace, X) = p.fft * X
\(p::DummyInplace, X) = p.fft \ X
length(p::DummyInplace) = length(p.fft)


text:
"""
    get_d_hsml( dx::Real, dy::Real,
                   hsml_inv::Real )

Computes the distance in 2D to the pixel center in units of the kernel support.
"""
@inline function get_d_hsml(dx::T, dy::T, hsml_inv::T) where T
    √( dx*dx + dy*dy ) * hsml_inv
end

"""
    function calculate_index(i::Integer, j::Integer, x_pixels::Integer)

Calculates the index of a flattened 2D image array.
"""
function calculate_index(i::T, j::T, x_pixels::T) where T
    return floor(T, i * x_pixels + j ) + 1
end

"""
    function calculate_weights_2D(  wk::Array{<:Real,1},
                                    iMin::Integer, iMax::Integer,
                                    jMin::Integer, jMax::Integer,
                                    x::Real, y::Real, hsml::Real, hsml_inv::Real,
                                    kernel::SPHKernel,
                                    x_pixels::Integer )

Calculates the kernel- and geometric weights of the pixels a particle contributes to.
"""
@fastmath function calculate_weights( wk::Vector{Float64},
                                      iMin::Integer, iMax::Integer,
                                      jMin::Integer, jMax::Integer,
                                      x::Real, y::Real, hsml::Real, hsml_inv::Real,
                                      kernel::SPHKernel,
                                      x_pixels::Integer )

    is_undersampled = false

    if hsml <= 1.0
        is_undersampled = true
    end

    distr_weight = 0.0
    n_distr_pix  = 0

    @inbounds for i = iMin:iMax
        x_dist, dx = get_x_dx(x, hsml, i)

        for j = jMin:jMax
            y_dist, dy = get_x_dx(y, hsml, j)

            # index in flattened 2D array
            idx = calculate_index(i, j, x_pixels)

            dxdy = dx * dy

            if is_undersampled

                wk[idx]       = dxdy
                distr_weight += dxdy
                n_distr_pix  += 1

            else # is_undersampled

                u = get_d_hsml(x_dist, y_dist, hsml_inv)

                # check if inside the kernel
                if u <= 1.0
                    # evaluate kernel
                    wk[idx]       = 𝒲₂(kernel, u, hsml_inv)
                    wk[idx]      *= dxdy
                    distr_weight += wk[idx]
                    n_distr_pix  += 1
                else
                    wk[idx] = 0.0
                end # u < 1.0

            end # is_undersampled

        end # j
    end # i

    return wk, n_distr_pix, distr_weight
end

"""
    get_quantities_2D( pos, weight, hsml,
                       rho, m, len2pix::T) where T

Helper function to convert quantities to pixel units and the correct data type.
"""
function get_quantities_2D( pos, weight, hsml,
                            rho, m, len2pix::T) where T

    hsml    *= T(len2pix)
    hsml_inv = T(1.0)/hsml
    area     = π * hsml^2

    rho *= T(1.0/(len2pix*len2pix*len2pix))
    dz   = m / rho / area

    return T.(pos), T(weight), hsml, hsml_inv, area, dz
end

"""
   sphMapping_2D( Pos, HSML,
                  M, Rho,
                  Bin_Q, Weights;
                  param::mappingParameters, kernel::SPHKernel,
                  show_progress::Bool=false )

Underlying function to map SPH data to a 2D grid.
"""
function sphMapping_2D( Pos, HSML,
                        M, Rho,
                        Bin_Q, Weights;
                        param::mappingParameters, kernel::SPHKernel,
                        show_progress::Bool=false )

    N = size(M,1)  # number of particles

    # max number of pixels over which the particle can be distributed
    N_distr = param.Npixels[1] * param.Npixels[2]

    image = zeros(Float64, N_distr, 2)

    if param.periodic
        k_start = 0
    else
        k_start = 7
    end

    # allocate arrays for weights
    wk = zeros(Float64, N_distr)

    if show_progress
        P = Progress(N)
    end

    # loop over all particles
    @inbounds for p = 1:N

        bin_q = Float64(Bin_Q[p])

        if bin_q == 0.0
            continue
        end

        _pos, weight, hsml, hsml_inv, area, dz = get_quantities_2D(Pos[:,p], Weights[p], HSML[p], Rho[p], M[p], param.len2pix)

        for k = k_start:7

            # simplify position quantities for performance
            x, y, z, skip_k = get_xyz( _pos, HSML[p], k, param)

            if skip_k
                continue
            end

            # calculate relevant pixel range
            iMin, iMax = pix_index_min_max( x, hsml, param.Npixels[1] )
            jMin, jMax = pix_index_min_max( y, hsml, param.Npixels[2] )

            wk, n_distr_pix, distr_weight = calculate_weights(wk, iMin, iMax, jMin, jMax,
                                                              x, y, hsml, hsml_inv,
                                                              kernel,
                                                              param.Npixels[1])

            # skip if the particle is not contained in the image
            # ( should only happen with periodic boundary conditions )
            if n_distr_pix == 0
                continue
            end

            weight_per_pix  = n_distr_pix / distr_weight
            kernel_norm     = area / n_distr_pix
            area_norm       = kernel_norm * weight_per_pix * weight * dz

            @inbounds for i = iMin:iMax, j = jMin:jMax
                idx = calculate_index(i, j, param.Npixels[1])

                image[idx,1], image[idx,2] = update_image( image[idx,1], image[idx,2],
                                                           wk[idx], bin_q, area_norm)

            end # i, j
        end # k

         # update for ProgressMeter
        if show_progress
            next!(P)
        end
    end # p

    return image

end # function


text:
# debug functionality

isdebug(group, mod=CUDA) =
    Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, group, mod) !== nothing


text:
using KTOptical
using Test

@testset "KTOptical.jl" begin
    # Write your tests here.
end


text:
### Day 15
## Timing is Everything
## Author: Thanasis Georgiou

# Hardcoded disks
disk1(t) = (t + 15) % 17 == 0
disk2(t) = (t + 2) % 3 == 0
disk3(t) = (t + 4) % 19 == 0
disk4(t) = (t + 2) % 13 == 0
disk5(t) = (t + 2) % 7 == 0
disk6(t) = (t + 0) % 5 == 0
disk7(t) = (t + 0) % 11 == 0

# Find time that works for problem 1
solution1 = 0
solution2 = 0
t = 0
while solution1 == 0 || solution2 == 0
  if solution1 == 0 && disk1(t + 1) && disk2(t + 2) && disk3(t + 3) && disk4(t + 4) && disk5(t + 5) && disk6(t + 6)
    solution1 = t
  end
  if solution2 == 0 && disk1(t + 1) && disk2(t + 2) && disk3(t + 3) && disk4(t + 4) && disk5(t + 5) && disk6(t + 6) && disk7(t + 7)
    solution2 = t
  end
  t += 1
end

println("Solutions: $solution1, $solution2")