text:
# detect if using SLURM
const IN_SLURM = "SLURM_JOBID" in keys(ENV)

using Distributed
IN_SLURM && using ClusterManagers

using ShaleDrillingModel
using Test
using StatsFuns
using Interpolations
using SparseArrays
using Calculus
using MarkovTransitionMatrices

# -------------------------------------------------
# price transition matrices
# -------------------------------------------------

nlogp = 15
nlogc = 11
nlogσ = 9

extrema_logp = [0.8510954,  2.405270]
extrema_logc = [5.0805023,  5.318393]
extrema_logσ = [-2.9558226, -1.542508]
sdlogσ = 0.18667531

Σpc = [0x1.4a1ecb9bcddfep-7  0x1.94e52859805cap-10;
       0x1.94e52859805cap-10 0x1.4cf39bf87b735p-9]

# make range for coefs
logp_space = range(extrema_logp[1] - log(2.0), stop=extrema_logp[2] + log(2.0), length=nlogp)
logc_space = range(extrema_logc[1] - log(2.0), stop=extrema_logc[2] + log(2.0), length=nlogc)
logσ_space = range(extrema_logσ[1] - 2*sdlogσ, stop=extrema_logσ[2] + 2*sdlogσ, length=nlogσ)

P_σ() = tauchen_1d(logσ_space, (x) -> x, sdlogσ)
P_pricecost() = tauchen_2d(Base.Iterators.product(logp_space, logc_space), (x) -> x, Σpc)
P_price() = tauchen_1d(logp_space, (x) -> x, Σpc[1,1])

# range for i,j block in matrix where block is n x n
blockrange(i,j,n) = ((i-1)*n+1):(i*n),  ((j-1)*n+1):(j*n)

function P_pricevol()
    Pσ = P_σ()
    P  = Matrix{eltype(logp_space)}(undef, nlogp*nlogσ, nlogp*nlogσ)
    for (j, logσj) in enumerate(logσ_space)
      for (i, logσi) in enumerate(logσ_space)
        P[blockrange(i,j,nlogp)...]        .= Pσ[i, j] .* tauchen_1d(logp_space, (x) -> x, exp(logσi)^2)
      end
    end
    return P
end
function P_pricecostvol()
    Pσ = P_σ()
    P = Matrix{eltype(logp_space)}(undef, nlogp*nlogc*nlogσ, nlogp*nlogc*nlogσ)
    for (j, logσj) in enumerate(logσ_space)
      for (i, logσi) in enumerate(logσ_space)
        P[blockrange(i,j,nlogp*nlogc)...] .= Pσ[i, j] .* tauchen_2d(Base.Iterators.product(logp_space, logc_space), (x) -> x, Σpc*exp(logσi))
      end
    end
    return P
end

minp = 1e-4
Πp      = MarkovTransitionMatrices.sparsify!(P_pricevol(),     minp)
Πpc     = MarkovTransitionMatrices.sparsify!(P_pricecostvol(), minp)
Πpconly = MarkovTransitionMatrices.sparsify!(P_pricecost(),    minp)

# -------------------------------------------------
# price transition matrices
# -------------------------------------------------

# some primitives
β = (1.02 / 1.125) ^ (1.0/4.0)  # real discount rate
royalty_rates = [0.125, 1.0/6.0, 0.1875, 0.20, 0.225, 0.25]
royalty_types = 1:length(royalty_rates)
geology_types = 1.3430409262656042:0.1925954901417719:5.194950729101042

# initial parameters
flowfuncname = :one_restr
θt = [-4.28566, -5.45746, -0.3, ] # ShaleDrillingModel.STARTING_log_ogip, ShaleDrillingModel.STARTING_σ_ψ,
σv = 1.0

θfull = vcat(θt, σv)

# problem sizes
nψ, dmx, nz, nv =  51, 3, size(Πp,1), 51
# wp = LeasedProblemContsDrill(dmx,4,5,3,2)
wp = LeasedProblem(dmx,4,5,3,2)
# wp = PerpetualProblem(dmx,4,5,3,2)
# lp = LeasedProblem(dmx,4,1,-1,0)

zspace, ψspace, dspace, d1space, vspace = (logp_space, logσ_space,), range(-4.5, stop=4.5, length=nψ), 0:dmx, 0:1, range(-3.0, stop=3.0, length=nv)

prim = dcdp_primitives(flowfuncname, β, wp, zspace, Πp, ψspace)
print_summary(prim)
tmpv = dcdp_tmpvars(prim)
evs = dcdp_Emax(prim)

## check sizes of models
ShaleDrillingModel.check_size(prim, evs)

include("BSplineTestFuns_runtests.jl")
include("state-space.jl")
# include("flow-payoffs.jl")

@testset  "testing flow gradients" begin
    let geoid = 2, roy = 0.2
        @test check_flowgrad(θt, σv, prim, geoid, roy)
    end
    @test check_dΠψ(σv, ψspace)
end

include("test_utility.jl")
include("test_transition.jl")

include("logsumexp3.jl")
include("vfit.jl")
include("vf_solve_terminal_and_infill.jl")
include("vf_solve_exploratory.jl")

zero!(tmpv)
let geoid = 2, roy = 0.25, itype = (geoid, roy,)
    solve_vf_all!(evs, tmpv, prim, θt, σv, itype, true)
end

include("vf_interpolation.jl")
include("parallel_solution.jl")
include("action_probabilities_new.jl")


text:
Metric(mat::T) where {N, T<:SMatrix{N, N, Float64}} = Metric{N, T}(Symmetric(mat))

function Metric(mat::AbstractMatrix{<:Real})
    N, M = size(mat)
    @assert N == M
    return Metric(SMatrix{N, N, Float64}(mat))
end

"""
    Metric(N::Int)

The Euclidean metric in `N` dimensions.
"""
Metric(N::Int) = Metric(SMatrix{N,N}(1.0I))

function Wedge(vectors::Vector{SVector{N, Float64}}) where N
    K = length(vectors)
    Wedge{N, K}(SVector{K, SVector{N, Float64}}(vectors))
end

"""
    Wedge(s::Simplex{N, K}) where {N, K}

The wedge product of all vectors emanating from the first vertex of a simplex.
"""
function Wedge(s::Simplex{N, K}) where {N, K}
    vectors = SVector{N, Float64}[p.coords - s.points[1].coords for p in s.points[2:end]]
    return Wedge(vectors)
end

export inner_product
"""
    inner_product(m::Metric{N}, v1::SVector{N, Float64}, v2::SVector{N, Float64}) where N
    inner_product(m::Metric{N}, w1::Wedge{N, K}, w2::Wedge{N, K}) where {N, K}

Compute the inner product of two vectors with respect to the given metric.
"""
inner_product(m::Metric{N}, v1::SVector{N, Float64}, v2::SVector{N, Float64}) where N =
    transpose(v1) * m.mat * v2

function inner_product(m::Metric{N}, w1::Wedge{N, K}, w2::Wedge{N, K}) where {N, K}
    if K == 0
        return 1.0
    else
        inner_products = hcat([[inner_product(m, v1, v2) for v1 in w1.vectors]
            for v2 in w2.vectors]...)
        return det(inner_products)
    end
end

export norm_square
"""
    norm_square(m::Metric{N}, v::SVector{N, Float64}) where N
    norm_square(m::Metric{N}, w::Wedge{N}) where N

Compute the inner product of a vector with itself with respect to the given metric. Note
that this value can be negative if the metric is not positive semi-definite.
"""
norm_square(m::Metric{N}, v::SVector{N, Float64}) where N = inner_product(m, v, v)

norm_square(m::Metric{N}, w::Wedge{N}) where N = inner_product(m, w, w)

import LinearAlgebra: norm
"""
    norm(m::Metric{N}, v::SVector{N, Float64}) where N
    norm(m::Metric{N}, w::Wedge{N}) where N

Compute the square root of the magnitude of the inner product of a vector with itself with
respect to the metric `m`.
"""
norm(m::Metric{N}, v::SVector{N, Float64}) where N = sqrt(abs(norm_square(m, v)))

norm(m::Metric{N}, w::Wedge{N}) where N = sqrt(abs(norm_square(m, w)))

export volume_square
"""
    volume_square(m::Metric{N}, s::Simplex{N, K}) where {N, K}

Compute the squared volume of a simplex with respect to the given metric. Note that this
value can be negative if the metric is not positive semi-definite.
"""
volume_square(m::Metric{N}, s::Simplex{N, K}) where {N, K} =
    norm_square(m, Wedge(s))/factorial(K-1)^2

export volume
"""
    volume(m::Metric{N}, s::Simplex{N, K}) where {N, K}

Compute the volume of a simplex with respect to the given metric.
"""
volume(m::Metric{N}, s::Simplex{N, K}) where {N, K} =
    sqrt(abs(norm_square(m, Wedge(s))))/factorial(K-1)

ext: jl


text:
using Test
using FillArrays
using HybridSystems
using Polyhedra

include("automata.jl")

@testset "Switched System" begin
    square = convexhull([1, 1], [-1, 1], [-1, -1], [1, -1])
    simplex = convexhull([0, 0], [0, 1], [1, 0])

    @testset "Switched System" begin
        s = discreteswitchedsystem([ones(2, 2), zeros(2, 2)])
        @test nstates(s) == 1
        @test states(s) == Base.OneTo(1)
        @test statedim(s, 1) == 2
        @test ntransitions(s) == 2
        @test eltype(transitions(s)) == transitiontype(s)
        @test transitions(s) == HybridSystems.OneStateTransition.([1, 2])
        @test source.(s, transitions(s)) == [1, 1]
        @test event.(s, transitions(s)) == [1, 2]
        @test target.(s, transitions(s)) == [1, 1]
        @test length(transitions(s)) == 2
        @test in_transitions(s, 1) == HybridSystems.OneStateTransition.([1, 2])
        @test out_transitions(s, 1) == HybridSystems.OneStateTransition.([1, 2])
        @testset "Switchings" begin
            @test length(collect(switchings(s, 2, 1, true))) == 4
            @test source.(s, switchings(s, 2, 1, true)) == ones(4)
            @test target.(s, switchings(s, 2, 1, true)) == ones(4)
            @test source.(s, switchings(s, 2, 1, false)) == ones(4)
            @test target.(s, switchings(s, 2, 1, false)) == ones(4)
            sw1 = first(switchings(s, 2, 1, true))
            sw2 = first(switchings(s, 2, 1, false))
            append!(sw1, sw2)
            @test sw1.seq == HybridSystems.OneStateTransition.(ones(Int, 4))
            prepend!(sw1, sw2)
            @test sw1.seq == HybridSystems.OneStateTransition.(ones(Int, 6))

            sw = HybridSystems.switchingsequence(s, 3)
            @test sw.len == 0
            prepend!(sw, sw2)
            @test sw.len == 2
            prepend!(sw, sw2)
            @test sw.len == 4
            @test sw.seq == HybridSystems.OneStateTransition.(ones(Int, 4))

            sw = HybridSystems.switchingsequence(s, 3)
            @test sw.len == 0
            append!(sw, sw2)
            @test sw.len == 2
            append!(sw, sw2)
            @test sw.len == 4
            @test sw.seq == HybridSystems.OneStateTransition.(ones(Int, 4))
        end
    end

    @testset "State Dependent Switched System" begin
        s = discreteswitchedsystem([ones(2, 2), zeros(2, 2)], [square, simplex])
        @test nstates(s) == 2
        @test states(s) == 1:2
        @test all(statedim.(s, states(s)) .== 2)
        @test stateset(s, 1) === square
        @test stateset(s, 2) === simplex
        #@test length(transitions(s)) == 4
        @test ntransitions(s) == 4
    end

    @testset "State Dependent Switched System" begin
        s = discreteswitchedsystem([ones(2, 2), zeros(2, 2)], Fill(square, 2))
        @test nstates(s) == 2
        @test states(s) == 1:2
        @test all(statedim.(s, states(s)) .== 2)
        @test stateset(s, 1) === square
        @test stateset(s, 2) === square
        #@test length(transitions(s)) == 4
        @test ntransitions(s) == 4
    end
end

@testset "Examples" begin
    @testset "Square" begin
        include("../examples/square.jl")
        s = square_example(Polyhedra.DefaultLibrary{Float64}())
        @test nstates(s) == 1
        @test states(s) == 1:1
        @test statedim(s, 1) == 2
        @test stateset(s, 1) isa Polyhedra.DefaultPolyhedron{Float64}
        @test inputdim(s, 1) == 0
    end
    @testset "Horizontal jump" begin
        include("../examples/horizontal_jump.jl")
        for shift in (false, true)
            s = horizontal_jump_example(Polyhedra.DefaultLibrary{Float64}(), shift)
            @test nstates(s) == 2
            @test states(s) == 1:2
            for u in states(s)
                @test mode(s, u) === s.modes[u]
                @test statedim(s, u) == 2
                @test stateset(s, u) isa Polyhedra.DefaultPolyhedron{Float64}
            end
            #@test length(transitions(s)) == 2
            @test ntransitions(s) == 2
            for t in transitions(s)
                @test stateset(s, t) isa Polyhedra.DefaultPolyhedron{Float64}
                @test inputset(s, t) isa Polyhedra.DefaultPolyhedron{Float64}
                @test target_mode(s, t) === s.modes[target(s, t)]
            end
        end
    end
    @testset "Cruise Control" begin
        include("../examples/cruise_control.jl")
        for sym in (false, true)
            s = cruise_control_example(7, 1, sym=sym)
            @test nstates(s) == 7
            @test transitiontype(s) == HybridSystems.LightTransition{LightGraphs.SimpleGraphs.SimpleEdge{Int64}}
            t = first(transitions(s))
            @test target(s, t) == 1
            @test source(s, t) == 1
            #@test sprint(show, s) == "Hybrid System with automaton HybridSystems.LightAutomaton{LightGraphs.SimpleGraphs.SimpleDiGraph{Int64},LightGraphs.SimpleGraphs.SimpleEdge{Int64}}({7, 21} directed simple Int64 graph, Dict(Edge 2 => 3=>2,Edge 6 => 5=>1,Edge 2 => 1=>1,Edge 4 => 4=>1,Edge 4 => 3=>2,Edge 3 => 6=>2,Edge 6 => 3=>2,Edge 4 => 6=>2,Edge 7 => 7=>1,Edge 6 => 6=>2,Edge 3 => 2=>1,Edge 5 => 4=>1,Edge 5 => 6=>2,Edge 7 => 6=>2,Edge 1 => 3=>2,Edge 5 => 3=>2,Edge 2 => 6=>2,Edge 1 => 1=>1,Edge 1 => 6=>2,Edge 3 => 3=>2,Edge 7 => 3=>2))"
            @test iszero(inputdim(s, 1))
            @test iszero(inputdim(s, 2))
        end
    end
end


text:
# categorical contains specific exceptions
# if we cannot map to categorical, we map to linear instead

struct HybridDiscretizer{N<:Real, D<:Integer} <: AbstractDiscretizer{N,D}
    cat :: CategoricalDiscretizer{N,D}
    lin :: LinearDiscretizer{N,D}
end

"""
a hybrid discretizers that maps a special missing value indicator (typically Inf or NaN)
to a discrete bin, but otherwise is a linear discretizer
"""
function datalineardiscretizer(binedges::Vector{Float64}, ::Type{D}=Int;
    missing_key::Float64=Inf,
    force_outliers_to_closest::Bool=DEFAULT_LIN_DISC_FORCE_OUTLIERS_TO_CLOSEST,
    ) where {D<:Integer}

    categorical = CategoricalDiscretizer(Dict{Float64, Int}(missing_key=>1))
    linear = LinearDiscretizer(binedges, D, force_outliers_to_closest=force_outliers_to_closest)
    HybridDiscretizer{Float64,D}(categorical, linear)
end

function supports_encoding(disc::HybridDiscretizer{N,D}, x::N) where {N<:Real,D<:Integer}
    if haskey(disc.cat.n2d, x) ||
        (disc.lin.force_outliers_to_closest && !isnan(x))
        return true
    else
        return disc.lin.binedges[1] ≤ x ≤ disc.lin.binedges[end]
    end
end
supports_decoding(disc::HybridDiscretizer{N,D}, d::D) where {N<:Real,D<:Integer} = 1 ≤ d ≤ nlabels(disc)

function encode(disc::HybridDiscretizer{N,D}, x::N) where {N,D<:Integer}

    if haskey(disc.cat.n2d, x)
        return encode(disc.cat, x) + disc.lin.nbins
    end

    convert(D, encode(disc.lin, x))
end
encode(disc::HybridDiscretizer{N,D}, x) where {N,D} = encode(disc, convert(N, x))::D
function encode(disc::HybridDiscretizer{N,D}, data::AbstractArray) where {N,D<:Integer}
    arr = Vector{D}(undef, length(data))
    for (i,x) in enumerate(data)
        arr[i] = encode(disc, x)
    end
    reshape(arr, size(data))
end

function decode(disc::HybridDiscretizer{N,D}, d::D, method::AbstractSampleMethod=SAMPLE_UNIFORM)::N where {N<:Real,D<:Integer}
    if d ≤ disc.lin.nbins
        retval = decode(disc.lin, d, method)
    else
        retval = decode(disc.cat,  d - disc.lin.nbins)
    end
    return retval
end
decode(disc::HybridDiscretizer{N,D}, d::E, method::AbstractSampleMethod=SAMPLE_UNIFORM) where {N<:Real,D<:Integer,E<:Integer} =
    decode(disc, convert(D, d), method)

function decode(disc::HybridDiscretizer{N,D}, data::AbstractArray{D}, ::AbstractSampleMethod=SAMPLE_UNIFORM) where {N,D<:Integer}
    arr = Vector{N}(undef, length(data))
    for (i,d) in enumerate(data)
        arr[i] = decode(disc, d)
    end
    reshape(arr, size(data))
end

Base.max(disc::HybridDiscretizer) = Base.max(disc.lin)
Base.min(disc::HybridDiscretizer) = Base.min(disc.lin)
Base.extrema(disc::HybridDiscretizer) = extrema(disc.lin)
Base.extrema(disc::HybridDiscretizer{N,D}, d::D) where {N<:AbstractFloat,D} = extrema(disc.lin, d)
Base.extrema(disc::HybridDiscretizer{N,D}, d::D) where {N<:Integer,D} = extrema(disc.lin, d)
totalwidth(disc::HybridDiscretizer) = totalwidth(disc.lin)

nlabels(disc::HybridDiscretizer) = disc.lin.nbins + nlabels(disc.cat)
bincenters(disc::HybridDiscretizer) = bincenters(disc.lin)
binwidth(disc::HybridDiscretizer{N,D}, d::D) where {N<:AbstractFloat,D} = binwidth(disc.lin, d)
binwidths(disc::HybridDiscretizer{N,D}) where {N<:AbstractFloat,D} = binwidths(disc.lin)


text:
@testset "AlgAssRelOrdIdl" begin
  Qx, x = FlintQQ["x"]
  f = x^2 - 10x - 8
  K, a = number_field(f, "a")
  KG = group_algebra(K, GrpAbFinGen([ 2 ]))

  @testset "Arithmetic" begin
    O = any_order(KG)
    I = 2*O
    J = 4*O

    @test I + J == I
    @test I*J == 8*O
    @test intersect(I, J) == J
    @test I^2 == J
    @test I^fmpz(2) == J

    @test norm(I) == 4*base_ring(O)
  end

  @testset "Locally free basis" begin
    Ky, y = K["y"]
    OK = maximal_order(K)
    g = y^2 + 1
    L, b = number_field(g, "b") # Gal(L/K) == C_2
    OL = maximal_order(L)

    LtoKG, KGtoL = Hecke._find_isomorphism(L, KG)
    basisOL = Vector{elem_type(KG)}()
    for i = 1:degree(L)
      for b in basis(pseudo_basis(OL, copy = false)[i][2])
        push!(basisOL, LtoKG(b*pseudo_basis(OL, copy = false)[i][1]))
      end
    end
    d = lcm([ denominator(b) for b in basisOL ])
    OKG = Order(KG, basis(KG))
    I = Hecke.ideal_from_lattice_gens(KG, OKG, [ d*b for b in basisOL ])

    p = prime_decomposition(OK, 3)[1][1]
    g = Hecke.locally_free_basis(I, p)
    OKGg = OKG*g
    mat_I = Hecke.coprime_bases(OKG, I, p)[4]
    mat_OKGg = Hecke.coprime_bases(OKG, OKGg, p)[4]
    t = det(mat_I*inv(mat_OKGg))
    @test valuation(t, p) == 0

    p = prime_decomposition(OK, 2)[1][1]
    @test_throws ErrorException  Hecke.locally_free_basis(I, p)
  end

  @testset "maximal integral ideal" begin
    K,z = cyclotomic_field(7)
    o = maximal_order(K)
    M = zeros(K,2,2,2)
    M[1,1,1] = one(K)
    M[1,2,2] = one(K)
    M[2,1,2] = one(K)
    M[2,2,1] = K(-12)
    E = Hecke.AlgAss(K,M)
    OE = maximal_order(E)
    @test isprime(numerator(norm(@inferred Hecke.maximal_integral_ideal(OE, 3*o, :left))))
  end

  @testset "factor" begin
    K,z = cyclotomic_field(7);
    E = Hecke.quaternion_algebra2(K, K(-1), K(-1))
    OE = maximal_order(E)
    list_ideals = [rand(E)*OE for i=1:10]
    list_factors = factor.(list_ideals)
    @test prod.(list_factors) == list_ideals
  end
end


text:
function act(T::TensorMap{ComplexSpace, 2, 2}, psi::TensorMap{ComplexSpace, 2, 1})
    chi = get_chi(psi)
    d = dims(domain(T))[2]
    fuse_tensor = isomorphism(ℂ^(chi*d), ℂ^d*ℂ^chi)

    @tensor Tpsi[-1, -2; -3] := T[3, -2, 2, 5] * psi[1, 2, 4] * fuse_tensor[-1, 3, 1] * fuse_tensor'[5, 4, -3]
    return Tpsi
end

function rrule(::typeof(act), T::TensorMap{ComplexSpace, 2, 2}, psi::TensorMap{ComplexSpace, 2, 1})
    chi = get_chi(psi)
    d = dims(domain(T))[2]
    fuse_tensor = isomorphism(ℂ^(chi*d), ℂ^d*ℂ^chi)

    @tensor Tpsi[-1, -2; -3] := T[3, -2, 2, 5] * psi[1, 2, 4] * fuse_tensor[-1, 3, 1] * fuse_tensor'[5, 4, -3]

    fwd = Tpsi
    function act_pushback(f̄wd)
        @tensor psi_pushback[-1, -2; -3] := f̄wd[1,4,2] * conj(fuse_tensor[1,3,-1]) * conj(fuse_tensor'[5,-3,2]) * conj(T[3,4,-2,5])
        return NoTangent(), NoTangent(), psi_pushback
    end
    return fwd, act_pushback
end

function transf_mat(psi::TensorMap{ComplexSpace, 2, 1}, phi::TensorMap{ComplexSpace, 2, 1})
    function lop(v::TensorMap{ComplexSpace, 1, 1})
        @tensor result[-2; -1] := psi'[1, -1, 3] * phi[-2, 3, 2] * v[2, 1]
        return result
    end
    return lop
end

function transf_mat_T(psi::TensorMap{ComplexSpace, 2, 1}, phi::TensorMap{ComplexSpace, 2, 1})
    function lop_T(v::TensorMap{ComplexSpace, 1, 1})
        @tensor result[-2; -1] := psi[1, 3, -1] * phi'[-2, 2, 3] * v[2, 1]
        return result
    end
    return lop_T
end

function transf_mat(A::TensorMap{ComplexSpace, 2, 1}, O::TensorMap{ComplexSpace, 2, 2}, B::TensorMap{ComplexSpace, 2, 1})
    function lop(v::TensorMap{ComplexSpace, 2, 1})
        @tensor v1[-1, -2; -3] := A[-1, 5, 1] * O[-2, 4, 5, 3] * B'[2, -3, 4] * v[1, 3, 2]
        return v1
    end
    return lop
end

function transf_mat_T(A::TensorMap{ComplexSpace, 2, 1}, O::TensorMap{ComplexSpace, 2, 2}, B::TensorMap{ComplexSpace, 2, 1})
    function lop_T(v::TensorMap{ComplexSpace, 2, 1})
        @tensor v1[-1, -2; -3] := A'[-1, 1, 5] * O'[5, -2, 3, 4] * B[2, 4, -3] * v[1, 3, 2]
        return v1
    end
    return lop_T
end

function ovlp(psi::TensorMap{ComplexSpace, 2, 1}, phi::TensorMap{ComplexSpace, 2, 1})
    chi_psi, chi_phi = get_chi(psi), get_chi(phi)
    v0 = TensorMap(rand, ComplexF64, ℂ^chi_phi, ℂ^chi_psi)
    w, _ = eigsolve(transf_mat(psi, phi), v0, 1)
    return w[1]
end

function rrule(::typeof(ovlp), psi::TensorMap{ComplexSpace, 2, 1}, phi::TensorMap{ComplexSpace, 2, 1})
    chi_psi, chi_phi = get_chi(psi), get_chi(phi)
    v0 = TensorMap(rand, ComplexF64, ℂ^chi_phi, ℂ^chi_psi)
    wr, vr = eigsolve(transf_mat(psi, phi), v0, 1)
    _, vl = eigsolve(transf_mat_T(psi, phi), v0, 1)
    fwd = wr[1]
    vl, vr = vl[1], vr[1]
    vr = vr / dot(vl, vr)

    function ovlp_pushback(f̄wd)
        @tensor psi_pushback[-2, -3; -1] := vl'[-2, 1] * phi[1, -3, 2] * vr[2, -1] * conj(f̄wd)
        @tensor phi_pushback[-1, -2; -3] := conj(vl'[1, -1]) * conj(psi'[2, 1, -2]) * conj(vr[-3, 2]) * f̄wd
        return NoTangent(), psi_pushback, phi_pushback
    end
    return fwd, ovlp_pushback
end

function nonherm_cost_func(T::TensorMap{ComplexSpace, 2, 2}, psi::TensorMap{ComplexSpace, 2, 1})
    Tpsi = act(T, psi)

    up = ovlp(Tpsi, Tpsi) * ovlp(psi, psi)
    dn = ovlp(psi, Tpsi) * ovlp(Tpsi, psi)

    return log(norm(up / dn))
end

function lambda_gamma(psi::TensorMap{ComplexSpace, 2, 1})
    chi = get_chi(psi)

    lop = transf_mat(psi, psi)
    lopT = transf_mat_T(psi, psi)

    ρ0 = TensorMap(rand, ComplexF64, ℂ^chi, ℂ^chi)
    wr, ρr = eigsolve(lop, ρ0, 1)
    _, ρl = eigsolve(lopT, ρ0, 1)
    ρl, ρr = ρl[1], ρr[1]
    norm_psi = wr[1] |> norm |> sqrt
    psi = psi / norm_psi

    _, s, v = tsvd(ρr)
    Y = sqrt(s) * v
    Yinv = v' * inv(sqrt(s))

    _, s, v = tsvd(ρl')
    X = sqrt(s) * v
    Xinv = v' * inv(sqrt(s))

    u, Λ, v = tsvd(X * Y')
    @tensor Γ[-1, -2; -3] :=
        v[-1, 1] * Yinv'[1, 3] * psi[3, -2, 4] * Xinv[4, 2] * u[2, -3]
    return norm_psi, Γ, Λ
end

function left_canonical(psi::TensorMap{ComplexSpace, 2, 1})

    norm_psi, Γ, Λ = lambda_gamma(psi)
    @tensor psi_L[-1, -2; -3] := Λ[-1, 1] * Γ[1, -2, -3]
    return norm_psi, psi_L

end

function iTEBD_truncate(psi::TensorMap{ComplexSpace, 2, 1}, chi::Integer)
    norm_psi, Γ, Λ = lambda_gamma(psi)
    chi0 = get_chi(psi)

    P = isometry(ℂ^chi0, ℂ^chi)
    @tensor psi_truncated[-1, -2; -3] :=
        P'[-1, 1] * Λ[1, 2] * Γ[2, -2, 3] * P[3, -3]
    return norm_psi, psi_truncated
end

function ln_fidelity(psi::TensorMap{ComplexSpace, 2, 1}, phi::TensorMap{ComplexSpace, 2, 1})
    up = ovlp(psi, phi) * ovlp(phi, psi)
    dn = ovlp(psi, psi) * ovlp(phi, phi)
    return (up/dn) |> norm |> log
end

function free_energy(T::TensorMap{ComplexSpace, 2, 2}, psi::TensorMap{ComplexSpace, 2, 1})
    Tpsi = act(T, psi)
    up = ovlp(psi, Tpsi)
    dn = ovlp(psi, psi)
    return (up/dn) |> norm |> log
end

function mps_add(psi::TensorMap{ComplexSpace, 2, 1}, phi::TensorMap{ComplexSpace, 2, 1})
    # https://github.com/Jutho/TensorKit.jl/issues/54
    # actually doesn't work for infinite MPS
    chi_psi, chi_phi = get_chi(psi), get_chi(phi)
    chi = chi_psi + chi_phi

    embedder1 = isometry(ℂ^chi, ℂ^chi_psi)
    embedder2 = TensorMap(zeros, ComplexF64, ℂ^chi, ℂ^chi_phi)
    embedder2.data[chi_psi+1:chi, :] = Matrix{ComplexF64}(I, chi_phi, chi_phi)

    @tensor psi_plus_phi[-1, -2; -3] := embedder1[-1, 1] * psi[1, -2, 2] * embedder1'[2, -3] +
             embedder2[-1, 1] * phi[1, -2, 2] * embedder2'[2, -3]

    return psi_plus_phi
end

function right_canonical_QR(psi::TensorMap{ComplexSpace, 2, 1}, tol::Float64=1e-15)

    chi = get_chi(psi)
    L0 = id(ℂ^chi)

    L, Q = rightorth(permute(psi, (1, ), (2, 3)))
    psi_R = permute(Q, (1, 2), (3, ))
    L = L / norm(L)
    δ = norm(L - L0)
    L0 = L

    ix= 0
    while δ > tol && ix < 200
        if ix >= 20 && ix % 10 == 0
            lop = transf_mat(psi, psi_R)
            _, vr = eigsolve(lop, L0, 1; tol=max(tol, δ/10))
            L = vr[1]'
            #print("--> ")
        end

        L, Q = rightorth(permute(psi * L, (1, ), (2, 3)))
        psi_R = permute(Q, (1, 2), (3, ))
        L = L / norm(L)

        δ = norm(L-L0)
        L0 = L

        #println(δ)
        ix += 1
    end

    #println(ix, " iterations")
    δ > tol && @warn "right_canonical_QR failed to converge. δ: $δ , tol: $tol"

    return L0', psi_R
end

function left_canonical_QR(psi::TensorMap{ComplexSpace, 2, 1}, tol::Float64=1e-15)
    chi = get_chi(psi)
    R0 = id(ℂ^chi)

    psi_L, R = leftorth(psi)
    R = R / norm(R)
    δ = norm(R - R0)
    R0 = R

    ix = 0
    while δ > tol && ix < 200
        if ix >= 20 && ix % 10 == 0
            lop_T = transf_mat_T(psi, psi_L)
            _, vl = eigsolve(lop_T, R0, 1; tol=max(tol, δ/10))
            R = vl[1]
            #print("--> ")
        end

        @tensor psi_tmp[-1, -2; -3] := R[-1, 1] * psi[1, -2, -3]
        psi_L, R = leftorth(psi_tmp)
        R = R / norm(R)

        δ = norm(R - R0)
        R0 = R

        #println(δ)
        ix += 1
    end

    #println(ix, " iterations")
    δ > tol && @warn "left_canonical_QR failed to converge. δ: $δ , tol: $tol"

    return R0, psi_L
end

function entanglement_spectrum(psi::TensorMap{ComplexSpace, 2, 1})
    _, psi = left_canonical_QR(psi)
    Y, _ = right_canonical_QR(psi)
    _, s, _ = tsvd(Y')
    s = diag(s.data) .^ 2
    return s ./ sum(s)
end

function entanglement_entropy(psi::TensorMap{ComplexSpace, 2, 1})
    s = entanglement_spectrum(psi)
    return -sum(s .* log.(s))
end

function expand(psi::TensorMap{ComplexSpace, 2, 1}, chi::Integer, perturb::Float64=1e-3)
    chi0, d = get_chi(psi), get_d(psi)
    if chi <= chi0
        @warn "chi not larger than current bond D, not expanded "
        return psi
    end

    phi_arr = perturb*rand(ComplexF64, chi, d, chi)
    phi_arr[1:chi0, :, 1:chi0] += toarray(psi)
    phi = arr_to_TensorMap(phi_arr)

    return phi
end

function tangent_map_tn(O::TensorMap{ComplexSpace, 2, 2}, AL::TensorMap{ComplexSpace, 2, 1}, AR::TensorMap{ComplexSpace, 2, 1})
    lop_R = transf_mat(AR, O, AR)
    lop_L = transf_mat_T(AL, O, AL)
    chi_mps = get_chi(AL)
    chi_mpo = dims(domain(O))[2]

    ER = TensorMap(rand, ComplexF64, ℂ^chi_mps*ℂ^chi_mpo, ℂ^chi_mps)
    EL = TensorMap(rand, ComplexF64, ℂ^chi_mps*ℂ^chi_mpo, ℂ^chi_mps)
    _, ER = eigsolve(lop_R, ER, 1)
    _, EL = eigsolve(lop_L, EL, 1)
    EL, ER = EL[1], ER[1]

    @tensor norm = ER[1,2,3]*EL'[3,1,2]
    ER = ER / norm

    @tensor map_AC[-1, -2, -3; -4, -5, -6] := EL'[-1, -4, 1] * O[1, -2, -5, 2] * ER[-6, 2, -3]
    @tensor map_C[-1, -2; -3, -4] := EL'[-1, -3, 1] * ER[-4, 1, -2]

    return map_AC, map_C
end

function tangent_map(O::TensorMap{ComplexSpace, 2, 2}, AL::TensorMap{ComplexSpace, 2, 1}, AR::TensorMap{ComplexSpace, 2, 1})
    lop_R = transf_mat(AR, O, AR)
    lop_L = transf_mat_T(AL, O, AL)
    chi_mps = get_chi(AL)
    chi_mpo = dims(domain(O))[2]

    ER = TensorMap(rand, ComplexF64, ℂ^chi_mps*ℂ^chi_mpo, ℂ^chi_mps)
    EL = TensorMap(rand, ComplexF64, ℂ^chi_mps*ℂ^chi_mpo, ℂ^chi_mps)
    _, ER = eigsolve(lop_R, ER, 1)
    _, EL = eigsolve(lop_L, EL, 1)
    EL, ER = EL[1], ER[1]

    @tensor norm = ER[1,2,3]*EL'[3,1,2]
    ER = ER / norm

    function map_AC(AC::TensorMap{ComplexSpace, 2, 1})
        @tensor updated_AC[-1, -2; -3] := AC[1, 3, 2] * EL'[-1, 1, 4] * ER[2, 5, -3] * O[4, -2, 3, 5]
        return updated_AC
    end
    function map_C(C::TensorMap{ComplexSpace, 1, 1})
        @tensor updated_C[-1; -2] := C[1, 2] * EL'[-1, 1, 3] * ER[2, 3, -2]
        return updated_C
    end
    return map_AC, map_C
end

"""
    calculate_ALR(AC::TensorMap{ComplexSpace, 2, 1}, C::TensorMap{ComplexSpace, 1, 1}) -> AL::TensorMap{ComplexSpace, 2, 1}, AR::TensorMap{ComplexSpace, 2, 1}, ϵL::Float64, ϵR::Float64

    From `AC` and `C`, calculate the AL and AR tensors. Calculate the convergence measures `ϵL` and `ϵR` as well.
"""
function calculate_ALR(AC::TensorMap{ComplexSpace, 2, 1}, C::TensorMap{ComplexSpace, 1, 1})
    UAC_l, RAC = leftorth(AC; alg=Polar())
    UC_l, RC = leftorth(C; alg=Polar())

    LAC, UAC_r = rightorth(permute(AC, (1,), (2,3)); alg=Polar())
    LC, UC_r = rightorth(C; alg=Polar())

    AL = UAC_l * UC_l'
    AR = permute(UC_r' * UAC_r, (1, 2), (3,))

    ϵL = norm(RAC - RC)
    ϵR = norm(LAC - LC)
    #ϵL = norm(AC - AL * C)
    #ϵR = norm(permute(AC, (1,), (2, 3)) - C * permute(AR, (1,), (2, 3)))

    return AL, AR, ϵL, ϵR
end

"""
    biorth_right_canonical(AU::TensorMap{ComplexSpace, 2, 1}, AD::TensorMap{ComplexSpace, 2, 1}) -> AU_R::TensorMap{ComplexSpace, 2, 1}, AD_R::TensorMap{ComplexSpace, 2, 1}

    right biorthogonalization of MPSes `AU` and `AD` into `AU_R` and `AD_R`
"""
function biorth_right_canonical(AU::TensorMap{ComplexSpace, 2, 1}, AD::TensorMap{ComplexSpace, 2, 1})
    chiU, chiD = get_chi(AU), get_chi(AD)

    lop = transf_mat(AU, AD)
    wR, CR = eigsolve(lop, TensorMap(rand, ComplexF64, ℂ^chiD, ℂ^chiU), 1)
    wR, CR = wR[1], CR[1]

    #UR, S, VR = tsvd(CR)
    #S = sqrt(S)

    #wR_abs, wR_arg = abs(wR), angle(wR)
    #Sinv = pinv(S)
    #@tensor AU_R[-1, -2; -3] := Sinv[-1, 1] * VR[1, 2] * AU[2, -2, 3] * VR'[3, 4] * S[4, -3] / sqrt(wR_abs)
    #@tensor AD_R[-1, -2; -3] := Sinv[-1, 1] * UR'[1, 2] * AD[2, -2, 3] * UR[3, 4] * S[4, -3] / sqrt(wR_abs) * exp(-1im*wR_arg)

    #return AU_R, AD_R

    @tensor AD_R[-1, -2; -3] := pinv(CR)[-1, 1] * AD[1, -2, 2] * CR[2, -3] / wR
    return AU, AD_R
end

"""
    biorth_left_canonical(AU::TensorMap{ComplexSpace, 2, 1}, AD::TensorMap{ComplexSpace, 2, 1}) -> AU_L::TensorMap{ComplexSpace, 2, 1}, AD_L::TensorMap{ComplexSpace, 2, 1}

    left biorthogonalization of MPSes `AU` and `AD` into `AU_L` and `AD_L`
"""
function biorth_left_canonical(AU::TensorMap{ComplexSpace, 2, 1}, AD::TensorMap{ComplexSpace, 2, 1})
    chiU, chiD = get_chi(AU), get_chi(AD)

    lopT = transf_mat_T(AU, AD)
    wL, CL = eigsolve(lopT, TensorMap(rand, ComplexF64, ℂ^chiD, ℂ^chiU), 1)
    wL, CL = wL[1], CL[1]

    #UL, S, VL = tsvd(CL')
    #S = sqrt(S)

    #wL_abs, wL_arg = abs(wL), angle(wL)
    #Sinv = pinv(S)

    #@tensor AU_L[-1, -2; -3] := S[-1, 1] * UL'[1, 2] * AU[2, -2, 3] * UL[3, 4] * Sinv[4, -3] / sqrt(wL_abs)
    #@tensor AD_L[-1, -2; -3] := S[-1, 1] * VL[1, 2] * AD[2, -2, 3] * VL'[3, 4] * Sinv[4, -3] / sqrt(wL_abs) * exp(1im*wL_arg)

    @tensor AD_L[-1, -2; -3] := CL'[-1, 1] * AD[1, -2, 2] * pinv(CL')[2, -3] / wL'
    return AU, AD_L
end

# FIXME. correct this.
function tangent_map(O::TensorMap{ComplexSpace, 2, 2}, AL::TensorMap{ComplexSpace, 2, 1}, AR::TensorMap{ComplexSpace, 2, 1}, BL::TensorMap{ComplexSpace, 2, 1}, BR::TensorMap{ComplexSpace, 2, 1})
    lop_R = transf_mat(AR, O, BR)
    lop_L = transf_mat_T(AL, O, BL)
    chi_mpsA, chi_mpsB = get_chi(AL), get_chi(BL)
    chi_mpo = dims(domain(O))[2]

    ER = TensorMap(rand, ComplexF64, ℂ^chi_mpsA*ℂ^chi_mpo, ℂ^chi_mpsB)
    EL = TensorMap(rand, ComplexF64, ℂ^chi_mpsA*ℂ^chi_mpo, ℂ^chi_mpsB)
    _, ER = eigsolve(lop_R, ER, 1)
    _, EL = eigsolve(lop_L, EL, 1)
    EL, ER = EL[1], ER[1]

    @tensor norm = ER[1,2,3]*EL'[3,1,2]
    ER = ER / norm

    function map_AC(AC::TensorMap{ComplexSpace, 2, 1})
        @tensor updated_AC[-1, -2; -3] := AC[1, 3, 2] * EL'[-1, 1, 4] * ER[2, 5, -3] * O[4, -2, 3, 5]
        return updated_AC
    end
    function map_C(C::TensorMap{ComplexSpace, 1, 1})
        @tensor updated_C[-1; -2] := C[1, 2] * EL'[-1, 1, 3] * ER[2, 3, -2]
        return updated_C
    end
    return map_AC, map_C
end


text:
@safetestset "OPTIONS requests" begin

  using Genie, HTTP

  route("/options", method = OPTIONS) do
    push!(params(:RESPONSE).headers, "X-Foo-Bar" => "Baz")
  end

  port = nothing
  port = rand(8500:8900)

  server = up(port)
  sleep(1)

  response = HTTP.request("OPTIONS", "http://localhost:$port") # unhandled, should get default response
  @test response.status == 200
  @test get(Dict(response.headers), "X-Foo-Bar", nothing) == nothing

  response = HTTP.request("OPTIONS", "http://localhost:$port/options") # handled
  @test response.status == 200
  @test get(Dict(response.headers), "X-Foo-Bar", nothing) == "Baz"

  down()
  sleep(1)
  server = nothing
  port = nothing
end


text:
npts = 64;

## METHOD 1

function cylinder(;r=1, h=1, n=8, k=1)
	V,FV = Lar.ring(0,r)([n,r])
	W,CW = Lar.extrudeSimplicial((V,FV), [h for i=1:k])
	V = W .- [0,0,h*k/2]
	FW = Lar.simplexFacets(CW)
	FV = Lar.simplexBoundary_3(CW,FW)
	EV = [[u,v] for (u,v) in sort(collect(Set(map(sort,cat(
					[[[u,v],[v,w],[w,u]] for (u,v,w) in FV]))))) if u!=v]
	return V,FV,EV
end

tube = cylinder(n=npts,h=2,k=2);

triple = Lar.Struct([ tube, Lar.r(pi/2,0,0), tube, Lar.r(0,pi/2,0), tube ]);
model = Lar.struct2lar(triple)
V,FV,EV = model
GL.VIEW([
	GL.GLFrame,
	GL.GLGrid( model[1:2]..., GL.COLORS[1], 0.5)
]);

model = CAGD.Model(V);
cFE = convert(Lar.ChainOp, Lar.coboundary_1(model.G, FV, EV));
CAGD.addModelCells!(model, 1, EV, signed = true);
CAGD.addModelCells!(model, 2, cFE);

A_model = CAGD.spatial_arrangement(model);

if todisplay
    V, CVs, FVs, EVs = Lar.pols2tria(A_model.G, A_model.T[1], A_model.T[2], A_model.T[3]);

    GL.VIEW(GL.GLExplode(V, FVs,  1.5,1.5,1.5,  99,1));
    GL.VIEW(GL.GLExplode(V, EVs,  1.5,1.5,1.5,  99,1));
    GL.VIEW(GL.GLExplode(V, CVs[1:end],  5,5,5,  99,0.5));
end;

## METHOD 2

# Build a cylinder

V,_ = Lar.rod()([npts, 1]);
EV = [
    [[2*i, (2*i+1)%2npts+1] for i = 1 : npts];          # horizontal upper edges
    [[2*i-1, (2*i)%2npts+1] for i = 1 : npts];          # horizontal lower edges
    [[2*i+1, 2*i+2] for i = 0 : npts-1];                # vertical edges
];
FE = [
    [i for i = 1 : npts],                                       # upper face
    [npts + i for i = 1 : npts],                                # lower face
    [[i, npts+i, 2*npts+i, 2*npts+i+1] for i = 1:npts-1]...,    # vertical faces
    [npts, 2*npts, 3*npts, 2*npts+1]
];

FV = Lar.cop2lar(Lar.lar2cop(FE) * Lar.lar2cop(EV));

# Build three cylinders

cyl = Lar.Struct([
    Lar.r(0,0,0),
    Lar.s(1.0, 1.0, 1.5),
    Lar.t(0,0,-1.5),
    (V,FV,EV)
]);
cyl = Lar.struct2lar(cyl);
tris = Lar.Struct([
    cyl,
    Lar.Struct([Lar.r(pi/2,0,0), cyl ]) ,
    Lar.Struct([Lar.r(0,pi/2,0), cyl ])
]);
V, FV, EV = Lar.struct2lar(tris);

if todisplay
    GL.VIEW([
        GL.GLFrame
        GL.GLGrid(V,EV,GL.COLORS[1],1)
    ]);
end;

# Model Generation

model = CAGD.Model(V);
cFE = convert(Lar.ChainOp, Lar.coboundary_1(model.G, FV, EV));
CAGD.addModelCells!(model, 1, EV, signed = true);
CAGD.addModelCells!(model, 2, cFE);

A_model = CAGD.spatial_arrangement(model);

if todisplay
    V, CVs, FVs, EVs = Lar.pols2tria(A_model.G, A_model.T[1], A_model.T[2], A_model.T[3])

    GL.VIEW(GL.GLExplode(V, FVs,  1.5,1.5,1.5,  99,1));
    GL.VIEW(GL.GLExplode(V, EVs,  1.5,1.5,1.5,  99,1));
    GL.VIEW(GL.GLExplode(V, CVs[1:end],  5,5,5,  99,0.5));
end;


text:
# Apply a log-transform on a Gamma-distributed variable to sample from it using a normal proposal

using Klara

plogtarget(p::Float64, v::Vector) = (v[1]-1)*p-exp(p)/v[2]+p

p = BasicContUnvParameter(:p, logtarget=plogtarget, nkeys=3)

model = likelihood_model([Constant(:k), Constant(:θ), p], isindexed=false)

mcsampler = MH()

mcrange = BasicMCRange(nsteps=100000, burnin=10000)

v0 = Dict(:k=>2., :θ=>1., :p=>log(10.))

job = BasicMCJob(model, mcsampler, mcrange, v0)

run(job)

chain = output(job)

mean(exp.(chain.value))

acceptance(chain, diagnostics=false)


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    AbstractNeighborhood

A neighborhood of points in a `N`-dimensional space
with coordinates of type `T`.
"""
abstract type AbstractNeighborhood{T,N} end

"""
    ndims(neigh)

Return the number of dimensions of neighborhood `neigh`.
"""
Base.ndims(::AbstractNeighborhood{T,N}) where {N,T} = N

"""
    coordtype(neigh)

Return the coordinate type of neighborhood `neigh`.
"""
coordtype(::AbstractNeighborhood{T,N}) where {N,T} = T

"""
    isneighbor(neigh, xₒ, x)

Tells whether or not the coordinates `x` are in the neighborhood
`neigh` centered at coordinates `xₒ`.
"""
isneighbor(neigh::AbstractNeighborhood,
           xₒ::AbstractVector, x::AbstractVector) = @error "not implemented"

"""
    volume(neigh)

Return the volume of neighborhood `neigh`.
"""
volume(neigh::AbstractNeighborhood) = @error "not implemented"

#------------------
# IMPLEMENTATIONS
#------------------
include("neighborhoods/ball.jl")
include("neighborhoods/cylinder.jl")


text:
n = 10
while n > 0
    println("Countdown: $n")
    global n = n - 1 # Loop introduces new scope
end

# println("Loop terminated.")

function countdown(n)
    while n > 0
        println("Countdown: $n")
        n = n - 1 # Within the same scope
    end
    println("Loop terminated.")
end

countdown(n)


text:
# adapted from https://gist.github.com/Keno/d0c2df947f67be543036238a0caeb1c6

let counter = 0
    global fast_sin, get_fast_sin_counter, reset_counter
    fast_sin(x::Float64) = (counter += 1; sin(x::Float64))
    get_fast_sin_counter() = counter
    reset_counter() = (counter = 0)
end

# implementation
# ==============

import Core:
    MethodInstance
const CC = Core.Compiler
import .CC:
    AbstractInterpreter,
    NativeInterpreter,
    InferenceState,
    OptimizationState,
    OptimizationParams,
    IRCode,
    optimize,
    argextype,
    widenconst
import Base:
    @invoke,
    get_world_counter,
    to_tuple_type
import Base.Meta: isexpr
using XCompiler

# implements very inefficient rewrite of `sin(::Float64)` -> `fast_sin(::Float64)`;
# ignores anything involved with inlining and discards static calls — this is just PoC
struct FastSinInterpreter <: XInterpreter
    native::NativeInterpreter
    cache_key::UInt
    function FastSinInterpreter(world = get_world_counter(); cache_key = 0, inlining=false, kwargs...)
        native = NativeInterpreter(world;
                                   inf_params = XCompiler.XInferenceParams(; kwargs...),
                                   opt_params = XCompiler.XOptimizationParams(; inlining, kwargs...),
                                   )
        cache_key = convert(UInt, cache_key)
        cache_key = hash(native.inf_params, cache_key)
        cache_key = hash(native.opt_params, cache_key)
        XCompiler.maybe_init_cache!(cache_key)

        return new(native, cache_key)
    end
end
XCompiler.native(interp::FastSinInterpreter)    = interp.native
XCompiler.cache_key(interp::FastSinInterpreter) = interp.cache_key

let
    tme = first(methods(fast_sin, (Float64,)))
    ttt = to_tuple_type((typeof(fast_sin), Float64,))

    function CC.optimize(interp::FastSinInterpreter, opt::OptimizationState, params::OptimizationParams, @nospecialize(result))
        @assert isnothing(opt.ir)

        linfo = opt.linfo
        if !(linfo.def === tme && linfo.specTypes === ttt)
            (; src, sptypes, slottypes) = opt
            for (i, x) in enumerate(src.code)
                if isexpr(x, :call) && length(x.args) == 2
                    ft = widenconst(argextype(x.args[1], src, sptypes, slottypes))
                    if ft === typeof(sin)
                        at = widenconst(argextype(x.args[2], src, sptypes, slottypes))
                        if at === Float64
                            src.code[i] = Expr(:call, GlobalRef(@__MODULE__, :fast_sin), x.args[2])
                        end
                    end
                end
            end
        end

        return optimize(native(interp), opt, params, result)
    end
end

# test
# ====

using Test

get_sin() = sin
function f(x, replace)
    reset_counter()
    @testset "simple" begin
        sin(x)
        @test get_fast_sin_counter() == (replace ? 1 : 0)
    end
    @testset "a bit complex but still inferred" begin
        get_sin()(x)
        @test get_fast_sin_counter() == (replace ? 2 : 0)
    end
    @testset "dynamic dispatch" begin
        get_sin()(Base.inferencebarrier(x))
        @test get_fast_sin_counter() == (replace ? 3 : 0) # fail, we can't hijack the dynamic dispatch
    end
end

@testset "testset" begin
    @testset "customized compilation" begin
        with(FastSinInterpreter()) do
            f(1.0, true)
        end
    end

    @testset "don't affect native code cache" begin
        f(1.0, false) # fail, since code cache can be inserted outside of `CC.code_cache`
    end
end


text:
using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra

@bench @testset "fiber native" begin
    nep = nep_gallery("nlevp_native_fiber")
    n = size(nep, 1)

    # An exact eigenvalue according (reported in NLEVP collection)
    sol_val= 7.139494306065948e-07;

    (λ,v)=quasinewton(nep,λ=7.14e-7,v=ones(n),
                      logger=displaylevel, armijo_factor=0.5,armijo_max=10)

    @test abs(λ-sol_val)<1e-10;

    # check that we "maintain" real arithmetic
    vv=real(v/v[1]);
    (λ1,v)=resinv(Float64,nep,λ=7.14e-7,v=vv,
                      logger=displaylevel)
    @test abs(λ-sol_val)<1e-10;

    @test eltype(v)==Float64

    @testset "Errors thrown" begin
        @test_throws MethodError nep_gallery("nlevp_native_fiber", 15)
        @test_throws MethodError nep_gallery("nlevp_native_fiber", t=15)
    end

end


text:
"""
Saves a global timestamp for a new run, so that all results from this run have the
same time stamp
"""
function new_run()
    open(joinpath(@__DIR__, "global_timestemp.jls"), "w") do io
        serialize(io, now())
    end
end

"""
Queries the time stamp of the last run
"""
function last_time_stamp()
    open(joinpath(@__DIR__, "global_timestemp.jls")) do io
        deserialize(io)
    end
end
"""
Database schema.
"""
immutable BenchResult
    name::String
    benchmark::BenchmarkTools.Trial
    N::Int
    typ::DataType
    device::String
    hardware::String
    codepath::String
    meandiffrence::Float64
    timestamp::DateTime
end

function Base.show(io::IO, b::BenchResult)
    println(io, "BenchResult:")
    println(io, "    ", b.name)
    print(io, "    ")
    show(io, b.benchmark)
    println(io)
    println(io, "    device: ", b.device)
    println(io, "    hardware: ", b.hardware)
    println(io, "    meandiff: ", b.meandiffrence)
    println(io, "    time: ", b.timestamp)
end

function BenchResult(
        x::BenchResult;
        name = x.name,
        benchmark = x.benchmark,
        N = x.N,
        typ = x.typ,
        device = x.device,
        hardware = x.hardware,
        codepath = x.codepath,
        meandiffrence = x.meandiffrence
    )
    BenchResult(
        name,
        benchmark,
        N,
        typ,
        device,
        hardware,
        codepath,
        meandiffrence
    )
end

function BenchResult(
        name::String,
        benchmark::BenchmarkTools.Trial,
        N::Int,
        typ::DataType,
        device::String,
        hardware::String,
        codepath::String,
        meandiff
    )
    BenchResult(
        name,
        benchmark,
        N,
        typ,
        device,
        hardware,
        codepath,
        meandiff,
        last_time_stamp()
    )
end

for field in fieldnames(BenchResult)
    constr_expr = Expr(:call, BenchResult)
    for field2 in fieldnames(BenchResult)
        if field2 == field
            push!(constr_expr.args, :(cval))
        else
            push!(constr_expr.args, :(x.$(field2)))
        end
    end
    @eval begin
        $field(x::BenchResult) = x.$field
        function $field(x::BenchResult, val)
            cval = convert($(fieldtype(BenchResult, field)), val)
            $constr_expr
        end
    end
end

const database_path = joinpath(@__DIR__, "..", "results", "data", "database.jld2")

get_database() = load(database_path)["database"]
update_database!(data::Vector{BenchResult}) = save(database_path, Dict("database" => data))

function append_data!(results::Vector{BenchResult})
    db = get_database()
    append!(db, results)
    update_database!(db)
    return
end


text:
#=
Copyright 2017 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

# Inlet and outlet boundary conditions are applied at the beginning and at the
# ending of each time step respectively
# (see [`solveModel`](godunov.html#solveModel)). Inlet (outlet) condition is
# applied by setting flow (pressure) value at the first (last) vessel node.
# [Inlet](boundary_conditions.html#inletCompatibility) and
# [outlet](boundary_conditions.html#outletCompatibility) compatibility
# relations are defined to compute the remaining quantities:
# `P`, `u`, `A`, `c`, and `Q`, `u`, `A`, `c` respectively.

# *function* __`setInletBC`__
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------- --------------------------------------------------------------
# `t`           `::Float` current time in seconds.
#
# `dt`          `::Float` $\Delta t$, current time step.
#
# `v`           `::Vessel` vessel data structure.
#
# `h`           `::Heart` heart data structure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# The flow (pressure) at the first node of the vessel (`Q[1]`) can be set by
# several
# functions. This function is selected by `BC_switch` from `h` data structure.
# The inlet compatibility relations are handled by
# [`inletCompatibility`](boundary_conditions.html#inletCompatibility).
# ----------------------------------------------------------------------------
#
function setInletBC(t :: Float64, dt :: Float64,
                    v :: Vessel,  h  :: Heart)

	if h.BC_switch == 1
		v.Q[1] = sinHeaviside(t, h)

	elseif h.BC_switch == 2
		v.Q[1] = gauss(t, h)

	elseif h.BC_switch == 3

  	if h.inlet_type == "Q"
			v.Q[1] = inputFromData(t, h)
  	else
			v.P[1] = inputFromData(t, h)
  	end

	end

	inletCompatibility(dt, v, h)
end

# `heart_data.BC_switch` = 1 sets a half sine flow waveform at the vessel
# inlet. Flow values at each time step are computed by the following function:

# *function* __`sinHeaviside`__ $\rightarrow$ `flow::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------- --------------------------------------------------------------
# `t`           `::Float` current time in seconds.
#
# `h`           `::Heart` heart data structure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# The flow waveform is always bigger than the initial flow rate assigned at
# the initialisation step, i.e. `initial_flow` ($Q_i$ in the picture). The
# sinusoid is multiplied by the desired maximum flow rate `flow_amplitude`
# ($A_Q$), and by the [`heaviside`](boundary_conditions.html#heaviside) step
# function. This is done in order to ensure only a positive inlet flow. We
# want also the waveform to be different than zero for the entire length of
# the systole, defined by `h.sys_T` variable. This is done by setting the
# `sin` argument as $\frac{2 \pi t}{T_s}$.
# In this way we can define a systolic impulse as narrow as wanted without
# worrying about the cardiac period which is completely unrelated.
# ----------------------------------------------------------------------------
#

#

#
# ----------------------------------------------------------------------------
# Returns:
# --------- ------------------------------------------------------------------
# `Q`       `::Float` flow rate value at current time.
# ----------------------------------------------------------------------------
#
function sinHeaviside(t :: Float64, h :: Heart)

	return h.flow_amplitude*(sin(2*pi*t/h.sys_T)*heaviside(-t+h.sys_T*0.5)) +
          h.initial_flow
end

# *function* __`heaviside`__ $\rightarrow$ `::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------- --------------------------------------------------------------
# `n`           `::Float`
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# $$
#   H(n) = \begin{cases}
#           0 \quad \text{if } n < 0, \\
#           1 \quad \text{if } n \geq 0.
#          \end{cases}
# $$
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Returns:
# ------------- --------------------------------------------------------------
# `H`           `::Float` either 1 or 0.
# ----------------------------------------------------------------------------
#
function heaviside(n :: Float64)

	if n < 0.

		return 0.

	else

		return 1.

	end

end

# `heart_data.BC_switch` = 2 sets a gaussian shaperd flow waveform at the
# vessel inlet. The bell shape has smoother transition at $T_s$ than
# [`sinHeaviside`](boundary_conditions.html#sinHeaviside). As a result, the
# numerical scheme produces less oscillations due to that discontinuity.
# Flow values at each time step are computed by the following function:

# *function* __`gauss`__ $\rightarrow$ `flow::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------- --------------------------------------------------------------
# `t`           `::Float` current time in seconds.
#
# `h`           `::Heart` heart data structure.
# ----------------------------------------------------------------------------
#
function gauss(t :: Float64, h :: Heart)

	# --------------------------------------------------------------------------
	# Functioning
	# --------------------------------------------------------------------------
	# `t_hat` variable tells us how many cardiac cycles have already been
	# simulated. This is calculated with `div(a, b)` which returns `a/b`,
	# truncating to an integer.
	# --------------------------------------------------------------------------
	t_hat = div(t,h.cardiac_T)
	# The gaussian bell should span the entire systolic period `h.sys_T` length.
	# For `gauss` to return a value different than zero, the current time `t`
	# must be within the current systolic period `h.sys_T` ($T_s$).
	if t < h.sys_T + h.cardiac_T*t_hat && t >= h.cardiac_T*t_hat
		# When this is confirmed, the inlet flow rate reads
		# $$
		#   Q = A_Q \exp{\left(- \frac{\left(t - \hat{t} T -
		#   \tfrac{T_s}{2}\right)^2}{2 \left(\tfrac{T_s}{8}\right)^2 }\right)}
		#   + Q_i.
		# $$
		# Variables are the same as those used in
		# [`sinHeaviside`](boundary_conditions.html#sinHeaviside).
		return h.flow_amplitude * exp( - (t-t_hat*h.cardiac_T -h.sys_T*0.5)^2 /
            (2*(h.sys_T*0.125)^2)) + initial_flow
	# Otherwise only the `initial_flow` is returned.
	#
	# --------------------------------------------------------------------------
	# Returns:
	# ---------- ---------------------------------------------------------------
	# `Q`        `::Float` inlet flow rate at current time.
	# --------------------------------------------------------------------------
	else
		return 0. + initial_flow

	end
end

# `heart_data.BC_switch` = 3 requires a flow waveform to be supplied by the
# user. The flow waveform is read by the following function:

# *function* __`inputFromData`__ $\rightarrow$ `flow::Float`
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------- --------------------------------------------------------------
# `t`           `::Float` current time in seconds.
#
# `h`           `::Heart` heart data structure.
# ----------------------------------------------------------------------------
#
function inputFromData(t :: Float64, h :: Heart)

	# --------------------------------------------------------------------------
	# Functioning
	# --------------------------------------------------------------------------
	# Inlet flow waveform is stored in `h` data structure inside `input_data`
	# matrix. `input_data` has two columns and as many rows as time steps. The
	# first column (`input_data[:,1]`) contains the time variable. The second
	# column (`input_data[:,2]`) contains the flow variable. The two arrays are
	# copied into two variables, `idt` and `idq` respectively.
	# --------------------------------------------------------------------------
	idt = h.input_data[:,1]
	idq = h.input_data[:,2]
	# `t_hat` variable tells us how many cardiac cycles have already been
	# simulated. This is calculated with `div(a, b)` which returns `a/b`,
	# truncating to an integer.
	t_hat = div(t,h.cardiac_T)
	# The current time is then updated by subtracting as many cardiac cycles as
	# those already passed. This is done in order to have a reference within the
	# cardiac cycle, i.e. in order to know at what point (*when*) in the cardiac
	# cycle the simulation is running.
	t -= t_hat*h.cardiac_T
	# The waveform is given in a discrete form and usually with a small amount
	# of points along the cardiac cycle. Usually, the numerical $\Delta t$ is
	# small enough to fall between two time steps reported in `idt`. This
	# loop skims all the entries in `idt` and finds the time interval containing
	# the current time `t`. The variable `idx` is initialised to zero for
	# simplicity.
	#
	#

	#

	idx = 0
	for i in 1:length(idt)
		if ((t >= idt[i]) && (t <= idt[i+1]))
			idx=i
			break
		end
	end
	# The flow `qu` is then computed at the current time by a linear
	# interpolation between the two known flows (`idq[idx]` and `idq[idx+1]`)
	# at the ends of the time interval.
	qu = idq[idx] + (t - idt[idx]) * (idq[idx+1] - idq[idx]) /
       (idt[idx+1] - idt[idx])
	# --------------------------------------------------------------------------
	# Returns:
	# ---------- ---------------------------------------------------------------
	# `qu`       `::Float` inlet flow rate at current time.
	# --------------------------------------------------------------------------
	return qu
end

# In the case of veins, the inlet is always coupled with three-element Windkessel.
# The inlet function is re-defined as follows

# *function* __`setInletBC`__
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------- --------------------------------------------------------------
# `v`           `::Vessel` vein data structure.
#
# `dt`          `::Float` $\Delta t$, current time step.
#
# `Pc`           `::Float` pressure across the peripheral compliance of the
#                three-element Windkessel.
#
# `R2`           `::Float` peripeheral resistance of the coupled three-element
#                Windkessel.
#
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# The flow at the first node of the vein (`Q[1]`) is set as the flow exiting
# the three-element Windkessel computed as
# $$
#   Q_{vein-in} = \frac{P_c}{R_2},
# $$
# and a small factor of `1e-10` is added to avoid numerical instabilities.
# Then, the inlet compatibility relations are computed.
# ----------------------------------------------------------------------------
function setInletBC(v  :: Vessel,  dt :: Float64,
                    Pc :: Float64, R2 :: Float64)
	v.Q[1] = Pc/R2 + 1e-10

	inletCompatibility(dt, v)
end

# The inlet compatibility relations are handled by
# [`inletCompatibility`](boundary_conditions.html#inletCompatibility).

# *function* __`inletCompatibility`__
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------- --------------------------------------------------------------
# `dt`          `::Float` $\Delta t$, current time step.
#
# `v`           `::Vessel` vessel data structure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# Compatibility relations are derived by using a technique called
# *extrapolation of characteristics* [^1].
# ----------------------------------------------------------------------------
#
function inletCompatibility(dt :: Float64, v :: Vessel, h :: Heart)

	# Riemann invariants are [computed](converter.html#riemannInvariants)
	# starting from variables at time `t`.
	# $$
	#   W_{11}^t = u - 4c, \quad W_{21}^t = u + 4c,
	# $$
	# are computed at the inlet node of the vessel.
	# `W12` and `W22` are the same quantities computed at the second node, i.e.
	# the closest node at the right hand side of the inlet node.
	W11, W21 = riemannInvariants(1, v)
	W12, W22 = riemannInvariants(2, v)

	# The two Riemann invariants within the two nodes are extrapolated with a
	# linear law. They reads
	# $$
	#   W_{11}^{t+1} = W_{11}^t + \left(W_{12}^t - W_{11}^t \right) (c_1 - u_1)
	#    \frac{\Delta t}{\Delta x},
	# $$
	# $$
	#   W_{21}^{t+1} = W_{21}^t - W_{11}^{t+1} + 2 \frac{Q_1}{A_1},
	# $$
	# where the 1 subscript indicates the inlet node.
	W11 += (W12-W11)*(v.c[1] - v.u[1])*dt/v.dx
	W21 = 2*v.Q[1]/v.A[1] - W11
	# Longitudinal velocity and wave velocity at the inlet node are retrieved
	# by [`rI2uc`](converter.html#rI2uc) (*Riemann invariants to u and c*)
	# function.
	v.u[1], v.c[1] = rI2uc(W11, W21)
	# `A` and `P` quantities follow from their definitions.

	if h.inlet_type == "Q"
		v.A[1] = v.Q[1]/v.u[1]
		v.P[1] = pressure(v.A[1], v.A0[1], v.beta[1], v.Pext)
	else
		v.A[1] = areaFromPressure(v.P[1], v.A0[1], v.beta[1], v.Pext)
		v.Q[1] = v.u[1]*v.A[1]
	end

end

function inletCompatibility(dt :: Float64, v :: Vessel)

	W11, W21 = riemannInvariants(1, v)
	W12, W22 = riemannInvariants(2, v)

	W11 += (W12-W11)*(v.c[1] - v.u[1])*dt/v.dx
	W21 = 2*v.Q[1]/v.A[1] - W11

	v.u[1], v.c[1] = rI2uc(W11, W21)

	v.A[1] = v.Q[1]/v.u[1]
	v.P[1] = pressure(v.A[1], v.A0[1], v.beta[1], v.Pext)

end

# Outlet boundary condition is applied to the last node of the involved
# vessel. Two boundary conditions are available. Either a reflection
# coefficient is provided or a three elements windkessel is coupled.

# *function* __`setOutletBC`__
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------- --------------------------------------------------------------
# `dt`          `::Float` $\Delta t$, current time step.
#
# `v`           `::Vessel` vessel data structure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# Outlet boundary conditions are selected in the `project.csv` file.
# ----------------------------------------------------------------------------
#
function setOutletBC(dt :: Float64, v :: Vessel)

	# A proximal resistance `R1` set to zero means that a reflection coefficient
	# was specified. Thus, the pressure `P` is first linearly extrapolated [^2]
	# as
	# $$
	#   P_M = 2 P_{M-1} - P_{M-2},
	# $$
	# where $M$ is the index of the vessel outlet node. Then, the compatibility
	# relations are applied by means of
	# [`outletCompatibility`](boundary_conditions.html#outletCompatibility)
	# function.
	if v.R1 == 0.

		v.P[end] = 2*v.P[end-1] - v.P[end-2]
		outletCompatibility(dt, v)
	# When all the parameters of the three element windkessel are specified,
	# [`wk3`](boundary_conditions.html#wk3) function is called. This function
	# solves lumped parameter equations and it also applies compatibility
	# relations.
	else
		wk3(dt, v)

	end

end

# Outlet compatibility relations compute all the quantities not directly
# assigned by the outlet boundary condition.

# *function* __`outletCompatibility`__
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------- --------------------------------------------------------------
# `dt`          `::Float` $\Delta t$, current time step.
#
# `v`           `::Vessel` vessel data structure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# Compatibility relations are derived by using a technique called
# *extrapolation of characteristics* as for the
# [inlet](boundary_conditions.html#inletCompatibility).
# ----------------------------------------------------------------------------
#
function outletCompatibility(dt :: Float64, v :: Vessel)

	W1M1, W2M1 = riemannInvariants(v.M-1, v)
	W1M, W2M   = riemannInvariants(v.M, v)

	# $$
	#   W_{2M}^{t+1} = W_2^t + \frac{W_{2 M-1}^t - W_{2M}^t}{\Delta x} (u_M^t
	#                   + c_M^t) \Delta t,
	# $$
	# $$
	#   W_{1M}^{t+1} = W_{1M}^0 - R_t \left(W_{2L}^{t+1} - W_{2L}^0 \right).
	# $$
	W2M += (W2M1 - W2M)*(v.u[end] + v.c[end])*dt/v.dx
	W1M = v.W1M0 - v.Rt * (W2M - v.W2M0)
	# Here the remaining quatities are computed by means of their definitions.
	v.u[end], v.c[end] = rI2uc(W1M, W2M)
	v.Q[end] = v.A[end]*v.u[end]

end

# The three element windkessel simulates the perfusion of downstream vessels.
# This 0D model is coupled by [`wk3`](boundary_conditions.html#wk3)
# function to 1D model terminal branches via the solution of a Riemann
# problem at the 0D/1D interface.
#
#

#

#
# $R_1$ is the proximal resistance, $R_2$ is the peripeheral reisistance,
# $C_c$ is the peripheral compliance, $P_c$ is the pressure across the
# peripheral compliance, $P_out$ is the pressure at the artery-vein interface,
# $P_e$ is the pressure at the 0D/1D interface.

# *function* __`wk3`__
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------- --------------------------------------------------------------
# `dt`          `::Float` $\Delta t$, current time step.
#
# `v`           `::Vessel` vessel data structure.
# ----------------------------------------------------------------------------
#
function wk3(dt :: Float64, v :: Vessel)

	# --------------------------------------------------------------------------
	# Functioning [^3]
	# --------------------------------------------------------------------------
	# At capillary level the pressure is assumed to be zero, i.e. $P_{out}=0.$
	# --------------------------------------------------------------------------
	Pout = 0.
	# The coupling is performed by assuming that an intermediate state
	# $(A^*,u^*)$ generates from $(A_l,u_l)$ (1D outlet) and $(A_r,u_r)$ (0D
	# inlet). This intermediate state must satisfy the windkessel equation
	# $$
	#   A^*u^* \left(1 + \frac{R_1}{R_2}\right) + C_c R_1 \frac{\partial
	#   (A^*u^*)}{\partial t} = \frac{P_e - P_{out}}{R_2} + C_c \frac{\partial
	#   P_e}{\partial t}.
	# $$
	Al = v.A[end]
	ul = v.u[end]
	# $P_c$ is computed at each time step from
	# $$
	#   C_c \frac{\partial P_c}{\partial t} = A^*u^* - \frac{P_c -
	#   P_{out}}{R_2},
	# $$
	# which is discretised numerically with a first-order scheme. $P_c$ is
	# is [initialised](initialise.html#initialiseVessel) to zero.
	v.Pc += dt/v.Cc * (Al*ul - (v.Pc-Pout)/v.R2)
	# @printf "%8.3e\n" v.Pc/v.R2*1e6
	# We consider $\beta$ and $A_0$ to be the same on both sides of the 0D/1D
	# interface. It yields the nonlinear equation
	# $$
	#   \mathcal{F}(A^*) = A^*R_1\left(u_l+4c_l\right)
	#   -4A^*R_1c^* - \frac{\beta}{A_0}\left(
	#   \sqrt{A^*}-\sqrt{A_0}\right) + P_c = 0,
	# $$
	# where $c_l$ and $c^*$ are the wave speeds calculated with $A_l$ and $A^*$,
	# respectively.
	As = Al

	# fun(As) = v.R1*(ul+4*sqrt(3*v.gamma[end]*sqrt(Al)*0.5))*As -
	# 	4*v.R1*sqrt(3*v.gamma[end]*sqrt(As)*0.5)*As -
	# 	(v.Pext + v.beta[end]*(sqrt(As/v.A0[end]) - 1)) + v.Pc

	# dfun(As) = v.R1*( ul +4*sqrt(1.5*v.gamma[end])*(sqrt(sqrt(Al)) - 1.25*sqrt(sqrt(As))) ) -
	# 	v.beta[end]*0.5/(sqrt(As*v.A0[end]))

	ssAl = sqrt(sqrt(Al))
	sgamma =  2*sqrt(6*v.gamma[end])
	sA0 = sqrt(v.A0[end])
	bA0 = v.beta[end]/sA0

	fun(As) = As*v.R1*(ul + sgamma * (ssAl - sqrt(sqrt(As)))) -
		(v.Pext + bA0*(sqrt(As) - sA0)) + v.Pc

	dfun(As) = v.R1*(ul + sgamma * (ssAl - 1.25*sqrt(sqrt(As)))) - bA0*0.5/sqrt(As)

	# $\mathcal{F}(A^*)=0$ is solved for $A^*$ with the Newton's method
	# implemented in [__Roots__](https://github.com/JuliaLang/Roots.jl) library.
	# $A^*$ is initialised equal to $A_l$.
	# As = newton(fun, dfun, As)
	As = newtone(fun, dfun, As)
	# Once $A^*$ is found, $u^*$ reads
	# $$
	#   u^* = \frac{P_e^* - P{out}}{A^*R_1},
	# $$
	# where $P_e^*$ is $P_e$ calculated with $A^*$.
	us = (pressure(As, v.A0[end], v.beta[end], v.Pext) - Pout)/(As*v.R1)

	v.A[end] = As
	v.u[end] = us

end

function newtone(f, df, x0)

	xn = x0 - f(x0)/df(x0)
	if abs(xn-x0)<= 1e-5
		return xn
	else
		newtone(f, df, xn)
	end
end

# *function* __`updateGhostCells`__
#
# ----------------------------------------------------------------------------
# Parameters:
# ------------- --------------------------------------------------------------
# `v`           `::Vessel` vessel data structure.
# ----------------------------------------------------------------------------
#
# ----------------------------------------------------------------------------
# Functioning
# ----------------------------------------------------------------------------
# Ghost cells are updated by assuming both ends as transmissive boundaries.
# ----------------------------------------------------------------------------
#
function updateGhostCells(v :: Vessel)

	v.U00A = v.A[1]
	v.U00Q = v.Q[1]
	v.U01A = v.A[2]
	v.U01Q = v.Q[2]

	v.UM1A = v.A[v.M]
	v.UM1Q = v.Q[v.M]
	v.UM2A = v.A[v.M-1]
	v.UM2Q = v.Q[v.M-1]

end

# `updateGhostCells` exploits julia's multiple dispatch property. This
# function can handle one or more than one vessels organised in an `Array`
# structure.
function updateGhostCells(vessels :: Array{Vessel, 1})

	for vessel in vessels
		updateGhostCells(vessel)
	end
#   if length(vessels) == 1
#     updateGhostCells(vessels[1])
#   else
#     pmap(updateGhostCells, vessels)
#   end
end

# ### References
#
# [^1]: Peiro, Joaquim, and Alessandro Veneziani.
# [*"Reduced models of the cardiovascular system
# ."*](http://link.springer.com/chapter/10.1007/978-88-470-1152-6_10#page-1)
# In Cardiovascular mathematics, pp. 347-394. Springer Milan, 2009.
#
# [^2]: Anderson, John David, and J. Wendt. *Computational fluid dynamics*.
# Vol. 206. New York: McGraw-Hill, 1995.
#
# [^3]: Fernandez, Miguel Angel, Vuk Milisic, and Alfio Quarteroni.
# [*"Analysis of a geometrical multiscale blood flow model based on
# the coupling of ODEs and hyperbolic
# PDEs."*](http://epubs.siam.org/doi/abs/10.1137/030602010)
# Multiscale Modeling & Simulation 4.1 (2005): 215-236.


text:
const TimeNS = typeof(time_ns())

struct ThreadLocalAccessRecord
    ptr::UInt
    time::TimeNS
    tag::UInt
end

const ACCESSES = Vector{ThreadLocalAccessRecord}[]
const LENGTHS = Int[]

function MemoryAccesses.init(n = 2^10)
    if length(ACCESSES) != Threads.nthreads() || length(LENGTHS) != 8 * Threads.nthreads()
        append!(empty!(ACCESSES), (Union{}[] for _ in 1:Threads.nthreads()))
        resize!(LENGTHS, 8 * Threads.nthreads())
    end
    fill!(LENGTHS, 0)
    for a in ACCESSES
        resize!(a, n)
        fill!(a, ThreadLocalAccessRecord(0, 0, 0))
    end
    init_reference()
    return
end

function MemoryAccesses.clear()
    fill!(LENGTHS, 0)
    for a in ACCESSES
        fill!(a, ThreadLocalAccessRecord(0, 0, 0))
    end
    clear_reference()
    return
end

"""
    MemoryAccesses.record(x, [tag = UInt(0)])

Record that `x` is accessed.

* If `x` is an array, `pointer(x, 1)` is recorded.
  (TODO: record the first and the last pointers?)
* Otherwise `x` is assumed to be convertable to an `UInt` (e.g., `Ptr`) whose
  value is recorded.

Optionally, arbitrary value `tag` that is convertable to an `UInt` can be passed
as the second argument.
"""
MemoryAccesses.record

MemoryAccesses.record(xs::AbstractArray, tag = UInt(0)) =
    MemoryAccesses.record(pointer(xs, 1), tag)

MemoryAccesses.record(ptr, tag = UInt(0)) = MemoryAccesses.record(UInt(ptr), tag)
function MemoryAccesses.record(ptr::UInt, tag = UInt(0))
    tid = Threads.threadid()
    i = LENGTHS[8*tid] += 1
    ACCESSES[tid][min(i, end)] = ThreadLocalAccessRecord(ptr, time_ns(), UInt(tag))
    return
end

Base.@propagate_inbounds accesslength(tid) = LENGTHS[8*tid]

MemoryAccesses.isfull() =
    any(length(a) < accesslength(tid) for (tid, a) in enumerate(ACCESSES))

function checkfull()
    MemoryAccesses.isfull() || return
    n = length(ACCESSES[1])
    @warn """
    Buffer is full. Call `MemoryAccesses.clear(n)` to increase the buffer size.
    Currently, `n = $n`.
    """
end

function MemoryAccesses.view()
    checkfull()
    accesses = [@view(a[1:min(end, accesslength(tid))]) for (tid, a) in enumerate(ACCESSES)]
    references = Iterators.flatten(REFERENCES)
    return (; accesses, references)
end

function MemoryAccesses.copy()
    v = MemoryAccesses.view()
    return (; accesses = map(copy, v.accesses), references = collect(v.references))
end


text:
export MADDPGManager

"""
    MADDPGManager(; agents::Dict{<:Any, <:Agent}, args...)
Multi-agent Deep Deterministic Policy Gradient(MADDPG) implemented in Julia. Here only works for simultaneous games whose action space is discrete.
See the paper https://arxiv.org/abs/1706.02275 for more details.

# Keyword arguments
- `agents::Dict{<:Any, <:NamedPolicy{<:Agent{<:DDPGPolicy, <:AbstractTrajectory}, <:Any}}`, here each agent collects its own information. While updating the policy, each `critic` will assemble all agents' trajectory to update its own network.
- `batch_size::Int`
- `update_freq::Int`
- `update_step::Int`, count the step.
- `rng::AbstractRNG`.
"""
mutable struct MADDPGManager{P<:DDPGPolicy, T<:AbstractTrajectory, N<:Any} <: AbstractPolicy
    agents::Dict{<:N, <:Agent{<:NamedPolicy{<:P, <:N}, <:T}}
    batch_size::Int
    update_freq::Int
    update_step::Int
    rng::AbstractRNG
end

# for simultaneous game with a discrete action space.
function (π::MADDPGManager)(env::AbstractEnv)
    while current_player(env) == chance_player(env)
        env |> legal_action_space |> rand |> env
    end
    Dict((player, ceil(agent.policy(env))) for (player, agent) in π.agents)
end

function (π::MADDPGManager)(stage::Union{PreEpisodeStage, PostActStage}, env::AbstractEnv)
    # only need to update trajectory.
    for (_, agent) in π.agents
        update!(agent.trajectory, agent.policy, env, stage)
    end
end

function (π::MADDPGManager)(stage::PreActStage, env::AbstractEnv, actions)
    # update each agent's trajectory.
    for (player, agent) in π.agents
        update!(agent.trajectory, agent.policy, env, stage, actions[player])
    end

    # update policy
    update!(π)
end

function (π::MADDPGManager)(stage::PostEpisodeStage, env::AbstractEnv)
    # collect state and a dummy action to each agent's trajectory here.
    for (_, agent) in π.agents
        update!(agent.trajectory, agent.policy, env, stage)
    end

    # update policy
    update!(π)
end

# update policy
function RLBase.update!(π::MADDPGManager)
    π.update_step += 1
    π.update_step % π.update_freq == 0 || return

    for (_, agent) in π.agents
        length(agent.trajectory) > agent.policy.policy.update_after || return
        length(agent.trajectory) > π.batch_size || return
    end

    # get training data
    temp_player = collect(keys(π.agents))[1]
    t = π.agents[temp_player].trajectory
    inds = rand(π.rng, 1:length(t), π.batch_size)
    batches = Dict((player, RLCore.fetch!(BatchSampler{SARTS}(π.batch_size), agent.trajectory, inds))
                for (player, agent) in π.agents)

    # get s, a, s′ for critic
    s = Flux.stack((batches[player][:state] for (player, _) in π.agents), 1)
    a = Flux.stack((batches[player][:action] for (player, _) in π.agents), 1)
    s′ = Flux.stack((batches[player][:next_state] for (player, _) in π.agents), 1)

    # for training behavior_actor
    mu_actions = Flux.stack(
        ((
            batches[player][:state] |> # get personal state information
            x -> send_to_device(device(agent.policy.policy.behavior_actor), x) |>
            agent.policy.policy.behavior_actor |> send_to_host
        ) for (player, agent) in π.agents), 1
    )
    # for training behavior_critic
    new_actions = Flux.stack(
        ((
            batches[player][:next_state] |> # get personal next_state information
            x -> send_to_device(device(agent.policy.policy.target_actor), x) |>
            agent.policy.policy.target_actor |> send_to_host
        ) for (player, agent) in π.agents), 1
    )

    for (player, agent) in π.agents
        p = agent.policy.policy # get DDPGPolicy struct
        A = p.behavior_actor
        C = p.behavior_critic
        Aₜ = p.target_actor
        Cₜ = p.target_critic

        γ = p.γ
        ρ = p.ρ

        _device(x) = send_to_device(device(A), x)

        # Note that here default A, C, Aₜ, Cₜ on the same device.
        s, a, s′ = _device((s, a, s′))
        mu_actions = _device(mu_actions)
        new_actions = _device(new_actions)
        r = _device(batches[player][:reward])
        t = _device(batches[player][:terminal])

        qₜ = Cₜ(vcat(s′, new_actions)) |> vec
        y = r .+ γ .* (1 .- t) .* qₜ

        gs1 = gradient(Flux.params(C)) do
            q = C(vcat(s, a)) |> vec
            loss = mean((y .- q) .^ 2)
            ignore() do
                p.critic_loss = loss
            end
            loss
        end

        update!(C, gs1)

        gs2 = gradient(Flux.params(A)) do
            loss = -mean(C(vcat(s, mu_actions)))
            ignore() do
                p.actor_loss = loss
            end
            loss
        end

        update!(A, gs2)

        # polyak averaging
        for (dest, src) in zip(Flux.params([Aₜ, Cₜ]), Flux.params([A, C]))
            dest .= ρ .* dest .+ (1 - ρ) .* src
        end
    end
end


text:
@testset "exp problems" begin
    @testset "exp1_test" begin
        #=
        primal
            min x + y + z
        s.t.
            y e^(x/y) <= z, y > 0 (i.e (x, y, z) are in the exponential primal cone) :w3, w_4, w_5
            x == 1 :w_1
            y == 2 :w_2

        dual
            max 2w_2 + w_1
        s.t.
            w_1 + w_3 == 1
            w_2 + w_4 == 1
            w_5 == 1
            (w_3, w_4, w_5) ∈ DualExponentialCone
        =#
        primal_model = exp1_test()
        dual_model, primal_dual_map = dual_model_and_map(primal_model)

        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 5
        list_of_cons = MOI.get(dual_model, MOI.ListOfConstraints())
        @test Set(list_of_cons) == Set(
            [
                (SAF{Float64}, MOI.EqualTo{Float64})
                (VVF, MOI.DualExponentialCone)
            ],
        )
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VVF,MOI.DualExponentialCone}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SAF{Float64},MOI.EqualTo{Float64}}(),
        ) == 3
        obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
        @test obj_type == SAF{Float64}
        obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
        @test MOI.get(dual_model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
        @test MOI.constant(obj) == 0.0
        @test MOI.coefficient.(obj.terms) == [2.0; 1.0]

        eq_con1_fun = MOI.get(
            dual_model,
            MOI.ConstraintFunction(),
            CI{SAF{Float64},MOI.EqualTo{Float64}}(2),
        )
        eq_con1_set = MOI.get(
            dual_model,
            MOI.ConstraintSet(),
            CI{SAF{Float64},MOI.EqualTo{Float64}}(2),
        )
        @test MOI.coefficient.(eq_con1_fun.terms) == [1.0; 1.0]
        @test MOI.constant.(eq_con1_fun) == 0.0
        @test MOI.constant(eq_con1_set) == 1.0
        eq_con2_fun = MOI.get(
            dual_model,
            MOI.ConstraintFunction(),
            CI{SAF{Float64},MOI.EqualTo{Float64}}(3),
        )
        eq_con2_set = MOI.get(
            dual_model,
            MOI.ConstraintSet(),
            CI{SAF{Float64},MOI.EqualTo{Float64}}(3),
        )
        @test MOI.coefficient.(eq_con2_fun.terms) == [1.0; 1.0]
        @test MOI.constant.(eq_con2_fun) == 0.0
        @test MOI.constant(eq_con2_set) == 1.0
        eq_con3_fun = MOI.get(
            dual_model,
            MOI.ConstraintFunction(),
            CI{SAF{Float64},MOI.EqualTo{Float64}}(4),
        )
        eq_con3_set = MOI.get(
            dual_model,
            MOI.ConstraintSet(),
            CI{SAF{Float64},MOI.EqualTo{Float64}}(4),
        )
        @test MOI.coefficient.(eq_con3_fun.terms) == [1.0]
        @test MOI.constant.(eq_con3_fun) == 0.0
        @test MOI.constant(eq_con3_set) == 1.0

        dual_exp_con = MOI.get(
            dual_model,
            MOI.ConstraintFunction(),
            CI{VVF,MOI.DualExponentialCone}(1),
        )
        @test dual_exp_con.variables == VI.(3:5)

        primal_con_dual_var = primal_dual_map.primal_con_dual_var
        @test primal_con_dual_var[CI{SAF{Float64},MOI.EqualTo{Float64}}(2)] ==
              [VI(1)]
        @test primal_con_dual_var[CI{SAF{Float64},MOI.EqualTo{Float64}}(3)] ==
              [VI(2)]
        @test primal_con_dual_var[CI{VVF,MOI.ExponentialCone}(1)] == VI.(3:5)

        primal_var_dual_con = primal_dual_map.primal_var_dual_con
        @test primal_var_dual_con[VI(1)] ==
              CI{SAF{Float64},MOI.EqualTo{Float64}}(2)
        @test primal_var_dual_con[VI(2)] ==
              CI{SAF{Float64},MOI.EqualTo{Float64}}(3)
        @test primal_var_dual_con[VI(3)] ==
              CI{SAF{Float64},MOI.EqualTo{Float64}}(4)
    end

    @testset "exp2_test" begin
        #=
        primal
            min x + y + z
        s.t.
            y e^(x/y) <= z, y > 0 (i.e (x, y, z) are in the exponential primal cone) :w3, w_4, w_5
            x == 1 :w_1
            y == 2 :w_2

        dual
            max 2w_2 + w_1
        s.t.
            w_1 + w_3 == 1
            w_2 + w_4 == 1
            w_5 == 1
            (w_3, w_4, w_5) ∈ DualExponentialCone
        =#
        primal_model = exp2_test()
        dual_model, primal_dual_map = dual_model_and_map(primal_model)

        @test MOI.get(dual_model, MOI.NumberOfVariables()) == 5
        list_of_cons = MOI.get(dual_model, MOI.ListOfConstraints())
        @test Set(list_of_cons) == Set(
            [
                (SAF{Float64}, MOI.EqualTo{Float64})
                (VVF, MOI.DualExponentialCone)
            ],
        )
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{VVF,MOI.DualExponentialCone}(),
        ) == 1
        @test MOI.get(
            dual_model,
            MOI.NumberOfConstraints{SAF{Float64},MOI.EqualTo{Float64}}(),
        ) == 3
        obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
        @test obj_type == SAF{Float64}
        obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
        @test MOI.get(dual_model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
        @test MOI.constant(obj) == 0.0
        @test MOI.coefficient.(obj.terms) == [2.0; 1.0]

        eq_con1_fun = MOI.get(
            dual_model,
            MOI.ConstraintFunction(),
            CI{SAF{Float64},MOI.EqualTo{Float64}}(2),
        )
        eq_con1_set = MOI.get(
            dual_model,
            MOI.ConstraintSet(),
            CI{SAF{Float64},MOI.EqualTo{Float64}}(2),
        )
        @test MOI.coefficient.(eq_con1_fun.terms) == [1.0; 1.0]
        @test MOI.constant.(eq_con1_fun) == 0.0
        @test MOI.constant(eq_con1_set) == 1.0
        eq_con2_fun = MOI.get(
            dual_model,
            MOI.ConstraintFunction(),
            CI{SAF{Float64},MOI.EqualTo{Float64}}(3),
        )
        eq_con2_set = MOI.get(
            dual_model,
            MOI.ConstraintSet(),
            CI{SAF{Float64},MOI.EqualTo{Float64}}(3),
        )
        @test MOI.coefficient.(eq_con2_fun.terms) == [1.0; 1.0]
        @test MOI.constant.(eq_con2_fun) == 0.0
        @test MOI.constant(eq_con2_set) == 1.0
        eq_con3_fun = MOI.get(
            dual_model,
            MOI.ConstraintFunction(),
            CI{SAF{Float64},MOI.EqualTo{Float64}}(4),
        )
        eq_con3_set = MOI.get(
            dual_model,
            MOI.ConstraintSet(),
            CI{SAF{Float64},MOI.EqualTo{Float64}}(4),
        )
        @test MOI.coefficient.(eq_con3_fun.terms) == [1.0]
        @test MOI.constant.(eq_con3_fun) == 0.0
        @test MOI.constant(eq_con3_set) == 1.0

        dual_exp_con = MOI.get(
            dual_model,
            MOI.ConstraintFunction(),
            CI{VVF,MOI.DualExponentialCone}(1),
        )
        @test dual_exp_con.variables == VI.(3:5)

        primal_con_dual_var = primal_dual_map.primal_con_dual_var
        @test primal_con_dual_var[CI{SAF{Float64},MOI.EqualTo{Float64}}(2)] ==
              [VI(1)]
        @test primal_con_dual_var[CI{SAF{Float64},MOI.EqualTo{Float64}}(3)] ==
              [VI(2)]
        @test primal_con_dual_var[CI{VAF{Float64},MOI.ExponentialCone}(1)] ==
              VI.(3:5)

        primal_var_dual_con = primal_dual_map.primal_var_dual_con
        @test primal_var_dual_con[VI(1)] ==
              CI{SAF{Float64},MOI.EqualTo{Float64}}(2)
        @test primal_var_dual_con[VI(2)] ==
              CI{SAF{Float64},MOI.EqualTo{Float64}}(3)
        @test primal_var_dual_con[VI(3)] ==
              CI{SAF{Float64},MOI.EqualTo{Float64}}(4)
    end
end


text:
# we are making an Octave model if this gets called, so lets load the Octave strategy -
include("./strategy/OctaveStrategy.jl")

function make_octave_model(path_to_model_file::String, path_to_output_dir::String; model_type::Symbol = :static)
    error("Ooops! Octave model generation has not yet been implemented.")
end

max_stars_repo_name: varnerlab/CFMG

max_issues_repo_name: varnerlab/CFMG

max_forks_repo_name: varnerlab/CFMG


text:
module Benchmark
using Lasso, GLMNet, BenchmarkLite, Distributions

# As in Friedman et al., generate random X and Y matrix such that SNR=3
function makeXY(ρ, nsamples, nfeatures)
    Σ = fill(ρ, nfeatures, nfeatures)
    Σ[diagind(Σ)] = 1
    X = rand(MvNormal(Σ), nsamples)'
    β = [(-1)^j*exp(-2*(j-1)/20) for j = 1:nfeatures]
    signal = (β'*Σ*β)[1]
    y = X*β + scale!(randn(nsamples), sqrt(signal/3))
    # @show uvar = var(y - X*(X\y))
    # @show evar = var(y - mean(y)) - uvar
    # @show evar/uvar
    (X, y)
end

type GLMNetOp{Dist,Naive} end
calc{Dist,Naive}(::GLMNetOp{Dist,Naive}, X, y) = glmnet(X, y, Dist(), naivealgorithm=Naive)
calc(::GLMNetOp{Binomial}, X, y) = glmnet(X, y, Binomial())

type LassoOp{Dist,Naive} end
calc{Dist,Naive}(::LassoOp{Dist,Naive}, X, y) = fit(LassoPath, X, y, Dist(), naivealgorithm=Naive, criterion=:coef)

type LassoBenchmark{Op} <: Proc end
Base.length(p::LassoBenchmark, n) = 0
Base.isvalid(p::LassoBenchmark, n) = true
Base.start(p::LassoBenchmark, n) = (gc(); inputs[n])
function Base.start{Naive}(p::LassoBenchmark{GLMNetOp{Binomial,Naive}}, n)
    X, y = inputs[n]
    yp = zeros(length(y), 2)
    yp[:, 1] = y .== 0
    yp[:, 2] = y .== 1
    gc()
    (X, yp)
end
Base.run{Op}(p::LassoBenchmark{Op}, n, s) =
    calc(Op(), s[1], s[2])
Base.done(p::LassoBenchmark, n, s) = true

Base.string{Dist,Naive}(p::LassoBenchmark{GLMNetOp{Dist,Naive}}) = "glmnet $(Dist.name.name) $(Naive ? "naive" : "cov")"
Base.string{Dist,Naive}(p::LassoBenchmark{LassoOp{Dist,Naive}}) = "Lasso.jl $(Dist.name.name) $(Naive ? "naive" : "cov")"

inputs = Dict{(Float64, Int, Int),(Matrix{Float64},Vector{Float64})}()
cfgs = vec([begin
                x = (ρ, N, p)
                inputs[x] = makeXY(ρ, N, p)
                x
            end for ρ in [0, 0.1, 0.2, 0.5, 0.9, 0.95],
               (N, p) in [(1000, 100), (5000, 100), (100, 200), (100, 500), (100, 1000), (100, 5000), #=(100, 20000)=#]])
rtable = run(Proc[LassoBenchmark{GLMNetOp{Normal,true}}(), LassoBenchmark{LassoOp{Normal,true}}(),
                  LassoBenchmark{GLMNetOp{Normal,false}}(), LassoBenchmark{LassoOp{Normal,false}}()], cfgs)
show(rtable; unit=:sec)

for (key, value) in inputs
    (X, y) = inputs[key]
    inputs[key] = (X, float64(rand(length(y)) .< 1./(1+exp(-y))))
end
rtable = run(Proc[LassoBenchmark{GLMNetOp{Binomial,true}}(), LassoBenchmark{LassoOp{Binomial,true}}()], cfgs)
show(rtable; unit=:sec)
end


text:
"""
    RaibertHopper <: Environment

    hopping robot, inspired by: "Dynamically Stable Legged Locomotion"
"""
struct RaibertHopper end

function raiberthopper(;
    representation=:minimal,
    timestep=0.05,
    gravity=[0.0; 0.0; -9.81],
    spring=0.0,
    damper=0.1,
    control_scaling=Diagonal(ones(3)),
    seed=1,
    contact_foot=true,
    contact_body=true,
    vis=Visualizer(),
    name=:robot,
    info=nothing,
    opts_step=SolverOptions(rtol=3.0e-4, btol=3.0e-4, undercut=1.5),
    opts_grad=SolverOptions(rtol=3.0e-4, btol=3.0e-4, undercut=1.5),
    T=Float64)

    mechanism = get_raiberthopper(;
        timestep,
        gravity,
        spring,
        damper,
        contact_foot,
        contact_body)
    initialize_raiberthopper!(mechanism)

    if representation == :minimal
        nx = minimal_dimension(mechanism)
    elseif representation == :maximal
        nx = maximal_dimension(mechanism)
    end

    nu = 3
    no = nx

    aspace = BoxSpace(nu,
        low=(-1.0e8 * ones(nu)),
        high=(1.0e8 * ones(nu)))
    ospace = BoxSpace(no,
        low=(-Inf * ones(no)),
        high=(Inf * ones(no)))

    rng = MersenneTwister(seed)

    z = get_maximal_state(mechanism)
    x = representation == :minimal ? maximal_to_minimal(mechanism, z) : z

    fx = zeros(nx, nx)
    fu = zeros(nx, nu)

    u_prev = zeros(nu)
    control_mask = raiberthopper_control_mask()
    build_robot(mechanism,
        vis=vis,
        name=name)

    TYPES = [RaibertHopper, T, typeof(mechanism), typeof(aspace), typeof(ospace), typeof(info)]
    env = Environment{TYPES...}(mechanism, representation, aspace, ospace,
        x, fx, fu,
        u_prev, control_mask' * control_scaling,
        nx, nu, no,
        info,
        [rng], vis,
        opts_step, opts_grad)

    return env
end

function raiberthopper_control_mask()
    [0 0 0 1 0 0 0;
	 0 0 0 0 1 0 0;
	 0 0 0 0 0 0 1]
end

function raiberthopper_nominal_max(;
    foot_radius=0.05)

    # initial state
    x1b1 = [0.0; 0.0; 0.5 + foot_radius]
    v1b1 = [0.0; 0.0; 0.0]
    q1b1 = [1.0; 0.0; 0.0; 0.0]
    ω1b1 = [0.0; 0.0; 0.0]
    z1b1 = [x1b1; v1b1; q1b1; ω1b1]

    x1b2 = [0.0; 0.0; 0.0 + foot_radius]
    v1b2 = [0.0; 0.0; 0.0]
    q1b2 = [1.0; 0.0; 0.0; 0.0]
    ω1b2 = [0.0; 0.0; 0.0]
    z1b2 = [x1b2; v1b2; q1b2; ω1b2]

    z1 = [z1b1; z1b2]
end

function raiberthopper_offset_max(x_shift, y_shift, z_shift)
    z = raiberthopper_nominal_max()
    shift = [x_shift; y_shift; z_shift]
    z[1:3] += shift
    z[13 .+ (1:3)] += shift
    return z
end

function visualize(env::Environment{RaibertHopper}, traj::Vector{Vector{T}};
    name=:robot,
    axes=false,
    grid=true) where T

	@assert size(traj[1]) == size(env.state)

    # convert to maximal representation
    z = [env.representation == :minimal ? minimal_to_maximal(env.mechanism, x) : x for x in traj]

    body_color = magenta
    foot_color = magenta
    env.mechanism.bodies[1].shape.color = body_color
    env.mechanism.bodies[2].shape.color = foot_color

    # build system
    build_robot(env.mechanism,
        vis=env.vis,
        name=name)

    n_leg = 100
    r_leg = 0.025
    for i = 1:n_leg
        setobject!(env.vis["leg$i"], GeometryBasics.Sphere(Point3f0(0),
            convert(Float32, r_leg)),
            MeshPhongMaterial(color=magenta))
    end

    # animate
    anim = MeshCat.Animation(convert(Int, floor(1.0 / env.mechanism.timestep)))
    for (t, x) in enumerate(z)
        x_body = x[1:3]
        x_foot = x[13 .+ (1:3)]
        dir = x_body - x_foot
        dir_norm = dir ./ norm(dir)

        MeshCat.atframe(anim, t) do
            set_robot(env.vis, env.mechanism, x,
                name=name)
            step = range(0.0, stop=norm(dir), length=n_leg)
            for i = 1:n_leg
                MeshCat.settransform!(env.vis["leg$i"], MeshCat.compose(MeshCat.Translation(step[i] .* dir_norm + x_foot), MeshCat.LinearMap(Dojo.RotY(0.0))))
            end
        end
    end
    MeshCat.setanimation!(env.vis, anim)

    set_camera!(env.vis, zoom=1.5)
    setvisible!(env.vis["/Axes"], axes)
    setvisible!(env.vis["/Grid"], grid)
end

function ghost(env::Environment{RaibertHopper}, traj::Vector{Vector{T}};
    timesteps=[t for t = 1:length(traj)],
    axes=false,
    grid=true,
    line=false) where T

    # convert to maximal representation
    z = [env.representation == :minimal ? minimal_to_maximal(env.mechanism, x) : x for x in traj]

    # color
    body_color = magenta
    foot_color = magenta
    env.mechanism.bodies[1].shape.color = body_color
    env.mechanism.bodies[2].shape.color = foot_color

    for t in timesteps
        # build system
        build_robot(env.vis, env.mechanism,
            name=Symbol("robot_$t"),
            color=(t == length(z) ? magenta : magenta_light))
        set_robot(env.vis, env.mechanism, z[t],
            name=Symbol("robot_$t"))

        x_body = z[t][1:3]
        x_foot = z[t][13 .+ (1:3)]

        leg = GeometryBasics.Cylinder(Point3f0(x_foot...), Point3f0(x_body...), convert(Float32, 0.025))
        setobject!(env.vis["leg_$t"], leg, MeshPhongMaterial(color=(t == length(z) ? magenta : magenta_light)))
    end
    setvisible!(env.vis[:robot], false)

    # body
    if line
        line_mat = LineBasicMaterial(color=color=RGBA(51.0 / 255.0, 1.0, 1.0, 1.0), linewidth=25.0)
        points = Vector{Point{3,Float64}}()
        for (i, xt) in enumerate(z)
            k = xt[1:3]
            push!(points, Point(k[1], k[2], k[3]))
        end
        setobject!(env.vis[:body_traj], MeshCat.Line(points, line_mat))

        # foot
        line_mat = LineBasicMaterial(color=color=RGBA(1.0, 153.0 / 255.0, 51.0 / 255.0, 1.0), linewidth=10.0)
        points = Vector{Point{3,Float64}}()
        for (i, xt) in enumerate(z)
            k = xt[13 .+ (1:3)]
            push!(points, Point(k[1], k[2], k[3]))
        end
        setobject!(env.vis[:foot_traj], MeshCat.Line(points, line_mat))
    end

    set_camera!(env.vis,
        zoom=1.5)
    setvisible!(env.vis["/Axes"], axes)
    setvisible!(env.vis["/Grid"], grid)

    open(env.vis)
end


text:
module TestBaselineMethods

include(joinpath("..", "fixture_learners.jl"))
using .FixtureLearners
fcp = FeatureClassification()
nfcp = NumericFeatureClassification()

using FactCheck

importall Orchestra.Transformers.BaselineMethods

facts("Baseline transformers") do
  context("Baseline learner does simple transforms") do
    bl = Baseline()
    instances = [
      1 1;
      1 2;
      2 2;
    ]
    labels = ["a"; "a"; "b"]
    fit!(bl, instances, labels)
    transformed = transform!(bl, instances)
    expected_transformed = ["a"; "a"; "a"]

    @fact transformed => expected_transformed
  end
  context("Identity returns instances as is") do
    id = Identity()
    fit!(id, nfcp.train_instances, nfcp.train_labels)
    transformed = transform!(id, nfcp.test_instances)
    expected_transformed = nfcp.test_instances

    @fact transformed => expected_transformed
  end
end

end # module


text:
# Synaptic kinetics
syn∞ = 1/(1 + exp((-35 - Vₘ)/5))
glut_kinetics = Gate(SteadyStateTau, :s,
                     syn∞,
                     (1 - syn∞)/(1/40), 1)

chol_kinetics = Gate(SteadyStateTau, :s,
                     syn∞,
                     (1 - syn∞)/(1/100), 1)

EGlut = Equilibrium{Leak}(-70mV, :Glut) # NOTE: -70mV reversal => IPSP
EChol = Equilibrium{Leak}(-80mV, :Chol) # Leak as alias for non-specific ion current

@named Glut = SynapticChannel(Leak, [glut_kinetics], EGlut);
@named Chol = SynapticChannel(Leak, [chol_kinetics], EChol);


text:
module TumorGrowth

print("Loading Packages... ")
import DataFrames: DataFrame
using DataFrames: DataFrameRow, nrow
using CSV
using Plots: palette, distinguishable_colors
using GLMakie, Makie
using LinearAlgebra
using Statistics
using StaticArrays
using StatsBase: countmap
using Random
using Distributions
using ProgressMeter
println("Done!")

(
    "simulation.jl",
    "plotting.jl",
    "pushing.jl",
    "pushing_animate.jl",
    "sampling.jl",
    "analysis.jl",
    "time_series.jl",
    "nonspatial.jl"
) .|> include

export data_import

"""
    data_import(path::String; delim="\t") :: DataFrame

Convenience function which loads a DataFrame from a '.csv' file and correctly
parses columns with elementtype 'Vector'.
"""
function data_import(path::String; delim="\t")
    data = DataFrame(CSV.File(path; delim=delim); copycols=false)
    arraycols = collect(first(data)) .|> val ->
        (val isa String || val isa InlineString) && all(occursin.(("[","]"), val))
    if any(arraycols)
        data[!,arraycols] .= data[!,arraycols] .|> Meta.parse .|> eval
    end
    for col in findall(.!(arraycols))
            data[!, col] .= data[!, col] |> Vector
    end
    return data
end

end  # module TumorGrowth

max_stars_repo_name: aangaji/TumorGrowth

max_issues_repo_name: aangaji/TumorGrowth

max_forks_repo_name: aangaji/TumorGrowth


text:
# CUSTOM NICE PRINTING FOR VARIOUS TYPES

# ABSTRACT TYPES

# single SITE

# all labels
function show(io::IO, s::S) where {L,D,S<:AbstractSite{L,D}}
    print(io, S, " @", point(s), ": ", label(s))
end
# string labels
function show(io::IO, s::S) where {L<:AbstractString,D,S<:AbstractSite{L,D}}
    print(io, S, " @", point(s), ": \"", label(s),"\"")
end
# symbol labels
function show(io::IO, s::S) where {L<:Symbol,D,S<:AbstractSite{L,D}}
    print(io, S, " @", point(s), ": :", label(s))
end

# single BOND

# all labels
function show(io::IO, b::B) where {L,N,B<:AbstractBond{L,N}}
    print(io, B, " ", from(b), "-->", to(b), " @", wrap(b), ": ", label(b))
end
# string labels
function show(io::IO, b::B) where {L<:AbstractString,N,B<:AbstractBond{L,N}}
    print(io, B, " ", from(b), "-->", to(b), " @", wrap(b), ": \"", label(b),"\"")
end
# symbol labels
function show(io::IO, b::B) where {L<:Symbol,N,B<:AbstractBond{L,N}}
    print(io, B, " ", from(b), "-->", to(b), " @", wrap(b), ": :", label(b))
end

# single UNITCELL
function show(io::IO, u::U) where {D,N,LS,LB,S<:AbstractSite{LS,D},B<:AbstractBond{LB,N},U<:AbstractUnitcell{S,B}}
    print(io, "Unitcell object\n--> type ", U, "\n--> ", length(sites(u)), " sites of type ", S, "\n--> ", length(bonds(u)), " bonds of type ", B)
end

# single LATTICE
function show(io::IO, la::L) where {D,N,LS,LB,U,S<:AbstractSite{LS,D},B<:AbstractBond{LB,N},L<:AbstractLattice{S,B,U}}
    print(io, "Lattice object\n--> type ", L, "\n--> ", length(sites(la)), " sites of type ", S, "\n--> ", length(bonds(la)), " bonds of type ", B, "\n--> unitcell of type ", U)
end


text:
using GemmKernels

using Pkg
Pkg.add(PackageSpec(name="XUnit", rev="9b756fcda72d813dbf017f8400d7c55251ef7d1b"))

using XUnit

import CUDA
import InteractiveUtils

@info "Julia details\n\n" * sprint(io->InteractiveUtils.versioninfo(io))
@info "CUDA details\n\n" * sprint(io->CUDA.versioninfo(io))

CUDA.allowscalar(false)

@testset runner=ParallelTestRunner() "GemmKernels.jl" begin
    include("tiling.jl")
    include("matmul.jl")
    include("blas.jl")
end


text:
###############################################################################
#
# An example of writing cell comments to a worksheet using XlsxWriter.
#
# For more advanced comment options see comments2.jl
#
# Original Python Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
# https://github.com/jmcnamara/XlsxWriter

using XlsxWriter

function test()

	wb = Workbook("comments1.xlsx")
	ws = add_worksheet!(wb)

	write!(ws, "A1", "Hello")
	write_comment!(ws, "A1", "This is a comment")

	close(wb)
	isfile("comments1.xlsx")

end

test()


text:
module RandomProjectionTree

using Printf

include("RPTree.jl")

end


text:
import Bijections

function part1(fname)
  unique_digits = (1, 4, 7, 8)
  unique_digit_lengths = (2, 4, 3, 7)
  n_unique_digits = 0

  open(fname, "r") do f
    data = readlines(f)
    for entry in data
      outputs = split(split(entry, '|')[2])
      length_outputs = [length(digit) for digit in outputs]
      n_unique_digits += sum([digit in unique_digit_lengths
                              for digit in length_outputs])
    end
  end
  println("1, 4, 7 and 8 appear in total $n_unique_digits times.")
end

function part2(fname)
  #=
  Strategy:
    - get 'cf', the 1
    - get 'acf', the 7
    - get 'bcdf', the 4
    - get 'abcdefg', the 8
    - get 'abcdfg', the 9, having the 4 as subset
    - get 'acdfg', the 3, having the 7 as subset
    - get 'abcefg', the 0, having the 1 as subset
    - get 'abdefg', the 6, the remaining digit of length 6.

    - get 'abdfg', the 5, having the 4 as subset
    - get 'acdeg', the 2, the remaining one
  =#
  open(fname, "r") do f
    data = readlines(f)
    sum_output = 0
    for entry in data
      signals, outputs = split.(split(entry, '|'))
      # Ordering of chars in signals and outputs needed for later comparisons.
      signals = join.(sort.(collect.(signals)))
      outputs = join.(sort.(collect.(outputs)))
      decoded = Bijections.Bijection{String, Int}()

      for (val, cond) in zip((1, 7, 4, 8, 9, 3, 0, 6, 2, 5), (
          x -> length(x) == 2, x -> length(x) == 3, x -> length(x) == 4,
          x -> length(x) == 7,
          x -> length(x)==6 && all([occursin(char, x) for char in decoded(4)]),
          x -> length(x)==5 && all([occursin(char, x) for char in decoded(7)]),
          x -> length(x)==6 && all([occursin(char, x) for char in decoded(1)]),
          x -> length(x) == 6,
          x -> occursin(decoded(8)[findfirst(x -> !occursin(x, decoded(9)), decoded(8))], x),
          x -> true
         ))
        decode_value!(decoded, signals, val, cond)
      end
      decoded_output = ""
      for output in outputs
        decoded_output *= string(decoded[output])
      end
      sum_output += parse(Int, decoded_output)
    end
    println("Sum of all outputs: $sum_output")
  end
end

function decode_value!(decoded, signals, value, condition)
  idx = findfirst(condition, signals)
  decoded[signals[idx]] = value
  deleteat!(signals, idx)
end

part1("test_input.txt")
part2("test_input.txt")

part1("input.txt")
part2("input.txt")