text:
######## Fractional Brownian Motion ########
"""
Fractional Brownian motion.

# Members
- hurst: the Hurst exponent
"""
struct FractionalBrownianMotion <: SSSIProcess
    hurst::Real

    function FractionalBrownianMotion(hurst::Real)
        0. < hurst < 1. || error("Hurst exponent must be bounded in 0 and 1.")
        new(hurst)
    end
end

ss_exponent(X::FractionalBrownianMotion) = X.hurst

@doc raw""" fBm_autocov(t,s,H)

Autocovariance function of standard fBm:
    $\frac 1 2 (|t|^{2H} + |s|^{2H} - |t-s|^{2H})$
"""
fBm_autocov = (t::Real,s::Real,H::Real) -> 1/2 * (abs(t)^(2H) + abs(s)^(2H) - abs(t-s)^(2H))

autocov(X::FractionalBrownianMotion, t::Real, s::Real) = fBm_autocov(t,s,X.hurst)

"""
    fBm_covmat(G1::AbstractVector{<:Real}, G2::AbstractVector{<:Real}, H::Real)

Compute the covariance matrix of fBm with hurst `H` between two grids `G1` and `G2`.
This gives the same result as `covmat(FractionalBrownianMotion(H), G1, G2)`.
"""
function fBm_covmat(G1::AbstractVector{<:Real}, G2::AbstractVector{<:Real}, H::Real)
    Σ = zeros(length(G1),length(G2))
    for (c,s) in enumerate(G2), (r,t) in enumerate(G1)
        Σ[r,c] = fBm_autocov(t, s, H)
    end
    return Σ
end

fBm_covmat(G, H) = Matrix(Symmetric(fBm_covmat(G, G, H)))

# Moving average kernels of fBm
"""
K_+ kernel
"""
function Kplus(x::Real, t::Real, H::Real)
    # @assert t>0
    p::Real = H-1/2
    v::Real = 0
    if x<0
        v = (t-x)^p - (-x)^p
    elseif 0<=x<t
        v = (t-x)^p
    else
        v = 0
    end
    return v
end

"""
K_- kernel
"""
function Kminus(x::Real, t::Real, H::Real)
    p::Real = H-1/2
    v::Real = 0
    if x<=0
        v = 0
    elseif 0<x<=t
        v = -(x)^p
    else
        v = (x-t)^p - (x)^p
    end
    return v
end

"""
K_+ + K_- kernel
"""
Kppm(x, t, H) = Kplus(x, t, H) + Kminus(x, t, H)

"""
K_+ - K_- kernel
"""
Kpmm(x, t, H) = Kplus(x, t, H) - Kminus(x, t, H)

######## Fractional Gaussian Noise ########

@doc raw"""
Fractional Gaussian noise (fGn) is the discrete time version of the continuous time differential process of a fBm: $ B(t+\delta) - B(t) $.  It is defined as $ B(n+l) - B(n) $, where `l` is the lag.

# Note
- The anti-causal convention adopted by the definition here has no impact since the process is stationary.
"""
struct FractionalGaussianNoise <: IncrementProcess{FractionalBrownianMotion}
    parent_process::FractionalBrownianMotion
    lag::Integer

    function FractionalGaussianNoise(hurst::Real, lag::Integer=1)
        lag >= 1 || error("Lag must be >= 1.")
        new(FractionalBrownianMotion(hurst), lag)
    end
end

ss_exponent(X::FractionalGaussianNoise) = X.parent_process.hurst
lag(X::FractionalGaussianNoise) = X.lag
step(X::FractionalGaussianNoise) = 1

@doc raw"""
    fGn_autocov(t::Real,H::Real,δ::Real)

Autocovariance function of stardard (continuous time) fGn: $ 1/2 (|t+δ|^{2H} + |t-δ|^{2H} - 2|t|^{2H}) $
"""
fGn_autocov = (t::Real,H::Real,δ::Real) -> 1/2 * (abs(t+δ)^(2H) + abs(t-δ)^(2H) - 2*abs(t)^(2H))

function autocov(X::FractionalGaussianNoise, n::Integer)
    return fGn_autocov(n, ss_exponent(X), lag(X))
end

"""
    fGn_covmat(G1::AbstractVector{<:Real}, G2::AbstractVector{<:Real}, H::Real, δ::Real)

Compute the covariance matrix of standard (continuous time) fractional Gaussian noise of Hurst exponent `H` and time lag `δ` between the grid `G1` and `G2`.

# Notes
The special case of discrete regular grids with integer lag `d` is equivalent to `covmat(FractionalGaussianNoise(H, d), 1:N, 1:M)`.
"""
function fGn_covmat(G1::AbstractVector{<:Real}, G2::AbstractVector{<:Real}, H::Real, δ::Real)
    Σ = zeros(length(G1),length(G2))
    for (c,s) in enumerate(G2), (r,t) in enumerate(G1)
        Σ[r,c] = fGn_autocov(t-s, H, δ)
    end
    return Σ
end

fGn_covmat(G::AbstractVector{<:Real}, H::Real, δ::Real) = Matrix(Symmetric(fGn_covmat(G, G, H, δ)))

fGn_covmat(N::Integer, H::Real, d::Integer) = fGn_covmat(1:N, H, d)

######## Fractional Wavelet noise (fWn) ########

@doc raw"""
A fractional Wavelet noise (fWn) is discrete time version of the high-pass filtered process of fBm:
$$
W(t) = \sum_{n=0}^{p-1} B(t-n\delta) \psi[n]
$$

It is defined as
$$
W[k] = \sum_{n=0}^{p-1} B(k-n) \psi[n]
$$

# Notes
- fGn is a special case of fWn.
```julia
d = 1  # lag for fGn
fGn = FracFin.FractionalGaussianNoise(0.7, d)

filter = vcat(1,zeros(d-1),-1)
fWn = FracFin.FractionalWaveletNoise(0.7, filter)

G = FracFin.covmat(fGn, 10)
W = FracFin.covmat(fWn, 10)
maximum(abs.(G-W))  # close to 0
```
"""
struct FractionalWaveletNoise <: FilteredProcess{DiscreteTime, FractionalBrownianMotion}
    parent_process::FractionalBrownianMotion
    filter::AbstractVector{<:Real}
    coeff::AbstractVector{<:Real}  # coefficients for the summation of |t + k|^(2H)
    supp::AbstractVector{<:Integer}  # support (index range) of the coefficients

    function FractionalWaveletNoise(hurst::Real, filter::AbstractVector{<:Real})
        @assert isapprox(sum(filter), 0; atol=1e-8)  "Filter must not contain DC."

        # whatever the mode of transform the support is always 1-length(filter):length(filter)-1, and the coefficients satisfy coeffs[k] = coeffs[-k]
        new(FractionalBrownianMotion(hurst), filter, fWn_cov_coeff(filter), 1-length(filter):length(filter)-1)
    end
end

ss_exponent(X::FractionalWaveletNoise) = X.parent_process.hurst
filter(X::FractionalWaveletNoise) = X.filter
step(X::FractionalWaveletNoise) = 1

"""
# Notes
- `fWn_cov_coeff(u,v)` is the reverse of `fWn_cov_coeff(v,u)`.
"""
function fWn_cov_coeff(u::AbstractVector{<:Real}, v::AbstractVector{<:Real})
    m = length(u)
    n = length(v)

    return -1 * native_conv(u, reverse(v))
end

fWn_cov_coeff(u::AbstractVector{<:Real}) = fWn_cov_coeff(u, u)

"""
Autocovariance function of a continuous time fWn.
"""
function fWn_autocov(X::FractionalWaveletNoise, t::Real, δ::Real; mode::Symbol=:causal)
    return 1/2 * sum(X.coeff .* abs.(t .+ δ * X.supp).^(2*ss_exponent(X)))
end

autocov(X::FractionalWaveletNoise, n::Integer) = fWn_autocov(X, n, 1)

"""
    covmat(X::FractionalWaveletNoise, G1::AbstractVector{<:Real}, G2::AbstractVector{<:Real}, δ::Real=1)

Compute the covariance matrix between the grid `G1` and `G2` of a continuous time fWn `X` with step `δ`.

# Notes
- Like `fGn_covmat()`, this function is equivalent to `covmat()` in case of regular grids.
"""
function covmat(X::FractionalWaveletNoise, G1::AbstractVector{<:Real}, G2::AbstractVector{<:Real}, δ::Real=1)
    Σ = zeros(length(G1),length(G2))
    for (c,t) in enumerate(G2), (r,s) in enumerate(G1)
        Σ[r,c] = fWn_autocov(X, s-t, δ)
    end

    return Σ
end

covmat(X::FractionalWaveletNoise, G::AbstractVector{<:Real}, δ::Real=1) = Matrix(Symmetric(covmat(X, G, G, δ)))

## Redundant functions
function fWn_covmat(G1::AbstractVector{<:Real}, G2::AbstractVector{<:Real}, filter::AbstractVector{<:Real}, H::Real, δ::Real)
    coeff = fWn_cov_coeff(filter)
    n = length(filter)
    Σ = zeros(length(G1),length(G2))

    for (c,t) in enumerate(G2), (r,s) in enumerate(G1)
        Σ[r,c] = 1/2 * sum(coeff .* abs.((s-t) .+ δ*(1-n:n-1)).^(2H))
    end

    return Σ
end

fWn_covmat(G::AbstractVector{<:Real}, f::AbstractVector{<:Real}, H::Real, δ::Real) = Matrix(Symmetric(fWn_covmat(G, G, f, H, δ)))

######## fWn bank (fWnb) ########

"""
A fractional Wavelet noise bank (fWnb) is a fBm filtered by a bank of filters, in other words, a collection of fWn.

# Notes
- This process is multivariate hence can not be properly handled in the scalar framework as fWn.
"""
struct FractionalWaveletNoiseBank <: AbstractRandomField # <: FilteredProcess{DiscreteTime, FractionalBrownianMotion}
    parent_process::FractionalBrownianMotion
    filters::AbstractVector
    coeffs::AbstractMatrix  # 2-tuple: pair coefficients for the summation of |t + k|^(2H) and the range of index
    supps::AbstractMatrix  # supports (index range) of the coefficients

    function FractionalWaveletNoiseBank(hurst::Real, filters::AbstractVector{<:AbstractVector{<:Real}}, mode::Symbol)
        @assert all(isapprox(sum(f), 0; atol=1e-8) for f in filters)  "Filters must not contain DC."

        # Gotcha! Double for-loops in array comprehension has reversed order than the usual double for-loop, e.g.
        # [(r,c) for r=1:3, c=1:2] is a 3-by-2 matrix
        coeffs = [fWn_cov_coeff(u,v) for u in filters, v in filters]  # the first for-loop variable `u` is the row dimension

        # Supports, following the same convention as in `convmask()`, but as correlation not convolution.
        supps = if mode in [:causal, :left]  # <-> :anticausal in `convmask`
            [(1-length(u):length(v)-1) for u in filters, v in filters]
        elseif mode in [:anticausal, :right]
            [(1-length(v):length(u)-1) for u in filters, v in filters]
        else
            [(1-floor(Int,length(v)/2)-ceil(Int,length(u)/2)):(floor(Int,length(u)/2)+ceil(Int,length(v)/2)-1) for u in filters, v in filters]
        end

        new(FractionalBrownianMotion(hurst), filters, coeffs, supps)
    end
end

ss_exponent(X::FractionalWaveletNoiseBank) = X.parent_process.hurst
filter(X::FractionalWaveletNoiseBank) = X.filters
step(X::FractionalWaveletNoiseBank) = 1

function fWn_autocov(A::AbstractMatrix{<:AbstractVector{<:Real}}, K::AbstractMatrix{<:AbstractVector{<:Real}}, H::Real, t::Real, δ::Real; partial::Bool=false, kwargs...)
    @assert size(A) == size(K) && size(A,1) == size(A,2)  "Mismatched dimensions."
    @assert 0 < H < 1  "Invalid Hurst exponent."

    J = size(A, 1)
    Σ = zeros(Real, (J,J))

    if partial
        for r=1:J
            Σ[r,r] = 1/2 * sum(A[r,r] .* abs.(t .+ δ*K[r,r]).^(2H))
        end
    else
        for c=1:J, r=1:J
            Σ[r,c] = 1/2 * sum(A[r,c] .* abs.(t .+ δ*K[r,c]).^(2H))
        end
    end
    return Σ
end

"""
Covariance of fWnb. This covariance is a matrix since fWnb is multivariate process.

# Args
- X: a fWnb process
- t: time
- δ: sampling step

# Notes
- filter at [r,c] is the reverse of that at [c,r]
- range of index at [r,c] is reverse of that at [c,r] times -1
By consequent the final matrix Σ(t) = Σ'(-t).
"""
function autocov(X::FractionalWaveletNoiseBank, t::Real, δ::Real; kwargs...)
    fWn_autocov(X.coeffs, X.supps, ss_exponent(X), t, δ; kwargs...)
end

# function autocov(X::FractionalWaveletNoiseBank, t::Real, δ::Real)
#     J = length(X.filters)
#     H = ss_exponent(X)
#     Σ = zeros(Real, (J,J))

#     for c=1:J, r=1:J
#         # X.coeffs[r,c] is (filter, range of index)
#         Σ[r,c] = 1/2 * sum(X.coeffs[r,c] .* abs.(t .+ δ*X.supps[r,c]).^(2H))
#     end
#     return Σ
# end

function fWn_covmat(A::AbstractMatrix{<:AbstractVector{<:Real}}, K::AbstractMatrix{<:AbstractVector{<:Real}}, H::Real, G1::AbstractVector{<:Real}, G2::AbstractVector{<:Real}; kwargs...)
    @assert size(A) == size(K) && size(A,1) == size(A,2)  "Dimension mismatches."

    J = size(A,1)
    l1, l2 = length(G1), length(G2)
    Σ = zeros(l1*J, l2*J)

    for c=0:l2-1,r=0:l1-1  # inner (row) iteration is on `r`
        Σ[(r*J+1):(r*J+J), (c*J+1):(c*J+J)] = fWn_autocov(A, K, H, G1[r+1]-G2[c+1], δ; kwargs...)
    end
    return Σ  # not forcing symmetry
end

function fWn_covmat(A::AbstractMatrix{<:AbstractVector{<:Real}}, K::AbstractMatrix{<:AbstractVector{<:Real}}, H::Real, G::AbstractVector{<:Real}; kwargs...)
    @assert size(A) == size(K) && size(A,1) == size(A,2)  "Dimension mismatches."

    J = size(A,1)
    l = length(G)
    Σ = zeros(l*J, l*J)

    for c=0:l-1,r=0:l-1
        # use the property:
        # autocov(x, t) = transpose(autocov(x, -t))
        Σ[(r*J+1):(r*J+J), (c*J+1):(c*J+J)] = (r>=c) ? fWn_autocov(A, K, H, G[r+1]-G[c+1], δ; kwargs...) : transpose(Σ[(c*J+1):(c*J+J), (r*J+1):(r*J+J)])
    end

    return Matrix(Symmetric(Σ))  # forcing symmetry
end

function fWn_covmat(A::AbstractMatrix{<:AbstractVector{<:Real}}, K::AbstractMatrix{<:AbstractVector{<:Real}}, H::Real, G::AbstractVector{<:Integer}; kwargs...)
    @assert size(A) == size(K) && size(A,1) == size(A,2)  "Dimension mismatches."

    J = size(A,1)
    l = length(G)
    Σ = zeros(l*J, l*J)

    Σs = [fWn_autocov(A, K, H, G[n]-G[1], 1; kwargs...) for n=1:length(G)]

    for c=0:l-1,r=0:l-1
        Σ[(r*J+1):(r*J+J), (c*J+1):(c*J+J)] = (r>=c) ? Σs[r-c+1] : transpose(Σs[c-r+1])
    end

    return Matrix(Symmetric(Σ))  # forcing symmetry
end

"""
Compute the covariance matrix of standard (continuous time) fractional Wavelet noise bank of Hurst exponent `H` and time lag `δ` between the grid `G1` and `G2`.
"""
function covmat(X::FractionalWaveletNoiseBank, G1::AbstractVector{<:Real}, G2::AbstractVector{<:Real}, δ::Real; kwargs...)
    return fWn_covmat(X.coeffs, X.supps, ss_exponent(X), G1, G2, δ; kwargs...)
end

function covmat(X::FractionalWaveletNoiseBank, G::AbstractVector{<:Real}, δ::Real; kwargs...)
    return fWn_covmat(X.coeffs, X.supps, ss_exponent(X), G, δ; kwargs...)
end

function covmat(X::FractionalWaveletNoiseBank, G::AbstractVector{<:Integer}; kwargs...)
    return fWn_covmat(X.coeffs, X.supps, ss_exponent(X), G; kwargs...)
end

# """
# Compute the covariance matrix of standard (continuous time) fractional Wavelet noise bank of Hurst exponent `H` and time lag `δ` between the grid `G1` and `G2`.
# """
# function covmat(X::FractionalWaveletNoiseBank, G1::AbstractVector{<:Real}, G2::AbstractVector{<:Real}, δ::Real)
#         J = length(X.filters)
#     l1, l2 = length(G1), length(G2)
#     Σ = zeros(l1*J, l2*J)

#     for c=0:l2-1,r=0:l1-1  # inner (row) iteration is on `r`
#         Σ[(r*J+1):(r*J+J), (c*J+1):(c*J+J)] = autocov(X, G1[r+1]-G2[c+1], δ)
#     end
#     return Σ  # not forcing symmetry
# end

# function covmat(X::FractionalWaveletNoiseBank, G::AbstractVector{<:Real}, δ::Real)
#     J = length(X.filters)
#     l = length(G)
#     Σ = zeros(l*J, l*J)

#     for c=0:l-1,r=0:l-1
#         # use the property:
#         # autocov(x, t) = transpose(autocov(x, -t))
#         Σ[(r*J+1):(r*J+J), (c*J+1):(c*J+J)] = (r>=c) ? autocov(X, G[r+1]-G[c+1], δ) : transpose(Σ[(c*J+1):(c*J+J), (r*J+1):(r*J+J)])
#     end

#     return Matrix(Symmetric(Σ))  # forcing symmetry
# end

# function covmat(X::FractionalWaveletNoiseBank, G::AbstractVector{<:Integer})
#     J = length(X.filters)
#     l = length(G)
#     Σ = zeros(l*J, l*J)

#     Σs = [autocov(X, G[n]-G[1], 1) for n=1:length(G)]

#     for c=0:l-1,r=0:l-1
#         Σ[(r*J+1):(r*J+J), (c*J+1):(c*J+J)] = (r>=c) ? Σs[r-c+1] : transpose(Σs[c-r+1])
#     end

#     return Matrix(Symmetric(Σ))  # forcing symmetry
# end

# covmat(X::FractionalWaveletNoiseBank, l::Integer) = covmat(X, 1:l)

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
@everywhere begin
    struct DynamicHistogram
        bins::Vector{Float64}
        start::Float64
        binwidth::Float64
        maxbins::Int
        lock::Threads.SpinLock
    end
    DynamicHistogram(start, binwidth, maxbins) =
        DynamicHistogram(zeros(1), start, binwidth, maxbins, Threads.SpinLock())
    function Base.:+(hist1::T, hist2::T) where {T<:DynamicHistogram}
        bins = zeros(max(length(hist1.bins), length(hist2.bins)))
        for i in 1:length(bins)
            bins[i] += length(hist1.bins) >= i ? hist1.bins[i] : 0.0
            bins[i] += length(hist2.bins) >= i ? hist2.bins[i] : 0.0
        end
        DynamicHistogram(bins, hist1.start, hist1.binwidth, hist1.maxbins, Threads.SpinLock())
    end
    function Base.push!(hist::DynamicHistogram, value::Float64)
        idx = 1
        binstart = hist.start
        lock(hist.lock) do
            while length(hist.bins) < hist.maxbins
                # Add bins if necessary
                if idx > length(hist.bins)
                    push!(hist.bins, 0.0)
                end
                if binstart <= value <= binstart+hist.binwidth
                    hist.bins[idx] += value
                    return
                end
                idx += 1
                binstart += hist.binwidth
            end
            throw(DomainError(value, "Exhausted allowed number of bins"))
        end
    end
end

@testset "Mutation" begin

@testset "@mutable" begin
    w = first(workers())
    x = remotecall_fetch(w) do
        Dagger.@mutable Ref{Int}()
    end
    @test fetch(Dagger.@spawn (x->x[] = myid())(x)) == w
    @test_throws_unwrap Dagger.ThunkFailedException fetch(Dagger.@spawn single=last(workers()) (x->x[] = myid())(x))
end # @testset "@mutable"

@testset "Shard" begin
    cs = Dagger.@shard Threads.Atomic{Int}(0)
    s = fetch(cs; raw=true)
    ctxprocs = Dagger.Sch.eager_context().procs
    for p in keys(s.chunks)
        @test p isa OSProc
        @test p in ctxprocs

        c = s.chunks[p]
        @test c.processor == p
        @test c.scope isa Dagger.ProcessScope
        @test c.scope.wid == p.pid

        @test fetch(c) isa Threads.Atomic{Int}
    end

    @testset "procs kwarg" begin
        procs = [OSProc(first(workers()))]
        cs = Dagger.@shard procs=procs Threads.Atomic{Int}(0)
        s = fetch(cs; raw=true)
        @test length(keys(s.chunks)) == 1
        p = first(keys(s.chunks))
        @test p isa Dagger.OSProc
        @test p.pid == first(workers())

        c = s.chunks[p]
        @test c.processor == p
        @test c.scope isa Dagger.ProcessScope
        @test c.scope.wid == p.pid
    end

    @testset "workers kwarg" begin
        cs = Dagger.@shard workers=[first(workers())] Threads.Atomic{Int}(0)
        s = fetch(cs; raw=true)
        @test length(keys(s.chunks)) == 1
        p = first(keys(s.chunks))
        @test p isa Dagger.OSProc
        @test p.pid == first(workers())

        c = s.chunks[p]
        @test c.processor == p
        @test c.scope isa Dagger.ProcessScope
        @test c.scope.wid == p.pid
    end

    @testset "per_thread kwarg" begin
        cs = Dagger.@shard per_thread=true Threads.Atomic{Int}(0)
        s = fetch(cs; raw=true)
        for p in keys(s.chunks)
            @test p isa Dagger.ThreadProc
            gp = Dagger.get_parent(p)
            @test gp in ctxprocs
            @test p in Dagger.get_processors(p)

            c = s.chunks[p]
            @test c.processor == p
            @test c.scope isa Dagger.ExactScope
            @test c.scope.processor == p
        end
    end

    # Can't mix procs and workers
    @test_throws ArgumentError Dagger.@shard procs=[OSProc] workers=[1] 1+1
    # Can't mix procs and per_thread=true
    @test_throws ArgumentError Dagger.@shard procs=[OSProc] per_thread=true 1+1

    # Need at least one processor
    @test_throws ArgumentError Dagger.@shard procs=[] 1+1
    @test_throws ArgumentError Dagger.@shard workers=[] 1+1

    @testset "Distributed Histogram" begin
        cs = Dagger.@shard DynamicHistogram(0.0, 1.0, 100)
        wait.([Dagger.@spawn push!(cs, 5.2) for i in 1:100])
        @test sum(sum(fetch.(map(identity, cs))).bins) ≈ 5.2*100
    end
end # @testset "Shard"

end # @testset "Mutation"

ext: jl


text:
# required by item/Calibration/mmle.jl
include("compute_likelihood.jl")
include("cubic_spline.jl")
include("rescale.jl")
include("convergence.jl")


text:
abstract type Mortality{R,W} <: CellRule{R,W} end

struct LinearMortality{R,W,D,T} <: Mortality{R,W}
    proportion::D
    threshold::T
end
LinearMortality{R,W}(d::D; threshold::T = 0.01) where {R,W,D,T} =
    LinearMortality{R,W}(d, threshold)

function DynamicGrids.applyrule(data, rule::LinearMortality, N, I)
    # potentially type asset N is occupancy, later
    N > zero(N) || return zero(N)
    N > rule.threshold || return zero(N)
    new = @fastmath N - rule.proportion * N
    new > zero(N) && return new
    return zero(N)
end

LinearMortality(layernames::T, decay::DT; threshold::TT = 0.01) where {T,DT,TT} = begin
    rules = Ruleset()
    for sym in layernames
        rules += LinearMortality{sym}(decay, threshold = threshold)
    end
    return rules
end


text:
# -*- coding: utf-8 -*-

using Base.Test

include("./wx_tests.jl")
include("./test_tk.jl")
w = test_tk("title日本語漢字表示申能")
sleep(1.5)
exists(w) && destroy(w)


text:
# module comparison

path,_ = splitdir(@__FILE__)
# using PyPlot
#Packages I haven't developed
import Interpolations
using PyPlot
close("all")
rc("figure", figsize=(3.5, 2.6))
rc("font", size=10.0)
rc("lines", linewidth=1.5)
rc("lines", markersize=3.0)
rc("legend", frameon=false)
rc("axes.spines", right=false, top=false)
rc("figure.subplot", left=0.19, bottom=0.18, top=0.97, right=0.95)
rc("axes", color_cycle=["348ABD", "A60628", "009E73", "7A68A6", "D55E00", "CC79A7"])
# plot_cycle=["#348ABD", "#A60628", "#009E73", "#7A68A6", "#D55E00", "#CC79A7"]
import JLD
import CSV
import Dierckx
import StatsBase

# Packages I played a very minor role
import PreComp
import Composites
import VTKtools

# Packages I played a minor role
import OptParams
import Atmosphere
import CCBlade #new CCBlade package w/o ND splined af data
import BeamFEA
import Snopt
import Gradients

# Packages I played a major role in
import BrentMin
import MotorPower
import Akima
# import VLM
include("oldvlm.jl")
import BPM
import AirfoilPrep

printiter = 1

include("$path/objective.jl")
include("$path/propstructmisc.jl")
include("$path/propwingmisc.jl")
include("$path/setup.jl")
include("$path/save_vtk.jl")

plyprops = plyproperties()

ctlparams = controlparams(
objective=:none, #:eta_total_cruise #:mass,
printfreq = 50,
savevtk = 500,
printobjective=true,
printdetailedoutput=true,
printconstraintviolations=true,
)

# Define general constraints
constraints = [
:p_materialfailure,
:p_localbuckling,
:w_materialfailure,
:w_localbuckling,

:totalmass,
:machtip,
:Lift,
:Drag,
:noise,
# :p_chordlamthick,
# :range,
:takeoff_dist,
:stall,
:capacity,
]

# Prop inputs
p_n_akima = 6
Dprop = 0.406
# p_radii =  1e-3*[17.5, 30, 45, 60, 75, 90, 105]/(Dprop/2)
p_radii = [0.205, .231, .400, .511, .615, .720, .789, .860, .930, .965, .999, 1.0]
p_chord_pts = Dprop/2*[1.56e-01, 1.56e-01, 1.46e-01, 1.47e-01, 1.52e-01, 1.54e-01, 1.52e-01, 1.38e-01, 1.10e-01, 9.13e-02, 6.87e-02, 6.0e-02]
p_twist_pts = pi/180*[5.27e+01, 5.27e+01, 4.39e+01, 3.85e+01, 3.35e+01, 2.91e+01, 2.69e+01, 2.47e+01, 2.35e+01, 2.25e+01, 2.19e+01, 2.09e+01] #+ p_pitch
J = 0.695  # prop advance ratio
velocity = 19.0
p_rpm = velocity/(J*Dprop)*60.0  # rotation rate
etaprop_center = [-0.35,0.35]

#check thrust coefficient make sure it is 0.168
twist_75 = interp1(p_radii, p_twist_pts, [0.75])[1]
println(twist_75*180/pi)
p_pitch = 30*pi/180 - twist_75
# println(p_pitch*180/pi)
p_twist_pts += p_pitch - 2.0*pi/180 +2.0*pi/180

#Wing Inputs
w_n_linear = 30
w_nondim_halfspany = [0.0,0.35,1.0]
w_halfspan = 2.58/2
b = w_halfspan*2
w_chord =[0.279, 0.279, 0.161]
w_sweep = [0.0,-5,0.0]*pi/180
w_aoa = 4.0*pi/180
cmac = 0.24

table_af_wing_blown = JLD.load("$path/airfoils/af_prop_E212.jld")
table_3D_wing_blown = table_af_wing_blown["NDtable"]
spline_3D_blown = AirfoilPrep.SplineND_from_tableND(table_3D_wing_blown)

table_af_wing_free = JLD.load("$path/airfoils/af_prop_E212.jld")
table_3D_wing_free = table_af_wing_free["NDtable"]
spline_3D_freestream = AirfoilPrep.SplineND_from_tableND(table_3D_wing_free)

# # Load the prop airfoil ND table
# E212_non_extrap = JLD.load("$(path)/airfoils/af_prop_E212.jld")
# E212_non_extrap = E212_non_extrap["NDtable"]
# omega = p_rpm*pi/30
# TSR = omega * Dprop/2 / velocity
#
# grid_alphas=[i for i in -180:1.0:180]
# Rtip = Dprop/2
# chord_prop_pt = p_chord_pts
# rpt = p_radii * Rtip
# r = linspace(p_radii[2],p_radii[end-1],p_n_akima) * Rtip
# chord_prop = Akima.interp_noderiv(rpt, chord_prop_pt, r)
# af = Array{Any}(p_n_akima)
# for i = 1:p_n_akima
#
#    r_over_R = r[i]
#    c_over_r = chord_prop[i]/r[i]
#    NDextrap3D_3Dtable = AirfoilPrep.NDTable_correction3D_extrap(E212_non_extrap,r_over_R,c_over_r,TSR;grid_alphas=grid_alphas)
#    afspl = AirfoilPrep.SplineND_from_tableND(NDextrap3D_3Dtable)
#    af[i] = afspl
# end
# JLD.save("$(path)/airfoils/af_epema.jld", "af", af)

af = JLD.load("$(path)/airfoils/af_epema.jld")
af = af["af"]
af = fill(af[1],p_n_akima)

af_xy = CSV.read("$(path)/airfoils/e212-il.csv",header = true, delim = ',',datarow = 10,allowmissing=:none)
xaf1 = reverse(af_xy[:,1])/100.0
yaf1 = reverse(af_xy[:,2])/100.0

# nameout = "$(path)/AirfoilData/e212.dat"
# write(nameout,"x\ty\n")
# open(nameout,"a") do x
#     for i = 1:length(xaf1)
#         write(x,"$(xaf1[i])\t$(yaf1[i])\n")
#     end
# end

#Filter airfoil to get points where strain is evaluated
n_strainpts = 10
xafstrain = zeros(n_strainpts)
yafstrain = zeros(n_strainpts)
step_simplify = round(Int,length(xaf1)/n_strainpts)
idx = 1
for i = 1:n_strainpts
    xafstrain[i] = xaf1[idx[1]]
    yafstrain[i] = yaf1[idx[1]]
    idx += step_simplify
end

# Fill airfoil data for each section
p_xaf = zeros(p_n_akima,length(xaf1))
p_yaf = zeros(p_n_akima,length(yaf1))
p_xafstrain = zeros(p_n_akima,length(xafstrain))
p_yafstrain = zeros(p_n_akima,length(yafstrain))
for i = 1:p_n_akima
    p_xaf[i,:] = copy(xaf1) #break links since using in multiple areas
    p_yaf[i,:] = copy(yaf1)
    p_xafstrain[i,:] = copy(xafstrain)
    p_yafstrain[i,:] = copy(yafstrain)
end

#TODO: different airfoil for wing
w_xaf = zeros(w_n_linear,length(xaf1))
w_yaf = zeros(w_n_linear,length(yaf1))
w_xafstrain = zeros(w_n_linear,length(xafstrain))
w_yafstrain = zeros(w_n_linear,length(yafstrain))
for i = 1:w_n_linear
    w_xaf[i,:] = copy(xaf1)
    w_yaf[i,:] = copy(yaf1)
    w_xafstrain[i,:] = copy(xafstrain)
    w_yafstrain[i,:] = copy(yafstrain)
end

altitude_cruise = 3000.0
designE = propwingdesign(;
    #Design Variables at this scope
    batt_mass = 0.0,
    w_aoa = w_aoa,
    p_rpm =p_rpm,
    p_pitch = 1.0,
    velocity = velocity,

    #Design Parameters at this scope
    constraints = constraints,
    af = af,
    p_xaf = p_xaf,
    p_yaf = p_yaf,
    p_xafstrain = p_xafstrain,
    p_yafstrain = p_yafstrain,
    w_xaf = w_xaf,
    w_yaf = w_yaf,
    w_xafstrain = w_xafstrain,
    w_yafstrain = w_yafstrain,
    spline_3D_freestream = spline_3D_freestream,
    spline_3D_blown = spline_3D_blown,

    # Design variables unique to the system function
    p_lam_tin = [ones(p_radii)*0.001*10;ones(p_radii)*0.001*85],
    p_orientation = [0.0,45.0],
    p_chord_pts = p_chord_pts,
    p_twist_pts = p_twist_pts,
    p_Rtip = Dprop/2,
    kv = 100.0,
    i0 = 2.0,
    w_chord_pts = w_chord,
    w_twist_pts = zeros(w_chord)+w_aoa,
    w_halfspan = w_halfspan,
    w_sweep_pts = w_sweep,
    w_dihedral_pts = zeros(w_chord),
    w_lam_tin = [
        ones(w_nondim_halfspany)*0.001*20;#t1 highmodulus_uni
        ones(w_nondim_halfspany)*0.001*20;#t2 highmodulus_weave
        ones(w_nondim_halfspany)*0.01*2;#t2 highmodulus_weave
        ones(w_nondim_halfspany)*0.001*20;#t2 highmodulus_weave
        ones(w_nondim_halfspany)*0.01*2;#t2 highmodulus_weave
        ],
    w_orientation = [0.0,45.0,45.0],
    w_airfoilthickness = zeros(w_chord), #Does nothing

    #Design variables unique to the system function
    etaprop_center = etaprop_center,
    p_n_akima = p_n_akima,
    p_radii =  p_radii,
    altitude = 0.0,
    numprops = 2,
    p_x_offset = ones(p_radii)*0.25,
    p_aerocenter = ones(p_radii)*0.275,
    p_usedmaterials = ["highmodulus_uni","highmodulus_weave"],
    p_webloc =[0.999],
    plyprops = plyprops,
    blades = 6,
    w_n_linear = w_n_linear,
    w_x_offset = ones(w_nondim_halfspany)*0.25,
    w_aerocenter = ones(w_nondim_halfspany)*0.27,
    w_nondim_halfspany = w_nondim_halfspany,
    w_usedmaterials = ["highmodulus_uni","highmodulus_weave","taylor_foam","highmodulus_weave","taylor_foam"],
    w_webloc = [0.5],
    prop_tilt = 0.0*pi/180,

    #Control parameters
    detailed = true,
    VTKfilename = "Epema.vtk",
    savefinalvtk = true,
    verification = false,
    plots = false,
    n_CP = 200,

)

#Call the PropWing System
detailed_outputE = system(ctlparams,designE)

println(detailed_outputE.CT)

exp_x = [4.61E-02, 1.39E-01, 2.44E-01, 2.87E-01, 4.00E-01, 4.38E-01, 5.62E-01, 6.87E-01, 8.09E-01]
# exp_x2 = [4.57E-02, 1.39E-01, 2.44E-01, 2.87E-01, 3.99E-01, 4.38E-01, 5.62E-01, 6.86E-01, 8.09E-01]
# exp_x3 = [4.61E-02, 1.39E-01, 2.44E-01, 2.88E-01, 4.00E-01, 4.38E-01, 5.62E-01, 6.86E-01, 8.08E-01]
exp_y1 = [5.12E-01, 5.79E-01, 7.54E-01, 6.86E-01, 3.31E-01, -5.93E-02, 3.20E-01, 3.15E-01, 2.83E-01]
exp_y2 = [4.97E-01, 5.63E-01, 7.09E-01, 5.03E-01, 1.44E-01, -1.01E-01, 3.06E-01, 3.03E-01, 2.70E-01]
exp_y3 = [5.27E-01, 5.97E-01, 7.98E-01, 8.70E-01, 5.14E-01, -1.47E-02, 3.31E-01, 3.29E-01, 2.95E-01]

other_vlm_x = [9.66E-03, 2.95E-02, 4.82E-02, 6.69E-02, 8.47E-02, 1.04E-01, 1.23E-01, 1.42E-01, 1.61E-01, 1.81E-01, 2.00E-01, 2.19E-01, 2.38E-01, 2.56E-01, 2.75E-01, 2.95E-01, 3.14E-01, 3.33E-01, 3.61E-01, 3.97E-01, 4.34E-01, 4.70E-01, 5.08E-01, 5.43E-01, 5.80E-01, 6.16E-01, 6.54E-01, 6.89E-01, 7.26E-01, 7.62E-01, 7.99E-01, 8.35E-01, 8.72E-01, 9.08E-01, 9.45E-01, 9.81E-01]
other_vlm_y = [5.18E-01, 5.20E-01, 5.16E-01, 5.19E-01, 5.22E-01, 5.30E-01, 5.41E-01, 5.55E-01, 5.73E-01, 6.02E-01, 6.75E-01, 8.13E-01, 8.74E-01, 8.83E-01, 8.44E-01, 7.54E-01, 5.29E-01, 4.39E-01, 4.19E-01, 1.67E-01, 1.12E-01, 1.63E-01, 2.70E-01, 2.95E-01, 3.04E-01, 3.06E-01, 3.03E-01, 2.97E-01, 2.89E-01, 2.78E-01, 2.65E-01, 2.50E-01, 2.32E-01, 2.09E-01, 1.78E-01, 1.29E-01]

# chord = [p.chord for p in detailed_outputE.panels]
# N = length(detailed_outputE.panels)
# CPy = zeros(N)
# for i = 1:N  # CP
#     rcp = VLM.control_point(detailed_outputE.panels[i])
#     CPy[i] = rcp[2]
# end

chord = detailed_outputE.panels.chord
CPy = detailed_outputE.panels.y

figname = "Epema"
figure(figname)
plot(CPy/designE.w_halfspan, detailed_outputE.cl.*chord/cmac)
errorbar(exp_x, exp_y1, yerr=[exp_y1-exp_y2, exp_y3-exp_y1], fmt="o")
plot(other_vlm_x, other_vlm_y, "k--")
xlim([0 , 1])
xlabel("Halfspan Location")
ylabel("Local Lift Coefficient")#L"c_l\ c / c_{mac}")
legend(["Our VLM", "Epema VLM", "Epema Exper."])
savefig("../MDAO_paper/figures/$figname.pdf", transparent=true)

# Error Calculation

spl_VLM0deg = Dierckx.Spline1D(CPy/designE.w_halfspan, detailed_outputE.cl.*chord/cmac)

error0 = 0.0
for i = 1:length(exp_x)
    error0 = error0+(spl_VLM0deg(exp_x[i])-exp_y1[i]).^2
end
sumofsquareserror0 = sqrt(error0/length(exp_x))
println("sumofsquareserror0 $sumofsquareserror0")


text:
include("0_setup.jl");

R"""
library(data.table)
library(magrittr)

system.time(lc <- dir($lc_csv_path, full.names = T) %>%
  sapply(dir, full.names = T) %>%
  lapply(fread, na.strings = NULL) %>%
  rbindlist)

feather::write_feather(lc, $lc_feather_path)
"""


text:
using TrussMorph
tm = TrussMorph
using Test
using Makie
using Colors
using Dates
using Printf

recompute = true
plot = true
node_dof = 3 # 2 for truss model, 3 for frames
full_node_dof = 3 # fixed for 2D cases
box_constr = undef
sym_var_ids = undef

dir = "/Users/yijiangh/Dropbox (MIT)/Course_Work/6.838_2019_spring/final_project/code/gh_validation"
# result_file_dir = joinpath(pwd(),"test","results")
result_file_dir = "/Users/yijiangh/Dropbox (MIT)/Course_Work/6.838_2019_spring/final_project/code/results"

# st_file_name  = "wiggle_truss_min_d.json"
# end_file_name  = "wiggle_truss_min_w.json"
# fp0 = joinpath(dir, st_file_name)
# fp1 = joinpath(dir, end_file_name)
# load_fp = joinpath(dir, "wiggle_truss_min_d_load_case.json")
# full_design_var_ids = vcat(collect(3:12), collect(15:24))
# full_design_var_ids = reshape(full_design_var_ids, (2, Int.(length(full_design_var_ids) / 2)))'
# design_var_ids = full_design_var_ids[:,2]

# st_file_name  = "funicular_arch.json"
# end_file_name  = "funicular_cable.json"
# fp0 = joinpath(dir, st_file_name)
# fp1 = joinpath(dir, end_file_name)
# load_fp = joinpath(dir, "funicular_arch_load_case.json")
# design_var_ids = collect(4:2:20)

st_file_name  = "2D_truss_4.json"
end_file_name  = "2D_truss_5.json"
# st_file_name  = "2D_truss_0.json"
# end_file_name  = "2D_truss_1.json"
fp0 = joinpath(dir, st_file_name)
fp1 = joinpath(dir, end_file_name)
load_fp = joinpath(dir, "2D_truss_0_load_case.json")
design_var_ids = [4, 6]
# design_var_ids = [3, 4, 5, 6]

# st_file_name  = "cm_truss_0.json"
# end_file_name  = "cm_truss_1.json"
# fp0 = joinpath(dir, st_file_name)
# fp1 = joinpath(dir, end_file_name)
# # load_fp = joinpath(dir, "cm_truss_2pts_uneven_load_case.json")
# # load_fp = joinpath(dir, "cm_truss_0_load_case.json")
# load_fp = joinpath(dir, "cm_truss_2pts_load_case.json")
# # design_var_ids = [4, 6, 8]
# design_var_ids = [3, 4, 6, 7, 8]
# # box_constr = [1.0 3; -5 5; -5 5; 3.0 5; -5 5]

path_disc = 20
parm_weight = 1.0
parm_smooth = 1e4
# parm_smooth = 1e4

# plot parameters
load_scale = 0.01
supp_scale = 0.1
line_width = 4.0

t0,_ = parse_truss_json(fp0)
t1,_ = parse_truss_json(fp1)
load = parse_load_json(load_fp, node_dof)

st_time = time()
morph_path, init_morph_path, opt_sE, opt_wE, opt_totE, init_sE, init_wE, init_totE, opt_sumE, init_sumE, parm =
    tm.compute_morph_path(t0, t1, load, node_dof, full_node_dof, design_var_ids, path_disc=path_disc,parm_smooth=parm_smooth, parm_weight=parm_weight, do_opt=recompute, box_constraint=box_constr)
@show time() - st_time

pure_st_file_name = SubString(st_file_name, 1:length(st_file_name)-length(".json"))
pure_end_file_name = SubString(end_file_name, 1:length(end_file_name)-length(".json"))
result_file_name = pure_st_file_name * "-" *  pure_end_file_name
f_file_dir = joinpath(result_file_dir, result_file_name)

if !recompute
    searchdir(path,key) = filter(x->occursin(key,x), readdir(path))

    opt_sE = Float64[]
    opt_wE = Float64[]
    opt_totE = Float64[]

    @assert(ispath(f_file_dir))
    result_jsons = searchdir(f_file_dir,string("-", path_disc+2, ".json"))
    # @show result_jsons
    morph_path = Matrix{Float64}[]
    morph_data = Dict()

    for json_fp in result_jsons
        full_json_fp = joinpath(f_file_dir, json_fp)

        parsed_truss, morph_data = parse_truss_json(full_json_fp, parse_morph=true)
        push!(morph_path, parsed_truss.X)
        push!(opt_wE, morph_data["weight_energy"])
        push!(opt_sE, morph_data["smoothness_energy"])
        push!(opt_totE, morph_data["tot_energy"])
        parm = [morph_data["smooth_parameter"], morph_data["weight_parameter"]]
        opt_sumE = morph_data["sum_energy"]
    end
    # @show opt_sumE
    # @show parsed_morph_path
else
    tm.save_morph_path_json(morph_path, result_file_dir, st_file_name, end_file_name, t0, opt_sE, opt_wE, opt_totE, opt_sumE, parm_smooth, parm_weight)
end

if plot
    plen = length(morph_path)
    r = ones(Float64, size(t0.T,1)) * 0.5
    color_array = Array{RGBAf0}(undef, plen)

    max_x = maximum([maximum(M[:,1]) for M in morph_path])
    min_x = minimum([minimum(M[:,1]) for M in morph_path])
    max_y = maximum([maximum(M[:,2]) for M in morph_path])
    min_y = minimum([minimum(M[:,2]) for M in morph_path])
    max_xlim = max_x - min_x
    max_ylim = max_y - min_y
    max_lim = max(max_xlim, max_ylim)
    limits = FRect(min_x, min_y,
                   max_xlim, max_ylim)

    scene = Scene()
    init_scene= Scene()
    tm.draw_load!(scene, t0.X, load, load_scale=max_ylim*load_scale, xaxis_label="x - OPT", plot_limits=limits)
    tm.draw_load!(init_scene, t0.X, load, load_scale=max_ylim*load_scale, xaxis_label="x - INIT(linear)", plot_limits=limits)

    for i=1:plen
        mcolor = Float32.((plen - i) / plen) * RGBAf0(1.0,0.0,0.0,0.1) + Float32.(i / plen) * RGBAf0(0.0,0.0,1.0,1)
        color_array[i] = mcolor

        tm.draw_truss!(scene, morph_path[i], t0.T, t0.S, r, supp_scale=supp_scale, color=mcolor, xaxis_label="x - OPT", line_width=line_width, plot_limits=limits)

        tm.draw_truss!(init_scene, init_morph_path[i], t0.T, t0.S, r, supp_scale=supp_scale, color=mcolor, xaxis_label="x - INIT(linear)", line_width=line_width, plot_limits=limits)
        # display(scene)
        # sleep(0.2)
    end

    sE_limits = FRect(1, 0, plen, max(maximum(opt_sE), maximum(init_sE)))
    wE_limits = FRect(1, 0, plen, max(maximum(opt_wE), maximum(init_wE)))
    totE_limits = FRect(1, 0, plen, max(maximum(opt_totE), maximum(init_totE)))

    scene_sE = barplot(opt_sE, color=color_array,
                       axis = (names = (axisnames = ("iter - OPT", "smoothness (neighbor) energy"),),
                              ),
                       limits=sE_limits,
                       )
    scene_wE = barplot(opt_wE, color=color_array,
                      axis = (names = (axisnames = ("iter - OPT", "weight energy"),),
                             ),
                      limits=wE_limits,
                      )
    scene_totalE = barplot(opt_totE, color=color_array,
                     axis = (
                        names = (axisnames = (string("iter - OPT total E: ", opt_sumE), string("total energy- pS:", parm[1], ", pW:", parm[2])),),
                            ),
                     limits=totE_limits,
                     )

    init_scene_sE = barplot(init_sE, color=color_array,
                       axis = (
                            names = (axisnames = ("iter INIT", "smoothness (neighbor) energy"),),
                              ),
                       limits=sE_limits,
                       )
    init_scene_wE = barplot(init_wE, color=color_array,
                      axis = (names = (axisnames = ("iter INIT", "weight energy"),),
                             ),
                      limits=wE_limits,
                      )
    init_scene_totalE = barplot(init_totE, color=color_array,
                     axis = (
                        names = (axisnames = (string("iter INIT -total E: ", init_sumE), string("total energy- pS:", parm[1], ", pW:", parm[2])),),
                            ),
                      limits=totE_limits,
                     )

    scene_final = hbox(
        vbox(scene_sE, scene_wE, scene_totalE),
        vbox(init_scene_sE, init_scene_wE, init_scene_totalE),
        vbox(scene, init_scene)
        )

    # scene_final = scene_sE

    Makie.save(joinpath(f_file_dir, result_file_name * "_" * string(Dates.now()) * ".png"), scene_final)

    # anim_sc = Scene()
    # display(anim_sc)

    # save morph gif
    # record(anim_sc, joinpath(f_file_dir, result_file_name * "_" * string(Dates.now()) * ".gif"), 1:plen) do i
    #     tm.draw_truss!(anim_sc, morph_path[i], t0.T, t0.S, r, supp_scale=supp_scale, color=color_array[i], xaxis_label="x - OPT", line_width)
    #     display(anim_sc)
    #     sleep(1)
    # end
    # t = Node(i)
    # color=color_array[lift(i->i,t)]

    # record(anim_sc, joinpath(f_file_dir, result_file_name * "_" * string(Dates.now()) * ".gif"), 1:plen) do i

    for i=1:plen
        anim_sc = Scene()
        # display(anim_sc)
        mcolor = Float32.((plen - i) / plen) * RGBAf0(1.0,0.0,0.0,1) + Float32.(i / plen) * RGBAf0(0.0,0.0,1.0,1)

        tm.draw_load!(anim_sc, morph_path[i], load, load_scale=max_ylim*load_scale, xaxis_label="x - OPT", plot_limits=limits)
        tm.draw_truss!(anim_sc, morph_path[i], t0.T, t0.S, r, supp_scale=supp_scale, xaxis_label="x - OPT", color = mcolor, line_width=line_width, plot_limits=limits)

        result_img_name = result_file_name * string("_", @sprintf("%03d",i),"-",plen) * ".png"
        Makie.save(joinpath(f_file_dir, result_img_name), anim_sc)
        display(anim_sc)
        sleep(0.1)
    end

    display(scene_final)
end


text:
const _UniformEncodingTables = Dict()
"""
    uniform_encoding_table(T, X::AbstractUnitRange, degree::Int)

Get the encoding table of `X` using uniform degree filter.

A uniform degree filter is defined as `f(x) = U(x) > degree ? P+1 : x`, where `P` is number
of bits of `x`, and `U(x)` is the circular bit transition of `x` with bit representation of
datatype `T`.

For example, `0b0000001` has ``2`` bit transitions thus is unchanged, while `0b00000101` has
``4`` bit transitions. The transition is count in circular sense, i.e., the highest and
lowest bits are also compared.

```jldoctest
julia> using LocalBinaryPatterns: uniform_encoding_table

julia> X = uniform_encoding_table(UInt8, 0:255, 2);

julia> X[0b0000001]
0x01

julia> X[0b0000101] # miscellaneous
0x09
```

This function is used to distinguish local binary patterns from texture patterns and
miscellaneous patterns. See [1] for more information.

!!! info "Runtime caches"
    For better performance, this function uses a runtime cache to store the lookup and
    shared among all callers. The result is expected to be used in read-only mode.

# References

- [1] T. Ojala, M. Pietikäinen, and T. Mäenpää, “A Generalized Local Binary Pattern Operator for Multiresolution Gray Scale and Rotation Invariant Texture Classification,” in _Advances in Pattern Recognition — ICAPR 2001, vol. 2013, S. Singh, N. Murshed, and W. Kropatsch, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg_, 2001, pp. 399–408. doi: 10.1007/3-540-44732-6_41.
"""
function uniform_encoding_table(::Type{T}, X::AbstractUnitRange, degree::Int) where T
    d = _UniformEncodingTables
    k = (T, degree, minimum(X), maximum(X))
    haskey(d, k) && return d[k]
    lookup = d[k] = _uniform_encoding_table(T, X, degree)
    return lookup
end
uniform_encoding_table(::Type{T}, nbits::Int, degree::Int) where T = uniform_encoding_table(T, 0:2^nbits-1, degree)

function _uniform_encoding_table(::Type{T}, X::AbstractUnitRange, degree::Int) where T<:Unsigned
    if !((maximum(X) <= typemax(T)) && (minimum(X) >= typemin(T)))
        throw(ArgumentError("range $(first(X)):$(last(X)) can't be uniquely represented by type $T."))
    end

    function _count_bit_transitions(x)
        # Eq. (10) in Ojala 2002
        count_ones(x << 7 & typemax(typeof(x))) + count_ones(x ⊻ (x >> 1))
    end
    # Eq. (9) for Ojala 2002
    lookup = map(X) do x
        n = _count_bit_transitions(T(x))
        T(ifelse(n > degree, 8sizeof(T)+1, x))
    end
    OffsetVector(lookup, X)
end


text:
#=##############################################################################
# DESCRIPTION
    Functions common to all grids that have the properties `_ndivsnodes`,
    `_ndivsnodes`, and `_override_vtkcelltype`.
# AUTHORSHIP
  * Author    : Eduardo J Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Jun 2018
  * License   : MIT License
=###############################################################################

get_ndivsnodes(self::GridExtentions) = self._ndivsnodes
get_ndivscells(self::GridExtentions) = self._ndivscells

function get_cellnodes(self::GridExtentions, i::Int64)
  return [get_node(self, n) for n in get_cell(self, i)]
end
function get_cellnodes(self::GridExtentions, coor::Array{Int64,1})
  return get_cellnodes(self, sub2ind(self._ndivscells, coor...))
end

"""
  `get_fieldval(grid, field_name, coor)`

Returns the value of node of coordinates `coor` (1-indexed) in the field
'field_name'.
"""
function get_fieldval(self::GridExtentions, field_name::String, coor::Array{Int64,1})

  if self.field[field_name]["entry_type"]=="node"
    return get_fieldval(self, field_name, Base._sub2ind(self._ndivsnodes, coor...)) #TODO: use LinearIndex instead of _sub2ind() (next line too)

  elseif self.field[field_name]["entry_type"]=="cell"
    return get_fieldval(self, field_name, Base._sub2ind(self._ndivscells, coor...))

  else
    error("Entry type $(self.field[field_name]["entry_type"]) not implemented.")
  end
end

"""
  `get_fieldval(grid, field_name, i)`

Returns the value of i-th node (1-indexed) in the field 'field_name'.
"""
function get_fieldval(self::GridExtentions, field_name::String, i::Int64)
  if !(field_name in keys(self.field))
    error("Field $field_name doesn't exist."*
          " Available fields: $(keys(self.field))")
  end

  return self.field[field_name]["field_data"][i]
end

"""
  `add_field(grid::Grid, field_name::String, field_type::String, field_data,
entry_type::String)`

Adds a field of data associated to each node.

NOTE: each data entry must be a single value if `field_type==scalar`, or a
      3-element array if `field_type==vector`.
"""
function add_field(grid::GridExtentions, field_name::String,
                            field_type::String, field_data, entry_type::String;
                            raise_warn::Bool=true)
  # Error cases
  if !(field_type in ["vector", "scalar"])
    error("Unkown field type $(field_type)")
  elseif !(entry_type in ["node", "cell", "system"])
    error("Unkown entry type $(entry_type)")
  elseif entry_type=="node" && size(field_data, 1)!=grid.nnodes
    error("Invalid node field size."*
            " Expected $(grid.nnodes), got $(size(field_data, 1))")
  elseif entry_type=="cell" && size(field_data, 1)!=grid.ncells
    error("Invalid cell field size."*
            " Expected $(grid.nnodes), got $(size(field_data, 1))")
  end

  if field_name in keys(grid.field) && raise_warn
    @warn("Overwritting field $field_name.")
  end

  grid.field[field_name] = Dict("field_name" => field_name,
                                "field_type" => field_type,
                                "entry_type" => entry_type,
                                "field_data" => field_data)
end

"""
  `calculate_field(grid::Grid, f, field_name::String, field_type::String)`

Evaluates the function `f` at each nodes and stores the values as a new field.

NOTE: f must return a single value if `field_type==scalar`, or a 3-element array
      if `field_type==vector`.
"""
function calculate_field(grid::GridExtentions, f, field_name::String,
                                        field_type::String, entry_type::String)
  if entry_type=="node"
    field_data = [f(get_node(grid, i)) for i in 1:grid.nnodes]
  elseif entry_type=="cell"
    field_data = [f(get_cell(grid, i)) for i in 1:grid.ncells]
  else
    error("Unkown entry type $(entry_type)")
  end
  add_field(grid, field_name, field_type, field_data, entry_type)
end

"""
  `save(grid, filename; args...)`

Outputs a vtk file of this grid. See generateVTK for a descrition of optional
arguments `args...`.
"""
function save(grid::GridExtentions, filename::String; O=nothing, Oaxis=nothing,
                                                                        args...)
  # Determines whether to add 0 to points for vtk file
  if grid.dims<=3
    aux1 = zeros(Float64, 3-grid.dims)
  else
    error("$(grid.dims)-dimensional grids can't be exported as VTKs!")
  end

  if O != nothing || Oaxis != nothing

    if grid.dims != 3
      error("Requested space transformation on grid of $(grid.dims) "*
            "dimensions. There is not such implementation yet!")
    end
    O = O==nothing ? zeros(3) : O
    Oaxis = Oaxis==nothing ? Array(1.0I, 3, 3) : Oaxis
    check_coord_sys(Oaxis)
    invOaxis = collect(Oaxis')

    points = [countertransform(get_node(grid, i), invOaxis, O) for i in 1:grid.nnodes]

  else

    points = [vcat(get_node(grid, i), aux1) for i in 1:grid.nnodes]

  end

  cells = [get_cell(grid, i) .- 1 for i in 1:grid.ncells]

  point_data = []
  cell_data = []
  for (key, val) in grid.field
    if val["entry_type"] == "node"
      push!(point_data, val)
    elseif val["entry_type"] == "cell"
      push!(cell_data, val)
    elseif val["entry_type"] == "system"
      nothing
    else
      error("Unkown entry type $(entry_type)")
    end
  end

  # Checks for quasi-dimensions (NDIVS[i]==0)
  qdims = length( [1 for ndiv in grid._ndivscells if ndiv==0] )
  dims = grid.dims - qdims

  ctype = grid._override_vtkcelltype

  return generateVTK(filename, points; cells=cells,
                      point_data=size(point_data,1)!=0 ? point_data : nothing,
                      cell_data=size(cell_data,1)!=0 ? cell_data : nothing,
                      _griddims=dims,
                      override_cell_type=ctype,
                      args...)
end


text:
#=
  Package: FishEBM
  File: types.jl
  Justin Angevaare, Devin Rose
  Type definitions for the structured stock-level and agent-level model components
  Created: May 2015
=#

"""
  Description: Assumptions regarding mortality, maturity, broodsize, and
    catchability of age specific adult agents.

  naturalmortality = Age specific survivorship (survivorship at carrying
    capacity if density depedence occurs).
  halfmature = Age at 50% mature (Binomial cdf assumed).
  broodsize = Age specific fecundity (i.e. mean quantity of eggs each spawning
    female will produce).
  fecunditycompensation = compensatory strength for changes in fecundity.
    Compensatory strength is a divisor of K which will result in a 68% change in
    fecundity - smaller values indicate lower compensation strength.
    Compensation function based on Normal CDF. Use NaN if compensation is
    assumed to not occur.
  maturitycompensation = compensatory strength for changes in age of sexual
    maturity. Compensatory strength is a divisor of K which will result in a 68%
    change in age of sexual maturity - smaller values indicate lower
    compensation strength. Compensation function based on Normal CDF. Use NaN if
    compensation is assumed to not occur.
  mortalitycompensation = compensatory strength for changes in age of sexual
    maturity. Compensatory strength is a divisor of K which will result in a 68%
    change in natural mortality - smaller values indicate lower compensation
    strength. Compensation function based on Normal CDF. Use NaN if compensation
    is assumed to not occur.

    Last Update: September 2016
"""
type AdultAssumptions
  naturalmortality::Vector{Float64}
  halfmature::Float64
  broodsize::Vector{Int64}
  fecunditycompensation::Float64
  maturitycompensation::Float64
  mortalitycompensation::Float64
  catchability::Vector{Float64}

  AdultAssumptions() = new()
  AdultAssumptions(naturalmortality, halfmature, broodsize, fecunditycompensation, maturitycompensation, mortalitycompensation, catchability) =
    new(naturalmortality, halfmature, broodsize, fecunditycompensation, maturitycompensation, mortalitycompensation, catchability)
end

"""
  Description: Assumptions regarding mortality, movement, and growth of
    autonomous agents.

    naturalmortality = Stage specific survivorship of each habitat type and stage
      of agents.
    extramortality = Stage specific survivorship of agents (each independent life
      stage) in an environment location influenced by anthropogenic effects.
    growth = Growth rate (in weeks) from one life stage to the next.
    movement = Movement weight matrices for each life stage.
    autonomy = Movement autonomy.

  Last Update: October 2016
"""
type AgentAssumptions
  naturalmortality::Array
  extramortality::Vector{Float64}
  growth::Vector{Int64}
  movement::Array
  autonomy::Vector{Float64}

  AgentAssumptions() = new()
  AgentAssumptions(naturalmortality, extramortality, growth, movement, autonomy) =
    new(naturalmortality, extramortality, growth, movement, autonomy)
end

"""
  Description: This "Environment agent" is used to track fish population
    dynamics during a simulation.

  Last Update: September 2016
"""
type EnviroAgent
   # locationID corresponds to the ID in the habitat array
  locationID::Int64

  # Number of alive agents alive in each cohort
  alive::Vector{Int64}

  # Independent mortality tracking by location
  killedNatural::Vector{Int64}
  killedExtra::Vector{Int64}
  harvest::Int64

  # The spawning week of each cohort during the simulation
  weekNum::Vector{Int64}

  EnviroAgent(locationID) = new(locationID, [0], [0,0,0,0], [0,0,0,0], 0, [0])
end

"""
  Description: A specialized type which contains layers of information to
    indicate spawning area, habitat type, and additional risks. Abstract matches areas of interest
    to habitat index number, and hashing matches index numbers to agent numbers. Hash contains agent
    number during run-time.

  Precondition: Location id should be specified as NaN when a valid location
    does not exist.

  Last Update: October 2016
"""
type EnvironmentAssumptions
  # Spawning map and a hashed list of areas
  spawning::Vector{Int64}
  spawningHash::Vector{Int64}

  # Habitat type array
  habitat::Array{Int64}

  # risk map and a hashed list of areas
  risk::Vector{Int64}
  riskHash::Vector{Int64}

  # Harvest map, hashed list and harvest zones
  harvest::Array{Int64}
  abstractHarvest::Vector{Int64}
  harvestZoneNum::Vector{Int64}
  harvestHash::Vector{Int64}

  # Constructors
  EnvironmentAssumptions() = new()
  EnvironmentAssumptions(spawning, spawningHash, habitat, risk, riskHash, harvest, abstractHarvest, harvestZoneNum, harvestHash) =
  new(spawning, spawningHash, habitat, risk, riskHash, harvest, abstractHarvest, harvestZoneNum, harvestHash)
end


text:
#===============================================================================

    construction.jl - LR, June 2020

    Routines to actually construct the Hamiltonian from a given many-body basis.

===============================================================================#
using SparseArrays, LinearAlgebra
include("typedefs.jl")

function get_particle_hops(state::SpinState, n_basis::Integer)::Array{Tuple{Int,Int,SpinState,Int},1}
    """ Returns possible single particle jumps from one state to another. This
        requires relatively little memory but using this speeds up the computation
        slightly and we don't have to re-compute things later (for additional
        terms, for instance).

        Number of entries per state: n_part*(n_basis-n_npart+1)
        Example: n_basis = 30 // n_part = 16 ---> 240 single particle hops

    """
    hops = Array{Tuple{Int,Int,SpinState,Int},1}()
    for i = 1:n_basis
        temp = annihilate(state, i)
        if !isnothing(temp)
            for j = 1:n_basis
                new_state = create(temp, j)
                if !isnothing(new_state)
                    push!(hops, (i,j,new_state,sign(state,i,j)))
                end
            end
        end
    end
    return hops
end

function find_n_basis(hilbert_space::Array{FullState,1})::Integer
    """ Finds the maximal index of the occupied single-partilce orbitals so that
        the loops are terminated at this level.

        Example:

                states = [|1000000>,|0000100>] -> find_n_basis(states) = 7

        Notes:
            -   Potentially this strategy could be produce a hit on performance,
                but it is also the most general one.
            -   TODO: this could be optimized if the lookup_table is assumed to
                be ordered - then the last entry would automatically give the
                correct answer -> would need OrderedDict. Is this optimization
                necessary?
    """
    n_basis::Integer = 0
    for state in hilbert_space
        s_up, s_down = f_to_s(state)

        for s in [s_up, s_down]
            m = max_orbital(s)
            if m > n_basis
                n_basis = m
            end
        end
    end
    return n_basis
end

function construct_hamiltonian(
        hilbert_space::Array{FullState,1};
        up_coeffs::Union{Nothing,OneBodyCoeffTensor}=nothing, # ↑ single-body coefficients
        down_coeffs::Union{Nothing,OneBodyCoeffTensor}=nothing, # ↓ single-body coefficients
        up_down_coeffs::Union{Nothing,TwoBodyCoeffTensor}=nothing, # ↑↓ interaction
        up_up_coeffs::Union{Nothing,TwoBodyCoeffTensor}=nothing, # ↑↑ interaction
        down_down_coeffs::Union{Nothing,TwoBodyCoeffTensor}=nothing # ↓↓ interaction
    )::Hamiltonian where {T <: Orbital}
    """ Loops through the entire Hilbert space and finds the Hamiltonian Matrix
        elements. Returns a Hamiltonian structure.
    """
    # These values are the entries of the Hamiltonian. Assuming that they wont
    # be larger than 2^32, we can use 32 bit integers here.
    row = Vector{SparseIndexType}()
    col = Vector{SparseIndexType}()
    data = Vector{DType}()

    # Find maximally occupied index.
    n_basis::Integer = find_n_basis(hilbert_space)
    lookup_table, inv_lookup_table = make_lookup_table(hilbert_space)

    # Loop over all states in the Hilbert space.
    for n = 1:length(lookup_table)
        # Split into single spin states.
        s_up, s_down = f_to_s(lookup_table[n])

        up_hops = get_particle_hops(s_up, n_basis)
        down_hops = get_particle_hops(s_down, n_basis)

        # ----------------------
        # One-body terms.

        if !isnothing(up_coeffs)
            for (l,i,new_up,sign) in up_hops
                new_state = get(inv_lookup_table, s_to_f(new_up, s_down), nothing)
                if !isnothing(new_state)
                    push!(col, n)
                    push!(row, new_state)
                    push!(data, -sign * up_coeffs[i,l])
                end
            end
        end # End  of ↑ section.

        if !isnothing(down_coeffs)
            for (l,i,new_down,sign) in down_hops
                new_state = get(inv_lookup_table, s_to_f(s_up, new_down), nothing)
                if !isnothing(new_state)
                    push!(col, n)
                    push!(row, new_state)
                    push!(data, -sign * down_coeffs[i,l])
                end
            end
        end # End  of ↓ section.

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

        # Two-body part.
        # Loop structure:
        #   k -> mu anihilator
        #   i -> mu creator
        #   l -> nu anihilator
        #   j -> nu creator

        # UP/DOWN loop (mu =  ↑, nu = ↓).
        if !isnothing(up_down_coeffs)
            for (k,i,new_up,sign_up) in up_hops
                for (l,j,new_down,sign_down) in down_hops
                    new_state = get(inv_lookup_table, s_to_f(new_up, new_down), nothing)
                    if !isnothing(new_state)
                        push!(col, n)
                        push!(row, new_state)
                        push!(data, sign_up*sign_down * up_down_coeffs[i,j,k,l])
                    end
                end
            end
        end # End  of ↑↓ section.

        # UP/UP loop  (mu =  ↑, nu = ↑).
        if !isnothing(up_up_coeffs)
            for (k,i,temp_up,temp_sign) in up_hops
                temp_hops = get_particle_hops(temp_up, n_basis)
                for (l,j,new_up,new_sign) in temp_hops
                    new_state = get(inv_lookup_table, s_to_f(new_up, s_down), nothing)
                    if !isnothing(new_state)
                        push!(col, n)
                        push!(row, new_state)
                        push!(data, temp_sign*new_sign * up_up_coeffs[i,j,k,l])
                    end
                end
            end
        end # End  of ↑↑ section.

        # DOWN/DOWN loop  (mu =  ↓, nu = ↓).
        if !isnothing(down_down_coeffs)
            for (k,i,temp_down,temp_sign) in down_hops
                temp_hops = get_particle_hops(temp_down, n_basis)
                for (l,j,new_down,new_sign) in temp_hops
                    new_state = get(inv_lookup_table, s_to_f(s_up, new_down), nothing)
                    if !isnothing(new_state)
                        push!(col, n)
                        push!(row, new_state)
                        push!(data, temp_sign*new_sign * down_down_coeffs[i,j,k,l])
                    end
                end
            end
        end # End  of ↓↓ section.

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

    end # End of state loop.

    # Return a Hamiltonian object.
    return Hamiltonian(row, col, data, length(lookup_table))
end


text:
@testset "maha" begin
    rng = MersenneTwister(123456)
    x = 2 * rand(rng)
    D_in = 3
    v1 = rand(rng, D_in)
    v2 = rand(rng, D_in)

    U = UpperTriangular(rand(rng, 3,3))
    P = Matrix(Cholesky(U, 'U', 0))
    @assert isposdef(P)

    k = MahalanobisKernel(P=P)

    @test kappa(k, x) == exp(-x)
    @test k(v1, v2) ≈ exp(-sqmahalanobis(v1, v2, P))
    @test kappa(ExponentialKernel(), x) == kappa(k, x)
    @test repr(k) == "Mahalanobis Kernel (size(P) = $(size(P)))"

    M1, M2 = rand(rng,3,2), rand(rng,3,2)

    function FiniteDifferences.to_vec(dist::SqMahalanobis)
        return vec(dist.qmat), x -> SqMahalanobis(reshape(x, size(dist.qmat)...))
    end
    a = rand()

    function test_mahakernel(U::UpperTriangular, v1::AbstractVector, v2::AbstractVector)
        return MahalanobisKernel(P=Array(U'*U))(v1, v2)
    end

    @test all(FiniteDifferences.j′vp(FDM, test_mahakernel, a, U, v1, v2)[1] .≈
        UpperTriangular(Zygote.pullback(test_mahakernel, U, v1, v2)[2](a)[1]))

    function test_sqmaha(U::UpperTriangular, v1::AbstractVector, v2::AbstractVector)
        return SqMahalanobis(Array(U'*U))(v1, v2)
    end

    @test all(FiniteDifferences.j′vp(FDM, test_sqmaha, a, U, v1, v2)[1] .≈
    UpperTriangular(Zygote.pullback(test_sqmaha, U, v1, v2)[2](a)[1]))

    # test_ADs(U -> MahalanobisKernel(P=Array(U' * U)), U, ADs=[:Zygote])
    @test_broken "Nothing passes (problem with Mahalanobis distance in Distances)"

    # Standardised tests.
    @testset "ColVecs" begin
        x0 = ColVecs(randn(D_in, 3))
        x1 = ColVecs(randn(D_in, 3))
        x2 = ColVecs(randn(D_in, 2))
        TestUtils.test_interface(k, x0, x1, x2)
    end
    @testset "RowVecs" begin
        x0 = RowVecs(randn(3, D_in))
        x1 = RowVecs(randn(3, D_in))
        x2 = RowVecs(randn(2, D_in))
        TestUtils.test_interface(k, x0, x1, x2)
    end
    test_params(k, (P,))
end


text:
##############################################################################
#
# Split string into tokens on whitespace
#
##############################################################################

tokenize{S <: Language, T <: AbstractString}(::Type{S}, s::T) = matchall(r"[^\s]+", s)


text:
import QMTK: forward

include("Forward_CPU.jl")
include("Forward_GPU.jl")


text:
using PolyChaos, GaussQuadrature, Test, LinearAlgebra

mytol = 1e-10
N = 5:10:100

names = ["hermite", "legendre"]
funs = [hermite, legendre, laguerre]

@testset "Compare $(names...) against GaussQuadrature" begin
    for i=1:length(names), n in N
        opq = OrthoPolyQ(names[i],n)
        mynw = nw(opq)
        n, w = funs[i](n)
        @test isapprox(norm(n-mynw[:,1],Inf),0;atol=mytol)
        @test isapprox(norm(w-mynw[:,2],Inf),0;atol=mytol)
    end
end

a = -0.5:0.1:2.0
b = copy(a)

@testset "Compare Jacobi against GaussQuadrature" begin
    for aa in a, bb in b, n in N
        opq = OrthoPolyQ("jacobi",n,Dict(:shape_a=>aa,:shape_b=>bb))
        mynw = nw(opq)
        n, w = jacobi(n,aa,bb)
        @test isapprox(norm(n-mynw[:,1],Inf),0;atol=mytol)
        @test isapprox(norm(w-mynw[:,2],Inf),0;atol=mytol)
    end
end

@testset "Compare Laguerre against GaussQuadrature" begin
    for n in N
        opq = OrthoPolyQ("laguerre",n,Dict(:shape=>0.,:rate=>1.))
        mynw = nw(opq)
        n, w = laguerre(n,0.)
        @test isapprox(norm(n-mynw[:,1],Inf),0;atol=mytol)
        @test isapprox(norm(w-mynw[:,2],Inf),0;atol=mytol)
    end
end


text:
#!/usr/bin/env julia

using DatagenCopulaBased
using NPZ
using Distributions
using RandomNumbers.Xorshifts

## tests on simulated sata with custom rmg

c = Gumbel_cop(2, 5.)
cg = Nested_Gumbel_cop([c], 1, 2.)

c = Frank_cop(2, 5.)
cf = Nested_Frank_cop([c], 1, 2.)

c = Clayton_cop(2, 5.)
cc = Nested_Clayton_cop([c], 1, 2.)

c = AMH_cop(2, .9)
ca = Nested_AMH_cop([c], 1, .3)

cch = Chain_of_Archimedeans([-0.9, 3], "clayton")

cmo = Marshall_Olkin_cop([1., 2., 3., 1., 2., 3., 4.])

cf = Frechet_cop(3, 0.5)

Σ = [1. 0.5 0.5; 0.5 1. 0.5; 0.5 0.5 1.]
ct = Student_cop(Σ, 1)

t = 100_000

r = Xoroshiro128Plus()
x = simulate_copula(t, cc; rng = r)
npzwrite("copula.npy", x)

Σ = [1. 0.8 0.6 0.6; 0.8 1. 0.6 0.6; 0.6 0.6 1. 0.6; 0.6 0.6 0.6 1.]
x = Array(transpose(rand(MvNormal(Σ),t)))
y = gcop2arch(x, ["clayton"=>[1,2,3]])
npzwrite("transform.npy", y)


text:
# Sample network and compute observables for observables stored as
# sparse matrices, therefore where we have a list of non-zero elements in the
# column.
function sample_network!(
  res::MCMCObsEvaluationCache,
  problem::ObservablesProblem,
  net::MatrixNet, σ, wholespace=false)

  σp = parent(res.σ)

  # The denominator of this state
  lnψ = net(σ)
  prob = wholespace ? exp(real(lnψ)) : 1.0
  res.Zave += prob
  # Iterate through all elements in row i_σ of the matrix computing
  # Since Julia has no CSR matrices, I iterate on columns of the
  # transpose matrix
  i_σ = index(σ)

  set_index!(row(σp), i_σ)
  set_index!(col(σp), i_σ)

  # TODO : this works only with NDM. Should generalize it to states...
  for (obs_id, Oᵗ) = enumerate(problem.ObservablesTransposed)
    O_loc = 0.0+0.0im
    for row_id = Oᵗ.colptr[i_σ]:(Oᵗ.colptr[i_σ+1]-1)
      # Find nonzero elements s by doing <i_sp|Oᵗ|i_σ>
      i_σp = Oᵗ.rowval[row_id]
      # BackConvert to int
      set_index!(col(σp), i_σp)
      # Compute the log(ψ(σ)/ψ(σ')), by only computing differences.
      log_ratio = net(σp) - lnψ
      O_loc += (Oᵗ.nzval[row_id]) * exp(log_ratio)
    end
    res.ObsAve[obs_id] += prob * O_loc
    push!(res.ObsVals[obs_id], prob * O_loc)
  end
  return res
end

# Sample network and compute observables for a KLocalOperator representation
# of observables without using lookuptables
#
function sample_network!(
  res::MCMCObsEvaluationCache,
  problem::ObservablesProblem{B,SM},
  net::MatrixNet, σ, wholespace=false) where {B, SM<:AbsLinearOperator}

  σp = parent(res.σ)

  # The denominator of this state
  lnψ = net(σ)
  prob = wholespace ? exp(real(lnψ)) : 1.0
  res.Zave += prob
  # Iterate through all elements in row i_σ of the matrix computing
  # Since Julia has no CSR matrices, I iterate on columns of the
  # transpose matrix
  i_σ = index(σ)

  set_index!(row(σp), i_σ)
  set_index!(col(σp), i_σ)

  # TODO : this works only with NDM. Should generalize it to states...
  for (obs_id, O) = enumerate(problem.ObservablesTransposed)
    O_loc = 0.0+0.0im
    #diffs_O = row_valdiff(O, row(σ.parent))
    for op=operators(O)
      r=local_index(row(σ.parent), sites(op))
      for (mel, changes)=op.op_conns[r]
        set_index!(col(σp), i_σ)
        for (site, val)=changes
          setat!(σp, site, val)
        end
        # Compute the log(ψ(σ)/ψ(σ')), by only computing differences.
        log_ratio = net(σp) - lnψ
        O_loc += conj(mel) * conj(exp(log_ratio)) # maybe a conj
      end
    end

    res.ObsAve[obs_id] += prob * O_loc
    push!(res.ObsVals[obs_id], prob * O_loc)
  end
  return res
end

function sample_network!(
  res::MCMCObsEvaluationCache,
  problem::ObservablesProblem{B,SM},
  net::MatrixNet, σ::LUState, wholespace=false) where {B, SM<:AbsLinearOperator}

  # The denominator of this state
  prob = wholespace ? exp(real(net(σ))) : 1.0
  res.Zave += prob
  #i_σ = index(σ)
  σs = state(σ)
  no_changes = changes(row(𝝝s))

  for (obs_id, O) = enumerate(problem.ObservablesTransposed)
    O_loc = 0.0+0.0im
    diffs_O = row_valdiff(O, row(σs))
    for (mel, changes)=diffs_O
      # Compute the log(ψ(σ)/ψ(σ')), by only computing differences.
      log_ratio = Δ_logψ(net, σ, changes, no_changes)
      O_loc += mel * exp(log_ratio)
    end

    res.ObsAve[obs_id] += prob * O_loc
    push!(res.ObsVals[obs_id], prob * O_loc)
  end
  return res
end


text: @time sumby_dict(id6, v1)


text:
function test_propagate(scenario)
    @unpack tx, rx, bfield, species, ground = scenario

    waveguide = HomogeneousWaveguide(bfield, species, ground)

    E, amp, phase = propagate(waveguide, tx, rx)
    @test eltype(E) == ComplexF64
    @test eltype(amp) == Float64
    @test eltype(phase) == Float64

    mesh = LMP.defaultmesh(tx.frequency)
    E2, amp2, phase2 = propagate(waveguide, tx, rx; mesh=mesh)
    @test E2 ≈ E    rtol=1e-3
    @test amp2 ≈ amp    rtol=1e-3
    @test phase2 ≈ phase    rtol=1e-3

    me = PhysicalModeEquation(tx.frequency, waveguide)
    modes = findmodes(me, mesh)
    E3, amp3, phase3 = propagate(waveguide, tx, rx; modes=modes)
    @test E3 ≈ E2    rtol=1e-3
    @test amp3 ≈ amp2    rtol=1e-3
    @test phase3 ≈ phase2    rtol=1e-3

    # mesh should be ignored
    E4, amp4, phase4 = propagate(waveguide, tx, rx; modes=modes, mesh=[1.0])
    @test E4 ≈ E3    rtol=1e-3
    @test amp4 ≈ amp3    rtol=1e-3
    @test phase4 ≈ phase3    rtol=1e-3

    # Are params being carried through?
    params = LMPParams(earthradius=6300e3)
    E5, amp5, phase5 = propagate(waveguide, tx, rx; params=params)
    @test !isapprox(E5, E, rtol=1e-3)
    @test !isapprox(amp5, amp, rtol=1e-3)
    @test !isapprox(phase5, phase, rtol=1e-3)

    # Mismatch between modes and modeterms
    E6, amp6, phase6 = propagate(waveguide, tx, rx; modes=modes, params=params)
    @test !isapprox(E6, E5, rtol=1e-3)
    @test !isapprox(amp6, amp5, rtol=1e-3)
    @test !isapprox(phase6, phase5, rtol=1e-3)
end

function test_propagate_segmented(scenario)
    @unpack tx, rx, bfield, species, ground, distances = scenario

    waveguide = SegmentedWaveguide([HomogeneousWaveguide(bfield[i], species[i], ground[i],
                                                         distances[i]) for i in 1:2])

    E, amp, phase = propagate(waveguide, tx, rx)
    @test eltype(E) == ComplexF64
    @test eltype(amp) == Float64
    @test eltype(phase) == Float64

    E1, amp1, phase1 = propagate(waveguide, tx, rx)
    @test E1 ≈ E    rtol=1e-3
    @test amp1 ≈ amp    rtol=1e-3
    @test phase1 ≈ phase   rtol=1e-3

    # Are params being carried through?
    params = LMPParams(earthradius=6300e3)
    E3, amp3, phase3 = propagate(waveguide, tx, rx; params=params)
    @test !isapprox(E3, E, rtol=1e-3)
    @test !isapprox(amp3, amp, rtol=1e-3)
    @test !isapprox(phase3, phase, rtol=1e-3)
end

function test_mcranges_segmented(scenario)
    # Test that single receiver location has correct phase without continuous Range
    # This also confirms that results are correct if mode conversion occurs at a
    # segment_range that is not also an output_range

    @unpack tx, rx, bfield, species, ground, distances = scenario

    waveguide = SegmentedWaveguide([HomogeneousWaveguide(bfield[i], species[i], ground[i],
                                                         distances[i]) for i in 1:2])

    _, ampref, phaseref = propagate(waveguide, tx, rx; unwrap=false)

    _, a1, p1 = propagate(waveguide, tx, GroundSampler(600e3, Fields.Ez))
    _, a2, p2 = propagate(waveguide, tx, GroundSampler(1400e3, Fields.Ez))
    _, a3, p3 = propagate(waveguide, tx, GroundSampler(1800e3, Fields.Ez))

    m1 = findfirst(isequal(600e3), rx.distance)
    m2 = findfirst(isequal(1400e3), rx.distance)
    m3 = findfirst(isequal(1800e3), rx.distance)

    @test a1 ≈ ampref[m1] atol=0.01
    @test a2 ≈ ampref[m2] atol=0.01
    @test a3 ≈ ampref[m3] atol=0.01
    @test p1 ≈ phaseref[m1] atol=1e-3
    @test p2 ≈ phaseref[m2] atol=1e-3
    @test p3 ≈ phaseref[m3] atol=1e-3
end

@testset "LongwaveModePropagator.jl" begin
    @info "Testing LongwaveModePropagator"

    @info "  Running:"

    # Just to save time, running with only one scenario
    @info "    Homogeneous ionospheres..."
    for scn in (resonant_scenario, interp_scenario)
        test_propagate(scn)
    end

    @info "    Segmented ionospheres..."
    for scn in (segmented_scenario,)
        test_propagate_segmented(scn)
    end

    # propagate(file,...) is tested in IO.jl
end


text:
using Core.Compiler: AbstractInterpreter, NativeInterpreter, InferenceState,
    OptimizationState, CodeInfo, CodeInstance, InferenceResult, WorldRange,
    IRCode, SSAValue, inlining_policy

struct InferredSource
    src::CodeInfo
    stmt_info::Vector{Any}
    effects::Effects
    rt::Any
    InferredSource(src::CodeInfo, stmt_info::Vector{Any}, effects, @nospecialize(rt)) =
        new(src, stmt_info, effects, rt)
end

struct OptimizedSource
    ir::IRCode
    src::CodeInfo
    isinlineable::Bool
    effects::Effects
end

maybe_create_optsource(@nospecialize(a), ::Effects) = a

function maybe_create_optsource(st::OptimizationState, effects::Effects)
    ir = st.ir
    ir === nothing && return st
    return OptimizedSource(ir, st.src, st.src.inlineable, effects)
end

const Remarks = Vector{Pair{Int, String}}

struct CthulhuInterpreter <: AbstractInterpreter
    native::AbstractInterpreter

    unopt::Dict{Union{MethodInstance, InferenceResult}, InferredSource}
    opt::Dict{MethodInstance, CodeInstance}

    remarks::Dict{MethodInstance, Remarks}
end

CthulhuInterpreter(interp::AbstractInterpreter=NativeInterpreter()) =
    CthulhuInterpreter(
        interp,
        Dict{MethodInstance, InferredSource}(),
        Dict{MethodInstance, CodeInstance}(),
        Dict{MethodInstance, Remarks}()
    )

import Core.Compiler: InferenceParams, OptimizationParams, get_world_counter,
    get_inference_cache, code_cache, lock_mi_inference, unlock_mi_inference, method_table
import Core.Compiler: InferenceState, WorldView
using Base: @invoke

Compiler.InferenceParams(interp::CthulhuInterpreter) = InferenceParams(interp.native)
Compiler.OptimizationParams(interp::CthulhuInterpreter) = OptimizationParams(interp.native)
Compiler.get_world_counter(interp::CthulhuInterpreter) = get_world_counter(interp.native)
Compiler.get_inference_cache(interp::CthulhuInterpreter) = get_inference_cache(interp.native)

# No need to do any locking since we're not putting our results into the runtime cache
Compiler.lock_mi_inference(interp::CthulhuInterpreter, mi::MethodInstance) = nothing
Compiler.unlock_mi_inference(interp::CthulhuInterpreter, mi::MethodInstance) = nothing
@static if v"1.8-beta2" <= VERSION < v"1.9-" || VERSION >= v"1.9.0-DEV.120"
Compiler.method_table(interp::CthulhuInterpreter) = method_table(interp.native)
else # if v"1.8-beta2" <= VERSION < v"1.9-" || VERSION >= v"1.9.0-DEV.120"
Compiler.method_table(interp::CthulhuInterpreter, sv::InferenceState) = method_table(interp.native, sv)
end # if v"1.8-beta2" <= VERSION < v"1.9-" || VERSION >= v"1.9.0-DEV.120"
struct CthulhuCache
    cache::Dict{MethodInstance, CodeInstance}
end

Compiler.code_cache(interp::CthulhuInterpreter) = WorldView(CthulhuCache(interp.opt), WorldRange(get_world_counter(interp)))
Compiler.get(wvc::WorldView{CthulhuCache}, mi::MethodInstance, default) = get(wvc.cache.cache, mi, default)
Compiler.haskey(wvc::WorldView{CthulhuCache}, mi::MethodInstance) = haskey(wvc.cache.cache, mi)
Compiler.setindex!(wvc::WorldView{CthulhuCache}, ci::CodeInstance, mi::MethodInstance) = setindex!(wvc.cache.cache, ci, mi)

Compiler.may_optimize(interp::CthulhuInterpreter) = true
Compiler.may_compress(interp::CthulhuInterpreter) = false
Compiler.may_discard_trees(interp::CthulhuInterpreter) = false
Compiler.verbose_stmt_info(interp::CthulhuInterpreter) = true

function Compiler.add_remark!(interp::CthulhuInterpreter, sv::InferenceState, msg)
    if !haskey(interp.remarks, sv.linfo)
        interp.remarks[sv.linfo] = Remarks()
    end
    push!(interp.remarks[sv.linfo], sv.currpc => msg)
end

function Compiler.finish(state::InferenceState, interp::CthulhuInterpreter)
    res = @invoke Compiler.finish(state::InferenceState, interp::AbstractInterpreter)
    key = Core.Compiler.any(state.result.overridden_by_const) ? state.result : state.linfo
    interp.unopt[key] = InferredSource(
        copy(state.src),
        copy(state.stmt_info),
        isdefined(Core.Compiler, :Effects) ? state.ipo_effects : nothing,
        state.result.result)
    return res
end

function Compiler.transform_result_for_cache(interp::CthulhuInterpreter, linfo::MethodInstance,
        valid_worlds::WorldRange, @nospecialize(inferred_result), ipo_effects::Effects = Effects())
    return maybe_create_optsource(inferred_result, ipo_effects)
end

# branch on https://github.com/JuliaLang/julia/pull/41328
@static if isdefined(Compiler, :is_stmt_inline)
function Compiler.inlining_policy(
    interp::CthulhuInterpreter, @nospecialize(src), stmt_flag::UInt8,
    mi::MethodInstance, argtypes::Vector{Any})
    @assert isa(src, OptimizedSource) || isnothing(src)
    if isa(src, OptimizedSource)
        if Compiler.is_stmt_inline(stmt_flag) || src.isinlineable
            return src.ir
        end
    else
        # the default inlining policy may try additional effor to find the source in a local cache
        return @invoke Compiler.inlining_policy(
            interp::AbstractInterpreter, nothing, stmt_flag::UInt8,
            mi::MethodInstance, argtypes::Vector{Any})
    end
    return nothing
end
else # @static if isdefined(Compiler, :is_stmt_inline)
function Compiler.inlining_policy(interp::CthulhuInterpreter)
    function (src)
        @assert isa(src, OptimizedSource)
        return src.isinlineable ? src.ir : nothing
    end
end
end # @static if isdefined(Compiler, :is_stmt_inline)

function Compiler.finish!(interp::CthulhuInterpreter, caller::InferenceResult)
    effects = EFFECTS_ENABLED ? caller.ipo_effects : nothing
    caller.src = maybe_create_optsource(caller.src, effects)
end


text:
using IDXParser
using Base.Test

# write your own tests here
@test 1 == 1

cwd = Base.dirname(Base.source_path(nothing))

labels = IDXParser.parseIDX("$(cwd)/t10k-labels-idx1-ubyte")

@test typeof(labels) == Array{UInt8, 1}
@test size(labels)   == (10000,)

images = IDXParser.parseIDX("$(cwd)/t10k-images-idx3-ubyte")

@test typeof(images) == Array{UInt8, 3}
@test size(images)   == (10000,28,28)


text:
function Acrobot()
    opts = SolverOptions(
        cost_tolerance_intermediate=1e-2,
        penalty_scaling = 1000.,
        penalty_initial = 0.001,
    )
    # model
    model = RobotZoo.Acrobot()
    n,m = RD.dims(model)

    # discretization
    tf = 5.0
    N = 101
    dt = tf/(N-1)

    # initial and final conditions
    x0 = @SVector [-pi/2, 0, 0, 0]
    xf = @SVector [+pi/2, 0, 0, 0]

    # objective
    Q = Diagonal(@SVector fill(1.0, n))
    R = Diagonal(@SVector fill(0.01, m))
    Qf = 100*Q
    obj = LQRObjective(Q*dt,R*dt,Qf,xf,N)

    # constraints
    conSet = ConstraintList(n,m,N)
    goal = GoalConstraint(xf)
    bnd  = BoundConstraint(n,m, u_min=-15, u_max=15)
    add_constraint!(conSet, goal, N:N)
    add_constraint!(conSet, bnd, 1:N-1)

    # initialization
    u0 = @SVector fill(0.0,m)

    # set up problem
    prob = Problem(model, obj, x0, tf, xf=xf, constraints=conSet)
    rollout!(prob)

    return prob, opts
end


text:
export node2tree, tree_findlabel, node_findlabel, node_findkey, node_find_leafkey, node_findkey_safe, share_labels
export node_clade, node_leavesclade, node_leavesclade_labels, tree_clade, tree_leavesclade, isclade, node_findroot
export lca, node_depth, node_divtime, node_ancestor_list, isancestor

###############################################################################################################
################################# Trees from nodes, finding labels, ... #######################################
###############################################################################################################

"""
	node2tree(root::TreeNode)

Create a `Tree` object from `root`. Keys are integers.
"""
function node2tree(root::TreeNode{T} where T)
	tree = Tree(root)
	key = 1
	leafkey = 1
	node2tree_addnode!(tree, root, key, leafkey)
	return tree
end

"""
	node2tree_addnode!(tree::Tree, node::TreeNode, key::Int64 ; addchildren = true)

Add existing `node::TreeNode` to `tree::Tree` using `key`. Recursively add children of `node` if `addchildren` is true.
Return new value of `key` to add further children to `tree`.
If `node` is a leaf node, also add it to `tree.leaves` with key `leafkey`. Return new value of `leafkey`.
"""
function node2tree_addnode!(tree::Tree, node::TreeNode, key::Int64, leafkey::Int64; addchildren = true)
	if in(key, keys(tree.nodes)) || in(leafkey, keys(tree.leaves))
		error("Trying to add node to an already existing key.")
	else
		tree.nodes[key] = node
		tree.lnodes[node.label] = node
		if node.isleaf
			tree.leaves[leafkey] = node
			tree.lleaves[node.label] = node
		end
		if addchildren
			ckey = key + 1
			node.isleaf ? cleafkey = leafkey + 1 : cleafkey = leafkey
			for c in node.child
				ckey, cleafkey = node2tree_addnode!(tree, c, ckey, cleafkey)
			end
		end
	end
	return ckey, cleafkey
end

"""
	name_nodes!(t::Tree)

Give a label to label-less nodes in `t`.
"""
function name_nodes!(t::Tree)
	name_nodes!(t.root, collect(keys(t.lnodes)))
end
function name_nodes!(r::TreeNode, labels ; i = 0)
	ii = i
	if r.isleaf && isempty(r.label)
		@warn "Label-less leaf node!"
	end
	if isempty(r.label)
		while in("NODE_$ii", labels)
			ii += 1
		end
		r.label = "NODE_$ii"
		ii += 1
		for c in r.child
			ii = name_nodes!(c, labels, i = ii)
		end
	end
	return ii
end

"""
	tree_findlabel(label::String, tree::Tree)

Find a sequence with label `label` in `tree`, and return a `(key, flag)` tuple. `key` is the key of the found sequence in `tree.nodes`. Defaults to `0` if no match is found. `flag` is `true` if a match is found, `false` otherwise.
"""
function tree_findlabel(label::String, tree::Tree)
	for k in keys(tree.nodes)
		if tree.nodes[k].label == label
			return k, true
		end
	end
	return 0, false
end

"""
	node_findlabel(label::String, root::TreeNode ; subtree = true)

Find label in tree defined by `root`. If `subtree`, only the children of `root` are searched. Otherwise, the whole tree is searched.

# Note
`subtree = false` is not yet implemented.
"""
function node_findlabel(label::String, root::TreeNode ; subtree = true)
	found, node = _node_findlabel(label, root)
	if !found
		@warn "Label $(label) was not found."
	end
	return node
end

"""
"""
function _node_findlabel(label::String, root::TreeNode)
	if root.label == label
		return true, root
	end
	for c in root.child
		found, out = _node_findlabel(label, c)
		if found
			return found, out
		end
	end
	return false, nothing
end

"""
	node_findkey(node, tree)

Find key corresponding to `node` in `tree`.
Return value is `nothing` if `node` is not found. --> not type stable.
"""
function node_findkey(node, tree)
	for i in keys(tree.nodes)
		if node == tree.nodes[i]
			return i
		end
	end
	return nothing
end

"""
	node_find_leafkey(node, tree)

Find leafkey corresponding to `node` in `tree`.
Return value is `nothing` if `node` is not found. --> not type stable.
"""
function node_find_leafkey(node, tree)
	if !node.isleaf
		return nothing
	end
	for i in keys(tree.leaves)
		if node == tree.leaves[i]
			return i
		end
	end
	return nothing
end

"""
	node_findkey_safe(node,  tree)

Type safe implementation of `node_findkey`. If nothing is found, an error is raised.
"""
function node_findkey_safe(node,tree)
	a = node_findkey(node,tree)
	if a == nothing
		error("Node $(node.label) was not found in tree.")
	end
	return a
end

"""
	share_labels(tree1, tree2)

Check if `tree1` and `tree2` share the same labels for leaf nodes.
"""
function share_labels(tree1, tree2)
	l1 = Set(l for l in keys(tree1.lleaves))
	l2 = Set(l for l in keys(tree2.lleaves))
	return l1 ==  l2
end

###############################################################################################################
################################################### Clades ####################################################
###############################################################################################################

"""
	node_clade(root::TreeNode)

Find and return clade corresponding to all descendants of `root`.
"""
function node_clade(root::TreeNode)
	if root.isleaf
		return [root]
	end
	clade = [root]
	for c in root.child
		append!(clade, node_clade(c))
	end
	return clade
end

"""
	node_clade_labels(root::TreeNode)

Find and return labels of nodes in clade corresponding to all descendants of `root`.
"""
function node_clade_labels(root::TreeNode)
	if root.isleaf
		return [root.label]
	end
	clade = [root.label]
	for c in root.child
		append!(clade, node_clade(c))
	end
	return clade
end

"""
	node_leavesclade(root::TreeNode)

Find and return clade corresponding to all descendants of `root` that are leaves.
"""
function node_leavesclade(root::TreeNode)
	if root.isleaf
		return [root]
	end
	clade = Array{TreeNode,1}(undef, 0)
	for c in root.child
		append!(clade, node_leavesclade(c))
	end
	return clade
end

"""
	node_leavesclade_labels(root::TreeNode)

Find and return labels of nodes in clade corresponding to all descendants of `root` that are leaves.
"""
function node_leavesclade_labels(root::TreeNode)
	if root.isleaf
		return [root.label]
	end
	clade = Array{String}(undef,0)
	for c in root.child
		append!(clade, node_leavesclade_labels(c))
	end
	return clade
end

"""
	tree_clade(tree::Tree, key)

Find and return keys of clade corresponding to all descendants of `tree.nodes[key]`.
"""
function tree_clade(tree::Tree, key)
	cl = node_clade(tree.nodes[key])
	out = map(x->node_findkey(x, tree), cl)
	return out
end

"""
	tree_leavesclade(tree::Tree, key)

Find and return leaves keys of clade corresponding to all leaves descendants of `tree.nodes[key]`.
"""
function tree_leavesclade(tree::Tree, key)
	cl = node_leavesclade(tree.nodes[key])
	return map(x->node_find_leafkey(x, tree), cl)
end

"""
	node_findroot(node::TreeNode ; maxdepth=1000)

Return root of the tree to which `node` belongs.
"""
function node_findroot(node::TreeNode ; maxdepth=1000)
	temp = node
	it = 0
	while !temp.isroot || it>maxdepth
		temp = temp.anc
		it += 1
	end
	if it>maxdepth
		@error "Could not find root after $maxdepth iterations."
	end
	return temp
end

"""
"""
function node_ancestor_list(node::TreeNode)
	list = [node.label]
	a = node
	while !a.isroot
		push!(list, a.anc.label)
		a = a.anc
	end
	return list
end

"""
	isclade(nodelist)

Check if `nodelist` is a clade. All nodes in `nodelist` should be leaves.
"""
function isclade(nodelist)
	out = true
	if !mapreduce(x->x.isleaf, *, nodelist, init=true)
		# verbose && println("F")
		out = false
	else
		claderoot = lca(nodelist)
		clade = node_leavesclade_labels(claderoot)
		# Now, checking if `clade` is the same as `nodelist`
		for c in clade
			flag = false
			for n in nodelist
				if n.label==c
					flag = true
					break
				end
			end
			if !flag
				out = false
				break
			end
		end
	end
	return out
end

"""
	find_clades(tree, label_list)

Find the set of clades for which leaves are exactly `label_list`.
"""
function find_clades(tree, label_list)
	r = lca([tree.lleaves[l] for l in label_list])
	cr = node_leavesclade_labels(r)
	if issetequal(cr, label_list) # `label_list` is the clade of `r`
		return [label_list]
	else # Go up from ll[1] and get the biggest possible clade
		ll = copy(label_list)
		return find_clades_!(ll, tree)
	end
end
function find_clades_!(list, tree)
	cl = []
	while !isempty(list)
		tmp = []
		r = tree.lleaves[list[1]]
		rflag = true
		idx = [1]
		while rflag # Go up from r as long as we can
			r = r.anc
			cr = node_leavesclade_labels(r)
			if !issubset(cr, list)
				rflag = false
			else
				idx = findall(x->in(x,cr),list) # Which elements of list are in cr
			end
		end
		push!(cl, list[idx])
		deleteat!(list, idx)
	end
	return cl
end

function find_clades!(r, nlist)

end

###############################################################################################################
######################################## LCA, divtime, ... ####################################################
###############################################################################################################
"""
	lca(i_node::TreeNode, j_node::TreeNode)

Find and return lowest common ancestor of `i_node` and `j_node`.
Idea is to go up in the tree in an asymmetric way on the side of the deeper node, until both are at equal distance from root. Then, problem is solved by going up in a symmetric way. (https://stackoverflow.com/questions/1484473/how-to-find-the-lowest-common-ancestor-of-two-nodes-in-any-binary-tree/6183069#6183069)
"""
function lca(i_node::TreeNode, j_node::TreeNode)

	if i_node.isroot
		return i_node
	elseif j_node.isroot
		return j_node
	end

	ii_node = i_node
	jj_node = j_node

	di = node_depth(ii_node)
	dj = node_depth(jj_node)
	while di != dj
		if di > dj
			ii_node = ii_node.anc
			di -=1
		else
			jj_node = jj_node.anc
			dj -=1
		end
	end
	while ii_node != jj_node
		ii_node = ii_node.anc
		jj_node = jj_node.anc
	end
	return ii_node
end

"""
	lca(nodelist)

Find the common ancestor of all nodes in `nodelist`. `nodelist` is an iterable collection of `TreeNode` objects.
"""
function lca(nodelist)
	# ca = nodelist[1]
	# Getting any element to start with
	ca = first(nodelist)
	for node in nodelist
		if !isancestor(ca, node)
			ca = lca(ca, node)
		end
	end
	return ca
end

"""
	isancestor(a:::TreeNode, node::TreeNode)

Check if `a` is an ancestor of `node`.
"""
function isancestor(a::TreeNode, node::TreeNode)
	if a==node
		return true
	else
		if node.isroot
			return false
		else
			return isancestor(a, node.anc)
		end
	end
end

"""
	node_depth(node::TreeNode)

Distance from `node` to root.
"""
function node_depth(node::TreeNode)
	d = 0
	_node = node
	while !_node.isroot
		_node = _node.anc
		d += 1
	end
	return d
end

"""
	node_divtime(i_node::TreeNode, j_node::TreeNode)

Compute divergence time between `i_node` and `j_node` by summing the `TreeNode.data.tau` values.
"""
function node_divtime(i_node::TreeNode, j_node::TreeNode)
	a_node = lca(i_node, j_node)
	tau = 0
	ii_node = i_node
	jj_node = j_node
	while ii_node != a_node
		tau += ii_node.data.tau
		ii_node = ii_node.anc
	end
	while jj_node != a_node
		tau += jj_node.data.tau
		jj_node = jj_node.anc
	end
	return tau
end

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility


text:
#=
@inline foldlargs(op, x) = x
@inline foldlargs(op, x1, x2, xs...) = foldlargs(op, op(x1, x2), xs...)

@inline tuplediff(xs, ys) = foldlargs((), xs...) do zs, x
    Base.@_inline_meta
    foldlargs(false, ys...) do found, y
        Base.@_inline_meta
        found ? true : y == x
    end ? zs : (zs..., x)
end
=#

@inline tuplediff(xs, ys) = foldlargs((), xs...) do zs, x
    foldlargs(false, ys...) do _, y
        y == x && reduced(zs)
    end |> ifunreduced() do _
        (zs..., x)
    end
end

@inline tupleunion(xs, ys) = (xs..., tuplediff(ys, xs)...)


text:
const excluded_constants = [
    "VK_TRUE",
    "VK_FALSE",
]

function ConstantDefinition(spec::SpecConstant)
    p = Dict(
        :category => :const,
        :name => remove_vk_prefix(spec.name),
        :value => spec.name,
    )
    ConstantDefinition(spec, p)
end

function include_constant(spec::SpecConstant)
    name = string(spec.name)
    startswith(name, "VK_") && !contains(name, r"(SPEC_VERSION|EXTENSION_NAME)") && name ∉ excluded_constants
end


text:
# Error weighted PPM

"""
    ewPPM{Char}( c::Int [,ϵ::Float64 = 1.0] )

Creates an Error Weighted PPM predictor with context length `c` and weighing paramter `ϵ`.

## Examples
```julia-repl
julia> p = ewPPM{Int}(3)
DiscretePredictors.ewPPM{Int}([*] (0)
, Int[], 3, 1.0, [1.0, 1.0, 1.0, 1.0], Dict{Int,Float64}[Dict{Int,Float64}(), Dict{Int,Float64}(), Dict{Int,Float64}(), Dict{Int,Float64}()])

julia> p = ewPPM{Int}(3, 0.95)
DiscretePredictors.ewPPM{Int}([*] (0)
, Int[], 3, 0.95, [1.0, 1.0, 1.0, 1.0], Dict{Int,Float64}[Dict{Int,Float64}(), Dict{Int,Float64}(), Dict{Int,Float64}(), Dict{Int,Float64}()])
```

Reference:
Pulliyakode, Saishankar Katri, and Sheetal Kalyani. "A modified ppm algorithm for online sequence prediction using short data records." IEEE Communications Letters 19.3 (2015): 423-426.
"""
mutable struct ewPPM{T} <: BasePredictor{T}
    model::Trie{T,Int}
    context::Vector{T}
    cxt_length::Int
    ϵ::Float64
    weights::Vector{Float64}
    last_prediction::Vector{Dict{T,Float64}}

    ewPPM{T}( _c::Int, _ϵ::Float64 = 1.0 ) where {T} = new( Trie{T,Int}(), Vector{T}(),
                                                _c, _ϵ,
                                                ones(Float64,_c+1), fill(Dict{T,Float64}(),_c+1) )
end

function add!( p::ewPPM{T}, sym::T ) where {T}
    p.model.value += 1;

    # Update weights
    error_vec   = Vector{Float64}(undef,p.cxt_length+1)
    for expert_idx = 1:p.cxt_length+1
        error_vec[expert_idx] = (get_best_symbol(p.last_prediction[expert_idx])==sym) ? 0.0 : 1.0
    end
    correction  = error_vec .- (1/(p.cxt_length+1))*sum(error_vec)
    for i = 1:p.cxt_length+1
        p.weights[i] = min(1,max(p.weights[i]-p.ϵ*correction[i],0))
    end

    # Create node path
    push!( p.context, sym );
    buffer  = p.context[1:end];
    while !isempty(buffer)
        p.model[buffer] = (haskey(p.model,buffer)) ? p.model[buffer]+1 : 1
        popfirst!( buffer )
    end

    if length(p.context) > p.cxt_length
        popfirst!( p.context )
    end

    nothing
end

function predict( p::ewPPM{T} ) where {T}
    # Clear away last predictions
    p.last_prediction   = fill(Dict{T,Float64}(),p.cxt_length+1)

    # Get prediction from 0-order expert
    symbols = Dict( k => (p.model[[k]]/p.model.value) for k ∈ keys(children(p.model,Vector{T}())) )

    # Save this for later
    p.last_prediction[1]    = deepcopy( symbols )

    # Apply 0-order expert weight
    for sym ∈ keys(symbols)
        symbols[sym] *= p.weights[1];
    end
    # println( "After Order-0: ", symbols );

    # Do prediction from all other higher orders
    no_of_active_experts    = min( p.model.value, p.cxt_length )
    for expert_idx  = 1:no_of_active_experts
        # Predict
        this_expert_prediction = predict_from_subcontext( p, p.context[1:expert_idx] )
        # Normalize probability out from experts - jugaad
        denom   = 0;
        for (sym,prob) ∈ this_expert_prediction denom += prob; end
        for sym ∈ keys(this_expert_prediction) this_expert_prediction[sym] /= denom;    end

        # Save
        p.last_prediction[expert_idx+1] = deepcopy( this_expert_prediction )
        # Apply weight to expert and add to final prediction
        for sym ∈ keys( symbols )
            symbols[sym] += (this_expert_prediction[sym]*p.weights[expert_idx+1])
        end
        # println( "After order-$expert_idx: ", symbols );
    end

    return symbols;
end

function info_string( p::ewPPM{T} ) where {T}
    return @sprintf( "ewPPM(%d,\$\\epsilon\$ = %3.2f)", p.cxt_length, p.ϵ )
end

function unique_string( p::ewPPM{T} ) where {T}
    return @sprintf( "ewPPM_%02d_%03d", p.cxt_length, trunc(Int,100*p.ϵ)  )
end

function predict_from_subcontext( p::ewPPM{T}, sub_cxt::Vector{T} ) where {T}
    # Create a dictionary with symbols
    symbols = Dict( k => (p.model[[k]]/p.model.value) for k ∈ keys(children(p.model,Vector{T}())) )

    buffer  = Vector{T}();
    for i = length(sub_cxt):-1:1
        pushfirst!( buffer, sub_cxt[i] )
        # println( "    Buffer: ", buffer, " -> ", p.model[buffer] );
        # Apply escape probability
        esc_prob = 1/p.model[buffer]
        for symbol in keys(symbols)
            symbols[symbol] *= esc_prob
        end
        # Get each child probability
        list_of_children    = children( p.model, buffer );
        for k ∈ keys(list_of_children)
            # println( "        ", k, " -> ", list_of_children[k].value )
            symbols[k] += list_of_children[k].value/p.model[buffer]
        end
    end

    return symbols
end


text:
module SoilPlantAirContinuum

using ..CanopyLayers: Canopy4RT, CanopyOpticals, CanopyRads, IncomingRadiation,
      LeafBios, RTCache, RTDimensions, SIF_fluxes!, SoilOpticals, SolarAngles,
      WaveLengths, big_leaf_partition, canopy_fluxes!, canopy_geometry!,
      canopy_matrices!, create_canopy_opticals, create_canopy_rads,
      create_canopy_rt, create_incoming_radiation, create_leaf_bios,
      create_rt_cache, create_rt_dims, create_wave_length, fluspect!,
      short_wave!, thermal_fluxes!
using ConstrainedRootSolvers: ReduceStepMethodND, SolutionToleranceND,
      find_peak
using DataFrames: DataFrame
using DocStringExtensions: TYPEDFIELDS
using ..Photosynthesis: AbstractPhotoModelParaSet, AirLayer, C3CLM, GCO₂Mode,
      Leaf, leaf_photosynthesis!, leaf_rd!, leaf_temperature_dependence!
using PkgUtility: GAS_R, GRAVITY, K_STEFAN, K_VON_KARMAN, M_H₂O, P_ATM, RT_25,
      T_0, T_25, YEAR_D, tinfo, ρ_H₂O
using ..PlantHydraulics: AbstractPlantOrganism, GrassLikeOrganism,
      PalmLikeOrganism, SteadyStateMode, TreeLikeOrganism, TreeSimple,
      create_grass, critical_flow, end_pressure, flow_profile!,
      pressure_profile!, roots_flow!, soil_p_25_swc, soil_swc,
      temperature_effects!
using ..StomataModels: AbstractStomatalModel, CanopyLayer, ESMBallBerry,
      EmpiricalStomatalModel, GswDrive, gas_exchange!, gsw_control!,
      prognostic_gsw!, stomatal_conductance, update_leaf_TP!
using UnPack: @unpack
using WaterPhysics: latent_heat_vapor, relative_diffusive_coefficient,
      saturation_vapor_pressure

# define constants
# TODO add unit conversion in PkgUtility
KG_2_MOL(FT)     = 1 / M_H₂O(FT);
KG_H_2_MOL_S(FT) = KG_2_MOL(FT) / 3600;

# export public types
export SPACMono, SPACSimple

# export public functions
export annual_profit, annual_simulation!, atmospheric_pressure,
       atmospheric_pressure_ratio, big_leaf_partition!, create_dataframe,
       gain_risk_map, initialize_spac_canopy!, layer_fluxes!, leaf_allocation!,
       leaf_gas_exchange!, leaf_gas_exchange_nonopt!, leaf_temperature,
       leaf_temperature_shaded, leaf_temperature_sunlit, optimize_flows!,
       optimize_hs!, optimize_leaf!, ppm_to_Pa, test_soil_from_psoil!,
       test_soil_from_swc!, update_Cab!, update_Kmax!, update_LAI!,
       update_VJR!, update_VJRWW!, update_Weibull!, vary_spac!, zenith_angle

include("types/container.jl" )
include("types/spacmono.jl"  )
include("types/spacsimple.jl")

include("bigleaf/gainriskmap.jl"   )
include("bigleaf/gasexchange.jl"   )
include("bigleaf/leafallocation.jl")
include("bigleaf/optimizeflow.jl"  )
include("bigleaf/optimizehs.jl"    )
include("bigleaf/optimizeleaf.jl"  )
include("bigleaf/partition.jl"     )
include("bigleaf/temperature.jl"   )
include("bigleaf/varytrait.jl"     )

include("layers/layer_fluxes.jl")
include("layers/initializert.jl")
include("layers/test_diurnal.jl")
include("layers/test_soil.jl"   )
include("layers/windspeed.jl"   )

include("planet/atmpressure.jl")
include("planet/solarangle.jl" )

include("simulation/annualprofit.jl"    )
include("simulation/annualsimulation.jl")
include("simulation/createdataframe.jl" )

include("soil/moisture.jl")

end # module

max_stars_repo_name: climate-machine/Land

max_issues_repo_name: climate-machine/Land

max_forks_repo_name: climate-machine/Land


text:
using PackVec
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

v = [1.0,2.0,3.0]

n = @unpack!(v, a::Real, b::Real, c::Real)

@test n == 3
@test a == v[1] == 1.0
@test b == v[2] == 2.0
@test c == v[3] == 3.0

n =  @pack!(v, a+=1, b+=2, c+=3)

@test n == 3
@test a == v[1] == 2.0
@test b == v[2] == 4.0
@test c == v[3] == 6.0

v = randn(1+4*4)
X = zeros(4,4)
n = @unpack!(v, a::Real, X)

@test n == 1+4*4
@test a == v[1]
@test X[1:4*4] == v[2:1+4*4]

n = @pack!(v, X, a)
@test n == 1+4*4
@test X[1:4*4] == v[1:4*4]
@test a == v[4*4+1]

v = randn(100)
n = @unpack!(v,
             D = Diagonal(zeros(3)),
             U = UpperTriangular(zeros(3,3)),
             L = LowerTriangular(zeros(3,3)),
             S = Symmetric(zeros(3,3)),
             C = LinAlg.Cholesky(zeros(3,3),'U'))

@test n == 3 + 6 + 6 + 6 + 6
@test D[1,1] == v[1]
@test D[3,3] == v[3]
@test U[1,1] == v[3+1]
@test U[3,3] == v[3+6]
@test L[1,1] == v[3+6+1]
@test L[3,3] == v[3+6+6]
@test S[1,1] == v[3+6+6+1]
@test S[3,3] == v[3+6+6+6]
@test C[:U][1,1] == v[3+6+6+6+1]
@test C[:U][3,3] == v[3+6+6+6+6]

vv = zeros(n)
nn = @pack!(vv, D, U, L, S, C)
@test n == nn
@test vv == v[1:n]

using StaticArrays
n = @unpack!(v, X::SVector{3}, Y::SMatrix{2,4})

@test n == 3 + 2*4
@test X[1] == v[1]
@test X[3] == v[3]
@test Y[1,1] == v[3+1]
@test Y[2,4] == v[3+2*4]

vv = zeros(n)
nn = @pack!(vv, X, Y)
@test n == nn
@test vv == v[1:n]

max_stars_repo_name: simonbyrne/PackVec