text:
abstract type BoundaryComponent{N} end

"""
    get_positions(boundary_component::BoundaryComponent{N})
        -> Vector{SVectorF{N}}

Return the positions of the boundary particles.
"""
function get_positions end

"""
    get_velocities(boundary_component::BoundaryComponent{N})
        -> Vector{SVectorF{N}}

Return the velocities of the boundary particles.
"""
function get_velocities end

"""
    get_velocity(boundary_component::BoundaryComponent{N}, i::Unsigned)
        -> SVectorF{N}

Return the velocity of the boundary particle with the given index.
"""
function get_velocity(boundary_component::BoundaryComponent, i::Unsigned)
    velocities = get_boundary_velocities(boundary_component)
    velocities[i]
end

"""
    get_pressures(boundary_component::BoundaryComponent)
        -> Vector{Float}

Return the virtual pressures of the boundary particles.
"""
function get_pressures end

"""
    get_pressure(boundary_component::BoundaryComponent, i::Unsigned)
        -> Float

Return the virtual pressure of the boundary particle with the given index.
"""
function get_pressure(boundary_component::BoundaryComponent, i::Unsigned)
    pressures = get_boundary_pressures(boundary_component)
    pressures[i]
end

"""
    compute_boundary_mass_density(
        boundary_component::BoundaryComponent,
        energy_component::EnergyComponent,
        eos::EquationOfState,
        idx_fluid::Unsigned,
        idx_boundary::Unsigned,
    )

Compute the virtual mass density of the boundary particle with index `idx_boundary`
due to its interaction with the fluid particle with index `idx_fluid`.
"""
function compute_mass_density end

"""
    updateboundaries!(
        boundary_component::BoundaryComponent{N},
        positions::AbstractVector{SVectorF{N}},
        velocities::Velocities{N},
        mass_component::MassComponent{MD,<:Kernel{N}},
        pressures::Pressures,
    )

Update the state of the boundary particles according to the given fluid particles.
"""
function updateboundaries! end

"""
    evolveboundaries!(boundary_component::BoundaryComponent, Δt::Number)

Evolves the boundary particles with the given time step.
"""
function evolveboundaries! end

include("boundary/none.jl")
include("boundary/wall.jl")

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
@testset "AnalyticVI" begin
    seed!(42)
    L = 3
    D = 10
    N = 20
    b = 5
    i = AnalyticVI()
    @test i isa AnalyticVI{Float64}
    @test AnalyticVI(; ϵ=0.0001f0) isa AnalyticVI{Float32}
    @test AnalyticSVI(10; ϵ=0.0001f0) isa AnalyticVI{Float32}
    @test repr(i) == "Analytic Variational Inference"
    @test AGP.ρ(i) == 1.0
    @test AGP.is_stochastic(i) == false

    i = AnalyticSVI(b)
    @test i isa AnalyticVI{Float64}
    AGP.set_ρ!(i, N / b)
    @test AGP.ρ(i) == N / b
    @test AGP.is_stochastic(i) == true
end

ext: jl


text:
"""
Main module for `GigaSOM.jl` - Huge-scale, high-performance flow cytometry clustering

The documentation is here: http://LCSB-BioCore.github.io/GigaSOM.jl
"""

module GigaSOM

using CSV
using DataFrames
using Distances
using Distributed
using DistributedData
using Distributions
using FCSFiles
using FileIO
using DistributedArrays
using NearestNeighbors
using Serialization
using StableRNGs

include("base/structs.jl")

include("base/dataops.jl")
include("base/trainutils.jl")

include("analysis/core.jl")
include("analysis/embedding.jl")

include("io/input.jl")
include("io/process.jl")
include("io/splitting.jl")

#core
export initGigaSOM, trainGigaSOM, mapToGigaSOM

#trainutils
export linearRadius, expRadius, gaussianKernel, bubbleKernel, thresholdKernel, distMatrix

#embedding
export embedGigaSOM

# structs
export Som

#io/input
export readFlowset,
    readFlowFrame,
    loadFCS,
    loadFCSHeader,
    getFCSSize,
    loadFCSSizes,
    loadFCSSet,
    selectFCSColumns,
    distributeFCSFileVector,
    distributeFileVector,
    getCSVSize,
    loadCSV,
    loadCSVSizes,
    loadCSVSet

#io/splitting
export slicesof, vcollectSlice, collectSlice

#io/process
export cleanNames!, getMetaData, getMarkerNames

#dataops (higher-level operations on data)
export dtransform_asinh

end # module


text:
const _parker_params = (
    A = [0.809, 0.330],
    T = [0.17046, 0.365],
    σ = [0.0563, 0.132],
    α = 1.050,
    β = 0.1685,
    s = 38.078,
    τ = 0.483,
)

"""
    aif_parker(t; params, hct = 0.42)

Returns the arterial input function defined in Parker et al. [@Parker2006].
The timepoints `t` are in minutes with the bolus arrival time defined as `t = 0`.
The model parameters from Ref [@Parker2006] will be used by default unless an optional `params` input is provided.
The hematocrit `hct` is used to convert the concentration in blood to plasma and has default value of 0.42.
"""
function aif_parker(t; params = _parker_params, hct = 0.42)
    @extract (A, T, σ, α, β, s, τ) params
    cb = [
        (A[1] / (σ[1] * sqrt(2 * pi))) * exp(-(t - T[1])^2 / (2 * (σ[1])^2)) +
        (A[2] / (σ[2] * sqrt(2 * pi))) * exp(-(t - T[2])^2 / (2 * (σ[2])^2)) +
        α * exp(-β * t) / (1.0 + exp(-s * (t - τ))) for t in t
    ]
    cb[t.<=0] .= 0
    ca = cb ./ (1 - hct)
    return ca
end

const _georgiou_params =
    (A = [0.37, 0.33, 10.06], m = [0.11, 1.17, 16.02], α = 5.26, β = 0.032, τ = 0.129)

"""
    aif_georgiou(t; params, hct = 0.35)

Returns the concentration in blood plasma using the model defined in Georgiou et al. [@Georgiou2018].
The timepoints `t` are in minutes with the bolus arrival defined as `t = 0`.
The model parameters from the paper will be used by default unless an optional `params` input is provided.
The hematocrit `hct` is used to convert the concentration in blood to plasma and has default value of 0.35.
"""
function aif_georgiou(t; params = _georgiou_params, hct = 0.35)
    @extract (A, m, α, β, τ) params
    # The AIF can be split into two parts
    exponential_part = zeros(size(t))
    gammavariate_part = zeros(size(t))

    for (i, t) in enumerate(t)
        exponential_part[i] = sum(A .* exp.(-m * t))
        if t > 7.5
            # Gamma variate part undefined for high N (i.e. t > 7.5 min), so fix its value
            gammavariate_part[i] = 3.035
            continue
        end
        N = fld(t, τ) # `fld` might be more efficient than `floor`
        for j = 0:N
            gammavariate_part[i] += gammavariate((j + 1) * α + j, β, t - j * τ)
        end
    end
    cb = exponential_part .* gammavariate_part
    cb[t.<=0] .= 0
    ca = cb ./ (1 - hct)
    return ca
end

function gammavariate(α, β, t)::Float64
    if t > 0
        gammavalue = (t^α * exp(-t / β)) / (β^(α + 1) * gamma(α + 1))
        if !isnan(gammavalue)
            return gammavalue
        end
    end
    return 0
end

"""
    aif_biexponential(t; params, hct = 0)

Returns an AIF defined by the biexponential model with the form ``D * (a_1 exp(-m_1 t) + a_2 exp(-m_2 t))``.
The `t` input is in minutes with the bolus arriving at 0.
The `params` can be a named tuple, e.g. `params = (D=0.1, a=[4.0, 4.78], m=[0.144, 0.011])`.
The `hct` is used to convert the concentration in blood to plasma.
The default value is set to to zero, i.e. no conversion takes place.
"""
function aif_biexponential(t; params, hct = 0)
    @extract (D, a, m) params
    cb = [D * (a[1] * exp(-m[1] * t) + a[2] * exp(-m[2] * t)) for t in t]
    cb[t.<=0] .= 0
    ca = cb ./ (1 - hct)
    return ca
end

"""
    aif_weinmann(t)

Returns an AIF defined by the biexponential model with parameters from
Tofts & Kermode [@Tofts1991] and Weinmann et al. [@Weinmann1984].
The `t` input is in minutes with the bolus arriving at 0.
The model describes concentration in plasma, therefore a hematocrit is not necessary.
"""
function aif_weinmann(t)
    params = (; D = 0.1, a = [3.99, 4.78], m = [0.144, 0.0111])
    return aif_biexponential(t; params, hct = 0)
end

"""
    aif_fritzhansen(t)

Returns an AIF defined by the biexponential model with parameters based on data from
Fritz-Hansen et al. [@Fritz-Hansen1996].
The model parameters were adopted from Whitcher & Schmid [@Whitcher2011]
The `t` input is in minutes with the bolus arriving at 0.
The model describes concentration in plasma, therefore a hematocrit is not necessary.
"""
function aif_fritzhansen(t; hct = 0)
    params = (; D = 1.0, a = [2.4, 0.62], m = [3.0, 0.016])
    return aif_biexponential(t; params, hct)
end

"""
    aif_orton1(t; params = (aB = 10.4, μB = 12.0, aG = 1.24, μG = 0.169), hct = 0.42)

Returns the concentration in blood plasma using the bi-exponential model defined in
Orton et al. [@Orton2008].
The timepoints `t` are in minutes with the bolus arriving at 0.
The model parameters from the paper will be used by default unless an optional `params` input is provided.
The `hct` is used to convert the concentration in blood to plasma and has default value of 0.42.
"""
function aif_orton1(t; params = (aB = 10.4, μB = 12.0, aG = 1.24, μG = 0.169), hct = 0.42)
    @extract (aB, μB, aG, μG) params
    AB = aB - aB * aG / (μB - μG)
    AG = aB * aG / (μB - μG)
    cb = [AB * exp(-μB * t) + AG * exp(-μG * t) for t in t]
    cb[t.<0] .= 0
    ca = cb ./ (1 - hct)
    return ca
end

"""
    aif_orton2(t; params = (aB = 344, μB = 20.2, aG = 1.24, μG = 0.172), hct = 0.42)

Returns the concentration in blood plasma using the 2nd model defined in
Orton et al. [@Orton2008].
The timepoints `t` are in minutes with the bolus arriving at 0.
The model parameters from the paper will be used by default unless an optional `params` input is provided.
The `hematocrit` is used to convert the concentration in blood to plasma and has default value of 0.42.
"""
function aif_orton2(t; params = (aB = 344, μB = 20.2, aG = 1.24, μG = 0.172), hct = 0.42)
    @extract (aB, μB, aG, μG) params
    AB = aB - aB * aG / (μB - μG)
    AG = aB * aG / (μB - μG)^2
    cb = [AB * t * exp(-μB * t) + AG * (exp(-μG * t) - exp(-μB * t)) for t in t]
    cb[t.<=0] .= 0
    ca = cb ./ (1 - hct)
    return ca
end

"""
    aif_orton3(t; params = (aB = 2.84, μB = 22.8, aG = 1.36, μG = 0.171), hct = 0.42)

Returns the concentration in blood plasma using the 3rd model defined in
Orton et al. [@Orton2008].
The timepoints `t` are in minutes with the bolus arriving at `t = 0`.
The model parameters from the paper will be used by default unless an optional `params` input is provided.
The `hct` is used to convert the concentration in blood to plasma and has default value of 0.42.
"""
function aif_orton3(t; params = (aB = 2.84, μB = 22.8, aG = 1.36, μG = 0.171), hct = 0.42)
    @extract (aB, μB, aG, μG) params
    tb = 2 * pi / μB
    cb = zeros(length(t))
    for (i, t) in enumerate(t)
        if t <= tb
            cb[i] = aB * (1 - cos(μB * t)) + aB * aG * _orton3_f(t, μG, μB)
        else
            cb[i] = aB * aG * _orton3_f(tb, μG, μB) * exp(-μG * (t - tb))
        end
    end
    cb[t.<=0] .= 0
    ca = cb ./ (1 - hct)
    return ca
end

function _orton3_f(t, α, μB)
    (1 / α) * (1 - exp(-α * t)) -
    (1 / (α^2 + μB^2)) * (α * cos(μB * t) + μB * sin(μB * t) - α * exp(-α * t))
end


text:
@enum dq_ref begin
    q = 1
    d = 2
end
@enum RI_ref begin
    R = 1
    I = 2
end
function dq_ri(δ::Float64)
    ## Uses the referenceframe of the Kundur page 852 of dq to RI
    dq_ri =  [sin(δ)  cos(δ);
             -cos(δ)  sin(δ)]
end

function ri_dq(δ::Float64)
    #Uses the reference frame of the Kundur page 852 of RI to dq
    ri_dq = [sin(δ) -cos(δ);
            cos(δ)  sin(δ)]
end


text:
using JSON

@testset "test matpower parser" begin
    @testset "30-bus case file" begin
        result = run_opf("../test/data/case30.m", ACPPowerModel, ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 204.96; atol = 1e-1)
    end

    @testset "30-bus case matpower data" begin
        data = PowerModels.parse_file("../test/data/case30.m")
        @test isa(JSON.json(data), String)

        result = run_opf(data, ACPPowerModel, ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 204.96; atol = 1e-1)
    end

    @testset "14-bus case file with bus names" begin
        data = PowerModels.parse_file("../test/data/case14.m")
        @test data["bus"]["1"]["bus_name"] == "Bus 1     HV"
        @test isa(JSON.json(data), String)
    end

    @testset "5-bus case file with pwl cost functions" begin
        data = PowerModels.parse_file("../test/data/case5_pwlc.m")
        @test data["gen"]["1"]["model"] == 1
        @test isa(JSON.json(data), String)
    end

    @testset "3-bus case file with hvdc lines" begin
        data = PowerModels.parse_file("../test/data/case3.m")
        @test length(data["dcline"]) > 0
        @test isa(JSON.json(data), String)
    end

    @testset "2-bus case file with spaces" begin
        result = run_pf("../test/data/case2.m", ACPPowerModel, ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 0.0; atol = 1e-1)
    end
end

@testset "test matpower data coercion" begin
    @testset "ACP Model" begin
        result = run_opf("../test/data/case14.m", ACPPowerModel, ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 8081.5; atol = 1e0)
        #@test result["status"] = bus_name
    end
    @testset "DC Model" begin
        result = run_opf("../test/data/case14.m", DCPPowerModel, ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 7642.6; atol = 1e0)
    end
    @testset "QC Model" begin
        result = run_opf("../test/data/case14.m", QCWRPowerModel, ipopt_solver)

        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 8075.1; atol = 1e0)
    end
end

@testset "test matpower extentions parser" begin
    @testset "3-bus extended constants" begin
        data = PowerModels.parse_file("../test/data/case3.m")

        @test data["const_int"] == 123
        @test data["const_float"] == 4.56
        @test data["const_str"] == "a string"
        @test isa(JSON.json(data), String)
    end

    @testset "3-bus extended matrix" begin
        data = PowerModels.parse_file("../test/data/case3.m")

        @test haskey(data, "areas")
        @test data["areas"]["1"]["col_1"] == 1
        @test data["areas"]["1"]["col_2"] == 1
        @test data["areas"]["2"]["col_1"] == 2
        @test data["areas"]["2"]["col_2"] == 3
        @test isa(JSON.json(data), String)
    end

    @testset "3-bus extended named matrix" begin
        data = PowerModels.parse_file("../test/data/case3.m")

        @test haskey(data, "areas_named")
        @test data["areas_named"]["1"]["area"] == 4
        @test data["areas_named"]["1"]["refbus"] == 5
        @test data["areas_named"]["2"]["area"] == 5
        @test data["areas_named"]["2"]["refbus"] == 6
        @test isa(JSON.json(data), String)
    end

    @testset "3-bus extended predefined matrix" begin
        data = PowerModels.parse_file("../test/data/case3.m")

        @test haskey(data, "areas_named")
        @test data["branch"]["1"]["rate_i"] == 50.2
        @test data["branch"]["1"]["rate_p"] == 45
        @test data["branch"]["2"]["rate_i"] == 36
        @test data["branch"]["2"]["rate_p"] == 60.1
        @test data["branch"]["3"]["rate_i"] == 12
        @test data["branch"]["3"]["rate_p"] == 30
        @test isa(JSON.json(data), String)
    end

    @testset "3-bus extended matrix from cell" begin
        data = PowerModels.parse_file("../test/data/case3.m")

        @test haskey(data, "areas_cells")
        @test data["areas_cells"]["1"]["col_1"] == "Area 1"
        @test data["areas_cells"]["1"]["col_2"] == 123
        @test data["areas_cells"]["1"]["col_4"] == "Slack \\\"Bus\\\" 1"
        @test data["areas_cells"]["1"]["col_5"] == 1.23
        @test data["areas_cells"]["2"]["col_1"] == "Area 2"
        @test data["areas_cells"]["2"]["col_2"] == 456
        @test data["areas_cells"]["2"]["col_4"] == "Slack Bus 3"
        @test data["areas_cells"]["2"]["col_5"] == 4.56
        @test isa(JSON.json(data), String)
    end

    @testset "3-bus extended named matrix from cell" begin
        data = PowerModels.parse_file("../test/data/case3.m")

        @test haskey(data, "areas_named_cells")
        @test data["areas_named_cells"]["1"]["area_name"] == "Area 1"
        @test data["areas_named_cells"]["1"]["area"] == 123
        @test data["areas_named_cells"]["1"]["area2"] == 987
        @test data["areas_named_cells"]["1"]["refbus_name"] == "Slack Bus 1"
        @test data["areas_named_cells"]["1"]["refbus"] == 1.23
        @test data["areas_named_cells"]["2"]["area_name"] == "Area 2"
        @test data["areas_named_cells"]["2"]["area"] == 456
        @test data["areas_named_cells"]["2"]["area2"] == 987
        @test data["areas_named_cells"]["2"]["refbus_name"] == "Slack Bus 3"
        @test data["areas_named_cells"]["2"]["refbus"] == 4.56
        @test isa(JSON.json(data), String)
    end

    @testset "3-bus extended predefined matrix from cell" begin
        data = PowerModels.parse_file("../test/data/case3.m")

        @test haskey(data, "areas_named")
        @test data["branch"]["1"]["name"] == "Branch 1"
        @test data["branch"]["1"]["number_id"] == 123
        @test data["branch"]["2"]["name"] == "Branch 2"
        @test data["branch"]["2"]["number_id"] == 456
        @test data["branch"]["3"]["name"] == "Branch 3"
        @test data["branch"]["3"]["number_id"] == 789
        @test isa(JSON.json(data), String)
    end

    @testset "3-bus tnep case" begin
        data = PowerModels.parse_file("../test/data/case3_tnep.m")

        @test haskey(data, "ne_branch")
        @test data["ne_branch"]["1"]["f_bus"] == 1
        @test data["ne_branch"]["1"]["construction_cost"] == 1
        @test isa(JSON.json(data), String)
    end

    @testset "`build_ref` for 3-bus tnep case" begin
        data = PowerModels.parse_file("../test/data/case3_tnep.m")
        ref = PowerModels.build_ref(data)

        @assert !(data["multinetwork"])
        ref = ref[:nw][0]

        @test haskey(data, "name")
        @test haskey(ref, :name)
        @test data["name"] == ref[:name]
    end
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

fst = WFST(["a","b","c"],["p","q","r"])
add_arc!(fst,1,1,"a","p",2)
add_arc!(fst,1,2,"b","q",3)
add_arc!(fst,2,2,"c","r",4)
initial!(fst,1)
final!(fst,2,5)

fst_plus = closure(fst; star=false)
println(fst_plus)

fst_star = closure(fst; star=true)
println(fst_star)
@test isinitial(fst_star,3)
@test isfinal(fst_star,3)
@test isfinal(fst_star,2)

W = typeofweight(fst)
ilabels=["a","b","c"]
o1,w1 = fst(ilabels)
@test o1 == ["p","q","r"]
@test w1 == W(2)*W(3)*W(4)*W(5)
o2,w2 = fst_star([ilabels;ilabels])
@test o2 == [o1;o1]
@test w2 == w1*w1
o3,w3 = fst_star([ilabels;ilabels])
@test o3 == o3
@test w2 == w3
o,w = fst([ilabels;ilabels])
@test w == zero(W)

# fst star accepts empty string and returns the same with weight one
o,w = fst_star(String[])
@test o == String[]
@test w == one(W)

# fst plus does not accept empty string
o,w = fst_plus(String[])
@test w == zero(W)


text:
function _maximise_nonlinear_through_budgeted(instance::CombinatorialInstance{T}, rewards::Dict{T, Float64},
                                              weights::Dict{T, Int}, max_weight::Int, nl_func::Function,
                                              solve_all_budgets_at_once::Bool) where T
  solutions = _maximise_nonlinear_through_budgeted_sub(instance, rewards, weights, max_weight, Val(solve_all_budgets_at_once))
  # TODO: Replace the mandatory parameter solve_all_budgets_at_once by a function defined on the solvers, so that each of them can indicate what they support? Then, if this parameter is not set, the most efficient implementation is called if available; otherwise, follow this parameter (and still warn if the required implementation is not available).
  # TODO: make it a trait of the underlying algorithm to get a default value. Do the same for ESCB2.

  best_solution = Int[]
  best_objective = -Inf
  for (budget, sol) in solutions
    # Ignore infeasible cases.
    if length(sol) == 0 || sol == [-1]
      continue
    end

    # Compute the maximum.
    f_x = nl_func(sol, budget)
    if f_x > best_objective
      best_solution = sol
      best_objective = f_x
    end
  end

  return best_solution, best_objective
end

function _maximise_nonlinear_through_budgeted_sub(instance::CombinatorialInstance{T}, rewards::Dict{T, Float64},
                                                  weights::Dict{T, Int}, max_weight::Int,
                                                  ::Val{true})::Dict{Int, Vector{T}} where T
  if ! applicable(solve_all_budgeted_linear, instance.solver, rewards, weights, max_weight)
    if applicable(solve_budgeted_linear, instance.solver, rewards, weights, max_weight)
      # @warn("The function solve_all_budgeted_linear is not defined for the solver $(typeof(instance.solver)).")
      return _maximise_nonlinear_through_budgeted_sub(instance, rewards, weights, max_weight, Val(false))
    else
      error("Neither solve_budgeted_linear nor solve_all_budgeted_linear are not defined for the solver $(typeof(instance.solver)).")
    end
  end

  return solve_all_budgeted_linear(instance.solver, rewards, weights, max_weight)
end

function _maximise_nonlinear_through_budgeted_sub(instance::CombinatorialInstance{T}, rewards::Dict{T, Float64},
                                                  weights::Dict{T, Int}, max_weight::Int,
                                                  ::Val{false})::Dict{Int, Vector{T}} where T
  if ! applicable(solve_budgeted_linear, instance.solver, rewards, weights, max_weight)
    if applicable(solve_all_budgeted_linear, instance.solver, rewards, weights, max_weight)
      # @warn("The function solve_budgeted_linear is not defined for the solver $(typeof(instance.solver)).")
      return _maximise_nonlinear_through_budgeted_sub(instance, rewards, weights, max_weight, Val(true))
    else
      error("Neither solve_budgeted_linear nor solve_all_budgeted_linear are not defined for the solver $(typeof(instance.solver)).")
    end
  end

  # Assumption: rewards is more complete than weights.
  @assert length(keys(rewards)) >= length(keys(weights))

  # Start the computation, memorising all solutions (one per value of the budget) along the way.
  solutions = Dict{Int, Vector{T}}()
  budget = 0
  while budget <= max_weight
    sol = solve_budgeted_linear(instance.solver, rewards, weights, budget)

    if length(sol) == 0 || sol == [-1]
      # Infeasible!
      for b in budget:max_weight
        solutions[b] = sol
      end
      break
    end

    # Feasible. Fill the dictionary as much as possible for this budget: if the constraint is not tight, the solution
    # is constant for all budgets between the value prescribed in the budget constraint until the obtained budget.
    sol_budget = sum(arm in keys(weights) ? weights[arm] : 0 for arm in keys(rewards) if arm in sol)
    for b in budget:sol_budget
      solutions[b] = sol
    end
    budget = sol_budget + 1
  end

  return solutions
end


text:
# Errors.

using Base.Libc: errno

"""
    e.g. constant_name(32; prefix="POLL") -> [:POLLNVAL]

Look up name for C-constant(s) with value `n`.
"""
function constant_name(n::Integer; prefix="", first=false)
    @nospecialize

    v = get(C.constants, n, Symbol[])
    if prefix != ""
        v = filter(x -> startswith(String(x), prefix), v)
    end
    if length(v) == 0
        string(n)
    elseif first || length(v) == 1
        string(v[1])
    else
        string("Maybe: ", join(("$(n)?" for n in v), ", "))
    end
end
constant_name(n; kw...) = constant_name(Int(n); kw...)

"""
Throw SystemError with `errno` info for failed C call.
"""
systemerror(p, errno::Cint=errno(); kw...) =
    Base.systemerror(p, errno; extrainfo=errname(errno))

errname(n) = constant_name(n; prefix="E")

"""
    @cerr [allow=C.EINTR] f(args...)

If `f` returns -1 throw SystemError with `errno` info.
"""
macro cerr(a, b=nothing)
    ex = b == nothing ? a : b
    allow = b == nothing ? b : a
    @require ex isa Expr && ex.head == :call
    @require allow == nothing || (
             allow isa Expr && allow.head == Symbol("=") &&
             allow.args[1] == :allow)
    f = ex.args[1]
    args = ex.args[2:end]
    r = gensym()
    condition = allow == nothing ?
                :($r == -1) :
                :($r == -1 && errno() ∉ $(allow.args[2]))
    esc(:(begin
        $r = $ex
        $condition && systemerror(dbstring($f, ($(args...),)))
        @assert $r >= -1
        $r
    end))
end

"""
    @czero f(args...)

If `f` returns non-zero throw SystemError with return value as `errno`.
"""
macro cerr0(ex)
    @require ex isa Expr && ex.head == :call
    f = ex.args[1]
    args = ex.args[2:end]
    r = gensym()
    esc(:(begin
        $r = $ex
        $r != 0 && systemerror(dbstring($f, ($(args...),)), $r)
        nothing
    end))
end

# End of file: errors.jl


text:
using PyPlot
if !@isdefined(CGS)
  include("../../../src/CGS.jl")
  using Main.CGS
end
include("../../../src/seager_mr.jl")

include("../../../src/loglinspace.jl")

#function mass_radius2(mass_ratio_total,mstar,smstar,flatchain,logplot)

nplanet = 7
# Give names to planets:
planet_name=["b","c","d","e","f","g","h"]
#cname = ["b","g","r","c","m","y","k"]
cname = ["C0","C1","C2","C3","C4","C5","C6","C7"]

r1 = 0.6
r2 = 1.2
nr = 200
if logplot
 m1 = 0.01
 m2 = 4.0
else
 m1 = 0.20
 m2 = 1.6
end
nm = 600
mrgrid = zeros(nplanet,nr,nm)

rgrid = zeros(nr)
mgrid = zeros(nm)
for i=1:nr
  rgrid[i] = (i-0.5)*(r2-r1)/nr+r1
end
for i=1:nm
  if logplot
    mgrid[i] = 10 .^((i-0.5)*(log10(m2)-log10(m1))/nm)*m1
  else
    mgrid[i] = (i-0.5)*(m2-m1)/nm+m1
  end
end

# Set up plotting window in Matplotlib (called by PyPlot):
#fig,axes = subplots()

# Number of posterior samples to select:
nsamp = 10000000
cmf = zeros(nplanet,nsamp)
#nsamp = 30000
npoints = [100,100,100,100,100,100,1000]
i0 = 20001
# Factor to convert from solar masses to Earth masses:
fac_mass = MSUN/MEARTH
# Factor to convert from solar radii to Earth radii:
fac_rad = RSUN/REARTH

# Average mass of star in solar masses:
#mstar = 0.09 # Mass of star
# Stand deviation of mass of star in solar masses:
#smstar  = 0.01 # Uncertainty on stellar mass

# Loop over planets:
#conf = [0.997, 0.95, 0.683]
conf = [0.95, 0.683]
nconf = length(conf)
level_mr = zeros(nplanet,nconf)
mavg = zeros(nplanet)
merror = zeros(nplanet)
ravg = zeros(nplanet)
rerror = zeros(nplanet)
# Photodynamic size:
pdsize = size(flatchain)[2]
# N-body HMC size:
nhmc = size(mass_ratio_total)[2]
# Use latest markov chain radius ratios:
# Set up vectors for samples:
msamp = zeros(nplanet,nsamp)
rsamp = zeros(nplanet,nsamp)
mstar_samp = zeros(nsamp)
rho_samp = zeros(nsamp)
rstar_samp = zeros(nsamp)
# Loop over number of samples:
for isamp=1:nsamp
  # Draw from posterior distribution:
  mplanet = mass_ratio_total[:,ceil(Int64,rand()*nhmc)]
  # Select a stellar mass:
  mstar_samp[isamp] = mstar+randn()*smstar
  # Compute the mass of planet in Earth mass units:
  msamp[:,isamp] = mplanet .*(mstar_samp[isamp]*fac_mass)
  # Select a density (relative to Solar):
  ipd = ceil(Int64,rand()*pdsize)
  rho_samp[isamp] = flatchain[15,ipd]
  # Compute the radius of the star in solar radii:
  rstar_samp[isamp] = (mstar_samp[isamp]/rho_samp[isamp])^(1/3)
  # Compute the radius of planet in Earth radii units:
  # Randomly select a radius ratio:
  rsamp[:,isamp] = flatchain[1:nplanet,ipd] .*(rstar_samp[isamp]*fac_rad)
  # Compute CMF from Zeng et al. (2016) formula:
  cmf[:,isamp] = (1.07 .-rsamp[:,isamp] .*msamp[:,isamp].^(-1/3.7))./0.21
  # Now compute density distribution:
  for j=1:nplanet
    if msamp[j,isamp] > m1 && msamp[j,isamp] < m2 && rsamp[j,isamp] > r1 && rsamp[j,isamp] < r2
      ir = ceil(Int64,(rsamp[j,isamp]-r1)/(r2-r1)*nr)
      if logplot
        im = ceil(Int64,(log(msamp[j,isamp])-log(m1))/(log(m2)-log(m1))*nm)
      else
        im = ceil(Int64,(msamp[j,isamp]-m1)/(m2-m1)*nm)
      end
      mrgrid[j,ir,im] +=1.0
    end
  end
end
for j=1:nplanet
  # Set levels for each planet:
  mrgrid[j,:,:] /= float(sum(mrgrid[j,:,:]))
  mrgrid_sort = sort(vec(mrgrid[j,:,:]))
  i0 = nm*nr
  mrgrid_cum = mrgrid_sort[i0]
  while mrgrid_cum < conf[1] && i0 > 1
    i0 -= 1
    for k=1:nconf
      if mrgrid_cum < conf[k] && (mrgrid_cum+mrgrid_sort[i0]) > conf[k]
        level_mr[j,k] = 0.5*(mrgrid_sort[i0]+mrgrid_sort[i0+1])
      end
    end
    mrgrid_cum += mrgrid_sort[i0]
  end
  println(j," Mass:   ",mean(msamp[j,:])," ",std(msamp[j,:])," fractional error: ",std(msamp[j,:])/mean(msamp[j,:])*100.0,"%")
  println(j," Radius: ",mean(rsamp[j,:])," ",std(rsamp[j,:])," fractional error: ",std(rsamp[j,:])/mean(rsamp[j,:])*100.0,"%")
  merror[j] = std(msamp[j,:])
  mavg[j] = mean(msamp[j,:])
  rerror[j] = std(rsamp[j,:])
  ravg[j] = mean(rsamp[j,:])
end
# Plot maximum likelihood values:
#mopt = readdlm("optimum_values.txt")
#ropt = [0.08581987985664834,0.08461786201718988,0.060328647716537474,0.07106532992445118, 0.08035101963995364, 0.08570949401350199, 0.057994985982196025, 53.833783736830554]
for j=1:nplanet
#  cs = axes.contour(mgrid,rgrid,mrgrid[j,:,:],levels = level_mr[j,:],colors=cname[j],linewidth=2.0,label=planet_name[j])
  cs = contour(mgrid,rgrid,mrgrid[j,:,:],levels = level_mr[j,:],colors=cname[j],linewidth=2.0)
#  plot(mopt[(j-1)*5+1]*fac_mass*mstar,ropt[j]*(mstar/ropt[8])^(1/3)*fac_rad,"o",label=planet_name[j],color=cname[j])
  plot(mavg[j],ravg[j],"o",label=planet_name[j],color=cname[j])
#  axes[:imshow](mrgrid[j,:,:],alpha=0.1)
end
# Make some plots:
#axes.set_xlabel(L"Mass $[M_\oplus]$")
xlabel(L"Mass $[M_\oplus]$")
#axes.set_ylabel(L"Radius $[R_\oplus]$")
ylabel(L"Radius $[R_\oplus]$")
#axes.axis([m1,m2,r1,r2])
axis([m1,m2,r1,r2])
text(1.0,1.00,"Earth")
text(0.85,0.93,"Venus")
if logplot
  #axes.set_xscale("log")
  set_xscale("log")
end
#mm = logarithmspace(-2.5,0.2,1000)
#plot(mm,seager_mr(mm,"H2O"),label="Water",linewidth=3,linestyle="--")
#plot(mm,seager_mr(mm,"Fe"),label="Iron",linewidth=3,linestyle="--")
#plot(mm,seager_mr(mm,"Earth"),label="Earth",linewidth=3,linestyle="--")
#plot(mm,seager_mr(mm,"MgSiO3"),label="Rock",linewidth=3,linestyle="--")
#axes[:legend](loc = "upper left",fontsize=10)
#axes.legend(loc = "upper left",fontsize=10)
#legend(loc = "upper left",fontsize=10)
#return merror
#end


text:
@testset "1143.longest-common-subsequence.jl" begin
    @test longest_common_subsequence("abcde", "ace") == 3
    @test longest_common_subsequence("abc", "abc") == 3
    @test longest_common_subsequence("abc", "def") == 0
end


text:
 #initialize queens variables and domains

    columns = [1, 2, 3, 4, 5, 6, 7, 8]
   rows = Dict(col => columns for col in columns)
   println("defining columns and rows")

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

mutable struct QueensConstraint
    variables::Vector{Int}

   function QueensConstraint(columns::Vector{Int})
   println("creating constraint")
        variables = columns
        return new(variables)
    end
end

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

#checks if the constraint is satisfied returns boolean
function satisfied(cons::QueensConstraint,assignment)
println("check if satisfied")
     for (q1c,q1r) in assignment # q1c = queen 1 column, q1r = queen 1 row
        for q2c in range(q1c + 1, stop=length(cons.variables)+1) # q2c = queen 2 column
            if q2c in keys(assignment)
                q2r = assignment[q2c] # q2r = queen 2 row

                if q1r == q2r # same row?
                  return false
                end

                if abs(q1r-q2r) == abs(q1c-q2c)  # same diagonal?
                 return false
                end
            end
        end
    end

   return true # no conflict
end

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

#csp struct for
struct CSP

    #the problem's fields
    variables::Vector
    domains::Dict{Int, Vector}
    constraints::Dict{Int, Vector{}}

    #constructor
    function CSP(vars, doms, cons=Dict())
    println("creating csp")
        variables = vars
        domains = doms
        constraints = cons

        for var in vars
            constraints[var] = Vector()
            if (!haskey(domains, var))
                error("Every variable should have a domain assigned to it.")
            end
        end
        return new(variables, domains, constraints)
    end
end

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

function add_constraint(csp:: CSP, constraint::QueensConstraint)

    for var in constraint.variables
    println("adding cons")
        if (!(var in csp.variables))
            error("Variable in constraint not in CSP")
        else
            push!(csp.constraints[var], constraint)
        end
    end
end

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

function consistent(csp:: CSP, variable, assignment)::Bool
println("check consistency")
    for constraint in csp.constraints[variable]
        if (!satisfied(constraint, assignment))
            return false
        end
        return true
    end
end

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

function backtracking_search(csp:: CSP,  assignment=Dict(), path=Dict())::Union{Dict,Nothing}
println("back tracking")
     # assignment is complete if every variable is assigned (our base case)
     if length(assignment) == length(csp.variables)
        return assignment
     end

    unassigned = [v for v in csp.variables if!(v in keys(assignment))]
#=
    # get all variables without assignments
    unassigned::Vector{Int} = []
    for v in csp.variables
        if (!haskey(assignment, v))
            push!(unassigned, v)
        end
    end

=#
    # get the every possible domain value of the first unassigned variable
    first = unassigned[1]
    # println(first)
    # pretty_print(csp.domains)
    for value in csp.domains[first]
        local_assignment = deepcopy(assignment)
        local_assignment[first] = value

        # if we're still consistent, we recurse (continue)
        if consistent(csp, first, local_assignment)
            # forward checking, prune future assignments that will be inconsistent
            for un in unassigned
                ass = deepcopy(local_assignment)
                for (i, val) in enumerate(csp.domains[un])
                    ass[un] = val
                    if un != first
                        if(!consistent(csp, un, ass))
                            deleteat!(csp.domains[un], i)
                        end
                    end
                end
            end

            path[first] = csp.domains
            # println("reduced")
            #out(csp.domains)
            result = backtracking_search(csp, local_assignment)

            #backtrack if nothing is found
            if result !== nothing
                out(path)
                return result
            end
        end
    end
    return nothing
end

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

function out(d::Dict, pre=1)

    for (k,v) in d
        if typeof(v) <: Dict
            s = "$(repr(k)) => "
            println(join(fill(" ", pre)) * s)
            out(v, pre+1+length(s))
        else
            println(join(fill(" ", pre)) * "$(repr(k)) => $(repr(v))")
        end
    end
    nothing
end
#-----------------------------------------------------------------

csp = CSP(columns,rows)
add_constraint(csp, QueensConstraint(columns))
solution = backtracking_search(csp)

if solution != Nothing
    print("solution is: ", solution)
else
    println("No solution found!")
end

max_stars_repo_path: Problem 1.jl

max_issues_repo_path: Problem 1.jl

max_forks_repo_path: Problem 1.jl


text:
#=
Implement a stack that has the following methods:

    push(val), which pushes an element onto the stack
    pop(), which pops off and returns the topmost element of the stack. If there are no elements in the stack, then it should throw an error or return null.
    max(), which returns the maximum value in the stack currently. If there are no elements in the stack, then it should throw an error or return null.

Each method should run in constant time.
=#

using Test
include("Solutions/problem43_implement_stack.jl")

push(2)
push(5)

@test stack == [2, 5]
@test max() == 5
@test pop() == 5
@test stack == [2]
@test pop() == 2
@test stack == []
@test isnothing(max())
@test isnothing(pop())

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
using Test, LinearAlgebra, SparseArrays, Random
using QDLDL
rng = Random.MersenneTwister(2401)

@testset "refactoring checks" begin

  m = 20
  n = 30

  A  = random_psd(n) + Diagonal(rand(n))
  B  = sprandn(m,n,0.2)
  C1 = -Diagonal(rand(m))
  C2 = -Diagonal(rand(m))
  M1 = [A B' ; B C1]
  M2 = [A B' ; B C2]
  M3 = [A B' ; B -pi*I]

  b = randn(m+n)
  F = qdldl(M1,perm = randperm(m+n))

  @test norm(M1\b - F\b) <= 1e-10

  #update the diagonal of M and F and try again
  update_diagonal!(F,(n+1):(m+n),diag(C2))
  @test norm(M2\b - F\b) <= 1e-10

  #update to a scalar and try again
  update_diagonal!(F,(n+1):(m+n),-pi)
  @test norm(M3\b - F\b) <= 1e-10

end


text:
using Brownian
using Base.Test

@test convert(FGN, FBM(0:0.5:10, 0.4)) == FGN(0.757858283255199,0.4)

p = FBM(0:1/2^10:1, 0.4)

# fBm using FFT approach

rand(p)
rand(p, fbm=false)
rand([p, p])

# fBm using Riemann-Liouville approach

rand(p, method=:rl)
rand(p, method=:rl, wtype=:improved)
rand([p, p], method=:rl)
rand([p, p], method=:rl, wtype=:improved)

# fBm using Cholesky approach

rand(p, method=:chol)
rand(p, fbm=false, method=:chol)
rand([p, p], method=:chol)

# Checking function for efficient Cholesky update

p = FBM(0:.1:.1*100, 0.4)
q = FBM(0:.1:.1*101, 0.4)

P = autocov(p)
Q = autocov(q)
c = chol(P)'
chol_update(convert(Array{Float64, 2}, c), Q)


text:
function strip_output_color(cur_string::AbstractString)
  replace(cur_string, r"\e[^m]*m", "")
end


text:
function type1totype0!(data::MIDIFile)
	if data.format != UInt8(1)
		error("Got type $(data.format); expecting type 1")
	end
	if dochannelsconflict(getprogramchangeevents(data))
		error("Conversion failed since different tracks
		patch different instruments to the same channel.
		Use getprogramchangeevents for inspection.")
	end
	for track in data.tracks
		toabsolutetime!(track)
	end
	for track in data.tracks[2:end]
		insertsorted!(data.tracks[1].events, track.events)
	end
	data.tracks = data.tracks[1:1]
	fromabsolutetime!(data.tracks[1])
	data.format = 0
	data
end

function type1totype0(data::MIDIFile)
	newdata = deepcopy(data)
	type1totype0!(newdata)
end

function type0totype1!(data::MIDIFile)
	if data.format != UInt8(0)
		error("Got type $(data.format); expecting type 0")
	end
	toabsolutetime!(data.tracks[1])
	push!(data.tracks, MIDITrack(Array{TrackEvent, 1}()))
	nofchannels = 0
	for event in data.tracks[1].events
		if !isa(event, MIDIEvent)
			push!(data.tracks[2].events, event)
		else
			channelnum = channelnumber(event)
			while channelnum > nofchannels
				push!(data.tracks, MIDITrack(Array{TrackEvent, 1}()))
				nofchannels += 1
			end
			push!(data.tracks[channelnum + 2].events, event)
		end
	end
	shift!(data.tracks)
	trackstoremove = []
	for i in 1:length(data.tracks)
		if length(data.tracks[i].events) == 0
			push!(trackstoremove, i)
		else
			fromabsolutetime!(data.tracks[i])
		end
	end
	deleteat!(data.tracks, trackstoremove)
	data.format = 1
	data
end

function type0totype1(data::MIDIFile)
	newdata = deepcopy(data)
	type0totype1!(newdata)
end

function insertsorted!(events1::Array{TrackEvent, 1},
					   events2::Array{TrackEvent, 1})
	map(x -> insertsorted!(events1, x), events2)
end

function insertsorted!(events1::Array{TrackEvent, 1},
					   event::TrackEvent)
	i = 0
	while i < length(events1) && events1[end - i].dT > event.dT
		i += 1
	end
	insert!(events1, length(events1) - i + 1, event)
end

function toabsolutetime!(track::MIDITrack)
	t = Int(0)
	for event in track.events
		t += event.dT
		event.dT = t
	end
end

function fromabsolutetime!(track::MIDITrack)
	t0 = t1 = Int(0)
	for event in track.events
		t1 = event.dT
		event.dT = event.dT - t0
		if event.dT < 0
			error("Negative deltas are not allowed. Please reorder your events.")
		end
		t0 = t1
	end
end

function getprogramchangeevents(data::MIDIFile)
	pgevents = []
	t = 0
	for track in data.tracks
		t += 1
		i = 0
		for event in track.events
			i += 1
			if event isa ProgramChangeEvent
				push!(pgevents, [t, i, event])
			end
		end
	end
	pgevents
end

function dochannelsconflict(pgevents)
	channels = Dict()
	for pgevent in pgevents
		newkey = pgevent[3].status
		if haskey(channels, newkey)
			if channels[newkey][1] != pgevent[1] &&
				channels[newkey][2] != pgevent[3].data[1]
				return true
			end
		else
			channels[newkey] = (pgevent[1], pgevent[3].data[1])
		end
	end
	false
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1ExternalDocumentation <: SwaggerModel
    description::Any # spec type: Union{ Nothing, String } # spec name: description
    url::Any # spec type: Union{ Nothing, String } # spec name: url

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1ExternalDocumentation(;description=nothing, url=nothing)
        o = new()
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1ExternalDocumentation, Symbol("description"), description)
        setfield!(o, Symbol("description"), description)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1ExternalDocumentation, Symbol("url"), url)
        setfield!(o, Symbol("url"), url)
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1ExternalDocumentation

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1ExternalDocumentation = Dict{Symbol,Symbol}(Symbol("description")=>Symbol("description"), Symbol("url")=>Symbol("url"))
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1ExternalDocumentation = Dict{Symbol,String}(Symbol("description")=>"String", Symbol("url")=>"String")
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1ExternalDocumentation }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1ExternalDocumentation))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1ExternalDocumentation }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1ExternalDocumentation[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1ExternalDocumentation }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1ExternalDocumentation[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1ExternalDocumentation)
    true
end

function validate_property(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1ExternalDocumentation }, name::Symbol, val)
end


text:
# NOTE:
#   Type fields
#     fieldnames(Stan.Sample)
#       num_samples, num_warmup, save_warmup, thin, adapt, algorithm

#     fieldnames(Stan.Hmc)
#       engine, metric, stepsize, stepsize_jitter

#     fieldnames(Stan.Adapt)
#       engaged, gamma, delta, kappa, t0, init_buffer, term_buffer, window

#   Ref
#     http://goedman.github.io/Stan.jl/latest/index.html#Types-1
sample(mf::T, ss::Stan.Sample) where {T<:Function} = sample(mf, ss.num_samples, ss.num_warmup, ss.save_warmup, ss.thin, ss.adapt, ss.alg)
sample(mf::T, num_samples::Int, num_warmup::Int, save_warmup::Bool, thin::Int, ss::Stan.Sample) where{T<:Function} =
  sample(mf, num_samples, num_warmup, save_warmup, thin, ss.adapt, ss.alg)

sample(mf::T, num_samples::Int, num_warmup::Int, save_warmup::Bool, thin::Int, adapt::Stan.Adapt, alg::Stan.Hmc) where {T<:Function} = begin
  if alg.stepsize_jitter != 0
    error("[Turing.sample] Turing does not support adding noise to stepsize yet.")
  end
  if adapt.engaged == false
    if isa(alg.engine, Stan.Static)   # hmc
      stepnum = Int(round(alg.engine.int_time / alg.stepsize))
      sample(mf, HMC(num_samples, alg.stepsize, stepnum); adapt_conf=adapt)
    elseif isa(alg.engine, Stan.Nuts) # error
      error("[Turing.sample] Stan.Nuts cannot be used with adapt.engaged set as false")
    end
  else
    if isa(alg.engine, Stan.Static)   # hmcda
      sample(mf, HMCDA(num_samples, num_warmup, adapt.delta, alg.engine.int_time); adapt_conf=adapt)
    elseif isa(alg.engine, Stan.Nuts) # nuts
      if alg.metric == Stan.dense_e
        sample(mf, NUTS(num_samples, num_warmup, adapt.delta); adapt_conf=adapt)
      else
        error("[Turing.sample] Turing does not support full covariance matrix for pre-conditioning yet.")
      end
    end
  end
end


text:
#Machinery to convert calls to random variables to scalls, e.g.
#normal(a,3) = scall(Normal, (), nocond, a, 3).
#scall has no optional args, so we can exploit type information properly.

as = map(i -> symbol(string("a", i)), 1:10)
#typed_as = map(a -> :($a::Union(Number, Array)), as)
macro entry(Dist, i)
    dist = symbol(string(lowercase(string(Dist))))
    esc(quote
            @sf function $dist($(as[1:i]...); dims::(Int...)=(), condition=nocond, whitened=default_whitened(condition, $Dist), sampler=default_sampler(condition, $Dist))
                scall($Dist, dims, condition, whitened, sampler, $(as[1:i]...))
            end
            $(Expr(:export, symbol(lowercase(string(Dist)))))
        end)
end
@entry(Normal, 2)
@entry(MvNormal, 2)
@entry(Chisq, 1)
@entry(Exponential, 1)
@entry(Bernoulli, 1)
@entry(Cauchy, 2)
@entry(Laplace, 2)
import Base.gamma
@entry(Gamma, 2)
@entry(Beta, 2)
@entry(Dirichlet, 1)
@entry(Categorical, 1)

default_whitened(::NoCond, _) = whitened
default_whitened(_, ::Type) = unwhitened
default_whitened(::NoCond, ::Type{Bernoulli})   = unwhitened
default_whitened(::NoCond, ::Type{Categorical}) = unwhitened

default_sampler(::NoCond, _) = langevinprior
default_sampler(_, ::Type) = conditioned
default_sampler(::NoCond, ::Type{Bernoulli})   = gibbs
default_sampler(::NoCond, ::Type{Categorical}) = gibbs

const VectorParams = Union(MvNormal, Dirichlet, Categorical)
const RealParams   = Union(Normal, LogNormal, Chisq, Exponential, Bernoulli, Cauchy, Laplace, Gamma, Beta)

const AN = Union(Array, Number)

#Conditioning functions.
scall(state::State, D, dims::(Int...), x::Union(Number, Array), ::UnWhitened, ::Conditioned, params...) = begin
    state.counter[end] += 1
    if isend(state, state.counter[end])
        push!(state, Sample(unwhitened, conditioned, x, logpdf(D, x, params...)))
    else
        state[state.counter[end]] = Sample(unwhitened, conditioned, x, logpdf(D, x, params...))
    end
    x
end

#Sampling functions
#Transform N(0, 1), (parameter independent)
@sf lognormal(μ, σ; condition=nocond, dims=()::(Int...)) =
    lognormal(μ, σ, dims, condition)
@sf lognormal(μ, σ, dims::(Int...), cond::NoCond) =
    exp(normal(μ, σ; dims=dims))
@sf lognormal(μ, σ, dims::(Int...), cond) =
    exp(normal(μ, σ; condition=log(cond), dims=dims))

#Transform N(0, 1), (parameter dependent)
macro direct_sampler(Dist, expr)
    name = Dist.args[1]
    params = Dist.args[2:end]
    esc(quote
        @sf scall(::Type{$name}, dims, ::NoCond, ::Whitened, sampler::Sampler, $(params...)) = $expr
    end)
end
@direct_sampler(Normal(μ, σ), μ .+ σ.*white_rand(maxsize(dims, μ, σ), sampler))
@direct_sampler(MvNormal(μ, Σ), μ + sqrtm(Σ)*white_rand((length(μ),), sampler))
@direct_sampler(Chisq(k::Int), squeeze(sum(white_rand(tuple(k,dims...), sampler.value).^2, 1), 1))
    #Only works for multiple samples, not multiple

#Transform cdf
normal_cdf(z) = 0.5*erfc(-z/√2)
@sf scall(D, dims, ::NoCond, ::Whitened, sampler::Sampler, params...) =
    quantile(D, normal_cdf(white_rand(maxsize(D, dims, params...), sampler)), params...)
#Non-differentiable fallback using Distributions.jl
import Distributions.quantile
quantile{D <: UnivariateDistribution}(d::Type{D}, p, params...) = begin
    quantile_(d, p, params...)
end
quantile_{D <: RealParams}(d::Type{D}, p, params::Number...) =
    quantile(d(params...), p)
import Base.broadcast
broadcast(d::DataType, argss...) =
    broadcast!((args...) -> d(args...), Array(d, maxsize(argss...)), argss...)
quantile_{D <: RealParams}(d::Type{D}, p, params...) = begin
    #Check for subtle bug where not enough gaussian random variables are drawn.
    @assert size(p) == maxsize(params...)
    broadcast(quantile, broadcast(d, params...), p)
end
quantile_{D <: VectorParams}(d::Type{D}, p, params...) =
    quantile(d(params...), p)
#Differentiable definitions.
quantile(::Type{Exponential}, p, scale) = -scale.*log1p(-p)
quantile(::Type{Bernoulli}, p, p0) = p.<p0
quantile(::Type{Cauchy}, p, location, scale) = location .- scale.*cospi(p)./sinpi(p)

#laplace_quantile(p::Real, location::Real, scale::Real) =
#    p < 0.5 ? location + scale*log(2.0*p) : location - scale*log(2.0*(1.0-p))
#laplace_quantile(p, location, scale) = broadcast(laplace_quantile, p, location, scale)
#@cdf_transform(laplace(location, scale), laplace_quantile(p, location, scale))
#
#mapfirst(f::Function, ps, args...) = broadcast(p -> f(p, args...), ps)
#categorical_quantile(p::Real, pv::Vector) = begin
#    i = 1
#    v = pv[1]
#    while v < p && i < length(pv)
#        i += 1
#        @inbounds v += pv[i]
#    end
#    i
#end
#@cdf_transform(categorical(pv), mapfirst(categorical_quantile, p, pv))
#

#Directly call Rmath quantile, because fallbacks are super-slow
quantile_gamma(p::Real, shape::Real, scale::Real) = begin
    ccall((:qgamma, "libRmath-julia"),
          Float64, (Float64, Float64, Float64, Int32, Int32),
          p, shape, scale, 1, 0)
end
const delta = 1E-6
import ReverseDiffOverload.diff_con
#@d3(quantile_gamma,
#    if 0.5<x
#        d*(quantile_gamma(x, y, z) - quantile_gamma(x-delta, y, z))/delta
#    else
#        d*(quantile_gamma(x+delta, y, z) - quantile_gamma(x, y, z))/delta
#    end,
#    d*(quantile_gamma(x, y+delta, z) - quantile_gamma(x, y, z))/delta,
#    d*(quantile_gamma(x, y, z+delta) - quantile_gamma(x, y, z))/delta)
quantile_gamma(x, y, z) = broadcast(quantile_gamma, x, y, z)
dx_quantile_gamma(x::Real, y::Real, z::Real) =
    if 0.5<x
        (quantile_gamma(x, y, z) - quantile_gamma(x-delta, y, z))/delta
    else
        (quantile_gamma(x+delta, y, z) - quantile_gamma(x, y, z))/delta
    end
dx_quantile_gamma(x, y, z) = broadcast(dx_quantile_gamma, x, y, z)
dy_quantile_gamma(x, y, z) = (quantile_gamma(x, y+delta, z) - quantile_gamma(x, y, z))/delta
dz_quantile_gamma(x, y, z) = (quantile_gamma(x, y, z+delta) - quantile_gamma(x, y, z))/delta
@d3(quantile_gamma, (res = d.*dx_quantile_gamma(x, y, z); isa(x, Real) ? sum(res) : res),
                    (res = d.*dy_quantile_gamma(x, y, z); isa(y, Real) ? sum(res) : res),
                    (res = d.*dz_quantile_gamma(x, y, z); isa(z, Real) ? sum(res) : res))
testdiff(quantile_gamma, 0.4, 5., 6.)
testdiff(quantile_gamma, 0.6, 5., 6.)
testdiff((x, y, z) -> sum(quantile_gamma(x, y, z)), [0.4, 0.6], 5., 6.)
testdiff((x, y, z) -> sum(quantile_gamma(x, y, z)), 0.4, [5., 6], 6.)

quantile(::Type{Gamma}, args::Real...) = quantile_gamma(args...)
quantile(::Type{Gamma}, args...) = broadcast(quantile_gamma, args...)

#Transform gamma
@direct_sampler(Beta(alpha, beta), begin
    X = gamma(alpha, 1; dims=dims)
    Y = gamma(beta, 1; dims=dims)
    X./(X+Y)
end)
@direct_sampler(Dirichlet(alpha::Vector), begin
    X = gamma(alpha, 1)
    X/sum(X)
end)

#MaxDims reports the dimensionality of the largest incoming parameter.
#Common cases for speed.
maxsize(as::Number...) = ()
maxsize(a::Array)            = size(a)
maxsize(r::Number,  a::Array) = size(a)
maxsize(a::Array,   r::Number) = size(a)
#General case.
maxdims(as::AN...) = max(map(ndims, as)...)
maxsize(as::AN...) = map(i -> max(map(a -> size(a, i), as)...), tuple(1:maxdims(as...)...))
#Include dims argument.
maxsize(dims::(),       as::Number...) = ()
maxsize(dims::(Int...), as::Number...) = dims
maxsize(dims::(),       as::AN...)     = maxsize(as...)
#Deal with distributions with vector arguments.
maxsize{D <: RealParams}(d::Type{D}, dims, as...) = maxsize(dims, as...)
maxsize{D <: VectorParams}(d::Type{D}, dims, as...) = dims

#UnWhitened sample
scall(state::State, d, dims::(Int...), condition::NoCond, ::UnWhitened, sampler::Sampler, params...) = begin
    state.counter[end] += 1
    isend_     = isend(state, state.counter[end])
    if isend_ || isa(state[state.counter[end]], Resample)
        dist = d(params...)
        val = rand(dist, dims...)
        sample = Sample(unwhitened, sampler, val, sum(logpdf(dist, val)))
        if isend_
            push!(state, sample)
        else#if isresample
            state[state.counter[end]] = sample
        end
    end
    #val
    state[state.counter[end]].value
    #state.counter += 1
    #if (state.endifs_before_resume > 0) ||
    #   (state.counter == length(state.trace) + 1)
    #    dist = d(params...)
    #    value = rand(dist, dims...)
    #    insert!(state.trace, state.counter,
    #            Sample(unwhitened, sampler, value, sum(logpdf(dist, value))))
    #elseif (state.trace[state.counter] == resample)
    #    dist = d(params...)
    #    value = rand(dist, dims...)
    #    state.trace[state.counter] = Sample(unwhitened, sampler, value, sum(logpdf(dist, value)))
    #end
    #state.trace[state.counter].logpdf = sum(logpdf(d(params...), state.trace[state.counter].value))
    #@assert size(state.trace[state.counter].value) == dims
    #state.trace[state.counter].value
end


text:
module OrthPolys

using SparseArrays
using LinearAlgebra: dot

import ACE

import ACE: evaluate!, evaluate_d!, fltype, read_dict, write_dict,
            alloc_B, alloc_dB, transform, transform_d, inv_transform,
            ACEBasis, ScalarACEBasis

using ACE.Transforms: DistanceTransform

import Base: ==

export transformed_jacobi

# this is a hack to prevent a weird compiler error that I don't understand
# at all yet
___f___(D::Dict) = (@show D)

function _fcut_(pl, tl, pr, tr, t)
   if (pl > 0 && t < tl) || (pr > 0 && t > tr)
      return zero(t)
   end
   return (t - tl)^pl * (t - tr)^pr
end

function _fcut_d_(pl, tl, pr, tr, t)
   if (pl > 0 && t < tl) || (pr > 0 && t > tr)
      return zero(t)
   end
   df = 0.0
   if pl > 0; df += pl * (t - tl)^(pl-1) * (t-tr )^pr ; end
   if pr  > 0; df += pr  * (t -  tr)^(pr -1) * (t-tl)^pl; end
   return df
end

@doc raw"""
`OrthPolyBasis:` defined a basis of orthonormal polynomials in terms of the
recursion coefficients. What is slightly unusual is that the polynomials have
an "envelope". This results in the recursion
```math
\begin{aligned}
   J_1(x) &= A_1 (x - x_l)^{p_l} (x - x_r)^{p_r} \\
   J_2 &= (A_2 x + B_2) J_1(x) \\
   J_{n} &= (A_n x + B_n) J_{n-1}(x) + C_n J_{n-2}(x)
\end{aligned}
```
Orthogonality is achieved with respect to a user-specified distribution, which
can be either continuous or discrete.

TODO: say more on the distribution!
"""
struct OrthPolyBasis{T} <: ScalarACEBasis
   # ----------------- the parameters for the cutoff function
   pl::Int        # cutoff power left
   tl::T          # cutoff left (transformed variable)
   pr::Int        # cutoff power right
   tr::T          # cutoff right (transformed variable)
   # ----------------- the recursion coefficients
   A::Vector{T}
   B::Vector{T}
   C::Vector{T}
   # ----------------- used only for construction ...
   #                   but useful to have since it defines the notion of orth.
   tdf::Vector{T}
   ww::Vector{T}
end

fltype(P::OrthPolyBasis{T}) where {T} = T
Base.length(P::OrthPolyBasis) = length(P.A)

==(J1::OrthPolyBasis, J2::OrthPolyBasis) =
      all( getfield(J1, sym) == getfield(J2, sym)
           for sym in (:pr, :tr, :pl, :tl, :A, :B, :C) )

write_dict(J::OrthPolyBasis{T}) where {T} = Dict(
      "__id__" => "ACE_OrthPolyBasis",
      "T" => write_dict(T),
      "pr" => J.pr,
      "tr" => J.tr,
      "pl" => J.pl,
      "tl" => J.tl,
      "A" => J.A,
      "B" => J.B,
      "C" => J.C
   )

OrthPolyBasis(D::Dict, T=read_dict(D["T"])) =
   OrthPolyBasis(
      D["pl"], D["tl"], D["pr"], D["tr"],
      Vector{T}(D["A"]), Vector{T}(D["B"]), Vector{T}(D["C"]),
      T[], T[]
   )

read_dict(::Val{:ACE_OrthPolyBasis}, D::Dict) = OrthPolyBasis(D)

# rand applied to a J will return a random transformed distance drawn from
# the measure w.r.t. which the polynomials were constructed.
# TODO: allow non-constant weights!
function ACE.rand_radial(J::OrthPolyBasis)
   @assert maximum(abs, diff(J.ww)) == 0
   return rand(J.tdf)
end

function OrthPolyBasis(N::Integer,
                       pcut::Integer,
                       tcut::T,
                       pin::Integer,
                       tin::T,
                       tdf::AbstractVector{T},
                       ww::AbstractVector{T} = ones(T, length(tdf))
                       ) where {T <: AbstractFloat}
   @assert pcut >= 0  && pin >= 0
   @assert N > 2

   if tcut < tin
      tl, tr = tcut, tin
      pl, pr = pcut, pin
   else
      tl, tr = tin, tcut
      pl, pr = pin, pcut
   end

   if minimum(tdf) < tl || maximum(tdf) > tr
      @warn("OrthoPolyBasis: t range outside [tl, tr]")
   end

   A = zeros(T, N)
   B = zeros(T, N)
   C = zeros(T, N)

   # normalise the weights s.t. <1, 1> = 1
   ww = ww ./ sum(ww)
   # define inner products
   dotw = (f1, f2) -> dot(f1, ww .* f2)

   # start the iteration
   _J1 = _fcut_.(pl, tl, pr, tr, tdf)
   a = sqrt( dotw(_J1, _J1) )
   A[1] = 1/a
   J1 = A[1] * _J1

   # a J2 = (t - b) J1
   b = dotw(tdf .* J1, J1)
   _J2 = (tdf .- b) .* J1
   a = sqrt( dotw(_J2, _J2) )
   A[2] = 1/a
   B[2] = -b / a
   J2 = (A[2] * tdf .+ B[2]) .* J1

   # keep the last two for the 3-term recursion
   Jprev = J2
   Jpprev = J1

   for n = 3:N
      # a Jn = (t - b) J_{n-1} - c J_{n-2}
      b = dotw(tdf .* Jprev, Jprev)
      c = dotw(tdf .* Jprev, Jpprev)
      _J = (tdf .- b) .* Jprev -c * Jpprev
      a = sqrt( dotw(_J, _J) )
      A[n] = 1/a
      B[n] = - b / a
      C[n] = - c / a
      Jprev, Jpprev = _J / a, Jprev
   end

   return OrthPolyBasis(pl, tl, pr, tr, A, B, C, collect(tdf), collect(ww))
end

alloc_B( J::OrthPolyBasis{T}) where {T} = zeros(T, length(J))
alloc_dB(J::OrthPolyBasis{T}) where {T} = zeros(T, length(J))
# alloc_B( J::OrthPolyBasis{T}, ::Integer) where {T} = zeros(T, length(J))
# alloc_dB(J::OrthPolyBasis{T}, ::Integer) where {T} = zeros(T, length(J))

# # TODO: revisit this to allow type genericity!!!
# alloc_B( J::OrthPolyBasis,
#          ::Union{SVector{3, TX}, AbstractVector{SVector{3, TX}}}) where {TX} =
#    zeros(TX, length(J))
# alloc_dB( J::OrthPolyBasis,
#          ::Union{SVector{3, TX}, AbstractVector{SVector{3, TX}}}) where {TX} =
#    zeros(TX, length(J))

# TODO -> should do a type promotion here???
alloc_B( J::OrthPolyBasis, ::TX) where {TX <: Number} = zeros(TX, length(J))
alloc_dB(J::OrthPolyBasis, ::TX) where {TX <: Number} = zeros(TX, length(J))

evaluate_P1(J::OrthPolyBasis, t) =
   J.A[1] * _fcut_(J.pl, J.tl, J.pr, J.tr, t)

function evaluate!(P, tmp, J::OrthPolyBasis, t; maxn=length(J))
   @assert length(P) >= maxn
   P[1] = evaluate_P1(J, t)
   if maxn == 1; return P; end
   P[2] = (J.A[2] * t + J.B[2]) * P[1]
   if maxn == 2; return P; end
   @inbounds for n = 3:maxn
      P[n] = (J.A[n] * t + J.B[n]) * P[n-1] + J.C[n] * P[n-2]
   end
   return P
end

function evaluate_d!(dP, tmp, J::OrthPolyBasis, t; maxn=length(J))
   @assert maxn <= length(dP)

   P1 = J.A[1] * _fcut_(J.pl, J.tl, J.pr, J.tr, t)
   dP[1] = J.A[1] * _fcut_d_(J.pl, J.tl, J.pr, J.tr, t)
   if maxn == 1; return dP; end

   α = J.A[2] * t + J.B[2]
   P2 = α * P1
   dP[2] = α * dP[1] + J.A[2] * P1
   if maxn == 2; return dP; end

   @inbounds for n = 3:maxn
      α = J.A[n] * t + J.B[n]
      P3 = α * P2 + J.C[n] * P1
      P2, P1 = P3, P2
      dP[n] = α * dP[n-1] + J.C[n] * dP[n-2] + J.A[n] * P1
   end
   return dP
end

function evaluate_ed!(P, dP, tmp, J::OrthPolyBasis, t; maxn=length(J))
   @assert maxn <= min(length(P), length(dP))

   P[1] = J.A[1] * _fcut_(J.pl, J.tl, J.pr, J.tr, t)
   dP[1] = J.A[1] * _fcut_d_(J.pl, J.tl, J.pr, J.tr, t)
   if maxn == 1; return dP; end

   α = J.A[2] * t + J.B[2]
   P[2] = α * P[1]
   dP[2] = α * dP[1] + J.A[2] * P[1]
   if maxn == 2; return dP; end

   @inbounds for n = 3:maxn
      α = J.A[n] * t + J.B[n]
      P[n] = α * P[n-1] + J.C[n] * P[n-2]
      dP[n] = α * dP[n-1] + J.C[n] * dP[n-2] + J.A[n] * P[n-1]
   end
   return dP
end

"""
`discrete_jacobi(N; pcut=0, tcut=1.0, pin=0, tin=-1.0, Nquad = 1000)`

A utility function to generate a jacobi-type basis
"""
function discrete_jacobi(N; pcut=0, tcut=1.0, pin=0, tin=-1.0, Nquad = 1000)
   tl, tr = minmax(tin, tcut)
   dt = (tr - tl) / Nquad
   tdf = range(tl + dt/2, tr - dt/2, length=Nquad)
   return OrthPolyBasis(N, pcut, tcut, pin, tin, tdf)
end

# ----------------------------------------------------------------
#   Transformed Polynomials Basis
# ----------------------------------------------------------------

struct TransformedPolys{T, TT, TJ} <: ScalarACEBasis
   J::TJ          # the actual basis
   trans::TT      # coordinate transform
   rl::T          # lower bound r
   ru::T          # upper bound r = rcut
end

==(J1::TransformedPolys, J2::TransformedPolys) = (
   (J1.J == J2.J) &&
   (J1.trans == J2.trans) &&
   (J1.rl == J2.rl) &&
   (J1.ru == J2.ru) )

TransformedPolys(J, trans, rl, ru) =
   TransformedPolys(J, trans, rl, ru)

write_dict(J::TransformedPolys) = Dict(
      "__id__" => "ACE_TransformedPolys",
      "J" => write_dict(J.J),
      "rl" => J.rl,
      "ru" => J.ru,
      "trans" => write_dict(J.trans)
   )

TransformedPolys(D::Dict) =
   TransformedPolys(
      read_dict(D["J"]),
      read_dict(D["trans"]),
      D["rl"],
      D["ru"]
   )

read_dict(::Val{:SHIPs_TransformedPolys}, D::Dict) =
   read_dict(Val{:ACE_TransformedPolys}(), D)

read_dict(::Val{:ACE_TransformedPolys}, D::Dict) = TransformedPolys(D)

Base.length(J::TransformedPolys) = length(J.J)
fltype(P::TransformedPolys{T}) where {T} = T

function ACE.rand_radial(J::TransformedPolys)
   t = ACE.rand_radial(J.J)
   return inv_transform(J.trans, t)
end

cutoff(J::TransformedPolys) = J.ru

alloc_B( J::TransformedPolys, args...) = alloc_B(J.J, args...)
alloc_dB(J::TransformedPolys) = alloc_dB(J.J)
alloc_dB(J::TransformedPolys, N::Integer) = alloc_dB(J.J)

function evaluate!(P, tmp, J::TransformedPolys, r; maxn=length(J))
   # transform coordinates
   t = transform(J.trans, r)
   # evaluate the actual polynomials
   evaluate!(P, nothing, J.J, t; maxn=maxn)
   return P
end

function evaluate_d!(dP, tmp, J::TransformedPolys, r; maxn=length(J))
   # transform coordinates
   t = transform(J.trans, r)
   dt = transform_d(J.trans, r)
   # evaluate the actual Jacobi polynomials + derivatives w.r.t. x
   evaluate_d!(dP, nothing, J.J, t, maxn=maxn)
   @. dP *= dt
   return dP
end

function evaluate_ed!(P, dP, tmp, J::TransformedPolys, r; maxn=length(J))
   # transform coordinates
   t = transform(J.trans, r)
   dt = transform_d(J.trans, r)
   # evaluate the actual Jacobi polynomials + derivatives w.r.t. x
   evaluate_d!(P, dP, nothing, J.J, t, maxn=maxn)
   @. dP *= dt
   return dP
end

"""
`transformed_jacobi(maxdeg, trans, rcut, rin = 0.0; kwargs...)` : construct
a `TransformPolys` basis with an inner polynomial basis of `OrthPolys` type.

* `maxdeg` : maximum degree
* `trans` : distance transform; normally `PolyTransform(...)`
* `rin, rcut` : inner and outer cutoff

**Keyword arguments:**

* `pcut = 2` : cutoff parameter
* `pin = 0` : inner cutoff parameter
* `Nquad = 1000` : number of quadrature points
"""
function transformed_jacobi(maxdeg::Integer,
                            trans::DistanceTransform,
                            rcut::Real, rin::Real = 0.0;
                            kwargs...)
   J =  discrete_jacobi(maxdeg; tcut = transform(trans, rcut),
                                tin = transform(trans, rin),
                                pcut = 2,
                                kwargs...)
   return TransformedPolys(J, trans, rin, rcut)
end

# ------------- MORE FUNCTIONALITY -------------

# include("products.jl");

end


text:
@enum Frame IncidentFrame TargetFrame

@enum Orientation FixedOrientation RandomOrientation

"""
A spherical particle.

$(FIELDS)
"""
Base.@kwdef struct Sphere
    "Radius of the sphere."
    r::Float64

    "x coordinate of the sphere."
    x::Float64

    "y coordinate of the sphere."
    y::Float64

    "z coordinate of the sphere."
    z::Float64

    "Complex refractive index of the sphere. Default is `1.0`."
    m::ComplexF64 = 1.0

    "Complex chiral factor of the sphere. Default is `0.0`."
    β::ComplexF64 = 0.0

    "The T matrix file specifying the scattering properties of the sphere. Default is `\"\"` (empty)."
    t_matrix::String = ""
end

"""
A horizontally infinite medium layer with certain thickness.

$(FIELDS)
"""
Base.@kwdef struct Layer
    "Thickness of the layer. Default is `0.0`."
    d::Float64 = 0.0

    "Complex refractive index of the layer. Default is `1.0`."
    m::ComplexF64 = 1.0

    "Complex chiral factor of the layer. Default is `0.0`."
    β::ComplexF64 = 0.0
end

"""
Configuration for an STMM run.

$(FIELDS)
"""
Base.@kwdef struct STMMConfig
    "The folder for storing intermediate results. Default is `\"\"`, which means a temporary folder will be created inside the current directory and be used as the working directory."
    working_directory::String = ""

    "The number of processes to be initiated. Default is `1`."
    number_processors::Int = 1

    "Name of the output file. Default is `\"mstm.out\"`."
    output_file::String = "mstm.out"

    "Name of the file for storing run information. Default is `\"mstm.run\"`."
    run_print_file::String = "mstm.run"

    "[**MSTM v3**] Rotate the target about the target coordinate frame using a z-y-z specified Euler angle. Default is `(0.0, 0.0, 0.0)`."
    target_euler_angles_deg::NTuple{3,Float64} = (0.0, 0.0, 0.0)

    """
    [**MSTM v3**, **MSTM v4**] Convergence criterion for determining the truncation order for the wave function expansions for each sphere. Default is `1e-6`.

    > Can also be set as a negative integer `-L` to force all sphere expansions to be truncated at `L`-th order.
    """
    mie_epsilon::Float64 = 1e-6

    "[**MSTM v3**, **MSTM v4**] Convergence criterion for estimating the maximum order of the T matrix for the cluster. Default is `1e-6`."
    translation_epsilon::Float64 = 1e-6

    "[**MSTM v3**, **MSTM v4**] Error criterion for the biconjugate gradient iterative solver. Default is `1e-8`."
    solution_epsilon::Float64 = 1e-8

    "[**MSTM v4**] The maximum truncation degree of the T matrix. This is needed to mostly to keep the code from crashing due to memory allocation limits. Default is `120`."
    max_t_matrix_order::Int = 120

    """
    [**MSTM v3**] Convergence criterion for T matrix solution. Default is `1e-7`.

    > MSTM v4 should also support this option, but it is not yet implemented.
    """
    t_matrix_epsilon::Float64 = 1e-7

    "[**MSTM v3**, **MSTM v4**] The maximum number of iterations attempted to reach a solution. Default is `5000`."
    max_iterations::Int64 = 5000

    """
    [**MSTM v3**, **MSTM v4**] `= false`, the near field translation matrices are calculated during iteration and
    not stored in memory; `= true`, near field matrices are calculated prior to the solution and stored in memory.
    Storing the matrices in memory results in some improvement in execution time, yet the amount of memory used
    vs. the resources available is not checked. Use this option at your own peril. Default is `false`.
    """
    store_translation_matrix::Bool = false

    """
    [**MSTM v3**] The number of processors to be used on parallel platform runs to calculate the
    expansion coefficients for the random orientation scattering matrix. The actual number of processors
    used during this step will be the minimum of the total number of processors for the run (the MPI size)
    and sm_number_processors. Default is `1`.
    """
    sm_number_processors::Int64 = 1

    """
    [**MSTM v3**] The threshold for the translation from near-field calculation to far-field calculation.

    - A value larger than the maximum sphere-to-sphere ``kr`` (i.e., a sufficiently large number, say 1e6) causes all spheres to belong to the near field set.
    - A value smaller than the minimum sphere-to-sphere ``kr`` (i.e., zero) causes all spheres to belong to the far field set.
    - A negative value enables the automatic criterion:

    ```math
    kr_{i-j}\\le\\left[\\frac{1}{2}(L_i+L_j)^2\\right]
    ```

    Default is `1e6`, which means all spheres are considered to be within each other's near field.
    """
    near_field_translation_distance::Float64 = 1e6

    """
    [**MSTM v3**] The maximum number of BCGM iterations for a given correction stage. Default is `20`.
    Only works when there are some sphere pairs belonging to the far field set.
    """
    iterations_per_correction::Int64 = 20

    "[**MSTM v3**] The minimum scattering polar angle. Default is `0.0`."
    θ_min::Float64 = 0.0

    "[**MSTM v3**] The maximum scattering polar angle. Default is `0.0`."
    θ_max::Float64 = 180.0

    "[**MSTM v3**] The minimum scattering azimuthal angle. Default is `0.0`."
    ϕ_min::Float64 = 0.0

    "[**MSTM v3**] The maximum scattering azimuthal angle. Default is `0.0`."
    ϕ_max::Float64 = 0.0

    "[**MSTM v3**, **MSTM v4**] The interval of θ for scattering matrix calculation. Default is `1.0`."
    Δθ::Float64 = 1.0

    """
    [**MSTM v3**, **MSTM v4**] This selects the reference coordinate frame upon which the scattering angles
    are based.

    - `IncidentFrame`, the frame is based on the incident field so that `θ = 0` is the forward scattering direction (i.e., the direction of incident field propagation) and `θ = β, φ = 180` is the target z axis direction.
    - `TargetFrame`, the frame is based on the target frame so that `θ = 0`` corresponds to the target z axis and `θ = β, φ = α`` is the incident field direction .false. Random orientation calculations force this option.
    """
    frame::Frame = TargetFrame

    "[**MSTM v3**, **MSTM v4**] Whether or not to normalize the scattering matrix. Default is `true`."
    normalize_scattering_matrix::Bool = true

    "[**MSTM v3**, **MSTM v4**] Whether or not to take azimuthal average for fixed orientation calculations. Default is `false`."
    azimuthal_average::Bool = false

    """
    [**MSTM v3**, **MSTM v4**] Dimensionless inverse width CB, at the focal point, of an incident Gaussian
    profile beam. Setting to `0`` selects plane wave incidence. The localized approximation used to
    represent the Gaussian beam is accurate for CB ≤ 0.2. Default is `0.0` (plane wave).
    """
    gaussian_beam_constant::Float64 = 0.0

    """
    [**MSTM v3**, **MSTM v4**] x, y, and z coordinates of the focal point of the Gaussian beam, relative to
    the origin defined by `spheres`. Default is `(0.0, 0.0, 0.0)`.
    """
    gaussian_beam_focal_point::NTuple{3,Float64} = (0.0, 0.0, 0.0)

    """
    [**MSTM v3**, **MSTM v4**]

    - Use `FixedOrientation` for a fixed orientation calculation.
    - Use `RandomOrientation` for a random orientation calculation.

    Default is `FixedOrientation`.
    """
    orientation::Orientation = FixedOrientation

    "[**MSTM v4**] Use Monte Carlo integration instead of analytical average. Default is `false`."
    use_monte_carlo_integration::Bool = false

    "[**MSTM v4**] Number of directions used for Monte Carlo integration. Works only when `use_monte_carlo_integration = true`. Default is `100`."
    number_incident_directions::Int = 100

    "The azimuthal angle (in degrees) of the incident wave. Default is `0.0`."
    α::Float64 = 0.0

    "The polar angle (in degrees) of the incident wave. Default is `0.0`."
    β::Float64 = 0.0

    "[**MSTM v4**] Whether or not to calculate the scattering matrix. Default is `true`."
    calculate_scattering_matrix::Bool = true

    """
    [**MSTM v3**] Whether or not to calculate the scattering oefficients. Default is `true`.
    When `= false`, the scattering coefficients will be read from the file specified by
    `scattering_coefficient_file`.
    """
    calculate_scattering_coefficients::Bool = true

    """
    [**MSTM v3**] File name for the scattering coefficient file. Default is `\"mstm.sca\"`.

    - When `calculate_scattering_coefficients = false`, the scattering coefficients will be read from this file.
    - When `calculate_scattering_coefficients = true`, the scattering coefficients will be written to this file.
    """
    scattering_coefficient_file::String = "mstm.scat"

    "[**MSTM v4**] Whether or not to use the FFT-based accelerated algorithm. Default is `false`."
    use_fft_translation::Bool = false

    """
    [**MSTM v4**]
    This determines the grid size ``d`` of the cubic lattice used in the algorithm. The
    formula is:

    ```math
    d=\\left(\\frac{V_S}{cN_{s,ext}}\\right)^{1/3}
    ```

    where ``V_S`` is the dimensionless total volume of the external spheres and ``c`` is `cell_volume_fraction`.
    The optimum value is usually where the total number of nodes is close to the total number of external
    spheres NS,ext. Setting `cell_volume_fraction = 0` will automate the setting of d using an estimate of
    the sphere volume fraction. Default is `0.0`. Works only when `use_fft_translation = true`.
    """
    cell_volume_fraction::Float64 = 0.0

    """
    [**MSTM v4**]
    Truncation degree for the node–based expansions of the scattered field. Set
    `node_order = −L`, where `L` is a positive integer, will set the node order to `ceil(d) + L`.
    Default is `-1`. Works only when `use_fft_translation = true`.
    """
    node_order::Int = -1

    """
    [**MSTM v4**] `= 0` will have the scattering matrix printed at a range of polar angles
    (θ), with angular increment given in degrees by `Δθ`. The scattering plane
    coincides with the incident azimuth plane. The limits on θ depend on the number of plane boundaries
    NB. scattering_map_model= 1 will print the scattering matrix in a pair (cosθ > 0, < 0) of 2D
    arrays, with coordinates `kx = sinθcosφ`, `ky = sinθsinφ`. The number of points on each axis is given
    by `scattering_map_dimension`. Default is `0`, and random orientation calculations force `0`.
    """
    scattering_map_model::Int = 0

    """
    [**MSTM v4**] The number of points on each axis. Default is `15`. Works only when `scattering_map_model = 1`.
    """
    scattering_map_dimension::Int = 15

    "[**MSTM v3**] Whether or not to print each iteration's error in the run print file. Default is `true`."
    track_iterations::Bool = true

    "[**MSTM v3**, **MSTM v4**] Whether or not to calculate the near field. Default is `false`."
    calculate_near_field::Bool = false

    "[**MSTM v3**] Coordinate of the near field plane. `1` = y-z, `2` = z-x, `3` = x-y. Default is `1`."
    near_field_plane_coord::Int64 = 1

    """
    [**MSTM v3**] Position of the near field plane.

    - When `near_field_plane_coord = 1`, this is the x-coordinate.
    - When `near_field_plane_coord = 2`, this is the y-coordinate.
    - When `near_field_plane_coord = 3`, this is the z-coordinate.

    Default is `0.0`.
    """
    near_field_plane_position::Float64 = 0.0

    """
    [**MSTM v3**] Rectangle region of the near field plane.

    - When passing `(ΔX, ΔY)`, the actual region would be `(X_min-ΔX, X_max+ΔX) x (Y_min-ΔY, Y_max+ΔY)`, where `X_min`, `X_max`,`Y_min` and `Y_max` are automatically calculated.
    - When passing `(X_min, Y_min, X_max, Y_max)`, the region would be used as it is.

    Default is `(0.0, 0.0)`.
    """
    near_field_plane_vertices::Union{NTuple{2,Float64},NTuple{4,Float64}} = (0.0, 0.0)

    "[**MSTM v4**] The minimum x coordinate of the near field calculation region. Default is `0.0`."
    near_field_x_min::Float64 = 0.0

    "[**MSTM v4**] The maximum x coordinate of the near field calculation region. Default is `0.0`."
    near_field_x_max::Float64 = 0.0

    "[**MSTM v4**] The minimum y coordinate of the near field calculation region. Default is `0.0`."
    near_field_y_min::Float64 = 0.0

    "[**MSTM v4**] The maximum y coordinate of the near field calculation region. Default is `0.0`."
    near_field_y_max::Float64 = 0.0

    "[**MSTM v4**] The minimum z coordinate of the near field calculation region. Default is `0.0`."
    near_field_z_min::Float64 = 0.0

    "[**MSTM v4**] The maximum z coordinate of the near field calculation region. Default is `0.0`."
    near_field_z_max::Float64 = 0.0

    "[**MSTM v3**, **MSTM v4**] The spacial step size of the grid points. Default is `0.1`."
    near_field_step_size::Float64 = 0.1

    """
    [**MSTM v4**] Controls where the incident field appears in the calculation results; `= 1`
    is the standard model, where the external field is = scattered + incident and the field
    inside the spheres is calculated from the internal field expansions; `≠ 1` has the external
    field due solely to the scattered field, and the field inside the particles is now
    internal-incident. Default is `1`.
    """
    near_field_calculation_model::Int = 1

    """
    [**MSTM v3**] A specific polarization state of the incident field is needed to calculate the near
    field. The field is taken to be linearly polarized, with a polarization angle of `γ` relative to the k–z
    plane. When `β = 0`, `γ` becomes the azimuth angle of the incident electric field vector relative to the
    cluster coordinate system. Default is `0.0`.
    """
    polarization_angle_deg::Float64 = 0.0

    """
    [**MSTM v4**] Whether or not to enable an accelerated algorithm for calculation of field value.
    Default is `true`.
    """
    store_surface_vector::Bool = true

    """
    [**MSTM v4**] Size of the grid used for re–expansion of fields resulting from plane
    boundary and periodic lattice interactions. Default is `5.0`. Experiment with this
    if your calculations are taking too long. Only works when `store_surface_vector = true`.
    """
    near_field_expansion_spacing::Float64 = 5.0

    """
    [**MSTM v4**] Expansion truncation order used for the grid re–expansion of the secondary
    field, per the above option. Should scale with near_field_expansion_spacing following a
    Mie–type criteria. The bigger the value, the more accurate the results, but the slower the
    calculations and the greater the memory demands. Default is `10`. Only works when
    `store_surface_vector = true`.
    """
    near_field_expansion_order::Int = 10

    "File name for the near field output. Default is `\"mstm.nf\"`."
    near_field_output_file::String = "mstm.nf"

    """
    [**MSTM v3**] The incident field component for the near field calculations -- for either the plane
    wave or Gaussian beam models –- is calculated using a single, regular VSWF expansion centered about
    the beam focal point. The plane wave epsilon is a convergence criterion for this expansion. Default is
    `1e-4`.
    """
    plane_wave_epsilon::Float64 = 1e-4

    """
    [**MSTM v3**] Option to calculate the T matrix or use an exisiting T matrix file. Default is `true`,
    which calculates the T matrix.
    """
    calculate_t_matrix::Bool = true

    """
    [**MSTM v3**]

    - When `calculate_t_matrix = true`, this is the file name for the T matrix output.
    - When `calculate_t_matrix = false`, this is the file name for reading the T matrix from.

    Default is `\"mstm.tm\"`.
    """
    t_matrix_file::String = "mstm.tm"

    """
    [**MSTM v3**, **MSTM v4**] Medium layers.

    - For **MSTM v3**, only the refractive index and chiral factor of the first layer will be used. Thickness is ignored.
    - For **MSTM v4**, the first layer locates at `z=0` and extends to `-∞`. The last layer extends to `+∞`. The thickness of the first and the last layer is ignored.
    """
    layers::Vector{Layer} = [Layer()]

    "[**MSTM v3**, **MSTM v4**] Spheres."
    spheres::Vector{Sphere} = Sphere[]

    "[**MSTM v4**] When `periodic = true`, the spheres will be considered to construct a periodic lattice in the x-y plane. Note that periodic lattice is incompatible with random orientation."
    periodic::Bool = false

    "[**MSTM v4**] When `periodic = true`, this option sets the size of the periodic lattice. The cell size will be automatically enlarged if it cannot contain all the spheres."
    cell_size::Tuple{Float64,Float64} = (0.0, 0.0)
end

"""
Helper function for building a configuration on basis of an existing configuration.

$(SIGNATURES)

"""
function STMMConfig(original::STMMConfig; kwargs...)
    dict = Dict(key => getfield(original, key) for key in propertynames(original))
    for (key, value) in kwargs
        dict[key] = value
    end
    return STMMConfig(; dict...)
end

export Orientation, FixedOrientation, RandomOrientation
export Frame, IncidentFrame, TargetFrame
export Sphere, Layer, STMMConfig


text:
columnbuttons = Observable{Any}(dom"div"())
plt = Observable{Any}()

function makebuttons(ph)
    prop = propertynames(ph)[5:end-3] #Removes name,x,y,z and ux,uy,uz
    propnm = [s for s=string.(prop)]
    buttons = button.(propnm)
    for (btn, key) in zip(reverse(buttons), reverse(prop))
        map!(t -> begin
            @manipulate for t0_ms = range(0,dur(seq),5)*1e3
                plot_phantom_map(ph, key; t0=t0_ms, darkmode)
            end
        end
        , plt, btn)
    end
    dom"div"(hbox(buttons))
end

map!(makebuttons, columnbuttons, pha_obs)

ui = dom"div"(vbox(columnbuttons, plt))
content!(w, "div#content", ui)


text:
@testset "Dual numbers" begin
  f(x,y) = x[1]^2+y[1]^2
  g(param) = Optim.minimum(optimize(x->f(x,param), -1,1))
  @test ForwardDiff.gradient(g, [1.0]) == [2.0]
end


text:
# Direct minimization of the energy

using Optim
using LineSearches

# This is all a bit annoying because our ψ is represented as ψ[k][G,n], and Optim accepts
# only dense arrays. We do a bit of back and forth using custom `pack` (ours -> optim's) and
# `unpack` (optim's -> ours) functions

# Orbitals inside each kblock must be kept orthogonal: the
# project_tangent and retract work per kblock
struct DMManifold <: Optim.Manifold
    Nk::Int
    unpack::Function
end
function Optim.project_tangent!(m::DMManifold, g, x)
    g_unpack = m.unpack(g)
    x_unpack = m.unpack(x)
    for ik = 1:m.Nk
        Optim.project_tangent!(Optim.Stiefel(),
                               g_unpack[ik],
                               x_unpack[ik])
    end
    g
end
function Optim.retract!(m::DMManifold, x)
    x_unpack = m.unpack(x)
    for ik = 1:m.Nk
        Optim.retract!(Optim.Stiefel(), x_unpack[ik])
    end
    x
end

# Array of preconditioners
struct DMPreconditioner
    Nk::Int
    Pks::Vector # Pks[ik] is the preconditioner for k-point ik
    unpack::Function
end
function LinearAlgebra.ldiv!(p, P::DMPreconditioner, d)
    p_unpack = P.unpack(p)
    d_unpack = P.unpack(d)
    for ik = 1:P.Nk
        ldiv!(p_unpack[ik], P.Pks[ik], d_unpack[ik])
    end
    p
end
function LinearAlgebra.dot(x, P::DMPreconditioner, y)
    x_unpack = P.unpack(x)
    y_unpack = P.unpack(y)
    sum(dot(x_unpack[ik], P.Pks[ik], y_unpack[ik])
        for ik = 1:P.Nk)
end
function precondprep!(P::DMPreconditioner, x)
    x_unpack = P.unpack(x)
    for ik = 1:P.Nk
        precondprep!(P.Pks[ik], x_unpack[ik])
    end
    P
end

"""
Computes the ground state by direct minimization. `kwargs...` are
passed to `Optim.Options()`. Note that the resulting ψ are not
necessarily eigenvectors of the Hamiltonian.
"""
direct_minimization(basis::PlaneWaveBasis; kwargs...) = direct_minimization(basis, nothing; kwargs...)
function direct_minimization(basis::PlaneWaveBasis{T}, ψ0;
                             prec_type=PreconditionerTPA,
                             optim_solver=Optim.LBFGS, tol=1e-6, kwargs...) where T
    if mpi_nprocs() > 1
        # need synchronization in Optim
        error("Direct minimization with MPI is not supported yet")
    end
    model = basis.model
    @assert model.temperature == 0  # temperature is not yet supported
    filled_occ = filled_occupation(model)
    n_spin = model.n_spin_components
    n_bands = div(model.n_electrons, n_spin * filled_occ, RoundUp)
    Nk = length(basis.kpoints)

    if ψ0 === nothing
        ψ0 = [ortho_qr(randn(Complex{T}, length(G_vectors(basis, kpt)), n_bands))
              for kpt in basis.kpoints]
    end
    occupation = [filled_occ * ones(T, n_bands) for ik = 1:Nk]

    # we need to copy the reinterpret array here to not raise errors in Optim.jl
    # TODO raise this issue in Optim.jl
    pack(ψ) = copy(reinterpret_real(pack_ψ(ψ)))
    unpack(x) = unpack_ψ(reinterpret_complex(x), size.(ψ0))

    # this will get updated along the iterations
    H = nothing
    energies = nothing
    ρ = nothing

    # computes energies and gradients
    function fg!(E, G, ψ)
        ψ = unpack(ψ)
        ρ = compute_density(basis, ψ, occupation)
        energies, H = energy_hamiltonian(basis, ψ, occupation; ρ=ρ)

        # The energy has terms like occ * <ψ|H|ψ>, so the gradient is 2occ Hψ
        if G !== nothing
            G = unpack(G)
            for ik = 1:Nk
                mul!(G[ik], H.blocks[ik], ψ[ik])
                G[ik] .*= 2*filled_occ
            end
        end
        energies.total
    end

    manif = DMManifold(Nk, unpack)

    Pks = [prec_type(basis, kpt) for kpt in basis.kpoints]
    P = DMPreconditioner(Nk, Pks, unpack)

    kwdict = Dict(kwargs)
    optim_options = Optim.Options(; allow_f_increases=true, show_trace=true,
                                  x_tol=pop!(kwdict, :x_tol, tol),
                                  f_tol=pop!(kwdict, :f_tol, -1),
                                  g_tol=pop!(kwdict, :g_tol, -1), kwdict...)
    res = Optim.optimize(Optim.only_fg!(fg!), pack(ψ0),
                         optim_solver(P=P, precondprep=precondprep!, manifold=manif,
                                      linesearch=LineSearches.BackTracking()),
                         optim_options)
    # return copy to ensure we have a plain array
    ψ = deepcopy(unpack(res.minimizer))

    # Final Rayleigh-Ritz (not strictly necessary, but sometimes useful)
    eigenvalues = []
    for ik = 1:Nk
        Hψk = H.blocks[ik] * ψ[ik]
        F = eigen(Hermitian(ψ[ik]'Hψk))
        push!(eigenvalues, F.values)
        ψ[ik] .= ψ[ik] * F.vectors
    end

    εF = nothing  # does not necessarily make sense here, as the
                  # Aufbau property might not even be true

    # We rely on the fact that the last point where fg! was called is the minimizer to
    # avoid recomputing at ψ
    (ham=H, basis=basis, energies=energies, converged=true,
     ρ=ρ, ψ=ψ, eigenvalues=eigenvalues, occupation=occupation, εF=εF, optim_res=res)
end


text:
__precompile__()
module SpecialFunctions

using Compat

if VERSION >= v"0.7.0-DEV.1760"
    depsfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
    if isfile(depsfile)
        include(depsfile)
    else
        error("SpecialFunctions is not properly installed. Please run " *
              "Pkg.build(\"SpecialFunctions\") and restart Julia.")
    end
else
    using Base.Math: openspecfun
end

if isdefined(Base, :airyai) && VERSION < v"0.7.0-DEV.986" #22763
    import Base: airyai, airyaix, airyaiprime, airyaiprimex,
                 airybi, airybix, airybiprime, airybiprimex,
                 besselh, besselhx, besseli, besselix, besselj, besselj0, besselj1,
                 besseljx, besselk, besselkx, bessely, bessely0, bessely1, besselyx,
                 hankelh1, hankelh1x, hankelh2, hankelh2x,
                 dawson, erf, erfc, erfcinv, erfcx, erfi, erfinv,
                 eta, digamma, invdigamma, polygamma, trigamma, zeta,
                 # deprecated
                 airy, airyx, airyprime
else
    export
        airyai,
        airyaiprime,
        airybi,
        airybiprime,
        airyaix,
        airyaiprimex,
        airybix,
        airybiprimex,
        besselh,
        besselhx,
        besseli,
        besselix,
        besselj,
        besselj0,
        besselj1,
        besseljx,
        besselk,
        besselkx,
        bessely,
        bessely0,
        bessely1,
        besselyx,
        dawson,
        erf,
        erfc,
        erfcinv,
        erfcx,
        erfi,
        erfinv,
        eta,
        digamma,
        invdigamma,
        polygamma,
        trigamma,
        hankelh1,
        hankelh1x,
        hankelh2,
        hankelh2x,
        zeta
end

export sinint,
       cosint

include("bessel.jl")
include("erf.jl")
include("sincosint.jl")
include("gamma.jl")
include("deprecated.jl")

end # module


text:
__precompile__()

module TMCM3110

export encode_command
function encode_command(m_address, n_command, n_type, n_motor, value)
  m_address = UInt8( m_address % (1<<8) )
  n_command = UInt8( n_command % (1<<8) )
  n_type    = UInt8( n_type % (1<<8)    )
  n_motor   = UInt8( n_motor % (1<<8)   )
  value     = Int32( value )
  values = [ parse(UInt8, bits(value)[ 1+(8*(i-1)) : 8+(8*(i-1)) ] , 2) for i in 1:4]

  checksum  = UInt8( (m_address + n_command + n_type + n_motor + sum(values)) % (1<<8) )

  tmcl_bytes = [m_address, n_command, n_type, n_motor, values..., checksum]

  return tmcl_bytes
end

export decode_reply
function decode_reply(reply)
  r_address  = UInt8( reply[1] )
  m_address  = UInt8( reply[2] )
  r_status   = UInt8( reply[3] )
  Int32(r_status) != 100 ? warn(STATUSCODES[Int(r_status)]) : nothing
  n_command  = UInt8( reply[4] )
  values     = [ UInt8(value) for value in reply[5:8] ]
  r_checksum = UInt8( reply[9] )
  checksum = UInt8( (r_address + m_address + r_status + n_command + sum(values)) % (1<<8) )
  value      = parse( UInt32, "$(bits(values[1]))$(bits(values[2]))$(bits(values[3]))$(bits(values[4]))" , 2 )
  value      = reinterpret(Int32, value)
  return value
end

export STATUSCODES
STATUSCODES = Dict( 100 => "Succesfully executed, no error",
                    101 => "Command loaded into TMCL program EEPROM",
                      1 => "Wrong Checksum",
                      2 => "Invalid command",
                      3 => "Wrong type",
                      4 => "Invalid value",
                      5 => "Configuration EEPROM locked",
                      6 => "Command not available" )

export COMMAND_NUMBERS
COMMAND_NUMBERS = Dict(  1 => "ROR",
                         2 => "ROL",
                         3 => "MST",
                         4 => "MVP",
                         5 => "SAP",
                         6 => "GAP",
                         7 => "STAP",
                         8 => "RSAP",
                         9 => "SGP",
                        10 => "GGP",
                        11 => "STGP",
                        12 => "RSGP",
                        13 => "RFS",
                        14 => "SIO",
                        15 => "GIO",
                        19 => "CALC",
                        20 => "COMP",
                        21 => "JC",
                        22 => "JA",
                        23 => "CSUB",
                        24 => "RSUB",
                        25 => "EI",
                        26 => "DI",
                        27 => "WAIT",
                        28 => "STOP",
                        30 => "SCO",
                        31 => "GCO",
                        32 => "CCO",
                        33 => "CALCX",
                        34 => "AAP",
                        35 => "AGP",
                        37 => "VECT",
                        38 => "RETI",
                        39 => "ACO"  )

export AXIS_PARAMETER
AXIS_PARAMETER = Dict(   0 => "target position",
                         1 => "actual position",
                         2 => "target speed",
                         3 => "actual speed",
                         4 => "max positioning speed",
                         5 => "max acceleration",
                         6 => "abs max current",
                         7 => "standby current",
                         8 => "target pos reached",
                         9 => "ref switch status",
                        10 => "right limit switch status",
                        11 => "left limit switch status",
                        12 => "right limit switch disable",
                        13 => "left limit switch disable",
                       130 => "minimum speed",
                       135 => "actual acceleration",
                       138 => "ramp mode",
                       140 => "microstep resolution",
                       141 => "ref switch tolerance",
                       149 => "soft stop flag",
                       153 => "ramp divisor",
                       154 => "pulse divisor",
                       160 => "step interpolation enable",
                       161 => "double step enable",
                       162 => "chopper blank time",
                       163 => "chopper mode",
                       164 => "chopper hysteresis dec",
                       165 => "chopper hysteresis end",
                       166 => "chopper hysteresis start",
                       167 => "chopper off time",
                       168 => "smartEnergy min current",
                       169 => "smartEnergy current downstep",
                       170 => "smartEnergy hysteresis",
                       171 => "smartEnergy current upstep",
                       172 => "smartEnergy hysteresis start",
                       173 => "stallGuard2 filter enable",
                       174 => "stallGuard2 threshold",
                       175 => "slope control high side",
                       176 => "slope control low side",
                       177 => "short protection disable",
                       178 => "short detection timer",
                       179 => "Vsense",
                       180 => "smartEnergy actual current",
                       181 => "stop on stall",
                       182 => "smartEnergy threshold speed",
                       183 => "smartEnergy slow run current",
                       193 => "ref. search mode",
                       194 => "ref. search speed",
                       195 => "ref. switch speed",
                       196 => "distance end switches",
                       204 => "freewheeling",
                       206 => "actual load value",
                       208 => "TMC262 errorflags",
                       209 => "encoder pos",
                       210 => "encoder prescaler",
                       212 => "encoder max deviation",
                       214 => "power down delay"    )

export INTERRUPT_VECTORS
INTERRUPT_VECTORS = Dict(  0 => "Timer 0",
                           1 => "Timer 1",
                           2 => "Timer 2",
                           3 => "Target position 0 reached",
                           4 => "Target position 1 reached",
                           5 => "Target position 2 reached",
                          15 => "stallGuard2 axis 0",
                          16 => "stallGuard2 axis 1",
                          17 => "stallGuard2 axis 2",
                          21 => "Deviation axis 0",
                          22 => "Deviation axis 1",
                          23 => "Deviation axis 2",
                          27 => "Left stop switch 0",
                          28 => "Right stop switch 0",
                          29 => "Left stop switch 1",
                          30 => "Right stop switch 1",
                          31 => "Left stop switch 2",
                          32 => "Right stop switch 2",
                          39 => "Input change 0",
                          40 => "Input change 1",
                          41 => "Input change 2",
                          42 => "Input change 3",
                          43 => "Input change 4",
                          44 => "Input change 5",
                          45 => "Input change 6",
                          46 => "Input change 7",
                         255 => "Global interrupts" )

end # module