text:
## ZeroOne

function classify(value::T, cutoff::Number) where {T<:Number}
    value >= cutoff ? T(1) : T(0)
end

function classify(value::Number, ::Type{LabelEnc.ZeroOne})
    classify(value, 0.5)
end

function classify(value::Number, lm::LabelEnc.ZeroOne{R}) where {R}
    R(classify(value, lm.cutoff))
end

## Margin

_sign(value::T) where {T} = ifelse(signbit(value), T(-1), T(1))::T

function classify(value::Number, ::Type{LabelEnc.MarginBased})
    _sign(value)
end

function classify(value::Number, lm::LabelEnc.MarginBased{R}) where {R}
    R(_sign(value))
end

## broadcast

function classify!(buffer::AbstractVector, values::AbstractVector, lm)
    buffer .= classify.(values, lm)
    buffer
end

function classify(values::AbstractVector{T}, cutoff::Number) where {T}
    classify.(values, cutoff)::Vector{T}
end

for KIND in (:(LabelEnc.MarginBased), :(LabelEnc.ZeroOne))
    @eval begin
        function classify(values::AbstractVector{T}, ::Type{($KIND)}) where {T}
            classify.(values, ($KIND))::Vector{T}
        end
        function classify(values::AbstractVector{T}, lm::L) where {T,L<:($KIND)}
            classify.(values, lm)::Vector{labeltype(L)}
        end
    end
end

## OneOfK

function classify(values::AbstractVector, ::Type{<:LabelEnc.OneOfK})
    argmax(values)
end

function classify(values::AbstractVector, lm::LabelEnc.OneOfK)
    classify(values, typeof(lm))
end

function classify!(buffer::T,
                   values::AbstractMatrix,
                   lm;
                   obsdim = LearnBase.default_obsdim(values)
                  ) where {T<:AbstractVector}
    classify!(buffer, values, lm, convert(LearnBase.ObsDimension,obsdim))::T
end

function classify(values::AbstractMatrix,
                  lm;
                  obsdim = LearnBase.default_obsdim(values))
    classify(values, lm, convert(LearnBase.ObsDimension,obsdim))
end

function classify!(buffer,
                   values::AbstractMatrix,
                   lm::LabelEnc.OneOfK,
                   obsdim::LearnBase.ObsDimension)
    classify!(buffer, values, typeof(lm), obsdim)
end

function classify(values::AbstractMatrix,
                  lm::LabelEnc.OneOfK,
                  obsdim::LearnBase.ObsDimension)
    classify(values, typeof(lm), obsdim)
end

function classify(values::AbstractMatrix,
                  ::Type{T},
                  ::Union{ObsDim.Last,ObsDim.Constant{2}}
                 ) where {T<:LabelEnc.OneOfK}
    K, N = size(values)
    buffer = Vector{Int}(undef, N)
    classify!(buffer, values, T, ObsDim.Last())
end

function classify!(buffer::AbstractVector,
                   values::AbstractMatrix{R},
                   ::Type{<:LabelEnc.OneOfK},
                   ::Union{ObsDim.Last,ObsDim.Constant{2}}
                  ) where {R<:Number}
    K, N = size(values)
    @assert length(buffer) == N
    @inbounds for n in 1:N
        imax = 0
        tmax = typemin(R)
        for k in 1:K
            tcur = values[k,n]
            if tcur > tmax
                imax = k
                tmax = tcur
            end
        end
        buffer[n] = imax
    end
    buffer
end

function classify(values::AbstractMatrix,
                  ::Type{T},
                  ::ObsDim.First
                 ) where {T<:LabelEnc.OneOfK}
    N, K = size(values)
    buffer = Vector{Int}(undef, N)
    classify!(buffer, values, T, ObsDim.First())
end

function classify!(buffer::AbstractVector,
                   values::AbstractMatrix{R},
                   ::Type{<:LabelEnc.OneOfK},
                   ::ObsDim.First
                  ) where {R<:Number}
    N, K = size(values)
    tmax = fill(typemin(R),N)
    @assert length(buffer) == N
    @inbounds for k in 1:K
        for n in 1:N
            tcur = values[n,k]
            if tcur > tmax[n]
                tmax[n] = tcur
                buffer[n] = k
            end
        end
    end
    buffer
end


text:
using Test
using Tables
using MLJBase
import MLJFlux
using CategoricalArrays
using ColorTypes
using Flux
import Random
import Random.seed!
using Statistics
import StatsBase
using StableRNGs

using ComputationalResources
using ComputationalResources: CPU1, CUDALibs

const RESOURCES = Any[CPU1(), CUDALibs()]
EXCLUDED_RESOURCE_TYPES = Any[]

MLJFlux.gpu_isdead() && push!(EXCLUDED_RESOURCE_TYPES, CUDALibs)

@info "MLJFlux supports these computational resources:\n$RESOURCES"
@info "Current test run to exclude resources with "*
    "these types, as unavailable:\n$EXCLUDED_RESOURCE_TYPES\n"*
    "Excluded tests marked as \"broken\"."

# alternative version of Short builder with no dropout; see
# https://github.com/FluxML/Flux.jl/issues/1372 and
# https://github.com/FluxML/Flux.jl/issues/1372
mutable struct Short2 <: MLJFlux.Builder
    n_hidden::Int     # if zero use geometric mean of input/output
    σ
end
Short2(; n_hidden=0, σ=Flux.sigmoid) = Short2(n_hidden, σ)
function MLJFlux.build(builder::Short2, rng, n, m)
    n_hidden =
        builder.n_hidden == 0 ? round(Int, sqrt(n*m)) : builder.n_hidden
    init = Flux.glorot_uniform(rng)
    return Flux.Chain(
        Flux.Dense(n, n_hidden, builder.σ, init=init),
        Flux.Dense(n_hidden, m, init=init))
end

seed!(123)

include("test_utils.jl")

@testset "penalizers" begin
    include("penalizers.jl")
end

@testset "core" begin
    include("core.jl")
end

@testset "builders" begin
    include("builders.jl")
end

@testset "mlj_model_interface" begin
    include("mlj_model_interface.jl")
end

@testset "regressor" begin
    include("regressor.jl")
end

@testset "classifier" begin
    include("classifier.jl")
end

@testset "image" begin
    include("image.jl")
end

@testset "integration" begin
    include("integration.jl")
end

ext: jl


text:
using Documenter, SialiaCLI

makedocs(;
    modules=[SialiaCLI],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/wardlem/SialiaCLI.jl/blob/{commit}{path}#L{line}",
    sitename="SialiaCLI.jl",
    authors="Mark Wardle, Potient",
    assets=String[],
)

deploydocs(;
    repo="github.com/wardlem/SialiaCLI.jl",
)


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, Distributions, Plots
	Random.seed!(0)
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing5.11"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	N = 10^5
	nMin, nStep, nMax = 10, 10, 200
	nn = Int(nMax/nStep)
	sampleSizes = nMin:nStep:nMax
	trueB = 5
	trueDist = Uniform(-2, trueB)

	MLEest(data) = maximum(data)
	MMest(data)  = mean(data) + sqrt(3)*std(data)
end

# ╔═╡ 5790308c-3af7-11eb-1ebf-b57c8c5e4061
res = Dict{Symbol,Array{Float64}}(
    ((sym) -> sym => Array{Float64}(undef,nn)).(
        [:MSeMLE,:MSeMM, :VarMLE,:VarMM,:BiasMLE,:BiasMM]))

# ╔═╡ 57b8a47c-3af7-11eb-3da1-aba055f63a52
for (i, n) in enumerate(sampleSizes)
    mleEst, mmEst = Array{Float64}(undef, N), Array{Float64}(undef, N)
    for j in 1:N
        sample    = rand(trueDist,n)
        mleEst[j] = MLEest(sample)
        mmEst[j]  = MMest(sample)
    end
    meanMLE, meanMM = mean(mleEst), mean(mmEst)
    varMLE, varMM = var(mleEst), var(mmEst)

    res[:MSeMLE][i] = varMLE + (meanMLE - trueB)^2
    res[:MSeMM][i] = varMM + (meanMM - trueB)^2
    res[:VarMLE][i] = varMLE
    res[:VarMM][i] = varMM
    res[:BiasMLE][i] = meanMLE - trueB
    res[:BiasMM][i] = meanMM - trueB
end

# ╔═╡ 57b9288e-3af7-11eb-2faf-3949ce5d7593
begin
	p1 = scatter(sampleSizes, [res[:MSeMLE] res[:MSeMM]], c=[:blue :red],
		label=["Mean sq.err (MLE)" "Mean sq.err (MM)"])
	p2 = scatter(sampleSizes, [res[:VarMLE] res[:VarMM]], c=[:blue :red],
		label=["Variance (MLE)" "Variance (MM)"])
	p3 = scatter(sampleSizes, [res[:BiasMLE] res[:BiasMM]], c=[:blue :red],
		label=["Bias (MLE)" "Bias (MM)"])

	plot(p1, p2, p3, ms=10, shape=:xcross, xlabel="n",
		layout=(1,3), size=(1200, 400))
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing5.11"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═5790308c-3af7-11eb-1ebf-b57c8c5e4061
# ╠═57b8a47c-3af7-11eb-3da1-aba055f63a52
# ╠═57b9288e-3af7-11eb-2faf-3949ce5d7593
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
module RestrictedDiscreteModelsTests

using Gridap.Arrays
using Gridap.ReferenceFEs
using Gridap.Geometry
using Test

domain = (0,2,0,2,0,2)
partition = (10,10,10)
oldmodel = CartesianDiscreteModel(domain,partition)

const R = 0.7

function is_in(coords)
  n = length(coords)
  x = (1/n)*sum(coords)
  d = x[1]^2 + x[2]^2 - R^2
  d < 0
end

oldgrid = get_grid(oldmodel)
oldcell_to_coods = get_cell_coordinates(oldgrid)
cell_to_mask = collect1d(apply(is_in,oldcell_to_coods))
cell_to_oldcell = findall(cell_to_mask)

labels = get_face_labeling(oldmodel)

ne = num_entities(labels)
fluid_entity = ne+1
solid_entity = ne+2
cell_to_entity = get_cell_entity(labels)
cell_to_entity .= fluid_entity
cell_to_entity[cell_to_oldcell] .= solid_entity
add_tag!(labels,"fluid",[fluid_entity])
add_tag!(labels,"solid",[solid_entity])

model = RestrictedDiscreteModel(oldmodel,cell_to_oldcell)
test_discrete_model(model)

trian = get_triangulation(model)
@test isa(trian,RestrictedTriangulation)

trian = Triangulation(model)
@test isa(trian,RestrictedTriangulation)

model = DiscreteModel(oldmodel,cell_to_oldcell)
test_discrete_model(model)
@test isa(model,RestrictedDiscreteModel)

model = DiscreteModel(oldmodel,cell_to_mask)
test_discrete_model(model)
@test isa(model,RestrictedDiscreteModel)

model_fluid = DiscreteModel(oldmodel,"fluid")
test_discrete_model(model_fluid)
@test isa(model_fluid,RestrictedDiscreteModel)

model_solid = DiscreteModel(oldmodel,"solid")
test_discrete_model(model_solid)
@test isa(model_solid,RestrictedDiscreteModel)

itrian = InterfaceTriangulation(model_fluid,model_solid)

#using Gridap.Visualization
#writevtk(itrian,"itrian",nsubcells=3,cellfields=["normal"=>get_normal_vector(itrian)])

#
#writevtk(model,"model")
#writevtk(trian,"trian")

oldtrian = Triangulation(oldmodel)
trian_s = SkeletonTriangulation(model)
trian_b = BoundaryTriangulation(model)

meas_K_b = cell_measure(trian_b,oldtrian)
meas_K_sl = cell_measure(trian_s.left,oldtrian)
meas_K_sr = cell_measure(trian_s.right,oldtrian)

oldcell_to_cell = zeros(Int,num_cells(oldmodel))
oldcell_to_cell[cell_to_oldcell] .= 1:length(cell_to_oldcell)

@test all( oldcell_to_cell[findall(meas_K_b .!= 0)] .!= 0 )
@test all( oldcell_to_cell[findall(meas_K_sl .!= 0)] .!= 0 )
@test all( oldcell_to_cell[findall(meas_K_sr .!= 0)] .!= 0 )

#using Gridap.Visualization
#writevtk(trian_b,"trian_b")
#writevtk(trian_s,"trian_s")
#writevtk(trian,"trian")
#writevtk(oldtrian,"oldtrian",celldata=["meas_K_b"=>meas_K_b,"meas_K_sl"=>meas_K_sl,"meas_K_sr"=>meas_K_sr])

end # module


text:
"""
    ScalarizeBridge{T, F, S}

Transforms a constraint `AbstractVectorFunction`-in-`vector_set_type(S)` where
`S <: LPCone{T}` to `F`-in-`S`.
"""
mutable struct ScalarizeBridge{T,F,S} <: AbstractBridge
    scalar_constraints::Vector{CI{F,S}}
    constants::Vector{T}
end
function bridge_constraint(
    ::Type{ScalarizeBridge{T,F,S}},
    model::MOI.ModelLike,
    f::MOI.AbstractVectorFunction,
    set::MOIU.VectorLinearSet,
) where {T,F,S}
    dimension = MOI.output_dimension(f)
    constants = MOI.constant(f, T)
    new_f = MOIU.scalarize(f, true)
    constraints = Vector{CI{F,S}}(undef, dimension)
    for i in 1:dimension
        constraints[i] = MOI.add_constraint(model, new_f[i], S(-constants[i]))
    end
    return ScalarizeBridge{T,F,S}(constraints, constants)
end

function MOI.supports_constraint(
    ::Type{ScalarizeBridge{T}},
    F::Type{<:MOI.AbstractVectorFunction},
    ::Type{<:MOIU.VectorLinearSet},
) where {T}
    # If `F` is `MOI.VectorAffineFunction{Complex{Float64}}`, `S` is `MOI.Zeros` and `T` is `Float64`,
    # it would create a set `MOI.EqualTo{Float64}` which is incorrect hence we say
    # we only support it if the coefficient type of `F` is `T`.
    return MOIU.is_coefficient_type(F, T)
end
function MOIB.added_constrained_variable_types(::Type{<:ScalarizeBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(
    ::Type{ScalarizeBridge{T,F,S}},
) where {T,F,S}
    return [(F, S)]
end
function concrete_bridge_type(
    ::Type{<:ScalarizeBridge{T}},
    F::Type{<:MOI.AbstractVectorFunction},
    S::Type{<:MOIU.VectorLinearSet},
) where {T}
    return ScalarizeBridge{T,MOIU.scalar_type(F),MOIU.scalar_set_type(S, T)}
end

# Attributes, Bridge acting as a model
function MOI.get(
    bridge::ScalarizeBridge{T,F,S},
    ::MOI.NumberOfConstraints{F,S},
) where {T,F,S}
    return length(bridge.scalar_constraints)
end
function MOI.get(
    bridge::ScalarizeBridge{T,F,S},
    ::MOI.ListOfConstraintIndices{F,S},
) where {T,F,S}
    return bridge.scalar_constraints
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::ScalarizeBridge)
    for ci in bridge.scalar_constraints
        MOI.delete(model, ci)
    end
end
function MOI.delete(
    model::MOI.ModelLike,
    bridge::ScalarizeBridge,
    i::IndexInVector,
)
    MOI.delete(model, bridge.scalar_constraints[i.value])
    deleteat!(bridge.scalar_constraints, i.value)
    deleteat!(bridge.constants, i.value)
    return
end

# Attributes, Bridge acting as a constraint
function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintFunction,
    bridge::ScalarizeBridge{T},
) where {T}
    func = MOIU.vectorize(MOI.get.(model, attr, bridge.scalar_constraints))
    if !(func isa MOI.VectorOfVariables)
        # `func` is in terms of bridged variables here while in
        # `bridge_constraint` it was in terms of the solver variables so
        # `MOI.constant(set)` might be different than `bridge.constants[i]`.
        for i in eachindex(bridge.scalar_constraints)
            set = MOI.get(
                model,
                MOI.ConstraintSet(),
                bridge.scalar_constraints[i],
            )
            func = MOIU.operate_output_index!(-, T, i, func, MOI.constant(set))
        end
    end
    return func
end
function MOI.get(
    ::MOI.ModelLike,
    ::MOI.ConstraintSet,
    bridge::ScalarizeBridge{T,F,S},
) where {T,F,S}
    return MOIU.vector_set_type(S)(length(bridge.scalar_constraints))
end

function MOI.supports(
    ::MOI.ModelLike,
    ::Union{MOI.ConstraintPrimalStart,MOI.ConstraintDualStart},
    ::Type{<:ScalarizeBridge},
)
    return true
end
function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimal,MOI.ConstraintPrimalStart},
    bridge::ScalarizeBridge,
)
    return MOI.get.(model, attr, bridge.scalar_constraints) .+ bridge.constants
end
function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimalStart,
    bridge::ScalarizeBridge,
    value,
)
    # TODO do no add constant if the primal status is a ray like in Vectorize
    return MOI.set.(
        model,
        attr,
        bridge.scalar_constraints,
        value .- bridge.constants,
    )
end
function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintDual,MOI.ConstraintDualStart},
    bridge::ScalarizeBridge,
)
    return MOI.get.(model, attr, bridge.scalar_constraints)
end
function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintDualStart,
    bridge::ScalarizeBridge,
    value,
)
    return MOI.set.(model, attr, bridge.scalar_constraints, value)
end
function MOI.modify(
    model::MOI.ModelLike,
    bridge::ScalarizeBridge{T,F,S},
    change::MOI.VectorConstantChange{T},
) where {T,F,S}
    bridge.constants = change.new_constant
    return MOI.set.(
        model,
        MOI.ConstraintSet(),
        bridge.scalar_constraints,
        S.(-change.new_constant),
    )
end
function MOI.modify(
    model::MOI.ModelLike,
    bridge::ScalarizeBridge,
    change::MOI.MultirowChange{T},
) where {T}
    for (index, value) in change.new_coefficients
        MOI.modify(
            model,
            bridge.scalar_constraints[index],
            MOI.ScalarCoefficientChange{T}(change.variable, value),
        )
    end
    return nothing
end
function MOI.set(
    model::MOI.ModelLike,
    ::MOI.ConstraintFunction,
    bridge::ScalarizeBridge{T,F,S},
    func,
) where {T,F,S}
    old_constants = bridge.constants
    bridge.constants = MOI.constant(func, T)
    new_func = MOIU.scalarize(func, true)
    MOI.set.(
        model,
        MOI.ConstraintFunction(),
        bridge.scalar_constraints,
        new_func,
    )
    for i in eachindex(bridge.constants)
        if bridge.constants[i] != old_constants[i]
            MOI.set(
                model,
                MOI.ConstraintSet(),
                bridge.scalar_constraints[i],
                S(-bridge.constants[i]),
            )
        end
    end
end

# TODO implement transform


text:
using Statistics, LinearAlgebra

"""
    sqrt_svd(A)

Returns the square root matrix by SVD
"""
function sqrt_svd(A::AbstractMatrix)
    F = svd(A)
    F.U * diagm(sqrt.(F.S)) * F.Vt
end

export RMSE

"""
    RMSE(e)

Returns the Root Mean Squared Error
"""
RMSE(e) = sqrt.(mean(e .^ 2))

"""
    normalise!( w )

Normalize the entries of a multidimensional array sum to 1.
"""
function normalise!(w)

    c = sum(w)
    # Set any zeros to one before dividing
    c += c == 0
    w ./= c

end

"""
    mk_stochastic!(w)

Ensure the matrix is stochastic, i.e.,
the sum over the last dimension is 1.
"""
function mk_stochastic!(w::Array{Float64,2})

    if first(size(w)) == 1
        normalise!(w)
    else
        # Copy the normaliser plane for each i.
        normaliser = sum(w, dims = 2)
        # Set zeros to 1 before dividing
        # This is valid since normaliser(i) = 0 iff T(i) = 0
        normaliser .+= normaliser .== 0
        w ./= normaliser
    end

end

"""
    sample_discrete(prob, r, c)

Sampling from a non-uniform distribution.
"""
function sample_discrete(prob, r, c)

    # this speedup is due to Peter Acklam
    cumprob = cumsum(prob)
    R = rand(r, c)
    M = zeros(Int64, (r, c))
    N = length(cumprob)
    for i = 1:N-1
        M .+= R .> cumprob[i]
    end
    M

end

"""
    inv_using_SVD(Mat, eigvalMax)

SVD decomposition of Matrix.
"""
function inv_using_SVD(Mat, eigvalMax)

    F = svd(Mat; full = true)
    eigval = cumsum(F.S) ./ sum(F.S)
    # search the optimal number of eigen values
    icut = findfirst(eigval .>= eigvalMax)

    U_1 = @view F.U[1:icut, 1:icut]
    V_1 = @view F.Vt'[1:icut, 1:icut]
    tmp1 = (V_1 ./ F.S[1:icut]') * U_1'

    if icut + 1 > length(eigval)
        tmp1
    else
        U_3 = @view F.U[icut+1:end, 1:icut]
        V_3 = @view F.Vt'[icut+1:end, 1:icut]
        tmp2 = (V_1 ./ F.S[1:icut]') * U_3'
        tmp3 = (V_3 ./ F.S[1:icut]') * U_1'
        tmp4 = (V_3 ./ F.S[1:icut]') * U_3'
        vcat(hcat(tmp1, tmp2), hcat(tmp3, tmp4))
    end

end

"""
    resample_multinomial( w )

Multinomial resampler.
"""
function resample_multinomial(w::Vector{Float64})

    m = length(w)
    q = cumsum(w)
    q[end] = 1.0 # Just in case...
    i = 1
    indx = Int64[]
    while i <= m
        sampl = rand()
        j = 1
        while q[j] < sampl
            j = j + 1
        end
        push!(indx, j)
        i = i + 1
    end
    indx
end

"""
    resample!( indx, w )

Multinomial resampler.
"""
function resample!(indx::Vector{Int64}, w::Vector{Float64})

    m = length(w)
    q = cumsum(w)
    i = 1
    while i <= m
        sampl = rand()
        j = 1
        while q[j] < sampl
            j = j + 1
        end
        indx[i] = j
        i = i + 1
    end
end

"""
    ensure_pos_sym(M; ϵ= 1e-8)

Ensure that matrix `M` is positive and symmetric to avoid numerical errors when numbers are small by doing `(M + M')/2 + ϵ*I`

reference : [StateSpaceModels.jl](https://github.com/LAMPSPUC/StateSpaceModels.jl)
"""
function ensure_pos_sym(M::Matrix{T}; ϵ::T = 1e-8) where {T<:AbstractFloat}
    return (M + M') ./ 2 + ϵ * I
end

function ensure_pos_sym!(M::Matrix{T}; ϵ::T = 1e-8) where {T<:AbstractFloat}
    M .= (M + M') ./ 2 + ϵ * I
end


text:
using Plots
using Test: @test_throws

export rect_sino

"""
`(sino, pos, ang) = rect_sino(sg, rects;
            oversample=1, xscale=1, yscale=1)`

Create sinogram projection of one or more rectangles.
Works for both parallel-beam geometry and for fan-beam geometry.

in
* `sg`                  sinogram geometry object from sino_geom()
* `rects`               [ne 6]
                        [centx centy widthx width y angle_degrees value]

options
* `oversample::Integer` oversampling factor for emulating "strips"
                        default 1: just 1 ray per detector element
* `xscale::Integer`     use -1 to flip in x (not recommended)
* `yscale::Integer`     use -1 to flip in y (not recommended)
* `how::Symbol`        `:fan` `:par` `:moj`

out
* `sino`                [nb na] sinogram
* `pos`                 [nb 1]  radial sample positions (or [nb na] is mojette)
* `ang`                 [na]    angular sample locations (in radians)
"""
function rect_sino(sg::MIRT_sino_geom,
    rects::AbstractArray{<:Real,2};
    oversample::Integer=1,
    xscale::Integer=1,
    yscale::Integer=1)
    how = sg.how
	size(rects, 2) != 6 && throw("6 parameters per rect")

    if how == :fan
        (sino, pos) = rect_sino_go(rects, sg.nb, sg.ds, sg.offset, sg.na,
                        sg.ar, sg.dso, sg.dod, sg.dfs, sg.source_offset,
                        xscale, yscale, oversample, 0)
    elseif how == :par
        (sino, pos) = rect_sino_go(rects, sg.nb, sg.dr, sg.offset, sg.na,
                        sg.ar, Inf, 1, 0, sg.source_offset, xscale, yscale, oversample, 0)
    elseif how == :moj
        (sino, pos) = rect_sino_go(rects, sg.nb, [], sg.offset, sg.na,
                        sg.ar, Inf, 1, 0, sg.source_offset, xscale, yscale, oversample, sg.d)
    else
        throw("sino geom $how not done")
    end

    return (sino, pos)
end

"""
`rect_sino_go()`
"""
function rect_sino_go(rects, nb, ds, offset_s, na, ang, dso, dod, dfs,
    source_offset, xscale, yscale, oversample, mojette)

    (pos, pos2) = rect_sino_pos([], nb, ds, offset_s, oversample, mojette, ang)

    sino = rect_sino_do(rects, pos2, ang[:]', xscale, yscale, dso, dod, dfs, source_offset)

    if oversample != 1
        sino = downsample2(sino, [oversample 1])
    end
    return (sino, pos)
end

"""
`rect_sino_pos()`

determine usual and fine "radial" sample positions
"""
function rect_sino_pos(pos, nb, ds, offset_s, nover, mojette, ang)
    wb = (nb - 1)/2 + offset_s
    if mojette != 0 # tricky mojette radial sampling
        # trick: ray_spacing aka ds comes from dx which is in mojette
        if !isempty(ds)
            throw("ds must be empty for mojette case")
        end

        dt = abs(mojette) * max.(abs.(cos.(ang)), abs.(sin.(ang)))' # [1, na]

        if !isempty(pos)
            throw("mojette requires empty 'pos'")
        end
        na = length(ang)
        pos_coarse = ((0:nb - 1) .- wb) * dt[:]' # [nb na]

        if nover > 1
            pos_fine = zeros(nover*nb, na)
            for ia in 1:na
                tmp = (-(over - 1):2:(nover - 1)) / (2*nover) * dt[ia]
                pos_fine[:, ia] = (tmp .+ pos_coarse[:,ia]')[:]
            end
        else
            pos_fine = pos_coarse
        end
    else # ordinary mojette sampling
        if isempty(pos)
            pos_coarse = ((0:(nb - 1)) .- wb) .* ds' # [nb]
        else
            pos_coarse = pos[:] # [nb 1]
            ds = pos[2] - pos[1]
            if any(abs(diff(pos) / ds - 1) > 1e-10)
				throw("uniform spacing required")
			end
        end
        if nover > 1
            # determine fine sampling positions
            pos_fine = (-(nover - 1):2:(nover-1)) / (2*nover) * ds
            pos_fine = pos_fine .+ pos_coarse' # [nover nb]
            pos_fine = pos_fine[:]
        else
            pos_fine = pos_coarse
        end
    end
    return (pos_coarse, pos_fine)
end

"""
`rect_sino_do()`

analytical line-integral projections rectangle
"""
function rect_sino_do(rects, pos, ang, xscale, yscale, dso, dod, dfs, source_offset)
	nb = size(pos, 1)
	na = maximum(size(ang))

	# effective radial and angular sample locations in parallel geometry
	if isinf(dso)
		if size(pos, 2) > 1 # mojette
			rads = pos
			angs = repeat(ang, nb, 1) # [nb na]
		else

			(rads, angs) = ndgrid(pos, ang') # [nb na]
		end
	else # fan
		if size(pos, 2) > 1 # mojette
			throw("mojette fan not supported")
		end
		dis_src_det = dso + dod

		if isinf(dfs) # flat detector
			gam = atan(pos / dis_src_det) # gamma
		else
			dis_foc_det = dfs + dis_src_det
			alf = pos / dis_foc_det
			gam = atan.(dis_foc_det * sin.(alf), dis_foc_det * cos.(alf) .- dfs) #gamma
		end

		rad = dso * sin.(gam) + source_offset * cos.(gam)
		rads = repeat(rad, 1, na) # [nb, na]

		angs = gam .+ ang  # [nb na] gamma + beta
	end

	sino = zeros(nb, na)

	cangs = cos.(angs)
	sangs = sin.(angs)

	#loop over rects
	#ticker reset
	ne = size(rects, 1)
	for ie in 1:ne
		#ticker(mfilename, ie, ne)
		rect = rects[ie, :] # [1, 6]

		cx = xscale * rect[1]
		cy = yscale * rect[2]
		wx = rect[3]
		wy = rect[4]
		if (wx <= 0) || (wy <= 0)
			throw("need positive rectangle sizes")
		end
		eang = rect[5] * (pi/180)
		if yscale == -1
			eang = -eang
		end
		if xscale == -1
			eang = pi - eang
		end
		val = rect[6]

		cos_rot = cangs .* cos(eang) + sangs .* sin(eang)
		sin_rot = sangs * cos(eang) - cangs * sin(eang)
		rp = sqrt.((wx * cos_rot).^2 + (wy * sin_rot).^2) # projected radius

		sp = cx .* cangs + cy .* sangs # radial shift
		dis = (rads - sp) ./ rp # scaled distance from center

		# projection angle after affine scaling and rotate
		abs_cos_ang_pi = wx * abs.(cos_rot) ./ rp
		abs_sin_ang_pi = wy * abs.(sin_rot) ./ rp

		# break points of the trapezoid
		len = 1 ./ max.(abs_cos_ang_pi, abs_sin_ang_pi)
		dmax = (abs_cos_ang_pi + abs_sin_ang_pi) / 2
		dbreak = abs.(abs_cos_ang_pi - abs_sin_ang_pi) / 2
		dmax_break = dmax - dbreak
		scale = val * wx * wy ./ rp .* len #lmax

		sino += scale .* trapezoid.(dis, -dmax, -dbreak, dbreak, dmax)
	end
	return sino
end

"""
'div0'
"""
function div0(x, y)
	if y == 0
		return 0
	else
		return x/y
	end
end

"""
"""
function trapezoid(t::Real, t1, t2, t3, t4)
	if t1 < t < t2
		return (t - t1)/(t2 - t1)
	elseif t2 <= t <= t3
		return 1
	elseif t3 < t < t4
		return (t4 - t)/(t4 - t3)
	else
		return 0
	end
end

"""
`rect_sino()`

shows doc strings
"""
function rect_sino()
	@doc rect_sino
end

"""
`rect_sino_test()`

internal test routine: standard sampling
"""
function rect_sino_test()
	down = 1
	ig = image_geom(nx=512, ny=504, fov=500)
	ig = ig.down(down)
	rect = [[60 10.5 38 27] .* ig.dx 0 1]
	# (xtrue, rect) = rect_im(ig, rect; oversample=3)

	sgf = sino_geom(:fan, nb = 888, na = 984, d = 1.0, orbit = 360,
			orbit_start = 0, offset = 0.75, dsd = 949, dod = 408, down=down) # fan
	sgp = sino_geom(:par, nb = 888, na = 984, down=down, d = 0.5, orbit = sgf.orbit,
			orbit_start = sgf.orbit_start, offset = 0.25) # parallel
	sgm = sino_geom(:moj, nb = 888, na = 984, down=down, d = 0.5, orbit = sgf.orbit,
			orbit_start = sgf.orbit_start, offset = 0.25) # mojette

	sino = sino_geom(:bad)
	@test_throws String rect_sino(sino, rect)

	sino_f = rect_sino(sgf, rect; oversample=1)
	sino_p = rect_sino(sgp, rect; oversample=1)
	sino_m = rect_sino(sgm, rect; oversample=2)
	r1 = rect_sino(sgf, rect; oversample=1, yscale=-1, xscale=-1)
	r2 = rect_sino(sgp, rect; oversample=1, yscale=-1, xscale=-1)
	r3 = rect_sino(sgm, rect; oversample=1, yscale=-1, xscale=-1)
	t = LinRange(-2, 5, 101)
	f = trapezoid.(t, 1, 2, 3, 4)
	div0(4, 2)
	div0(1, 0)
end

"""
`rect_sino_show`
"""
function rect_sino_show()
	down = 1
	ig = image_geom(nx=512, ny=504, fov=500)
	ig = ig.down(down)
	rect = [[60 10.5 38 27] .* ig.dx 0 1]

	sgf = sino_geom(:fan, nb = 888, na = 984, d = 1.0, orbit = 360,
			orbit_start = 0, offset = 0.75, dsd = 949, dod = 408, down=down) # fan
	sgp = sino_geom(:par, nb = 888, na = 984, down=down, d = 0.5, orbit = sgf.orbit,
			orbit_start = sgf.orbit_start, offset = 0.25) # parallel

	#sg = sino_geom(:ge1, nb=888, na=984, down=down)
	#sino = rect_sino(sg, rect; oversample=10)

	xtrue = rect_im(ig, rect; oversample=3)
	sino_f = rect_sino(sgf, rect; oversample=1)[1]
	sino_p = rect_sino(sgp, rect; oversample=1)[1]

	x1 = jim(sino_f, title="fan")
	x2 = jim(sino_p, title="parallel")
	x3 = jim(xtrue, title="xtrue")
	#x4 = jim(sino, title="ya")

	plot(x1, x2, x3)
end

"""
`rect_sino(:show)`

`rect_sino(:test)`
"""
function rect_sino(test::Symbol)
	if test == :show
		return rect_sino_show()
	end
	rect_sino()
	rect_sino_test()
	rect_sino(:show)
	true
end

#rect_sino(:test)
#rect_sino(:show)


text:
struct LaserPolarization{C}
    ξx::C
    ξy::C
end


text:
using HiRedis
using Base.Test

using Logging

Logging.configure(level=INFO)

start_session("127.0.0.1", 6379)

info("select db 2 for testing, flushdb")
selectdb(2)
flushdb()

info("test incr")
incr("fooincr")
incr("fooincr")
fooincr = incr("fooincr")
@test fooincr == 3

info("test set and get")
kvset("foo", "bar")
@test kvget("foo") == "bar"

info("test hset and hget")
hset("customer:1", "first_name", "Mark")
@test hget("customer:1", "first_name") == "Mark"

info("test hmset and hmget")
hmset("customer:1", "first_name", "Alice", "last_name", "Bean")
reply = hmget("customer:1", "first_name", "last_name")
@test reply[2] == "Bean"

info("test hgetall")
reply = hgetall("customer:1")
@test reply["last_name"] == "Bean"

info("test sadd and smembers")
sadd("account_keyset:1", 1)
sadd("account_keyset:1", 2)
sadd("account_keyset:1", 2)
sadd("account_keyset:1", 3)
@test length(smembers("account_keyset:1")) == 3

info("pipeline commands")
sadd("pipeline_test:1", "foo", pipeline=true)
sadd("pipeline_test:1", "bar", pipeline=true)
sadd("pipeline_test:1", "bar", pipeline=true)
sadd("pipeline_test:1", "baz", pipeline=true)
smembers("pipeline_test:1", pipeline=true)
replies = get_reply()
@test length(replies[5]) == 3

info("test pipelining using the @pipeline macro")
@pipeline begin
    incr("maz")
    incr("maz")
    incr("maz")
    kvget("maz")
end
replies = get_reply()
@test int(replies[4]) == 3

info("flushdb, select db 1")
flushdb()
selectdb(1)

end_session()


text:
mutable struct Financial
    inflation::Float64
    apr::Float64
    term::Int64
    downpayment::Float64
    XhR::Float64
end

max_stars_repo_name: arnavgautam/PVCR

max_issues_repo_name: arnavgautam/PVCR

max_forks_repo_name: arnavgautam/PVCR


text:
module SuperSub

supers = split("⁰¹²³⁴⁵⁶⁷⁸⁹", "")
subs = split("₀₁₂₃₄₅₆₇₈₉", "")
superscript(i::Integer) = join(map(i -> supers[i+1], reverse(digits(i))), "")
superscript(r::Rational) = "$(superscript(num(r))).$(superscript(den(r)))"
subscript(i::Integer) = join(map(i -> subs[i+1], reverse(digits(i))), "")
subscript(r::Rational) = "$(subscript(num(r))).$(subscript(den(r)))"

export superscript, subscript

end # module


text:
###
### **mag.h** -- fixed-precision unsigned floating-point numbers for bounds
###

### Types, macros and constants

### Memory management
arbcall"void mag_init(mag_t x)"
arbcall"void mag_clear(mag_t x)"
arbcall"void mag_swap(mag_t x, mag_t y)"
#ni arbcall"mag_ptr _mag_vec_init(slong n)"
#ni arbcall"void _mag_vec_clear(mag_ptr v, slong n)"
arbcall"slong mag_allocated_bytes(const mag_t x)"

### Special values
arbcall"void mag_zero(mag_t res)"
arbcall"void mag_one(mag_t res)"
arbcall"void mag_inf(mag_t res)"
arbcall"int mag_is_special(const mag_t x)"
arbcall"int mag_is_zero(const mag_t x)"
arbcall"int mag_is_inf(const mag_t x)"
arbcall"int mag_is_finite(const mag_t x)"

### Assignment and conversions
arbcall"void mag_init_set(mag_t res, const mag_t x)"
arbcall"void mag_set(mag_t res, const mag_t x)"
arbcall"void mag_set_d(mag_t res, double x)"
#ns arbcall"void mag_set_fmpr(mag_t res, const fmpr_t x)"
arbcall"void mag_set_ui(mag_t res, ulong x)"
#ni arbcall"void mag_set_fmpz(mag_t res, const fmpz_t x)"
arbcall"void mag_set_d_lower(mag_t res, double x)"
arbcall"void mag_set_ui_lower(mag_t res, ulong x)"
#ni arbcall"void mag_set_fmpz_lower(mag_t res, const fmpz_t x)"
#ni arbcall"void mag_set_d_2exp_fmpz(mag_t res, double x, const fmpz_t y)"
#ni arbcall"void mag_set_fmpz_2exp_fmpz(mag_t res, const fmpz_t x, const fmpz_t y)"
arbcall"void mag_set_ui_2exp_si(mag_t res, ulong x, slong y)"
#ni arbcall"void mag_set_d_2exp_fmpz_lower(mag_t res, double x, const fmpz_t y)"
#ni arbcall"void mag_set_fmpz_2exp_fmpz_lower(mag_t res, const fmpz_t x, const fmpz_t y)"
arbcall"double mag_get_d(const mag_t x)"
arbcall"double mag_get_d_log2_approx(const mag_t x)"
#ns arbcall"void mag_get_fmpr(fmpr_t res, const mag_t x)"
#ni arbcall"void mag_get_fmpq(fmpq_t res, const mag_t x)"
#ni arbcall"void mag_get_fmpz(fmpz_t res, const mag_t x)"
#ni arbcall"void mag_get_fmpz_lower(fmpz_t res, const mag_t x)"

### Comparisons
arbcall"int mag_equal(const mag_t x, const mag_t y)"
arbcall"int mag_cmp(const mag_t x, const mag_t y)"
arbcall"int mag_cmp_2exp_si(const mag_t x, slong y)"
arbcall"void mag_min(mag_t res, const mag_t x, const mag_t y)"
arbcall"void mag_max(mag_t res, const mag_t x, const mag_t y)"

### Input and output
#mo arbcall"void mag_print(const mag_t x)" # clashes with Base.print
#ns arbcall"void mag_fprint(FILE * file, const mag_t x)"
arbcall"char * mag_dump_str(const mag_t x)"
arbcall"int mag_load_str(mag_t x, const char * str)"
#ns arbcall"int mag_dump_file(FILE * stream, const mag_t x)"
#ns arbcall"int mag_load_file(mag_t x, FILE * stream)"

### Random generation
#ns arbcall"void mag_randtest(mag_t res, flint_rand_t state, slong expbits)"
#ns arbcall"void mag_randtest_special(mag_t res, flint_rand_t state, slong expbits)"

### Arithmetic
arbcall"void mag_add(mag_t res, const mag_t x, const mag_t y)"
arbcall"void mag_add_ui(mag_t res, const mag_t x, ulong y)"
arbcall"void mag_add_lower(mag_t res, const mag_t x, const mag_t y)"
arbcall"void mag_add_ui_lower(mag_t res, const mag_t x, ulong y)"
#ni arbcall"void mag_add_2exp_fmpz(mag_t res, const mag_t x, const fmpz_t e)"
arbcall"void mag_add_ui_2exp_si(mag_t res, const mag_t x, ulong y, slong e)"
arbcall"void mag_sub(mag_t res, const mag_t x, const mag_t y)"
arbcall"void mag_sub_lower(mag_t res, const mag_t x, const mag_t y)"
arbcall"void mag_mul_2exp_si(mag_t res, const mag_t x, slong y)"
#ni arbcall"void mag_mul_2exp_fmpz(mag_t res, const mag_t x, const fmpz_t y)"
arbcall"void mag_mul(mag_t res, const mag_t x, const mag_t y)"
arbcall"void mag_mul_ui(mag_t res, const mag_t x, ulong y)"
#ni arbcall"void mag_mul_fmpz(mag_t res, const mag_t x, const fmpz_t y)"
arbcall"void mag_mul_lower(mag_t res, const mag_t x, const mag_t y)"
arbcall"void mag_mul_ui_lower(mag_t res, const mag_t x, ulong y)"
#ni arbcall"void mag_mul_fmpz_lower(mag_t res, const mag_t x, const fmpz_t y)"
arbcall"void mag_addmul(mag_t z, const mag_t x, const mag_t y)"
arbcall"void mag_div(mag_t res, const mag_t x, const mag_t y)"
arbcall"void mag_div_ui(mag_t res, const mag_t x, ulong y)"
#ni arbcall"void mag_div_fmpz(mag_t res, const mag_t x, const fmpz_t y)"
arbcall"void mag_div_lower(mag_t res, const mag_t x, const mag_t y)"
arbcall"void mag_inv(mag_t res, const mag_t x)"
arbcall"void mag_inv_lower(mag_t res, const mag_t x)"

### Fast, unsafe arithmetic
arbcall"void mag_fast_init_set(mag_t x, const mag_t y)"
arbcall"void mag_fast_zero(mag_t res)"
arbcall"int mag_fast_is_zero(const mag_t x)"
arbcall"void mag_fast_mul(mag_t res, const mag_t x, const mag_t y)"
arbcall"void mag_fast_addmul(mag_t z, const mag_t x, const mag_t y)"
arbcall"void mag_fast_add_2exp_si(mag_t res, const mag_t x, slong e)"
arbcall"void mag_fast_mul_2exp_si(mag_t res, const mag_t x, slong e)"

### Powers and logarithms
arbcall"void mag_pow_ui(mag_t res, const mag_t x, ulong e)"
#ni arbcall"void mag_pow_fmpz(mag_t res, const mag_t x, const fmpz_t e)"
arbcall"void mag_pow_ui_lower(mag_t res, const mag_t x, ulong e)"
#ni arbcall"void mag_pow_fmpz_lower(mag_t res, const mag_t x, const fmpz_t e)"
arbcall"void mag_sqrt(mag_t res, const mag_t x)"
arbcall"void mag_sqrt_lower(mag_t res, const mag_t x)"
arbcall"void mag_rsqrt(mag_t res, const mag_t x)"
arbcall"void mag_rsqrt_lower(mag_t res, const mag_t x)"
arbcall"void mag_hypot(mag_t res, const mag_t x, const mag_t y)"
arbcall"void mag_root(mag_t res, const mag_t x, ulong n)"
arbcall"void mag_log(mag_t res, const mag_t x)"
arbcall"void mag_log_lower(mag_t res, const mag_t x)"
arbcall"void mag_neg_log(mag_t res, const mag_t x)"
arbcall"void mag_neg_log_lower(mag_t res, const mag_t x)"
arbcall"void mag_log_ui(mag_t res, ulong n)"
arbcall"void mag_log1p(mag_t res, const mag_t x)"
arbcall"void mag_exp(mag_t res, const mag_t x)"
arbcall"void mag_exp_lower(mag_t res, const mag_t x)"
arbcall"void mag_expinv(mag_t res, const mag_t x)"
arbcall"void mag_expinv_lower(mag_t res, const mag_t x)"
arbcall"void mag_expm1(mag_t res, const mag_t x)"
arbcall"void mag_exp_tail(mag_t res, const mag_t x, ulong N)"
arbcall"void mag_binpow_uiui(mag_t res, ulong m, ulong n)"
arbcall"void mag_geom_series(mag_t res, const mag_t x, ulong N)"

### Special functions
arbcall"void mag_const_pi(mag_t res)"
arbcall"void mag_const_pi_lower(mag_t res)"
arbcall"void mag_atan(mag_t res, const mag_t x)"
arbcall"void mag_atan_lower(mag_t res, const mag_t x)"
arbcall"void mag_cosh(mag_t res, const mag_t x)"
arbcall"void mag_cosh_lower(mag_t res, const mag_t x)"
arbcall"void mag_sinh(mag_t res, const mag_t x)"
arbcall"void mag_sinh_lower(mag_t res, const mag_t x)"
arbcall"void mag_fac_ui(mag_t res, ulong n)"
arbcall"void mag_rfac_ui(mag_t res, ulong n)"
arbcall"void mag_bin_uiui(mag_t res, ulong n, ulong k)"
arbcall"void mag_bernoulli_div_fac_ui(mag_t res, ulong n)"
arbcall"void mag_polylog_tail(mag_t res, const mag_t z, slong s, ulong d, ulong N)"
arbcall"void mag_hurwitz_zeta_uiui(mag_t res, ulong s, ulong a)"


text:
@testset "warnings" begin
    @testset "exclusivity warning" begin
        withenv("JULIA_DEBUG" => "all", "SUPPRESS_THREADINGUTILITIES_WARNING" => "0", "JULIA_EXCLUSIVE" => "") do
            @test (@test_logs (:debug, "The JULIA_EXCLUSIVE environment variable is not set to 1. Therefore, the kernel is allowed to move threads to different cores. We recommend that you set the JULIA_EXCLUSIVE environment variable to 1. To suppress this warning, set the SUPPRESS_THREADINGUTILITIES_WARNING environment variable to 1.") match_mode=:any ThreadingUtilities._print_exclusivity_warning()) isa Nothing
        end
        withenv("SUPPRESS_THREADINGUTILITIES_WARNING" => "1", "JULIA_EXCLUSIVE" => "") do
            @test ThreadingUtilities._print_exclusivity_warning() isa Nothing
        end
        withenv("JULIA_DEBUG" => "all", "SUPPRESS_THREADINGUTILITIES_WARNING" => "0", "JULIA_EXCLUSIVE" => "0") do
            @test (@test_logs (:debug, "The JULIA_EXCLUSIVE environment variable is not set to 1. Therefore, the kernel is allowed to move threads to different cores. We recommend that you set the JULIA_EXCLUSIVE environment variable to 1. To suppress this warning, set the SUPPRESS_THREADINGUTILITIES_WARNING environment variable to 1.") match_mode=:any ThreadingUtilities._print_exclusivity_warning()) isa Nothing
        end
        withenv("SUPPRESS_THREADINGUTILITIES_WARNING" => "1", "JULIA_EXCLUSIVE" => "0") do
            @test ThreadingUtilities._print_exclusivity_warning() isa Nothing
        end
    end

    @testset "utility functions" begin
        @test ThreadingUtilities._string_to_bool("true")
        @test ThreadingUtilities._string_to_bool("1")
        @test !ThreadingUtilities._string_to_bool("false")
        @test !ThreadingUtilities._string_to_bool("0")
        @test !ThreadingUtilities._string_to_bool("")
    end
end


text:
import Juniper
using Nonconvex, LinearAlgebra, Test

f(x::AbstractVector) = sqrt(x[2])
g(x::AbstractVector, a, b) = (a*x[1] + b)^3 - x[2]

@testset "First order - $first_order" for first_order in [true, false]
    options = JuniperIpoptOptions(first_order = first_order)
    @testset "Simple constraints" begin
        m = Model(f)
        addvar!(m, [0.0, 0.0], [10.0, 10.0])
        add_ineq_constraint!(m, x -> g(x, 2, 0))
        add_ineq_constraint!(m, x -> g(x, -1, 1))

        alg = JuniperIpoptAlg()
        r1 = Nonconvex.optimize(m, alg, [1.234, 2.345], options = options)
        global s1 = sum(r1.minimizer)
        @test abs(r1.minimum - sqrt(8/27)) < 1e-6
        @test norm(r1.minimizer - [1/3, 8/27]) < 1e-6

        setinteger!(m, 1, true)
        r2 = Nonconvex.optimize(m, alg, [1.234, 2.345], options = options)
        global s2 = sum(r2.minimizer)
        @test r2.minimizer[1] - round(Int, r2.minimizer[1]) ≈ 0 atol = 1e-7

        setinteger!(m, 1, false)
        setinteger!(m, 2, true)
        r3 = Nonconvex.optimize(m, alg, [1.234, 2.345], options = options)
        global s3 = sum(r3.minimizer)
        @test r3.minimizer[2] - round(Int, r3.minimizer[2]) ≈ 0 atol = 1e-7
    end

    @testset "Equality constraints" begin
        m = Model(f)
        addvar!(m, [0.0, 0.0], [10.0, 10.0])
        add_ineq_constraint!(m, x -> g(x, 2, 0))
        add_ineq_constraint!(m, x -> g(x, -1, 1))
        add_eq_constraint!(m, x -> sum(x) - s1)
        alg = JuniperIpoptAlg()
        r1 = Nonconvex.optimize(m, alg, [1.234, 2.345], options = options)
        @test abs(r1.minimum - sqrt(8/27)) < 1e-6
        @test norm(r1.minimizer - [1/3, 8/27]) < 1e-6

        m = Model(f)
        addvar!(m, [0.0, 0.0], [10.0, 10.0])
        add_ineq_constraint!(m, x -> g(x, 2, 0))
        add_ineq_constraint!(m, x -> g(x, -1, 1))
        add_eq_constraint!(m, x -> sum(x) - s2)
        setinteger!(m, 1, true)
        r2 = Nonconvex.optimize(m, alg, [1.234, 2.345], options = options)
        @test r2.minimizer[1] - round(Int, r2.minimizer[1]) ≈ 0 atol = 1e-7

        m = Model(f)
        addvar!(m, [0.0, 0.0], [10.0, 10.0])
        add_ineq_constraint!(m, x -> g(x, 2, 0))
        add_ineq_constraint!(m, x -> g(x, -1, 1))
        add_eq_constraint!(m, x -> sum(x) - s3)
        setinteger!(m, 2, true)
        r3 = Nonconvex.optimize(m, alg, [1.234, 2.345], options = options)
        @test r3.minimizer[2] - round(Int, r3.minimizer[2]) ≈ 0 atol = 1e-7
    end

    @testset "Block constraints" begin
        m = Model(f)
        addvar!(m, [0.0, 0.0], [10.0, 10.0])
        add_ineq_constraint!(m, FunctionWrapper(x -> [g(x, 2, 0), g(x, -1, 1)], 2))

        alg = JuniperIpoptAlg()
        r1 = Nonconvex.optimize(m, alg, [1.234, 2.345], options = options)
        @test abs(r1.minimum - sqrt(8/27)) < 1e-6
        @test norm(r1.minimizer - [1/3, 8/27]) < 1e-6

        setinteger!(m, 1, true)
        r2 = Nonconvex.optimize(m, alg, [1.234, 2.345], options = options)
        @test r2.minimizer[1] - round(Int, r2.minimizer[1]) ≈ 0 atol = 1e-7

        setinteger!(m, 1, false)
        setinteger!(m, 2, true)
        r3 = Nonconvex.optimize(m, alg, [1.234, 2.345], options = options)
        @test r3.minimizer[2] - round(Int, r3.minimizer[2]) ≈ 0 atol = 1e-7
    end

    @testset "Infinite bounds" begin
        @testset "Infinite upper bound" begin
            m = Model(f)
            addvar!(m, [0.0, 0.0], [Inf, Inf])
            add_ineq_constraint!(m, x -> g(x, 2, 0))
            add_ineq_constraint!(m, x -> g(x, -1, 1))

            alg = JuniperIpoptAlg()
            r1 = Nonconvex.optimize(m, alg, [1.234, 2.345], options = options)
            @test abs(r1.minimum - sqrt(8/27)) < 1e-6
            @test norm(r1.minimizer - [1/3, 8/27]) < 1e-6

            setinteger!(m, 1, true)
            r2 = Nonconvex.optimize(m, alg, [1.234, 2.345], options = options)
            @test r2.minimizer[1] - round(Int, r2.minimizer[1]) ≈ 0 atol = 1e-7

            setinteger!(m, 1, false)
            setinteger!(m, 2, true)
            r3 = Nonconvex.optimize(m, alg, [1.234, 2.345], options = options)
            @test r3.minimizer[2] - round(Int, r3.minimizer[2]) ≈ 0 atol = 1e-7
        end
        @testset "Infinite lower bound" begin
            m = Model(f)
            addvar!(m, [-Inf, -Inf], [10, 10])
            add_ineq_constraint!(m, x -> g(x, 2, 0))
            add_ineq_constraint!(m, x -> g(x, -1, 1))

            alg = JuniperIpoptAlg()
            r = Nonconvex.optimize(m, alg, [1.234, 2.345], options = options)
            @test abs(r.minimum - sqrt(8/27)) < 1e-6
            @test norm(r.minimizer - [1/3, 8/27]) < 1e-6
        end
        @testset "Infinite upper and lower bound" begin
            m = Model(f)
            addvar!(m, [-Inf, -Inf], [Inf, Inf])
            add_ineq_constraint!(m, x -> g(x, 2, 0))
            add_ineq_constraint!(m, x -> g(x, -1, 1))

            alg = JuniperIpoptAlg()
            r = Nonconvex.optimize(m, alg, [1.234, 2.345], options = options)
            @test abs(r.minimum - sqrt(8/27)) < 1e-6
            @test norm(r.minimizer - [1/3, 8/27]) < 1e-6
        end
    end
end


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

using CompilerTools
using Base.Test

## Tests for CompilerTools.LivenessAnalysis

function test_liveness_2(x::Int64, y::Int64, z::Int64)
	 A = rand(x, y)
	 B = rand(y, z)
	 C = zeros(x, z)
	 for i = 1:x
	     for j = 1:y
	     	 A[i,j] += 1;
	     end
	     for j = 1:z
	     	 for k = 1:y
		     C[i,j] += A[i,k] * B[k,j] + 1
		 end
	     end
	  end
end

ast_lv_2 = code_typed(test_liveness_2, (Int64,Int64,Int64))[1]

linfo, body = CompilerTools.LambdaHandling.lambdaToLambdaVarInfo(ast_lv_2)

#CompilerTools.CFGs.set_debug_level(4)
#CompilerTools.LivenessAnalysis.set_debug_level(3)
#CompilerTools.LambdaHandling.set_debug_level(3)
lives_2 = CompilerTools.LivenessAnalysis.from_lambda(ast_lv_2)
#=
println( length(lives_2.basic_blocks[lives_2.cfg.basic_blocks[6]].live_in))
println( length(lives_2.basic_blocks[lives_2.cfg.basic_blocks[6]].live_out))
println( length(lives_2.basic_blocks[lives_2.cfg.basic_blocks[6]].def))
println( length(lives_2.basic_blocks[lives_2.cfg.basic_blocks[6]].use))

println( length(lives_2.basic_blocks[lives_2.cfg.basic_blocks[14]].live_in))
println( length(lives_2.basic_blocks[lives_2.cfg.basic_blocks[14]].live_out))
println( length(lives_2.basic_blocks[lives_2.cfg.basic_blocks[14]].def))
println( length(lives_2.basic_blocks[lives_2.cfg.basic_blocks[14]].use))
=#

#@test length(lives_2.basic_blocks[lives_2.cfg.basic_blocks[6]].live_in) == 10
#@test length(lives_2.basic_blocks[lives_2.cfg.basic_blocks[6]].live_out) == 10
#@test length(lives_2.basic_blocks[lives_2.cfg.basic_blocks[6]].def) == 6
#@test length(lives_2.basic_blocks[lives_2.cfg.basic_blocks[6]].use) == 3

#@test length(lives_2.basic_blocks[lives_2.cfg.basic_blocks[14]].live_in) == 13
#@test length(lives_2.basic_blocks[lives_2.cfg.basic_blocks[14]].live_out) == 13
#@test length(lives_2.basic_blocks[lives_2.cfg.basic_blocks[14]].def) == 8
#@test length(lives_2.basic_blocks[lives_2.cfg.basic_blocks[14]].use) == 6

#@test length(lives_2.cfg.basic_blocks) == 26

deps = CompilerTools.TransitiveDependence.computeDependencies(lives_2)

#for entry in deps
#    vd = CompilerTools.LambdaHandling.getVarDef(entry[1], linfo)
#    println(entry[1], " ", entry[2], " ", vd.name)
#end


text:
using Documenter, Stan

DOC_ROOT = joinpath(dirname(pathof(Stan)), "..", "docs")
DocDir =  joinpath(DOC_ROOT, "src")

page_list = Array{Pair{String, Any}, 1}();
append!(page_list, [Pair("Intro", "INTRO.md")]);
append!(page_list, [Pair("Installation", "INSTALLATION.md")]);
append!(page_list, [Pair("Walkthrough", "WALKTHROUGH.md")]);
append!(page_list, [Pair("Walkthrough2", "WALKTHROUGH2.md")]);
append!(page_list, [Pair("Versions", "VERSIONS.md")]);

makedocs(
  format = Documenter.HTML(prettyurls = haskey(ENV, "GITHUB_ACTIONS")),
  root = DOC_ROOT,
  modules = Module[],
  sitename = "Stan.jl",
  authors = "Rob J Goedman",
  pages = page_list,
)

deploydocs(
  root = DOC_ROOT,
  repo = "github.com/StanJulia/Stan.jl.git",
  devbranch = "master",
  push_preview = true,
)


text:
"""
from Arb docs:

A variable of type arf_t holds an arbitrary-precision binary floating-point number:
that is, a rational number of the form x⋅2y where x,y∈Z and x is odd, or one of
the special values zero, plus infinity, minus infinity, or NaN (not-a-number).
There is currently no support for negative zero, unsigned infinity, or a NaN with a payload.
""" ArbFloat

mutable struct ArbFloat{P} <: AbstractFloat    # P is the precision in bits
    exp::Int        # fmpz         exponent of 2 (2^exp)
    size::UInt      # mp_size_t    nwords and sign (lsb holds sign of significand)
    d1::UInt        # significand  unsigned, immediate value or the initial span
    d2::UInt        #   (d1, d2)   the final part indicating the significand, or 0

    function ArbFloat{P}() where {P}
        z = new{P}()
        ccall(@libarb(arf_init), Cvoid, (Ref{ArbFloat},), z)
        finalizer(arf_clear, z)
        return z
    end
end

arf_clear(x::ArbFloat{P}) where {P} = ccall(@libarb(arf_clear), Cvoid, (Ref{ArbFloat},), x)

ArbFloat{P}(x::ArbFloat{P}) where {P} = x
ArbFloat(x::ArbFloat{P}) where {P} = x

ArbFloat{P}(x::Missing) where {P} = missing
ArbFloat(x::Missing) = missing

@inline sign_bit(x::ArbFloat{P}) where {P} = isodd(x.size)

ArbFloat(x, prec::Int) = prec>=MINIMUM_PRECISION ? ArbFloat{workingbits(prec)}(x) : throw(DomainError("bit precision $prec < $MINIMUM_PRECISION"))

function ArbFloat(x::T; bits::Int=0, digits::Int=0, base::Int=iszero(bits) ? 10 : 2) where {T<:Number}
    if base === 10
        digits = digits > 0 ? bits4digits(digits) : (bits > 0 ? bits : DEFAULT_PRECISION.x)
    elseif base === 2
        digits = bits > 0 ? bits : (digits > 0 ? digits : DEFAULT_PRECISION.x)
    else
        throw(ErrorException("base expects 2 or 10"))
    end
    ArbFloat(x, digits)
end

swap(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = ccall(@libarb(arf_swap), Cvoid, (Ref{ArbFloat}, Ref{ArbFloat}), x, y)

function copy(x::ArbFloat{P}) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arf_set), Cvoid, (Ref{ArbFloat}, Ref{ArbFloat}), z, x)
    return z
end

function copy(x::ArbFloat{P}, bitprecision::Int, roundingmode::RoundingMode) where {P}
    z = ArbFloat{P}()
    rounding = match_rounding_mode(roundingmode)
    rounddir = ccall(@libarb(arf_set_round), Cint, (Ref{ArbFloat}, Ref{ArbFloat}, Clong, Cint), z, x, bitprecision, rounding)
    return z
end

copy(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} = copy(x, P, roundingmode)
copy(x::ArbFloat{P}, bitprecision::Int) where {P} = copy(x, bitprecision, RoundNearest)

function ArbFloat{P}(x::Int64) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arf_set_si), Cvoid, (Ref{ArbFloat}, Clong), z, x)
    return z
end
ArbFloat{P}(x::T) where {P, T<:Union{Int8, Int16, Int32}} = ArbFloat{P}(Int64(x))

function ArbFloat{P}(x::UInt64) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arf_set_ui), Cvoid, (Ref{ArbFloat}, Culong), z, x)
    return z
end
ArbFloat{P}(x::T) where {P, T<:Union{UInt8, UInt16, UInt32}} = ArbFloat{P}(UInt64(x))

function ArbFloat{P}(x::Float64) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arf_set_d), Cvoid, (Ref{ArbFloat}, Cdouble), z, x)
    return z
end
ArbFloat{P}(x::T) where {P, T<:Union{Float16, Float32}} = ArbFloat{P}(Float64(x))

function ArbFloat{P}(x::BigFloat) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arf_set_mpfr), Cvoid, (Ref{ArbFloat}, Ref{BigFloat}), z, x)
    return z
end
ArbFloat{P}(x::BigInt) where {P} = ArbFloat{P}(BigFloat(x))
ArbFloat{P}(x::Rational{T}) where {P, T<:Signed} = ArbFloat{P}(BigFloat(x))

function ArbFloat{P}(x::Irrational{S}) where {P,S}
    prec = precision(BigFloat)
    newprec = max(prec, P + 32)
    setprecision(BigFloat, newprec)
    y = BigFloat(x)
    z = ArbFloat{P}(y)
    setprecision(BigFloat, prec)
    return z
end

Int64(x::ArbFloat{P}) where {P} = Int64(x, RoundNearest)
function Int64(x::ArbFloat{P}, roundingmode::RoundingMode) where {P}
    rounding = match_rounding_mode(roundingmode)
    z = ccall(@libarb(arf_get_si), Clong, (Ref{ArbFloat}, Cint), x, rounding)
    return z
end
Int32(x::ArbFloat{P}) where {P} = Int32(Int64(x))
Int32(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} = Int32(Int64(x), roundingmode)
Int16(x::ArbFloat{P}) where {P} = Int16(Int64(x))
Int16(x::ArbFloat{P}, roundingmode::RoundingMode) where {P} = Int16(Int64(x), roundingmode)

BigFloat(x::ArbFloat{P}) where {P} = BigFloat(x, RoundNearest)
function BigFloat(x::ArbFloat{P}, roundingmode::RoundingMode) where {P}
    rounding = match_rounding_mode(roundingmode)
    z = BigFloat(0, workingprecision(x))
    roundingdir = ccall(@libarb(arf_get_mpfr), Cint, (Ref{BigFloat}, Ref{ArbFloat}, Cint), z, x, rounding)
    return z
end
BigFloat(x::ArbFloat{P}, bitprecision::Int) where {P} = BigFloat(x, bitprecision, RoundNearest)
function BigFloat(x::ArbFloat{P}, bitprecision::Int, roundingmode::RoundingMode) where {P}
    rounding = match_rounding_mode(roundingmode)
    z = BigFloat(0, bitprecision)
    roundingdir = ccall(@libarb(arf_get_mpfr), Cint, (Ref{BigFloat}, Ref{ArbFloat}, Cint), z, x, rounding)
    return z
end

BigInt(x::ArbFloat{P}) where {P} = BigInt(trunc(BigFloat(x)))

function Base.Integer(x::ArbFloat{P}) where {P}
    if isinteger(x)
       abs(x) <= typemax(Int64) ? Int64(x) : BigInt(x)
    else
       throw(InexactError("$x"))
    end
end

for (F,A) in ((:floor, :arf_floor), (:ceil, :arf_ceil))
    @eval begin
        function $F(x::ArbFloat{P}) where {P}
            z = ArbFloat{P}()
            ccall(@libarb($A), Cvoid, (Ref{ArbFloat}, Ref{ArbFloat}), z, x)
            return z
        end
        $F(::Type{T}, x::ArbFloat{P}) where {P, T} = T($F(x))
    end
end

trunc(x::ArbFloat{P}) where {P} = signbit(x) ? ceil(x) : floor(x)
trunc(::Type{T}, x::ArbFloat{P}) where {P, T} = T(trunc(x))

midpoint(x::ArbFloat{P}) where {P} = x
radius(x::ArbFloat{P}) where {P} = zero(ArbFloat{P})


text:
raw = read(stdin, String)
raw = convert(Array{Any}, split(raw, ".\n", keepempty=false))

# raw[i] = [left, [[count 1, color 1], ...]]
for (i, r) in enumerate(raw)
    raw[i] = convert(Array{Any}, map(strip, split(r, "contain")))

    # Left
    raw[i][begin] = raw[i][begin][begin:end - length(" bags")]

    # Right
    raw[i][2] = convert(Array{Any}, map(strip, split(raw[i][2], ',')))

    if raw[i][2][begin] == "no other bags"
        raw[i][2] = []
        continue
    end

    for (j, rule) in enumerate(raw[i][2])
        first = findfirst(' ', rule)
        last = findlast(' ', rule)

        count, color = parse(Int, strip(rule[1:first - 1])), strip(rule[first + 1:last - 1])
        raw[i][2][j] = [count, color]
    end
end

# Every color containing a shiny gold
shiny = Set(["shiny gold"])
lastcount = nothing
while lastcount != length(shiny)
    global lastcount = length(shiny)

    for (left, right) in raw
        # If an item in right contains a shiny gold bag, left
        # contains a shiny gold bag
        any((x) -> x[2] ∈ shiny, right) && (push!(shiny, left))
    end
end

# Without shiny gold
println("Part 1 : $(length(shiny) - 1)")

name2prod = Dict(left => right for (left, right) in raw)
name2count = Dict()

# DP approach (memoization)
function count(color)
    global name2count

    if haskey(name2count, color)
        return name2count[color]
    end

    c = 1
    # For each right item
    for (number, name) in name2prod[color]
        c += number * count(name)
    end

    name2count[color] = c

    return c
end

println("Part 2 : $(count("shiny gold") - 1)")


text:
# # RealNVP network
# For the definition of this network and concepts of normalizing flow,
# please refer this realnvp blog: https://lilianweng.github.io/lil-log/2018/10/13/flow-based-deep-generative-models.html,
# and the pytorch notebook: https://github.com/GiggleLiu/marburg/blob/master/solutions/realnvp.ipynb

using NiLang, NiLang.AD
using LinearAlgebra
using DelimitedFiles
using Plots

# `include` the optimizer, you can find it under the `Adam.jl` file in the `examples/` folder.
include(NiLang.project_relative_path("examples", "Adam.jl"))

# ## Model definition
# First, define the single layer transformation and its behavior under `GVar` - the gradient wrapper.
struct RealNVPLayer{T}
    ## transform network
    W1::Matrix{T}
    b1::Vector{T}
    W2::Matrix{T}
    b2::Vector{T}
    y1::Vector{T}
    y1a::Vector{T}

    ## scaling network
    sW1::Matrix{T}
    sb1::Vector{T}
    sW2::Matrix{T}
    sb2::Vector{T}
    sy1::Vector{T}
    sy1a::Vector{T}
end

"""collect parameters in the `layer` into a vector `out`."""
function collect_params!(out, layer::RealNVPLayer)
    k=0
    for field in [:W1, :b1, :W2, :b2, :sW1, :sb1, :sW2, :sb2]
        v = getfield(layer, field)
        nv = length(v)
        out[k+1:k+nv] .= vec(v)
        k += nv
    end
    return out
end

"""dispatch vectorized parameters `out` into the `layer`."""
function dispatch_params!(layer::RealNVPLayer, out)
    k=0
    for field in [:W1, :b1, :W2, :b2, :sW1, :sb1, :sW2, :sb2]
        v = getfield(layer, field)
        nv = length(v)
        vec(v) .= out[k+1:k+nv]
        k += nv
    end
    return out
end

function nparameters(n::RealNVPLayer)
    sum(x->length(getfield(n, x)), [:W1, :b1, :W2, :b2, :sW1, :sb1, :sW2, :sb2])
end

# Then, we define `network` and how to access the parameters.
const RealNVP{T} = Vector{RealNVPLayer{T}}

nparameters(n::RealNVP) = sum(nparameters, n)

function collect_params(n::RealNVP{T}) where T
    out = zeros(T, nparameters(n))
    k = 0
    for layer in n
        np = nparameters(layer)
        collect_params!(view(out, k+1:k+np), layer)
        k += np
    end
    return out
end

function dispatch_params!(network::RealNVP, out)
    k = 0
    for layer in network
        np = nparameters(layer)
        dispatch_params!(layer, view(out, k+1:k+np))
        k += np
    end
    return network
end

function random_realnvp(nparams::Int, nhidden::Int, nhidden_s::Int, nlayer::Int; scale=0.1)
    random_realnvp(Float64, nparams, nhidden, nhidden_s::Int, nlayer; scale=scale)
end

function random_realnvp(::Type{T}, nparams::Int, nhidden::Int, nhidden_s::Int, nlayer::Int; scale=0.1) where T
    nin = nparams÷2
    scale = T(scale)
    y1 = zeros(T, nhidden)
    sy1 = zeros(T, nhidden_s)
    RealNVPLayer{T}[RealNVPLayer(
            randn(T, nhidden, nin)*scale, randn(T, nhidden)*scale,
            randn(T, nin, nhidden)*scale, randn(T, nin)*scale, y1, zero(y1),
            randn(T, nhidden_s, nin)*scale, randn(T, nhidden_s)*scale,
            randn(T, nin, nhidden_s)*scale, randn(T, nin)*scale, sy1, zero(sy1),
            ) for _ = 1:nlayer]
end

# ## Loss function
#
# In each layer, we use the information in `x` to update `y!`.
# During computing, we use to vector type ancillas `y1` and `y1a`,
# both of them can be uncomputed at the end of the function.

@i function onelayer!(x::AbstractVector{T}, layer::RealNVPLayer{T},
                y!::AbstractVector{T}, logjacobian!::T; islast) where T
    @routine @invcheckoff begin
        ## scale network
        scale ← zero(y!)
        ytemp2 ← zero(y!)
        i_affine!(layer.sy1, layer.sW1, layer.sb1, x)
        @inbounds for i=1:length(layer.sy1)
            if (layer.sy1[i] > 0, ~)
                layer.sy1a[i] += layer.sy1[i]
            end
        end
        i_affine!(scale, layer.sW2, layer.sb2, layer.sy1a)

        ## transform network
        i_affine!(layer.y1, layer.W1, layer.b1, x)
        ## relu
        @inbounds for i=1:length(layer.y1)
            if (layer.y1[i] > 0, ~)
                layer.y1a[i] += layer.y1[i]
            end
        end
    end
    ## inplace multiply exp of scale! -- dangerous
    @inbounds @invcheckoff for i=1:length(scale)
        @routine begin
            expscale ← zero(T)
            tanhscale ← zero(T)
            if (islast, ~)
                tanhscale += tanh(scale[i])
            else
                tanhscale += scale[i]
            end
            expscale += exp(tanhscale)
        end
        logjacobian! += tanhscale
        ## inplace multiply!!!
        temp ← zero(T)
        temp += y![i] * expscale
        SWAP(temp, y![i])
        temp -= y![i] / expscale
        temp → zero(T)
        ~@routine
    end

    ## affine the transform layer
    i_affine!(y!, layer.W2, layer.b2, layer.y1a)
    ~@routine
    ## clean up accumulated rounding error, since this memory is reused.
    @safe layer.y1 .= zero(T)
    @safe layer.sy1 .= zero(T)
end

# A realnvp network always transforms inputs reversibly.
# We update one half of `x!` a time, so that input and output memory space do not clash.
@i function realnvp!(x!::AbstractVector{T}, network::RealNVP{T}, logjacobian!) where T
    @invcheckoff for i=1:length(network)
        np ← length(x!)
        if (i%2==0, ~)
            @inbounds onelayer!(x! |> subarray(np÷2+1:np), network[i], x! |> subarray(1:np÷2), logjacobian!; islast=i==length(network))
        else
            @inbounds onelayer!(x! |> subarray(1:np÷2), network[i], x! |> subarray(np÷2+1:np), logjacobian!; islast=i==length(network))
        end
    end
end

# How to obtain the log-probability of a data.

@i function logp!(out!::T, x!::AbstractVector{T}, network::RealNVP{T}) where T
    (~realnvp!)(x!, network, out!)
    @invcheckoff for i = 1:length(x!)
        @routine begin
            xsq ← zero(T)
            @inbounds xsq += x![i]^2
        end
        out! -= 0.5 * xsq
        ~@routine
    end
end

# The negative-log-likelihood loss function

@i function nll_loss!(out!::T, cum!::T, xs!::Matrix{T}, network::RealNVP{T}) where T
    @invcheckoff for i=1:size(xs!, 2)
        @inbounds logp!(cum!, xs! |> subarray(:,i), network)
    end
    out! -= cum!/(@const size(xs!, 2))
end

# ## Training

function train(x_data, model; num_epochs = 800)
    num_vars = size(x_data, 1)
    params = collect_params(model)
    optimizer = Adam(; lr=0.01)
    for epoch = 1:num_epochs
        loss, a, b, c = nll_loss!(0.0, 0.0, copy(x_data), model)
        if epoch % 50 == 1
            println("epoch = $epoch, loss = $loss")
            display(showmodel(x_data, model))
        end
        _, _, _, gmodel = (~nll_loss!)(GVar(loss, 1.0), GVar(a), GVar(b), GVar(c))
        g = grad.(collect_params(gmodel))
        update!(params, grad.(collect_params(gmodel)), optimizer)
        dispatch_params!(model, params)
    end
    return model
end

function showmodel(x_data, model; nsamples=2000)
    scatter(x_data[1,1:nsamples], x_data[2,1:nsamples]; xlims=(-5,5), ylims=(-5,5))
    zs = randn(2, nsamples)
    for i=1:nsamples
        realnvp!(view(zs, :, i), model, 0.0)
    end
    scatter!(zs[1,:], zs[2,:])
end

# you can find the training data in `examples/` folder
x_data = Matrix(readdlm(NiLang.project_relative_path("examples", "train.dat"))')

import Random; Random.seed!(22)
model = random_realnvp(Float64, size(x_data, 1), 10, 10, 4; scale=0.1)

# Before training, the distribution looks like
# ![before](../../asset/nice_before.png)
model = train(x_data, model; num_epochs=800)

# After training, the distribution looks like
# ![before](../../asset/realnvp_after.png)


text:
struct VelocityFields <: FieldSet
    u::FaceFieldX
    v::FaceFieldY
    w::FaceFieldZ
end

struct TracerFields <: FieldSet
    T::CellField
    S::CellField
end

struct PressureFields <: FieldSet
    pHY′::CellField
    pNHS::CellField
end

struct SourceTerms <: FieldSet
    Gu::FaceFieldX
    Gv::FaceFieldY
    Gw::FaceFieldZ
    GT::CellField
    GS::CellField
end

struct StepperTemporaryFields <: FieldSet
    fC1::CellField
    fCC1::CellField
    fCC2::CellField
end

function VelocityFields(arch::Architecture, grid::Grid)
    u = FaceFieldX(arch, grid)
    v = FaceFieldY(arch, grid)
    w = FaceFieldZ(arch, grid)
    VelocityFields(u, v, w)
end

function TracerFields(arch::Architecture, grid::Grid)
    θ = CellField(arch, grid)  # Temperature θ to avoid conflict with type T.
    S = CellField(arch, grid)
    TracerFields(θ, S)
end

function PressureFields(arch::Architecture, grid::Grid)
    pHY′ = CellField(arch, grid)
    pNHS = CellField(arch, grid)
    PressureFields(pHY′, pNHS)
end

function SourceTerms(arch::Architecture, grid::Grid)
    Gu = FaceFieldX(arch, grid)
    Gv = FaceFieldY(arch, grid)
    Gw = FaceFieldZ(arch, grid)
    GT = CellField(arch, grid)
    GS = CellField(arch, grid)
    SourceTerms(Gu, Gv, Gw, GT, GS)
end

function StepperTemporaryFields(arch::Architecture, grid::Grid)
    fC1 = CellField(arch, grid)

    # Forcing Float64 for these fields as it's needed by the Poisson solver.
    fCC1 = CellField(Complex{Float64}, arch, grid)
    fCC2 = CellField(Complex{Float64}, arch, grid)
    StepperTemporaryFields(fC1, fCC1, fCC2)
end


text:
using DifferentialEquations, Test

f = function (r, yp, y, p, tres)
    r[1]  = -0.04*y[1] + 1.0e4*y[2]*y[3]
    r[2]  = -r[1] - 3.0e7*y[2]*y[2] - yp[2]
    r[1] -=  yp[1]
    r[3]  =  y[1] + y[2] + y[3] - 1.0
end
u0 = [1.0, 0, 0]
du0 = [-0.04, 0.04, 0.0]
prob_dae_resrob = DAEProblem(f,du0,u0,(0.0,100000.0))

prob = prob_dae_resrob
sol = solve(prob)

@test typeof(sol.alg) <: DifferentialEquations.Sundials.IDA

max_stars_repo_name: ChrisRackauckas/DiffEq

max_issues_repo_name: ChrisRackauckas/DiffEq


text:
module PageTest

using Base.Test
using PageDown

pages = ["1", "2", "3", "4", "5"]
page = Page(pages)

@test 1 == page.start
@test 1 == page.step

next!(page)
@test 2 == page.start

prev!(page)
@test 1 == page.start

prev!(page)
@test 1 == page.start

last!(page)
@test page.count == page.start

first!(page)
@test 1 == page.start

@test (page,["1"]) == proper(page)

# coverage
page = Page(pages, 1, 10)
@test_throws ArgumentError go!(page, 10)
last!(page)
@open 1
Base.Markdown.term(STDOUT, @first)

pages = ["1", "2", "3", "4", "5"]
page = Page(pages, 1, 1, Dict())

option = Dict(:IMG=>"IMG", :imgcmd=>"ls", :imgcat=>"", :imgfolder=>"", :opencmd=>"echo")
using Base.Markdown
pages = [md"""
         # img
          - IMG sample.png
         """]
page = Page(pages, option)
@next

pages = [md"""
         # github
          - http://github.com
         """]
page = Page(pages, option)
@open 1

end


text:
push!(LOAD_PATH, "/Users/sean/GitRepos/ParticleMapper/ParticleGridMapper.jl")
push!(LOAD_PATH, "/Users/sean/GitRepos/Chia-Yu_Hu_ParticleMapper/OctreeBH.jl")

using ParticleGridMapper
using StaticArrays
using PyPlot

const N = 3 #spatial dimension
const T = Float64
const Npart = 1000 #number of particles
const BOXSIZE = 1.0

using Random
Random.seed!(1114)

function get_spherical_cloud(Ngas, rmax)
    X = Vector{SVector{N,T}}(undef, Ngas)
    for i in 1:Ngas
        r3d = rmax * 2 #s.t. the while-loop will be evaluated at least once
        while r3d > rmax
            m_frac = rand() * 0.99;      #m_frac = 0 ~ 0.99
            r3d = rmax * m_frac;
        end
        phi       = ( 2. * rand() - 1. ) * pi
        cos_theta = ( 2. * rand() - 1. )
        x = r3d * sqrt( 1. - cos_theta^2 ) * cos(phi)
        y = r3d * sqrt( 1. - cos_theta^2 ) * sin(phi)
        z = r3d * cos_theta
        X[i] = @SVector [x,y,z]
    end
    return X
end

X = get_spherical_cloud(Npart, 0.5*BOXSIZE)
topnode_length = @SVector(ones(N)) * BOXSIZE  #actual length of tree

center = @SVector [0.,0.,0.]

Mtot = 1.0
mass = ones(Npart) .* (Mtot / Npart) #particle mass

radius = ParticleGridMapper.norm.(X./BOXSIZE)
rho = 2*Mtot/(4*pi) .* radius.^(-2)
volume = mass ./ rho

Nngb = 32
hsml = ( ( Nngb * (Mtot / Npart) / (4*pi/3) ) ./ rho).^(1/3)

part = [DataP2G{N,T}(SVector(X[i]), i, hsml[i], mass[i], 0.0, T[]) for i in eachindex(X)]

#build the tree
treeNgb = buildtree(part, center, topnode_length);

boxsizes = @SVector(ones(N)) * BOXSIZE  #for periodic B.C.

nx = ny = 128
ix,iy = 1,2

xmin = (-0.5, -0.5, -0.5).* BOXSIZE
xmax = (0.5, 0.5, 0.5).* BOXSIZE

println("Cartensian NGP...")
image_cart_NGP = map_particle_to_3Dgrid_NGP_thread(rho, mass, X, (nx,ny,nx), xmin, xmax);
dimlos = filter!(x->(x!=ix)&&(x!=iy), [1,2,3])[1]
image_cart_NGP=dropdims(sum(image_cart_NGP, dims=dimlos),dims=dimlos);

println("Cartensian SPH...")
image_cart_SPH = map_particle_to_2Dgrid_loopP_thread(rho, mass./rho, X, hsml, xmin, xmax,
    ngrids=(nx,ny,nx), xaxis=ix, yaxis=iy, column=true, pbc = (true, true, true));

max_depth = 10
max_depth_true = get_max_tree_depth(treeNgb)
@show max_depth, max_depth_true

treeAMR = deepcopy(treeNgb)
set_max_depth_AMR!(treeAMR, max_depth)
balance_all_level!(treeAMR)

println("AMR SPH...")
@time map_particle_to_AMRgrid_SPH!(treeAMR, rho, volume, X, hsml, boxsizes, treeNgb=treeNgb)
@time map_particle_to_AMRgrid_SPH!(treeAMR, rho, volume, X, hsml, boxsizes)
rhoAMR_SPH = get_AMRfield(treeAMR)
image_SPH = project_AMRgrid_to_image(nx, ny, ix, iy, treeAMR, center, boxsizes);

#treeNgb = buildtree(part, center, topnode_length);
#treeAMR = deepcopy(treeNgb)
#set_max_depth_AMR!(treeAMR, max_depth)
#balance_all_level!(treeAMR)
#
#@time for i in 1:length(X)
#    map_particle_to_AMRgrid_loopP_recursive!(treeAMR, rho[i], volume[i], hsml[i], X[i], boxsizes, 100, BOXSIZE)
#end
#rhoAMR_SPH_2 = get_AMRfield(treeAMR, max_depth=max_depth)
#image_SPH_2 = project_AMRgrid_to_image_thread(nx, ny, ix, iy, treeAMR, center, boxsizes);
#@show all(rhoAMR_SPH_2.≈rhoAMR_SPH)
#@show all(image_SPH_2.≈image_SPH)

println("AMR MFM...")
map_particle_to_AMRgrid_MFM!(treeAMR, rho, X, hsml, boxsizes)
rhoAMR_MFM = get_AMRfield(treeAMR)
image_MFM = project_AMRgrid_to_image(nx, ny, ix, iy, treeAMR, center, boxsizes);

println("AMR NGP...")
map_particle_to_AMRgrid_NGP!(treeAMR, rho)
rhoAMR_NGP = get_AMRfield(treeAMR)
image_NGP = project_AMRgrid_to_image(nx, ny, ix, iy, treeAMR, center, boxsizes);

clf()
fig, ax = subplots(2, 3, figsize=(15, 10))
ax[1].set_title("particle distribution")
ax[1].plot(getindex.(X,1), getindex.(X,2), ".", ms=3, color="red")
ax[1].axis([-0.5*BOXSIZE, 0.5*BOXSIZE, -0.5*BOXSIZE, 0.5*BOXSIZE])
ax[2].set_title("AMR MFM")
ax[2].imshow(log10.(image_MFM'), origin="lower", vmin=-1, vmax=1.5, interpolation="none", extent=(-0.5*BOXSIZE, 0.5*BOXSIZE, -0.5*BOXSIZE, 0.5*BOXSIZE))
ax[3].set_title("Cartesian NGP")
ax[3].imshow(log10.(image_cart_NGP'), origin="lower", vmin=-1, vmax=1.5, interpolation="none", extent=(-0.5*BOXSIZE, 0.5*BOXSIZE, -0.5*BOXSIZE, 0.5*BOXSIZE))
ax[4].set_title("AMR NGP")
ax[4].imshow(log10.(image_NGP'), origin="lower", vmin=-1, vmax=1.5, interpolation="none", extent=(-0.5*BOXSIZE, 0.5*BOXSIZE, -0.5*BOXSIZE, 0.5*BOXSIZE))
ax[5].set_title("Cartesian SPH")
ax[5].imshow(log10.(image_cart_SPH'), origin="lower", vmin=-1, vmax=1.5, interpolation="none", extent=(-0.5*BOXSIZE, 0.5*BOXSIZE, -0.5*BOXSIZE, 0.5*BOXSIZE))
ax[6].set_title("AMR SPH")
ax[6].imshow(log10.(image_SPH'), origin="lower", vmin=-1, vmax=1.5, interpolation="none", extent=(-0.5*BOXSIZE, 0.5*BOXSIZE, -0.5*BOXSIZE, 0.5*BOXSIZE))

fig.tight_layout()
savefig("compare_MFM_SPH_NGP_cloud.png")


text:
# TODO: this distribution may need clean-up
"""
    PoissonBinomial(p)

A *Poisson-binomial distribution* describes the number of successes in a sequence of independent trials, wherein each trial has a different success rate.
It is parameterized by a vector `p` (of length ``K``), where ``K`` is the total number of trials and `p[i]` corresponds to the probability of success of the `i`th trial.

```math
P(X = k) = \\sum\\limits_{A\\in F_k} \\prod\\limits_{i\\in A} p[i] \\prod\\limits_{j\\in A^c} (1-p[j]), \\quad \\text{ for } k = 0,1,2,\\ldots,K,
```

where ``F_k`` is the set of all subsets of ``k`` integers that can be selected from ``\\{1,2,3,...,K\\}``.

```julia
PoissonBinomial(p)   # Poisson Binomial distribution with success rate vector p

params(d)            # Get the parameters, i.e. (p,)
succprob(d)          # Get the vector of success rates, i.e. p
failprob(d)          # Get the vector of failure rates, i.e. 1-p
```

External links:

* [Poisson-binomial distribution on Wikipedia](http://en.wikipedia.org/wiki/Poisson_binomial_distribution)

"""
struct PoissonBinomial{T<:Real} <: DiscreteUnivariateDistribution
    p::Vector{T}
    pmf::Vector{T}

    function PoissonBinomial{T}(p::AbstractArray) where {T <: Real}
        pb = poissonbinomial_pdf_fft(p)
        @assert isprobvec(pb)
        new{T}(p, pb)
    end
end

function PoissonBinomial(p::AbstractArray{T}) where {T <: Real}
    for i in eachindex(p)
        @check_args(PoissonBinomial, 0 <= p[i] <= 1)
    end
    PoissonBinomial{T}(p)
end

PoissonBinomial(p::AbstractArray{T}, ::NoArgCheck) where {T<:Real} = PoissonBinomial{T}(p)

@distr_support PoissonBinomial 0 length(d.p)

#### Conversions

function PoissonBinomial(::Type{PoissonBinomial{T}}, p::Vector{S}) where {T, S}
    return PoissonBinomial(Vector{T}(p))
end
function PoissonBinomial(::Type{PoissonBinomial{T}}, d::PoissonBinomial{S}) where {T, S}
    return PoissonBinomial(Vector{T}(d.p), NoArgCheck())
end

#### Parameters

ntrials(d::PoissonBinomial) = length(d.p)
succprob(d::PoissonBinomial) = d.p
failprob(d::PoissonBinomial{T}) where {T} = one(T) .- d.p

params(d::PoissonBinomial) = (d.p,)
@inline partype(::PoissonBinomial{T}) where {T} = T

#### Properties

mean(d::PoissonBinomial) = sum(succprob(d))
var(d::PoissonBinomial) = sum(succprob(d) .* failprob(d))

function skewness(d::PoissonBinomial{T}) where {T}
    v = zero(T)
    s = zero(T)
    p,  = params(d)
    for i in eachindex(p)
        v += p[i] * (one(T) - p[i])
        s += p[i] * (one(T) - p[i]) * (one(T) - T(2) * p[i])
    end
    return s / sqrt(v) / v
end

function kurtosis(d::PoissonBinomial{T}) where {T}
    v = zero(T)
    s = zero(T)
    p,  = params(d)
    for i in eachindex(p)
        v += p[i] * (one(T) - p[i])
        s += p[i] * (one(T) - p[i]) * (one(T) - T(6) * (one(T) - p[i]) * p[i])
    end
    s / v / v
end

entropy(d::PoissonBinomial) = entropy(Categorical(d.pmf))
median(d::PoissonBinomial) = median(Categorical(d.pmf)) - 1
mode(d::PoissonBinomial) = argmax(d.pmf) - 1
modes(d::PoissonBinomial) = [x  - 1 for x in modes(Categorical(d.pmf))]

#### Evaluation

quantile(d::PoissonBinomial, x::Float64) = quantile(Categorical(d.pmf), x) - 1

function mgf(d::PoissonBinomial{T}, t::Real) where {T}
    p,  = params(d)
    prod(one(T) .- p .+ p .* exp(t))
end

function cf(d::PoissonBinomial{T}, t::Real) where {T}
    p,  = params(d)
    prod(one(T) .- p .+ p .* cis(t))
end

pdf(d::PoissonBinomial, k::Integer) = insupport(d, k) ? d.pmf[k+1] : 0
function logpdf(d::PoissonBinomial{T}, k::Int) where T<:Real
    insupport(d, k) ? log(d.pmf[k + 1]) : -T(Inf)
end

# Computes the pdf of a poisson-binomial random variable using
# fast fourier transform
#
#     Hong, Y. (2013).
#     On computing the distribution function for the Poisson binomial
#     distribution. Computational Statistics and Data Analysis, 59, 41–51.
#
function poissonbinomial_pdf_fft(p::AbstractArray)
    n = length(p)
    ω = 2 / (n + 1)

    x = Vector{Complex{Float64}}(undef, n+1)
    lmax = ceil(Int, n/2)
    x[1] = 1/(n + 1)
    for l=1:lmax
        logz = 0.
        argz = 0.
        for j=1:n
            zjl = 1 - p[j] + p[j] * cospi(ω*l) + im * p[j] * sinpi(ω * l)
            logz += log(abs(zjl))
            argz += atan(imag(zjl), real(zjl))
        end
        dl = exp(logz)
        x[l + 1] = dl * cos(argz) / (n + 1) + dl * sin(argz) * im / (n + 1)
        if n + 1 - l > l
            x[n + 1 - l + 1] = conj(x[l + 1])
        end
    end
    [max(0, real(xi)) for xi in _dft(x)]
end

# A simple implementation of a DFT to avoid introducing a dependency
# on an external FFT package just for this one distribution
function _dft(x::Vector{T}) where T
    n = length(x)
    y = zeros(complex(float(T)), n)
    @inbounds for j = 0:n-1, k = 0:n-1
        y[k+1] += x[j+1] * cis(-π * float(T)(2 * mod(j * k, n)) / n)
    end
    return y
end

#### Sampling

sampler(d::PoissonBinomial) = PoissBinAliasSampler(d)


text:
# Autogenerated wrapper script for libsharp2_jll for i686-w64-mingw32
export libsharp2

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("libsharp2")
JLLWrappers.@declare_library_product(libsharp2, "libsharp2-0.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libsharp2,
        "bin\\libsharp2-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import Test: @test, @testset, @test_throws, @test_broken
import Luna: Maths, Grid, Fields
import Dierckx
import HCubature: hquadrature
import Random: seed!
import FFTW

@testset "Derivatives" begin
    f(x) = @. 4x^3 + 3x^2 + 2x + 1

    @test isapprox(Maths.derivative(f, 1, 1), 12+6+2)
    @test isapprox(Maths.derivative(f, 1, 2), 24+6)
    @test isapprox(Maths.derivative(f, 1, 3), 24)

    e(x) = @. exp(x)

    x = [1, 2, 3, 4, 5]
    @test isapprox(Maths.derivative(e, 1, 5), exp(1), rtol=1e-6)
    @test isapprox(Maths.derivative.(e, x, 5), exp.(x), rtol=1e-6)

    @test isapprox(Maths.derivative(x -> exp.(2x), 1, 1), 2*exp(2))
    @test isapprox(Maths.derivative(x -> exp.(2x), 1, 2), 4*exp(2))
    @test isapprox(Maths.derivative(x -> exp.(-x.^2), 0, 1), 0, atol=1e-14)
end

@testset "Moments" begin
    x = collect(range(-10, stop=10, length=513))
    y = Maths.gauss.(x, 1, x0=1)
    @test Maths.moment(x, y) ≈ 1
    @test Maths.moment(x, y, 2) ≈ 2
    @test Maths.rms_width(x, y) ≈ 1

    x0 = [-2.5, -1.0, 0.0, 1.0, 2.5]
    σ = [0.1, 0.5, 1.0, 1.5, 1.5]
    y = zeros(length(x), length(x0))
    for ii = 1:length(x0)
        y[:, ii] = Maths.gauss.(x, σ[ii], x0=x0[ii])
    end
    @test_throws DomainError Maths.moment(x, y, dim=2)
    @test all(isapprox.(transpose(Maths.moment(x, y, dim=1)), x0, atol=1e-5))
    @test all(isapprox.(transpose(Maths.rms_width(x, y, dim=1)), σ, atol=1e-5))

    yt = transpose(y)
    @test_throws DomainError Maths.moment(x, yt, dim=1)
    xm = Maths.moment(x, yt, dim=2)
    @test all(isapprox.(xm, x0, atol=1e-5))
end

@testset "Fourier" begin
    t = collect(range(-10, stop=10, length=513))
    Et = Maths.gauss.(t, fwhm=4).*cos.(4*t)
    EtA = Maths.hilbert(Et)
    @test maximum(abs.(EtA)) ≈ 1
    @test all(isapprox.(real(EtA), Et, atol=1e-9))

    hilbert! = Maths.plan_hilbert!(Et)
    out = complex(Et)
    hilbert!(out, Et)
    @test all(out .≈ EtA)

    hilbert = Maths.plan_hilbert(Et)
    out = hilbert(Et)
    @test all(out .≈ EtA)

    t = collect(range(-10, stop=10, length=512))
    Et = Maths.gauss.(t, fwhm=4).*cos.(4*t)
    to, Eto = Maths.oversample(t, Et, factor=4)
    @test 4*size(Et)[1] == size(Eto)[1]
    @test all(isapprox.(Eto[1:4:end], Et, rtol=1e-6))

    Etc = Maths.gauss.(t, fwhm=4).*exp.(1im*4*t)
    to, Etco = Maths.oversample(t, Etc, factor=4)
    @test 4*size(Etc)[1] == size(Etco)[1]
    @test all(isapprox.(Etco[1:4:end], Etc, rtol=1e-6))
end

@testset "integration" begin
    x = collect(range(0, stop=8π, length=2^14)).*1e-15
    y = cos.(x.*1e15)
    yi = sin.(x.*1e15)./1e15
    yic = similar(yi)
    yic2 = copy(y)
    Maths.cumtrapz!(yic, y, x)
    Maths.cumtrapz!(yic2, x[2]-x[1])
    @test isapprox(yi, yic, rtol=1e-6)
    @test isapprox(yi, yic2, rtol=1e-6)

    ω = [1e15, 2e15]'
    y = cos.(x.*ω)
    yi =  sin.(x.*ω)./ω
    yic = similar(y)
    yic2 = copy(y)
    Maths.cumtrapz!(yic, y, x)
    Maths.cumtrapz!(yic2, x[2]-x[1])
    @test isapprox(yi, yic, rtol=1e-6)
    @test isapprox(yi, yic2, rtol=1e-6)
end

@testset "series" begin
    sumfunc(x, n) = x + 1/factorial(n)
    e, succ, steps = Maths.aitken_accelerate(sumfunc, 0, rtol=1e-10)
    e2, succ, steps = Maths.converge_series(sumfunc, 0, rtol=1e-10)
    @test isapprox(e, exp(1), rtol=1e-10)
    @test isapprox(e, e2, rtol=1e-10)
    sumfunc(x, n) = x + 1/2^n
    o, succ, steps = Maths.aitken_accelerate(sumfunc, 0, n0=1, rtol=1e-10)
    @test isapprox(o, 1, rtol=1e-10)
    serfunc(x, n) = (x + 2/x)/2
    sqrt2, succ, steps = Maths.aitken_accelerate(serfunc, 1, rtol=1e-10)
    @test isapprox(sqrt2, sqrt(2), rtol=1e-10)
end

@testset "windows" begin
    x = collect(range(-10, stop=10, length=2048))
    pl = Maths.planck_taper(x, -5, -4, 7, 8)
    @test all(pl[x .< -5] .== 0)
    @test all(pl[-4 .< x .< 7] .== 1)
    @test all(pl[8 .< x] .== 0)
end

@testset "fwhm" begin
    fw = 0.2
    x = collect(range(-1, stop=1, length=2048))
    δx = x[2] - x[1]
    y = Maths.gauss.(x, fwhm=fw)
    @test isapprox(Maths.fwhm(x, y), fw, rtol=1e-4)
    @test isapprox(Maths.fwhm(x, y, method=:spline), fw, rtol=1e-4)
    @test Maths.fwhm(x, y, method=:spline) == Maths.fwhm(x, y, method=:spline, minmax=:max)
    @test abs(Maths.fwhm(x, y, method=:nearest) - fw) < δx

    fw = 0.1
    sep = 0.5
    y = Maths.gauss.(x, fwhm=fw, x0=-sep/2) + Maths.gauss.(x, fwhm=fw, x0=sep/2)
    @test isapprox(Maths.fwhm(x, y, minmax=:min), fw, rtol=1e-4)
    @test isapprox(Maths.fwhm(x, y, method=:spline, minmax=:min), fw, rtol=1e-4)
    @test isapprox(Maths.fwhm(x, y, minmax=:max), fw+sep, rtol=1e-4)
    @test isapprox(Maths.fwhm(x, y, method=:spline, minmax=:max), fw+sep, rtol=1e-4)

    fw = 1:10
    x = collect(range(-50, 50; length=2^14))
    y = zeros((length(x), length(fw)))
    for (idx, fwi) in enumerate(fw)
        y[:, idx] .= Maths.gauss.(x; fwhm=fwi)
    end
    @test all([isapprox(f, fwi; rtol=1e-4) for (f, fwi) in zip(Maths.fwhm(x, y; dim=1), fw)])
    @test all([isapprox(f, fwi; rtol=1e-4) for (f, fwi) in zip(Maths.fwhm(x, y'; dim=2), fw)])

    hw = 1
    f(x) = Maths.gauss.(x, fwhm=2*hw)
    @test Maths.hwhm(f, 0) ≈ hw
    @test Maths.hwhm(f, 0; direction=:bwd) ≈ hw

    # Test that FWHM max still works with noise-burst pulses
    x = collect(range(-0.5, stop=1.5, length=2^14))
    y = zero(x)
    seed!(123)
    for i = 1:10000
        y .+= rand()*Maths.gauss.(x, x0=rand(), fwhm=rand()/1000)
    end
    m = 1.1*maximum(y)
    y .+= m*Maths.gauss.(x, x0=-0.05, fwhm=0.05)
    y .+= m*Maths.gauss.(x, x0=1.05, fwhm=0.05)
    @test isapprox(Maths.fwhm(x, y, minmax=:max), 0.05 + 1.1, rtol=1e-4)
    @test isapprox(Maths.fwhm(x, y, minmax=:max, method=:spline), 0.05 + 1.1, rtol=1e-4)
end

@testset "CSpline" begin
    import Random: shuffle
    x = range(0.0, 2π, length=100)
    y = sin.(x)
    spl = Maths.CSpline(x, y);
    fslow(x0) = x0 <= spl.x[1] ? 2 :
                x0 >= spl.x[end] ? length(spl.x) :
                findfirst(x -> x>x0, spl.x)
    ff = Maths.FastFinder(x)
    @test_throws ErrorException Maths.FastFinder(x[end:-1:1])
    @test_throws ErrorException Maths.FastFinder(shuffle(x))
    @test_throws ErrorException Maths.FastFinder(vcat(x[1], x))
    @test all(abs.(spl.(x) .- y) .< 5e-18)
    x2 = range(0.0, 2π, length=300)
    idcs = spl.ifun.(x2)
    idcs_slow = fslow.(x2)
    idcs_ff = ff.(x2)
    idcs_ff_bw = ff.(x2[end:-1:1])
    @test idcs == idcs_slow
    @test idcs_ff == idcs_slow
    @test idcs_ff_bw == idcs_slow[end:-1:1]
    for i = 1:10
        x2r = shuffle(x2)
        @test ff.(x2r) == fslow.(x2r)
    end
    # Create new FastFinder, immediately index backwards - does this still work?
    ff = Maths.FastFinder(x)
    @test ff.(x2[end:-1:1]) == idcs_slow[end:-1:1]
    # Accuracy
    @test maximum(spl.(x2) - sin.(x2)) < 5e-8
    @test abs(Maths.derivative(spl, 1.3, 1) - cos(1.3)) < 1.7e-7
    @test maximum(cos.(x2) - Maths.derivative.(spl, x2, 1)) < 2.1e-6
    # Extrapolation
    ff = Maths.FastFinder(x)
    x3 = range(-0.5, 2π+0.5, length=200)
    @test ff.(x3[end:-1:1]) == fslow.(x3[end:-1:1])
    @test ff.(x3) == fslow.(x3)
    spl = Maths.CSpline(x, y; bounds_error=true)
    @test_throws DomainError spl(-0.5)
    @test_throws DomainError spl(10)
end

@testset "BSpline" begin
    x = range(0.0, 2π, length=100)
    y = sin.(x)
    spl = Maths.BSpline(x, y)
    @test all(abs.(spl.(x) .- y) .< 3e-16)
    x2 = range(0.0, 2π, length=300)
    @test maximum(spl.(x2) - sin.(x2)) < 5e-8
    # these use the actual spline derivative
    @test abs(Maths.derivative(spl, 1.3, 1) - cos(1.3)) < 1.7e-7
    @test maximum(cos.(x2) - Maths.derivative.(spl, x2, 1)) < 2.1e-6
    # test second derivative
    @test maximum(-sin.(x2) - Maths.derivative.(spl, x2, 2)) < 2.0e-4
    # test direct finite differences
    @test abs(invoke(Maths.derivative, Tuple{Any,Any,Integer}, spl, 1.3, 1) - cos(1.3)) < 1.7e-7
    @test maximum(cos.(x2) .- invoke.(Maths.derivative, Tuple{Any,Any,Integer}, spl, x2, 1)) < 2.1e-6
    # test roots
    yr = x.^2 .- 1.0
    splr = Maths.BSpline(x, yr)
    @test Maths.roots(splr) == [1.0]
    # test complex
    yi = sin.(x .+ π/6)
    yc = complex.(y, yi)
    splc = Maths.BSpline(x, complex.(y, yi))
    @test all(abs.(splc.(x) .- yc) .< 5e-16)
    @test maximum(abs.(splc.(x2) .- complex.(sin.(x2), sin.(x2 .+ π/6)))) < 2.6e-7
    @test abs(Maths.derivative(splc, 1.3, 1) - complex(cos(1.3), cos(1.3 + π/6))) < 2.5e-7
    @test abs(Maths.derivative(splc, 1.3, 2) - complex(-sin(1.3), -sin(1.3 + π/6))) < 2.5e-3
    # test Julia evaluation vs original Dierckx
    @test all(spl.(x2) .== spl.rspl.(x2))
    # test full spline Derivatives
    spl1 = Maths.differentiate_spline(spl, 1)
    @test maximum(abs.(cos.(x2) .- spl1.(x2))) < 5.1e-6
    spl2 = Maths.differentiate_spline(spl, 2)
    @test isapprox(-sin.(x2), spl2.(x2),  rtol=3e-5)
    spl3 = Maths.differentiate_spline(spl, 3)
    @test isapprox(-cos.(x2), spl3.(x2),  rtol=8e-4)
    # test other orders
    spl_k1 = Maths.BSpline(x, y, order=1)
    @test all(abs.(spl_k1.(x) .- y) .< 3e-16)
    @test maximum(spl_k1.(x2) - sin.(x2)) < 6e-4
    spl_k2 = Maths.BSpline(x, y, order=2)
    @test all(abs.(spl_k2.(x) .- y) .< 3.4e-16)
    @test maximum(spl_k2.(x2) - sin.(x2)) < 2e-5
    spl_k4 = Maths.BSpline(x, y, order=4)
    @test all(abs.(spl_k4.(x) .- y) .< 4.5e-16)
    @test maximum(spl_k4.(x2) - sin.(x2)) < 2e-8
    spl_k5 = Maths.BSpline(x, y, order=5)
    @test all(abs.(spl_k5.(x) .- y) .< 4.5e-16)
    @test maximum(spl_k5.(x2) - sin.(x2)) < 1.1e-10
    # check 5th order spline derivatives
    @test abs(Maths.derivative(spl_k5, 1.3, 1) - cos(1.3)) < 2e-11
    @test abs(Maths.derivative(spl_k5, 1.3, 2) + sin(1.3)) < 1.9e-8
    @test abs(Maths.derivative(spl_k5, 1.3, 3) + cos(1.3)) < 1.6e-7
    @test abs(Maths.derivative(spl_k5, 1.3, 4) - sin(1.3)) < 2e-4
    @test abs(Maths.derivative(spl_k5, 1.3, 5) - cos(1.3)) < 2e-3
end

@testset "randgauss" begin
    import Statistics: std, mean
    seed!(1234)
    x = Maths.randgauss(1, 0.5, 1000000)
    @test isapprox(std(x), 0.5, rtol=1e-3)
    @test isapprox(mean(x), 1, rtol=1e-3)
    seed!(1234)
    x = Maths.randgauss(10, 0.1, 1000000)
    @test isapprox(std(x), 0.1, rtol=1e-3)
    @test isapprox(mean(x), 10, rtol=1e-3)
    seed!(1234)
    x = Maths.randgauss(-1, 0.5, 1000000)
    @test isapprox(std(x), 0.5, rtol=1e-3)
    @test isapprox(mean(x), -1, rtol=1e-3)
    seed!(1234)
    x = Maths.randgauss(1, 0.5, (1000, 1000))
    @test isapprox(std(x), 0.5, rtol=1e-3)
    @test isapprox(mean(x), 1, rtol=1e-3)
end

@testset "linterp" begin
    x = range(0.0, 2π, length=100)
    y = sin.(x)
    l = Maths.LinTerp(x, y)
    @test all(abs.(l.(x) .- y) .< 3e-16)
    x2 = range(0.0, 2π, length=300)
    @test isapprox(l.(x2), sin.(x2), rtol=4e-4)
    @test Maths.linterp(0.5, 0.0, 0.0, 1.0, 1.0) ≈ 0.5
    @test Maths.linterp(1.5, 0.0, 0.0, 1.0, 1.0) ≈ 1.5
    @test Maths.linterp(-0.5, 0.0, 0.0, 1.0, 1.0) ≈ -0.5
end

@testset "gaussnorm" begin
    for power = 2:2:10
        @test Maths.gaussnorm(1, power=power) ≈ hquadrature(
            x -> Maths.gauss(x, 1; power=power), -10, 10)[1]
        @test Maths.gaussnorm(fwhm=1, power=power) ≈ hquadrature(
            x -> Maths.gauss(x; fwhm=1, power=power), -10, 10)[1]
    end
end

@testset "findpeaks" begin
    grid = Grid.RealGrid(1.0, 800e-9, (160e-9, 3000e-9), 10e-12)
    dt = grid.t[2] - grid.t[1]
    x = Array{Float64}(undef, length(grid.t))
    FT = FFTW.plan_rfft(x, 1)
    Eω = FT * x
    fill!(Eω, 0.0)
    positions = (-367e-15, -10e-15, 100e-15, 589e-15)
    widths = (30e-15, 3e-15, 100e-15, 200e-15)
    powers = (1e3, 1e4, 1e2, 2e3)
    for i in 1:length(positions)
        field = Fields.GaussField(λ0=800e-9, τfwhm=widths[i], power=powers[i],
                                  ϕ=[2π*i/length(widths), positions[i]])
        Eω .+= field(grid, FT)
    end
    Et = FT \ Eω
    It = Fields.It(Et, grid)
    pks = Maths.findpeaks(grid.t, It, threshold=10.0, filterfw=false)
    for i in 1:length(positions)
        @test isapprox(pks[i].position, positions[i], atol=dt*1.2)
        @test isapprox(pks[i].fw, widths[i], rtol=1e-7, atol=dt*1.2)
        @test isapprox(pks[i].peak, powers[i], rtol=1e-2)
    end
end


text:
# ---------------------- #
# Piecewise linear Basis #
# ---------------------- #

function Basis(::Lin, breaks::Vector, evennum::Int=0)
    n = length(breaks)  # 28
    !(issorted(breaks)) && error("breaks should be increasing")

    if evennum != 0
        if length(breaks) == 2
            breaks = linspace(breaks[1], breaks[2], evennum)
        else
            if length(breaks) < 2
                error("breaks must have at least 2 elements")
            end

            if any(abs(diff(diff(breaks))) .> 5e-15*mean(abs(breaks)))
                error("Breaks not evenly spaced")
            end
            evennum = length(breaks)
        end
    end
    n = length(breaks)
    a = breaks[1]
    b = breaks[end]
    Basis(Lin(), n, a, b, LinParams(breaks, evennum))
end

# define methods for LinParams type
Basis(p::LinParams) = Basis(Lin(), p.breaks, p.evennum)

nodes(p::LinParams) = p.breaks

function derivative_op(p::LinParams, order::Int=1)
    breaks, evennum = p.breaks, p.evennum

    newbreaks = breaks
    n = length(breaks)
    D = Array(SparseMatrixCSC{Float64, Int}, abs(order))

    for i in 1:order
        d = 1./diff(newbreaks)
        d = sparse([1:n-1; 1:n-1], [1:n-1; 2:n], [-d; d], n-1, n)
        if i > 1
            D[i] = d*D[i-1]
        else
            D[1] = d
        end
        newbreaks = (newbreaks[1:end-1]+newbreaks[2:end])/2
        n = n-1
    end

    for i in -1:-1:order
        newbreaks=[dot([3, -1], newbreaks[1:2]);
                   (newbreaks[1:end-1]+newbreaks[2:end]);
                   dot([-1, 3], newbreaks[end-1:end])]/2
        d = diff(newbreaks)'
        n = n+1
        d = tril(repmat(d, n, 1), -1)
        if i<-1
            D[-i] = d*D[-i-1]
        else
            D[1] = d
        end
        #adjustment to make value at original left endpoint equal 0
        if evennum > 0
            temp = evalbase(LinParams(newbreaks, length(newbreaks)),
                            breaks[1], 0)*D[-i]
        else
            temp = evalbase(LinParams(newbreaks, 0), breaks[1], 0)*D[-i]
        end
        D[-i] = D[-i]-repmat(temp, length(newbreaks), 1)
    end

    params = LinParams(newbreaks, evennum > 0 ? evennum : 0)

    return D, params

end

function _prep_evalbase(p::LinParams, x::Real)
    n = length(p.breaks)
    if p.evennum != 0
        _ind = fix((x-p.breaks[1])*((n-1)/(p.breaks[end]-p.breaks[1])))
        ind = clamp(_ind+1, 1, n-1)
    else
        ind = lookup(p.breaks, x, 3)
    end
    1, n, [ind]
end

function _prep_evalbase(p::LinParams, x::AbstractArray)
    m = size(x, 1)
    n = length(p.breaks)

    # Determine the maximum index of
    #   the breakpoints that are less than or equal to x,
    #   (if x=b use the index of the next to last breakpoint).
    if p.evennum != 0
        ind = fix((x-p.breaks[1]).*((n-1)./(p.breaks[end]-p.breaks[1]))) + 1
        clamp!(ind, 1, n-1)
    else
        ind = lookup(p.breaks, x, 3)
    end

    return m, n, ind
end

function evalbase(::Type{SparseMatrixCSC}, p::LinParams,
                  x::Union{Real,AbstractArray}=nodes(p), order::Int=0)
    # 46-49
    if order != 0
        D, params = derivative_op(p, order)
        B = evalbase(params, x, 0) * D[end]
        return B
    end

    m, n, ind = _prep_evalbase(p, x)
    z = similar(x)
    for i in 1:length(x)
        z[i] = (x[i]-p.breaks[ind[i]])/(p.breaks[ind[i]+1]-p.breaks[ind[i]])
    end

    return sparse(vcat(1:m, 1:m), vcat(ind, ind+1), vcat(1-z, z), m, n)
end

function evalbase(::Type{SplineSparse}, p::LinParams,
                  x::Union{Real,AbstractArray}=nodes(p), order::Int=0)
    # 46-49
    if order != 0
        error("derivatives un-supported right now")
    end

    m, n, ind = _prep_evalbase(p, x)

    z = Array(eltype(x), 2*length(x))
    for i in 1:length(x)
        ix = 2i
        z[ix] = (x[i]-p.breaks[ind[i]])/(p.breaks[ind[i]+1]-p.breaks[ind[i]])
        z[ix-1] = 1 - z[ix]
    end

    return SplineSparse(1, 2, n, z, ind)
end

evalbase(p::LinParams, x::Union{Real,AbstractArray}=nodes(p), order::Int=0) =
    evalbase(SparseMatrixCSC, p, x, order)

function evalbase(p::LinParams, x, order::AbstractArray{Int})
    out = Array(SparseMatrixCSC{Float64,Int}, size(order))

    for I in eachindex(order)
        out[I] = evalbase(p, x, order[I])
    end
    return out
end


text:
@everywhere begin
    using DrWatson
    import JSON
    using Logging
    using HDF5
    using PWS
end

f = ARGS[1]
dict = JSON.parsefile(projectdir("_research", "tmp", f))
@info "Read file" file = projectdir("_research", "tmp", f)

run_name = dict["run_name"]
duration = dict["duration"]
num_responses = dict["num_responses"]

mean_L = dict["mean_L"]
num_receptors = dict["num_receptors"]
Y_tot = dict["Y_tot"]

LR_timescale = dict["LR_timescale"]
Y_timescale = dict["Y_timescale"]

dtimes = collect(0.0:0.04:duration)

system_fn = () -> PWS.chemotaxis_system(
    mean_L = mean_L,
    num_receptors = num_receptors,
    Y_tot = Y_tot,
    LR_timescale = LR_timescale,
    Y_timescale = Y_timescale,
    dtimes = dtimes,
    Y_ratio = 1/6,
    LR_ratio = 0.5
)

algorithm = SMCEstimate(dict["smc_samples"])

mi = PWS.run_parallel(system_fn, algorithm, num_responses)
result = Dict("Samples" => mi, "DiscreteTimes" => dtimes)

function DrWatson._wsave(filename, result::Dict)
    h5open(filename, "w") do file
        PWS.write_hdf5!(file, result)
    end
end

save_dict = Dict("Duration" => duration, "TauLR" => LR_timescale, "TauY" => Y_timescale)

filename = savename(save_dict, "hdf5")
local_path = datadir("chemotaxis", run_name, filename)
tagsave(local_path, merge(dict, result), storepatch=false)
@info "Saved to" filename

# upload to SUN storage

# include("smbclient.jl")
# sun_path = joinpath(sun_home, "data", "chemotaxis", run_name, filename)
# mkpath(dirname(sun_path))
# cp(local_path, sun_path)


text:
#=
This file is auto-generated. Do not edit.
=#

mutable struct HydroDispatch <: HydroGen
    name::String
    available::Bool
    bus::Bus
    activepower::Float64
    reactivepower::Float64
    tech::TechHydro
    op_cost::TwoPartCost
    _forecasts::InfrastructureSystems.Forecasts
    internal::InfrastructureSystemsInternal
end

function HydroDispatch(name, available, bus, activepower, reactivepower, tech, op_cost, _forecasts=InfrastructureSystems.Forecasts(), )
    HydroDispatch(name, available, bus, activepower, reactivepower, tech, op_cost, _forecasts, InfrastructureSystemsInternal())
end

function HydroDispatch(; name, available, bus, activepower, reactivepower, tech, op_cost, _forecasts=InfrastructureSystems.Forecasts(), )
    HydroDispatch(name, available, bus, activepower, reactivepower, tech, op_cost, _forecasts, )
end

# Constructor for demo purposes; non-functional.

function HydroDispatch(::Nothing)
    HydroDispatch(;
        name="init",
        available=false,
        bus=Bus(nothing),
        activepower=0.0,
        reactivepower=0.0,
        tech=TechHydro(nothing),
        op_cost=TwoPartCost(nothing),
        _forecasts=InfrastructureSystems.Forecasts(),
    )
end

"""Get HydroDispatch name."""
get_name(value::HydroDispatch) = value.name
"""Get HydroDispatch available."""
get_available(value::HydroDispatch) = value.available
"""Get HydroDispatch bus."""
get_bus(value::HydroDispatch) = value.bus
"""Get HydroDispatch activepower."""
get_activepower(value::HydroDispatch) = value.activepower
"""Get HydroDispatch reactivepower."""
get_reactivepower(value::HydroDispatch) = value.reactivepower
"""Get HydroDispatch tech."""
get_tech(value::HydroDispatch) = value.tech
"""Get HydroDispatch op_cost."""
get_op_cost(value::HydroDispatch) = value.op_cost
"""Get HydroDispatch _forecasts."""
get__forecasts(value::HydroDispatch) = value._forecasts
"""Get HydroDispatch internal."""
get_internal(value::HydroDispatch) = value.internal


text:
include(find_in_path("PatternDispatch.jl"))

module TestOps
using PatternDispatch, PatternDispatch.Graph

unordered(x,y) = !(x >= y || y >= x)

@assert (@qpat x) == (@qpat x)
@assert (@qpat x) == (@qpat y)
@assert unordered((@qpat x::Int), (@qpat x::String))
@assert (@qpat x) > (@qpat x::Int)
@assert (@qpat x::Int::Real) == (@qpat x::Int)
@assert (@qpat x) > (@qpat x::Int)

@assert (@qpat x) > (@qpat (x,y))
@assert unordered((@qpat x::Int), (@qpat (x,y)))
@assert (@qpat (x,y)::Number) == nullpat
@assert unordered((@qpat (x::Int,y)), (@qpat (x,y::Int)))
@assert (@qpat (x::Number,y)) > (@qpat (x::Int,y::Int))

@assert (@qpat 1::Float64) == nullpat
@assert (@qpat x::Any) == (@qpat x)
@assert (@qpat (x::None,z)) == nullpat
@assert (@qpat (x,y)::Tuple) == (@qpat (x,y))

end


text:
@testset "debug" begin

# not really a test, since it's not verified.  but if you run this you
# should see appropriate output to stdout.

@test ParserCombinator.truncate("1234567890", 10) == "1234567890"
@test ParserCombinator.truncate("1234567890",  9) == "1234...90"
@test ParserCombinator.truncate("1234567890",  8) == "123...90"
@test ParserCombinator.truncate("1234567890",  7) == "123...0"
@test ParserCombinator.truncate("1234567890",  6) == "12...0"
@test ParserCombinator.truncate("1234567890",  5) == "1...0"

println("one level")
parse_dbg("ab", Trace(Dot()))

println("multiple")
parse_dbg("ab", Equal("a") + Trace(Dot()[0:end]) + Equal("b"))

grammar = p"\d+" + Eos()
debug, task = make(Debug, "123abc", grammar; delegate=NoCache)
@test_throws ParserException once(task)
@test debug.max_iter == 4

end