text:
using IRTools.Inner: Slot

function logprob!(ir)
  vars = []
  vartuple = argument!(ir)
  target = Slot(:target)
  pushfirst!(ir, :($target = 0.0))
  for (v, st) in ir
    if iscall(st.expr, rand)
      T = exprtype(ir, v)
      ir[v] = Expr(:call, getindex, vartuple, length(vars)+1)
      prob = insertafter!(ir, v, Expr(:call, logpdf, st.expr.args[2], v))
      target′ = insertafter!(ir, prob, Expr(:call, +, target, prob))
      insertafter!(ir, target′, :($target = $target′))
      push!(vars, (v, T))
    elseif iscall(st.expr, observe)
      prob = insert!(ir, v, Expr(:call, ifelse, st.expr.args[2], 0.0, -Inf))
      target′ = insert!(ir, v, Expr(:call, +, target, prob))
      insert!(ir, v, :($target = $target′))
      delete!(ir, v)
    end
  end
  return!(ir, push!(ir, Expr(:call, tuple, returnvalue(blocks(ir)[end]), target)))
  ssa!(ir)
  return ir, vars
end

logprob(ir) = logprob!(copy(ir))


text:
# This file is part of Kpax3. License is MIT.

function readposteriork(fileroot::String)
  fp = string(fileroot, "_posterior_k.csv")
  d = DelimitedFiles.readdlm(fp, ',', String, header=true)

  len = size(d[1], 1)

  k = zeros(Int, len)
  pk = zeros(Float64, len)

  for i in 1:len
    k[i] = parse(Int, d[1][i, 1])
    pk[i] = parse(Float64, d[1][i, 2])
  end

  (k, pk)
end

function readposteriorP(fileroot::String)
  fp = string(fileroot, "_posterior_R.csv")
  d = DelimitedFiles.readdlm(fp, ',', String, header=true)

  len = size(d[1], 1)
  n = convert(Int, (1 + sqrt(1 + 8 * len)) / 2)

  id = fill("", n)
  P = ones(Float64, n, n)

  id[1] = d[1][1, 1]
  idx = 0
  for j in 1:(n - 1)
    id[j + 1] = d[1][j, 2]

    for i in (j + 1):n
      idx += 1
      P[i, j] = P[j, i] = parse(Float64, d[1][idx, 3])
    end
  end

  (id, P)
end

function readposteriorC(fileroot::String)
  fp = string(fileroot, "_posterior_C.csv")
  d = DelimitedFiles.readdlm(fp, ',', String, header=true)

  len = size(d[1], 1)

  site = zeros(Int, len)
  aa = fill("", len)
  freq = zeros(Float64, len)
  C = zeros(Float64, 3, len)

  for b in 1:len
    site[b] = parse(Int, d[1][b, 1])
    aa[b] = d[1][b, 2]
    freq[b] = parse(Float64, d[1][b, 3])
    C[1, b] = parse(Float64, d[1][b, 4])
    C[2, b] = parse(Float64, d[1][b, 5])
    C[3, b] = parse(Float64, d[1][b, 6])
  end

  (site, aa, freq, C)
end

ext: jl


text:
"""
operations
"""

function Base.:+(x::AbstractNumberSymbolic{Tv}) where Tv
    SymbolicExpression{Tv}(x.params, :+, [x], [])
end

function Base.:-(x::AbstractNumberSymbolic{Tv}) where Tv
    SymbolicExpression{Tv}(x.params, :-, [x], [])
end

for op in [:+, :-, :*, :/, :^]
    @eval function Base.$op(x::AbstractNumberSymbolic{Tx}, y::AbstractNumberSymbolic{Ty}) where {Tx<:Number,Ty<:Number}
        s = union(x.params, y.params)
        Tv = promote_type(Tx,Ty)
        SymbolicExpression{promote_type(Tx,Ty)}(s, $(Expr(:quote, op)), [x, y], [])
    end

    @eval function Base.$op(x::Tx, y::AbstractNumberSymbolic{Ty}) where {Tx<:Number,Ty<:Number}
        Tv = promote_type(Tx,Ty)
        SymbolicExpression{promote_type(Tx,Ty)}(y.params, $(Expr(:quote, op)), [x, y], [])
    end

    @eval function Base.$op(x::AbstractNumberSymbolic{Tx}, y::Ty) where {Tx<:Number,Ty<:Number}
        Tv = promote_type(Tx,Ty)
        SymbolicExpression{promote_type(Tx,Ty)}(x.params, $(Expr(:quote, op)), [x, y], [])
    end
end

function Base.exp(x::AbstractNumberSymbolic{Tv}) where Tv
    SymbolicExpression{Tv}(x.params, :exp, [x], [])
end

function Base.log(x::AbstractNumberSymbolic{Tv}) where Tv
    SymbolicExpression{Tv}(x.params, :log, [x], [])
end

function Base.sqrt(x::AbstractNumberSymbolic{Tv}) where Tv
    SymbolicExpression{Tv}(x.params, :sqrt, [x], [])
end

###

function Base.sum(x::Vector{<:AbstractNumberSymbolic{T}}) where {T<:Number}
    sum(convert(AbstractVectorSymbolic{T}, x))
end

function Base.sum(x::AbstractVectorSymbolic{T}) where {T<:Number}
    SymbolicExpression{T}(x.params, :sum, [x], [])
end

function dot(x::Vector{<:AbstractNumberSymbolic{Tx}}, y::Vector{<:AbstractNumberSymbolic{Ty}}) where {Tx<:Number,Ty<:Number}
    dot(convert(AbstractVectorSymbolic{Tx}, x), convert(AbstractVectorSymbolic{Ty}, y))
end

function dot(x::AbstractVectorSymbolic{Tx}, y::Vector{<:AbstractNumberSymbolic{Ty}}) where {Tx<:Number,Ty<:Number}
    dot(x, convert(AbstractVectorSymbolic{Ty}, y))
end

function dot(x::Vector{<:AbstractNumberSymbolic{Tx}}, y::AbstractVectorSymbolic{Ty}) where {Tx<:Number,Ty<:Number}
    dot(convert(AbstractVectorSymbolic{Tx}, x), y)
end

function dot(x::AbstractVectorSymbolic{Tx}, y::AbstractVectorSymbolic{Ty}) where {Tx<:Number,Ty<:Number}
    Tv = promote_type(Tx,Ty)
    s = union(x.params, y.params)
    SymbolicExpression{Tv}(s, :dot, AbstractVectorSymbolic{Tv}[x, y], [])
end

for op in [:+, :-]
    @eval function Base.$op(x::AbstractVectorSymbolic{Tx}, y::AbstractVectorSymbolic{Ty}) where {Tx<:Number,Ty<:Number}
        s = union(x.params, y.params)
        Tv = promote_type(Tx,Ty)
        SymbolicVectorExpression{promote_type(Tx,Ty)}(s, $(Expr(:quote, op)), [x, y], x.dim, [])
    end

    @eval function Base.$op(x::Tx, y::AbstractVectorSymbolic{Ty}) where {Tx<:Number,Ty<:Number}
        Tv = promote_type(Tx,Ty)
        SymbolicVectorExpression{promote_type(Tx,Ty)}(y.params, $(Expr(:quote, op)), [x, y], x.dim, [])
    end

    @eval function Base.$op(x::AbstractVectorSymbolic{Tx}, y::Ty) where {Tx<:Number,Ty<:Number}
        Tv = promote_type(Tx,Ty)
        SymbolicVectorExpression{promote_type(Tx,Ty)}(x.params, $(Expr(:quote, op)), [x, y], x.dim, [])
    end
end


text:
abstract type HVRuleModel <: MixingRule end

struct HVRule{γ} <: HVRuleModel
    components::Array{String,1}
    activity::γ
end

@registermodel HVRule
export HVRule
function HVRule(components::Vector{String}; activity = Wilson, userlocations::Vector{String}=String[],activity_userlocations::Vector{String}=String[], verbose::Bool=false)
    init_activity = activity(components;userlocations = activity_userlocations,verbose)

    model = HVRule(components, init_activity)
    return model
end

function mixing_rule(model::RKModel,V,T,z,mixing_model::HVRuleModel,α,a,b,c)
    n = sum(z)
    invn2 = (one(n)/n)^2
    b̄ = dot(z,Symmetric(b),z) * invn2
    c̄ = dot(z,c)/n
    ā = b̄*((sum(z[i]*a[i,i]*α[i]/b[i,i] for i ∈ @comps))/n-excess_gibbs_free_energy(mixing_model.activity,1e5,T,z)/n/log(2))
    return ā,b̄,c̄
end

function mixing_rule(model::PRModel,V,T,z,mixing_model::HVRuleModel,α,a,b,c)
    λ = 1/(2*√(2))*log((2+√(2))/(2-√(2)))
    n = sum(z)
    invn2 = (one(n)/n)^2
    b̄ = dot(z,Symmetric(b),z) * invn2
    c̄ = dot(z,c)/n
    ā = b̄*((sum(z[i]*a[i,i]*α[i]/b[i,i] for i ∈ @comps))/n-excess_gibbs_free_energy(mixing_model.activity,1e5,T,z)/n/λ)
    return ā,b̄,c̄
end


text:
export Standardize

export standardize, standardize!, transform

# various transforms

### Standardization

immutable Standardize
    dim::Int
    mean::Vector{Float64}
    scale::Vector{Float64}

    function Standardize(d::Int, m::Vector{Float64}, s::Vector{Float64})
        Base.depwarn("Standardize is deprecated. Please use mean_and_std and zscore in StatsBase instead.", :Standardize)
        lenm = length(m)
        lens = length(s)
        lenm == d || lenm == 0 || throw(DimensionMismatch("Inconsistent dimensions."))
        lens == d || lens == 0 || throw(DimensionMismatch("Inconsistent dimensions."))
        new(d, m, s)
    end
end

indim(t::Standardize) = t.dim
outdim(t::Standardize) = t.dim

function transform!{YT<:Real,XT<:Real}(y::DenseArray{YT,1}, t::Standardize, x::DenseArray{XT,1})
    d = t.dim
    length(x) == length(y) == d || throw(DimensionMismatch("Inconsistent dimensions."))

    m = t.mean
    s = t.scale

    if isempty(m)
        if isempty(s)
            if !is(x, y)
                copy!(y, x)
            end
        else
            for i = 1:d
                @inbounds y[i] = x[i] * s[i]
            end
        end
    else
        if isempty(s)
            for i = 1:d
                @inbounds y[i] = x[i] - m[i]
            end
        else
            for i = 1:d
                @inbounds y[i] = s[i] * (x[i] - m[i])
            end
        end
    end
    return y
end

function transform!{YT<:Real,XT<:Real}(y::DenseArray{YT,2}, t::Standardize, x::DenseArray{XT,2})
    d = t.dim
    size(x,1) == size(y,1) == d || throw(DimensionMismatch("Inconsistent dimensions."))
    n = size(x,2)
    size(y,2) == n || throw(DimensionMismatch("Inconsistent dimensions."))

    m = t.mean
    s = t.scale

    if isempty(m)
        if isempty(s)
            if !is(x, y)
                copy!(y, x)
            end
        else
            for j = 1:n
                xj = view(x, :, j)
                yj = view(y, :, j)
                for i = 1:d
                    @inbounds yj[i] = xj[i] * s[i]
                end
            end
        end
    else
        if isempty(s)
            for j = 1:n
                xj = view(x, :, j)
                yj = view(y, :, j)
                for i = 1:d
                    @inbounds yj[i] = xj[i] - m[i]
                end
            end
        else
            for j = 1:n
                xj = view(x, :, j)
                yj = view(y, :, j)
                for i = 1:d
                    @inbounds yj[i] = s[i] * (xj[i] - m[i])
                end
            end
        end
    end
    return y
end

transform!{T<:AbstractFloat}(t::Standardize, x::DenseArray{T,1}) = transform!(x, t, x)
transform!{T<:AbstractFloat}(t::Standardize, x::DenseArray{T,2}) = transform!(x, t, x)

transform{T<:Real}(t::Standardize, x::DenseArray{T,1}) = transform!(Array(Float64, size(x)), t, x)
transform{T<:Real}(t::Standardize, x::DenseArray{T,2}) = transform!(Array(Float64, size(x)), t, x)

# estimate a standardize transform

function estimate{T<:Real}(::Type{Standardize}, X::DenseArray{T,2}; center::Bool=true, scale::Bool=true)
    d, n = size(X)
    n >= 2 || error("X must contain at least two columns.")

    m = Array(Float64, ifelse(center, d, 0))
    s = Array(Float64, ifelse(scale, d, 0))

    if center
        fill!(m, 0.0)
        for j = 1:n
            xj = view(X, :, j)
            for i = 1:d
                @inbounds m[i] += xj[i]
            end
        end
        scale!(m, 1.0 / n)
    end

    if scale
        fill!(s, 0.0)
        if center
            for j = 1:n
                xj = view(X, :, j)
                for i = 1:d
                    @inbounds s[i] += abs2(xj[i] - m[i])
                end
            end
        else
            for j = 1:n
                xj = view(X, :, j)
                for i = 1:d
                    @inbounds s[i] += abs2(xj[i])
                end
            end
        end
        for i = 1:d
            @inbounds s[i] = sqrt((n - 1) / s[i])
        end
    end

    return Standardize(d, m, s)
end

# standardize

function standardize{T<:Real}(X::DenseArray{T,2}; center::Bool=true, scale::Bool=true)
    t = estimate(Standardize, X; center=center, scale=scale)
    Y = transform(t, X)
    return (Y, t)
end

function standardize!{T<:AbstractFloat}(X::DenseArray{T,2}; center::Bool=true, scale::Bool=true)
    t = estimate(Standardize, X; center=center, scale=scale)
    Y = transform!(t, X)
    return (Y, t)
end


text:
using Test, Statistics, LinearAlgebra, Random
# make stable linear dynamics
Random.seed!(0)
eye = DifferentialDynamicProgramming.eye
costs = map(1:10) do MCiteration
    h    = .01  # time step
    n    = 10   # state dimension
    m    = 2    # control dimension
    A    = randn(n,n)
    A    = A-A' # skew-symmetric = pure imaginary eigenvalues
    A    = exp(h*A)        # discrete time
    B    = h*randn(n,m)

    # quadratic costs
    Q    = h*eye(n)
    R    = .1*h*eye(m)

    # control limits
    lims = []# ones(m,1)*[-1 1]*.6

    T    = 1000             # horizon
    x0   = ones(n,1)        # initial state
    u0   = .1*randn(m,T)    # initial controls

    # optimization problem
    N    = T+1
    fx   = A
    fu   = B
    cxx  = Q
    cxu  = zeros(size(B))
    cuu  = R

    # Specify dynamics functions
    function lin_dyn_df(x,u,Q,R)
        u[isnan.(u)] .= 0
        cx  = Q*x
        cu  = R*u
        fxx=fxu=fuu = []
        return fx,fu,fxx,fxu,fuu,cx,cu,cxx,cxu,cuu
    end
    function lin_dyn_f(x,u,A,B)
        u[isnan.(u)] .= 0
        xnew = A*x + B*u
        return xnew
    end

    function lin_dyn_cost(x,u,Q)
        c = 0.5*sum(x.*(Q*x)) + 0.5*sum(u.*(R*u))
        return c
    end

    f(x,u,i)     = lin_dyn_f(x,u,A,B)
    costfun(x,u) = lin_dyn_cost(x,u,Q)
    df(x,u)      = lin_dyn_df(x,u,Q,R)
    # plotFn(x)  = plot(squeeze(x,2)')

    # run the optimization
    @time x, u, L, Vx, Vxx, cost, otrace = iLQG(f,costfun,df, x0, u0, lims=lims, verbosity=3);
    # using Plots
    # plot(x', title="States", subplot=1, layout=(3,1), show=true)
    # plot!(u', title="Control signals", subplot=2, show=true)
    # plot!(cost, title="Cost", subplot=3, show=true)

    sum(cost)
end

@test maximum(costs) < 25 # This should be the case most of the times
@test mean(costs) < 10 # This should be the case most of the times
@test minimum(costs) < 5 # This should be the case most of the times


text:
using ParticleFilters
using Test
using Statistics

p = [8, 10, 12]
# particles in 2D state space
p2d = [
    [8, 8],
    [10, 10],
    [12, 12],
]
w = [4.0, 1.0, 1.0]

@testset "beliefs" begin

    @testset "unweighted" begin
        b = ParticleCollection(p)
        @inferred mean(b)
        @inferred cov(b)
        @test mean(b) ≈ 10
        @test cov(b) ≈ 8/3
        @test var(b) ≈ 8/3

        b = ParticleCollection(p2d)
        @inferred mean(b)
        @inferred cov(b)
        @test mean(b) ≈ [10, 10]
        @test cov(b) ≈ fill(8/3, 2, 2)
        @test var(b) ≈ fill(8/3, 2)
    end

    @testset "weighted" begin
        true_cov = cov([8, 8, 8, 8, 10, 12], corrected=false)
        b = WeightedParticleBelief(p, w)
        @inferred mean(b)
        @inferred cov(b)
        @test mean(b) ≈ 9
        @test cov(b) ≈ true_cov
        @test var(b) ≈ true_cov

        b = WeightedParticleBelief(p2d, w)
        @inferred mean(b)
        @inferred cov(b)
        @test mean(b) ≈ [9, 9]
        @test cov(b) ≈ fill(true_cov, 2, 2)
        @test var(b) ≈ fill(true_cov, 2)
    end
end


text:
import ..Continuations: step!, new_sweep!
using ..Continuations: AbstractContinuationSolver, ContinuationOptions,
    SweepSetup

abstract type BifurcationCache{PC} <: AbstractContinuationCache{PC} end

struct BifurcationSolver{
        R <: ContinuationSolver,
        P <: BifurcationProblem,
        C <: BifurcationCache,
        S <: BifurcationSolution,
        } <: AbstractContinuationSolver
    super::R
    prob::P
    opts::ContinuationOptions
    cache::C
    sol::S
end

function BifurcationSolver(prob::BifurcationProblem,
                           opts::ContinuationOptions)
    super = ContinuationSolver(prob, opts)
    cache = BifurcationCache(prob, super.cache)
    sol = BifurcationSolution(super.sol, sweeptype(prob, super))
    return BifurcationSolver(super, prob, opts, cache, sol)
end

TimeKind(::Type{<: BifurcationSolver{_, P}}) where {_, P} = TimeKind(P)

function step!(solver::BifurcationSolver)
    step!(solver.super)
    # calling [[../continuations/solver.jl::step!]]

    analyze!(solver.cache, solver.opts)
    record!(solver.sol, solver.cache)
    check_sweep_length(solver.sol.sweeps[end])
    post_step!(solver)  # TODO: remove
end

# This exists only for quick-and-dirty Bogdanov-Takens point handling
# in codim2 solver.  See: [[../codim2/solver.jl::post_step!]]
function post_step!(::Any) end

function new_sweep!(solver::BifurcationSolver, setup::SweepSetup)
    new_sweep!(solver.super, setup)
    # calling [[../continuations/solver.jl::new_sweep!]]

    allocate_sweep!(solver.sol, as(solver, ContinuationSolver))

    for u in setup.past_points
        re_analyze!(solver, u)
    end
    re_analyze!(solver, setup.u0)
    check_sweep_length(solver.sol.sweeps[end])
end
# TODO: `new_sweep!(solver.super, setup)` sets up `solver.super.cache`
# but `re_analyze!(solver, u)` rewrites the cache.  It causes no
# problem at the moment since the last `re_analyze!(solver, setup.u0)`
# set the cache back to the first state.

function allocate_sweep!(sol::BifurcationSolution, solver::ContinuationSolver)
    super = as(sol, ContinuationSolution)
    sweep = BifurcationSweep(super.sweeps[end], solver)
    push!(sol.sweeps, sweep)
end

function push_special_interval!(sweep::BifurcationSweep,
                                cache::BifurcationCache)
    push_special_interval!(sweep,
                           cache.point_type,
                           as(cache, ContinuationCache).J)
end

function push_special_interval!(sweep::BifurcationSweep,
                                point_type,
                                J1)
    super = as(sweep, ContinuationSweep)
    point = SpecialPointInterval(
        timekind(sweep),
        point_type,
        length(sweep) - 1,  # point_index
        super.u[end - 1],   # u0
        super.u[end],       # u1
        J1,
        WeakRef(sweep),
    )
    push!(sweep.special_intervals, point)
end
# `push_special_interval!` is relying on that it is called **after**
# the new point (`u1` in above) is `push!`ed to `super.u`.  With this
# assumption, `point_index` is same as the index used for
# `SimpleBifurcationInterval` in `push_point!` (see:
# [[../continuations/solution.jl::push_point!]]).  Note that `super.u`
# is updated via `push_point!` via `step!(::ContinuationSolver)`
# before `push_special_interval!` is called via `record!`.

"""
    is_special(point_type) :: Bool
"""
is_special(point_type::T) where T = point_type != regular_point(T)

function record!(sol::BifurcationSolution, cache::BifurcationCache)
    super = as(cache, ContinuationCache)
    sweep = sol.sweeps[end]
    push!(sweep.jacobians, copy(super.J))
    push!(sweep.eigvals, copy(cache.eigvals))

    if is_special(cache.point_type)
        push_special_interval!(sweep, cache)
    end
end


text:
@doc raw"""
    NelderMead(M, F [, p])
perform a nelder mead minimization problem for the cost funciton `F` on the
manifold `M`. If the initial population `p` is not given, a random set of
points is chosen.

This algorithm is adapted from the Euclidean Nelder-Mead method, see
[https://en.wikipedia.org/wiki/Nelder–Mead_method](https://en.wikipedia.org/wiki/Nelder–Mead_method)
and
[http://www.optimization-online.org/DB_FILE/2007/08/1742.pdf](http://www.optimization-online.org/DB_FILE/2007/08/1742.pdf).

# Input

* `M` – a manifold $\mathcal M$
* `F` – a cost function $F\colon\mathcal M\to\mathbb R$ to minimize
* `population` – (n+1 `random_point(M)`) an initial population of $n+1$ points, where $n$
  is the dimension of the manifold `M`.

# Optional

* `stopping_criterion` – ([`StopAfterIteration`](@ref)`(2000)`) a [`StoppingCriterion`](@ref)
* `retraction` – (`exp`) a `retraction(M,x,ξ)` to use.
* `α` – (`1.`) reflection parameter ($\alpha > 0$)
* `γ` – (`2.`) expansion parameter ($\gamma$)
* `ρ` – (`1/2`) contraction parameter, $0 < \rho \leq \frac{1}{2}$,
* `σ` – (`1/2`) shrink coefficient, $0 < \sigma \leq 1$

and the ones that are passed to [`decorate_options`](@ref) for decorators.

# Output
* either `x` the last iterate or the complete options depending on the optional
  keyword `return_options`, which is false by default (hence then only `x` is
  returned).
"""
function NelderMead(M::MT,
    F::Function,
    population = [random_point(M) for i=1:(manifoldDimension(M)+1) ];
    stopping_criterion::StoppingCriterion = StopAfterIteration(200000),
    α = 1., γ = 2., ρ=1/2, σ = 1/2,
    return_options=false,
    kwargs... #collect rest
  ) where {MT <: Manifold,T}
    p = CostProblem(M,F)
    o = NelderMeadOptions(population, stopping_criterion;
    α = α, γ = γ, ρ = ρ, σ = σ)
    o = decorate_options(o; kwargs...)
    resultO = solve(p,o)
    if return_options
        return resultO
    else
        return get_solver_result(resultO)
    end
end
#
# Solver functions
#
function initialize_solver!(p::P,o::O) where {P <: CostProblem, O <: NelderMeadOptions}
    # init cost and x
    o.costs = get_cost.(Ref(p), o.population )
    o.x = o.population[argmin(o.costs)] # select min
end
function step_solver!(p::P,o::O,iter) where {P <: CostProblem, O <: NelderMeadOptions}
    m = mean(p.M, o.population)
    ind = sortperm(o.costs) # reordering for cost and p, i.e. minimizer is at ind[1]
    ξ =log( p.M, m, o.population[last(ind)])
    # reflect last
    xr = exp(p.M, m, - o.α * ξ )
    Costr = get_cost(p,xr)
    # is it better than the worst but not better than the best?
    if Costr >= o.costs[first(ind)] && Costr < o.costs[last(ind)]
        # store as last
        o.population[last(ind)] = xr
        o.costs[last(ind)] = Costr
    end
    # --- Expansion ---
    if Costr < o.costs[first(ind)] # reflected is better than fist -> expand
        xe = exp(p.M, m, - o.γ * o.α * ξ)
        Coste = get_cost(p,xe)
        if Coste < Costr # expanded successful
            o.population[last(ind)] = xe
            o.costs[last(ind)] = Coste
        else # expansion failed but xr is still quite good -> store
            o.population[last(ind)] = xr
            o.costs[last(ind)] = Costr
        end
    end
    # --- Contraction ---
    if Costr > o.costs[ind[end-1]] # even worse than second worst
        if Costr < o.costs[last(ind)] # but at least better tham last
            # outside contraction
            xc = exp(p.M, m, - o.ρ*ξ)
            Costc = get_cost(p,xc)
            if Costc < Costr # better than reflected -> store as last
                o.population[last(ind)] = xr
                o.costs[last(ind)] = Costr
            end
        else # even worse than last -> inside contraction
            # outside contraction
            xc = exp(p.M, m, o.ρ*ξ)
            Costc = get_cost(p,xc)
            if Costc < o.costs[last(ind)] # better than last ? -> store
                o.population[last(ind)] = xr
                o.costs[last(ind)] = Costr
            end
        end
    end
    # --- Shrink ---
    for i=2:length(ind)
        o.population[ind[i]] = shortest_geodesic(p.M, o.population[ind[1]], o.population[ind[i]], o.σ)
        # update cost
        o.costs[ind[i]] = get_cost(p, o.population[ind[i]])
    end
    # store best
    o.x = o.population[ argmin(o.costs) ]
end
get_solver_result(p::P,o::O) where {P <: CostProblem, O <: NelderMeadOptions} = o.x


text:
# Hardy-Weinberg: A1A1: (1-p2)^2, A1A2: 2p2(1-p2), A2A2: p2^2
# ADDITIVE_MODEL=1 (A1A1: 0, A1A2: 1, A2A2: 2)
@inline function convert(::Type{T}, a2freq::AbstractFloat, ::Val{1}) where T <: AbstractFloat
    T(2a2freq)
end
# DOMINANT_MODEL=2 (A1A1: 0, A1A2: 1, A2A2: 1)
@inline function convert(::Type{T}, a2freq::AbstractFloat, ::Val{2}) where T <: AbstractFloat
    T(a2freq * (2 - a2freq))
end
# RECESSIVE_MODEL=3 (A1A1: 0, A1A2: 0, A2A2: 1)
@inline function convert(::Type{T}, a2freq::AbstractFloat, ::Val{3}) where T <: AbstractFloat
    T(abs2(a2freq))
end

"""
    Base.copyto!(v, s, P, Q; model=ADDITIVE_MODEL, center=false, scale=false)

Copy SnpArray `s` to a numeric vector or matrix `v`. Impute missing genotypes
according to ADMIXTURE estimates `P` (population A2 allele frequencies of SNPs)
and `Q` (population fractions of individuals). Columns of `P` should match columns of
`s`. Columns of `Q` should match rows of `s`. Note if the inferred minor allele
frequency from `P` and `Q` is ≤0.01, that genotype will be imputed according to
the inferred minor allele frequency.

# Positional arguments
- `v`: output vector or array.
- `s`: SnpArray.
- `P`: `K x S` array of A2 popluation allele frequencies.
- `Q`: `K x N` array of individual population fractions.

# Keyword arguments
- `model::Union{Val{1}, Val{2}, Val{3}}=ADDITIVE_MODEL`: `ADDITIVE_MODEL` (default), `DOMINANT_MODEL`, or `RECESSIVE_MODEL`.
- `center::Bool=false`: center each genotype dosage by its mean.
- `scale::Bool=false`: scale each genotype dosage by its standard deviation.
"""
function Base.copyto!(
    v      :: AbstractVecOrMat{T1},
    s      :: AbstractSnpArray,
    P      :: AbstractMatrix{T2},
    Q      :: AbstractMatrix{T2};
    model  :: Union{Val{1}, Val{2}, Val{3}} = ADDITIVE_MODEL,
    center :: Bool = false,
    scale  :: Bool = false
    ) where {T1 <: AbstractFloat, T2 <: AbstractFloat}
    m, n, K = size(s, 1), size(s, 2), size(P, 1)
    # center, scale, impute
    @inbounds for j in 1:n, i in 1:m
        vij = SnpArrays.convert(T1, s[i, j], model)
        a2freq = T2(0)
        for k in 1:K
            a2freq += Q[k, i] * P[k, j]
        end
        μij = SnpArrays.convert(T1, a2freq, model)
        σij = model == ADDITIVE_MODEL ? sqrt(μij * (1 - μij / 2)) : sqrt(μij * (1 - μij))
        v[i, j] = (isnan(vij) || a2freq ≤ 0.01 || a2freq ≥ 0.99) ? μij : vij
        center && (v[i, j] -= μij)
        scale  && (v[i, j] /= σij)
    end
    v
end

function Base.convert(
    ::Type{T},
    s::AbstractSnpArray,
    P::AbstractMatrix,
    Q::AbstractMatrix;
    kwargs...) where T <: Array
    T(s, P, Q; kwargs...)
end
Array{T,N}(s::AbstractSnpArray, P::AbstractMatrix, Q::AbstractMatrix; kwargs...) where {T,N} =
copyto!(Array{T,N}(undef, size(s)), s, P, Q; kwargs...)

function grm_admixture(
    s::AbstractSnpArray,
    P::AbstractMatrix,
    Q::AbstractMatrix,
    ::Type{T} = Float64
    ) where T <: AbstractFloat
    m, n = size(s)
    # convert genotype
    tic = time()
    G = Mmap.mmap(Matrix{T}, m, n) # slightly faster than G = Matrix{T}(undef, m, n)
    Base.copyto!(G, s, P, Q, model=ADDITIVE_MODEL, center=true, scale=true)
    @printf("convert genotype: %.2f seconds\n", time() - tic)
    # GG'
    tic = time()
    Φ = G * transpose(G) # m x m
    @printf("Φ = GG': %.2f seconds\n", time() - tic)
    # convert G to {0,1} matrix
    tic = time()
    @inbounds for (idx, g) in enumerate(G)
        G[idx] = ifelse(iszero(g), T(0), T(1))
        # iszero(g) || (G[idx] = T(1))
    end
    @printf("convert G to {0,1} matrix: %.2f seconds\n", time() - tic)
    # Sij = # SNPs observed in both individuals i and j
    tic = time()
    S = G * transpose(G)
    @printf("S = GG': %.2f seconds\n", time() - tic)
    # Φ /= 2S
    @inbounds for j in 1:m, i in 1:j
        Φ[i, j] /= 2S[i, j]
    end
    copytri!(Φ, 'U')
end # function grm_admixture


text:
abstract type BtgOptions end

"""
Options in btg model
TODOs:
    should have a default/suggested rangeθ and rangeλ for users
    could add MLE/bayesian options here
"""
mutable struct Options<:BtgOptions
    # model options
    transform_type::String
    kernel_type::String
    quadrature_type::Dict{String,String} # specify the type of discretization for different BTG model parameters
    quadrature_size::Dict{String,Int} # specify the grid size for different types of discretization
    # model parameter options
    parameter_range::Dict{String,Array{T, 2}} where T<:Real
    parameter_prior::Dict{String,priorType}
    parameter_dimθ::Int # dimension of theta, 1 or d, i.e. single lengthscale or dimension of input

    # prediction options
    confidence_level::T where T<:Real

    function Options(rangeθ::Array{T, 2}, rangeλ::Array{T, 2}, trainingData::AbstractTrainingData, parameter_dimθ::Int = 1) where T<:Real
        @assert size(rangeλ) == (1, 2) "Should provide 1-dimensional and two-end range of θ"
        @assert size(rangeθ)[2] == 2 "Should provide two-end range of θ"
        if parameter_dimθ > 1
            d = getDimension(trainingData)
            @assert d == parameter_dimθ "Multi-dimensional parameter θ should match dimension of input data"
            # @assert d == parameter_dimθ
        end
        if parameter_dimθ > 1 && size(rangeθ)[1] == 1
            rangeθs = repeat(rangeθ, parameter_dimθ, 1) # adjust rangeθs to match parameter_dimθ
        elseif parameter_dimθ == 1 && size(rangeθ)[1] > 1
            @warn "Given multiple lengthscale range, BTG adjusting parameter_dimθ to match."
            parameter_dimθ = size(rangeθ)[1]
            rangeθs = rangeθ
        else
            rangeθs = rangeθ
        end
        transform_type = "BoxCox" # else: "IdentityTransform", "ShiftedBoxCox"
        kernel_type = "Gaussian"
        quadrature_type = Dict("θ" => "Gaussian", "λ" => "Gaussian") # recommended combo: ["SparseGrid", "GaussianQuadrature"], ["SparseCarlo", "SparseCarlo"], ["QuasiMonteCarlo", "QuasiMonteCarlo"]
        quadrature_size = Dict("Gaussian" => 12, "MonteCarlo" => 400) # for "SparseCarlo" and "SparseGrid", depends on input dimension.
        parameter_range = Dict("θ" => rangeθs, "λ" => rangeλ)
        parameter_prior = Dict("θ" => inverseUniform(rangeθ), "λ" => Uniform(rangeλ))
        confidence_level = .95
        new(transform_type, kernel_type, quadrature_type, quadrature_size, parameter_range, parameter_prior, parameter_dimθ, confidence_level)
    end
end

function print(O::BtgOptions)
    println("\n\n=============== BTG OPTIONS =============== ")
    println("   transform type:       $(O.transform_type)")
    println("   Kernel type:          $(O.kernel_type)")
    println("   quadrature type:      $(O.quadrature_type)")
    println("   quadrature size:      $(O.quadrature_size)")
    println("   Confidence level:     $(O.confidence_level)")
    println("   Parameter range:      ")
    for key in keys(O.parameter_range)
        println("                         $(key): $(O.parameter_range[key])")
    end
    println("   Parameter prior:      ")
    for key in keys(O.parameter_prior)
        println("                         $(key): $(O.parameter_prior[key])")
    end
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using Joy
using Base.Test

orig = readlines("file")

@testset "Numerical argument parser" begin
    @test Joy.parse_n(['a','b', 'c', '2', '3']) == 23
    @test Joy.parse_n(['a','1', 'b', 'c', '3']) == 3
    @test Joy.parse_n(['a','1', 'b', '1', 'c']) == 1
    @test Joy.parse_n(['a']) == 1
    @test Joy.parse_n([]) == 1
end

@testset "chr2ind range" begin
    s = "aåäöasdf"
    @test Joy.chr2ind(s, 2:4) == 2:1:7
end

@testset "Find symbol" begin
    @test Joy.findsymbol("Lörem ipsum dolor sit amet", 'i', 1) == 7
    @test Joy.findsymbol("Lörem ipsum dolor sit amet", 'i', 2) == 20
end

buffer = Joy.self
Joy.open(buffer, "file")

@testset "Movements, keybindings, delete" begin
    Joy.handle_raw(buffer, 'j')
    @test Joy.line(buffer) == orig[2]
    Joy.replay(buffer, "2j")
    @test Joy.line(buffer) == orig[4]
    Joy.replay(buffer, "dd")
    @test Joy.line(buffer) == orig[5]
    Joy.replay(buffer, "gg")
    Joy.replay(buffer, "i\r\r\e")
    Joy.replay(buffer, "kia\e")
    @test Joy.line(buffer) == "a"
    cmd = ":settext(self, map(z->\"\$z \"*line(self, z), 1:height(self)))\r"
    Joy.replay(buffer, cmd)
    Joy.replay(buffer, "12gg")
    @test Joy.line(buffer)[1:2] == "12"
end

@testset "Word handling" begin
    @test Joy.nextword_pos("Lörem ipsum dolor sit amet", 2) == 7
    @test Joy.nextword_pos(" Lörem ipsum dolor sit amet") == 2
    @test Joy.nextword_pos("Lörem ipsum dolor sit amet", 4) == 19
    @test Joy.nextword_pos("(", 2) == -1
    Joy.open(buffer, "file")
    Joy.replay(buffer, "99k99h")
    Joy.replay(buffer, "wwww")
    @test Joy.line(buffer)[Joy.x(buffer):Joy.x(buffer)+2] == "sit"
    Joy.replay(buffer, "50hW")
    @test Joy.line(buffer)[Joy.x(buffer):Joy.x(buffer)+4] == "ipsum"
    Joy.replay(buffer, "b")
    @test Joy.line(buffer)[Joy.x(buffer):Joy.x(buffer)+4] == "Lorem"
    Joy.replay(buffer, "wj10w")
    @test Joy.line(buffer) == ""
    Joy.replay(buffer, "e")
    @test Joy.line(buffer)[Joy.x(buffer)] == 't'
    Joy.replay(buffer, "3b")
    @test Joy.line(buffer)[Joy.x(buffer)] == '.'
end

@testset "Delete mode" begin
    Joy.replay(buffer, "99k99h")
    Joy.replay(buffer, "dw")
    @test buffer.text[1][1] == orig[1][2]
    Joy.replay(buffer, "2dw")
    @test buffer.text[1][1] == orig[1][14]
    # Not perfect, should override this
    Joy.replay(buffer, "dfc")
    @test buffer.text[1][1] == 'c'
end

@testset "Yank/Paste" begin
    Joy.replay(buffer, "99k99h")
    Joy.settext(buffer, ["c", "c"])
    Joy.paste(buffer, [1,1], "a\nb")
    @test buffer.text == ["ca", "b", "c"]
    Joy.settext(buffer, ["Lorem", "ipsum"])
    Joy.replay(buffer, "99k99h")
    Joy.replay(buffer, "\"ayw")
    @test buffer.registers['a'] == "Lorem\n"
    Joy.settext(buffer, ["a", "b"])
    Joy.replay(buffer, "99k99h")
    Joy.replay(buffer, "yyp")
    @test buffer.text == ["a", "a", "b"]
    Joy.settext(buffer, ["a", "b"])
    Joy.replay(buffer, "99k99h")
    Joy.replay(buffer, "ddp")
    Joy.settext(buffer, ["b", "a"])
end

@testset "Search" begin
    Joy.replay(buffer, "99k99h")
    Joy.settext(buffer, ["Lorem ipsum", "Lorem2 ipsum2"])
    Joy.replay(buffer, "/ipsum\r")
    @test Joy.pos(buffer) == [1,7]
    Joy.replay(buffer, "/ipsum2\r")
    @test Joy.pos(buffer) == [2,8]
    Joy.replay(buffer, "99k99h")
    Joy.replay(buffer, "d/Lorem2\r")
    @test buffer.text == ["Lorem2 ipsum2"]
    Joy.replay(buffer, "y/ipsum2\r")
    @test buffer.registers['"'] == "Lorem2 "
    Joy.settext(buffer, ["Lorem ipsum", "Lorem2 ipsum2"])
    Joy.replay(buffer, "99k99h")
    Joy.replay(buffer, "/ipsum\r")
    Joy.replay(buffer, "n")
    @test Joy.pos(buffer) == [2,8]
end

@testset "Undo" begin
    Joy.settext(buffer, ["Lorem ipsum", "ipsum lorem"])
    Joy.replay(buffer, "99k99h")
    Joy.replay(buffer, "ddu")
    @test buffer.text == ["Lorem ipsum", "ipsum lorem"]
    Joy.replay(buffer, "\x12")
    @test buffer.text == ["ipsum lorem"]
end


text:
using CSV
using DelimitedFiles
using Plots
include("loadfunctions.jl")
function runmodelprecipitationzones_all_output(Potential_Evaporation::Array{Float64,1}, Precipitation_All_Zones::Array{Array{Float64,2},1}, Temperature_Elevation_Catchment::Array{Float64,2}, Inputs_All_Zones::Array{Array{HRU_Input,1},1}, Storages_All_Zones::Array{Array{Storages,1},1}, SlowStorage::Float64, parameters::Array{Parameters,1}, slow_parameters::Slow_Paramters, Area_Zones::Array{Float64,1}, Area_Zones_Percent::Array{Float64,1}, Elevation_Percentage::Array{Array{Float64,1},1}, Elevation_Zone_Catchment::Array{Float64,1}, ID_Prec_Zones::Array{Int64,1}, Nr_Elevationbands_All_Zones::Array{Int64,1}, Elevations_Each_Precipitation_Zone)
        Total_Discharge = zeros(length(Precipitation_All_Zones[1][:,1]))
        Total_Snowstorage = zeros(length(Precipitation_All_Zones[1][:,1]))
        Total_GWstorage = zeros(length(Precipitation_All_Zones[1][:,1]))
        #Total_Snow_Elevations = Array{Float64,2}[]
        count = zeros(length(Precipitation_All_Zones[1][:,1]), length(Elevation_Zone_Catchment))
        # Snow_Overall_Objective_Function = 0
        # Snow_Elevations_All = zeros(length(Precipitation_All_Zones[1][:,1]), length(Elevation_Zone_Catchment))
        # Snow_Extend_All = zeros(length(observed_snow_cover[1][:,1]), length(Elevation_Zone_Catchment))
        # Snow_Extend_All_Observed = zeros(length(observed_snow_cover[1][:,1]), length(Elevation_Zone_Catchment))
        # percentage = zeros(length(Precipitation_All_Zones[1][:,1]), length(Elevation_Zone_Catchment))
        # percentage_Snow_Extend = zeros(length(observed_snow_cover[1][:,1]), length(Elevation_Zone_Catchment))
        # percentage_soil = zeros(length(Precipitation_All_Zones[1][:,1]), 4)
        # Soilstorage_All = zeros(length(Precipitation_All_Zones[1][:,1]), 4)
        # Faststorage_All = zeros(length(Precipitation_All_Zones[1][:,1]), 4)
        for i in 1: length(ID_Prec_Zones)
                # take the storages and input of the specific precipitation zone
                Inputs_HRUs = Inputs_All_Zones[i]
                Storages_HRUs = Storages_All_Zones[i]
                # run the model for the specific precipitation zone
                Discharge, Snow_Extend, GWstorage, Snowstorage, Snow_Elevations, Soilstorage, Faststorage = runmodel_alloutput(Area_Zones[i], Potential_Evaporation, Precipitation_All_Zones[i], Temperature_Elevation_Catchment,
                        Inputs_HRUs[1], Inputs_HRUs[2], Inputs_HRUs[3], Inputs_HRUs[4],
                        Storages_HRUs[1], Storages_HRUs[2], Storages_HRUs[3], Storages_HRUs[4], SlowStorage,
                        parameters[1], parameters[2], parameters[3], parameters[4], slow_parameters, Nr_Elevationbands_All_Zones[i], Elevation_Percentage[i])
                # sum up the discharge of all precipitation zones
                Total_Discharge += Discharge
                Total_GWstorage += GWstorage * Area_Zones_Percent[i]
                Total_Snowstorage += Snowstorage * Area_Zones_Percent[i]
        end
        # calculate the mean difference over all precipitation zones
        # Snow_Elevations_All = Snow_Elevations_All ./ percentage
        # Snow_Extend_All = Snow_Extend_All ./percentage_Snow_Extend
        # Snow_Extend_All_Observed = Snow_Extend_All_Observed ./percentage_Snow_Extend
        # Soilstorage_All = Soilstorage_All ./ percentage_soil
        # Faststorage_All = Faststorage_All ./ percentage_soil
        #println("works too")
        return Total_Discharge::Array{Float64,1}, Total_GWstorage::Array{Float64,1}, Total_Snowstorage::Array{Float64,1}
end

function run_projections_feistritz(path_to_projection, path_to_best_parameter, startyear, endyear, period, spinup)
        local_path = "/home/sarah/"
        # ------------ CATCHMENT SPECIFIC INPUTS----------------
        ID_Prec_Zones = [109967]
        # size of the area of precipitation zones
        Area_Zones = [115496400.]
        Area_Catchment = sum(Area_Zones)
        Area_Zones_Percent = Area_Zones / Area_Catchment
        Snow_Threshold = 600
        Height_Threshold = 4000

        Mean_Elevation_Catchment = 900 # in reality 917
        # two last entries of array are height of temp measurement
        Elevations_Catchment = Elevations(200.0, 400.0, 1600.0, 488., 488.)
        Sunhours_Vienna = [8.83, 10.26, 11.95, 13.75, 15.28, 16.11, 15.75, 14.36, 12.63, 10.9, 9.28, 8.43]
        # where to skip to in data file of precipitation measurements
        Skipto = [24]
        # get the areal percentage of all elevation zones in the HRUs in the precipitation zones
        Areas_HRUs =  CSV.read(local_path*"HBVModel/Feistritz/HBV_Area_Elevation.csv", skipto=2, decimal='.', delim = ',')
        # get the percentage of each HRU of the precipitation zone
        Percentage_HRU = CSV.read(local_path*"HBVModel/Feistritz/HRU_Prec_Zones.csv", header=[1], decimal='.', delim = ',')
        Elevation_Catchment = convert(Vector, Areas_HRUs[2:end,1])

        # timeperiod for which model should be run (look if timeseries of data has same length)
        Timeseries = readdlm(path_to_projection*"pr_model_timeseries.txt")
        Timeseries = Date.(Timeseries, Dates.DateFormat("y,m,d"))
        if endyear <= Dates.year(Timeseries[end])
                startyear =  endyear - 29 - spinup
                indexstart_Proj = findfirst(x-> x == startyear, Dates.year.(Timeseries))[1]
                indexend_Proj = findlast(x-> x == endyear, Dates.year.(Timeseries))[1]
        else
                endyear = Dates.year(Timeseries[end])
                startyear = endyear - 29 - spinup # -3 for the spinup time
                indexend_Proj = length(Timeseries)
                indexstart_Proj = findfirst(x-> x == startyear, Dates.year.(Timeseries))[1]
                print(Timeseries[end])
        end
        println(startyear, " ", endyear, "\n")
        indexstart_Proj = findfirst(x-> x == startyear, Dates.year.(Timeseries))[1]
        indexend_Proj = findlast(x-> x == endyear, Dates.year.(Timeseries))[1]
        Timeseries = Timeseries[indexstart_Proj:indexend_Proj]

        #------------ TEMPERATURE AND POT. EVAPORATION CALCULATIONS ---------------------
        #Temperature is the same in whole catchment
        # Temperature Measurements are taken at Maria Luggau
        Projections_Temperature = readdlm(path_to_projection*"tas_10510_sim1.txt", ',')
        Temperature_Daily = Projections_Temperature[indexstart_Proj:indexend_Proj] ./ 10
        Temperature_Daily = Temperature_Daily[:,1]

        Elevation_Zone_Catchment, Temperature_Elevation_Catchment, Total_Elevationbands_Catchment = gettemperatureatelevation(Elevations_Catchment, Temperature_Daily)
        # get the temperature data at the mean elevation to calculate the mean potential evaporation
        Temperature_Mean_Elevation = Temperature_Elevation_Catchment[:,findfirst(x-> x==Mean_Elevation_Catchment, Elevation_Zone_Catchment)]
        Potential_Evaporation = getEpot_Daily_thornthwaite(Temperature_Mean_Elevation, Timeseries, Sunhours_Vienna)

        # ------------- LOAD PRECIPITATION DATA OF EACH PRECIPITATION ZONE ----------------------
        # get elevations at which precipitation was measured in each precipitation zone
        Elevations_109967= Elevations(200., 400., 1600., 563.,488.)
        Elevations_All_Zones = [Elevations_109967]

        #get the total discharge
        Total_Discharge = zeros(length(Temperature_Daily))
        Inputs_All_Zones = Array{HRU_Input, 1}[]
        Storages_All_Zones = Array{Storages, 1}[]
        Precipitation_All_Zones = Array{Float64, 2}[]
        Precipitation_Gradient = 0.0
        Elevation_Percentage = Array{Float64, 1}[]
        Nr_Elevationbands_All_Zones = Int64[]
        Elevations_Each_Precipitation_Zone = Array{Float64, 1}[]

        for i in 1: length(ID_Prec_Zones)
                Precipitation_Zone = readdlm(path_to_projection*"pr_"*string(ID_Prec_Zones[i])*"_sim1.txt", ',')
                Precipitation_Zone = Precipitation_Zone[indexstart_Proj:indexend_Proj] ./ 10
                Elevation_HRUs, Precipitation, Nr_Elevationbands = getprecipitationatelevation(Elevations_All_Zones[i], Precipitation_Gradient, Precipitation_Zone)
                push!(Precipitation_All_Zones, Precipitation)
                push!(Nr_Elevationbands_All_Zones, Nr_Elevationbands)
                push!(Elevations_Each_Precipitation_Zone, Elevation_HRUs)

                index_HRU = (findall(x -> x==ID_Prec_Zones[i], Areas_HRUs[1,2:end]))
                # for each precipitation zone get the relevant areal extentd
                Current_Areas_HRUs = convert(Matrix, Areas_HRUs[2: end, index_HRU])
                # the elevations of each HRU have to be known in order to get the right temperature data for each elevation
                Area_Bare_Elevations, Bare_Elevation_Count = getelevationsperHRU(Current_Areas_HRUs[:,1], Elevation_Catchment, Elevation_HRUs)
                Area_Forest_Elevations, Forest_Elevation_Count = getelevationsperHRU(Current_Areas_HRUs[:,2], Elevation_Catchment, Elevation_HRUs)
                Area_Grass_Elevations, Grass_Elevation_Count = getelevationsperHRU(Current_Areas_HRUs[:,3], Elevation_Catchment, Elevation_HRUs)

                Area_Rip_Elevations, Rip_Elevation_Count = getelevationsperHRU(Current_Areas_HRUs[:,4], Elevation_Catchment, Elevation_HRUs)
                #print(Bare_Elevation_Count, Forest_Elevation_Count, Grass_Elevation_Count, Rip_Elevation_Count)
                println((Area_Bare_Elevations), " ", Bare_Elevation_Count,"\n")
                println((Area_Forest_Elevations), " ", Forest_Elevation_Count,"\n")
                Area_Bare_Elevations = [0.0]
                Bare_Elevation_Count = [1]
                @assert 0.999 <= sum(Area_Bare_Elevations) <= 1.0001 || sum(Area_Bare_Elevations) == 0

                @assert 0.999 <= sum(Area_Forest_Elevations) <= 1.0001
                @assert 0.999 <= sum(Area_Grass_Elevations) <= 1.0001
                @assert 0.999 <= sum(Area_Rip_Elevations) <= 1.0001

                Area = Area_Zones[i]
                Current_Percentage_HRU = Percentage_HRU[:,1 + i]/Area
                # calculate percentage of elevations
                Perc_Elevation = zeros(Total_Elevationbands_Catchment)
                for j in 1 : Total_Elevationbands_Catchment
                        for h in 1:4
                                Perc_Elevation[j] += Current_Areas_HRUs[j,h] * Current_Percentage_HRU[h]
                        end
                end
                Perc_Elevation = Perc_Elevation[(findall(x -> x!= 0, Perc_Elevation))]
                @assert 0.99 <= sum(Perc_Elevation) <= 1.01
                push!(Elevation_Percentage, Perc_Elevation)
                # calculate the inputs once for every precipitation zone because they will stay the same during the Monte Carlo Sampling
                bare_input = HRU_Input(Area_Bare_Elevations, Current_Percentage_HRU[1],zeros(length(Bare_Elevation_Count)) , Bare_Elevation_Count, length(Bare_Elevation_Count), (Elevations_All_Zones[i].Min_elevation + 100, Elevations_All_Zones[i].Max_elevation - 100), (Snow_Threshold, Height_Threshold), 0, [0], 0, [0], 0, 0)
                forest_input = HRU_Input(Area_Forest_Elevations, Current_Percentage_HRU[2], zeros(length(Forest_Elevation_Count)) , Forest_Elevation_Count, length(Forest_Elevation_Count), (Elevations_All_Zones[i].Min_elevation + 100, Elevations_All_Zones[i].Max_elevation - 100), (Snow_Threshold, Height_Threshold), 0, [0], 0, [0],  0, 0)
                grass_input = HRU_Input(Area_Grass_Elevations, Current_Percentage_HRU[3], zeros(length(Grass_Elevation_Count)) , Grass_Elevation_Count,length(Grass_Elevation_Count), (Elevations_All_Zones[i].Min_elevation + 100, Elevations_All_Zones[i].Max_elevation - 100), (Snow_Threshold, Height_Threshold), 0, [0], 0, [0],  0, 0)
                rip_input = HRU_Input(Area_Rip_Elevations, Current_Percentage_HRU[4], zeros(length(Rip_Elevation_Count)) , Rip_Elevation_Count, length(Rip_Elevation_Count), (Elevations_All_Zones[i].Min_elevation + 100, Elevations_All_Zones[i].Max_elevation - 100),(Snow_Threshold, Height_Threshold), 0, [0], 0, [0],  0, 0)

                all_inputs = [bare_input, forest_input, grass_input, rip_input]
                #print(typeof(all_inputs))
                push!(Inputs_All_Zones, all_inputs)

                bare_storage = Storages(0, zeros(length(Bare_Elevation_Count)), zeros(length(Bare_Elevation_Count)), zeros(length(Bare_Elevation_Count)), 0)
                forest_storage = Storages(0, zeros(length(Forest_Elevation_Count)), zeros(length(Forest_Elevation_Count)), zeros(length(Bare_Elevation_Count)), 0)
                grass_storage = Storages(0, zeros(length(Grass_Elevation_Count)), zeros(length(Grass_Elevation_Count)), zeros(length(Bare_Elevation_Count)), 0)
                rip_storage = Storages(0, zeros(length(Rip_Elevation_Count)), zeros(length(Rip_Elevation_Count)), zeros(length(Bare_Elevation_Count)), 0)

                all_storages = [bare_storage, forest_storage, grass_storage, rip_storage]
                push!(Storages_All_Zones, all_storages)
        end
        # ---------------- CALCULATE OBSERVED OBJECTIVE FUNCTIONS -------------------------------------
        # calculate the sum of precipitation of all precipitation zones to calculate objective functions
        #Total_Precipitation = Precipitation_All_Zones[1][:,1]*Area_Zones_Percent[1] + Precipitation_All_Zones[2][:,1]*Area_Zones_Percent[2] + Precipitation_All_Zones[3][:,1]*Area_Zones_Percent[3]
        Total_Precipitation = Precipitation_All_Zones[1][:,1]
        index_spinup = findfirst(x -> Dates.year(x) == startyear + spinup, Timeseries)
        #print("index",index_spinup,"\n")
        # evaluations chouls alsways contain whole year
        index_lastdate = findlast(x -> Dates.year(x) == endyear, Timeseries)
        print("index",typeof(index_lastdate),typeof(index_spinup),"\n")
        Timeseries_Obj = Timeseries[index_spinup: end]
        # ---------------- START MONTE CARLO SAMPLING ------------------------
        GWStorage = 70.0
        All_Discharge = zeros(length(Timeseries_Obj))
        All_Snowstorage = zeros(length(Timeseries_Obj))
        All_Snowmelt = zeros(length(Timeseries_Obj))
        All_Snow_Cover = transpose(length(Elevation_Zone_Catchment))
        # get the parameter sets of the calibrations
        best_calibrations = readdlm(path_to_best_parameter, ',')
        parameters_best_calibrations = best_calibrations[:,10:29]
        for n in 1 : 1:size(parameters_best_calibrations)[1]
                Current_Inputs_All_Zones = deepcopy(Inputs_All_Zones)
                Current_Storages_All_Zones = deepcopy(Storages_All_Zones)
                Current_GWStorage = deepcopy(GWStorage)
                # use parameter sets of the calibration as input
                beta_Bare, beta_Forest, beta_Grass, beta_Rip, Ce, Interceptioncapacity_Forest, Interceptioncapacity_Grass, Interceptioncapacity_Rip, Kf_Rip, Kf, Ks, Meltfactor, Mm, Ratio_Pref, Ratio_Riparian, Soilstoaragecapacity_Bare, Soilstoaragecapacity_Forest, Soilstoaragecapacity_Grass, Soilstoaragecapacity_Rip, Temp_Thresh = parameters_best_calibrations[n, :]
                bare_parameters = Parameters(beta_Bare, Ce, 0, 0.0, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Bare, Temp_Thresh)
                forest_parameters = Parameters(beta_Forest, Ce, 0, Interceptioncapacity_Forest, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Forest, Temp_Thresh)
                grass_parameters = Parameters(beta_Grass, Ce, 0, Interceptioncapacity_Grass, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Grass, Temp_Thresh)
                rip_parameters = Parameters(beta_Rip, Ce, 0.0, Interceptioncapacity_Rip, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Rip, Temp_Thresh)
                slow_parameters = Slow_Paramters(Ks, Ratio_Riparian)

                parameters = [bare_parameters, forest_parameters, grass_parameters, rip_parameters]
                parameters_array = parameters_best_calibrations[n, :]

                # Discharge, Snow_Cover, Snow_Melt = runmodelprecipitationzones_future(Potential_Evaporation, Precipitation_All_Zones, Temperature_Elevation_Catchment, Current_Inputs_All_Zones, Current_Storages_All_Zones, Current_GWStorage, parameters, slow_parameters, Area_Zones, Area_Zones_Percent, Elevation_Percentage, Elevation_Zone_Catchment, ID_Prec_Zones, Nr_Elevationbands_All_Zones, Elevations_Each_Precipitation_Zone)
                # All_Discharge = hcat(All_Discharge, Discharge[index_spinup:index_lastdate])
                # #All_Snowstorage = hcat(All_Snowstorage, Snow_Storage[index_spinup:index_lastdate])
                # All_Snowmelt = hcat(All_Snowmelt, Snow_Melt[index_spinup:index_lastdate])
                # All_Snow_Cover = vcat(All_Snow_Cover, Snow_Cover[:,index_spinup:index_lastdate])
                Discharge, GWstorage, Snowstorage = runmodelprecipitationzones_all_output(Potential_Evaporation, Precipitation_All_Zones, Temperature_Elevation_Catchment, Current_Inputs_All_Zones, Current_Storages_All_Zones, Current_GWStorage, parameters, slow_parameters, Area_Zones, Area_Zones_Percent, Elevation_Percentage, Elevation_Zone_Catchment, ID_Prec_Zones, Nr_Elevationbands_All_Zones, Elevations_Each_Precipitation_Zone)
                All_Snowstorage = hcat(All_Snowstorage, Snowstorage[index_spinup: index_lastdate])
        end
        #All_Goodness = transpose(All_Goodness[:, 2:end])
        #All_Discharge = transpose(All_Discharge[:, 2:end])
        All_Snowstorage = transpose(All_Snowstorage[:,2:end])
        #All_Snowmelt = transpose(All_Snowmelt[:,2:end])
        #All_Snow_Cover = All_Snow_Cover[2:end,:]
        #save the results for the projections
        #writedlm(path_to_projection*"100_model_results_05_10.csv", All_Goodness, ',')
        #writedlm(path_to_projection*"300_model_results_discharge_"*period*".csv", All_Discharge, ',')
        writedlm(path_to_projection*"300_model_results_snow_storage_"*period*".csv", All_Snowstorage, ',')
        # writedlm(path_to_projection*"300_model_results_snow_melt_"*period*".csv", All_Snowmelt, ',')
        # writedlm(path_to_projection*"300_model_results_snow_cover_"*period*".csv", All_Snow_Cover, ',')
        # writedlm(path_to_projection*"results_epot_"*period*".csv", Potential_Evaporation[index_spinup: index_lastdate], ',')
        # writedlm(path_to_projection*"results_precipitation_"*period*".csv", Total_Precipitation[index_spinup: index_lastdate], ',')
        #return All_Discharge
end
path = "/home/sarah/Master/Thesis/Data/Projektionen/new_station_data_rcp45/rcp45/"
# 14 different projections
Name_Projections = readdir(path)
# run the model for all projections using the best 100 parameter sets
# for (i, name) in enumerate(Name_Projections)
#         print(i)
#         name = Name_Projections[i]
#         Discharge_present = run_projections_gailtal(path*name*"/Gailtal/", "/home/sarah/Master/Thesis/Calibrations/Gailtal_less_dates/Gailtal_Parameterfit_All_less_dates_best_proj.csv", 1978, 2010, "past_2010")
#         Discharge_future = run_projections_gailtal(path*name*"/Gailtal/", "/home/sarah/Master/Thesis/Calibrations/Gailtal_less_dates/Gailtal_Parameterfit_All_less_dates_best_proj.csv", 2068, 2100, "future_2100")
#         #print(size(Discharge_present))
# end

for (i, name) in enumerate(Name_Projections)
        print(i)
        name = Name_Projections[i]
        Discharge_present = run_projections_feistritz(path*name*"/Pitten/", "/home/sarah/Master/Thesis/Calibrations/Feistritz_less_dates/Feistritz_Parameterfit_All_less_dates_unique_best_300.csv", 1981, 2010, "past_2010",3)
        Discharge_future = run_projections_feistritz(path*name*"/Pitten/", "/home/sarah/Master/Thesis/Calibrations/Feistritz_less_dates/Feistritz_Parameterfit_All_less_dates_unique_best_300.csv", 2071, 2100, "future_2100", 10)
end

path = "/home/sarah/Master/Thesis/Data/Projektionen/new_station_data_rcp85/rcp85/"
# 14 different projections
Name_Projections = readdir(path)

for (i, name) in enumerate(Name_Projections)
        print(i)
        name = Name_Projections[i]
        Discharge_present = run_projections_feistritz(path*name*"/Pitten/", "/home/sarah/Master/Thesis/Calibrations/Feistritz_less_dates/Feistritz_Parameterfit_All_less_dates_unique_best_300.csv", 1981, 2010, "past_2010",3)
        Discharge_future = run_projections_feistritz(path*name*"/Pitten/", "/home/sarah/Master/Thesis/Calibrations/Feistritz_less_dates/Feistritz_Parameterfit_All_less_dates_unique_best_300.csv", 2071, 2100, "future_2100", 10)
end

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
function generators_range(line::Ref{UInt}, ns::Namespace, arg1::AetherlangObject, arg2::AetherlangObject)::AetherlangObject
    AetherlangObject(arg1.dataref:arg2.dataref, Ref{String}(AETH_BUILTIN_TYPES.collection))
end
function generators_range(line::Ref{UInt}, ns::Namespace, arg1::AetherlangObject, arg2::AetherlangObject, arg3::AetherlangObject)::AetherlangObject
    AetherlangObject(arg1.dataref:arg2.dataref:arg3.dataref, Ref{String}(AETH_BUILTIN_TYPES.collection))
end
function generators_gencoll(line::Ref{UInt}, ns::Namespace, exp::AetherlangObject, var::AetherlangObject, coll::AetherlangObject)::AetherlangObject
    v::Vector = []
    nns::Namespace = copy(ns)
    for i in 1:length(coll.dataref)
        nns[var.dataref] = AetherlangObject(coll.dataref[i])
        push!(v, aetherlang_eval!(exp.dataref, nns, Ref{UInt}(0)).dataref)
    end
    AetherlangObject(v)
end
function generators_gencoll(line::Ref{UInt}, ns::Namespace, exp::AetherlangObject, var::AetherlangObject, coll::AetherlangObject, cond::AetherlangObject)::AetherlangObject
    v::Vector = []
    nns::Namespace = copy(ns)
    for i in 1:length(coll.dataref)
        nns[var.dataref] = AetherlangObject(coll.dataref[i])
        if aetherlang_eval!(cond.dataref, nns, Ref{UInt}(0)).type[] !== AETH_BUILTIN_TYPES.void
            push!(v, aetherlang_eval!(exp.dataref, nns, Ref{UInt}(0)).dataref)
        end
    end
    AetherlangObject(v)
end
function generators_gencoll(line::Ref{UInt}, ns::Namespace, exp::AetherlangObject, var::AetherlangObject, coll::AetherlangObject, cond::AetherlangObject, otherwise_exp::AetherlangObject)::AetherlangObject
    v::Vector = []
    nns::Namespace = copy(ns)
    for i in 1:length(coll.dataref)
        nns[var.dataref] = AetherlangObject(coll.dataref[i])
        if aetherlang_eval!(cond.dataref, nns, Ref{UInt}(0)).type[] !== AETH_BUILTIN_TYPES.void
            push!(v, aetherlang_eval!(exp.dataref, nns, Ref{UInt}(0)).dataref)
        else
            push!(v, aetherlang_eval!(otherwise_exp.dataref, nns, Ref{UInt}(0)).dataref)
        end
    end
    AetherlangObject(v)
end

generators_namespace_modify = Namespace(
    "range" => AetherlangObject(generators_range),
    "gencoll" => AetherlangObject(generators_gencoll)
)
generators_namespace_modify["rng"] = generators_namespace_modify["range"]

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
using SeisNoise
using Test


text:
# ## [Пример №2, краткий, стандартные параметры](@id de_2)
#
# В этом примере:
# - получаем стандартные входные параметры,
# - решается прямая задача,
# - генерируются экспериментальные данные.
# Только в отличии от [Подробного примера №1](@ref de_1)
# будем пользоваться функциями-сокращениями.
#
# Этот файл будет использоваться в тестировании.
#
# Задаем параметры с помощью стандартной функции
# [`NonLinearReactionAdvectionDiffusionWithFrontData.dparams`](@ref)
using NonLinearReactionAdvectionDiffusionWithFrontData

a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀ = NonLinearReactionAdvectionDiffusionWithFrontData.dparams();
##
u, XX, TP = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ);
nothing #hide

# Генерируем экспериментальные данные функцией
# [`NonLinearReactionAdvectionDiffusionWithFrontData.generate_obs_data`](@ref).
ϕl, ϕr, ϕ, f1_data, f2_data = generate_obs_data(u, Xₙ, N, Tₘ, M, qₙ, ulₘ, urₘ);
nothing # hide


text:
export AbstractNLS # TODO: do not export
export parameter_size, residue_size
export eval_r, eval_r_J
export eval_nls_fobj, eval_nls_∇fobj, eval_nls_∇∇fobj

using LinearAlgebra: dot, mul!
using LinearAlgebra.BLAS: BlasFloat, syrk!, gemv!

@doc raw"""
```julia
abstract type AbstractNLS end
```

Defines an abstract non-linear least squares problem (NLS). In our
context such problem is essentially a differentiable function ``\mathbf{r}``:

```math
\mathbf{r}: \theta\in\mathbb{R}^{n_θ}\mapsto \mathbf{r}(\mathbf{\theta})\in\mathbb{R}^{n_S}
```
where:
- ``\mathbf{r}(\mathbf{θ})∈\mathbb{R}^{n_S}`` is the residue vector,
- ``\mathbf{θ}∈\mathbb{R}^{n_θ}`` is the parameter vector to be optimized

The objective function to minimize is:

```math
f(θ)=\frac{1}{2}\| \mathbf{r}(θ) \|^2
```

The classical approach uses a linear approximation of ``\mathbf{r}``:
```math
\mathbf{r}(\mathbf{θ}+δ\mathbf{θ})\approx \mathbf{r}(\mathbf{θ}) + \mathbf{J}(\mathbf{θ})\cdot δ\mathbf{θ}
```
where ``\mathbf{J}`` is the Jacobian:
```math
\mathbf{J}_{i,j}=\partial_j r^i(\mathbf{θ}),\ i\in[1,n_S],\ j\in[1,n_θ]
```
This leads to
```math
f(\mathbf{θ}+δ\mathbf{θ})\approx f(\mathbf{θ}) + \langle \nabla f, δ\mathbf{θ} \rangle + \frac{1}{2}  \langle \nabla^2 f \cdot δ\mathbf{θ},  δ\mathbf{θ} \rangle
```

Where the gradient ``\nabla f`` is ``\mathbf{J}^t \mathbf{r}`` and the
(approximate) Hessian ``\nabla^2 f`` is ``\mathbf{J}^t \mathbf{J}``.

To implement such model, you must define the following functions:
- [`parameter_size`](@ref) : returns ``n_θ``
- [`residue_size`](@ref) : returns ``n_S``
- [`eval_r`](@ref) : compute ``\mathbf{r}``
- [`eval_r_J`](@ref) : compute ``(\mathbf{r}, \mathbf{J})``
"""
abstract type AbstractNLS end

# ================================================================
# Interface...
# ================================================================
#
@doc raw"""
    parameter_size(nls::AbstractNLS)

Return the dimension ``n_θ`` of the parameter vector ``θ``.
"""
parameter_size(nls::AbstractNLS) = @assert(false,"To implement")

@doc raw"""
    sample_size(nls::AbstractNLS)

Return the dimension ``n_S`` of the residue vector ``r``.
"""
residue_size(nls::AbstractNLS) = @assert(false,"To implement")

@doc raw"""
```julia
eval_r(nls::AbstractNLS,
        θ::AbstractVector) -> r
```

Compte the residual vector ``\mathbf{r}``
"""
eval_r(nls::AbstractNLS,
       θ::AbstractVector) = @assert(false,"To implement")

@doc raw"""
```julia
eval_r_J(nls::AbstractNLS,θ::AbstractVector) -> (r,J)
```

Compute the residual the vector ``\mathbf{r}`` and its Jacobian ``\mathbf{J}``
"""
eval_r_J(nls::AbstractNLS,
         θ::AbstractVector) = @assert(false,"To implement")

# ================================================================
# Convenience functions...
# ================================================================
#

@doc raw"""
```julia
eval_nls_fobj(r::AbstractVector{T}) -> f(θ)
```

Compute ``f(θ)=\frac{1}{2}\| \mathbf{r}(\mathbf{θ}) \|^2``
"""
eval_nls_fobj(r::AbstractVector) = dot(r,r)/2

@doc raw"""
```julia
eval_nls_∇fobj(r,J) -> ∇fobj
```

Compute the gradient: ``\nabla f(\mathbf{θ}) = \mathbf{J}^t\mathbf{r}``
"""
function eval_nls_∇fobj(r::AbstractVector, J::AbstractMatrix)
    J'*r
end

@doc raw"""
```julia
eval_nls_∇∇fobj(J) -> ∇∇fobj
```

Compute the (approximate) Hessian: ``\nabla^2 f(\mathbf{θ}) = \mathbf{J}^t\mathbf{J}``
"""
function eval_nls_∇∇fobj(J::AbstractMatrix)
    Symmetric(J'*J)
end


text:
module Tel

using JLD
using XGBoost
using Logging
@Logging.configure(level=DEBUG)

const data_dir = "../data"
const result_dir = "../result"

include("preprocess.jl")
include("model.jl")
include("eval.jl")
include("utils.jl")
include("reformat.jl")

end


text:
module GeneralizedSVD

import Base: show,getproperty
using LinearAlgebra
using Base: promote_typeof

export GenSVD, gsvd, gsvd!

include("preproc.jl")
include("householderqr.jl")
include("csd.jl")

# Immutable composite type GenSVD
struct GenSVD{T,S} <: Factorization{T}
    U::S
    V::S
    Q::S
    alpha::Vector
    beta::Vector
    k::Int
    l::Int
    R1::S

    # Inner constructor
    function GenSVD{T,S}(U::AbstractMatrix{T}, V::AbstractMatrix{T},
                       Q::AbstractMatrix{T}, alpha::Vector,
                       beta::Vector, k::Int, l::Int,
                       R1::AbstractMatrix{T}) where {T,S}
        new(U, V, Q, alpha, beta, k, l, R1)
    end
end

# Outer constructor
function GenSVD(U::AbstractMatrix{T}, V::AbstractMatrix{T},
              Q::AbstractMatrix{T}, alpha::Vector,
              beta::Vector, k::Int, l::Int,
              R1::AbstractMatrix{T}) where T
    GenSVD{T,typeof(U)}(U, V, Q, alpha, beta, k, l, R1)
end

# iteration for destructuring into components
Base.iterate(S::GenSVD) = (S.U, Val(:V))
Base.iterate(S::GenSVD, ::Val{:V}) = (S.V, Val(:Q))
Base.iterate(S::GenSVD, ::Val{:Q}) = (S.Q, Val(:C))
Base.iterate(S::GenSVD, ::Val{:C}) = (S.C, Val(:S))
Base.iterate(S::GenSVD, ::Val{:S}) = (S.S, Val(:R))
Base.iterate(S::GenSVD, ::Val{:R}) = (S.R, Val(:done))
Base.iterate(S::GenSVD, ::Val{:done}) = nothing

# Explicitely construct C, S, vals and R
@inline function getproperty(F::GenSVD{T}, d::Symbol) where T
    Fa = getfield(F, :alpha)
    Fb = getfield(F, :beta)
    Fk = getfield(F, :k)
    Fl = getfield(F, :l)
    FU = getfield(F, :U)
    FV = getfield(F, :V)
    FQ = getfield(F, :Q)
    FR = getfield(F, :R1)
    if d === :vals
        return Fa[1:Fk + Fl] ./ Fb[1:Fk + Fl]
    elseif d === :C
        m = size(FU, 1)
        if m - Fk - Fl >= 0
            return [Matrix{T}(I, Fk, Fk)  zeros(T, Fk, Fl)            ;
                    zeros(T, Fl, Fk)      Diagonal(Fa[Fk + 1:Fk + Fl]);
                    zeros(T, m - Fk - Fl, Fk + Fl)                    ]
        else
            return [Matrix{T}(I, m, Fk) [zeros(T, Fk, m - Fk); Diagonal(Fa[Fk + 1:m])] zeros(T, m, Fk + Fl - m)]
        end
    elseif d === :S
        m = size(FU, 1)
        p = size(FV, 1)
        if m - Fk - Fl >= 0
            return [zeros(T, Fl, Fk) Diagonal(Fb[Fk + 1:Fk + Fl]); zeros(T, p - Fl, Fk + Fl)]
        else
            return [zeros(T, p, Fk) [Diagonal(Fb[Fk + 1:m]); zeros(T, Fk + p - m, m - Fk)] [zeros(T, m - Fk, Fk + Fl - m); Matrix{T}(I, Fk + p - m, Fk + Fl - m)]]
        end
    elseif d === :R
        n = size(FQ, 1)
        return [zeros(T, Fk + Fl, n - Fk - Fl) FR]
    else
        getfield(F, d)
    end
end

Base.propertynames(F::GenSVD) =
    (:vals, :C, :S, :R, fieldnames(typeof(F))...)

# show() is used for neat display
function show(io::IO, mime::MIME{Symbol("text/plain")},
              F::GenSVD{<:Any,<:AbstractArray})
    summary(io, F); println(io)
    println(io, "U factor:")
    show(io, mime, F.U)
    println(io, "\n V factor:")
    show(io, mime, F.V)
    println(io, "\n Q factor:")
    show(io, mime, F.Q)
    println(io, "\n C factor:")
    show(io, mime, F.C)
    println(io, "\n S factor:")
    show(io, mime, F.S)
    println(io, "\n R factor:")
    show(io, mime, F.R)
    println(io, "\n Generalized singular values:")
    show(io, mime, F.vals)
    println(io, "\n k factor:")
    show(io, mime, F.k)
    println(io, "\n l factor:")
    show(io, mime, F.l)
    println(io, "\n numerical rank of the matrix [A; B]:")
    show(io, mime, F.k + F.l)
end

# promotion type to use for eigenvalues of a Matrix{T}
eigtype(T) = promote_type(Float32, typeof(zero(T)/sqrt(abs2(one(T)))))

copy_oftype(A::AbstractArray{T}, ::Type{T}) where {T} = copy(A)
copy_oftype(A::AbstractArray{T,N}, ::Type{S}) where {T,N,S} = convert(AbstractArray{S,N}, A)

"""
    gsvd(A, B) -> GenSVD

Compute the generalized SVD of an m-by-n matrix `A` and a p-by-n matrix `B`,
returning a `GenSVD` factorization object `F` such that

```math
     A = F.U * F.C * F.R * F.Q',    B = F.V * F.S * F.R * F.Q'
```

# Examples
```jldoctest
julia> A = randn(3, 2); B = randn(4, 2);

julia> F = gsvd(A, B);

julia> [A; B] ≈ [F.U*F.C; F.V*F.S]*F.R*F.Q'
true
```
"""
function gsvd(A::StridedMatrix{T}, B::StridedMatrix{T}) where T<:AbstractFloat
    return gsvd!(copy(A),copy(B))
end

function gsvd(A::StridedMatrix{TA}, B::StridedMatrix{TB}) where {TA,TB}
    S = promote_type(eigtype(TA),TB)
    return gsvd!(copy_oftype(A, S), copy_oftype(B, S))
end

"""
    gsvd!(A, B) -> GenSVD
`gsvd!` is the same as [`gsvd`](@ref), but modifies the arguments
`A` and `B` in-place, instead of making copies.
"""
function gsvd!(A::StridedMatrix{T}, B::StridedMatrix{T}) where T<:AbstractFloat
    m, n = size(A)
    p = size(B)[1]

    if n != size(B)[2]
        throw(ArgumentError("A's column number must equal that of B"))
    end

    # Step 1:
    # preprocess A, B
    U, V, Q, k, l, A, B = preproc!(A, B)

    # Early exit if A23 doesn't exist
    if l == 0
        return GenSVD(U, V, Q, ones(T, k), zeros(T, k), k, 0, A[1:k+l,n-k-l+1:n])
    end

    if m == k
        return GenSVD(U, V, Q, [ones(T, k); zeros(T, l)], [zeros(T, k); ones(T, l)], k, l, [A[1:k,n-k-l+1:n];B[1:l,n-k-l+1:n]])
    end

    # A23 is upper triangular
    if m-k-l >= 0
        @views A23 = A[k+1:k+l, n-l+1:n]
    # A23 is upper trapezoidal
    else
        @views A23 = A[k+1:m, n-l+1:n]
    end
    @views B13 = B[1:l, n-l+1:n]

    # Step 2:
    # QR decomposition of A23 and B13
    Q1, Q2, R23 = householderqr(A23, B13)

    # Step 3:
    # CSD of Q1, Q2
    U1, V1, Z1, cosine, sine = csd2by1!(Q1, Q2)

    # Step 4:
    # set alpha, beta
    alpha = ones(T, k+l)
    @views alpha[k+1:k+l] = cosine
    beta = zeros(T, k+l)
    @views beta[k+1:k+l] = sine

    # Step 5:
    # update U
    t = min(m, k+l)
    U[1:m,k+1:t] = @views U[1:m,k+1:t] * U1

    # Step 6:
    # update V
    V[1:p,1:l] = @views V[1:p,1:l] * V1

    # Step 7:
    # set W
    W = Z1' * R23

    # Step 8:
    # compute RQ decomposition of W
    W_ = copy(W)
    W, tau = LAPACK.gerqf!(W)
    Q3 = LAPACK.orgrq!(W, tau, length(tau))

    # Step 9:
    # update R13(a)
    @views A13 = A[1:k, n-l+1:n]*Q3'
    # form R1
    R1 = zeros(T, k+l, k+l)
    @views R1[1:k, 1:k+1] = A[1:k,n-k-l+1:n-l+1]
    @views R1[1:k, k+1:k+l] = A13
    @views R1[k+1:k+l, k+1:k+l]= W_ * Q3'

    # Step 10:
    # update Q
    @views Q[1:n,n-l+1:n] = Q[1:n,n-l+1:n] * Q3'

    GenSVD(U, V, Q, alpha, beta, k, l, triu!(R1))
end

end


text:
function test_reciprocal_all()
    print(rpad("Math.Series.1/x...", 50, ' '))

    reciprocal = Brainstorm._Math._Series._Reciprocal.reciprocal

    setprecision(4_000_000) do
        @test all([reciprocal(d, x) - trunc(BigInt, (1/big(x))*(big(10)^d)) for x in (1_000_003, 1_000_033, 1_000_037), d in (1, 10, 100, 1_000, 10_000, 100_000, 1_000_000)] .== 0)
    end

    println("PASS")
end


text:
abstract type AbstactBoolean <: BasicDiscreteFunction end

"""
Boolean functions of the type B^n -> B
"""
struct BooleanFunction{T} <: AbstactBoolean
    m :: Vector{T}
end

"""
Boolean map of type B^n -> B^m
"""
struct BooleanMap{T} <: AbstactBoolean
    m :: Matrix{T}
end

const BoolFun = BooleanFunction
const BoolMap = BooleanMap

size(f :: AbstactBoolean) = size(f.m)
arity(f :: BoolFun) = Int(log2(size(f)[1]))
arity(f :: BoolMap) = Int(log2(size(f)[2]))
nfuns(f :: BoolFun) = 1
nfuns(f :: BoolMap) = size(f)[1]

domain(f :: AbstactBoolean) = BCube(arity(f))
codomain(f :: AbstactBoolean) = BCube(nfuns(f))

# standard constructor from domain (and range)
BoolFun(bc :: BCube) = BoolFun(zeros(Int, length(bc)))
BoolMap(dom :: BCube, rng :: BCube) = BoolMap(zeros(Int, (ndims(rng), length(dom))))

# construct from arities - number
BoolFun(arity :: Int) = BoolFun(BCube(arity))
BoolMap(m :: Int, n :: Int) = BoolMap(BCube(m), BCube(n))

# construct from values
BoolFun(tup :: NTuple) = BoolFun(collect(tup))

# apply to BoolVec
(f :: BoolFun)(v :: BoolVec) = BoolVec(f.m[value(v) + 1], 1)
(f :: BoolMap)(v :: BoolVec) = BoolVec(f.m[:, value(v) + 1])
# apply to Iterable
function (f :: AbstactBoolean)(it)
    x = BoolVec(it)
    return f(x)
end
# apply to bunch of numbers
(fs :: AbstactBoolean)(v...) = fs(tuple(v...))

setindex!(f :: BoolFun, v, k :: BoolVec) = (f.m[value(k) + 1] = v)
function setindex!(f :: BoolMap, v, k :: BoolVec)
    col = value(k) + 1
    for row in 1 : nfuns(f)
        f.m[row, col] = v[row]
    end
    return nothing
end
setindex!(f :: AbstactBoolean, v, k) = (f[BoolVec(k)] = v)
setindex!(f :: AbstactBoolean, v, k...) = (f[tuple(k...)] = v)

function direct(fs :: Union{NTuple{N, T}, Vector{T}}) where T <: BoolFun where N
    dom = domain(fs[1])
    dim_dom = ndims(dom)
    dim_rng = length(fs)
    fprod = BoolMap(dim_dom, dim_rng)
    for x in dom
        y = Tuple(f(x) for f in fs)
        fprod[x] = y
    end
    return fprod
end

mat(f :: AbstactBoolean) = f.m
proj(f :: BoolMap, i) = BoolFun(f.m[i, :])
getindex(f :: BoolMap, i) = proj(f, i)
length(f :: BoolMap) = nfuns(f)
iterate(f :: BoolMap) = (f[1], 2)
iterate(f :: BoolMap, i) = (i <= length(f)) ? (f[i], i+1) : nothing

max_stars_repo_name: kirtsar/DiscreteFunctions

max_issues_repo_name: kirtsar/DiscreteFunctions


text:
module Conductor

using ModelingToolkit,
      Catalyst,
      Unitful,
      Unitful.DefaultSymbols,
      InteractiveUtils,
      IfElse,
      Symbolics,
      SymbolicUtils,
      Setfield,
      MacroTools

const MTK = ModelingToolkit

import Symbolics:
    get_variables,
    Symbolic,
    value,
    tosymbol,
    VariableDefaultValue,
    wrap

import ModelingToolkit:
    toparam,
    isparameter,
    Equation,
    defaults,
    AbstractSystem,
    get_eqs,
    get_states,
    get_observed,
    get_defaults,
    get_ps,
    get_systems,
    _merge,
    renamespace,
    hasdefault,
    getdefault,
    setdefault,
    AbstractTimeDependentSystem,
    AbstractSystem,
    independent_variables,
    get_variables!

import ModelingToolkit.SciMLBase: parameterless_type

import Unitful:
    Time,
    Voltage,
    Current,
    Molarity,
    ElectricalConductance

import SymbolicUtils: FnType, Rule
import Unitful: mV, mS, cm, µF, mF, µm, pA, nA, mA, µA, ms, mM, µM
import Base: show, display

export Gate, AlphaBeta, SteadyStateTau, IonChannel, PassiveChannel, SynapticChannel
export EquilibriumPotential, Equilibrium, Equilibria, MembranePotential, IonCurrent
export AuxConversion, D, Network
export Simulation, Concentration, IonConcentration
export @named
export Calcium, Sodium, Potassium, Chloride, Cation, Anion, Leak, Ion, NonIonic
export t
export CompartmentSystem, ConductanceSystem
export output, get_output, timeconstant, steadystate, forward_rate, reverse_rate, hasexponent, exponent
export Sphere, Cylinder, Point, area, radius, height

const ℱ = Unitful.q*Unitful.Na # Faraday's constant
const t = let name = :t; only(@parameters $name) end
const D = Differential(t)
const ExprValues = Union{Expr,Symbol,Number}  # For use in macros

# Metadata IDs
struct ConductorUnits end # temporary shim until we implement MTK's unit checking

isfunction(ex::ExprValues) = try return eval(ex) isa Function catch; return false end

function extract_symbols(ex::ExprValues, out::Vector{Symbol}=[])
    if ~isfunction(ex) && isa(ex, Symbol)
        union!(out, [ex])
    end
    return ex
end

# Basic symbols
function MembranePotential(V0 = -60mV; dynamic = true, name::Symbol = :Vₘ)
    # remove when we have proper unit checks
    V0_val = ustrip(Float64, mV, V0)
    return dynamic ? only(@variables $name(t) = V0_val) : only(@parameters $name = V0_val)
end

@enum Location Outside Inside

# Custom Unitful.jl quantities
@derived_dimension SpecificConductance 𝐈^2*𝐋^-4*𝐌^-1*𝐓^3 # conductance per unit area
@derived_dimension SpecificCapacitance 𝐈^2*𝐋^-4*𝐌^-1*𝐓^4 # capacitance per unit area
@derived_dimension ConductancePerFarad 𝐓^-1 # S/F cancels out to 1/s; perhaps find a better abstract type?

#const TimeF64 = Quantity{Float64, 𝐓, U} where U

include("ions.jl")
include("gates.jl")
include("channels.jl")
include("compartments.jl")
include("networks.jl")
include("io.jl")

#=
function Simulation(network; time::Time, system = false)
    t_val = ustrip(Float64, ms, time)
    simplified = structural_simplify(network)
    if system
        return simplified
    else
        @info repr("text/plain", simplified)
        return ODAEProblem(simplified, [], (0., t_val), [])
    end
end
=#

function Simulation(neuron::CompartmentSystem; time::Time, return_system = false)
    odesys = convert(ODESystem, neuron)
    t_val = ustrip(Float64, ms, time)
    simplified = structural_simplify(odesys)
    if return_system
        return simplified
    else
        @info repr("text/plain", simplified)
        return ODAEProblem(simplified, [], (0., t_val), [])
    end
end

end # module


text:
# Affine coupling layer from Dinh et al. (2017)
# Includes 1x1 convolution from in Putzky and Welling (2019)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

export CouplingLayerGlow, CouplingLayerGlow3D

"""
    CL = CouplingLayerGlow(C::Conv1x1, RB::ResidualBlock; logdet=false)

or

    CL = CouplingLayerGlow(n_in, n_hidden; k1=3, k2=1, p1=1, p2=0, s1=1, s2=1, logdet=false, ndims=2) (2D)

    CL = CouplingLayerGlow(n_in, n_hidden; k1=3, k2=1, p1=1, p2=0, s1=1, s2=1, logdet=false, ndims=3) (3D)

    CL = CouplingLayerGlow3D(n_in, n_hidden; k1=3, k2=1, p1=1, p2=0, s1=1, s2=1, logdet=false) (3D)

 Create a Real NVP-style invertible coupling layer based on 1x1 convolutions and a residual block.

 *Input*:

 - `C::Conv1x1`: 1x1 convolution layer

 - `RB::ResidualBlock`: residual block layer consisting of 3 convolutional layers with ReLU activations.

 - `logdet`: bool to indicate whether to compte the logdet of the layer

 or

 - `n_in`, `n_hidden`: number of input and hidden channels

 - `k1`, `k2`: kernel size of convolutions in residual block. `k1` is the kernel of the first and third
    operator, `k2` is the kernel size of the second operator.

 - `p1`, `p2`: padding for the first and third convolution (`p1`) and the second convolution (`p2`)

 - `s1`, `s2`: stride for the first and third convolution (`s1`) and the second convolution (`s2`)

 - `ndims` : number of dimensions

 *Output*:

 - `CL`: Invertible Real NVP coupling layer.

 *Usage:*

 - Forward mode: `Y, logdet = CL.forward(X)`    (if constructed with `logdet=true`)

 - Inverse mode: `X = CL.inverse(Y)`

 - Backward mode: `ΔX, X = CL.backward(ΔY, Y)`

 *Trainable parameters:*

 - None in `CL` itself

 - Trainable parameters in residual block `CL.RB` and 1x1 convolution layer `CL.C`

 See also: [`Conv1x1`](@ref), [`ResidualBlock`](@ref), [`get_params`](@ref), [`clear_grad!`](@ref)
"""
struct CouplingLayerGlow <: NeuralNetLayer
    C::Conv1x1
    RB::Union{ResidualBlock, FluxBlock}
    logdet::Bool
end

@Flux.functor CouplingLayerGlow

# Constructor from 1x1 convolution and residual block
function CouplingLayerGlow(C::Conv1x1, RB::ResidualBlock; logdet=false)
    RB.fan == false && throw("Set ResidualBlock.fan == true")
    return CouplingLayerGlow(C, RB, logdet)
end

# Constructor from 1x1 convolution and residual Flux block
CouplingLayerGlow(C::Conv1x1, RB::FluxBlock; logdet=false) = CouplingLayerGlow(C, RB, logdet)

# Constructor from input dimensions
function CouplingLayerGlow(n_in::Int64, n_hidden::Int64; k1=3, k2=1, p1=1, p2=0, s1=1, s2=1, logdet=false, ndims=2)

    # 1x1 Convolution and residual block for invertible layer
    C = Conv1x1(n_in)
    RB = ResidualBlock(Int(n_in/2), n_hidden; k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2, fan=true, ndims=ndims)

    return CouplingLayerGlow(C, RB, logdet)
end

CouplingLayerGlow3D(args...;kw...) = CouplingLayerGlow(args...; kw..., ndims=3)

# Forward pass: Input X, Output Y
function forward(X::AbstractArray{Float32, 4}, L::CouplingLayerGlow)

    # Get dimensions
    k = Int(L.C.k/2)

    X_ = L.C.forward(X)
    X1, X2 = tensor_split(X_)

    Y2 = copy(X2)
    logS_T = L.RB.forward(X2)
    S = Sigmoid(logS_T[:,:,1:k,:])
    T = logS_T[:, :, k+1:end, :]
    Y1 = S.*X1 + T
    Y = tensor_cat(Y1, Y2)

    L.logdet == true ? (return Y, glow_logdet_forward(S)) : (return Y)
end

# Inverse pass: Input Y, Output X
function inverse(Y::AbstractArray{Float32, 4}, L::CouplingLayerGlow; save=false)

    # Get dimensions
    k = Int(L.C.k/2)
    Y1, Y2 = tensor_split(Y)

    X2 = copy(Y2)
    logS_T = L.RB.forward(X2)
    S = Sigmoid(logS_T[:,:,1:k,:])
    T = logS_T[:, :, k+1:end, :]
    X1 = (Y1 - T) ./ (S .+ eps(1f0)) # add epsilon to avoid division by 0
    X_ = tensor_cat(X1, X2)
    X = L.C.inverse(X_)

    save == true ? (return X, X1, X2, S) : (return X)
end

# Backward pass: Input (ΔY, Y), Output (ΔX, X)
function backward(ΔY::AbstractArray{Float32, 4}, Y::AbstractArray{Float32, 4}, L::CouplingLayerGlow; set_grad::Bool=true)

    # Recompute forward state
    k = Int(L.C.k/2)
    X, X1, X2, S = inverse(Y, L; save=true)

    # Backpropagate residual
    ΔY1, ΔY2 = tensor_split(ΔY)
    ΔT = copy(ΔY1)
    ΔS = ΔY1 .* X1
    if L.logdet
        set_grad ? (ΔS -= glow_logdet_backward(S)) : (ΔS_ = glow_logdet_backward(S))
    end

    ΔX1 = ΔY1 .* S
    if set_grad
        ΔX2 = L.RB.backward(cat(SigmoidGrad(ΔS, S), ΔT; dims=3), X2) + ΔY2
    else
        ΔX2, Δθrb = L.RB.backward(cat(SigmoidGrad(ΔS, S), ΔT; dims=3), X2; set_grad=set_grad)
        _, ∇logdet = L.RB.backward(cat(SigmoidGrad(ΔS_, S), 0f0.*ΔT; dims=3), X2; set_grad=set_grad)
        ΔX2 += ΔY2
    end
    ΔX_ = tensor_cat(ΔX1, ΔX2)
    if set_grad
        ΔX = L.C.inverse((ΔX_, tensor_cat(X1, X2)))[1]
    else
        ΔX, Δθc = L.C.inverse((ΔX_, tensor_cat(X1, X2)); set_grad=set_grad)[1:2]
        Δθ = cat(Δθc, Δθrb; dims=1)
    end

    if set_grad
        return ΔX, X
    else
        L.logdet ? (return ΔX, Δθ, X, cat(0f0*Δθ[1:3], ∇logdet; dims=1)) : (return ΔX, Δθ, X)
    end
end

## Jacobian-related functions

function jacobian(ΔX::AbstractArray{Float32, 4}, Δθ::Array{Parameter, 1}, X, L::CouplingLayerGlow)

    # Get dimensions
    k = Int(L.C.k/2)

    ΔX_, X_ = L.C.jacobian(ΔX, Δθ[1:3], X)
    X1, X2 = tensor_split(X_)
    ΔX1, ΔX2 = tensor_split(ΔX_)

    Y2 = copy(X2)
    ΔY2 = copy(ΔX2)
    ΔlogS_T, logS_T = L.RB.jacobian(ΔX2, Δθ[4:end], X2)
    S = Sigmoid(logS_T[:,:,1:k,:])
    ΔS = SigmoidGrad(ΔlogS_T[:,:,1:k,:], nothing; x=logS_T[:,:,1:k,:])
    T = logS_T[:, :, k+1:end, :]
    ΔT = ΔlogS_T[:, :, k+1:end, :]
    Y1 = S.*X1 + T
    ΔY1 = ΔS.*X1 + S.*ΔX1 + ΔT
    Y = tensor_cat(Y1, Y2)
    ΔY = tensor_cat(ΔY1, ΔY2)

    # Gauss-Newton approximation of logdet terms
    JΔθ = L.RB.jacobian(zeros(Float32, size(ΔX2)), Δθ[4:end], X2)[1][:, :, 1:k, :]
    GNΔθ = cat(0f0*Δθ[1:3], -L.RB.adjointJacobian(tensor_cat(SigmoidGrad(JΔθ, S), zeros(Float32, size(S))), X2)[2]; dims=1)

    L.logdet ? (return ΔY, Y, glow_logdet_forward(S), GNΔθ) : (return ΔY, Y)
end

function adjointJacobian(ΔY, Y, L::CouplingLayerGlow)
    return backward(ΔY, Y, L; set_grad=false)
end

## Other utils

# Clear gradients
function clear_grad!(L::CouplingLayerGlow)
    clear_grad!(L.C)
    clear_grad!(L.RB)
end

# Get parameters
function get_params(L::CouplingLayerGlow)
    p1 = get_params(L.C)
    p2 = get_params(L.RB)
    return cat(p1, p2; dims=1)
end

# Logdet (correct?)
glow_logdet_forward(S) = sum(log.(abs.(S))) / size(S, 4)
glow_logdet_backward(S) = 1f0./ S / size(S, 4)


text:
export @set, @optic, @reset, @modify, @delete, @insert
using MacroTools

"""
    @set assignment

Return a modified copy of deeply nested objects.

# Example
```jldoctest
julia> using Accessors

julia> struct T;a;b end

julia> t = T(1,2)
T(1, 2)

julia> @set t.a = 5
T(5, 2)

julia> t
T(1, 2)

julia> t = @set t.a = T(2,2)
T(T(2, 2), 2)

julia> @set t.a.b = 3
T(T(2, 3), 2)
```
Supports the same syntax as [`@optic`](@ref). See also [`@reset`](@ref).
"""
macro set(ex)
    setmacro(identity, ex, overwrite=false)
end

"""
    @reset assignment

Shortcut for `obj = @set obj...`.

# Example
```jldoctest
julia> using Accessors

julia> t = (a=1,)
(a = 1,)

julia> @reset t.a=2
(a = 2,)

julia> t
(a = 2,)
```
Supports the same syntax as [`@optic`](@ref). See also [`@set`](@ref).
"""
macro reset(ex)
    setmacro(identity, ex, overwrite=true)
end

"""

    @modify(f, obj_optic)

Define an optic and call [`modify`](@ref) on it.
```jldoctest
julia> using Accessors

julia> xs = (1,2,3);

julia> ys = @modify(xs |> Elements() |> If(isodd)) do x
           x + 1
       end
(2, 2, 4)
```
Supports the same syntax as [`@optic`](@ref). See also [`@set`](@ref).
"""
macro modify(f, obj_optic)
    modifymacro(identity, f, obj_optic)
end

"""

    @delete obj_optic

Define an optic and call [`delete`](@ref) on it.
```jldoctest
julia> using Accessors

julia> xs = (1,2,3);

julia> ys = @delete xs[2]
(1, 3)
```
Supports the same syntax as [`@optic`](@ref). See also [`@set`](@ref).
"""
macro delete(ex)
    obj, optic = parse_obj_optic(ex)
    Expr(:call, delete, obj, optic)
end

"""
    @insert assignment

Return a modified copy of deeply nested objects.

# Example
```jldoctest
julia> using Accessors

julia> t = (a=1, b=2);

julia> @insert t.c = 5
(a = 1, b = 2, c = 5)

julia> t
(a = 1, b = 2)
```

Supports the same syntax as [`@optic`](@ref). See also [`@set`](@ref).
"""
macro insert(ex)
    insertmacro(identity, ex, overwrite=false)
end

"""
    modifymacro(optictransform, f, obj_optic)

This function can be used to create a customized variant of [`@modify`](@ref).
See also [`opticmacro`](@ref), [`setmacro`](@ref).
"""

function modifymacro(optictransform, f, obj_optic)
    f = esc(f)
    obj, optic = parse_obj_optic(obj_optic)
    :(($modify)($f, $obj, $(optictransform)($optic)))
end

foldtree(op, init, x) = op(init, x)
foldtree(op, init, ex::Expr) =
    op(foldl((acc, x) -> foldtree(op, acc, x), ex.args; init=init), ex)

need_dynamic_optic(ex) =
    foldtree(false, ex) do yes, x
        yes || x === :end || (x === :begin) || x === :_
    end

replace_underscore(ex, to) = postwalk(x -> x === :_ ? to : x, ex)

function lower_index(collection::Symbol, index, dim)
    if isexpr(index, :call)
        return Expr(:call, lower_index.(collection, index.args, dim)...)
    elseif index === :end
        if dim === nothing
            return :($(Base.lastindex)($collection))
        else
            return :($(Base.lastindex)($collection, $dim))
        end
    elseif (index === :begin)
        if dim === nothing
            return :($(Base.firstindex)($collection))
        else
            return :($(Base.firstindex)($collection, $dim))
        end
    end
    return index
end

function parse_obj_optics(ex)
    if @capture(ex, (front_ |> back_))
        obj, frontoptic = parse_obj_optics(front)
        backoptic = try
            # allow e.g. obj |> first |> _.a.b
            obj_back, backoptic = parse_obj_optics(back)
            if obj_back == esc(:_)
                backoptic
            else
                (esc(back),)
            end
        catch ArgumentError
            backoptic = (esc(back),)
        end
        return obj, tuple(frontoptic..., backoptic...)
    elseif @capture(ex, front_[indices__])
        obj, frontoptic = parse_obj_optics(front)
        if any(need_dynamic_optic, indices)
            @gensym collection
            indices = replace_underscore.(indices, collection)
            dims = length(indices) == 1 ? nothing : 1:length(indices)
            lindices = esc.(lower_index.(collection, indices, dims))
            optic = :($DynamicIndexLens($(esc(collection)) -> ($(lindices...),)))
        else
            index = esc(Expr(:tuple, indices...))
            optic = :($IndexLens($index))
        end
    elseif @capture(ex, front_.property_)
        property isa Union{Symbol,String} || throw(ArgumentError(
            string("Error while parsing :($ex). Second argument to `getproperty` can only be",
                   "a `Symbol` or `String` literal, received `$property` instead.")
        ))
        obj, frontoptic = parse_obj_optics(front)
        optic = :($PropertyLens{$(QuoteNode(property))}())
    elseif @capture(ex, f_(front_))
        obj, frontoptic = parse_obj_optics(front)
        optic = esc(f) # function optic
    else
        obj = esc(ex)
        return obj, ()
    end
    return (obj, tuple(frontoptic..., optic))
end

"""
    opticcompose([optic₁, [optic₂, [optic₃, ...]]])

Compose `optic₁`, `optic₂` etc. There is one subtle point here:
While the two composition orders `(optic₁ ⨟ optic₂) ⨟ optic₃` and `optic₁ ⨟ (optic₂ ⨟ optic₃)` have equivalent semantics, their performance may not be the same.

The `opticcompose` function tries to use a composition order, that the compiler likes. The composition order is therefore not part of the stable API.
"""
opticcompose() = identity
opticcompose(args...) = opcompose(args...)

function parse_obj_optic(ex)
    obj, optics = parse_obj_optics(ex)
    optic = Expr(:call, opticcompose, optics...)
    obj, optic
end

function get_update_op(sym::Symbol)
    s = String(sym)
    if !endswith(s, '=') || isdefined(Base, sym)
        # 'x +=' etc. is actually 'x = x +', and so '+=' isn't defined in Base.
        # '>=' however is a function, and not an assignment operator.
        msg = "Operation $sym doesn't look like an assignment"
        throw(ArgumentError(msg))
    end
    Symbol(s[1:end-1])
end

struct _UpdateOp{OP,V}
    op::OP
    val::V
end
(u::_UpdateOp)(x) = u.op(x, u.val)

"""
    setmacro(optictransform, ex::Expr; overwrite::Bool=false)

This function can be used to create a customized variant of [`@set`](@ref).
It works by applying `optictransform` to the optic that is used in the customized `@set` macro
at runtime.

# Example
```julia
function mytransform(optic::Lens)::Lens
    ...
end
macro myset(ex)
    setmacro(mytransform, ex)
end
```
See also [`opticmacro`](@ref).
"""
function setmacro(optictransform, ex::Expr; overwrite::Bool=false)
    @assert ex.head isa Symbol
    @assert length(ex.args) == 2
    ref, val = ex.args
    obj, optic = parse_obj_optic(ref)
    dst = overwrite ? obj : gensym("_")
    val = esc(val)
    ret = if ex.head == :(=)
        :($dst = $set($obj, ($optictransform)($optic), $val))
    else
        op = get_update_op(ex.head)
        f = :($_UpdateOp($op,$val))
        :($dst = $modify($f, $obj, ($optictransform)($optic)))
    end
    ret
end

"""
    insertmacro(optictransform, ex::Expr; overwrite::Bool=false)

This function can be used to create a customized variant of [`@insert`](@ref).
It works by applying `optictransform` to the optic that is used in the customized `@insert` macro at runtime.

# Example
```julia
function mytransform(optic::Lens)::Lens
    ...
end
macro myinsert(ex)
    insertmacro(mytransform, ex)
end
```
See also [`opticmacro`](@ref),  [`setmacro`](@ref).
"""
function insertmacro(optictransform, ex::Expr; overwrite::Bool=false)
    @assert ex.head isa Symbol
    @assert length(ex.args) == 2
    @assert ex.head == :(=)
    ref, val = ex.args
    obj, optic = parse_obj_optic(ref)
    dst = overwrite ? obj : gensym("_")
    val = esc(val)
    :($dst = $insert($obj, ($optictransform)($optic), $val))
end

"""
    @optic

Construct an optic from property access and similar.

# Example

```jldoctest
julia> using Accessors

julia> struct T;a;b;end

julia> t = T("A1", T(T("A3", "B3"), "B2"))
T("A1", T(T("A3", "B3"), "B2"))

julia> l = @optic _.b.a.b
(@optic _.b) ∘ (@optic _.a) ∘ (@optic _.b)

julia> l(t)
"B3"

julia> set(t, l, 100)
T("A1", T(T("A3", 100), "B2"))

julia> t = ("one", "two")
("one", "two")

julia> set(t, (@optic _[1]), "1")
("1", "two")
```

See also [`@set`](@ref).
"""
macro optic(ex)
    opticmacro(identity, ex)
end

"""
    opticmacro(optictransform, ex::Expr)

This function can be used to create a customized variant of [`@optic`](@ref).
It works by applying `optictransform` to the created optic at runtime.
```julia
# new_optic = mytransform(optic)
macro myoptic(ex)
    opticmacro(mytransform, ex)
end
```
See also [`setmacro`](@ref).
"""
function opticmacro(optictransform, ex)
    obj, optic = parse_obj_optic(ex)
    if obj != esc(:_)
        msg = """Cannot parse optic $ex. Lens expressions must start with _, got $obj instead."""
        throw(ArgumentError(msg))
    end
    :($(optictransform)($optic))
end

_show(io::IO, optic::PropertyLens{field}) where {field} = print(io, "(@optic _.$field)")
_show(io::IO, optic::IndexLens) = print(io, "(@optic _[", join(repr.(optic.indices), ", "), "])")
Base.show(io::IO, optic::Union{IndexLens, PropertyLens}) = _show(io, optic)
Base.show(io::IO, ::MIME"text/plain", optic::Union{IndexLens, PropertyLens}) = _show(io, optic)

# debugging
show_composition_order(optic) = (show_composition_order(stdout, optic); println())
show_composition_order(io::IO, optic) = show(io, optic)
function show_composition_order(io::IO, optic::ComposedOptic)
    print(io, "(")
    show_composition_order(io, optic.outer)
    print(io, " ∘  ")
    show_composition_order(io, optic.inner)
    print(io, ")")
end


text:
"""
    ArrayPool{T}

Helps to maintain the pool of reusable arrays of different sizes
and reduce the burden on garbage collection.

# Type parameters
 * `T`: the type of array elements

"""
struct ArrayPool{T}
    len_pools_lock::Threads.SpinLock        # lock for thread-safe access to the pool
    len_pools::Dict{Int, Vector{Vector{T}}} # pools of vectors of different lengths

    ArrayPool{T}() where T =
        new{T}(Threads.SpinLock(), Dict{Int, Vector{Vector{T}}}())
end

"""
    acquire!(pool::ArrayPool{T}, size) where T -> Array{T}

Gets an array of a specific size from the pool.
The acquired array should be returned back to the pool using [`release!`](@ref).
The `size` could be either an integer or a tuple of integers.
"""
function acquire!(pool::ArrayPool{T}, len::Integer) where T
    lock(pool.len_pools_lock)
    len_pool = haskey(pool.len_pools, len) ?
               pool.len_pools[len] :
               get!(() -> Vector{Vector{T}}(), pool.len_pools, len)
    res = isempty(len_pool) ? Vector{T}(undef, len) : pop!(len_pool)
    unlock(pool.len_pools_lock)
    return res
end

acquire!(pool::ArrayPool, size::Tuple) =
    reshape(acquire!(pool, prod(size)), size)

"""
    release!(pool::ArrayPool{T}, arr::Array{T}) where T

Releases the array previously obtained by [`acquire!`](@ref) back into the pool.
"""
function release!(pool::ArrayPool{T}, arr::Array{T}) where T
    len = length(arr)
    len_pool = get(pool.len_pools, len, nothing)
    if len_pool !== nothing
        #@info "release($(size(arr))) ($(length(len_pool)))"
        (length(len_pool) <= 100) || error("Overflow of $len-sized vectors pool")
        lock(pool.len_pools_lock)
        push!(len_pool, vec(arr))
        unlock(pool.len_pools_lock)
    else
        throw(DimensionMismatch("No $len-element arrays were acquired before"))
    end
    return pool
end


text:
"""\nFunction that creates Utils struct from model definition\n# Arguments\n- `mdl::Model`: mutable struct Model.\n- `method_cube::String="CPU"`: the method for the utils cube. Try "CPU", "GPU","Mmap","Function".\n# Examples\n```julia-repl\njulia> execute_cube(my_model, "CPU")\n```\n"""
function execute_cube(mdl::Model,method_cube::String="CPU")

    (method_cube == "CPU") && ((cubefun, precompilation) = (cube_CPU, true))
    (method_cube == "GPU") && ((cubefun, precompilation) = (cube_GPU, true))
    (method_cube == "Mmap") && ((cubefun, precompilation) = (cube_Mmap, true))
    (method_cube == "Function") && ((cubefun, precompilation) = (cube_function, false))

    ( ! @isdefined cubefun ) && error("1.Method provided for Cube is incorrect. Check for spelling errors.")

    # precompilation = true

    # if (method_cube == "Function")

    #     precompiletion = false

    # end

    if precompilation == true

        (method_cube == "CPU") && (qb = cube_CPU_qb(mdl))
        (method_cube == "GPU") && (qb = cube_GPU_qb(mdl))
        (method_cube == "Mmap") && (qb = cube_Mmap_qb(mdl))
        return( Utils(cubefun,qb) )

    end

    return( Utils(cubefun) )

end


text:
"""
    $FUNCTIONNAME(pc::PlayerCost, x::SVector{nx}, u::SVector{nu}, t::AbstractFloat)

A convencience implementation of the cost quadraticization.
"""
function _quadraticize_ad(pc::PlayerCost, g::GeneralGame, x::SVector, u::SVector,
                          t::AbstractFloat)
    @warn "You are using the fallback quadraticization using ForwardDiff.
    Consider implementing a custom `quadraticize` for your `ControlSystem`
    type." maxlog=1

    x_cost = x->pc(g, x, u, t)
    u_cost = u->pc(g, x, u, t)
    # we can compute the gradient and the hessian in one go
    diff_x = DiffResults.HessianResult(x)
    diff_x = ForwardDiff.hessian!(diff_x, x_cost, x)
    # the linear state component of the cost is the gradient in x
    l = DiffResults.gradient(diff_x)
    # the quadratic state component of the cost is the hessian in x
    Q = DiffResults.hessian(diff_x)

    diff_u = DiffResults.HessianResult(u)
    diff_u = ForwardDiff.hessian!(diff_u, u_cost, u)
    # the linear control cost is the gradient in u
    r = DiffResults.gradient(diff_u)
    # the quadratic control component o the cost is the hessian in u
    R = DiffResults.hessian(diff_u)

    return QuadraticPlayerCost(l, Q, r, R)
end

quadraticize!(qcache::QuadCache, pc::PlayerCost, g::GeneralGame, x::SVector,
              u::SVector, t::AbstractFloat) = _quadraticize_ad(pc, g, x, u, t)


text:
using Optim, HTTP, GLM, LinearAlgebra, Random, Statistics, DataFrames, CSV, FreqTables
using Pkg
Pkg.add("Distributions")
using Pkg
Pkg.add("ForwardDiff")

    #Question1
    using DataFrames
using CSV
using HTTP
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2021/master/ProblemSets/PS4-mixture/nlsw88t.csv"
df = CSV.read(HTTP.get(url).body, DataFrame)
X = [df.age df.white df.collgrad]
Z = hcat(df.elnwage1, df.elnwage2, df.elnwage3, df.elnwage4,
df.elnwage5, df.elnwage6, df.elnwage7, df.elnwage8)
y = df.occ_code
#Using Q1 from Pset3
function mlogit_with_Z(theta, X, Z, y)

    alpha = theta[1:end-1]
    gamma = theta[end]
    K = size(X,2)
    J = length(unique(y))
    N = length(y)
    bigY = zeros(N,J)
    for j=1:J
        bigY[:,j] = y.==j
    end
    bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

    T = promote_type(eltype(X),eltype(theta))
    num   = zeros(T,N,J)
    dem   = zeros(T,N)
    for j=1:J
        num[:,j] = exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)
        dem .+= num[:,j]
    end

    P = num./repeat(dem,1,J)

    loglike = -sum( bigY.*log.(P) )

    return loglike
end
startvals = [2*rand(7*size(X,2)).-1; .1]
td = TwiceDifferentiable(theta -> mlogit_with_Z(theta, X, Z, y), startvals; autodiff = :forward)
# run the optimizer
theta_hat_optim_ad = optimize(td, startvals, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
theta_hat_mle_ad = theta_hat_optim_ad.minimizer
# evaluate the Hessian at the estimates
H  = Optim.hessian!(td, theta_hat_mle_ad)
theta_hat_mle_ad_se = sqrt.(diag(inv(H)))
println([theta_hat_mle_ad theta_hat_mle_ad_se]
#Question 2
#Interpretation of gamma#
#Question 3
n.
using Distributions
include("lgwt.jl") # make sure the function gets read in
# define distribution
d = Normal(0,1) # mean=0, standard deviation=1
e.
# get quadrature nodes and weights for 7 grid points
nodes, weights = lgwt(7,-4,4)
# now compute the integral over the density and verify it's 1
sum(weights.*pdf.(d,nodes))
# now compute the expectation and verify it's 0
sum(weights.*nodes.*pdf.(d,nodes))
#computing integrals using the above
d = Normal(0,2)
nodes, weights = lgwt(7,-5,5)
sum(weights.*(nodes.).*pdf.(d,nodes))

d = Normal(0,2)
nodes, weights = lgwt(10,-5,5)
sum(weights.*(nodes.).*pdf.(d,nodes))

#Question4


text:
# get the baseline of an ECG using median filters

"""
    getECGBaseline(ecg::Vector{<:Real}, samplerate::Real)

Get the baseline of an ECG signal for baseline correction. Source: Advances in Cardiac Signal Processing - Acharya, U.R. and Suri, J. and Spaan, J.A.E. and Krishnan, S.M. and Technologies, B.
    - ISBN: 9783540366751 page: 58f. adaption by Jan F. Kraemer

# Args:

* ecg: ECG signal
* samplerate: Sampling rate
"""
function getECGBaseline(ecg::Vector{<:Real}, samplerate::Real)
    samplerate > 10 || throw(DomainError("Samplerate needs to be positive and larger than 10 Hz, is: $samplerate"))

    # Baseline drift via median filter
    # 1. get rid of QRS
    baseline = medfilt1(ecg, n = Int(ceil(0.25*samplerate)))
    # 2. get rid of P and T waves
    baseline = medfilt1(baseline, n = Int(ceil(0.75*samplerate)))

    # Lowpass filter to avoid sudden changes
    low = Lowpass(5, fs = samplerate)
    cheby = Chebyshev2(5, 10)
    lowChebyFilt = digitalfilter(low, cheby)
    baseline = filtfilt(lowChebyFilt, baseline)

    return baseline
end


text:
using CxxWrap: ConstCxxRef
using CxxWrap.StdLib: StdString, SharedPtr, UniquePtr, WeakPtr
using CxxWrap.CxxWrapCore:
    CxxBool,
    CxxChar,
    CxxWchar,
    CxxSignedChar,
    CxxUChar,
    CxxShort,
    CxxUShort,
    CxxInt,
    CxxUInt,
    CxxLong,
    CxxULong,
    CxxLongLong,
    CxxULongLong,
    SmartPointer