text:
export IntrinsicSubspaceHomotopy, set_subspaces!

struct LinearSubspaceGeodesicInfo
    Q::Matrix{ComplexF64}
    Q_cos::Matrix{ComplexF64}
    Θ::Vector{Float64}
    U::Matrix{ComplexF64}
    γ1::Matrix{ComplexF64}
end

function LinearSubspaceGeodesicInfo(start::LinearSubspace, target::LinearSubspace)
    Q, Θ, U = geodesic_svd(target, start)
    Q_cos = target.intrinsic.Y * U
    γ1 = similar(Q_cos)
    n, k = size(γ1)
    for j = 1:k
        s, c = sincos(Θ[j])
        for i = 1:n
            γ1[i, j] = Q_cos[i, j] * c + Q[i, j] * s
        end
    end
    LinearSubspaceGeodesicInfo(Q, Q_cos, Θ, U, γ1)
end

const GEODESIC_LRU = LRU{
    Tuple{LinearSubspace{ComplexF64},LinearSubspace{ComplexF64}},
    LinearSubspaceGeodesicInfo,
}(
    maxsize = 128,
)

"""
    IntrinsicSubspaceHomotopy(F::System, V::LinearSubspace, W::LinearSubspace)
    IntrinsicSubspaceHomotopy(F::AbstractSystem, V::LinearSubspace, W::LinearSubspace)

Creates a homotopy ``H(x,t) = F(γ(t)x)`` where ``γ(t)`` is a family of affine subspaces
such that ``γ(1) = V`` and ``γ(0) = W``.
Here ``γ(t)`` is the geodesic between `V` and `W` in the affine Grassmanian, i.e.,
it is the curve of minimal length connecting `V` and `W`.
See also [`LinearSubspace`](@ref) and [`geodesic`](@ref) and the references therein.
"""
Base.@kwdef mutable struct IntrinsicSubspaceHomotopy{S<:AbstractSystem} <: AbstractHomotopy
    system::S

    start::LinearSubspace{ComplexF64}
    target::LinearSubspace{ComplexF64}
    geodesic::LinearSubspaceGeodesicInfo

    J::Matrix{ComplexF64}
    x::Vector{ComplexF64}
    ẋ::Vector{ComplexF64}
    x_high::Vector{ComplexDF64}
    t_cache::Base.RefValue{ComplexF64}
    # For AD
    taylor_t_cache::Base.RefValue{ComplexF64}
    taylor_γ::NTuple{4,Matrix{ComplexF64}}
    v::Vector{ComplexF64}
    tx³::TaylorVector{4,ComplexF64}
    tx²::TaylorVector{3,ComplexF64}
    tx¹::TaylorVector{2,ComplexF64}
end

## Implementation details

# The computation is performed using (implicit) Stiefel coordinates. If `size(F) = m, n`
# and `V` and `W` are affine subspaces of dimension `k, then the computation is performed in
# the affine Grassmanian Graff(k,n) embedded in the Grassmanian(k+1,n+1) using the system
# ``[F(γ(t)v); (γ(t)v)[end] - 1]``. Here the `(γ(t)v)[k+1] - 1` condition ensures that
# we are in the affine Grassmanian and that `(γ(t)v)[1:k]` is the correct value in \C^n.

IntrinsicSubspaceHomotopy(
    F::ModelKit.System,
    start::LinearSubspace,
    target::LinearSubspace;
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
) = IntrinsicSubspaceHomotopy(fixed(F; compile = compile), start, target)

function IntrinsicSubspaceHomotopy(
    system::AbstractSystem,
    start::LinearSubspace,
    target::LinearSubspace;
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
)
    IntrinsicSubspaceHomotopy(
        system,
        convert(LinearSubspace{ComplexF64}, start),
        convert(LinearSubspace{ComplexF64}, target),
    )
end

function IntrinsicSubspaceHomotopy(
    system::AbstractSystem,
    start::LinearSubspace{ComplexF64},
    target::LinearSubspace{ComplexF64},
)
    geodesic = get!(GEODESIC_LRU, (start, target)) do
        LinearSubspaceGeodesicInfo(start, target)
    end
    Q = geodesic.Q
    tx³ = TaylorVector{4}(ComplexF64, size(Q, 1))
    IntrinsicSubspaceHomotopy(
        system = system,
        start = start,
        target = target,
        geodesic = LinearSubspaceGeodesicInfo(start, target),
        J = zeros(ComplexF64, size(system) .+ (1, 1)),
        x = zeros(ComplexF64, size(Q, 1)),
        ẋ = zeros(ComplexF64, size(Q, 1)),
        x_high = zeros(ComplexDF64, size(Q, 1)),
        t_cache = Ref(complex(NaN, NaN)),
        taylor_t_cache = Ref(complex(NaN, NaN)),
        taylor_γ = tuple((similar(Q) for i = 0:3)...),
        v = zeros(ComplexF64, size(Q, 2)),
        tx³ = tx³,
        tx² = TaylorVector{3}(tx³),
        tx¹ = TaylorVector{2}(tx³),
    )
end
Base.size(H::IntrinsicSubspaceHomotopy) = (size(H.system)[1] + 1, dim(H.start) + 1)

"""
    set_subspaces!(H::IntrinsicSubspaceHomotopy, start::LinearSubspace, target::LinearSubspace)

Update the homotopy `H` to track from the affine subspace `start` to `target`.
"""
function set_subspaces!(
    H::IntrinsicSubspaceHomotopy,
    start::LinearSubspace,
    target::LinearSubspace,
)
    H.start = start
    H.target = target
    H.geodesic = get!(GEODESIC_LRU, (start, target)) do
        LinearSubspaceGeodesicInfo(start, target)
    end
    H
end
start_parameters!(H::IntrinsicSubspaceHomotopy, p::LinearSubspace) =
    set_subspaces!(H, p, H.target)
target_parameters!(H::IntrinsicSubspaceHomotopy, q::LinearSubspace) =
    set_subspaces!(H, H.starts, q)
parameters!(H::IntrinsicSubspaceHomotopy, p::LinearSubspace, q::LinearSubspace) =
    set_subspaces!(H, p, q)

function γ!(H::IntrinsicSubspaceHomotopy, t::Number)
    H.t_cache[] != t || return first(H.taylor_γ)
    if isreal(t)
        _γ!(H, real(t))
    else
        _γ!(H, t)
    end
    H.t_cache[] = t

    first(H.taylor_γ)
end
@inline function _γ!(H::IntrinsicSubspaceHomotopy, t::Number)
    @unpack Q, Q_cos, Θ = H.geodesic
    γ = first(H.taylor_γ)
    n, k = size(γ)
    @inbounds for j = 1:k
        Θⱼ = Θ[j]
        s, c = sincos(t * Θⱼ)
        for i = 1:n
            γ[i, j] = Q_cos[i, j] * c + Q[i, j] * s
        end
    end
    γ
end

γ̇!(H::IntrinsicSubspaceHomotopy, t::Number) = isreal(t) ? _γ̇!(H, real(t)) : _γ̇!(H, t)
@inline function _γ̇!(H::IntrinsicSubspaceHomotopy, t::Number)
    @unpack Q, Q_cos, Θ = H.geodesic
    _, γ̇ = H.taylor_γ
    n, k = size(γ̇)
    @inbounds for j = 1:k
        Θⱼ = Θ[j]
        s, c = sincos(t * Θⱼ)
        ċ = -s * Θⱼ
        ṡ = c * Θⱼ
        for i = 1:n
            γ̇[i, j] = Q_cos[i, j] * ċ + Q[i, j] * ṡ
        end
    end
    γ̇
end

function set_solution!(u::Vector, H::IntrinsicSubspaceHomotopy, x::AbstractVector, t)
    (length(x) == length(H.x) - 1) ||
        throw(ArgumentError("Cannot set solution. Expected extrinsic coordinates."))

    set_solution!(view(H.x, 1:length(x)), H.system, x)
    H.x[end] = 1

    if isone(t)
        LA.mul!(u, H.geodesic.γ1', H.x)
    elseif iszero(t)
        LA.mul!(u, H.geodesic.Q_cos', H.x)
    else
        LA.mul!(u, γ!(H, t)', H.x)
    end
end

function get_solution(H::IntrinsicSubspaceHomotopy, u::AbstractVector, t)
    if isone(t)
        (@view H.geodesic.γ1[1:end-1, :]) * u
    elseif iszero(t)
        (@view H.geodesic.Q_cos[1:end-1, :]) * u
    else
        γ = γ!(H, t)
        (@view γ[1:end-1, :]) * u
    end
end

function ModelKit.evaluate!(u, H::IntrinsicSubspaceHomotopy, v::Vector{ComplexDF64}, t)
    γ = γ!(H, t)
    n = first(size(H.system))
    LA.mul!(H.x_high, γ, v)
    evaluate!(u, H.system, H.x_high)
    u[n+1] = H.x_high[end] - 1.0
end

function ModelKit.evaluate!(u, H::IntrinsicSubspaceHomotopy, v::AbstractVector, t)
    γ = γ!(H, t)
    n = first(size(H.system))
    LA.mul!(H.x, γ, v)
    evaluate!(u, H.system, H.x)
    u[n+1] = H.x[end] - 1.0
    u
end

function ModelKit.evaluate_and_jacobian!(
    u,
    U,
    H::IntrinsicSubspaceHomotopy,
    v::AbstractVector,
    t,
)
    γ = γ!(H, t)

    LA.mul!(H.x, γ, v)
    evaluate_and_jacobian!(u, H.J, H.system, H.x)
    LA.mul!(U, H.J, γ)

    n = first(size(H.system))
    u[n+1] = H.x[end] - 1

    m = length(v)
    for j = 1:m
        U[n+1, j] = γ[end, j]
    end

    nothing
end

function ModelKit.taylor!(u, ::Val{1}, H::IntrinsicSubspaceHomotopy, v, t)
    γ = γ!(H, t)
    γ̇ = γ̇!(H, t)

    # apply chain rule
    #    d/dt [F(γ(t)v); (γ(t)v)[end] - 1] = [J_F(γ(t)v)* γ̇(t)*v;  (γ̇(t)v)[end]]
    LA.mul!(H.x, γ, v)
    LA.mul!(H.ẋ, γ̇, v)
    evaluate_and_jacobian!(u, H.J, H.system, H.x)
    LA.mul!(u, H.J, H.ẋ)
    M = size(H, 1)
    u[M] = H.ẋ[end]

    u
end

function _taylor_γ!(H::IntrinsicSubspaceHomotopy, t::Number)
    @unpack geodesic, taylor_γ = H
    @unpack Q, Q_cos, Θ, U = geodesic

    γ, γ¹, γ², γ³ = taylor_γ
    n, k = size(γ)
    @inbounds for j = 1:k
        Θⱼ = Θ[j]
        s, c = sincos(t * Θⱼ)
        c¹ = -s * Θⱼ
        s¹ = c * Θⱼ
        Θⱼ_2 = 0.5 * Θⱼ^2
        s² = -s * Θⱼ_2
        c² = -c * Θⱼ_2
        Θⱼ_3 = Θⱼ_2 * Θⱼ / 3
        s³ = -c * Θⱼ_3
        c³ = s * Θⱼ_3
        for i = 1:n
            γ[i, j] = Q_cos[i, j] * c + Q[i, j] * s
            γ¹[i, j] = Q_cos[i, j] * c¹ + Q[i, j] * s¹
            γ²[i, j] = Q_cos[i, j] * c² + Q[i, j] * s²
            γ³[i, j] = Q_cos[i, j] * c³ + Q[i, j] * s³
        end
    end

end

function taylor_γ!(H::IntrinsicSubspaceHomotopy, t::Number)
    H.taylor_t_cache[] != t || return H.taylor_γ

    if isreal(t)
        _taylor_γ!(H, real(t))
    else
        _taylor_γ!(H, t)
    end
    H.taylor_t_cache[] = t

    H.taylor_γ
end

function ModelKit.taylor!(
    u,
    v::Val{2},
    H::IntrinsicSubspaceHomotopy,
    tv::TaylorVector,
    t,
    incr::Bool,
)
    γ, γ¹, γ², γ³ = taylor_γ!(H, t)
    x, x¹, x² = vectors(H.tx²)
    v, v¹ = vectors(tv)

    if incr
        x .= H.x
        x¹ .= H.ẋ
    else
        LA.mul!(x, γ, v)
        LA.mul!(x¹, γ¹, v)
    end

    H.v .= v¹
    LA.mul!(H.x, γ, v¹)
    x¹ .+= H.x

    LA.mul!(H.x, γ¹, H.v)
    H.v .= v
    LA.mul!(H.x, γ², H.v, true, true)
    x² .= H.x

    taylor!(u, Val(2), H.system, H.tx²)

    n = first(size(H.system))
    u[n+1] = x²[end]

    u
end

function ModelKit.taylor!(
    u,
    v::Val{3},
    H::IntrinsicSubspaceHomotopy,
    tv::TaylorVector,
    t,
    incr::Bool,
)
    γ, γ¹, γ², γ³ = taylor_γ!(H, t)
    x, x¹, x², x³ = vectors(H.tx³)
    v, v¹, v² = vectors(tv)

    if !incr
        LA.mul!(x, γ, v)
        LA.mul!(x¹, γ¹, v)
        LA.mul!(x¹, γ, v¹, true, true)
        LA.mul!(x², γ², v)
        LA.mul!(x², γ¹, v¹, true, true)
    end

    H.v .= v²
    LA.mul!(H.x, γ, H.v)
    x² .+= H.x

    LA.mul!(H.x, γ¹, H.v)
    H.v .= v
    LA.mul!(H.x, γ³, H.v, true, true)
    H.v .= v¹
    LA.mul!(H.x, γ², H.v, true, true)
    x³ .= H.x

    taylor!(u, Val(3), H.system, H.tx³)

    n = first(size(H.system))
    u[n+1] = x³[end]

    u
end


text:
#
module MultTest

export X
export ptest43, ptest51, ptest52, ptest53

import Polynomials
using Polynomials

# Standard polynomial
X = poly([0.0])
#some test polynomials used by Zeng

# polynomial with roots 1:4 to the power of k
# k = 1..8
function ptest46(k::Int)
    ((X-1)^4 * (X-2)^3 * (X-3)^2 * (X-4) ) ^k
end

# polynomial with non-integer roots
function ptest51(dig::Int)
   p = (X - 10/11)^5 * (X - 20/11)^5 * (X - 30/11)^5
   Poly(map(i->round(p.a[i], dig), 1:length(p.a)))
end

# polynomial with adjacent roots close to 1 with given multiplicity
# epsilon = 0.1, 0.01, 0.001 ...
function ptest52(epsilon::Number)
    (X - 1 + epsilon)^20 * (X - 1)^20 * (X + 0.5)^5
end

# polynomial with given roots (including all conjugate roots) to the power of 2^k
function ptest53(k::Int)
	c = [0.5+im, -1.0+0.2im, -0.1+im, +0.8+0.6im, -0.7+0.7im,
		 1.4, -0.4+0.9im, 0.9, -0.8+0.3im, 0.3+0.8im, 0.6+0.4im]

    fX = Poly([1.0])
    for i = 1:length(c)
        a = c[i]
        if isreal(a)
            fX = fX * (X - a)
        else
            fX = fX * ((X-real(a))^2 + imag(a)^2)
        end
    end

	for i = 1:k
        fX = fX * fX
    end
    fX
end

end

import MultTest
using MultTest

ext: jl


text:
module TestRDA
    using Base.Test
    using DataArrays
    using DataFrames

    # R code generating test .rdas
        # df = data.frame(num=c(1.1, 2.2))
        # save(df, file='minimal.rda')

        # df['int'] = c(1L, 2L)
        # df['logi'] = c(T, F)
        # df['chr'] = c('ab', 'c')
        # df['factor'] = factor(df[['chr']])
        # #utf=c('Ж', '∰')) R handles it, read_rda doesn't.
        # save(df, file='types.rda')

        # df[2, ] = NA
        # df['chr'] = NULL  # NA characters breaking read_rda
        # save(df, file='NAs.rda')

    df = DataFrame(num = [1.1, 2.2])
    @assert isequal(DataFrame(read_rda("test/data/RDA/minimal.rda")["df"]), df)

    df["int"] = Int32[1, 2]
    df["logi"] = [true, false]
    df["chr"] = ["ab", "c"]
    df["factor"] = pool(df["chr"])
    @assert isequal(DataFrame(read_rda("test/data/RDA/types.rda")["df"]), df)

    df[2, :] = NA
    df = df[:, ["num", "int", "logi", "factor"]]  # (NA) chr breaks read_rda
    @assert isequal(DataFrame(read_rda("test/data/RDA/NAs.rda")["df"]), df)
end


text:
@testset "Lennard-Jones Unit Tests" begin
    ϵ = 1.0u"eV"
    σ = 0.25u"bohr"
    rcutoff = 2.0u"Å"
    species = [:Ar, :H]
    p = LennardJones(ϵ, σ, rcutoff, species)

    @test p isa EmpiricalPotential{NamedTuple{(:ϵ, :σ)},NamedTuple{(:rcutoff,)}}

    @test get_rcutoff(p) == austrip(rcutoff)
    @test get_species(p) == (:Ar, :H)

    @test get_parameters(p) == (; ϵ=austrip(ϵ), σ=austrip(σ))
    @test set_parameters(p, (; ϵ=2.0, σ=3.0)) == LennardJones(2.0, 3.0, austrip(rcutoff), species)
    @test serialize_parameters(p) == [austrip(ϵ), austrip(σ)]
    @test deserialize_parameters(p, [2.0, 3.0]) == LennardJones(2.0, 3.0, austrip(rcutoff), species)

    @test get_hyperparameters(p) == (; rcutoff=austrip(rcutoff))
    @test set_hyperparameters(p, (; rcutoff=1.0)) == LennardJones(austrip(ϵ), austrip(σ), 1.0, species)
    @test serialize_hyperparameters(p) == [austrip(rcutoff)]
    @test deserialize_hyperparameters(p, [1.0]) == LennardJones(austrip(ϵ), austrip(σ), 1.0, species)

    r = @SVector[1.0, 1.0, 1.0]
    R = norm(r)

    @test potential_energy(R, p) isa Float64
    @test force(R, r, p) isa SVector{3,Float64}
end


text:
using NewtonsMethod
using Test

@testset "NewtonsMethod.jl" begin
   testf(x)= (x-1)^3
    testf_prime(x)=3*(x-1)^2
    testf1(x)= exp(3*x)-1
    testf1_prime(x)=3*exp(3*x)

    #test for functions with derivative

    root1, conv, iter  = NewtonsMethod.newtonroot(testf,testf_prime; x_0 = 0.0)
    root2, conv, iter  = NewtonsMethod.newtonroot(testf1,testf1_prime; x_0 = 0.0)

    @test abs(root1 - 1.0)<1e-6
    @test abs(root2 - 0.0)<1e-6

    #do the same for functions without derivative

    root1, conv, iter  = NewtonsMethod.newtonroot(testf; x_0 = 0.0)
    root2, conv, iter  = NewtonsMethod.newtonroot(testf1; x_0 = 0.0)

    @test abs(root1 - 1.0) < 1e-6
    @test abs(root2 - 0.0) < 1e-6

    #Test if maxiter works
    root1, conv, iter  = NewtonsMethod.newtonroot(testf, testf_prime; x_0 = 0.0,  maxiter = 5)
    # if the function converged, return number of iteration, else return nothing
    @test iter == nothing ||iter <= 5
    #Test different tolerance, increasing it to 1e-10
    root1, conv, iter  = NewtonsMethod.newtonroot(testf, testf_prime; x_0=0.0,  tol = 1e-10)
    @test abs(root1-1.0)<1e-7

    #Test lack of convergence

    test_nc(x) = x^2+1

    rootnc, conv, iter  = NewtonsMethod.newtonroot(test_nc; x_0 = 1.0)
    @test rootnc==nothing

    x0 = BigFloat(0.79)
    root1, conv, iter  = NewtonsMethod.newtonroot(testf; x_0 = x0)
    root2, conv, iter  = NewtonsMethod.newtonroot(testf1; x_0 = x0)

    @test abs(root1 - 1.0)<1e-6
    @test abs(root2 - 0.0)<1e-6

end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_antdsunburst

"""
    ''_antdsunburst(;kwargs...)

An AntdSunburst component.

Keyword arguments:
- `id` (String; optional)
- `annotations` (optional)
- `appendPadding` (Real | Array of Reals; optional)
- `autoFit` (Bool; optional)
- `className` (String; optional)
- `color` (optional): . color has the following type: String | Array of Strings | lists containing elements 'func'.
Those elements have the following types:
  - `func` (String; optional)
- `colorField` (String; optional)
- `data` (Dict; required)
- `drilldown` (optional): . drilldown has the following type: lists containing elements 'enabled', 'breadCrumb'.
Those elements have the following types:
  - `enabled` (Bool; optional)
  - `breadCrumb` (optional): . breadCrumb has the following type: lists containing elements 'rootText', 'dividerText', 'textStyle', 'activeTextStyle', 'position'.
Those elements have the following types:
  - `rootText` (String; optional)
  - `dividerText` (String; optional)
  - `textStyle` (optional)
  - `activeTextStyle` (optional)
  - `position` (a value equal to: 'top-left', 'bottom-left'; optional)
- `height` (Real; optional)
- `hierarchyConfig` (optional): . hierarchyConfig has the following type: lists containing elements 'field', 'ignoreParentValue'.
Those elements have the following types:
  - `field` (String; optional)
  - `ignoreParentValue` (Bool; optional)
- `innerRadius` (Real; optional)
- `label` (optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `locale` (a value equal to: 'zh-CN', 'en-US'; optional)
- `meta` (optional)
- `padding` (Real | Array of Reals | a value equal to: 'auto'; optional)
- `radius` (Real; optional)
- `rawFields` (Array of Strings; optional)
- `reflect` (Bool; optional)
- `renderer` (a value equal to: 'canvas', 'svg'; optional)
- `style` (Dict; optional)
- `sunburstStyle` (optional): . sunburstStyle has the following type: lists containing elements 'func'.
Those elements have the following types:
  - `func` (String; optional)
- `tooltip` (optional)
- `width` (Real; optional)
"""
function ''_antdsunburst(; kwargs...)
        available_props = Symbol[:id, :annotations, :appendPadding, :autoFit, :className, :color, :colorField, :data, :drilldown, :height, :hierarchyConfig, :innerRadius, :label, :loading_state, :locale, :meta, :padding, :radius, :rawFields, :reflect, :renderer, :style, :sunburstStyle, :tooltip, :width]
        wild_props = Symbol[]
        return Component("''_antdsunburst", "AntdSunburst", "feffery_antd_charts", available_props, wild_props; kwargs...)
end

max_stars_repo_name: CNFeffery/feffery-antd-charts

max_issues_repo_name: CNFeffery/feffery-antd-charts

max_forks_repo_name: CNFeffery/feffery-antd-charts


text:
using IRTools, Test
using IRTools: IR, @dynamo

function f(x)
    N = 4
    for i1 in 1:3    # single loop works without issue
        for i2 in 1:N   # needs to be a variable `1:4` works fine
        end
    end
    0.0 # same error with `x`
end

@test @code_ir(f(1)) isa IR

function f(a, b)
    u = 1
    while true
        if true
        end
    end
    f(u)
end

@test @code_ir(f(1, 2)) isa IR

# issue 30
@dynamo function foo(a...)
    ir = IR(a...)
    return ir
end

mylog2(x) = ccall((:log2, Base.Math.libm), Float64, (Float64,), x)

@test foo(mylog2, 3.3) === mylog2(3.3)

bar(a, b) = a > b ? a : b

function foo(a, b)
  bar(a, b)
end

let
  ir1 = @code_ir foo(1, 2)
  ir2 = @code_ir bar(1, 2)
  ir = IRTools.inline(ir1, IRTools.var(4), ir2)
  f = IRTools.func(ir)
  @test f(nothing, 2, 3) == 3
  @test f(nothing, 3, 2) == 3
end

function foo1(x)
  while x < 10
    x += 1
  end
  return x
end

function foo2(x)
  if x < 100
    x + foo1(x)
  else
    x
  end
end

let
  ir = @code_ir foo2(1)
  ir2 = @code_ir foo1(1)
  ir3 = IRTools.inline(ir, IRTools.var(4), ir2)
  @test IRTools.func(ir3)(nothing, 2) == 12
  @test IRTools.func(ir3)(nothing, 101) == 101
end

function foo1(x)
  while x < 10
    x += 1
  end
  return x
end

function foo2(x)
  foo1(x)
  x > 0 ? x : 0
end

let
  ir = @code_ir foo2(1)
  ir2 = @code_ir foo1(1)
  ir3 = IRTools.inline(ir, IRTools.var(3), ir2)
  @test IRTools.func(ir3)(nothing, 2) == foo2(2)
  @test IRTools.func(ir3)(nothing, -2) == foo2(-2)
end


text:
using Dates: Hour, Minute, Second, Millisecond, days, hour, minute, second, millisecond

"""
    localtime(::ZonedDateTime) -> DateTime

Creates a local or civil `DateTime` from the given `ZonedDateTime`. For example the
`2014-05-30T08:11:24-04:00` would return `2014-05-30T08:11:24`.
"""
localtime(zdt::ZonedDateTime) = zdt.utc_datetime + zdt.zone.offset

"""
    utc(::ZonedDateTime) -> DateTime

Creates a utc `DateTime` from the given `ZonedDateTime`. For example the
`2014-05-30T08:11:24-04:00` would return `2014-05-30T12:11:24`.
"""
utc(zdt::ZonedDateTime) = zdt.utc_datetime

"""
    timezone(::ZonedDateTime) -> TimeZone

Returns the `TimeZone` used by the `ZonedDateTime`.
"""
timezone(zdt::ZonedDateTime) = zdt.timezone

Dates.days(zdt::ZonedDateTime) = days(localtime(zdt))

for period in (:Hour, :Minute, :Second, :Millisecond)
    accessor = Symbol(lowercase(string(period)))
    @eval begin
        Dates.$accessor(zdt::ZonedDateTime) = $accessor(localtime(zdt))
        Dates.$period(zdt::ZonedDateTime) = $period($accessor(zdt))
    end
end

Base.eps(::ZonedDateTime) = Millisecond(1)


text:
module Topologies

export Topology, TopologyIndex, coordinates, is_valid, index_type

abstract type Topology end
abstract type TopologyIndex end

# TODO: specify types

function coordinates end

function is_valid end

function index_type end

include("square_greed_topologies.jl")

include("hex_greed_topologies.jl")

end

max_stars_repo_name: Tiendil/pcg-julia

max_issues_repo_name: Tiendil/pcg-julia

max_forks_repo_name: Tiendil/pcg-julia


text:
#= Read and write Pajek files

Most of things specified in the format (like node and edge attributes) are supported.
See [1](http://www.pfeffer.at/txt2pajek/txt2pajek.pdf), and
[2](http://vlado.fmf.uni-lj.si/pub/networks/pajek/doc/draweps.htm)
for the reference. =#

const PAJEK_ATTRIBUTES = Dict("ic"=>:color, "c"=>:color, "w"=>:width, "l"=>:label,
    "bw"=>:border_width, "fos"=>:font_size, "lr"=>:label_angle, "s"=>:arrow_size)

const PAJEK_ATTR_DEFAULTS = Dict(:label => "", :shape => :circle,
    :color => "Grey", :weight => 1., :width => 1.)

const PAJEK_SHAPES = Dict("box" => :square, "triangle" => :triangle,
    "ellipse" => :circle, "diamond" => :diamond)

function _skip_to_header(io::IO)
    for line in eachline(io)
        line[1] == '*' && return line
    end
end

function _is_directed_lookahead(io::IO)
    pos = position(io)
    for line in eachline(io)
        if line[1] == '*' && contains(line, "Arcs")
            seek(io, pos)
            return true
        end
    end
    seek(io, pos)
    return false
end

function _parse_nodes(io::IO, g::Graph, header::String, parse_attrs::Bool)
    n = parse(Int, split(header)[end])
    addnodes!(g, n)
    if !parse_attrs
        return _skip_to_header(io)
    else
        for line in eachline(io)
            line[1] == '*' && return line

            tokens = split(line)
            ind = parse(Int, shift!(tokens))
            attrs = Dict{Symbol,Any}(:label => String(strip(shift!(tokens), '"')))

            i = 1
            while length(tokens) > 0
                token = shift!(tokens)
                if i < 4 && !isalpha(token[1])  # coordX, coordY, value, and/or shape
                    attrs[(:x, :y, :size)[i]] = float(token)
                elseif token in ("box", "triangle", "ellipse", "diamond")
                    attrs[:shape] = PAJEK_SHAPES[token]
                elseif token == "ic"            # color
                    if length(tokens) > 0
                        attrs[:color] = String(shift!(tokens))
                    end
                elseif length(tokens) > 0
                    attrs[Symbol(token)] = float(shift!(tokens))
                end
                i += 1
            end

            for (k, v) = attrs
                if !hasnodeattr(g, k)
                    default = get(PAJEK_ATTR_DEFAULTS, k, 0.)
                    addnodeattr!(g, k, default)
                end
                g.nodeattrs[k, ind] = v
            end
        end
    end
end

function _parse_edges(io::IO, g::Graph, directed::Bool, parse_attrs::Bool)
    for line in eachline(io)
        line[1] == '*' && return line

        tokens = split(line)
        x, y = parse(Int, shift!(tokens)), parse(Int, shift!(tokens))
        attrs = Dict{Symbol,Any}()

        if parse_attrs
            i = 1
            while length(tokens) > 0
                t = shift!(tokens)
                if i == 1 && isnumber(t[1])
                    attrs[:weight] = float(t)
                elseif length(tokens) > 0
                    if t == "l"
                        attrs[:label] = String(strip(shift!(tokens), '"'))
                    elseif t == "w"
                        attrs[:width] = float(shift!(tokens))
                    elseif t == "c"
                        attrs[:color] = String(shift!(tokens))
                    else
                        t_ = shift!(tokens)
                        attrs[Symbol(t)] = t_[1] == '"' ? String(strip(t_, '"')) : float(t_)
                    end
                end
                i += 1
            end
            for k = keys(attrs)
                if !hasedgeattr(g, k)
                    default = get(PAJEK_ATTR_DEFAULTS, k, 0.)
                    addedgeattr!(g, k, default)
                end
            end
        end

        addedge!(g, x, y; attrs...)
        isdirected(g) && !directed && addedge!(g, y, x; attrs...)
    end
end

"""
    read_pajek(io::IO[; parse_attrs=true])

Read graph from a Pajek-formatted file.

Since Pajek format is not very well documented, there may be errors parsing
some files. Most of the errors will come from parsing the node and edge attributes.
To suppress it, pass `parse_attrs=false` to this function, which also
dramatically speeds up the parsing.
"""
function read_pajek(io::IO; parse_attrs=true)
    g = Graph(directed=_is_directed_lookahead(io))
    header = _skip_to_header(io)
    while !eof(io)
        header = lowercase(header)
        if contains(header, "vertices")
            header = _parse_nodes(io, g, header, parse_attrs)
        elseif contains(header, "arcs")
            header = _parse_edges(io, g, true, parse_attrs)
        else # "edges"
            header = _parse_edges(io, g, false, parse_attrs)
        end
    end
    return g
end

"""
    writepajek(filename, g::Graph)

Serialize the `Graph` instance `g` in Pajek format.
"""
function writepajek(filename, g)
    f = open(filename, "w")
    # Nodes
    println(f, "*Vertices $(nodecount(g))")
    if hasnodeattr(g, :label)
        for (i, l) = enumerate(g[:, :label])
            println(f, i, '"', l, '"')
        end
    end
    # Edges
    println(f, isdirected(g) ? "*Arcs" : "*Edges")
    for e = edges(g)
        println(f, e.source, ' ', e.target)
    end
    close(f)
end


text:
function moi_sdplib(optimizer, path; verbose = false, test = false, scalar = false)

    if verbose
        println("running: $(path)")
    end
    MOI.empty!(optimizer)
    if test
        @test MOI.is_empty(optimizer)
    end

    n, m, F, c = sdplib_data(path)

    nvars = sympackedlen(n)

    X = MOI.add_variables(optimizer, nvars)
    vov = MOI.VectorOfVariables(X)
    cX = MOI.add_constraint(optimizer, vov, MOI.PositiveSemidefiniteConeTriangle(n))

    Xsq = Matrix{MOI.VariableIndex}(undef, n,n)
    ivech!(Xsq, X)
    Xsq = Matrix(Symmetric(Xsq,:U))

    # Objective function
    objf_t = [MOI.ScalarAffineTerm(F[0][idx...], Xsq[idx...])
        for idx in zip(findnz(F[0])[1:end-1]...)]
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(objf_t, 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    # Linear equality constraints
    for k in 1:m
        ctr_k = [MOI.ScalarAffineTerm(F[k][idx...], Xsq[idx...])
            for idx in zip(findnz(F[k])[1:end-1]...)]
        if scalar
            MOI.add_constraint(optimizer, MOI.ScalarAffineFunction(ctr_k, 0.0), MOI.EqualTo(c[k]))
        else
            MOI.add_constraint(optimizer,
                MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1], ctr_k), [-c[k]]), MOI.Zeros(1))
        end
    end

    MOI.optimize!(optimizer)

    objval = MOI.get(optimizer, MOI.ObjectiveValue())

    stime = -1.0
    try
        stime = MOI.get(optimizer, MOI.SolveTime())
    catch
        println("could not query time")
    end

    Xsq_s = MOI.get.(optimizer, MOI.VariablePrimal(), Xsq)
    minus_rank = length([eig for eig in eigen(Xsq_s).values if eig < -1e-4])
    if test
        @test minus_rank == 0
    end
    # @test tr(F[0] * Xsq_s) - obj < 1e-1
    # for i in 1:m
    #     @test abs(tr(F[i] * Xsq_s)-c[i]) < 1e-1
    # end

    verbose && sdplib_eval(F,c,n,m,Xsq_s)

    rank = -1
    status = 0
    if MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
        status = 1
    end
    return (objval, stime, rank, status)
end


text:
@testset "type qualifiers" begin
	c``

	@eval c"""
	struct CconstStruct {
		int i;
		int const ci;

		struct {
			int i;
			int const ci;
		} s;

		struct {
			int i;
			int const ci;
		} const cs;

		union {
			int i;
			int const ci;
		} u;

		union {
			int i;
			int const ci;
		} const cu;

		struct {
			int i;
			int const ci;
		} as[1];

		struct {
			int i;
			int const ci;
		} const acs[1];

		struct {
			int si;
			int const sci;
		};

		union {
			int ui;
			int const uci;
		};
	};
	"""
	@test sizeof(c"struct CconstStruct") == sizeof(Cint)*(2*6+3)

	x = c"struct CconstStruct"()
	@test x.i isa Cint
	@test x.ci isa Cint
	@test x.s isa Cstruct
	@test x.s.i isa Cint
	@test x.s.ci isa Cint
	@test x.cs isa Cstruct
	@test x.cs.i isa Cint
	@test x.cs.ci isa Cint
	@test x.u isa Cunion
	@test x.u.i isa Cint
	@test x.u.ci isa Cint
	@test x.cu isa Cunion
	@test x.cu.i isa Cint
	@test x.cu.ci isa Cint
	@test x.as isa Carray{<:Cstruct, 1}
	@test x.as[1] isa Cstruct
	@test x.as[1].i isa Cint
	@test x.as[1].ci isa Cint
	@test x.acs isa Carray{<:Cstruct, 1}
	@test x.acs[1] isa Cstruct
	@test x.acs[1].i isa Cint
	@test x.acs[1].ci isa Cint
	@test x.si isa Cint
	@test x.sci isa Cint
	@test x.ui isa Cint
	@test x.uci isa Cint

	ptr = Libc.malloc(x)
	@test eltype(ptr.i) <: Cint
	@test eltype(ptr.ci) <: Cconst{Cint}
	@test eltype(ptr.s) <: Cstruct
	@test eltype(ptr.s.i) <: Cint
	@test eltype(ptr.s.ci) <: Cconst{Cint}
	@test eltype(ptr.cs) <: Cconst{<:Cstruct}
	@test eltype(ptr.cs.i) <: Cconst{Cint}
	@test eltype(ptr.cs.ci) <: Cconst{Cint}
	@test eltype(ptr.u) <: Cunion
	@test eltype(ptr.u.i) <: Cint
	@test eltype(ptr.u.ci) <: Cconst{Cint}
	@test eltype(ptr.cu) <: Cconst{<:Cunion}
	@test eltype(ptr.cu.i) <: Cconst{Cint}
	@test eltype(ptr.cu.ci) <: Cconst{Cint}
	@test eltype(ptr.as) <: Carray{<:Cstruct, 1}
	@test eltype(ptr.as[1]) <: Cstruct
	@test eltype(ptr.as[1].i) <: Cint
	@test eltype(ptr.as[1].ci) <: Cconst{Cint}
	@test eltype(ptr.acs) <: Carray{<:Cconst{<:Cstruct}, 1}
	@test eltype(ptr.acs[1]) <: Cconst{<:Cstruct}
	@test eltype(ptr.acs[1].i) <: Cconst{Cint}
	@test eltype(ptr.acs[1].ci) <: Cconst{Cint}
	@test eltype(ptr.si) <: Cint
	@test eltype(ptr.sci) <: Cconst{Cint}
	@test eltype(ptr.ui) <: Cint
	@test eltype(ptr.uci) <: Cconst{Cint}

	cptr = reinterpret(Cptr{Cconst{eltype(ptr)}}, ptr)
	@test eltype(cptr.i) <: Cconst{Cint}
	@test eltype(cptr.ci) <: Cconst{Cint}
	@test eltype(cptr.s) <: Cconst{<:Cstruct}
	@test eltype(cptr.s.i) <: Cconst{Cint}
	@test eltype(cptr.s.ci) <: Cconst{Cint}
	@test eltype(cptr.cs) <: Cconst{<:Cstruct}
	@test eltype(cptr.cs.i) <: Cconst{Cint}
	@test eltype(cptr.cs.ci) <: Cconst{Cint}
	@test eltype(cptr.u) <: Cconst{<:Cunion}
	@test eltype(cptr.u.i) <: Cconst{Cint}
	@test eltype(cptr.u.ci) <: Cconst{Cint}
	@test eltype(cptr.cu) <: Cconst{<:Cunion}
	@test eltype(cptr.cu.i) <: Cconst{Cint}
	@test eltype(cptr.cu.ci) <: Cconst{Cint}
	@test eltype(cptr.as) <: Cconst{<:Carray{<:Cstruct, 1}}
	@test eltype(cptr.as[1]) <: Cconst{<:Cstruct}
	@test eltype(cptr.as[1].i) <: Cconst{Cint}
	@test eltype(cptr.as[1].ci) <: Cconst{Cint}
	@test eltype(cptr.acs) <: Cconst{<:Carray{<:Cconst{<:Cstruct}, 1}}
	@test eltype(cptr.acs[1]) <: Cconst{<:Cstruct}
	@test eltype(cptr.acs[1].i) <: Cconst{Cint}
	@test eltype(cptr.acs[1].ci) <: Cconst{Cint}
	@test eltype(cptr.si) <: Cconst{Cint}
	@test eltype(cptr.sci) <: Cconst{Cint}
	@test eltype(cptr.ui) <: Cconst{Cint}
	@test eltype(cptr.uci) <: Cconst{Cint}

	@test ptr.i+1 == ptr.ci
	@test ptr.s+1 == ptr.cs
	@test ptr.i+6 == ptr.u.i
	@test ptr.u.i-6 == ptr.i

	f(p) =
		p.i[] + p.ci[] +
		p.s.i[] + p.s.ci[] +
		p.cs.i[] + p.cs.ci[] +
		p.u.i[] + p.u.ci[] +
		p.cu.i[] + p.cu.ci[] +
		p.as[1].i[] + p.as[1].ci[] +
		p.acs[1].i[] + p.acs[1].ci[] +
		p.si[] + p.sci[] +
		p.ui[] + p.uci[]
	@test f(ptr) == 0
	@test @allocated(f(ptr)) == 0

	@test ptr.i[] == 0
	ptr.i[] = 1
	@test ptr.i[] == 1

	@test ptr.ci[] == 0
	@test_throws ErrorException ptr.ci[] = 1
	@test ptr.ci[] == 0

	@test ptr.s[] == eltype(ptr.s)()
	@test ptr.s.i[] == 0
	@test ptr.s.ci[] == 0
	ptr.s[] = eltype(ptr.s)(i = 2, ci = 3)

	@test ptr.s.i[] == 2
	ptr.s.i[] = 22
	@test ptr.s.i[] == 22
	@test_throws ErrorException ptr.s[].i = 222
	@test ptr.s.i[] == 22

	@test ptr.s.ci[] == 3
	@test_throws ErrorException ptr.s.ci[] = 33
	@test ptr.s.ci[] == 3
	@test_throws ErrorException ptr.s[].ci = 333
	@test ptr.s.ci[] == 3

	@test ptr.cs[] == unqualifiedtype(eltype(ptr.cs))()
	@test ptr.cs.i[] == 0
	@test ptr.cs.ci[] == 0
	@test_throws ErrorException ptr.cs[] = unqualifiedtype(eltype(ptr.cs))(i = 4, ci = 5)

	@test ptr.cs.i[] == 0
	@test_throws ErrorException ptr.cs.i[] = 44
	@test ptr.cs.i[] == 0
	@test_throws ErrorException ptr.cs[].i = 444
	@test ptr.cs.i[] == 0

	@test ptr.cs.ci[] == 0
	@test_throws ErrorException ptr.cs.ci[] = 55
	@test ptr.cs.ci[] == 0
	@test_throws ErrorException ptr.cs[].ci = 555
	@test ptr.cs.ci[] == 0

	@test ptr.u[] == eltype(ptr.u)()

	@test ptr.u.i[] == 0
	ptr.u[] = eltype(ptr.u)(i = 6)
	@test ptr.u.i[] == 6
	ptr.u.i[] = 66
	@test ptr.u.i[] == 66
	ptr.u.i[] = 666
	@test ptr.u.i[] == 666

	@test ptr.u.ci[] == 666
	ptr.u[] = eltype(ptr.u)(ci = 7)
	@test ptr.u.ci[] == 7
	@test_throws ErrorException ptr.u.ci[] = 77
	@test ptr.u.ci[] == 7
	@test_throws ErrorException ptr.u[].ci = 777
	@test ptr.u.ci[] == 7

	@test ptr.cu[] == unqualifiedtype(eltype(ptr.cu))()

	@test ptr.cu.i[] == 0
	@test_throws ErrorException ptr.cu[] = unqualifiedtype(eltype(ptr.cu))(i = 8)
	@test ptr.cu.i[] == 0
	@test_throws ErrorException ptr.cu.i[] = 88
	@test ptr.cu.i[] == 0
	@test_throws ErrorException ptr.cu.i[] = 888
	@test ptr.cu.i[] == 0

	@test ptr.cu.ci[] == 0
	@test_throws ErrorException ptr.cu[] = unqualifiedtype(eltype(ptr.cu))(ci = 9)
	@test ptr.cu.ci[] == 0
	@test_throws ErrorException ptr.cu.ci[] = 99
	@test ptr.cu.ci[] == 0
	@test_throws ErrorException ptr.cu[].ci = 999
	@test ptr.cu.ci[] == 0

	@test ptr.si[] == 0
	ptr.si[] = 1
	@test ptr.si[] == 1

	@test ptr.sci[] == 0
	@test_throws ErrorException ptr.sci[] = 2
	@test ptr.sci[] == 0

	@test ptr.ui[] == 0
	ptr.ui[] = 5
	@test ptr.ui[] == 5

	@test ptr.uci[] == 5
	@test_throws ErrorException ptr.uci[] = 6
	@test ptr.uci[] == 5

	@test f(ptr) == 51
	@test f(cptr) == 51
	@test_throws ErrorException cptr[] = x
	ptr[] = x
	@test f(cptr) == 0
	@test @allocated(f(cptr)) == 0

	@test cptr.i[] == 0
	@test_throws ErrorException cptr.i[] = 1
	@test cptr.i[] == 0

	@test cptr.s.i[] == 0
	@test_throws ErrorException cptr.s.i[] = 2
	@test cptr.s.i[] == 0

	@test cptr.u.i[] == 0
	@test_throws ErrorException cptr.u.i[] = 3
	@test cptr.u.i[] == 0

	@test cptr.si[] == 0
	@test_throws ErrorException cptr.si[] = 4
	@test cptr.si[] == 0

	@test cptr.ui[] == 0
	@test_throws ErrorException cptr.ui[] = 5
	@test cptr.ui[] == 0

	@test f(cptr) == 0
	@test @allocated(f(cptr)) == 0

	@eval c"""
	struct SelfRef {
		struct SelfRef *ptr;
	};
	"""
	@test sizeof(c"struct SelfRef") == sizeof(Cptr)

	Libc.free(ptr)
end


text:
# BBI Header
# ==========

const BED_MAGIC = 0x8789F2EB
const WIG_MAGIC = 0x888FFC26

# Supplemental Table 5.
struct Header
    magic::UInt32
    version::UInt16
    zoom_levels::UInt16
    chromosome_tree_offset::UInt64
    full_data_offset::UInt64
    full_index_offset::UInt64
    field_count::UInt16
    defined_field_count::UInt16
    auto_sql_offset::UInt64
    total_summary_offset::UInt64
    uncompress_buf_size::UInt32
    reserved::UInt64
end

const HEADER_SIZE = 64

function Base.read(io::IO, ::Type{Header})
    return Header(
        read(io, UInt32), read(io, UInt16), read(io, UInt16),
        read(io, UInt64), read(io, UInt64), read(io, UInt64),
        read(io, UInt16), read(io, UInt16), read(io, UInt64),
        read(io, UInt64), read(io, UInt32), read(io, UInt64))
end

function Base.write(stream::IO, header::Header)
    return write(
        stream,
        header.magic, header.version, header.zoom_levels,
        header.chromosome_tree_offset, header.full_data_offset,
        header.full_index_offset, header.field_count,
        header.defined_field_count, header.auto_sql_offset,
        header.total_summary_offset, header.uncompress_buf_size,
        header.reserved)
end


text:
module ClimateBase
export NCDataset

include("core/coredefs.jl")
include("core/prettyprint.jl")
include("core/nc_io.jl")
include("core/aggregation.jl")
include("core/interpolation.jl")

include("physical_dimensions/spatial.jl")
include("physical_dimensions/spatial_equalarea.jl")
include("physical_dimensions/temporal.jl")

include("climate/solar.jl")
include("climate/albedo.jl")

include("tsa/continuation.jl")
include("tsa/decomposition.jl")

end # module


text:
@doc raw"""
    function replay(f, t :: T) where T <: Trace

Given a stochastic function `f` and a trace `t`, makes `sample` calls behave as though
they had sampled the values in `t` at the corresponding addresses.

Returns a tuple `(t_new, g)`, where `t_new` is a trace and `g` is a function.
The function signature of `g` is the same as that of `f` with the first argument removed; that
is, if `f(t :: Trace, params...)`, then `g(params...)`. Computation is delayed, so each of the
latent nodes in `t_new` has `interpretation = REPLAYED`.
Calling `g(params...)` executes the computation and each latent node in `t_new` reverts to
its original interpretation.
"""
function replay(f, t :: T) where T <: Trace
    t_new = deepcopy(t)
    interpret_latent!(t_new, REPLAYED)
    g(params...) = f(t_new, params...)
    (t_new, g)
end

@doc raw"""
    function replace(f, r :: Dict)

Given a mapping `r` from addresses to distribution-like (currently `Distributions`
or `Array{Any, 1}`s), replaces the current distributions at that set of addresses
with this set of distributions. Returns a function `g` that has return signature
`(t :: Trace, rval)` where `rval` is a return value of `f`.
"""
function replace(f, r :: Dict)
    function g(t :: Trace, params...)
        rval = f(t, params...)
        t = replace(t, r)
        (t, rval)
    end
    g
end

@doc raw"""
    function replace(t :: Trace, r :: Dict)

Given a mapping `r` from addresses to distribution-like (currently `Distributions`
or `Array{Any, 1}`s), replaces the current distributions at that set of addresses
with this set of distributions. Returns the modified trace.
"""
function replace(t :: T, r :: Dict) where T <: Trace
    for (a, d) in r
        old_n = deepcopy(t[a])
        t[a] = node(old_n.value, old_n.address, d, old_n.observed, old_n.interpretation)
    end
    t
end

@doc raw"""
    function rewrite(f, t :: T, r :: Dict) where T <: Trace

Rewrites the history of the trace to make it appear as if the values in the trace
were sampled at the addresses in the keys of `r` from the corresponding distributions in the values of `r`.
Returns a function with call signature `g(params...)` that returns `(t :: Trace, rval)`, where `rval`
is the return type of `f`.
"""
function rewrite(f, t :: T, r :: Dict) where T <: Trace
    function g(params...)
        new_t, h = replay(f, t)
        rval = h(params...)
        new_t = replace(new_t, r)
        (new_t, rval)
    end
    g
end

@doc raw"""
    function block(f, t :: T, addresses) where T <: Trace

Given a stochastic function `f`, a trace `t`, and an iterable of addresses, converts traced
randomness into untraced randomness.

Returns a tuple `(t_new, g)`, where `t_new` is a trace and `g` is a function.
The function signature of `g` is the same as that of `f` with the first argument removed; that
is, if `f(t :: Trace, params...)`, then `g(params...)`. Computation is delayed, so each of
the latent nodes in `t_new` has `interpretation = BLOCKED`. Calling `g(params...)` executes the
computation and each latent node in `t_new` with an address in `addresses` is removed.
"""
function block(f, t :: T, addresses) where T <: Trace
    t_new = deepcopy(t)
    for a in addresses
        t_new[a].interpretation = BLOCKED
    end
    g(params...) = f(t_new, params...)
    (t_new, g)
end

function block(f, addresses)
    function g(t::Trace, params...)
        r = f(t, params...)
        for a in addresses
            delete!(t.trace, a)
        end
        logprob!(t)
        r
    end
    g
end

@doc raw"""
    block(f :: F, t :: Trace) where F <: Function

Converts all traced randomness into untraced randomness.
"""
block(f, t :: Trace) = block(f, t, keys(t))

@doc raw"""
    function update(f, r :: SamplingResults{I}) where I <: InferenceType

Given a stochastic function `f` and a `SamplingResults` `r`, update the prior
predictive to the posterior predictive by jointly replacing all latent sample sites with the joint empirical
posterior. Returns a stochastic function `g` with the same call signature as `f`. This function will modify in
place the trace passed into it as the first argument.
"""
function update(f, r :: SamplingResults{I}) where I <: InferenceType
    function g(t :: Trace, params...)
        t_sampled = StatsBase.sample(r.traces)  # maintains correct dims throughout
        kts = keys(t_sampled)
        marginals = Dict()
        for a in kts
            sampled_not_obs = t_sampled[a].observed == false
            not_in_t = !(a in keys(t))
            in_t_and_not_obs = !not_in_t && t[a].observed == false
            if (sampled_not_obs && not_in_t) || in_t_and_not_obs
                sample(t, a, r)
                marginals[a] = r[a]
            end
        end
        h = rewrite(f, t, marginals)
        r = h(params...)
    end
    g
end

@doc raw"""
    function condition(f, evidence :: Dict)

Condition a trace modified by `f` on `evidence`, which maps addresses to observed
evidence associated with that address. Returns a function with call signature
identical to that of `f` and return signature `(t :: Trace, rtype)` where `rtype`
is the return type of `f`.
"""
function condition(f, evidence :: Dict; bare=false)
    function g(t :: Trace, params...)
        for (a, e) in evidence
            t[a] = node(e, a, [CONDITIONED], true, CONDITIONED)
        end
        r = f(t, params...)
        !bare ? (t, r) : r
    end
    g
end

export block, replay, update, replace, rewrite, condition


text:
include("../src/Bf.jl")
using .BfInterpreter

bfArray = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8]
bfArray2 = [0.212046, 0.244306, 0.377881, 0.830958, 0.64464, 0.719155, 0.0548315, 0.721247, 0.21178, 0.533973, 0.708894, 0.657724, 0.33828, 0.233056, 0.389256, 0.341476, 0.108992, 0.569987, 0.773708, 0.501025, 0.676659, 0.199096, 0.509212, 0.840806, 0.367272, 0.751362, 0.19773, 0.379025, 0.349044, 0.808335, 0.426569, 0.657044, 0.199403, 0.704683, 0.316259, 0.833642, 0.729846, 0.180223, 0.428097, 0.353757, 0.45328, 0.565585, 0.955719, 0.437635, 0.111651, 0.61881, 0.677909, 0.114505, 0.226029, 0.266834, 0.272665, 0.472271, 0.905649, 0.0814154, 0.195086, 0.458175, 0.215172, 0.963086, 0.348278, 0.353836, 0.982223, 0.363623, 0.508142, 0.933824, 0.0441713, 0.808518, 0.825324, 0.741406, 0.428842, 0.858324, 0.260967, 0.621218, 0.28133, 0.814533, 0.893156, 0.159597, 0.127285, 0.912122, 0.681966, 0.0892357, 0.570726, 0.211872, 0.304449, 0.489168, 0.478844, 0.150225, 0.838168, 0.151734, 0.227039, 0.221976, 0.465813, 0.400767, 0.352967, 0.709896, 0.18062, 0.847185, 0.341455, 0.789339, 0.0218322, 0.233557, 0.704715, 0.739691, 0.109016, 0.416451, 0.333344, 0.687655, 0.933704, 0.275334, 0.459957, 0.467693, 0.968968, 0.811969, 0.0491648, 0.592472, 0.0321572, 0.21054, 0.676775, 0.584325, 0.563272, 0.271314, 0.449049, 0.062233, 0.0905708, 0.785225, 0.580624, 0.591054, 0.329226, 0.974927, 0.146002, 0.00222847, 0.196426, 0.589131, 0.05884, 0.434156, 0.530672, 0.225588, 0.856623, 0.157457, 0.352767, 0.643674, 0.714834, 0.610502, 0.0992828, 0.163109, 0.66382, 0.752591, 0.567825, 0.893828, 0.964615, 0.211538]

#triche = "++++++++++[>+>+++>+++++++>++++++++++<<<<-]>>>++.>+.+++++++..+++.<<++++++++++++++.------------.>+++++++++++++++.>.+++.------.--------.<<.+."
triche = ",[>++>+<+<-]>.++-++++++++++++++++-+-+[+++..><.[++>.,.-].,[[+,,<++]+>.]+.>+.<,<,,"
#bft2 = bfType(triche)

#bft = bfType(bfArray2)
#@show join(bft.bfcode, "")

#b = runnable!(bft)

#println(" bf code is runnable : $b")
if true
    #try
         #output= execute(bft2, ['t', 'o', 'i', 'm', 'i', 'd','t', 'o', 'i', 'm', 'i', 'd', 't', 'o', 'i', 'm', 'i', 'd'])
         instructionsSet = BfInterpreter.getInstructionsDict()
         output = execute(triche , UInt8[18], instructionsSet)
         #output = execute(bft2, [Char(10), Char(1)])
         #output = execute(bft2)
         #@show Int(output[1][1])
         if length(output) < 40
             @show join(output, "")
         else
             @show join(output[1:40], "")
         end
         println()
     #=catch y
         if isa(y, InexactError)
             println("WARNING : cell is corrupted : negative value : \n $y")
         end
         println("toto $y")
     #end =#
 else
     println(" This code is not runnable")
end


text:
###############################################################################
#
# Containers for Simulation results
#
###############################################################################
"""
    mutable struct SPACContainer1L{FT}

Struct that contains 1-layer gas exchange information.

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct SPACContainer1L{FT<:AbstractFloat}
    "Mean gross photosynthetic rate `[μmol m⁻² s⁻¹]`"
    ag::FT = FT(0)
    "Mean net photosynthetic rate `[μmol m⁻² s⁻¹]`"
    an::FT = FT(0)
    "Leaf internal CO₂ partial pressure `[Pa]`"
    c ::FT = FT(0)
    "Flow rate per basal area `[mol s⁻¹ m⁻²]`"
    e ::FT = FT(0)
    "Leaf diffusive conductance to H₂O `[mol m⁻² s⁻¹]`"
    gh::FT = FT(0)
    "Xylem end pressure `[MPa]`"
    p ::FT = FT(0)
    "Leaf temperature `[K]`"
    t ::FT = FT(0)
end

"""
    mutable struct SPACContainer2L{FT}

Struct that contains 2-layer gas exchange information.

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct SPACContainer2L{FT<:AbstractFloat}
    # gas exchange information
    cont_sh::SPACContainer1L{FT} = SPACContainer1L{FT}();
    cont_sl::SPACContainer1L{FT} = SPACContainer1L{FT}();

    # partition information
    "Shaded layer fraction"
    frac_sh::FT = FT(0)
    "Sunlit layer fraction"
    frac_sl::FT = FT(0)
    "Shaded layer leaf area `[m²]`"
    la_sh  ::FT = FT(0)
    "Sunlit layer leaf area `[m²]`"
    la_sl  ::FT = FT(0)
    "Shaded layer LAI"
    lai_sh ::FT = FT(0)
    "Sunlit layer LAI"
    lai_sl ::FT = FT(0)
    "Shaded layer PAR `[μmol m⁻² s⁻¹]`"
    par_sh ::FT = FT(0)
    "Sunlit layer PAR `[μmol m⁻² s⁻¹]`"
    par_sl ::FT = FT(0)
    "Shaded layer absorbed energy `[W m⁻²]`"
    rad_sh ::FT = FT(0)
    "Sunlit layer absorbed energy `[W m⁻²]`"
    rad_sl ::FT = FT(0)
end


text:
struct FokkerPlanckOperator{M,N,ℳ,𝒩,RT,CT,HT,MT,WT} <: CollisionOperator{M,N,ℳ,𝒩,RT,CT}
    grid::Grid2d{M,N,RT}
    ft::FourierTransform{ℳ,𝒩,RT,CT}
    D::Vector{Matrix{CT}}
    Δ⁻¹::Matrix{CT}

    hfunc::HT
    mfunc::MT
    wfunc::WT

    Du::Vector{SharedArray{RT,2}}
    Dh::Vector{SharedArray{RT,2}}
    m::SharedArray{RT,2}

    Dû::Vector{SharedArray{CT,2}}
    Dĥ::Vector{SharedArray{CT,2}}
    m̂::SharedArray{CT,2}

    𝔽::Vector{SharedArray{RT,2}}
    𝔻::Matrix{SharedArray{RT,2}}

    function FokkerPlanckOperator{M,N,ℳ,𝒩,RT,CT,HT,MT,WT}(grid::Grid2d{M,N,RT}, ft::FourierTransform{ℳ,𝒩,RT,CT},
                                                           hfunc::HT, mfunc::MT, wfunc::WT) where {M,N,ℳ,𝒩,RT,CT,HT,MT,WT}
        D   = get_gradient(ft)
        Δ⁻¹ = get_inverse_laplacian(ft, sign=-1)

        Du = [SharedArray{RT}((M,N)), SharedArray{RT}((M,N))]
        Dh = [SharedArray{RT}((M,N)), SharedArray{RT}((M,N))]
        m  = SharedArray{RT}((M,N))

        Dû = [SharedArray{CT}((ℳ,𝒩)), SharedArray{CT}((ℳ,𝒩))]
        Dĥ = [SharedArray{CT}((ℳ,𝒩)), SharedArray{CT}((ℳ,𝒩))]
        m̂  = SharedArray{CT}((ℳ,𝒩))

        𝔽 = [SharedArray{RT}((M,N)), SharedArray{RT}((M,N))]
        𝔻 = [SharedArray{RT}((M,N)) for k ∈ 1:2, l ∈ 1:2]

        new(grid, ft, D, Δ⁻¹, hfunc, mfunc, wfunc, Du, Dh, m, Dû, Dĥ, m̂, 𝔽, 𝔻)
    end
end

function FokkerPlanckOperator(grid::Grid2d{M,N,RT}, ft::FourierTransform{ℳ,𝒩,RT,CT};
                              hfunc=hfunc_fokker_planck!, mfunc=mfunc_fokker_planck!, wfunc=wfunc_fokker_planck!) where {M,N,ℳ,𝒩,RT,CT}
    op = FokkerPlanckOperator{M,N,ℳ,𝒩,RT,CT, typeof(hfunc), typeof(mfunc), typeof(wfunc)}(grid, ft, hfunc, mfunc, wfunc)
    nworkers() > 1 ? initialize_workers(op, ft.ℳcut, ft.𝒩cut) : nothing
    return op
end

function initialize_workers(op::FokkerPlanckOperator{M,N,ℳ,𝒩,RT,CT}, ℳcut=0, 𝒩cut=0) where {M,N,ℳ,𝒩,RT,CT}
    # # add parallel workers
    # if haskey(ENV, "JULIA_NUM_THREADS")
    #     nw = parse(Int, ENV["JULIA_NUM_THREADS"])
    # else
    #     nw = 1
    # end
    #
    # addprocs(nw)

    local gDh = op.Dh
    local gDu = op.Du
    local gm  = op.m

    local gDû = op.Dû
    local gm̂  = op.m̂

    local g𝔽  = op.𝔽
    local g𝔻  = op.𝔻

    local gwfunc  = op.wfunc

    local wn::Int = div(N, nworkers())

    # @eval @everywhere Dh = $gDh
    # @eval @everywhere Du = $gDu
    # @eval @everywhere m  = $gm

    @sync for (p,w) in enumerate(workers())
        @spawnat w begin

            global Dh = gDh
            global Du = gDu
            global m  = gm

            global Dû = gDû
            global m̂  = gm̂

            global 𝔽  = g𝔽
            global 𝔻  = g𝔻

            global wfunc = gwfunc

            global j1 = wn*(p-1) + 1
            global j2 = wn*p

            global gr = Grid2d(M,N)
            global ft = FourierTransform(gr, ℳcut=ℳcut, 𝒩cut=𝒩cut)

            # global lg = [zeros(RT,M,N), zeros(RT,M,N)]
            # global lw = Array{Array{RT,2}}(2,2)
            # global lŵ = Array{Array{CT,2}}(2,2)
            #
            # for l in 1:size(lw,2)
            #     for k in 1:size(lw,1)
            #         lw[k,l] = zeros(RT,M,N)
            #     end
            # end
            #
            # for l in 1:size(lŵ,2)
            #     for k in 1:size(lŵ,1)
            #         lŵ[k,l] = zeros(CT,ℳ,𝒩)
            #     end
            # end

        end
    end
end

@generated function collision_operator!(op::FokkerPlanckOperator{M,N,ℳ,𝒩,RT,CT,HT,MT,WT}, u::Matrix{RT}, divJ::Matrix{RT}) where {M,N,ℳ,𝒩,RT,CT,HT,MT,WT}

    local ϕ::Matrix{RT}
    local h::Matrix{RT}
    local m::Matrix{RT}

    local Du::Matrix{RT}
    local Dh::Matrix{RT}

    local û::Matrix{CT}
    local ϕ̂::Matrix{CT}
    local ĥ::Matrix{CT}
    local m̂::Matrix{CT}

    local Dû::Matrix{CT}
    local Dĥ::Matrix{CT}

    local J::Vector{Matrix{RT}}
    local Ĵ::Vector{Matrix{CT}}

    local divĴ::Matrix{CT}

    ϕ  = zeros(RT,M,N)
    h  = zeros(RT,M,N)
    m  = zeros(RT,M,N)

    Du = zeros(RT,M,N)
    Dh = zeros(RT,M,N)

    û  = zeros(CT,ℳ,𝒩)
    ϕ̂  = zeros(CT,ℳ,𝒩)
    ĥ  = zeros(CT,ℳ,𝒩)
    m̂  = zeros(CT,ℳ,𝒩)

    Dû = zeros(CT,ℳ,𝒩)
    Dĥ = zeros(CT,ℳ,𝒩)

    J  = [zeros(RT,M,N), zeros(RT,M,N)]
    Ĵ  = [zeros(CT,ℳ,𝒩), zeros(CT,ℳ,𝒩)]

    divĴ = zeros(CT,ℳ,𝒩)

    quote
        @assert size(u) == size(divJ) == (M,N)

        frfft!(op.ft, u, $û)

        apply_operator!(op.Δ⁻¹, $û, $ϕ̂)
        irfft!(op.ft, $ϕ̂, $ϕ)

        op.mfunc(u, $m, op.grid)
        frfft!(op.ft, $m, $m̂)

        op.m .= $m
        op.m̂ .= $m̂

        op.hfunc(u, $ϕ, $h, op.grid)
        frfft!(op.ft, $h, $ĥ)

        apply_operator!(op.D, $û, op.Dû)
        apply_operator!(op.D, $ĥ, op.Dĥ)

        # FFT does not seem to work on Shared Array!
        #
        # for k in 1:length(op.Du)
        #     irfft!(op.ft, op.Dû[k], op.Du[k])
        # end
        #
        # for k in 1:length(op.Dh)
        #     irfft!(op.ft, op.Dĥ[k], op.Dh[k])
        # end

        for k in 1:length(op.Du)
            $Dû .= op.Dû[k]
            irfft!(op.ft, $Dû, $Du)
            op.Du[k] .= $Du
        end

        for k in 1:length(op.Dh)
            $Dĥ .= op.Dĥ[k]
            irfft!(op.ft, $Dĥ, $Dh)
            op.Dh[k] .= $Dh
        end

        compute_convolution!(op)

        for k in 1:length($J)
            @inbounds for j in 1:size($J[k],2)
                for i in 1:size($J[k],1)
                    $J[k][i,j] = op.𝔽[k][i,j] * op.m[i,j]
                end
            end

            for l in 1:length($J)
                @inbounds for j in 1:size($J[k],2)
                    for i in 1:size($J[k],1)
                        $J[k][i,j] += op.𝔻[k,l][i,j] * op.Du[l][i,j]
                    end
                end
            end

            prfft!(op.ft, $J[k], $Ĵ[k])
        end

        apply_operator!(op.D, $Ĵ, $divĴ)

        irfft!(op.ft, $divĴ, divJ)

        divJ
    end
end

function compute_convolution!(op::FokkerPlanckOperator{M,N,ℳ,𝒩,RT,CT}) where {M,N,ℳ,𝒩,RT,CT}
    if nworkers() > 1
        # @sync @parallel for j = 1:N
        #     # convolution_kernel_fourier!(gr, ft, j, j, wfunc, Dh, Dû, m̂, 𝔽, 𝔻)
        #     convolution_kernel_trapezoidal!(gr, j, j, wfunc, Dh, Du, m, 𝔽, 𝔻)
        # end

        @sync for w in workers()
            @spawnat w begin
                # convolution_kernel_fourier!(gr, ft, j1, j2, wfunc, Dh, Dû, m̂, 𝔽, 𝔻)
                convolution_kernel_trapezoidal!(gr, j1, j2, wfunc, Dh, Du, m, 𝔽, 𝔻)
            end
        end

        # @sync for w in workers()
        #     # @async remotecall_wait(convolution_kernel_fourier!, w, gr, ft, j1, j2, wfunc, Dh, Dû, m̂, 𝔽, 𝔻)
        #     @async remotecall_wait(convolution_kernel_trapezoidal!, w, gr, j1, j2, wfunc, Dh, Du, m, 𝔽, 𝔻)
        # end

    else
        # convolution_kernel_fourier!(op.grid, op.ft, 1, N, op.wfunc, op.Dh, op.Dû, op.m̂, op.𝔽, op.𝔻)
        convolution_kernel_trapezoidal!(op.grid, 1, N, op.wfunc, op.Dh, op.Du, op.m, op.𝔽, op.𝔻)
    end
end

@generated function convolution_kernel_fourier!(gr::Grid2d{M,N}, ft::FourierTransform{ℳ,𝒩,RT,CT}, j1::Int, j2::Int, wfunc::Function,
                                                Dh::Vector{SharedArray{RT,2}}, Dû::Vector{SharedArray{CT,2}}, m̂::SharedArray{CT,2},
                                                𝔽::Vector{SharedArray{RT,2}}, 𝔻::Matrix{SharedArray{RT,2}}) where {M,N,ℳ,𝒩,RT,CT}

    local g::Vector{Matrix{RT}} = [zeros(RT,M,N), zeros(RT,M,N)]
    local w::Matrix{Matrix{RT}} = Array{Array{RT,2}}(2,2)
    local ŵ::Matrix{Matrix{CT}} = Array{Array{CT,2}}(2,2)

    for l in 1:size(w,2)
        for k in 1:size(w,1)
            w[k,l] = zeros(RT,M,N)
        end
    end

    for l in 1:size(ŵ,2)
        for k in 1:size(ŵ,1)
            ŵ[k,l] = zeros(CT,ℳ,𝒩)
        end
    end

    quote
        @assert 1 ≤ j1 ≤ j2 ≤ N
        @assert length($g) == length(Dh) == length(Dû)
        @assert length($g) == length(𝔽) == size(𝔻,1) == size(𝔻,2)

        for j in j1:j2
            for i in 1:M

                for k in 1:length($g)
                    $g[k] .= Dh[k][i,j] .- Dh[k]
                end

                wfunc(gr, i, j, $g, $w)

                for l in 1:size($w,2)
                    for k in 1:size($w,1)
                        prfft!(ft, $w[k,l], $ŵ[k,l])
                    end
                end

                for k in 1:length(𝔽)
                    𝔽[k][i,j] = 0
                    for l in 1:length(𝔽)
                        𝔽[k][i,j] -= real(fourier_quadrature($ŵ[k,l], Dû[l], ft.μ, gr))
                    end
                end

                for l in 1:size(𝔻,2)
                    for k in 1:size(𝔻,1)
                        𝔻[k,l][i,j] = real(fourier_quadrature($ŵ[k,l], m̂, ft.μ, gr))
                    end
                end

            end
        end
    end
end

@generated function convolution_kernel_trapezoidal!(gr::Grid2d{M,N}, j1::Int, j2::Int, wfunc::Function,
                                                    Dh::Vector{SharedArray{RT,2}}, Du::Vector{SharedArray{RT,2}}, m::SharedArray{RT,2},
                                                    𝔽::Vector{SharedArray{RT,2}}, 𝔻::Matrix{SharedArray{RT,2}}) where {M,N,RT}

    local g::Vector{Matrix{RT}} = [zeros(RT,M,N), zeros(RT,M,N)]
    local w::Matrix{Matrix{RT}} = [zeros(RT,M,N) for k ∈ 1:2, l ∈ 1:2]

    quote
        @assert 1 ≤ j1 ≤ j2 ≤ N
        @assert length($g) == length(Dh) == length(Du)
        @assert length($g) == length(𝔽) == size(𝔻,1) == size(𝔻,2)

        for j in j1:j2
            for i in 1:M

                for k in 1:length($g)
                    $g[k] .= Dh[k][i,j] .- Dh[k]
                end

                wfunc(gr, i, j, $g, $w)

                for k in 1:length(𝔽)
                    𝔽[k][i,j] = 0
                    for l in 1:length(𝔽)
                        𝔽[k][i,j] -= trapezoidal_quadrature($w[k,l], Du[l], gr)
                    end
                end

                for l in 1:size(𝔻,2)
                    for k in 1:size(𝔻,1)
                        𝔻[k,l][i,j] = trapezoidal_quadrature($w[k,l], m, gr)
                    end
                end

            end
        end
    end
end

function mfunc_fokker_planck!(u::Matrix{RT}, m::Union{Array{RT,2},SharedArray{RT,2}}, grid::Grid2d) where {RT}
    m .= 1
end

function hfunc_fokker_planck!(u::Matrix{RT}, ϕ::Matrix{RT}, h::Union{Array{RT,2},SharedArray{RT,2}}, grid::Grid2d) where {RT}
    h .= ϕ
end

@generated function wfunc_fokker_planck!(grid::Grid2d{M,N,RT}, i::Int, j::Int, g::Vector{Matrix{RT}}, w::Matrix{Matrix{RT}}) where {M,N,RT}
    # local g¹::Matrix{RT} = zeros(RT,M,N)
    local g²::Matrix{RT} = zeros(RT,M,N)

    quote
        @assert size(w,1) == size(w,2) == length(g)
        @assert size(w[1,1],1) == size(g[1],1) == M
        @assert size(w[1,1],2) == size(g[1],2) == N

        for l in 1:size(w,2)
            for k in 1:size(w,1)
                w[k,l] .= 0
            end
        end

        $g² .= 0

        for k in 1:size(g,1)
            @inbounds for j in 1:size($g²,2)
                for i in 1:size($g²,1)
                    $g²[i,j] += g[k][i,j]^2
                end
            end
        end

        # $g¹ .= sqrt.($g²)

        for k in 1:size(w,1)
            @inbounds for j in 1:size(w[k,k],2)
                for i in 1:size(w[k,k],1)
                    w[k,k][i,j] = $g²[i,j]
                end
            end
        end

        for l in 1:size(w,2)
            for k in 1:size(w,1)
                @inbounds for j in 1:size(w[k,l],2)
                    for i in 1:size(w[k,l],1)
                        w[k,l][i,j] -= g[k][i,j] * g[l][i,j]
                    end
                end
            end
        end

        w
    end
end


text:
# solve the minimum length Hamiltonian problem; in this algorithm, the
# initial and final stops are *fixed* and the feasible region consists of alternating
# intermediate stops
struct MIP_IF <: HamiltonianSolver end

min_path(T, stops, ::MIP_IF) = _min_length_path_mip(T, stops)

function _min_length_path_mip(T::AbstractMatrix, stops::Vector{<:AbstractStop}; optimizer)
     stop_names = [name(s) for s in stops]
     _min_length_path_mip(T, stop_names, optimizer=optimizer)
end

# solve the minimum length Hamiltonian problem using a MIP formulation
# *where the first and last stops are fixed*
function _min_length_path_mip(T::AbstractMatrix, stops::Vector{<:String}; optimizer)
    M = view(T, stops, stops)

    model = Model(optimizer)

    # number of stops
    N = length(stops)

    # binary decision variables: 1 if x_{st} stop s goes to stop t and 0 otherwise
    @variable(model, x[1:N, 1:N], Bin)

    # cost function
    @objective(model, Min, sum(x[s, t] * M[s, t] for s=1:N, t=1:N))

    for s = 1:N  # loop over source stop
        # no edge reaches the initial stop
        @constraint(model, x[s, 1] == 0)

        # no self-loops
        @constraint(model, x[s, s] == 0)

        # each stop arrives to one and only one stop (unless the last one)
        if s < N
            @constraint(model, sum(x[s, t] for t in 2:N) == 1)
        end
    end

    # can't jump directly from first to last stop
    @constraint(model, x[1, N] == 0)

    for t = 1:N  # loop over target stop
        # no edge after the final stop
        @constraint(model, x[N, t] == 0)

        # each stop receives one and only one edge (unless the first stop)
        if t > 1
           @constraint(model, sum(x[s, t] for s in 1:N) == 1)
        end
    end

    # no repetitions (not required in ppple)
    for  s= 1:N, t = 1:N
        @constraint(model, x[s, t] + x[t, s] <= 1)
    end

    optimize!(model)

    while !is_mip_solved(model, x)
        optimize!(model)
    end

    minvalue = JuMP.objective_value(model)
    Mpath = JuMP.getvalue.(x)

    result = similar(stops)
    result[1] = stops[1]
    row = 1
    for i in 1:N-1
        idx = findfirst(==(1,), view(Mpath, row, :))
        @assert !isnothing(idx)
        result[i+1] = stops[idx]
        row = idx
    end
    return minvalue, result
end

# Check whether the Hamiltonian problem is solved,
# which requires that there is a unique chain of length N
# if there is no such chain => we add it as a constraint to the problem
#
# For instance, if there are 5 stops we may find a chain of length 3:
#
# 3 -> 4 -> 3
# 1 -> 2 -> 5 : x[1, 2] + x[2, 5] = 2
# The idea is to add a new constraint that excludes the chain 1 -> 2 -> 5
function is_mip_solved(model, x)
    minvalue = JuMP.objective_value(model)
    Mpath = JuMP.getvalue.(x)

    # check if the problem is solved <=> the path starting at 1 covers all stops
    N = size(Mpath, 1)
    chain_idx = Vector{Int}()
    push!(chain_idx, 1)
    while true
        idx = findfirst(==(1,), view(Mpath, chain_idx[end], :))
        push!(chain_idx, idx)
        idx == N && break
    end

    if length(chain_idx) < N
        @constraint(model, sum(x[chain_idx, chain_idx]) <= length(chain_idx) - 2)
        return false
    end
    return true
end

# solve the minimum length Hamiltonian problem; in this algorithm, the
# initial stop is *fixed*, and optimization is based on the remaining m-1 stops
struct MIP_I <: HamiltonianSolver end

# TODO add MIP_I naive implementation is to iteratively call MIP_IF for each
# of the remaining n - 1 stops


text:
mutable struct State{T, D, L, H, P}
       Qu :: T
       Qv :: T
       Qθ :: T
       Qs :: T
       Qb :: T
        K :: D
        ℓ :: L
        h :: H # boundary layer depth
    plume :: P # plumes
end

function State(grid, mixing_length, boundary_layer_depth, nonlocal_flux; T=Float64)
    mixing_length = instantiate_mixing_length(mixing_length)
    boundary_layer_depth = 4.0 #instantiate_boundary_layer_depth(boundary_layer_depth)
    plume = instantiate_plume(nonlocal_flux, grid)
    K = CellField(grid)
    return State((zero(T) for i=1:5)..., K, mixing_length, boundary_layer_depth, plume)
end

"""
    update_state!(model)

Update the top flux conditions and mixing depth for `model`
and store in `model.state`.
"""
function update_state!(m)
    m.state.Qu = getbc(m, m.bcs.U.top)
    m.state.Qv = getbc(m, m.bcs.V.top)
    m.state.Qθ = getbc(m, m.bcs.T.top)
    m.state.Qs = getbc(m, m.bcs.S.top)
    m.state.Qb = m.constants.g * (m.constants.α * m.state.Qθ - m.constants.β * m.state.Qs)

    zero_out_negative_tke!(m.solution.e)
    update_near_wall_tke!(m)
    update_boundary_layer_depth!(m)

    update_mixing_length!(m)
    update_diffusivity!(m)
    update_nonlocal_flux!(m)

    return nothing
end

# Fallbacks
instantiate_plume(nonlocal_flux, grid) = nothing #(T=nothing, S=nothing, e=nothing, W²=nothing)
update_nonlocal_flux!(model) = nothing
update_near_wall_tke!(m) = nothing
instantiate_mixing_length(mixing_length_parameters) = nothing
instantiate_boundary_layer_depth(boundary_layer_depth_parameters) = nothing
update_mixing_length!(m) = nothing
update_boundary_layer_depth!(m) = nothing

@inline function zero_out_negative_tke!(e)
    for i in eachindex(e)
        @inbounds e[i] = ifelse(e[i] < 0, zero(eltype(e)), e[i])
    end
    return nothing
end

function update_diffusivity!(m)
    for i in 1:m.grid.N
        @inbounds m.state.K[i] = diffusivity_mixing_length(m, i) * sqrt_e(m, i)
    end

    # Neumann condition on eddy diffusivity over boundary
    @inbounds m.state.K[0] = m.state.K[1]
    @inbounds m.state.K[m.grid.N+1] = m.state.K[m.grid.N]

    return nothing
end


text:
# --- PRIVATE METHODS BELOW HERE ----------------------------------------------------------------- #
# --- PRIVATE METHODS ABOVE HERE ----------------------------------------------------------------- #

# --- PUBLIC METHODS BELOW HERE ------------------------------------------------------------------ #
function lsm_price(contractSet::Set{VLAbstractAsset}, underlyingPriceTable::Array{Float64,2};
    riskFreeRate::Float64 = 0.015, decisionLogic::Function=_american_decision_logic)::VLResult

    # main -
    try

        # what is the size of the price table?
        (number_of_rows, number_of_cols) = size(underlyingPriceTable)

    catch error
        return VLResult(error)
    end
end
# --- PUBLIC METHODS ABOVE HERE ------------------------------------------------------------------ #


text:
using LinearAlgebra
using Plots
using ForwardDiff: HessianConfig,JacobianConfig,hessian!,jacobian!,hessian,jacobian

"""
Some struct to make things tidy
"""
struct Target
    neglog_density::Function#the neglog_density function
    d::Integer#Dimension of the problem
    hessian_cfg::Union{HessianConfig}
    hessian_result::Array{<:Real,2}
    u::AbstractArray{<:Real,1}#normalization parameter for modifiedCholesky
end

"""
The Soft absolute value function, with scalar inputs
"""
function sabs(x::Real,u::Real)
    log2 = log(2)
    temp = (x*log2/u)
    ret = log(exp(temp)+exp(-temp))*u/log2
    return ret
end

"""
Implementation of Modified Cholesky Decomposition routine
taken from :
`Modified Cholesky Riemann Manifold Hamiltonian Monte Carlo:
exploiting sparsity for fast sampling of high-dimensional targets`
    by Tore Selland Kleppe. DOI 10.1007/s11222-017-9763-5
"""
function modifiedCholesky(A::Symmetric{Float64,Array{Float64,2}},u::AbstractArray{<:Real,1},K::Integer=1)
    size_A = size(A)
    d = size_A[1]
    Id = Matrix{Float64}(I,size_A)
    L̃ = Id
    D = diag(A)
    for j in 1:size_A[1]
        if j>1
            L̃[j,1:j-1] = L̃[j,1:j-1]./D[1:j-1]

            if j<d
                L̃[j+1:d,j] = A[j+1:d,j]
                L̃[j+1:d,j] = L̃[j+1:d,j] - L̃[j+1:d,1:j-1]*L̃[j,1:j-1]
            end
        else
            if j<d
                L̃[j+1:d,j] = A[j+1:d,j]
            end
        end
        if j>K
            D[j] = sabs(D[j],u[j])
        end
        if j<d
            D[j+1:d] = D[j+1:d] - L̃[j+1:d,j].^2/D[j]
        end
    end
    L = LowerTriangular(L̃*diagm(D))
    return L
end

"""
the metric matrix function, evaluated for a specific negative log target probability distribution
"""
function metric(x::AbstractArray{<:Real,1},target::Target)
    H = Symmetric(hessian!(target.hessian_result,target.neglog_density,x,target.hessian_cfg))
    L = modifiedCholesky(H,target.u)
    return L
end

"""
Hamiltonian of a Reimannian Manifold HMC, see eq 1 of
`Modified Cholesky Riemann Manifold Hamiltonian Monte Carlo:
exploiting sparsity for fast sampling of high-dimensional targets`
    by Tore Selland Kleppe. DOI 10.1007/s11222-017-9763-5
"""
# function hamiltonian(x::AbstractArray{<:Real,1},p::AbstractArray{<:Real,1},u::AbstractArray{<:Real,1},
#     neglog_target_pi,hessian_tape::HessianTape,hessian_result::Array{<:Real,2})
function hamiltonian(x::AbstractArray{<:Real,1},p::AbstractArray{<:Real,1},target::Target)

    L = metric(x,target)
    # L = L̃*diagm(D)
    p̂ = L\p
    ham = target.neglog_density(x) + 0.5p̂'*p̂ + logdet(L)
    return ham
end

struct Hamiltonian
    target::Target
    fun::Function
    jacobian!::Function
    jacobian_tape::Union{JacobianTape,CompiledTape}
    jacobian_result::Array{<:Real,1}
    Hamiltonian(target) = begin
        fun(x,p) = hamiltonian(x,p,target)
        x,p = randn(target.d),randn(target.d)
        jacobian_tape = JacobianTape(fun,(x,p))
        jacobian_result = similar(x)
    end
end

"""
Explicit Leap Frog algorithm: using 10.1103/PhysRevE.94.043303
"""
function generalized_leap_frog(ϵ::Real,ℓ::Integer,ω::Real
                            ,x::AbstractArray{<:Real,1},p::AbstractArray{<:Real,1},
                            target::Target)
    c = cos(2ω*ϵ)
    s = sin(2ω*ϵ)

    x̃ = copy(x)
    p̃ = copy(p)

    for i in 1:ℓ
        #2
        jacobian!(target.jacobian_result,jacobian_tape,(x,p̃))
        p -= 0.5ϵ*(target.jacobian_result)

        #3
        L = metric(x,u,target.hessian_tape,target.hessian_result)
        x̃ += 0.5ϵ*(L'\L\p̃) #this is more efficient

        #4
        jacobian!(target.jacobian_result,target.jacobian_tape,(x̃,p))
        p̃ -= 0.5ϵ*(target.jacobian_result)

        #5
        L = metric(x̃,u,target.hessian_tape,target.hessian_result)
        x += 0.5ϵ*(L'\L\p) #this is more efficient

        #6 no need to compute c = cos(2ωϵ); s = sin(2ωϵ) every timestep
        #7
        x = (x+x̃ + c*(x-x̃) + s*(p-p̃))/2

        #8
        p = (p+p̃ - s*(x-x̃) + c*(p-p̃))/2

        #9
        x̃ = (x+x̃ - c*(x-x̃) - s*(p-p̃))/2

        #10
        p̃ = (p+p̃ + s*(x-x̃) - c*(p-p̃))/2

        #11
        jacobian!(target.jacobian_result,target.jacobian_tape,(x̃,p))
        p̃ -= 0.5ϵ*(target.jacobian_result)

        #12
        L = metric(x̃,u,target.hessian_tape,target.hessian_result)
        x += 0.5ϵ*(L'\L\p) #this is more efficient

        #13
        jacobian!(target.jacobian_result,target.jacobian_tape,(x,p̃))
        p -= 0.5ϵ*(target.jacobian_result)

        #14
        L = metric(x,u,target.hessian_tape,target.hessian_result)
        x̃ += 0.5ϵ*(L'\L\p̃) #this is more efficient
    end

    return x,p
end

"""
Example of target density
"""
function funnel_neglog(x::AbstractArray{<:Real,1})
    nLP = (x[1]^2/(2exp(x[2])) + x[2]/2 + x[2]^2/18)
    # nLP = (x^2/(2exp(y)) + y/2 + y^2/18)
    return nLP
end

n = 2
x = randn(n)
p = randn(n)
w = randn(n)
# hessConfig = HessianConfig(x)
hessian_result = similar(x,n,n)

ϵ = 0.5n^(-0.25)
ℓ = 1.5/ϵ

target_pi = funnel_neglog
hess_cfg = HessianConfig(target_pi,x)

hessian!(hessian_result,target_pi,x,hess_cfg)
u = eigmin(hessian_result)ones(n)
target = Target(funnel_neglog,n,hess_cfg,hessian_result,abs.(u))
metric(x,target)
hamiltonian(x,p,target)

target_hamiltonian(x::AbstractArray{<:Real,1}) = hamiltonian(x::AbstractArray{<:Real,1},p::AbstractArray{<:Real,1},target)
jac_cfg = JacobianConfig(target_hamiltonian,x)
jacobian_result = similar(x)
jacobian!(jacobian_result,target_hamiltonian,x)

# #test function
# function negLogPosterior(x)
#     return sqrt(sum(x.^2))
# end

# """
# Julia does not have a meshgrid function?
# """
# function meshgrid(x, y)
#     return (repeat(x',(length(y)),1), repeat(y,1,(length(x))))
# end

#

#
# u = [1:n...]/1.
#
# L̃,D = modifiedCholesky(Symmetric(H),u,1)
#
# x = [-5:0.1:30...]
# y=[-5:0.1:5...]
# X,Y = meshgrid(x,y)
# Z = @. exp(-(X^2/(2exp(Y)) + Y/2 + Y^2/18))
# heatmap(x,y,Z,fill=(true,:greys_r))

max_stars_repo_path: Julia files (legacy)/RMHMC_f.jl

max_issues_repo_path: Julia files (legacy)/RMHMC_f.jl

max_forks_repo_path: Julia files (legacy)/RMHMC_f.jl