text:
@testset "972.equal-rational-numbers.jl" begin
    @test is_rational_equal("0.(52)", "0.5(25)") == true
    @test is_rational_equal("0.1666(6)", "0.166(66)") == true
    @test is_rational_equal("0.9(9)", "1.") == true
end


text:
function create_plan(img, u, v, padfac, m)
    dx,dy = pixelsizes(img)
    uv = Matrix{eltype(u)}(undef, 2, length(u))
    uv[1,:] = u*dx
    uv[2,:] = v*dy
    dx,dy = pixelsizes(img)
    ny,nx = size(img)
    nnx = nextpow(2, padfac*nx)
    nny = nextpow(2, padfac*ny)
    nsx = nnx÷2-nx÷2
    nsy = nny÷2-ny÷2
    pimg = PaddedView(zero(eltype(img)), img.im,
                      (1:nnx, 1:nny),
                      (nsx+1:nsx+nx, nsy+1:nsy+ny)
                     )
    plan = plan_nfft(uv, size(pimg), m)
    phases = @. exp(1im*π*(u*dx + v*dy))*dx*dy#*exp(1im*π*(u*x + v*y))*dx*dy
    return plan, phases, pimg
end

Base.@kwdef struct NFFT <: FourierTransform
    padfac::Int = 2
    m::Int = 10
end

struct NFFTCache{C} <: AbstractCache
    plan::C
    padfac::Int
    m::Int
end

ext: jl


text:
using .JuMP: @constraint, Model
import .JuMP

"""
    add_constraint!(model::JuMP.Model, X, Set)

Add constraint `xᵢ ∈ Set` where `xᵢ` is a column of `X` (or if `X` is a vector `xᵢ=X`).

### Input
- `model` -- `JuMP.Model` type
- `X` -- variable vector `X[1:n]` or variable matrix `X`[1:n, 1:N]` where `n` is the number of states and
`N` is the number of time steps
- `Set` -- lazy set or vector type

### Output
Nothing.

### Note
Only works if `X` is a `Array{VariableRef,1}` or `VariableRef` variable types,
`DenseAxisArrays` are not yet supported yet.

The following types for the constraint `Set` are supported:
- polyhedral sets →  Inequality constraint `Fxᵢ ≤ g` where `F, g = constraints_list(P)`.
- Vector          → Equality constraint `xᵢ==Set`
- Singleton       → Equality constraint `xᵢ==element(Set)`
- Universe        → No constraint
"""
function add_constraint!(model::JuMP.Model, X, Set) end

function add_constraint!(model::JuMP.Model, X, A::AbstractVector{<:Real})
    _throw_dim_mismatch(X, A)

    for j=1:size(X,2)
        JuMP.@constraint(model, X[:, j] .== A)
    end
    nothing
end

function add_constraint!(model::JuMP.Model, X, S::Singleton)
    _throw_dim_mismatch(X, S)

    for j=1:size(X,2)
        JuMP.@constraint(model, X .== element(S))
    end
    nothing
end

function add_constraint!(model::JuMP.Model, X, U::Universe)
    _throw_dim_mismatch(X, U)
    # do nothing
    nothing
end

function add_constraint!(model::JuMP.Model, X, P::LazySet)
    @assert applicable(constraints_list, P) "the number of constraint of P needs "*
                 "to be finite; try overapproximating P with an `HPolytope` first"
    _throw_dim_mismatch(X, P)

    F, g = tosimplehrep(P)
    isempty(P) && return nothing
    for j=1:size(X,2)
        JuMP.@constraint(model, F*X[:, j] .<= g)
    end
    nothing
end

_throw_dim_mismatch(X, P::LazySet) =  InvariantSets.dim(P) != size(X, 1) &&
    throw(ArgumentError("the dimension of the constraint Set $(InvariantSets.dim(P)) "*
                        "doesn't match with the rows of the variable array with size $(size(X))"))
_throw_dim_mismatch(X, A::AbstractArray) =  length(A) != size(X, 1) &&
    throw(ArgumentError("the dimension of the constraint Set $(length(A)) "*
                        "doesn't match with the rows of the variable array with size $(size(X))"))

max_stars_repo_name: mforets/InvariantSets

max_issues_repo_name: mforets/InvariantSets


text:
@testset "environment variables" begin
    config = DockerConfig(;
        base_image = "julia:latest",
        env = Dict{String, String}("MY_ENVIRONMENT_VARIABLE" => "hello_world")
    )

    with_container() do container
        code = """
        println("This was a success.")
        if !haskey(ENV, "MY_ENVIRONMENT_VARIABLE")
            throw(ErrorException("environment variable does not exist"))
        end
        if ENV["MY_ENVIRONMENT_VARIABLE"] != "hello_world"
            throw(ErrorException("environment variable has the wrong value"))
        end
        """
        @test success(container, config, `julia -e $(code)`)
    end
end


text:
using LazySets.Approximations: AbstractDirections

"""
    LGG09{N, AM, TN<:AbstractDirections} <: AbstractContinuousPost

Implementation of Girard - Le Guernic algorithm for reachability analysis
using support functions.

## Fields

- `δ`            -- step-size of the discretization
- `approx_model` -- (optional, default: `Forward`) approximation model;
                    see `Notes` below for possible options
- `template`     -- (alias: `dirs`) struct that holds the directions (either lazily or concretely)
                    for each support function evaluation defining the template
- `vars`         -- (optional, default: all variables are computed) an integer or a vector of integers
                    specifying the variables of interest to automatically construct a template
                    using canonical directions; requires that `n` (or `dim`) is specified as well
- `static`       -- (optional, default: `false`) if `true`, use statically sized arrays
- `threaded`     -- (optional, default: `false`) if `true`, use multi-threading
                    to compute different template directions in parallel
- `sparse`       -- (optional, default: `false`) if `true`, convert the matrix exponential
                    obtained after discretization to a sparse matrix
- `cache`        -- (optional, default: `true`) if `true`, use a cache for intermediate
                    computations in the set recurrence loop
- `vars`         -- (optional, default: `missing`) used to specify the variables instead
                    of passing the template

## Notes

The type fields are:

- `N`        -- number type of the step-size
- `AM`       -- type of the approximation model
- `TN`       -- type of the abstract directions that define the template

The flag `threaded=true` specifies that the support functions
along different directions should be computed in parallel.
See the section on [Multi-threading](@ref) for details on how to setup
the number of threads.

## References

This is an implementation of the algorithm from [[LGG09]](@ref).

These methods are described at length in the dissertation [[LG09]](@ref).
"""
struct LGG09{N, AM, VN, TN<:AbstractDirections{N, VN}, S, VA} <: AbstractContinuousPost
    δ::N
    approx_model::AM
    template::TN
    static::S
    threaded::Bool
    sparse::Bool
    cache::Bool
    vars::VA
end

# convenience constructor using symbols
function LGG09(; δ::N,
               approx_model::AM=Forward(sih=:concrete, exp=:base, setops=:lazy),
               template=nothing,
               dirs=nothing, # alias for template
               vars=missing, # shortcut to specify variables of interest
               n=missing,    # required to construct the template from vars
               dim=missing,  # alias for n
               static::Bool=false,
               threaded::Bool=false,
               sparse::Bool=false,
               cache::Bool=true) where {N, AM}

    if !ismissing(vars)
        if !ismissing(n)
            _n = n
        elseif !ismissing(dim)
            _n = dim
        else
            throw(ArgumentError("the ambient dimension, `n` (or `dim`), should be specified"))
        end
        directions = _get_template(vars, _n)
    elseif !isnothing(dirs)
        directions = _get_template(dirs)
    elseif !isnothing(template)
        directions = _get_template(template)
    else
        throw(ArgumentError("either `vars`, `dirs` or `template` should be specified"))
    end

    return LGG09(δ, approx_model, directions, Val(static), threaded, sparse, cache, vars)
end

_get_template(template::AbstractDirections) = template
_get_template(template::AbstractVector{N}) where {N<:Number} = CustomDirections([template])
_get_template(template::AbstractVector{VT}) where {N<:Number, VT<:AbstractVector{N}} = CustomDirections(template)

_get_template(vars::Int, n::Int) = _get_template((vars,), n)

function _get_template(vars::VecOrTuple, n::Int)
    m = length(vars)
    @assert m ≤ n "the number of variables should not exceed `n`"
    dirs = Vector{SingleEntryVector{Float64}}(undef, 2*m)

    @inbounds for (i, vi) in enumerate(vars)
        d₊ = SingleEntryVector(vi, n, 1.0)
        d₋ = SingleEntryVector(vi, n, -1.0)
        dirs[i] = d₊
        dirs[i+m] = d₋
    end
    return CustomDirections(dirs, n)
end

step_size(alg::LGG09) = alg.δ
numtype(::LGG09{N}) where {N} = N
setrep(alg::LGG09) = setrep(alg.template)

function rsetrep(alg::LGG09{N, AM, VN, TN}) where {N, AM, VN, TN}
    SN = SubArray{N, 1, Matrix{N}, Tuple{Base.Slice{Base.OneTo{Int}}, Int}, true}
    return TemplateReachSet{N, VN, TN, SN}
end

include("reach_homog.jl")
include("reach_inhomog.jl")
include("post.jl")


text:
using Random,Distributions,Polynomials
using Statistics,QuadGK,FFTW
using Optim,Roots

Random.seed!(150000)

#ZSB 零息债券定价参数 gama,b1,b2
function gamma(κ_r::Float64,σ_r::Float64)
    return √(κ_r^2+2*σ_r^2)
end

function b1(κ_r::Float64,θ_r::Float64,σ_r::Float64,r0::Float64,T::Float64)
    g=gamma(κ_r,σ_r)
    result=((2*g*exp((κ_r+g)*T/2))/(2*g+(κ_r+g)*(exp(g*T)-1)))^(2*κ_r*θ_r/σ_r^2)

    return result
end

function b2(κ_r::Float64,θ_r::Float64,σ_r::Float64,r0::Float64,T::Float64)
    g=gamma(κ_r,σ_r)
    result=2*(exp(g*T)-1)/(2*g+(κ_r+g)*(exp(g*T)-1))
    return result
end

function B(κ_r::Float64,θ_r::Float64,σ_r::Float64,r0::Float64,T::Float64)

    b_1=b1(κ_r,θ_r,σ_r,r0,T)
    b_2=b2(κ_r,θ_r,σ_r,r0,T)
    result=b_1*exp(-b_2*r0)
    return result
end
###################end###############

#H93模型特征函数值
function H93_char_func(u,T::Float64,r::Float64,κ_v::Float64,θ_v::Float64,σ_v::Float64,ρ::Float64,v0::Float64)
    c1=κ_v*θ_v
    c2=-√((ρ*σ_v*u*1im-κ_v)^2-σ_v^2*(-u*1im-u^2))
    c3=(κ_v-ρ*σ_v*u*1im+c2)/(κ_v-ρ*σ_v*u*1im-c2)

    H1=(r*u*1im*T+(c1/σ_v^2)*((κ_v-ρ*σ_v*u*1im+c2)*T-2*log((1-c3*exp(c2*T))/(1-c3))))
    H2=((κ_v-ρ*σ_v*u*1im+c2)/σ_v^2*((1-exp(c2*T))/(1-c3*exp(c2*T))))

    char_fun_value=exp(H1+H2*v0)

    return char_fun_value
end

#M76模型特征函数值
function M76_char_func(u,T::Float64,λ::Float64,μ::Float64,δ::Float64)
        ω=-λ*(exp(μ+0.5*(δ^2))-1)
        tmp1=λ*(map(x->exp(x)-1,1im*u*μ .-0.5*(u.^2)*δ^2))
        value=map(x->exp(x*T),(1im*u*ω .+tmp1))
        return value
end

#BCC模型特征函数值
function BCC_char_func(u,T::Float64,r::Float64,κ_v::Float64,θ_v::Float64,σ_v::Float64,ρ::Float64,v0::Float64,λ::Float64,μ::Float64,δ::Float64)
    BCC1=H93_char_func(u,T,r,κ_v,θ_v,σ_v,ρ,v0)
    BCC2=M76_char_func(u,T,λ,μ,δ)

    return BCC1*BCC2
end

#H93模型内在价值
function H93_int_func(u,S0,K,T,r,κ_v,θ_v,σ_v,ρ,v0)
    char_fun_value=H93_char_func(u-1im*0.5,T,r,κ_v,θ_v,σ_v,ρ,v0)
    int_func_value=1/(u^2+0.25)*real(exp(1im*u*log(S0/K))*char_fun_value)
    return int_func_value
end

#BCC模型内在价值函数值
function BCC_int_func(u,S0,K,T,r,κ_v,θ_v,σ_v,ρ,v0,λ,μ,δ)
    char_fun_value=BCC_char_func(u-1im*0.5,T,r,κ_v,θ_v,σ_v,ρ,v0,λ,μ,δ)
    int_func_value=1/(u^2+0.25)*real(exp(1im*u*log(S0/K))*char_fun_value)

    return int_func_value
end

#H93模型欧式看涨期权值
function H93_call_func(S0::Float64,K::Float64,T::Float64,r::Float64,κ_v::Float64,θ_v::Float64,σ_v::Float64,ρ::Float64,v0)
    int_func_value=quadgk(x->H93_int_func(x,S0,K,T,r,κ_v,θ_v,σ_v,ρ,v0),0,Inf)[1]
    print("h93_int:",int_func_value)
    print("temp1:",exp(-r*T))
    print("temp2:",√(S0*K)/π*int_func_value)

    call_value=max(0,S0-exp(-r*T)*√(S0*K)/π*int_func_value)
    return call_value
end

#BCC模型欧式看涨期权值
function BCC_call_value(S0::Float64,K::Float64,T::Float64,r::Float64,κ_v::Float64,θ_v::Float64,σ_v::Float64,ρ::Float64,v0::Float64,λ::Float64,μ::Float64,δ::Float64)
    int_value=quadgk(x->BCC_int_func(x,S0,K,T,r,κ_v,θ_v,σ_v,ρ,v0,λ,μ,δ),0,Inf)[1]
    call_value=max(0,S0-exp(-r*T)*√(S0*K)/π*int_value)

    return call_value
end
#
# κ_v,θ_v,σ_v,v0,T=1.5,0.02,0.15,0.01,1.0
# ρ=0.1
# λ = 0.25
# μ = -0.2
# δ = 0.1
# σ = √(v0)
# S0 = 100.0
# K = 100.0
# T = 1.0
# r = 0.05
#
#
# # B_res=BCC_call_value(S0,K,T,r,κ_v,θ_v,σ_v,ρ,v0,λ,μ,δ)
# h93=H93_call_func(S0,K,T,r,κ_v,θ_v,σ_v,ρ,v0)
# # h93=π
# print(h93)
# u=2
# h93=H93_char_func(u,T,r,κ_v,θ_v,σ_v,ρ,v0)
# h93_int=H93_int_func(u,S0,K,T,r,κ_v,θ_v,σ_v,ρ,v0)
# print("char",h93)
# print("int",h93_int)


text:
module KLDivergence

using ..DeterministicVI: ElboArgs, VariationalParams, init_thread_pool!
using ...Model: CanonicalParams, ids, prior, NUM_SOURCE_TYPES, NUM_COLOR_COMPONENTS
using ...SensitiveFloats: SensitiveFloat, add_sources_sf!
using Compat
using ForwardDiff, ReverseDiff, DiffBase

# Calculate the Kullback-Leibler divergences between pairs of well-known parameteric
# distributions, and derivatives with respect to the parameters of the first distributions.

######################
# Beta Distributions #
######################

"""
Returns the KL divergence between a pair of beta distributions

Args:
- α₁: the shape parameter for the first beta distribution
- β₁: the scale parameter for the first beta distribution
- α₂: the shape parameter for the second beta distribution
- β₂: the scale parameter for the second beta distribution
"""
function beta_kl(α₁, β₁, α₂, β₂)
    α₁_plus_β₁ = α₁ + β₁
    α₂_plus_β₂ = α₂ + β₂
    α₁_minus_α₂ = α₁ - α₂
    β₁_minus_β₂ = β₁ - β₂
    log_term = lgamma(α₁_plus_β₁) - lgamma(α₁) - lgamma(β₁)
    log_term -= lgamma(α₂_plus_β₂) - lgamma(α₂) - lgamma(β₂)
    apart_term = α₁_minus_α₂ * digamma(α₁) + β₁_minus_β₂ * digamma(β₁)
    together_term = -(α₁_minus_α₂ + β₁_minus_β₂) * digamma(α₁_plus_β₁)
    kl = log_term + apart_term + together_term
    assert(!isnan(kl))
    return kl
end

#############################
# Categorical Distributions #
#############################

"""
Returns the KL divergence between a pair of categorical distributions
"""
function categorical_kl(p₁, p₂)
    kl = sum(a * (log(a) - log(b)) for (a, b) in zip(p₁, p₂))
    assert(!isnan(kl))
    return kl
end

##########################
# Gaussian Distributions #
##########################

"""
Returns the KL divergence between a pair of univariate Gaussian distributions

Args:
- μ₁: the mean for the first Gaussian distribution
- var₁: the variance for the first Gaussian distribution
- μ₂: the mean for the second Gaussian distribution
- var₂: the variance for the second Gaussian distribution
"""
gaussian_kl(μ₁, var₁, μ₂, var₂) = .5 * (log(var₂) - log(var₁) + (var₁ + (μ₁ - μ₂)^2) / var₂ - 1)

#####################################
# Multivariate Normal Distributions #
#####################################

"""
Returns the KL divergence between a pair of multivariate normal distributions,
the first having a diagonal covariance matrix

Args:
- μ₁: the mean for the first normal distribution
- var₁: the variance parameters for the first normal distribution
- μ₂: the mean for the second normal distribution
- Σ₂: the covariance matrix for the second normal distribution
"""
function diagmvn_mvn_kl(μ₁, var₁, μ₂, Σ₂, inv_Σ₂ = inv(Σ₂), logdet_Σ₂ = logdet(Σ₂))
    μ₂_minus_μ₁ = μ₂ - μ₁
    kl = sum(diag(inv_Σ₂) .* var₁) - length(μ₂)
    kl += dot(μ₂_minus_μ₁, inv_Σ₂ * μ₂_minus_μ₁)
    kl += logdet_Σ₂ - sum(log.(var₁))
    assert(!isnan(kl))
    return 0.5 * kl
end

###########################################
# Subtracting KL divergences from sources #
###########################################

kl_source_a(vs) = categorical_kl(vs[ids.is_star], prior.is_star)

function kl_source_r(vs)
    kl = zero(eltype(vs))
    for i in 1:NUM_SOURCE_TYPES
        kl += vs[ids.is_star[i]] * gaussian_kl(vs[ids.flux_loc[i]], vs[ids.flux_scale[i]],
                                            prior.flux_mean[i], prior.flux_var[i])
    end
    assert(!isnan(kl))
    return kl
end

function kl_source_k(vs)
    kl = zero(eltype(vs))
    for i in 1:NUM_SOURCE_TYPES
        kl += vs[ids.is_star[i]] * categorical_kl(vs[ids.k[:, i]], prior.k[:, i])
    end
    assert(!isnan(kl))
    return kl
end

function kl_source_c(vs)
    kl = zero(eltype(vs))
    for i in 1:NUM_SOURCE_TYPES
        μ₁, var₁ = vs[ids.color_mean[:, i]], vs[ids.color_var[:, i]]
        a = vs[ids.is_star[i]]
        for d in 1:NUM_COLOR_COMPONENTS
            μ₂, Σ₂ = prior.color_mean[:, d, i], prior.color_cov[:, :, d, i]
            kl += a * vs[ids.k[d, i]] * diagmvn_mvn_kl(μ₁, var₁, μ₂, Σ₂)
        end
    end
    assert(!isnan(kl))
    return kl
end

function source_e_log_prob(vs)
    x = vs[ids.gal_radius_px]
    μ = prior.gal_radius_px_mean
    σ² = prior.gal_radius_px_var
    kl = -0.5 * (log(2pi) + log(σ²) + (x - μ)^2 / σ²)
    assert(!isnan(kl))
    return kl
end

"""
Subtract the KL divergences for a single source.
"""
function subtract_kl(vs)
    kl = zero(eltype(vs))
    kl -= kl_source_a(vs)
    kl -= kl_source_k(vs)
    kl -= kl_source_r(vs)
    kl -= kl_source_c(vs)

    # negative log probability is the kl divergence between a
    # variational distribution that is a point mass and the prior
    kl += source_e_log_prob(vs)
    return kl
end

###################
# Differentiation #
###################

using ForwardDiff: Dual, JacobianConfig, pickchunksize
using ReverseDiff: TrackedArray, TrackedReal, CompiledTape, GradientTape

const KLGradientTape{T} = GradientTape{typeof(subtract_kl),
                                       TrackedArray{T,T,1,Vector{T},Vector{T}},
                                       TrackedReal{T,T,Void}}

struct KLHelper{N,T}
    gradient_tape::CompiledTape{:kl_gradient,KLGradientTape{T}}
    nested_gradient_tape::CompiledTape{:kl_nested_gradient,KLGradientTape{Dual{N,T}}}
    dual_buffer::Vector{Dual{N,T}}
    jacobian_config::JacobianConfig{N,T,Vector{Dual{N,T}}}
end

const PARAM_LENGTH = length(CanonicalParams)

function KLHelper(::Type{Dual{N,T}}) where {N,T}
    dual_buffer = zeros(Dual{N,T}, PARAM_LENGTH)
    jacobian_config = JacobianConfig{N}(rand(T, PARAM_LENGTH))
    gradient_tape = CompiledTape{:kl_gradient}(GradientTape(subtract_kl, rand(T, PARAM_LENGTH)))
    nested_gradient_tape = CompiledTape{:kl_nested_gradient}(GradientTape(subtract_kl, rand(Dual{N,T}, PARAM_LENGTH)))
    ReverseDiff.compile(gradient_tape)
    ReverseDiff.compile(nested_gradient_tape)
    return KLHelper{N,T}(gradient_tape, nested_gradient_tape, dual_buffer, jacobian_config)
end

function kl_gradient!(out, x, helper::KLHelper)
    return ReverseDiff.gradient!(out, helper.gradient_tape, x)
end

function kl_hessian!(out, x, helper::KLHelper)
    f = t -> ReverseDiff.gradient!(helper.dual_buffer, helper.nested_gradient_tape, t)
    return ForwardDiff.jacobian!(out, f, x, helper.jacobian_config)
end

###############
# Entry Point #
###############

function subtract_kl_source!(kl_source::SensitiveFloat, result::DiffBase.DiffResult,
                             vs, helper::KLHelper)
    if kl_source.has_gradient
        kl_gradient!(result, vs, helper)
        kl_source.v[] = DiffBase.value(result)
    else
        kl_source.v[] = subtract_kl(vs)
    end
    if kl_source.has_hessian
        kl_hessian!(kl_source.h, vs, helper)
    end
    return kl_source
end

function subtract_kl_all_sources!(ea::ElboArgs,
                                  vp::VariationalParams{T},
                                  accum::SensitiveFloat,
                                  kl_source::SensitiveFloat = get_kl_source(T),
                                  helper::KLHelper = get_kl_helper(T)) where {T}
    result = DiffBase.DiffResult(zero(T), kl_source.d)
    for sa in 1:length(ea.active_sources)
        subtract_kl_source!(kl_source, result, vp[ea.active_sources[sa]], helper)
        add_sources_sf!(accum, kl_source, sa)
    end
    return accum
end

############
# __init__ #
############

const CHUNK_SIZE = pickchunksize(PARAM_LENGTH)

const KL_HELPER_FLOAT_POOL = Vector{KLHelper{CHUNK_SIZE,Float64}}()
const KL_HELPER_DUAL_POOL = Vector{KLHelper{CHUNK_SIZE,Dual{1,Float64}}}()

const KL_SOURCE_FLOAT_POOL = Vector{SensitiveFloat{Float64}}()
const KL_SOURCE_DUAL_POOL = Vector{SensitiveFloat{Dual{1,Float64}}}()

get_kl_helper(::Type{Float64}) = KL_HELPER_FLOAT_POOL[Base.Threads.threadid()]
get_kl_helper(::Type{Dual{1,Float64}}) = KL_HELPER_DUAL_POOL[Base.Threads.threadid()]

get_kl_source(::Type{Float64}) = KL_SOURCE_FLOAT_POOL[Base.Threads.threadid()]
get_kl_source(::Type{Dual{1,Float64}}) = KL_SOURCE_DUAL_POOL[Base.Threads.threadid()]

function __init__()
    init_thread_pool!(KL_HELPER_FLOAT_POOL, () -> KLHelper(Dual{CHUNK_SIZE,Float64}))
    init_thread_pool!(KL_HELPER_DUAL_POOL, () -> KLHelper(Dual{CHUNK_SIZE,Dual{1,Float64}}))
    init_thread_pool!(KL_SOURCE_FLOAT_POOL, () -> SensitiveFloat{Float64}(PARAM_LENGTH, 1, true, true))
    init_thread_pool!(KL_SOURCE_DUAL_POOL, () -> SensitiveFloat{Dual{1,Float64}}(PARAM_LENGTH, 1, true, false))
end

# explicitly call this for use with compiled system image
__init__()

end # module


text:
using Oceananigans:
    CPU, GPU, AbstractGrid, AbstractPoissonSolver,
    BC, Periodic, ModelBoundaryConditions

using Oceananigans.Grids: RegularCartesianGrid

# PoissonBCs are named XYZ, where each of X, Y, and Z is either
# 'P' (for Periodic) or 'N' (for Neumann).
struct PPN <: PoissonBCs end
struct PNN <: PoissonBCs end

# Not yet supported:
#struct PPP <: PoissonBCs end
#struct NNN <: PoissonBCs end
#struct NPP <: PoissonBCs end
#struct PNP <: PoissonBCs end
#struct NPN <: PoissonBCs end
#struct NNP <: PoissonBCs end

poisson_bc_symbol(::BC) = :N
poisson_bc_symbol(::BC{<:Periodic}) = :P

"""
    PoissonBCs(bcs)

Returns the boundary conditions for the Poisson solver corresponding
to the model boundary conditions `bcs`.
"""
function PoissonBCs(bcs)
    # We assume that bounary conditions on all fields are
    # consistent with the boundary conditions on one side of ``u``.
    x = poisson_bc_symbol(bcs.u.x.left)
    y = poisson_bc_symbol(bcs.u.y.left)
    z = poisson_bc_symbol(bcs.u.z.left)

    return eval(Expr(:call, Symbol(x, y, z)))
end

PoissonBCs(model_bcs::ModelBoundaryConditions) = PoissonBCs(model_bcs.solution)

PoissonSolver(::CPU, pbcs::PoissonBCs, grid::AbstractGrid) = PoissonSolverCPU(pbcs, grid)
PoissonSolver(::GPU, pbcs::PoissonBCs, grid::AbstractGrid) = PoissonSolverGPU(pbcs, grid)

unpack_grid(grid::AbstractGrid) = grid.Nx, grid.Ny, grid.Nz, grid.Lx, grid.Ly, grid.Lz

"""
    ω(M, k)

Return the `M`th root of unity raised to the `k`th power.
"""
@inline ω(M, k) = exp(-2im*π*k/M)

"""
    λi(grid::AbstractGrid, ::PoissonBCs)

Return an Nx×1×1 array of eigenvalues satisfying the discrete form of Poisson's
equation with periodic boundary conditions in the x-dimension on `grid`.
"""
function λi(grid::AbstractGrid, ::PoissonBCs)
    Nx, Ny, Nz, Lx, Ly, Lz = unpack_grid(grid)
    is = reshape(1:Nx, Nx, 1, 1)
    @. (2sin((is-1)*π/Nx) / (Lx/Nx))^2
end

"""
    λj(grid::AbstractGrid, ::PPN)

Return an 1×Ny×1 array of eigenvalues satisfying the discrete form of Poisson's
equation with periodic boundary conditions in the y-dimension on `grid`.
"""
function λj(grid::AbstractGrid, ::PPN)
    Nx, Ny, Nz, Lx, Ly, Lz = unpack_grid(grid)
    js = reshape(1:Ny, 1, Ny, 1)
    @. (2sin((js-1)*π/Ny) / (Ly/Ny))^2
end

"""
    λj(grid::AbstractGrid, ::PNN)

Return an 1×Ny×1 array of eigenvalues satisfying the discrete form of Poisson's
equation with staggered Neumann boundary conditions in the y-dimension on `grid`.
"""
function λj(grid::AbstractGrid, ::PNN)
    Nx, Ny, Nz, Lx, Ly, Lz = unpack_grid(grid)
    js = reshape(1:Ny, 1, Ny, 1)
    @. (2sin((js-1)*π/(2Ny)) / (Ly/Ny))^2
end

"""
    λk(grid::AbstractGrid, ::PoissonBCs)

Return an 1×1×Nz array of eigenvalues satisfying the discrete form of Poisson's
equation with staggered Neumann boundary conditions in the y-dimension on `grid`.
"""
function λk(grid::AbstractGrid, ::PoissonBCs)
    Nx, Ny, Nz, Lx, Ly, Lz = unpack_grid(grid)
    ks = reshape(1:Nz, 1, 1, Nz)
    @. (2sin((ks-1)*π/(2Nz)) / (Lz/Nz))^2
end


text:
@testset "PREBUILD" begin

  for f in readdir(joinpath(@__DIR__, "test_bin"))
    file = joinpath(@__DIR__, "test_bin", f)
    x1,   x2 = Pickle.Torch.THload(file)
    px1, px2 = libtorch.load(file)
    @test x1[2] ≈ px1[2].numpy()
    @test x2[2] ≈ px2[2].numpy()
  end

end


text:
#=

T intersection
    ________________________________________________________
        <---- seg 2     main street <----- seg 1
        -----> seg 3                 -----> seg 4
    _________________            ___________________________
                    | seg . seg  |
                    |  5  .  6   |
                    |     .      |
                    |     .      |
                    |     .      |

=#

# State type
struct PedCarMDPState
    crash::Bool
    ego::VehicleState
    ped::VehicleState
    car::VehicleState
    route::SVector{2, LaneTag}
end

function Base.:(==)(v1::VehicleState, v2::VehicleState)
    return v1.posG == v2.posG && v1.v == v2.v
end
function Base.isapprox(v1::VehicleState, v2::VehicleState)
    return v1.posG ≈ v2.posG && v1.v ≈ v2.v
end
Base.hash(veh::VehicleState, h::UInt) = hash(veh.posF, hash(veh.v, h))

# copy b to a
function Base.copyto!(a::PedCarMDPState, b::PedCarMDPState)
    a.crash = b.crash
    a.ego = b.ego
    a.ped = b.ped
    a.car = b.car
    a.route = b.route
end

function Base.hash(s::PedCarMDPState, h::UInt64 = zero(UInt64))
    return hash(s.crash, hash(s.ego, hash(s.ped, hash(s.car, hash(s.route, h)))))
end

function Base.:(==)(a::PedCarMDPState, b::PedCarMDPState)
    return a.crash == b.crash && a.ego == b.ego && a.ped == b.ped && a.car == b.car && a.route == b.route
end

# Action type
const PedCarMDPAction = UrbanAction
const CarRoute = SVector{2, LaneTag}
const CarState = Tuple{VehicleState, CarRoute}
const OFF_ROUTE = SVector{2, LaneTag}(LaneTag(0,0), LaneTag(0,0))
const N_NEXT = 12
const N_CAR_START = 12 # for v_res=2.0, had to hardcode for performance
const N_PED_START = 12 # for v_ped_res = 1.0
const CarTransitionDict = Dict{Tuple{VehicleState, CarRoute, LonAccelDirection}, SparseCat{Vector{Tuple{VehicleState, CarRoute}}, Vector{Float64}}}
const PedTransitionDict = Dict{Tuple{VehicleState, ConstantSpeedDawdling}, SparseCat{Vector{VehicleState}, Vector{Float64}}}
const EgoTransitionDict = Dict{Tuple{VehicleState, LonAccelDirection}, SparseCat{Vector{VehicleState}, Vector{Float64}}}

"""
    PedCarMDP
Implementation of a driving scenario involving three agents (an ego vehicle, another car and a pedestrian) at a T-shape intersections.
"""
@with_kw mutable struct PedCarMDP <: MDP{PedCarMDPState, PedCarMDPAction}
    env::UrbanEnv = UrbanEnv(params=UrbanParams(nlanes_main=1,
                crosswalk_pos =  [VecSE2(6, 0., pi/2), VecSE2(-6, 0., pi/2), VecSE2(0., -5., 0.)],
                crosswalk_length =  [14.0, 14., 14.0],
                crosswalk_width = [4.0, 4.0, 3.1],
                stop_line = 22.0))
    ΔT::Float64 = 0.5
    pos_res::Float64 = 2.0
    vel_res::Float64 = 2.0
    vel_ped_res::Float64 = 1.0
    car_action_space::Vector{Float64} = collect(-9.0:1.0:4.0)
    ped_action_space::Vector{Float64} = collect(0.0:1.0:2.0)
    car_models::Dict{SVector{2, LaneTag}, DriverModel} = get_car_models(env, get_ttc_model)
    car_type::VehicleDef = VehicleDef()
    ego_type::VehicleDef = VehicleDef()
    ped_type::VehicleDef = VehicleDef(AgentClass.PEDESTRIAN, 1.0, 1.0)
    a_noise::Float64 = 1.0
    v_noise::Float64 = 1.0
    ped_birth::Float64 = 0.7
    car_birth::Float64 = 0.7
    ego_start::Float64 = env.params.stop_line - ego_type.length/2
    ego_goal::LaneTag = LaneTag(2,1)
    off_grid::VehicleState = VehicleState(VecSE2(UrbanEnv().params.x_min+VehicleDef().length/2, env.params.y_min+VehicleDef().width/2, 0), Frenet(env.roadway[LaneTag(5,1)], 25.1, -26.5, pi/2), 0.)
    collision_cost::Float64 = -1.0
    action_cost::Float64 = 0.
    goal_reward::Float64 = 1.
    γ::Float64 = 0.95
    _ped_grid::Dict{LaneTag, RectangleGrid{3}} = init_ped_grid(env, pos_res, vel_ped_res)
    _car_grid::Dict{LaneTag, RectangleGrid{2}} = init_car_grid(env, pos_res, vel_res)
    _l_grid::Dict{LaneTag, RectangleGrid{1}} = init_l_grid(env, pos_res)
    _v_grid::RectangleGrid{1} = init_v_grid(env, vel_res)
    _car_transition_dict::CarTransitionDict = CarTransitionDict()
    _ped_transition_dict::PedTransitionDict = PedTransitionDict()
    _ego_transition_dict::EgoTransitionDict = EgoTransitionDict()
    _collision_checker::Dict{Tuple{Vararg{VehicleState, 3}}, Bool} = Dict{Tuple{Vararg{VehicleState, 3}}, Bool}()
end


text:
using Test, BenchmarkTools

input = 9798

xs = 1:300
ys = (1:300)'

rack_id = xs .+ 10 # Find the fuel cell's rack ID, which is its X coordinate plus 10.
# Begin with a power level of the rack ID times the Y coordinate.
# Increase the power level by the value of the grid serial number (your puzzle input).
# Set the power level to itself multiplied by the rack ID.
# Keep only the hundreds digit of the power level (so 12345 becomes 3; numbers with no hundreds digit become 0).
# Subtract 5 from the power level.

function star1(input, sizes=1:300)
    xs = (1:300) .+ 10
    ys = (1:300)'

    power = ((((xs .* ys .+ input) .* xs) .÷ 100) .% 10) .- 5

    maxpower = typemin(Int)
    maxxysz = (0,0,0)
    for sz in sizes
        for x in 1:(300-sz+1)
            for y in 1:(300-sz+1)
                pow = sum(power[xx,yy] for xx in x:(x+sz-1), yy in y:(y+sz-1))
                maxpower = max(maxpower, pow)
                maxxysz = maxpower == pow ? (x,y,sz) : maxxysz
            end
        end
    end
    return maxxysz, maxpower
end

star1(18, 3)

@test star1(18) == ((33,45), 29)
@test star1(42) == ((21,61), 30)

xy, maxp = star1(input)


text:
using DiffEqParamEstim, OrdinaryDiffEq, StochasticDiffEq, ParameterizedFunctions,
      DiffEqBase, RecursiveArrayTools
using Test

pf_func = function (du,u,p,t)
  du[1] = p[1] * u[1] - p[2] * u[1]*u[2]
  du[2] = -3 * u[2] + u[1]*u[2]
end

u0 = [1.0;1.0]
tspan = (0.0,10.0)
p = [1.5,1.0]
prob = ODEProblem(pf_func,u0,tspan,p)
sol = solve(prob,Tsit5())

t = collect(range(0, stop=10, length=200))
randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)

monte_prob = EnsembleProblem(prob)
obj = build_loss_objective(monte_prob,Tsit5(),L2Loss(t,data),maxiters=10000,
                                                      verbose=false,num_monte=5)

import Optim
result = Optim.optimize(obj, [1.3,0.8], Optim.BFGS())
@test result.minimizer ≈ [1.5,1.0] atol=3e-1

pg_func = function (du,u,p,t)
  du[1] = 1e-6u[1]
  du[2] = 1e-6u[2]
end
prob = SDEProblem(pf_func,pg_func,u0,tspan,p)
sol = solve(prob,SRIW1())

monte_prob = EnsembleProblem(prob)

# Too stochastic for CI
#=
srand(200)
obj = build_loss_objective(monte_prob,SRIW1(),L2Loss(t,data),maxiters=1000,
                           verbose=false,verbose_opt=false,verbose_steps=1,num_monte=50)

result = Optim.optimize(obj, [1.4,0.95], Optim.BFGS())
@test result.minimizer ≈ [1.5,1.0] atol=3e-1
=#


text:
# Utility functions for processing data and seismograms

"""
    apply_split!(n, e, phi, dt)

Apply a splitting operator to two traces.  `n` is assumed to be oriented north,
and `e` east.
"""
function apply_split!(n::T, e::T, phi::Number, dt::Number) where {T<:Seis.Trace}
    rotate_traces!(n, e, -phi)
    # Move the fast trace forward
    shift_trace!(n, -dt)
    rotate_traces!(n, e, phi)
    nothing
end

"""
    array_shift!(a, n, pad=zero(eltype(a))) -> a

Shift the values within an array back by `n` places, and replace
the shifted values with `pad`.  This function
works similarly to `circshift`, but rather
than elements shifting back round to the other end of the array when
shifted off one end, replace the shifted elements with `pad`.
"""
function array_shift!(a, n, pad=zero(eltype(a)))
    n == 0 && return a
    N = length(a)
    # Just zero out trace if shifted the whole way
    if abs(n) >= N
        a .= pad
        return a
    end
    if n > 0
        @inbounds for i in N:-1:(n+1)
            a[i] = a[i-n]
        end
        @inbounds a[1:n] .= pad
    elseif n < 0
        @inbounds for i in 1:(N+n)
            a[i] = a[i-n]
        end
        @inbounds a[end+n+1:end] .= pad
    end
    a
end

"""
    check_trace_order(t1, t2) -> a, b

Return the traces `t1` and `t2` as `a` and `b`, where `b` is always clockwise of `a`.
"""
check_trace_order(t1, t2) = Seis.angle_difference(t1.sta.azi, t2.sta.azi) > 0 ? (t1, t2) : (t2, t1)

"""
    copy_trace!(to::T, from::T) where {T<:Seis.Trace}

Copy the `Seis.Trace` trace `a` into `b` without allocating.
"""
@inline function copy_trace!(to::T, from::T) where {T<:Seis.Trace}
    to.t .= from.t
    nothing
end

"""
    minloc2(a) -> imin, jmin

Return an 2-tuple which contains the location `(imin,jmin)` of the minimum value of a
2-dimensional array.
"""
function minloc2(A)
    imin = jmin = 0
    min = zero(eltype(A))
    for j in 1:size(A, 2)
        for i in 1:size(A, 1)
            if i == j == 1
                imin = jmin = 1
                min = A[i,j]
            else
                if A[i,j] < min
                    imin = i
                    jmin = j
                    min = A[i,j]
                end
            end
        end
    end
    imin, jmin
end

"""
    maxabsloc2(A) -> imax, jmax

Return the indices of a two-dimensional array `A` where its maximum
absolute value occurs.
"""
function maxabsloc2(A)
    imin = jmin = 0
    maxval = zero(eltype(A))
    for j in 1:size(A, 2)
        for i in 1:size(A, 1)
            val = abs(A[i,j])
            if i == j == 1
                imin = jmin = 1
                maxval = val
            else
                if val > maxval
                    imin = i
                    jmin = j
                    maxval = val
                end
            end
        end
    end
    imin, jmin
end

"""
    rotate_traces!(s1, s2, phi)

Rotate the pair of traces `s1` and `s2` clockwise through `phi`°.  Particle motion
appears to move anticlockwise.
"""
function rotate_traces!(s1, s2, phi)
    # Minimal version of Seis.rotate_through! which does not update headers
    sinp, cosp = sincos(deg2rad(phi))
    @inbounds for i in 1:Seis.nsamples(s1)
        s1.t[i], s2.t[i] = cosp*s1.t[i] - sinp*s2.t[i], sinp*s1.t[i] + cosp*s2.t[i]
    end
end

"""
    shift_trace!(t, dt) -> t

Shift a trace `t` backwards in time by `dt` seconds.
"""
function shift_trace!(t, dt)
    T = eltype(t.t)
    n = round(Int, dt/t.delta)
    if n == 0
        return t
    end
    array_shift!(t.t, n)
end

max_stars_repo_name: anowacki/SeisSplit

max_issues_repo_name: anowacki/SeisSplit

max_forks_repo_name: anowacki/SeisSplit


text:
@testset "PSNR" begin
    @info "test: PSNR"

    iqi = PSNR()
    sz_img_3 = (3, 3, 3)

    # Gray image
    type_list = generate_test_types([Bool, Float32, N0f8], [Gray])
    A = [1.0 1.0 1.0; 1.0 1.0 1.0; 0.0 0.0 0.0]
    B = [1.0 1.0 1.0; 0.0 0.0 0.0; 1.0 1.0 1.0]
    for T in type_list
        test_ndarray(iqi, sz_img_3, T)

        a = A .|> T
        b = B .|> T

        # scalar output
        @test assess_psnr(a, b) == assess(PSNR(), a, b)
        @test assess_psnr(a, b) == assess_psnr(a, b, 1.0)
        @test isinf(assess_psnr(a, a))

        # vector output
        @test size(assess_psnr(a, b, [1.0])) == (1,)
        @test mean(assess_psnr(a, b, [1.0])) == assess_psnr(a, b, 1.0)
        @test all(isinf.(assess_psnr(a, a, [1.0])))

        # FIXME: the result of Bool type is not strictly equal to others
        eltype(T) <: Bool && continue
        test_numeric(iqi, a, b, T)
        test_numeric(iqi, channelview(a), channelview(b), T)
    end
    test_cross_type(iqi, A, B, type_list)

    # RGB image
    type_list = generate_test_types([Float32, N0f8], [RGB])
    A = [RGB(0.0, 0.0, 0.0) RGB(0.0, 1.0, 0.0) RGB(0.0, 1.0, 1.0)
        RGB(0.0, 0.0, 1.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 1.0, 0.0)
        RGB(1.0, 1.0, 1.0) RGB(1.0, 0.0, 1.0) RGB(0.0, 0.0, 0.0)]
    B = [RGB(0.0, 0.0, 0.0) RGB(0.0, 0.0, 1.0) RGB(1.0, 1.0, 1.0)
        RGB(0.0, 1.0, 0.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 0.0, 1.0)
        RGB(0.0, 1.0, 1.0) RGB(1.0, 1.0, 0.0) RGB(0.0, 0.0, 0.0)]
    for T in type_list
        test_ndarray(iqi, sz_img_3, T)

        a = A .|> T
        b = B .|> T

        # scalar output
        @test assess_psnr(a, b) == assess(PSNR(), a, b) == PSNR()(a, b)
        @test assess_psnr(a, b) == assess_psnr(a, b, 1.0) == PSNR()(a, b, 1.0)
        @test assess_psnr(a, b) == assess_psnr(channelview(a), channelview(b))
        @test isinf(assess_psnr(a, a))

        # vector output
        @test_throws ArgumentError assess_psnr(a, b, [1.0])
        @test assess_psnr(a, b, [1.0, 1.0, 1.0]) == assess(PSNR(), a, b, [1.0, 1.0, 1.0]) == PSNR()(a, b, [1.0, 1.0, 1.0])
        @test size(assess_psnr(a, b, [1.0, 1.0, 1.0])) == (3,)
        @test mean(assess_psnr(a, b, [1.0, 1.0, 1.0])) != assess_psnr(a, b) # generally they doesn't equal
        @test all(isinf.(assess_psnr(a, a, [1.0, 1.0, 1.0])))

        test_numeric(iqi, a, b, T)
        test_numeric(iqi, channelview(a), channelview(b), T; filename="references/PSNR_2d_RGB")
    end
    type_list = generate_test_types([Float32, N0f8], [RGB, BGR])
    test_cross_type(iqi, A, B, type_list)
    @test isapprox(assess_psnr(Lab.(A), B), assess_psnr(A, B); rtol=1e-5)

    # general Color3 images that doesn't have peakval inferred
    type_list = generate_test_types([Float32], [Lab, HSV])
    A = [RGB(0.0, 0.0, 0.0) RGB(0.0, 1.0, 0.0) RGB(0.0, 1.0, 1.0)
        RGB(0.0, 0.0, 1.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 1.0, 0.0)
        RGB(1.0, 1.0, 1.0) RGB(1.0, 0.0, 1.0) RGB(0.0, 0.0, 0.0)]
    B = [RGB(0.0, 0.0, 0.0) RGB(0.0, 0.0, 1.0) RGB(1.0, 1.0, 1.0)
        RGB(0.0, 1.0, 0.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 0.0, 1.0)
        RGB(0.0, 1.0, 1.0) RGB(1.0, 1.0, 0.0) RGB(0.0, 0.0, 0.0)]
    for T in type_list
        a = A .|> T
        b = B .|> T

        @test assess_psnr(a, b, [1.0, 1.0, 1.0]) == assess(PSNR(), a, b, [1.0, 1.0, 1.0]) == PSNR()(a, b, [1.0, 1.0, 1.0])
        @test all(isinf.(assess_psnr(A, A, [1.0, 1.0, 1.0])))
    end
end


text:
using MPAlign

n = 200
λ = 3
s = 0.8
PG = Pair_ER(n,λ,s)

@time mp_fin = run_bp(PG,15,verbose=true)

M_loglr = create_matrix_lr(PG,mp_fin)

ov1, ov2, v1,v2 = eval_M(PG, M_loglr)
println(ov1," | ",ov2)


text:
import Base.length
import Base.getindex

immutable MinusPlusOneVector
	field  :: Vector{Float64}
	numPlus :: Int64
	function MinusPlusOneVector{T}(inputfield::Vector{T})
		floatfield=ones(size(inputfield))
		numPlus=0
		for i = 1 : Base.size(inputfield)[1]
			if inputfield[i]==one(T)
				numPlus+=1
			elseif inputfield[i]==-one(T)
				floatfield[i] = -1
			else
				throw(ArgumentError("Not 1 or -1"))
			end
		end
		new(floatfield,numPlus)
	end
	function MinusPlusOneVector(floatfield::Vector{Float64})
		numPlus=0
		for i in floatfield
			if i == one(Float64)
				numPlus+=1
			elseif i == -one(Float64)
			else
				throw(ArgumentError("Not 1 or -1"))
			end
		end
		new(floatfield,numPlus)
	end
end

function MinusPlusOneVector{T}(fieldOfT::Vector{T},onesCollection::Set{T})
	floatfield=ones(size(fieldOfT))
	for i in 1 : size(fieldOfT)[1]
		if ! in(fieldOfT[i],onesCollection)
			floatfield[i]=-1
		end
	end
	MinusPlusOneVector(floatfield)
end

# These are used because they get imported by default !
function getindex(a::MinusPlusOneVector, b::UnitRange{Int64})
	getindex(a.field,b)
end
function getindex(a::MinusPlusOneVector, b::Int64)
	getindex(a.field,b)
end
function getindex(a::MinusPlusOneVector, b::Array{Int64,1})
	getindex(a.field,b)
end

function length(a::MinusPlusOneVector)
	length(a.field)
end


text:
module Runtests

using ModiaMedia
#using  Test

#@testset "Test ModiaMedia" begin
   Medium = getMedium("ConstantPropertyLiquidWater")

   println("\n... Medium = ", Medium)

   state  = setState_pT(Medium, 1e5, 300.0)
   d      = density(Medium,state)
   h      = specificEnthalpy(Medium,state)

   println("\nd = ", d)
   println("h = ", h)
   ModiaMedia.standardPlot(Medium; figure=1)

   Medium2 = getMedium("SimpleAir")

   println("\n... Medium2 = ", Medium2)

   state2  = setState_pT(Medium2, 1e5, 300.0)
   d2      = density(Medium2,state2)
   h2      = specificEnthalpy(Medium2,state2)

   println("\nd = ", d2)
   println("h   = ", h2)
   ModiaMedia.standardPlot(Medium2; figure=2)

   Medium3 = getMedium("N2")

   println("\n... Medium3 = ", Medium3)

   state3 = setState_pT(Medium3, 1e5, 300.0)
   d3     = density(Medium3,state3)
   h3     = specificEnthalpy(Medium3,state3)

   println("\nd = ", d3)
   println("h = ", h3)
   ModiaMedia.standardPlot(Medium3; figure=3)

#end

end


text:
module TimeSeriesInterface

import Base: +, -, *

using Dates
using RecipesBase
using Statistics

include("timeseries.jl")
include("forecast.jl")
include("models.jl")
include("FileFormats/FileFormats.jl")

end


text:
isdefined(Base, :__precompile__) && __precompile__()
@info "PowerGraphics.jl loads Plots.jl. Precompile might take a while"
module PowerGraphics

export plot_demand
export plot_dataframe
export plot_pgdata
export plot_results
export plot_fuel
export plot_demand!
export plot_dataframe!
export plot_pgdata!
export plot_results!
export plot_fuel!
export make_fuel_dictionary
export report
export get_generation_data
export get_load_data
export get_service_data
export categorize_data

#I/O Imports
import Dates
import TimeSeries
import Reexport
import Requires
import Colors
import DataFrames
import YAML
Reexport.@reexport using Plots
import DataStructures: OrderedDict, SortedDict
import PowerSystems
import InfrastructureSystems
import PowerSimulations
import InteractiveUtils

const PSY = PowerSystems
const IS = InfrastructureSystems
const PSI = PowerSimulations

include("definitions.jl")
include("plot_data.jl")
include("fuel_results.jl")
include("plot_recipes.jl")
include("plotly_recipes.jl")
include("make_report.jl")
include("call_plots.jl")

function __init__()
    Requires.@require Weave = "44d3d7a6-8a23-5bf8-98c5-b353f8df5ec9" include(
        "make_report.jl",
    )
end

end #module

max_stars_repo_name: NREL-SIIP/InfrastructureSystemsPlots

max_issues_repo_name: NREL-SIIP/InfrastructureSystemsPlots

max_forks_repo_name: NREL-SIIP/InfrastructureSystemsPlots


text:
export GitPatch, lines, nlines, DiffHunk, hunks, nchanges, delta

type GitPatch
    ptr::Ptr{Void}

    function GitPatch(ptr::Ptr{Void})
        @assert ptr != C_NULL
        p = new(ptr)
        finalizer(p, free!)
        return p
    end
end

function free!(p::GitPatch)
    if p.ptr != C_NULL
        ccall((:git_patch_free, libgit2), Void, (Ptr{Void},), p.ptr)
        p.ptr = C_NULL
    end
end

Base.convert(::Type{Ptr{Void}}, p::GitPatch) = p.ptr

Base.diff(repo::GitRepo, blob::GitBlob, other::Nothing, opts=nothing) = begin
    old_path_ptr, new_path_ptr = zero(Ptr{Uint8}), zero(Ptr{Uint8})
    if opts != nothing
        if get(opts, :old_path, nothing) != nothing
            old_path = convert(Ptr{Uint8}, opts[:old_path]::String)
        end
        if get(opts, :new_path, nothing) != nothing
            new_path = convert(Ptr{Uint8}, opts[:new_path]::String)
        end
    end
    gopts = parse_git_diff_options(opts)
    patch_ptr = Ptr{Void}[0]
    @check ccall((:git_patch_from_blobs, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Uint8}, Ptr{Void}, Ptr{Uint8}, Ptr{DiffOptionsStruct}),
                 patch_ptr, blob, old_path_ptr, C_NULL, new_path_ptr, &gopts)
    return GitPatch(patch_ptr[1])
end

Base.diff(repo::GitRepo, blob::GitBlob, other::GitBlob, opts::MaybeDict=nothing) = begin
    old_path_ptr, new_path_ptr = zero(Ptr{Uint8}), zero(Ptr{Uint8})
    if opts != nothing
        if get(opts, :old_path, nothing) != nothing
            old_path_ptr = convert(Ptr{Uint8}, opts[:old_path]::String)
        end
        if get(opts, :new_path, nothing) != nothing
            new_path_ptr = convert(Ptr{Uint8}, opts[:new_path]::String)
        end
    end
    gopts = parse_git_diff_options(opts)
    patch_ptr = Ptr{Void}[0]
    @check ccall((:git_patch_from_blobs, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Uint8}, Ptr{Void}, Ptr{Uint8}, Ptr{DiffOptionsStruct}),
                 patch_ptr, blob, old_path_ptr, other, new_path_ptr, &gopts)
    return GitPatch(patch_ptr[1])
end

Base.diff(repo::GitRepo, blob::GitBlob, other::String, opts::MaybeDict=nothing) = begin
    old_path_ptr, new_path_ptr = zero(Ptr{Uint8}), zero(Ptr{Uint8})
    if opts != nothing
        if get(opts, :old_path, nothing) != nothing
            old_path_ptr = convert(Ptr{Uint8}, opts[:old_path]::String)
        end
        if get(opts, :new_path, nothing) != nothing
            new_path_ptr = convert(Ptr{Uint8}, opts[:new_path]::String)
        end
    end
    gopts = parse_git_diff_options(opts)
    patch_ptr = Ptr{Void}[0]
    @check ccall((:git_patch_from_blob_and_buffer, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Uint8},
                  Ptr{Uint8}, Csize_t, Ptr{Uint8}, Ptr{DiffOptionsStruct}),
                 patch_ptr, blob, old_path_ptr, buffer, length(buffer), new_path_ptr, &gopts)
    return GitPatch(patch_ptr[1])
end

type GitPatchStat
    adds::Int
    dels::Int
end

Base.stat(p::GitPatch) = begin
    adds = Csize_t[0]
    dels = Csize_t[0]
    @check ccall((:git_patch_line_stats, libgit2), Cint,
                 (Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Void}),
                 C_NULL, adds, dels, p)
    return GitPatchStat(int(adds[1]), int(dels[1]))
end

nchanges(p::GitPatch) = (s = stat(p); s.adds + s.dels)

type DiffHunk
    patch::GitPatch
    header::UTF8String
    line_count::Int
    hunk_index::Int
    old_start::Int
    old_lines::Int
    new_start::Int
    new_lines::Int

    function DiffHunk(p::GitPatch, ptr::Ptr{DiffHunkStruct}, idx::Integer, lc::Integer)
        @assert ptr != C_NULL
        h = unsafe_load(ptr)::DiffHunkStruct
        head_arr = zeros(Uint8, 128)
        #TODO: get rid of this ugly hack
        for i=1:128
            head_arr[i] = getfield(h.header, symbol("c$i"))
        end
        return new(p,
                   utf8(bytestring(convert(Ptr{Uint8}, head_arr))),
                   lc,
                   idx,
                   h.old_start,
                   h.old_lines,
                   h.new_start,
                   h.new_lines)
    end
end

function hunks(p::GitPatch)
    nhunks = ccall((:git_patch_num_hunks, libgit2), Csize_t, (Ptr{Void},), p)
    if nhunks == 0
        return nothing
    end
    err = zero(Cint)
    hunk_ptr  = Ptr{DiffHunkStruct}[0]
    lines_ptr = Csize_t[0]
    hs = DiffHunk[]
    for i in 1:nhunks
        err = ccall((:git_patch_get_hunk, libgit2), Cint,
                    (Ptr{Ptr{DiffHunkStruct}}, Ptr{Csize_t}, Ptr{Void}, Csize_t),
                    hunk_ptr, lines_ptr, p, i-1)
        err == GitErrorConst.GIT_OK || break
        push!(hs, DiffHunk(p, hunk_ptr[1], i, lines_ptr[1]))
    end
    if err != GitErrorConst.GIT_OK
        throw(GitError(err))
    end
    return hs
end

function line_origin_to_symbol(o)
    o == GitConst.DIFF_LINE_CONTEXT   && return :context
    o == GitConst.DIFF_LINE_ADDITION  && return :addition
    o == GitConst.DIFF_LINE_DELETION  && return :deletion
    o == GitConst.DIFF_LINE_ADD_EOFNL && return :eof_newline_added
    o == GitConst.DIFF_LINE_DEL_EOFNL && return :eof_newline_removed
    o == GitConst.DIFF_LINE_CONTEXT_EOFNL && return :eof_no_newline
    throw(ArgumentError("Unknown line origin constant $o"))
end

type DiffLine
    hunk::DiffHunk
    line_origin::Symbol
    content::ByteString
    old_lineno::Int
    new_lineno::Int
    content_offset::Union(Nothing, Int)

    function DiffLine(h::DiffHunk, ptr::Ptr{DiffLineStruct})
        @assert ptr != C_NULL
        l = unsafe_load(ptr)::DiffLineStruct
        c = Array(Uint8, l.content_len)
        for i in 1:l.content_len
            c[i] = unsafe_load(l.content, i)
        end
        return new(h,
                   line_origin_to_symbol(l.origin),
                   bytestring(convert(Ptr{Uint8}, c), l.content_len),
                   l.old_lineno,
                   l.new_lineno,
                   l.content_offset == -1 ? nothing : l.content_offset)
    end
end

function lines(h::DiffHunk)
    err::Cint = 0
    ls = DiffLine[]
    lc = h.line_count
    hi = h.hunk_index
    line_ptr = Ptr{DiffLineStruct}[0]
    for i in 1:lc
        err = ccall((:git_patch_get_line_in_hunk, libgit2), Cint,
                    (Ptr{Ptr{DiffLineStruct}}, Ptr{Void}, Csize_t, Csize_t),
                    line_ptr, h.patch, hi-1, i-1)
        err == GitErrorConst.GIT_OK || break
        push!(ls, DiffLine(h, line_ptr[1]))
    end
    if err != GitErrorConst.GIT_OK
        throw(GitError(err))
    end
    return ls
end

function cb_patch_print(delta_ptr::Ptr{Void}, hunk_ptr::Ptr{Void},
                        line_ptr::Ptr{DiffLineStruct}, payload::Ptr{Void})
    l = unsafe_load(line_ptr)::DiffLineStruct
    s = unsafe_pointer_to_objref(payload)::Vector{Uint8}
    add_origin = false
    if l.origin == GitConst.DIFF_LINE_CONTEXT ||
       l.origin == GitConst.DIFF_LINE_ADDITION ||
       l.origin == GitConst.DIFF_LINE_DELETION
       add_origin = true
    end
    prev_len = length(s)
    if add_origin
        resize!(s, prev_len + l.content_len + 1)
        s[prev_len + 1] = l.origin
        for i in 1:l.content_len
            s[prev_len + i + 1] = unsafe_load(l.content, i)
        end
    else
        resize!(s, prev_len + l.content_len)
        for i in 1:l.content_len
            s[prev_len + i] = unsafe_load(l.content, i)
        end
    end
    return GitErrorConst.GIT_OK
end

const c_cb_patch_print = cfunction(cb_patch_print, Cint,
                                  (Ptr{Void}, Ptr{Void}, Ptr{DiffLineStruct}, Ptr{Void}))

Base.string(p::GitPatch) = begin
    s = Uint8[]
    @check ccall((:git_patch_print, libgit2), Cint,
                 (Ptr{Void}, Ptr{Void}, Any), p, c_cb_patch_print, &s)
    return utf8(bytestring(convert(Ptr{Uint8}, s)))
end

nhunks(p::GitPatch) = int(ccall((:git_patch_num_hunks, libgit2), Cint, (Ptr{Void},), p))

function delta(p::GitPatch)
    delta_ptr = ccall((:git_patch_get_delta, libgit2), Ptr{DiffDeltaStruct}, (Ptr{Void},), p)
    return DiffDelta(delta_ptr)
end

function line_stats(p::GitPatch)
    lines = Csize_t[0]
    additions = Csize_t[0]
    deletions = Csize_t[0]
    @check ccall((:git_patch_line_stats, libgit2), Cint,
                 (Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Csize_t}, Ptr{Void}),
                 lines, additions, deletions, p)
    return (int(lines[1]), int(additions[1]), int(deletions[1]))
end

nlines(p::GitPatch) = reduce(+, 0, line_stats(p))


text:
export DEBUG_METADATA_VERSION, strip_debuginfo!

DEBUG_METADATA_VERSION() = API.LLVMDebugMetadataVersion()

strip_debuginfo!(mod::Module) = API.LLVMStripModuleDebugInfo(mod)

if version() >= v"8.0"
set_subprogram!(func::Function, sp::Metadata) = API.LLVMSetSubprogram(func, sp)
get_subprogram(func::Function) = Metadata(API.LLVMGetSubprogram(func))
end


text:
function div(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:IEEEFloat}
    (!isfinite(x) || isnan(y)) && return nan(DoubleFloat{T})
    !isfinite(y) && return zero(DoubleFloat{T})
    return trunc(x / y)
end

function fld(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:IEEEFloat}
    (!isfinite(x) || isnan(y)) && return nan(DoubleFloat{T})
    !isfinite(y) && return zero(DoubleFloat{T})
    return floor(x / y)
end

function cld(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:IEEEFloat}
    (!isfinite(x) || isnan(y)) && return nan(DoubleFloat{T})
    !isfinite(y) && return zero(DoubleFloat{T})
    return ceil(x / y)
end

function rem(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:IEEEFloat}
    (!isfinite(x) || isnan(y)) && return nan(DoubleFloat{T})
    !isfinite(y) && return x
    return x - div(x,y) * y
end

function mod(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:IEEEFloat}
    (!isfinite(x) || isnan(y)) && return nan(DoubleFloat{T})
    !isfinite(y) && return x
    return x - fld(x,y) * y
end

function divrem(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:IEEEFloat}
    (!isfinite(x) || isnan(y)) && return nan(DoubleFloat{T})
    !isfinite(y) && return zero(DoubleFloat{T})
    dv = trunc(x / y)
    rm = x - dv * y
    return dv, rm
end

function fldmod(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:IEEEFloat}
    (!isfinite(x) || isnan(y)) && return nan(DoubleFloat{T})
    !isfinite(y) && return zero(DoubleFloat{T})
    fr = floor(x / y)
    md = x - fr * y
    return fr, md
end


text:
import Reactive: set_value!

Base.getindex(x::Node) = value(x)
function Base.setindex!(x::Node{T}, value) where T
    set_value!(x, convert(T, value))
    force_update!()
    push!(x, value)
end

# I don't want to use map anymore, it's so ambigious, especially to newcomers.
# TODO should this become it's own function?
const lift = Reactive.map

to_value(x) = value(x)
to_value(x::Nothing) = x

to_node(::Type{T1}, x::Node{T2}, name = :node) where {T1, T2} = signal_convert(Node{T1}, x, name)
to_node(x::T, name = :node) where T = to_node(T, x)
to_node(::Type{T}, x, name = :node) where T = to_node(T, Node(T, x, name = string(name)))

signal_convert(::Type{Signal{T1}}, x::Signal{T1}, name = :node) where T1 = x
signal_convert(::Type{Signal{T1}}, x::Signal{T2}, name = :node) where {T1, T2} = map(x-> convert(T1, x), x, typ = T1, name = string(name))
signal_convert(::Type{Signal{T1}}, x::T2, name = :node) where {T1, T2} = Node(T1, convert(T1, x), name = string(name))
signal_convert(t, x, name = :node) = x

node(name, node) = Node(node, name = string(name))
node(name, node::Node) = map(identity, node, name = string(name))

function close_all_nodes(any::T) where T
    for field in fieldnames(T)
        value = getfield(any, field)
        (value isa Node) && close(value, true)
    end
end

function disconnect!(s::Node)
    unpreserve(s)#; empty!(s.actions)
    s.parents = (); close(s, false)
    return
end

function get_children(x::Signal)
    filter!(x-> x != nothing, map(x-> x.value, Reactive.nodes[Reactive.edges[x.id]]))
end
function children_with(action::AT, signal::Signal, rest::Signal...) where AT
    c = filter!(get_children(signal)) do x
        any(y-> isa(y, AT), x.actions)
    end
    isempty(c) && return false, Signal(nothing)
    node = first(c) # what if more?
    has_node = all(rest) do signal
        node in get_children(signal)
    end
    has_node, node
end

"""
    map_once(closure, inputs::Signal....)::Signal

Like Reactive.foreach, in the sense that it will be preserved even if no reference is kept.
The difference is, that you can call map once multiple times with the same closure and it will
close the old result Signal and register a new one instead.

```
function test(s1::Signal)
    s3 = map_once(x-> (println("1 ", x); x), s1)
    s3 = map_once(x-> (println("2 ", x); x), s1)

end
test(Signal(1), Signal(2))
>

"""
function map_once(
        f, input::Signal, inputsrest::Signal...;
        init = f(map(value, (input,inputsrest...))...),
        typ = typeof(init), name = Reactive.auto_name!("map", input, inputsrest...)
    )
    output = Signal(typ, init, (input, inputsrest...); name = name)
    action_func = function ()
        Reactive.set_value!(output, f(value.((input, inputsrest...))...))
    end
    # The action func should be unique regarding the types of inputs + f!
    # so we can use it to figure out if map_once was already called with this f + input types.
    # At least for closures, the name should be unique once per function body,
    # so exactly the kind of uniqueness we need!
    has, prev_output = children_with(action_func, input, inputsrest...)
    # if all input nodes have the action_func, then `node` must be the result
    # of a previous call to map_once, meaning we're about to replace this node!
    # clean up previously connected signal!
    has && disconnect!(prev_output)
    # add the action to the new output
    Reactive.add_action!(action_func, output)
    preserve(output)
end

#
# call_count = Ref(1)
# io = IOBuffer()
# noyield_println(x...) = println(io, x...)
# function example1(s1)
#     # only works with closure, which will have it's unique,
#     # which stays consistent per call
#     cc = call_count[]
#     s2 = map_once(v-> noyield_println("func 1 call $cc :", v), s1)
#     push!(s1, 22)
#     # this is a unique closure, so it won't disconnect s2
#     s3 = map_once(v-> noyield_println("func 2 call $cc :", v), s1)
#     push!(s1, 33)
#     call_count[] = call_count[] + 1
#     s2, s3
# end
# s1 = Signal(1)
# call_count[] = 1
# s2 = example1(s1)
# yield() # give a chance to run event loop
# println(String(take!(io)))
# func 1 call 1 :1
# func 2 call 1 :1
# func 1 call 1 :22
# func 2 call 1 :22
# func 1 call 1 :33
# func 2 call 1 :33
# s2 = example1(s1)
# func 1 call 2 :33 # first call to maponce
# func 2 call 2 :33 # first call to maponce
# func 2 call 1 :22 # first call to push!, func 2 from call 1 is still registered
# func 1 call 2 :22 # first call to push!, call 2
# func 1 call 2 :33 # second call to push!, call 1 is gone now!
# func 2 call 2 :33
#
# function no_closure(v)
#     println("func 1 ", v)
# end
# function example2(s1)
#     # println is a function which will stay the same
#     s2 = map_once(no_closure, s1)
#     push!(s1, 22)
#     # so it will globally replace any map_once(println, call, s1)
#     # (need to be those exact signals)
#     s3 = map_once(no_closure, s1)
#     push!(s1, 33)
#     s2
# end
#
# s1 = Signal(1)
# s2 = example2(s1)
#
# function example3(s1, call)
#     # normal behaviour doesn't disconnect (with map instead of foreach it will disconnect )
#     # whenever gc thinks it should
#     s2 = foreach(println, call, s1)
#     push!(s1, 22)
#     s2
# end
# s1 = Signal(1)
# c = Signal("call 1 value: ")
# s2 = example1(c, s1)
# call 1 value: 1 # from init 1
# call 1 value: 1 # from init 2
# call 1 value: 22 # from first push
# call 1 value: 33 # from second push call 1
# call 1 value: 33 # from second push call 2
#
# s3 = example1(c, s1)
#     call 1 value: 22
#     call 1 value: 22
#     call 1 value: 22
#     call 1 value: 22
#     call 1 value: 22
#     call 1 value: 22
#
# s1 = Signal(1)
# c  = Signal("call 1, value ")
# s2 = example2(s1, c)
# >call 1, value 1
# >call 1, value 22
# s3 = example2(s1, 2)
# >call 2, value 22
# >call 1, value 22 # old signal still firing
# >call 2, value 22
# # same as above
# s1 = Signal(1)
# c = Signal("call 1, value ")
# s2 = example3(s1, c)
# s3 = example3(s1, c)
# >call 2, value 22
# >call 1, value 22 # old signal still firing
# >call 2, value 22
# get_children(c)[1] == s2
# get_children(s1)[1] == s2
# s2.parents  == (c, s1)
# c.id
# s1.id
# s2.parents[2].id
#
# s1 = Signal(1)
# s2 = Signal("sada")
# s1.id
# s2.id
# s3 = map(string, s1, s2)
# s3.parents == (s1, s2)


text:
module Poincare

using GPUBenchmarks, BenchmarkTools, GPUArrays
import CUDAnative

const cu = CUDAnative
description = """
Poincare section of a chaotic neuronal network.
The domination of OpenCL in this benchmark might be due to a better use of vector intrinsics in Transpiler.jl, but needs some
more investigations.
Result of calculation:
![](https://github.com/JuliaGPU/GPUBenchmarks.jl/blob/master/results/plots/poincare_result.png?raw=true)
"""

# Original poincare implementation by https://github.com/RainerEngelken
# GPU version by Simon Danisch

function poincare_inner{N}(rv, result::CUDAnative.CuDeviceArray, c, π, ::Val{N}, n)
    # find next spiking neuron
    ϕ₁, ϕ₂, ϕ₃ = rv[1], rv[2], rv[3]
    πh = π / 2f0
    π2 = π * 2f0
    for unused = 1:N
        if ϕ₁ > ϕ₂
            if ϕ₁ > ϕ₃
                # first neuron is spiking
                dt = πh - ϕ₁
                # evolve phases till next spike time
                ϕ₁ = -πh
                ϕ₂ = cu.atan(cu.tan(ϕ₂ + dt) - c)
                ϕ₃ += dt
                # save state of neuron 2 and 3
                x = Cuint(cu.max(cu.ceil(((ϕ₂ + πh) / π) * (Float32(n) - 1f0)) + 1f0, 1f0))
                y = Cuint(cu.max(cu.ceil(((ϕ₃ + πh) / π) * (Float32(n) - 1f0)) + 1f0, 1f0))
                i1d = GPUArrays.gpu_sub2ind((n, n), (x, y)) # convert to linear index
                accum = result[i1d]
                # this is unsafe, since it could read + write from different threads, but good enough for the stochastic kind of process we're doing
                result[i1d] = accum + 1f0
                continue
            end
        else
            if ϕ₂ > ϕ₃
                # second neuron is spiking
                dt = πh - ϕ₂
                # evolve phases till next spike time
                ϕ₁ += dt
                ϕ₂ = -πh
                ϕ₃ = cu.atan(cu.tan(ϕ₃ + dt) - c)
                continue
            end
        end
        # third neuron is spikinga
        dt = πh - ϕ₃
        # evolve phases till next spike time
        ϕ₁ += dt
        ϕ₂ = cu.atan(cu.tan(ϕ₂ + dt) - c)
        ϕ₃ = -πh
    end
    return
end

function poincare_inner{N}(rv, result, c, π, ::Val{N}, n)
    # find next spiking neuron
    ϕ₁, ϕ₂, ϕ₃ = rv[1], rv[2], rv[3]
    πh = π / 2f0
    π2 = π * 2f0
    for unused = 1:N
        if ϕ₁ > ϕ₂
            if ϕ₁ > ϕ₃
                # first neuron is spiking
                dt = πh - ϕ₁
                # evolve phases till next spike time
                ϕ₁ = -πh
                ϕ₂ = atan(tan(ϕ₂ + dt) - c)
                ϕ₃ += dt
                # save state of neuron 2 and 3
                x = Cuint(max(round(((ϕ₂ + πh) / π) * (Float32(n) - 1f0)) + 1f0, 1f0))
                y = Cuint(max(round(((ϕ₃ + πh) / π) * (Float32(n) - 1f0)) + 1f0, 1f0))
                i1d = GPUArrays.gpu_sub2ind((n, n), (x, y)) # convert to linear index
                accum = result[i1d]
                # this is unsafe, since it could read + write from different threads, but good enough for the stochastic kind of process we're doing
                result[i1d] = accum + 1f0
                continue
            end
        else
            if ϕ₂ > ϕ₃
                # second neuron is spiking
                dt = πh - ϕ₂
                # evolve phases till next spike time
                ϕ₁ += dt
                ϕ₂ = -πh
                ϕ₃ = atan(tan(ϕ₃ + dt) - c)
                continue
            end
        end
        # third neuron is spikinga
        dt = πh - ϕ₃
        # evolve phases till next spike time
        ϕ₁ += dt
        ϕ₂ = atan(tan(ϕ₂ + dt) - c)
        ϕ₃ = -πh
    end
    return
end

function poincare_inner(n, seeds::GPUArray, result, c, π, val::Val{N}) where N
    foreach(poincare_inner, seeds, result, c, Float32(pi), val, n)
end
function poincare_inner(n, seeds::Array, result, c, π, val::Val{N}) where N
    for rv in seeds
        poincare_inner(rv, result[], c, π, val, n)
    end
end

function makeresult(bench, N, device, hardware, mdiff)
    BenchResult(
        "Poincare",
        bench,
        N,
        Float32,
        string(device),
        hardware,
        @__FILE__,
        mdiff
    )
end

function execute(device)
    results = BenchResult[]
    hardware, AT = init(device)
    c = 1f0; divisor = 2^11
    srand(2)
    for i = 3:9
        N = 10^i
        ND = Cuint(1024)
        result = AT(zeros(Float32, ND, ND))
        _n = div(N, divisor)
        seeds = AT([ntuple(i-> rand(Float32), Val{3}) for x in 1:divisor])
        bench = @benchmark begin
            $(poincare_inner)($ND, $seeds, $(Base.RefValue(result)), $c, $(Float32(pi)), $(Val{_n}()))
            synchronize($(result)) # synchronize for the benchmark
        end
        jl_results = zeros(Float32, ND, ND)
        result = AT(jl_results)
        _n = div(N, divisor)
        jl_seeds = [ntuple(i-> rand(Float32), Val{3}) for x in 1:divisor]
        seeds = AT(jl_seeds)
        poincare_inner(ND, seeds, Base.RefValue(result), c, Float32(pi), Val{_n}())
        poincare_inner(ND, jl_seeds, Base.RefValue(jl_results), c, Float32(pi), Val{_n}())
        mdiff = meandifference(jl_results, result)
        @show mdiff
        push!(results, makeresult(bench, N, device, hardware, mdiff))
    end
    results
end

# using FileIO, Interpolations, Colors, GPUBenchmarks, GPUArrays, ColorVectorSpace, FixedPointNumbers
# device = :opencl
# hardware, AT = init(device)
# N = 10^9
# ND = Cuint(512)
# result = AT(zeros(Float32, ND, ND))
# divisor = 2^11
# _n = div(N, divisor)
# seeds = AT([ntuple(i-> rand(Float32), Val{3}) for x in 1:divisor])
# _n = div(N, divisor)
# c = 1f0
# seeds = AT([ntuple(i-> rand(Float32), Val{3}) for x in 1:divisor])
# poincare_inner(ND, seeds, Base.RefValue(result), c, Float32(pi), Val{_n}())
# synchronize(result)
# cn = 100
# cmap = interpolate(colormap("Oranges", cn), BSpline(Linear()), OnCell());
# img_color = map(Array(result)) do val
#     val = val / 2000f0
#     val = clamp(val, 0f0, 1f0);
#     idx = (val * (cn - 1)) + 1.0
#     RGB{N0f8}(cmap[idx])
# end
# #save as an image
# save(GPUBenchmarks.dir("results", "plots", "poincare_result.png"), img_color)

end


text:
"""
Example modelling a linear beam element

Formulation:

K=-E*A*[dw_i!dx*du_j!dx]
"""

# local coordinates:
local_interp=InterpolationFunction(
    :(a0+a1*xi), # interpolation function expression
    [:xi], # local coordinate system degrees of freedom
    [:a0, :a1], # linear coefficients for regression
    [
        [0] # conditions to impose with regression: absolute value (order 0 derivative in the first axis)
    ];
    maximum_derivative=1
)

# absolute coordinates
abs_interp=InterpolationFunction(
    :(a0+a1*x), # interpolation function expression
    [:x], # absolute coordinate system degrees of freedom
    [:a0, :a1], # linear coefficients for regression
    [
        [0] # conditions to impose with regression: absolute value (order 0 derivative)
    ];
    maximum_derivative=1
)

# defining an integration domain:
dmn=Domain(
    local_interp, # local coordinate system
    [
        [0.0, 1.0] # boundaries for local coordinates: xi between 0 and 1
    ],
    [
        [0.0], # positions at which to impose values for interpolation:
        [1.0] # xi=0, xi=1
    ]
)

# defining a residual for within the beam:
domain_resd=Residual(
    [
        ( # variable w (weight), interpolated by InterpolationFunction abs_interp,
            :w, abs_interp, [
                (:dw!dx, [1]) # accompanied by its first order derivative
            ]
        )
        ( # variable u (displacement), interpolated by InterpolationFunction abs_interp,
            :u, abs_interp, [
                (:du!dx, [1]) # accompanied by its first order derivative
            ]
        )
    ],
    :(-E*(dw!dx*du!dx)),
    3; # Gauss-Legendre quadrature order
    extra_args=[:E, :A] # additional argument: elasticity
)

dmn_resd_func=get_volume_residual_function(domain_resd)

# defining a specific finite element:
fin=Finel(
    [
        [0.0],
        [1.0] # defined by corner points at x=0 and x=1 (corresponding to xi=0 and xi=1 in local coordinates: see Domain)
    ],
    dmn
)

# obtaining elasticity matrix:
K=zeros(Float64, 2, 2)

ws=zeros(Float64, 2)
us=zeros(Float64, 2)

E=1.0
A=1.0

for i=1:2
    ws[i]=1.0

    for j=1:2
        us[j]=1.0

        variable_vals=[ws, us]
        K[i, j]=dmn_resd_func(variable_vals, fin, E, A)

        us[j]=0.0
    end

    ws[i]=0.0
end

@assert isapprox(
    K,
    [
        -1.0 1.0;
        1.0 -1.0
    ]
)


text:
using Test
using HMCResearchRandomWalks.Arlett

@testset "Nontrivial PARAMETERS_arlett.jl" begin
    @test isapprox(Arlett.SENSOR_CENTER_MAX_X, 12.5)
    @test isapprox(Arlett.SENSOR_INNER_ADJ_MIN_X, 32.5)
    @test isapprox(Arlett.SENSOR_INNER_ADJ_MAX_X, 57.5)
    @test isapprox(Arlett.SENSOR_OUTER_ADJ_MIN_X, 77.5)
    @test isapprox(Arlett.SENSOR_OUTER_ADJ_MAX_X, 102.5)
    @test isapprox(Arlett.ENZYME_MAX_Y, 1.65)
    @test isapprox(Arlett.SPAWN_LEFT_X, -13.499)
    @test isapprox(Arlett.SPAWN_RIGHT_X, 13.499)
    @test isapprox(Arlett.SPAWN_ENZYME_MAX_Y, 1.649)
    @test isapprox(Arlett.MARGIN_OF_COLLISION, 0.2)
    @test isapprox(Arlett.ESCAPE_X, 122.5)
    @test isapprox(Arlett.ESCAPE_Y, 122.5)
    @test isapprox(Arlett.SAFE_MIN_Y, 1.85)
    @test isapprox(Arlett.SAFE_MAX_Y, 122.3)
    @test isapprox(Arlett.SAFE_MAX_X, 122.3)
    #TODO: Enzyme-Water Corners
end

@testset "spawnrandompoint() in Arlett" begin
    leftwallspawn = Arlett.spawnrandompoint(0.0)
    @test isapprox(leftwallspawn[1], -13.499)
    @test isapprox(leftwallspawn[2], 1.501)

    leftcornerspawn = Arlett.spawnrandompoint(0.149/27.296)
    @test isapprox(leftcornerspawn[1], -13.499)
    @test isapprox(leftcornerspawn[2], 1.649)

    rightwallspawn = Arlett.spawnrandompoint(1.0)
    @test isapprox(rightwallspawn[1], 13.499)
    @test isapprox(rightwallspawn[2], 1.501)

    rightcornerspawn = Arlett.spawnrandompoint(27.147/27.296)
    @test isapprox(rightcornerspawn[1], 13.499)
    @test isapprox(leftcornerspawn[2], 1.649)

    centerspawn = Arlett.spawnrandompoint(0.5)
    @test isapprox(centerspawn[1], 0)
    @test isapprox(centerspawn[2], 1.649)

    @test Arlett.inenz(leftwallspawn)
    @test Arlett.inenz(leftcornerspawn)
    @test Arlett.inenz(rightcornerspawn)
    @test Arlett.inenz(rightwallspawn)
    @test Arlett.inenz(centerspawn)
    @test !Arlett.inwalls(leftwallspawn[1], leftwallspawn[2])
    @test !Arlett.inwalls(rightwallspawn[1], rightwallspawn[2])

    randspawn = Arlett.spawnrandompoint()
    @test Arlett.inenz(randspawn)
    @test !Arlett.inwalls(randspawn[1], randspawn[2])
    @test isapprox(abs(randspawn[1]), 13.499) || isapprox(randspawn[2], 1.649)
end

@testset "whereoutsidespawn in Arlett" begin
    @test Arlett.whereoutsidespawn([-14, 1.6]) == "W"
    @test Arlett.whereoutsidespawn([-14, 1.7]) == "NW"
    @test Arlett.whereoutsidespawn([0, 1.7]) == "N"
    @test Arlett.whereoutsidespawn([14, 1.7]) == "NE"
    @test Arlett.whereoutsidespawn([14, 1.6]) == "E"
    @test Arlett.whereoutsidespawn([0, 1.6]) == "whereoutsidespawn Error"
end

"testing functions for locations"
function rand_x() return Arlett.ESCAPE_X * rand() * 2 - Arlett.ESCAPE_X end
function rand_y() return Arlett.ESCAPE_Y * rand() - 0.5 * Arlett.ESCAPE_Y end

@testset "flow_arlett.jl" begin
    set_FLOW_BIAS(true)
    coordinates_without_flow = [
        [-90, 0.1], [-45, 0.1], [0, 0.1], [45, 0.1], [90, 0.1], # in well, in ppd layer
        [-90, 1.4], [-45, 1.4], [45, 1.4], [90, 1.4],           # in well, in water
        [0, 1.4], [0, 1.6]                                      # in well, in enzymatic layer
    ]
    for coordinate in coordinates_without_flow
        @test Arlett.flow_arlett(coordinate) == 0
    end
    y_vals = [
        1.7,
        5,
        25,
        50,
        100
    ]
    correct_flow = [
        4.58262333 * 10^-5
        1.34607810 * 10^-4
        6.67729077 * 10^-4
        1.32218322 * 10^-3
        2.59126672 * 10^-3
    ]
    correct_speed = [
        16.92441308
        49.71297043
        246.60378821
        488.30491661
        956.999193998
    ]
    for i in 1:5
        x = rand_x() # flow is x invariant
        @test isapprox(Arlett.flow_arlett([x, y_vals[i]]), correct_flow[i])
        @test isapprox(Arlett.getspeed(y_vals[i]), correct_speed[i])
    end
end

@testset "locations_arlett.jl" begin
    @test !Arlett.inwalls(rand_x(), 1.51)

    @test Arlett.inwater([rand_x(), 1.66])
    @test !Arlett.inwater([0.0, 1.65])

    @test Arlett.inenz([0.0, 0.151])
    @test Arlett.inenz([0.0, 1.65])
    @test Arlett.inenz([-13.4, 1.6])

    @test !Arlett.inoverflowenz(1.0, 0.15)
    @test !Arlett.inoverflowenz(1.0, 1.5)
    @test Arlett.inoverflowenz(1.0, 1.6)
    @test Arlett.inoverflowenz(13.4, 1.6)
    @test !Arlett.inoverflowenz(1.0, 1.66)
    @test !Arlett.inoverflowenz(-13.6, 1.6)

    @test Arlett.incentersensorx(12.49*rand())
    @test Arlett.ininnersensorx(12.49*rand() + 32.5)
    @test Arlett.ininnersensorx(12.49*rand() - 57.5)
    @test Arlett.inoutersensorx(12.49*rand() + 77.5)
    @test Arlett.inoutersensorx(12.49*rand() - 102.5)

    @test Arlett.insensor([rand_x(), -1.0])
    @test Arlett.insensor([rand_x(), 0.0])
    @test !Arlett.insensor([rand_x(), 0.1])

    for x in [-90.0, -45.0, 0.0, 45.0, 90.0]
        @test !Arlett.inwalls(x, rand_y())
        @test Arlett.inwalls(x + 14, 1.5)
        @test !Arlett.inwalls(x + 14, 1.51)

        if x != 0.0
            @test Arlett.inwater([x, 0.151])
            @test Arlett.inwater([x, 1.65])
            @test !Arlett.inenz([x, 0.151])
            @test !Arlett.inenz([x, 1.65])
        end
        @test !Arlett.inwater([x, 0.15])
        @test Arlett.inppd([x, 0.15])
        @test !Arlett.inppd([x, 0.151])
    end
end

@testset "sensorcases() in locations_arlett.jl" begin
    @test Arlett.sensorcases(-90.0) == "left outer sensor"
    @test Arlett.sensorcases(-45.0) == "left inner sensor"
    @test Arlett.sensorcases(0.0) == "center sensor"
    @test Arlett.sensorcases(45.0) == "right inner sensor"
    @test Arlett.sensorcases(90.0) == "right outer sensor"
end

# @testset "wallcases() in locations_arlett.jl" begin
#
# end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
using Revise
using UnivariateGaussianProcess

# An example usage of the package.
function run_gp()
    obj_func = sin
    bounds = (-π, π)
    cov_func = (x1,x2) -> matern1(x1, x2; σ=2)
    acq_func = ucb

    noise = 1e-2
    init_mean = (x) -> x
    x_points = [-π, 0., π]
    y_points = obj_func.(x_points)
    points = DATA_POINTS(x_points, y_points)

    return gaussian_process(
        obj_func,
        bounds,
        cov_func,
        acq_func;
        noise,
        init_mean,
        evals = 6,
        points,
        plot = true,
        plot_only_final = false,
    )
end


text:
using DashBootstrapComponents

tab1_content = dbc_card(
    dbc_cardbody([
        html_p("This is tab 1!", className = "card-text"),
        dbc_button("Click here", color = "success"),
    ]),
    className = "mt-3",
);

tab2_content = dbc_card(
    dbc_cardbody([
        html_p("This is tab 2!", className = "card-text"),
        dbc_button("Don't click here", color = "danger"),
    ]),
    className = "mt-3",
);

tabs = dbc_tabs([
    dbc_tab(tab1_content, label = "Tab 1"),
    dbc_tab(tab2_content, label = "Tab 2"),
    dbc_tab("This tab's content is never seen", label = "Tab 3", disabled = true),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
### Credit to djsegal
### https://discourse.julialang.org/t/is-there-a-faster-bisection-root-solver-that-uses-atol/12658/21

# f, cur_a, cur_b, cur_f_a, cur_f_b  = ap, x1, x0, fx1, fx0;
# cur_flag, abstol = true, cbrt(eps())
# cur_c, cur_f_c, cur_d, cur_f_d, bad_streak = NaN,NaN,NaN,NaN,0

@inline check_approx(x, y, atol, rtol) = norm(x-y) <= atol + rtol*max(norm(x), norm(y))
@inline check_approx(x, y, atol) = norm(x-y) <= atol

function custom_bisection(f, cur_a::T, cur_b::T, cur_f_a::T, cur_f_b::T, abstol=sqrt(eps(max(cur_a,cur_b))) , cur_flag = true) where {T}
  # cur_s = zero(T)
  local cur_s::T
  cur_d, cur_f_d = T(NaN), T(NaN)

  cur_c = middle(cur_a, cur_b)
  cur_f_c = f(cur_c)
  # for iter = 1:100
  for iteration ∈ 1:100
    init_a, init_b = cur_a, cur_b
    init_f_a, init_f_b = cur_f_a, cur_f_b

    if abs(cur_f_a) < abs(cur_f_b)
      cur_a, cur_b = cur_b, cur_a
      cur_f_a, cur_f_b = cur_f_b, cur_f_a
    end

    check_approx(cur_f_a, 0.0, abstol) && return cur_a
    check_approx(cur_f_b, 0.0, abstol) && return cur_b
    check_approx(cur_f_c, 0.0, abstol) && return cur_c
    check_approx(cur_f_d, 0.0, abstol) && return cur_d

    check_approx(cur_a, cur_b, abstol) && return cur_c
    check_approx(cur_f_a, cur_f_b, abstol) && return NaN

    is_good_a = isfinite(cur_f_a)
    is_good_b = isfinite(cur_f_b)
    is_good_c = isfinite(cur_f_c)
    is_good_d = isfinite(cur_f_d)

    ( is_good_a && is_good_b && cur_f_a * cur_f_b > 0 ) && return NaN

    ( is_good_a || is_good_b ) || return NaN

    if !(is_good_a && is_good_b && is_good_c)
      if is_good_d && ( check_approx(cur_c, cur_a, abstol) || check_approx(cur_c, cur_b, abstol) )
        cur_g = cur_d
        cur_f_g = cur_f_d
      else
        cur_g = cur_c
        cur_f_g = cur_f_c
      end

      if is_good_a
        cur_root = custom_bisection(f, cur_a, cur_b, cur_f_a, cur_f_b, abstol, cur_flag)
        isnan(cur_root) &&
          ( cur_root = custom_bisection(f, cur_a, cur_b, cur_f_a, cur_f_b, abstol, cur_flag) )
        return cur_root
      else
        cur_root = custom_bisection(f, cur_a, cur_b, cur_f_a, cur_f_b, abstol, cur_flag)
        isnan(cur_root) &&
          ( cur_root = custom_bisection(f, cur_a, cur_b, cur_f_a, cur_f_b, abstol, cur_flag) )
        return cur_root
      end
    end

    # cur_d = NaN
    # cur_s = zero(T)

    if !check_approx(cur_f_a, cur_f_c, abstol) && !check_approx(cur_f_b, cur_f_c, abstol)
      cur_term_1 = cur_f_b / ( cur_f_a - cur_f_c )
      cur_term_2 = cur_f_c / ( cur_f_a - cur_f_b )
      cur_term_3 = cur_f_a / ( cur_f_b - cur_f_c )

      cur_s = cur_a * cur_term_1 * cur_term_2
      cur_s -= cur_b * cur_term_2 * cur_term_3
      cur_s += cur_c * cur_term_3 * cur_term_1
    else
      cur_s = cur_b - cur_f_b * ( cur_b - cur_a ) / ( cur_f_b - cur_f_a )
    end

    s_bound = ( T(3) * cur_a + cur_b ) * T(0.25)
    cur_bool = ( cur_s < cur_b && cur_s < s_bound  ) || ( cur_s > cur_b && cur_s > s_bound )

    if cur_flag
      cur_diff = abs( cur_b - cur_c )
    else
      cur_diff = abs( cur_d - cur_c )
    end

    cur_bool |= abs( cur_s - cur_b ) >= ( cur_diff * T(0.5) )
    cur_bool |= check_approx(cur_diff, 0.0, abstol)

    cur_flag = cur_bool
    cur_bool && ( cur_s = middle(cur_a, cur_b) )

    cur_d = cur_c
    cur_c = cur_b

    cur_f_d = cur_f_c
    cur_f_c = cur_f_b

    cur_f_s = f(cur_s)

    if cur_f_a * cur_f_s < zero(T)
      cur_b = cur_s
      cur_f_b = cur_f_s
    else
      cur_a = cur_s
      cur_f_a = cur_f_s
    end

    # is_bad_f_a = check_approx(cur_f_a, init_f_a, abstol)
    # is_bad_f_b = check_approx(cur_f_b, init_f_b, abstol)
    #
    # is_bad_f_a |= abs(cur_f_a) > abs(init_f_a)
    # is_bad_f_b |= abs(cur_f_b) > abs(init_f_b)
    #
    # if ( is_bad_f_a && is_bad_f_b )
    #   bad_streak += 1
    # else
    #   bad_streak = 0
    # end
    #
    # ( bad_streak > 4 ) && return NaN

    # is_unchanged = check_approx(cur_a, init_a, atol=2*eps())
    # is_unchanged &= check_approx(cur_b, init_b, atol=2*eps())
    #
    # if is_unchanged
    #   cur_c = middle(cur_a, cur_b)
    #   cur_f_c = f(cur_c)
    #   cur_flag = true
    #   # cur_d, cur_f_d = NaN, NaN
    # end
  end
  cur_s
end#cur_c::Number=NaN, cur_f_c::Number=NaN, cur_d::Number=NaN, cur_f_d::Number=NaN, bad_streak

function custom_bisection_old(f, cur_a::T, cur_b::T, cur_f_a::T, cur_f_b::T, cur_flag::Bool=true, abstol::Number=sqrt(eps(max(cur_a,cur_b))), cur_c::T=NaN, cur_f_c::T=NaN, cur_d::T=NaN, cur_f_d::T=NaN, bad_streak::Integer=0) where {T}
  init_a = cur_a
  init_b = cur_b

  # @show cur_a, cur_b

  init_f_a = cur_f_a
  init_f_b = cur_f_b

  isnan(cur_c) && ( cur_c = middle(cur_a, cur_b) )
  isnan(cur_f_c) && ( cur_f_c = f(cur_c) )

  if abs(cur_f_a) < abs(cur_f_b)
    cur_a, cur_b = cur_b, cur_a
    cur_f_a, cur_f_b = cur_f_b, cur_f_a
  end

  check_approx(cur_f_a, 0.0, abstol) && return cur_a
  check_approx(cur_f_b, 0.0, abstol) && return cur_b
  check_approx(cur_f_c, 0.0, abstol) && return cur_c
  check_approx(cur_f_d, 0.0, abstol) && return cur_d

  check_approx(cur_a, cur_b, abstol) && return cur_c
  check_approx(cur_f_a, cur_f_b, abstol) && return NaN

  is_good_a = isfinite(cur_f_a)
  is_good_b = isfinite(cur_f_b)
  is_good_c = isfinite(cur_f_c)
  is_good_d = isfinite(cur_f_d)

  ( is_good_a && is_good_b && cur_f_a * cur_f_b > 0 ) && return NaN

  ( is_good_a || is_good_b ) || return NaN

  if !(is_good_a && is_good_b && is_good_c)
    if is_good_d && ( check_approx(cur_c, cur_a, abstol) || check_approx(cur_c, cur_b, abstol) )
      cur_g = cur_d
      cur_f_g = cur_f_d
    else
      cur_g = cur_c
      cur_f_g = cur_f_c
    end

    # if is_bad_a
    #   cur_root = custom_bisection_old(f, cur_g, cur_b, cur_f_c, cur_f_b, abstol=abstol)
    #   isnan(cur_root) &&
    #     ( cur_root = custom_bisection_old(f, cur_a, cur_g, cur_f_a, cur_f_c, abstol=abstol) )
    #   return cur_root
    # else
    #   cur_root = custom_bisection_old(f, cur_a, cur_g, cur_f_a, cur_f_c, abstol=abstol)
    #   isnan(cur_root) &&
    #     ( cur_root = custom_bisection_old(f, cur_g, cur_b, cur_f_c, cur_f_b, abstol=abstol) )
    #   return cur_root
    # end
    if is_good_a
      cur_root = custom_bisection_old(f, cur_a, cur_b, cur_f_a, cur_f_b, cur_flag, abstol)
      isnan(cur_root) &&
        ( cur_root = custom_bisection_old(f, cur_a, cur_b, cur_f_a, cur_f_b, cur_flag, abstol) )
      return cur_root
    else
      cur_root = custom_bisection_old(f, cur_a, cur_b, cur_f_a, cur_f_b, cur_flag, abstol)
      isnan(cur_root) &&
        ( cur_root = custom_bisection_old(f, cur_a, cur_b, cur_f_a, cur_f_b, cur_flag, abstol) )
      return cur_root
    end
  end

  # cur_d = NaN
  # cur_s = zero(T)
  local cur_s::T

  if !check_approx(cur_f_a, cur_f_c, abstol) && !check_approx(cur_f_b, cur_f_c, abstol)
    cur_term_1 = cur_f_b / ( cur_f_a - cur_f_c )
    cur_term_2 = cur_f_c / ( cur_f_a - cur_f_b )
    cur_term_3 = cur_f_a / ( cur_f_b - cur_f_c )

    cur_s = cur_a * cur_term_1 * cur_term_2
    cur_s -= cur_b * cur_term_2 * cur_term_3
    cur_s += cur_c * cur_term_3 * cur_term_1
  else
    # cur_term = ( cur_b - cur_a )
    # cur_term /= ( cur_f_b - cur_f_a )
    #
    # cur_s = cur_b
    # cur_s -= cur_f_b * cur_term

    cur_s = cur_b - cur_f_b * ( cur_b - cur_a ) / ( cur_f_b - cur_f_a )
  end

  s_bound = ( T(3) * cur_a + cur_b ) * T(0.25)
  cur_bool = ( cur_s < cur_b && cur_s < s_bound  ) || ( cur_s > cur_b && cur_s > s_bound )

  if cur_flag
    cur_diff = abs( cur_b - cur_c )
  else
    cur_diff = abs( cur_d - cur_c )
  end

  cur_bool |= abs( cur_s - cur_b ) >= ( cur_diff * T(0.5) )
  cur_bool |= check_approx(cur_diff, 0.0, abstol)

  cur_flag = cur_bool
  cur_bool && ( cur_s = middle(cur_a, cur_b) )

  cur_d = cur_c
  cur_c = cur_b

  cur_f_d = cur_f_c
  cur_f_c = cur_f_b

  cur_f_s = f(cur_s)

  if cur_f_a * cur_f_s < zero(T)
    cur_b = cur_s
    cur_f_b = cur_f_s
  else
    cur_a = cur_s
    cur_f_a = cur_f_s
  end

  is_bad_f_a = check_approx(cur_f_a, init_f_a, abstol)
  is_bad_f_b = check_approx(cur_f_b, init_f_b, abstol)

  is_bad_f_a |= abs(cur_f_a) > abs(init_f_a)
  is_bad_f_b |= abs(cur_f_b) > abs(init_f_b)

  if ( is_bad_f_a && is_bad_f_b )
    bad_streak += 1
  else
    bad_streak = 0
  end

  ( bad_streak > 4 ) && return NaN

  is_unchanged = check_approx(cur_a, init_a, 2*eps())
  is_unchanged &= check_approx(cur_b, init_b, 2*eps())

  if is_unchanged
    cur_c = middle(cur_a, cur_b)
    cur_f_c = f(cur_c)

    # cur_root = custom_bisection_old(f, cur_a, cur_c, cur_f_a, cur_f_c, abstol=abstol, bad_streak=bad_streak)
    cur_root = custom_bisection_old(f, cur_a, cur_b, cur_f_a, cur_f_b, cur_flag, abstol, NaN, NaN, NaN, NaN, bad_streak)
    isnan(cur_root) || return cur_root

    # cur_root = custom_bisection_old(f, cur_c, cur_b, cur_f_c, cur_f_b, abstol=abstol, bad_streak=bad_streak)
    # return cur_root
    return custom_bisection_old(f, cur_a, cur_b, cur_f_a, cur_f_b, cur_flag, abstol, NaN, NaN, NaN, NaN, bad_streak)
  end

  # custom_bisection_old(
  #   f, cur_a, cur_b, cur_f_a, cur_f_b,
  #   cur_flag, abstol=abstol,
  #   cur_c=cur_c, cur_f_c=cur_f_c,
  #   cur_d=cur_d, cur_f_d=cur_f_d,
  #   bad_streak=bad_streak
  # )
  custom_bisection_old(
    f, cur_a, cur_b, cur_f_a, cur_f_b,
    cur_flag, abstol,
    cur_c, cur_f_c,
    cur_d, cur_f_d,
    bad_streak
  )
end#cur_c::Number=NaN, cur_f_c::Number=NaN, cur_d::Number=NaN, cur_f_d::Number=NaN, bad_streak


text:
module CartesianMeshes

using StaticArrays

export AbstractCartesianMesh, CartesianMesh, PeriodicCartesianMesh
export nb_dims, nb_cells, all_cells
export nb_inner_faces, inner_faces
export nb_boundary_faces, boundary_faces
export cell_centers, cells_centers, cell_center, face_center
export cells_next_to_inner_face, cell_next_to_boundary_face
export cell_volume, face_area, normal_vector

# Tools
struct Half{N} <: Real
    n::N
end
import Base.Int, Base.convert, Base.+, Base.-, Base.<, Base.==, Base.print
Base.Int(h::Half{Int}) = Base.div(h.n, 2)
Base.convert(::Type{Int}, h::Half{Int}) = Int(h)
Base.convert(::Type{<:Number}, h::Half) = h.n/2
Base.:+(a::Half, b::Half) = Half(a.n+b.n)
Base.:-(a::Half) = Half(-a.n)
Base.:-(a::Half, b::Half) = Half(a.n-b.n)
Base.:(==)(a::Half, b::Half) = a.n == b.n
Base.:<(a::Half, b::Half) = a.n < b.n
is_int(h::Half{Int}) = mod(h.n, 2) == 0
Base.print(io::IO, h::Half) = print(io, "$(h.n)/2")

cell_centers(mesh) = map(i -> cell_center(mesh, i), all_cells(mesh))
cells_centers(mesh) = cell_centers(mesh)

const FaceIndex{N} = NTuple{N, Half{Int}}
const CellIndex{N} = CartesianIndex{N}

#################################################################
abstract type AbstractCartesianMesh{D, L} end
#################################################################

#########################################################
struct CartesianMesh{D, L} <: AbstractCartesianMesh{D, L}
    x_min::SVector{D, L}
    x_max::SVector{D, L}
    nb_cells::SVector{D, Int64}
end
CartesianMesh(x_min::Number, x_max::Number, nb_cells::Int) = CartesianMesh{1, typeof(x_min)}(SVector(x_min), SVector(x_max), (nb_cells,))
CartesianMesh(nb_cells::Int) = CartesianMesh{1, Float64}(SVector(0.0), SVector(1.0), (nb_cells,))
CartesianMesh(nx::Int, ny::Int) = CartesianMesh{2, Float64}(SVector(0.0, 0.0), SVector(1.0, 1.0), (nx, ny))
CartesianMesh(nx::Int, ny::Int, nz::Int) = CartesianMesh{3, Float64}(SVector(0.0, 0.0, 0.0), SVector(1.0, 1.0, 1.0), (nx, ny, nz))

CartesianMesh(x_min::NTuple{N}, x_max::NTuple{N}, nb_cells::NTuple{N, Int64}) where N = CartesianMesh{N, eltype(x_min)}(SVector(x_min...), SVector(x_max...), SVector(nb_cells...))

#################################################################
struct PeriodicCartesianMesh{D, L} <: AbstractCartesianMesh{D, L}
    x_min::SVector{D, L}
    x_max::SVector{D, L}
    nb_cells::SVector{D, Int64}
end
PeriodicCartesianMesh(x_min::Number, x_max::Number, nb_cells::Int) = PeriodicCartesianMesh{1, typeof(x_min)}(SVector(x_min), SVector(x_max), (nb_cells,))
PeriodicCartesianMesh(nb_cells::Int) = PeriodicCartesianMesh{1, Float64}(SVector(0.0), SVector(1.0), (nb_cells,))
PeriodicCartesianMesh(nx::Int, ny::Int) = PeriodicCartesianMesh{2, Float64}(SVector(0.0, 0.0), SVector(1.0, 1.0), (nx, ny))
PeriodicCartesianMesh(nx::Int, ny::Int, nz::Int) = PeriodicCartesianMesh{3, Float64}(SVector(0.0, 0.0, 0.0), SVector(1.0, 1.0, 1.0), (nx, ny, nz))

PeriodicCartesianMesh(x_min::NTuple{N}, x_max::NTuple{N}, nb_cells::NTuple{N, Int64}) where N = PeriodicCartesianMesh{N, eltype(x_min)}(SVector(x_min...), SVector(x_max...), SVector(nb_cells...))

nb_dims(mesh::AbstractCartesianMesh{D}) where D = D

dx(mesh::AbstractCartesianMesh) = @. (mesh.x_max - mesh.x_min)/mesh.nb_cells

nb_cells(mesh::AbstractCartesianMesh) = prod(mesh.nb_cells)
all_cells(mesh::AbstractCartesianMesh) = CartesianIndices(Tuple(mesh.nb_cells))

# INNER FACES
nb_inner_faces(mesh::CartesianMesh{1}) = mesh.nb_cells[1] - 1
inner_faces(mesh::CartesianMesh{1}) = ((Half(n),) for n in 3:2:2*mesh.nb_cells[1])
nb_boundary_faces(mesh::CartesianMesh{1}) = 2
boundary_faces(mesh::CartesianMesh{1}) = ((Half(1),), (Half(2*mesh.nb_cells[1]+1),))

nb_inner_faces(mesh::PeriodicCartesianMesh{1}) = mesh.nb_cells[1]
inner_faces(mesh::PeriodicCartesianMesh{1}) = ((Half(n),) for n in 1:2:2*mesh.nb_cells[1])
nb_boundary_faces(mesh::PeriodicCartesianMesh{1}) = 0
boundary_faces(mesh::PeriodicCartesianMesh{1}) = Tuple([])

const _left = (-Half(1), Half(0), Half(0))
const _right = (Half(1), Half(0), Half(0))
const _down = (Half(0), -Half(1), Half(0))
const _up = (Half(0), Half(1), Half(0))
const _back = (Half(0), Half(0), -Half(1))
const _front = (Half(0), Half(0), Half(1))
_bottom_cells(mesh::CartesianMesh{2}) = ((n, m) = mesh.nb_cells; CartesianIndices((1:n, 1:1)))
_top_cells(mesh::CartesianMesh{2}) = ((n, m) = mesh.nb_cells; CartesianIndices((1:n, m:m)))
_left_cells(mesh::CartesianMesh{2}) = ((n, m) = mesh.nb_cells; CartesianIndices((1:1, 1:m)))
_right_cells(mesh::CartesianMesh{2}) = ((n, m) = mesh.nb_cells; CartesianIndices((n:n, 1:m)))

_face(dir, cell::CellIndex{N}) where N = Half.(2 .* Tuple(cell)) .+ dir[1:N]

nb_inner_faces(mesh::CartesianMesh{2}) = 2 * mesh.nb_cells[1] * mesh.nb_cells[2] - mesh.nb_cells[1] - mesh.nb_cells[2]
inner_faces(mesh::CartesianMesh{2}) = ((n, m) = mesh.nb_cells; (_face(dir, cell) for (dir, cells) in ((_right, CartesianIndices((1:(n-1), 1:m))), (_up, CartesianIndices((1:n, 1:(m-1))))) for cell in cells))
nb_boundary_faces(mesh::CartesianMesh{2}) = 2 * (mesh.nb_cells[1] + mesh.nb_cells[2])
boundary_faces(mesh::CartesianMesh{2}) = (_face(dir, cell) for (dir, cells) in ((_down, _bottom_cells(mesh)), (_right, _right_cells(mesh)), (_up, _top_cells(mesh)), (_left, _left_cells(mesh))) for cell in cells)

nb_inner_faces(mesh::PeriodicCartesianMesh) = 2 * prod(mesh.nb_cells)
inner_faces(mesh::PeriodicCartesianMesh{N}) where N = (_face(dir, cell) for cell in CartesianIndices(Tuple(mesh.nb_cells)) for dir in (_left, _down, _back)[1:N])
nb_boundary_faces(mesh::PeriodicCartesianMesh) = 0
boundary_faces(mesh::PeriodicCartesianMesh) = Tuple([])

# CELLS NEXT TO FACE
_direction(i_face::FaceIndex{1}) = 1
_direction(i_face::FaceIndex{2}) = !is_int(i_face[1]) ? 1 : 2  # 1 = normale horizontale, 2 = normale verticale
_direction(i_face::FaceIndex{3}) = !is_int(i_face[1]) ? 1 : !is_int(i_face[2]) ? 2 : 3

_dir_step(i_face::FaceIndex{1}) = (1,)
_dir_step(i_face::FaceIndex{2}) = _direction(i_face) == 1 ? (1, 0) : (0, 1)
_dir_step(i_face::FaceIndex{3}) = _direction(i_face) == 1 ? (1, 0, 0) : _direction(i_face) == 2 ? (0, 1, 0) : (0, 0, 1)

_cells_next(i_face::FaceIndex{N}) where N = (Int.(i_face), Int.(i_face) .+ _dir_step(i_face))
_cells_next(i_face::Half{Int}) = _cells_next((i_face,))

cells_next_to_inner_face(::CartesianMesh, i_face) = CellIndex.(_cells_next(i_face))
function cells_next_to_inner_face(mesh::PeriodicCartesianMesh{N}, i_face::FaceIndex{N}) where N
    if Half(1) in i_face
        i_cell_1, i_cell_2 = _cells_next(i_face)
        return (CellIndex(mod1.(i_cell_1, Tuple(mesh.nb_cells))),
                CellIndex(mod1.(i_cell_2, Tuple(mesh.nb_cells))))
    else
        return CellIndex.(_cells_next(i_face))
    end
end

cell_next_to_boundary_face(mesh::CartesianMesh{1}, i_face) = i_face == (Half(1),) ? CellIndex(1) : CellIndex(mesh.nb_cells[1])

function cell_next_to_boundary_face(mesh::CartesianMesh{2}, i_face)
    if i_face[1] == Half(1)  # left boundary
        return CellIndex(Int(i_face[1]) + 1, Int(i_face[2]))
    elseif i_face[2] == Half(1)  # bottom boundary
        return CellIndex(Int(i_face[1]), Int(i_face[2]) + 1)
    elseif _direction(i_face) == 1  # right boundary
        return CellIndex(Int(i_face[1]), Int(i_face[2]))
    elseif _direction(i_face) == 2  # top boundary
        return CellIndex(Int(i_face[1]), Int(i_face[2]))
    else
        error()
    end
end

# GEOMETRY
cell_center(mesh::AbstractCartesianMesh{1}, i_cell) = (i_cell[1] - 0.5) * dx(mesh)[1]
cell_center(mesh::AbstractCartesianMesh{N}, i_cell) where N = (SVector{N}(Tuple(i_cell)...) .- 0.5) .* dx(mesh)
face_center(mesh::AbstractCartesianMesh{1}, i_face) = (i_face[1] - 0.5) * dx(mesh)[1]
face_center(mesh::AbstractCartesianMesh{N}, i_face) where N = (SVector{N, Float64}(Tuple(i_face)...) .- 0.5) .* dx(mesh)

cell_volume(mesh::AbstractCartesianMesh, i_cell) = prod(dx(mesh))

face_area(mesh::AbstractCartesianMesh{1, T}, i_face) where T = one(T)
face_area(mesh::AbstractCartesianMesh{2}, i_face) = _direction(i_face) == 1 ? dx(mesh)[2] : dx(mesh)[1]
function face_area(mesh::AbstractCartesianMesh{3}, i_face)
    Δx, Δy, Δz = dx(mesh)
    if _direction(i_face) == 1; return Δy*Δz
    elseif _direction(i_face) == 2; return Δx*Δz
    elseif _direction(i_face) == 3; return Δx*Δy
    end
end

normal_vector(mesh::CartesianMesh{1, T}, i_face) where T = i_face == (Half(1),) ? -one(T) : one(T)
normal_vector(mesh::PeriodicCartesianMesh{1, T}, i_face) where T = one(T)

_normal(::Type{T}, i_face) where T = map(x -> x == 1 ? one(T) : zero(T)/oneunit(T), _dir_step(i_face))
# Division by oneunit is required for compatibility with Unitful
function normal_vector(mesh::PeriodicCartesianMesh{N, T}, i_face) where {N, T}
    SVector(_normal(T, i_face)...)
end
function normal_vector(mesh::AbstractCartesianMesh{N, T}, i_face) where {N, T}
    n = SVector(_normal(T, i_face)...)
    if Half(1) in i_face
        return -n
    else
        return n
    end
end

include("./plot_recipes.jl")

end # module