text:
struct RejectionExact <: AbstractRejectionExact end

function solve(problem::PDMPProblem, Flow::Function; verbose::Bool = false, save_rejected = false, ind_save_d = -1:1, ind_save_c = -1:1, n_jumps = Inf64, save_positions = (false, true), save_rate = false, finalizer = finalize_dummy)
	verbose && println("#"^30)
	verbose && printstyled(color=:red,"--> Start Rejection method\n")

	# initialise the problem. If I call twice this function, it should give the same result...
	init!(problem)

	# we declare the characteristics for convenience
	caract = problem.caract
	ratecache = caract.ratecache
	simjptimes = problem.simjptimes

	ti, tf = problem.tspan
	# it is faster to pre-allocate arrays and fill it at run time
	n_jumps  += 0 # to hold initial vector
	n_reject  = 0  # to hold the number of rejects
	nsteps  = 1
	npoints = 2 # number of points for ODE integration

	xc0 = caract.xc
	xd0 = caract.xd

	# Set up initial variables
	t = ti

	X0 = copy(xc0)
	Xd = copy(xd0)
	res_ode = zeros(2, length(X0))

	X0, _, Xd, _, xc_hist, xd_hist, res_ode, ind_save_d, ind_save_c = allocate_arrays(ti, xc0, xd0, n_jumps; rejection = true)

	tp = [ti, tf]		   # vector to hold the time interval over which to integrate the flow

	#variables for rejection algorithm
	reject = true
	lambda_star = 0.0 # this is the bound for the rejection method
	ppf = caract.R(ratecache.rate, X0, Xd, caract.parms, t, true)

	δt = simjptimes.tstop_extended

	while (t < tf) && (nsteps < n_jumps)
		verbose && println("--> step : ",nsteps," / ", n_jumps)
		reject = true
		while reject && (nsteps < n_jumps)
			tp .= [t, min(tf, t + δt / ppf[2]) ] #mettre un lambda_star?
			Flow(res_ode, X0, Xd, tp)

			@inbounds for ii in eachindex(X0)
				X0[ii] = res_ode[end, ii]
			end
			verbose && println("----> δt = ", δt, ", t∈", tp, ", dt = ", tp[2]-tp[1], ", xc = ", X0)

			t = tp[end]
			ppf = caract.R(ratecache.rate, X0, Xd, caract.parms, t, true)
			@assert ppf[1] <= ppf[2] "(Rejection algorithm) Your bound on the total rate is wrong, $ppf"
			if t == tf
				reject = false
			else
				reject = rand() < 1 - ppf[1] / ppf[2]
			end
			δt = -log(rand())
			if reject
				n_reject += 1
			end
		end

		# there is a jump!
		ppf = caract.R(ratecache.rate, X0, Xd, caract.parms, t, false)

		if (t < tf)
			verbose && println("----> Jump!, ratio = ", ppf[1] / ppf[2], ", xd = ", Xd)
			# make a jump
			ev = pfsample(ratecache.rate)

			# we perform the jump
			affect!(caract.pdmpjump, ev, X0, Xd, caract.parms, t)
		end

		nsteps += 1
		pushTime!(problem, t)
		push!(xc_hist, X0[ind_save_c])
		push!(xd_hist, Xd[ind_save_d])
		save_rate && push!(problem.rate_hist, sum(ratecache.rate))
		finalizer(ratecache.rate, caract.xc, caract.xd, caract.parms, t)
	end
	if verbose println("--> Done") end
	if verbose println("--> xd = ",xd_hist[:,1:nsteps]) end
	return PDMPResult(problem.time, xc_hist, xd_hist, problem.rate_hist, save_positions, nsteps, n_reject)
end

function solve(problem::PDMPProblem, algo::Rejection{Tode}; reltol = 1e-7, abstol = 1e-9, kwargs...) where {Tode <: Symbol}
	ode = algo.ode
	@assert ode in [:cvode, :lsoda, :adams, :bdf]

	caract = problem.caract

	# define the ODE flow
	if ode == :cvode || ode == :bdf
		Flow0 = (X0_,Xd,tp_) -> Sundials.cvode(  (tt,x,xdot) -> caract.F(xdot,x,Xd,caract.parms,tt), X0_, tp_, abstol = abstol, reltol = reltol, integrator = :BDF)
	elseif	ode == :adams
		Flow0 = (X0_,Xd,tp_) -> Sundials.cvode(  (tt,x,xdot) -> caract.F(xdot,x,Xd,caract.parms,tt), X0_, tp_, abstol = abstol, reltol = reltol, integrator = :Adams)
	elseif ode == :lsoda
		Flow0 = (X0_,Xd,tp_) -> LSODA.lsoda((tt,x,xdot,data) -> caract.F(xdot,x,Xd,caract.parms,tt), X0_, tp_, abstol = abstol, reltol = reltol)
	end

	Flow = (out,X0_,Xd,tp_) -> (out .= Flow0(X0_,Xd,tp_))

	return solve(problem, Flow; kwargs...)
end

function solve(problem::PDMPProblem, algo::Talgo; kwargs...) where {Talgo <: AbstractRejectionExact}
	Flow = (res_ode, X0, Xd, tp) -> problem.caract.F(res_ode, X0, Xd, problem.caract.parms, tp)
	solve(problem, Flow; kwargs...)
end


text:
using BinaryProvider
include("compile.jl")

const verbose = ("--verbose" in ARGS)
const prefix  = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["libextraf"], :libextraf),
]

src_url  = "https://mkatase.github.io/github-hosted/data/libmakef.tar.gz"
src_hash = "d33fef76bf8026112b7d441a4038225d679f305eb23fe2dfb5b249ce270eafbe"

src_path = joinpath(prefix, "downloads", "libmakef.tar.gz")

libname  = "libextraf"

if !isfile(src_path) || !verify(src_path, src_hash; verbose=verbose)
    println("Go to compile")
    compile(libname, src_url, src_hash, prefix=prefix, verbose=verbose)
end

write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)
println("end of build.jl script")

ext: jl


text:
using YAML, DataFrames, CSV, Plots
using Statistics
using Streamfall

HERE = @__DIR__
DATA_PATH = joinpath(HERE, "../../test/data/hymod/")

# Load and generate stream network
network = YAML.load_file(joinpath(DATA_PATH, "hymod_network.yml"))
sn = create_network("HyMod Network", network)

# Load climate data
date_format = "YYYY-mm-dd"
obs_data = CSV.File(joinpath(DATA_PATH, "leaf_river_data.csv"),
                    comment="#",
                    dateformat=date_format) |> DataFrame

hist_streamflow = obs_data[:, "leaf_river_outflow"]
climate_data = obs_data[:, ["Date", "leaf_river_P", "leaf_river_ET"]]
climate = Climate(climate_data, "_P", "_ET")

# This will set node parameters to the optimal values found
metric = (obs, sim) -> 1.0 - Streamfall.NNSE(obs, sim)
calibrate!(sn, climate, hist_streamflow; metric=metric, MaxTime=90.0)

# Save calibrated network spec to file
Streamfall.save_network_spec(sn, "hymod_example_calibrated.yml")

# Run the model as a single system
run_basin!(sn, climate)

# Get node
nid, node = sn["leaf_river"]

@info "Mean flow (ft^3/s)" mean(node.outflow)

# Get performance and plot
obs = hist_streamflow[:, "leaf_river_outflow"]
@info "RMSE" Streamfall.RMSE(obs, node.outflow)
@info "NSE" Streamfall.NSE(obs, node.outflow)

plot(obs)
plot!(node.outflow)


text:
module JLBOX_MODULETest
include("$(pwd())/test/helper.jl")
reload("$(pwd())/src/JLBOX_MODULE.jl")
using JLBOX_MODULE

facts("JLBOX_MODULE") do
    @fact isempty(lintfile("$(pwd())/src/JLBOX_MODULE.jl", returnMsgs=true)) => true
    @fact isempty(lintfile("$(pwd())/test/JLBOX_MODULE_test.jl", returnMsgs=true)) => true
end

end # module JLBOX_MODULETest

max_stars_repo_name: jbnicolai/jlbox

max_issues_repo_name: jbnicolai/jlbox

max_forks_repo_name: jbnicolai/jlbox


text:
struct MaximalCoupling{U1<:Distribution, U2<:Distribution}
    p::U1
    q::U2
    function MaximalCoupling{U1, U2}(p::U1, q::U2) where {U1, U2}
        length(p) == length(q) ||
            throw(DimensionMismatch("Coupled distributions of different dimension."))
        eltype(p) == eltype(q) ||
            throw(ArgumentError("Coupled distributions of different types."))
        new(p, q)
    end
end
MaximalCoupling(p::U1, q::U2) where {U1, U2} = MaximalCoupling{U1, U2}(p, q)

length(coup::MaximalCoupling) = length(coup.p)
eltype(coup::MaximalCoupling) = Union{eltype(coup.p), eltype(coup.q)}


text:
## functions for generic programming within SymPy
## if a generic method should apply to a sympy variable, it would
## be defined here.

## Iterator for Sym
import Base.iterate
iterate(x::Sym) = (x.__pyobject__, 0)
iterate(x::Sym, state) = nothing

Base.isless(a::Sym, b::Sym) = sympy.Lt(a,b) == True
Base.isless(a::Sym, b::Number) = sympy.Lt(promote(a,b)...) == True
Base.isless(a::Number, b::Sym) = sympy.Lt(promote(a,b)...) == True

Base.isequal(a::Sym, b::Sym) = Eq(a,b) == True
Base.isequal(a::Sym, b::Number) = Eq(promote(a,b)...) == True
Base.isequal(a::Number, b::Sym) = Eq(promote(a,b)...) == True

# Floating point bits
Base.eps(::Type{Sym}) = zero(Sym)
Base.eps(::Sym) = zero(Sym)
Base.signbit(x::Sym) = x < 0
Base.copysign(x::Sym, y::Number) = abs(x) * sign(y)
Base.flipsign(x::Sym, y) = signbit(y) ? -x : x
Base.typemax(::Type{Sym}) = oo
Base.typemin(::Type{Sym}) = -oo

Base.fld(x::SymbolicObject, y) = floor(x/y)
Base.cld(x::SymbolicObject, y) = ceil(x/y)
Base.mod(x::SymbolicObject, args...)= Mod(x, args...)
#Base.mod1
#Base.mod2pi
#Base.fldmod

function Base.round(x::Sym; kwargs...)
    length(free_symbols(x)) > 0 && throw(ArgumentError("can't round a symbolic expression"))
    round(N(x); kwargs...)
end

function Base.trunc(x::Sym; kwargs...)
    length(free_symbols(x)) > 0 && throw(ArgumentError("can't truncate a symbolic expression"))
    trunc(N(x); kwargs...)
end

# check on type of number
# these are boolean: true/false; not tru/false/nothing,as in SymPy
Base.isfinite(x::Sym) = !is_(:infinite, x)
Base.isinf(x::Sym) = is_(:infinite, x)
Base.isnan(x::Sym) = x == sympy.nan
Base.isinteger(x::Sym) = is_(:integer, x)
Base.iseven(x::Sym) = is_(:even, x)
Base.isreal(x::Sym) = is_(:real, x)
Base.isodd(x::Sym) = is_(:odd, x)

## zero and one (zeros?)
Base.zero(x::Sym) = Sym(0)
Base.zero(::Type{Sym}) = Sym(0)

Base.one(x::Sym) = Sym(1)
Base.one(::Type{Sym}) = Sym(1)

## float, complex, real, image,
Base.float(x::Sym) = _float(N(x))
_float(x::Sym) = throw(ArgumentError("variable must have no free symbols"))
_float(x) = float(x)
Base.Float64(x::Sym) = _Float64(N(x))
_Float64(x::Sym) = throw(ArgumentError("variable must have no free symbols"))
_Float64(x) = Float64(x)

Base.complex(::Type{Sym}) = Sym
Base.complex(r::Sym) = real(r) + imag(r) * im
function Base.complex(r::Sym, i)
    isreal(r) || throw(ArgumentError("r and i must not be complex"))
    isreal(i) || throw(ArgumentError("r and i must not be complex"))
    N(r) + N(i) * im
end
Base.complex(xs::AbstractArray{Sym}) = complex.(xs) # why is this in base?

Base.conj(x::SymbolicObject) = x.conjugate()
function Base.transpose(f::Sym)::Sym
    if pycall_hasproperty(PyObject(f), :transpose)
        f.transpose()
    else
        f
    end
end

Base.real(::Type{Sym}) = Sym
Base.real(x::Sym) = _real(N(x))
_real(x::Sym) = sympy.:re(x)
_real(x) = real(x)

Base.imag(x::Sym) = _imag(N(x))
_imag(x::Sym) = sympy.im(x)
_imag(x) = imag(x)

# sympy.div for poly division
Base.div(x::Sym, y::Union{Sym,Number}) = convert(Sym, sympy.floor(x/convert(Sym,y)))
Base.rem(x::Sym, y::Union{Sym,Number}) = x-Sym(y)*Sym(sympy.floor.(x/y))

Base.denominator(x::SymbolicObject) = denom(x)
Base.numerator(x::SymbolicObject) = numer(x)


text:
abstract time_each_annual
abstract layers_each_time

type Layer_1st <: layers_each_time
  start :: Array{Int64,1}
  count :: Array{Int64,1}
end

SurfaceLayer = Layer_1st([1,1,0,1], [0,0,1,1])

function select_layer( dim , L :: layers_each_time)
  error("No function for a pure abstract layers_each_time")
end

function select_layer ( dim , L :: Layer_1st )
  idices = find(L.start .<=0 )
  start = copy(L.start)
  for jdx in idices
    start[jdx] = start[jdx]+dim[jdx]
  end

  idices = find(L.count .<=0 )
  count = copy(L.count)
  for jdx in idices
    count[jdx] = count[jdx]+dim[jdx]
  end
  return (start, count)
end
export select_layer, SurfaceLayer, Layer_1st, layers_each_time

#= abstract time_each_annual =#

function select_months( y , T :: time_each_annual)
  error("No function for a pure abstract time_each_annual")
end

type Annual_Months <: time_each_annual
  seasonidx :: Array{Int64,2}
end

function select_months( y  , T :: Annual_Months)
  y1 = y[:,:,:,T.seasonidx[:]]
  y2 = reshape(y1, (
                    size(y1,1),size(y1,2),size(y1,3),
                    size(T.seasonidx,1),size(T.seasonidx,2)))
  return(squeeze(mean(y2,4),4))
end
Seasons_each_Year = Annual_Months(reshape([1:12],3,4))
Year_each_Year = Annual_Months(reshape([1:12],12,1))
Months_each_Year = Annual_Months(reshape([1:12],1,12))

export time_each_annual, select_months
export Annual_Months
export Seasons_each_Year, Year_each_Year, Months_each_Year


text:
"""
    Calculates ECI Position of SV

    $SIGNATURES
    ´sat_state´: satellite state, combining decoded data, code- and carrierphase

    This function calculates the position of the SV in ECI coorinates.
    The implementation follows IS-GPS-200K
"""
function sat_position_ECI(sat_state::SatelliteState)
    is_sat_healthy_and_decodable(sat_state.decoder_state) || throw(BadData("SV not decoded properly"))

    t = calc_corrected_time(sat_state)
    tₖ = t - sat_state.decoder_state.data.t_oe
    tₖ = check_crossover(tₖ)

    # Excentric and true anomaly
    tₛᵥ = calc_uncorrected_time(sat_state)
    dtₛᵥ = code_phase_offset(sat_state.decoder_state, tₛᵥ)
    Eₖ = calc_eccentric_anomaly(sat_state.decoder_state, t, dtₛᵥ)

    e  = sat_state.decoder_state.data.e
    vₖ =  2 * atan(sqrt((1 + e) / (1 - e)) * tan(Eₖ / 2))
    # Argument of latitude
    Φₖ = vₖ + sat_state.decoder_state.data.ω

    # Second harmonic perturbations
    δuₖ = sat_state.decoder_state.data.C_us * sin(2 * Φₖ) + sat_state.decoder_state.data.C_uc * cos(2 * Φₖ);
    δrₖ = sat_state.decoder_state.data.C_rs * sin(2 * Φₖ) + sat_state.decoder_state.data.C_rc * cos(2 * Φₖ);
    δiₖ = sat_state.decoder_state.data.C_is * sin(2 * Φₖ) + sat_state.decoder_state.data.C_ic * cos(2 * Φₖ);

    # Corrected orbit parameters
    uₖ = Φₖ + δuₖ;
    rₖ = sat_state.decoder_state.data.sqrt_A^2 * (1 - e * cos(Eₖ)) + δrₖ;
    iₖ = sat_state.decoder_state.data.i_0 + δiₖ + sat_state.decoder_state.data.IDOT * tₖ;

    # Positions in orbital plane
    xₖs = rₖ * cos(uₖ);
    yₖs = rₖ * sin(uₖ);

    # Corrected longitude of ascending node
    Ωₖ = sat_state.decoder_state.data.Ω_0 + sat_state.decoder_state.data.Ω_dot * tₖ;

    # ECI coordinates
    sinΩₖ = sin(Ωₖ);
    cosΩₖ = cos(Ωₖ);
    cosiₖ = cos(iₖ);
    siniₖ = sin(iₖ);

    xₖ = xₖs * cosΩₖ - yₖs * cosiₖ * sinΩₖ;
    yₖ = xₖs * sinΩₖ + yₖs * cosiₖ * cosΩₖ;
    zₖ = yₖs * siniₖ;

    return [xₖ, yₖ, zₖ]
end

"""
    Calculates ECEF Position by converting ECI position

    $SIGNATURES
    ´sat_state´: satellite state, combining decoded data, code- and carrierphase

    This function calculates the position of the SV in ECI coorinates and
    converts this position into ECEF coordinates.
    The implementation follows IS-GPS-200K
"""
function sat_position_ECI_2_ECEF(sat_state::SatelliteState)
    pos_ECI = sat_position_ECI(sat_state)

    tₖ = calc_corrected_time(sat_state)
    tₖ = tₖ - sat_state.decoder_state.data.t_oe
    tₖ = check_crossover(tₖ)
    θ = sat_state.decoder_state.constants.Ω_dot_e * (tₖ + sat_state.decoder_state.data.t_oe)
    cosθ  = cos(θ)
    sinθ  = sin(θ)
    R = [cosθ sinθ; -sinθ cosθ]
    pos_ECI[1:2] = R * pos_ECI[1:2]
    return pos_ECI
end

"""
    Calculates ECEF Position of SV

    $SIGNATURES
    ´sat_state´: satellite state, combining decoded data, code- and carrierphase

    This function calculates the position of the SV in ECEF coorinates.
    The implementation follows IS-GPS-200K Table 20-IV.
"""
function sat_position_ECEF(sat_state::SatelliteState)

    is_sat_healthy_and_decodable(sat_state.decoder_state) || throw(BadData("SV not decoded properly"))
    F = sat_state.decoder_state.constants.F
    e = sat_state.decoder_state.data.e
    μ = sat_state.decoder_state.constants.μ
    Ω_dot_e = sat_state.decoder_state.constants.Ω_dot_e

    A = sat_state.decoder_state.data.sqrt_A^2
    n₀ = sqrt(μ / (A^3))
    n = n₀ + sat_state.decoder_state.data.Δn

    t = calc_corrected_time(sat_state)
    tₖ = t - sat_state.decoder_state.data.t_oe
    tₖ = check_crossover(tₖ)

    # Excentric and true anomaly
    t_sv = calc_uncorrected_time(sat_state)
    dt_sv = code_phase_offset(sat_state.decoder_state, t_sv)
    E = calc_eccentric_anomaly(sat_state.decoder_state, t, dt_sv)

    vₖ = 2 * atan(sqrt((1 + e) / (1 - e)) * tan(E / 2))

    Φₖ = vₖ + sat_state.decoder_state.data.ω

    δuₖ = sat_state.decoder_state.data.C_us * sin(2 * Φₖ) + sat_state.decoder_state.data.C_uc * cos(2 * Φₖ)
    δrₖ = sat_state.decoder_state.data.C_rs * sin(2 * Φₖ) + sat_state.decoder_state.data.C_rc * cos(2 * Φₖ)
    δiₖ = sat_state.decoder_state.data.C_is * sin(2 * Φₖ) + sat_state.decoder_state.data.C_ic * cos(2 * Φₖ)

    uₖ = Φₖ + δuₖ
    rₖ = A * (1 - e * cos(E)) + δrₖ
    iₖ = sat_state.decoder_state.data.i_0 + δiₖ + sat_state.decoder_state.data.IDOT * tₖ

    xₖs = rₖ * cos(uₖ)
    yₖs = rₖ * sin(uₖ)

    Ωₖ = sat_state.decoder_state.data.Ω_0 +
        (sat_state.decoder_state.data.Ω_dot - Ω_dot_e) * tₖ -
        sat_state.decoder_state.data.t_oe * Ω_dot_e

    xₖ = xₖs * cos(Ωₖ) - yₖs * cos(iₖ) * sin(Ωₖ)
    yₖ = xₖs * sin(Ωₖ) + yₖs * cos(iₖ) * cos(Ωₖ)
    zₖ = yₖs * sin(iₖ)

    position = [xₖ, yₖ, zₖ]
    return position
end

"""
    Checks if satellite data is correctly tranmitted
    $SIGNATURES
    ´decoder_state´: Decoder

    Checks if satellite data contains the needed information and
    if errors during decoding occured
"""
function is_sat_healthy_and_decodable(decoder_state::GNSSDecoderState)
    decoder_state.data.svhealth == "000000" ? decoder_state.subframes_decoded[1:3] == [1,1,1] : false
end


text:
using DataInterpolations, Test
using QuadGK
using DataInterpolations: integral

function test_integral(func, tspan, name::String)
  t1 = minimum(tspan)
  t2 = maximum(tspan)
  @testset "$name" begin
    qint, err = quadgk(func, t1, t2)
    aint = integral(func, t1, t2)
    @test isapprox(qint, aint, atol=1e-8)
  end
end

# Linear Interpolation
u = 2.0collect(1:10)
t = 1.0collect(1:10)
A = LinearInterpolation(u,t)

test_integral(A, t, "Linear Interpolation (Vector)")

# Quadratic Interpolation
u = [1.0, 4.0, 9.0, 16.0]
t = [1.0, 2.0, 3.0, 4.0]
A = QuadraticInterpolation(u,t)

test_integral(A, t, "Quadratic Interpolation (Vector)")

# QuadraticSpline Interpolation
u = [0.0, 1.0, 3.0]
t = [-1.0, 0.0, 1.0]

A = QuadraticSpline(u,t)

test_integral(A, t, "Quadratic Interpolation (Vector)")

# CubicSpline Interpolation
u = [0.0, 1.0, 3.0]
t = [-1.0, 0.0, 1.0]

A = CubicSpline(u,t)

test_integral(A, t, "Cubic Spline Interpolation (Vector)")


text:
"""
using Pkg; Pkg.add("LBFGSB"); Pkg.add("FFTW"); Pkg.add("ASE")
import Conda
run(`conda create -n conda_jl python conda`)
ENV["CONDA_JL_HOME"] = "/home/schinavro/anaconda3/envs/conda_jl"

ENV["PYTHON"] = "/home/schinavro/anaconda3/bin/python"
Pkg.build("PyCall")
"""

using NPZ
import MPI
MPI.Init()

comm = MPI.COMM_WORLD
nprc = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)
root = 0

println("I am rank " * string(rank) * " of nprc " * string(nprc))
indict = npzread(ARGS[1])
if rank == root
    println(size(indict["ak"]))
end
# indict = MPI.bcast(indict, root, comm)

t0 = indict["epoch"]
Nk = indict["Nk"]
D = indict["D"]
init = indict["init"]
fin = indict["fin"]

ak = indict["ak"]

dir = (fin - init) / t0
ν = π/t0

x(θ, t::Real) = init + dir .* t + 2 * sum((θ .* sin.(ν*t * [1:1:Nk;])), dims=1) / (2*(Nk + 1))   # (1 x D)
ẋ(θ, t::Real) = dir .+ 2 * ν * sum((θ .*[1:1:Nk;].* cos.(ν * t * [1:1:Nk;])), dims=1) / (2 * (Nk + 1))
ẍ(θ, t::Real) = -2 * ν * ν * sum((θ .* ([1:1:Nk;].^2) .* sin.(ν * t * [1:1:Nk;])), dims=1) / (2 * (Nk + 1))
#ak = (rand(Nk, D) .- 0.5) .* 0.5 .* LinRange(1, 0, Nk) .* LinRange(1, 0, Nk)

A = [-200., -100., -170, 15]
a = [-1, -1, -6.5, 0.7]
b = [0., 0., 11., 0.6]
c = [-10, -10, -6.5, 0.7]
x0 = [1, 0, -0.5, -1]
y0 = [0, 0.5, 1.5, 1]              # (4, )

"""
Kinetic
"""
function T(θ, t)
    𝐯 = ẋ(θ, t)
    return 0.5 * 1 * sum(𝐯 .* 𝐯)
end
"""
Muller Brown potential
"""
function V(θ, t)
    # 1 x D
    𝑥, 𝑦 = x(θ, t)
    # (4, )
    Vk = A .* exp.(a.*(𝑥 .- x0).^2 + b.*(𝑥.-x0).*(𝑦.-y0) + c .* (𝑦.-y0).^2)
    return sum(Vk) / 100
end
"""
Muller Brown gradient
"""
function ∇V(θ, t)
    # 1 x D
    𝑥, 𝑦 = x(θ, t)
    # (4, )
    𝑥_x0 = 𝑥 .- x0
    𝑦_y0 = 𝑦 .- y0
    # (4, )
    Vk = @. A * exp(a * (𝑥_x0)^2 + b*(𝑥_x0)*(𝑦_y0) + c * (𝑦_y0)^2)
    # (1)
    Fx = sum(@. Vk * (2 * a * 𝑥_x0 + b * 𝑦_y0))
    Fy = sum(@. Vk * (b * 𝑥_x0 + 2 * c * 𝑦_y0))
    # (2, )
    return [Fx, Fy] / 100
end
"""
Muller Brown hessian
"""
function ∇²V(θ, t)
    𝑥, 𝑦 = x(θ, t)            # 1 x D
    # (4, )
    𝑥_x0 = 𝑥 .- x0
    𝑦_y0 = 𝑦 .- y0
    # (4, )
    Vk = @. A * exp(a * (𝑥_x0)^2 + b*(𝑥_x0)*(𝑦_y0) + c * (𝑦_y0)^2)
    dx = @. (2 * a * 𝑥_x0 + b * 𝑦_y0)
    dy = @. (b * 𝑥_x0 + 2 * c * 𝑦_y0)
    # (4, ) -> (1)
    Hxx = sum(@. Vk * (2 * a + dx * dx))
    Hxy = sum(@. Vk * (b + dx * dy))
    Hyy = sum(@. Vk * (2 * c + dy * dy))
    return [Hxx Hxy; Hxy Hyy] / 100
end
"""
Lagrangian
"""
function ADMD(θ, t)
    _T = T(θ, t)
    _V = V(θ, t)
    return _T - _V + muE * (_T+_V - Et) * (_T+_V - Et)
end

"""
Feynman-Kac
"""
muE = 200
Et = -0.47
β = 3
γ = 8 * 3
𝒟 = 1 / (β*γ)

function OnsagerMachlup(θ, t)
    _T = T(θ, t)
    _V = V(θ, t)

    # F = ∇V(θ, t)
    # ℍ = ∇²V(θ, t)
    # 𝐯 = vec(ẋ(θ, t))
    # ℍ𝐯 = ℍ * 𝐯
    # return β^2 * 𝒟 / 4 * (F'F) - 𝒟 * β / 2 * (ℍ𝐯'ℍ𝐯) + muE * (_T+_V - Et) * (_T+_V - Et)
    F = ∇V(θ, t)
    𝐯 = vec(ẋ(θ, t))
    Fv = F * 𝐯
    return β^2 * 𝒟 / 4 * (F'F) - 𝒟 * β / 2 * sqrt(Fv'Fv) + muE * (_T+_V - Et) * (_T+_V - Et)
end

function OnsagerMachlup1(θ, t)
    dV = ∇V(θ, t)
    𝐚 = vec(ẍ(θ, t))
    # ∇V𝐚 = dV .* 𝐚
    #- 2 * (∇V𝐚'∇V𝐚)
    ℍ = ∇²V(θ, t)
    𝐯 = vec(ẋ(θ, t))
    ℍ𝐯 = ℍ * 𝐯
    # 2 * (ℍ𝐯'ℍ𝐯)
    return γ*(𝐚'𝐚) + (dV'dV)/2γ - (ℍ𝐯'𝐚)
end

function OnsagerMachlup4(θ, t)
    dV = ∇V(θ, t)
    𝐚 = vec(ẍ(θ, t))
    # ∇V𝐚 = dV .* 𝐚
    #- 2 * (∇V𝐚'∇V𝐚)
    ℍ = ∇²V(θ, t)
    𝐯 = vec(ẋ(θ, t))
    ℍ𝐯 = ℍ * 𝐯
    # 2 * (ℍ𝐯'ℍ𝐯)
    return (dV'dV) - (ℍ𝐯'ℍ𝐯)
end

function OnsagerMachlup2(θ, t)
    dV = ∇V(θ, t)
    # 𝐯 = vec(ẋ(θ, t))
    𝐚 = vec(ẍ(θ, t))
    ∇V𝐚 = dV .* 𝐚
    # ℍ = ∇²V(θ, t)
    # ℍ𝐯 = ℍ * 𝐚
    # 2 * (ℍ𝐯'ℍ𝐯)
    return (𝐚'𝐚) + (dV'dV) - 2 * (∇V𝐚'∇V𝐚)
end

function OnsagerMachlup2(θ, t)
    dV = ∇V(θ, t)
    dV_ = ∇V(θ, t + 0.001)
    # 𝐯 = vec(ẋ(θ, t))
    𝐚 = vec(ẍ(θ, t))
    ∇V𝐚 = dV .* 𝐚
    # ℍ = ∇²V(θ, t)
    # ℍ𝐯 = ℍ * 𝐚
    # 2 * (ℍ𝐯'ℍ𝐯)
    return (𝐚'𝐚) + (dV'dV) - 2 * (∇V𝐚'∇V𝐚)
end

using QuadGK
"""
Action
"""
function S(θ)
    dt0 = t0 / nprc
    a = rank * dt0
    b = (rank + 1) * dt0
    _S = quadgk(t -> L(θ, t), a, b, rtol=1e-8)[1]
    _S = MPI.Reduce(_S, (S1, S2) -> S1 + S2, root, comm)
    return _S
end
"""
Onsager Machlup action
"""
function S2(θ)
    dt0 = t0 / nprc
    a = rank * dt0
    b = (rank + 1) * dt0
    _S = quadgk(t -> L(θ, t), a, b, rtol=1e-8)[1]
    _S = MPI.Reduce(_S, (S1, S2) -> S1 + S2, root, comm)
    return _S
end

using ForwardDiff

function dS2(g, θ)
    g = vec(ForwardDiff.gradient(S, θ))
    println("dS: ", max(abs.(g)...))
end

function dS(g, θ)
    Nk, D = size(θ)
    n = Nk * D
    # a = div(Nk * D, nprc) * rank + 1
    # b = div(Nk * D, nprc) * (rank + 1)
    # obj.g = MPI.Gather(dSi, root)
    function 𝒮(θ)
        return quadgk(t -> L(θ, t), 0, t0, rtol=1e-8)[1]
    end
    g[1:n] = vec(ForwardDiff.gradient(𝒮, θ))
    println("dS: ", max(abs.(g)...))
end

if indict["L"] == 1
    L = ADMD
elseif indict["L"] == 2
    L = OnsagerMachlup
else
    L = OnsagerMachlup
end

# using Optim
module_mpi = include("./lbfgsb_mpi.jl")
optimizer = module_mpi.L_BFGS_B(Nk * D, 40)
result, resk = optimizer(S, dS, ak, comm, m=20)

# result = optimize(S, ak, BFGS(),
#                       Optim.Options(g_tol = 1e-16, iterations=1000),
#                       autodiff=:forward)
# resk = Optim.minimizer(result)

if rank == root
    # println(map(t -> x(resk, t), LinRange(0, t0, Nk +2)))
    npzwrite("result.npz", Dict("resk" => resk))
end

MPI.Barrier(comm)
MPI.Finalize()

max_stars_repo_name: schinavro/taps

max_issues_repo_name: schinavro/taps

max_forks_repo_name: schinavro/taps


text:
######################################################################
# finfields.jl: Elliptic curves over finite fields
######################################################################

export frobeniustrace, cardinality, frobeniuspolynomial, randXZ, has_montgomery, torsionpoint, torsionXZ, card_over_extension
export base_extend

######################################################################
# Trace, cardinality and Frobenius polynomial
######################################################################

function frobeniustrace(E::EllipticCurve{T}) where T<:FinFieldElem
	#Call PARI
	throw(NotImplementedError("Link to PARI to compute cardinalities of elliptic curves over finite fields"))
end

function cardinality(E::EllipticCurve{T}) where T<:FinFieldElem
	t = frobeniustrace(E)
	q = order(base_ring(E))
	return q + 1 - t
end

function frobeniuspolynomial(E::EllipticCurve{T}, Card) where T<:FinFieldElem
	q = order(base_ring(E))
	t = q + 1 - Card
	A, X = PolynomialRing(Nemo.ZZ, "X")
	return X^2 - t * X + q
end

######################################################################
# Random points on elliptic curves
######################################################################

function rand(E::EllipticCurve{T}) where T<:FinFieldElem
	K = base_ring(E)
	A, Y = PolynomialRing(K, "Y")
	x = rand(K)
	poly = equation(E, x, Y)
	bool, y = any_root(poly)
	while bool == false
		x = rand(K)
		poly = equation(E, x, Y)
		(bool, y) = any_root(poly)
	end
	return Point(x, y, Nemo.one(K), E)
end

function randXZ(E::Montgomery{T}) where T<:FinFieldElem
	P = rand(E)
	return XZPoint(P)
end

######################################################################
# Torsion points
######################################################################

function torsionpoint(E::EllipticCurve{T}, l, Card) where T<:FinFieldElem
	cofactor = div(Card, l)
	@assert cofactor * l == Card
	P = rand(E)
	Q = cofactor * P
	while isinfinity(Q)
		P = rand(E)
		Q = cofactor * P
	end
	isinfinity(l * Q) || throw(ArgumentError("Given curve has no such rational torsion points"))
	return Q
end

function torsionXZ(E::Montgomery{T}, l, Card) where T<:FinFieldElem
	cofactor = div(Card, l)
	@assert cofactor * l == Card
	P = randXZ(E)
	Q = cofactor * P
	while isinfinity(Q)
		P = randXZ(E)
		Q = cofactor * P
	end
	isinfinity(l * Q) || throw(ArgumentError("Given curve has no such torsion rational points"))
	return Q
end

######################################################################
# Base extensions
######################################################################

function base_extend(E::AbstractWeierstrass{T}, K::FinField) where T<:FinFieldElem
	CurveType = curvetype(E)
	K1 = base_ring(E)
	p1 = characteristic(K1)
	p = characteristic(K)
	a1, a2, a3, a4, a6 = a_invariants(E)
	((degree(K1) == 1) & (p1 == p)) || throw(ArgumentError("Invalid field extension"))
	return CurveType(convert(a1, K), convert(a2, K), convert(a3, K), convert(a4, K), convert(a6, K))
end

function base_extend(E::Montgomery{T}, K::FinField) where T<:FinFieldElem
	K1 = base_ring(E)
	p1 = characteristic(K1)
	p = characteristic(K)
	((degree(K1) == 1) & (p1 == p)) || throw(ArgumentError("Invalid field extension"))
	return Montgomery(convert(E.A, K), convert(E.B, K))
end

function card_over_extension(Card::Integer, p::Integer, r::Integer) where T
    t = p + 1 - Card
    W = Base.power_by_squaring([[0, -p] [1, t]], r) * [2, t]
    return p^r + 1 - W[1]
end

######################################################################
# Existence of a Montgomery model
######################################################################

"""
Decide whether a given curve in short Weierstrass form has a Mongomery rational model.
"""
function has_montgomery(E::ShortWeierstrass{T}) where T<:FinFieldElem
	K = base_ring(E)
	A, X = PolynomialRing(K, "X")
	_, _, _, a, b = a_invariants(E)
	poly = X^3 + a * X + b
	r = roots(poly)
	n = length(r)
	for i = 1:n
		alpha, _ = r[i]
		(test, beta) = issquare(3*alpha^2 + E.a)
		test && return (true, Montgomery(3 * alpha // beta, 1 // beta))
	end
	return (false, Montgomery(K(1), K(1)))
end


text:
# __precompile__()

module CaesarLCMTypes

using LCMCore
using StaticArrays: SVector, SMatrix

export
  apriltag_detections_t,
  apriltag_t,
  image_t,
  image_metadata_t,
  point_cloud_t,
  pose_node_t,
  pose_pose_nh_t,
  pose_pose_xyh_nh_t,
  pose_pose_xyh_t,
  prior_zpr_t,
  example_t,
  generic_variable_t,
  generic_factor_t,
  brookstone_supertype_t

include("apriltag_t.jl")
include("apriltag_detections_t.jl")
include("image_metadata_t.jl")
include("image_t.jl")
include("point_cloud_t.jl")
include("pose_node_t.jl")
include("pose_pose_nh_t.jl")
include("pose_pose_xyh_nh_t.jl")
include("pose_pose_xyh_t.jl")
include("prior_zpr_t.jl")
include("example_t.jl")

# generic types for robots of any variety -- can add separate optimizations.
include("generic_variable_t.jl")
include("generic_factor_t.jl")

# supertypes
include("brookstone_supertype_t.jl")

end


text:
""" Simplicial sets in one, two, and three dimensions.

For the time being, this module provides data structures only for [delta
sets](https://en.wikipedia.org/wiki/Delta_set), also known as [semi-simplicial
sets](https://ncatlab.org/nlab/show/semi-simplicial+set). These include the face
maps but not the degeneracy maps of a simplicial set. In the future we may add
support for simplicial sets. The analogy to keep in mind is that graphs are to
semi-simpicial sets as reflexive graphs are to simplicial sets.

Also provided are the fundamental operators on simplicial sets used in nearly
all geometric applications, namely the boundary and coboundary (discrete
exterior derivative) operators. For additional operators, see the
`DiscreteExteriorCalculus` module.
"""
module SimplicialSets
export Simplex, V, E, Tri, SimplexChain, VChain, EChain, TriChain,
  SimplexForm, VForm, EForm, TriForm,
  HasDeltaSet, HasDeltaSet1D, HasDeltaSet2D,
  AbstractDeltaSet1D, DeltaSet1D, OrientedDeltaSet1D, EmbeddedDeltaSet1D,
  AbstractDeltaSet2D, DeltaSet2D, OrientedDeltaSet2D, EmbeddedDeltaSet2D,
  ∂, boundary, coface, d, coboundary, exterior_derivative,
  simplices, nsimplices, point, volume,
  orientation, set_orientation!, orient!, orient_component!,
  src, tgt, nv, ne, vertices, edges, has_vertex, has_edge, edge_vertices,
  add_vertex!, add_vertices!, add_edge!, add_edges!,
  add_sorted_edge!, add_sorted_edges!,
  triangle_edges, triangle_vertices, ntriangles, triangles,
  add_triangle!, glue_triangle!, glue_sorted_triangle!

using LinearAlgebra: det
using SparseArrays
using StaticArrays: @SVector, SVector, SMatrix

using Catlab, Catlab.CategoricalAlgebra, Catlab.Graphs
import Catlab.Graphs: src, tgt, nv, ne, vertices, edges, has_vertex, has_edge,
  add_vertex!, add_vertices!, add_edge!, add_edges!
using ..ArrayUtils

# 0-D simplicial sets
#####################

@present DeltaCategory0D(FreeSchema) begin
  V::Ob
end

""" Abstract type for C-sets that contain a delta set of some dimension.

This dimension could be zero, in which case the delta set consists only of
vertices (0-simplices).
"""
@abstract_acset_type HasDeltaSet

vertices(s::HasDeltaSet) = parts(s, :V)
nv(s::HasDeltaSet) = nparts(s, :V)
nsimplices(::Type{Val{0}}, s::HasDeltaSet) = nv(s)

has_vertex(s::HasDeltaSet, v) = has_part(s, :V, v)
add_vertex!(s::HasDeltaSet; kw...) = add_part!(s, :V; kw...)
add_vertices!(s::HasDeltaSet, n::Int; kw...) = add_parts!(s, :V, n; kw...)

# 1D simplicial sets
####################

@present DeltaCategory1D <: DeltaCategory0D begin
  E::Ob
  (∂v0, ∂v1)::Hom(E, V) # (∂₁(0), ∂₁(1))
end

""" Abstract type for C-sets that contain a one-dimensional delta set.
"""
@abstract_acset_type HasDeltaSet1D <: HasDeltaSet

""" Abstract type for one-dimensional delta sets, aka semi-simplicial sets.
"""
@abstract_acset_type AbstractDeltaSet1D <: HasDeltaSet1D

""" A one-dimensional delta set, aka semi-simplicial set.

Delta sets in 1D are isomorphic to graphs (in the category theorist's sense).
The source and target of an edge can be accessed using the face maps [`∂`](@ref)
(simplicial terminology) or `src` and `tgt` maps (graph-theoretic terminology).
More generally, this type implements the graphs interface in `Catlab.Graphs`.
"""
@acset_type DeltaSet1D(DeltaCategory1D, index=[:∂v0,:∂v1]) <: AbstractDeltaSet1D

edges(s::HasDeltaSet1D) = parts(s, :E)
edges(s::HasDeltaSet1D, src::Int, tgt::Int) =
  (e for e in coface(1,1,s,src) if ∂(1,0,s,e) == tgt)

ne(s::HasDeltaSet1D) = nparts(s, :E)
nsimplices(::Type{Val{1}}, s::HasDeltaSet1D) = ne(s)

has_edge(s::HasDeltaSet1D, e) = has_part(s, :E, e)
has_edge(s::HasDeltaSet1D, src::Int, tgt::Int) =
  has_vertex(s, src) && any(e -> ∂(1,0,s,e) == tgt, coface(1,1,s,src))

src(s::HasDeltaSet1D, args...) = subpart(s, args..., :∂v1)
tgt(s::HasDeltaSet1D, args...) = subpart(s, args..., :∂v0)
face(::Type{Val{(1,0)}}, s::HasDeltaSet1D, args...) = subpart(s, args..., :∂v0)
face(::Type{Val{(1,1)}}, s::HasDeltaSet1D, args...) = subpart(s, args..., :∂v1)

coface(::Type{Val{(1,0)}}, s::HasDeltaSet1D, args...) = incident(s, args..., :∂v0)
coface(::Type{Val{(1,1)}}, s::HasDeltaSet1D, args...) = incident(s, args..., :∂v1)

""" Boundary vertices of an edge.
"""
edge_vertices(s::HasDeltaSet1D, e...) = SVector(∂(1,0,s,e...), ∂(1,1,s,e...))

add_edge!(s::HasDeltaSet1D, src::Int, tgt::Int; kw...) =
  add_part!(s, :E; ∂v1=src, ∂v0=tgt, kw...)

function add_edges!(s::HasDeltaSet1D, srcs::AbstractVector{Int},
                    tgts::AbstractVector{Int}; kw...)
  @assert (n = length(srcs)) == length(tgts)
  add_parts!(s, :E, n; ∂v1=srcs, ∂v0=tgts, kw...)
end

""" Add edge to simplicial set, respecting the order of the vertex IDs.
"""
add_sorted_edge!(s::HasDeltaSet1D, v₀::Int, v₁::Int; kw...) =
  add_edge!(s, min(v₀, v₁), max(v₀, v₁); kw...)

""" Add edges to simplicial set, respecting the order of the vertex IDs.
"""
function add_sorted_edges!(s::HasDeltaSet1D, vs₀::AbstractVector{Int},
                           vs₁::AbstractVector{Int}; kw...)
  add_edges!(s, min.(vs₀, vs₁), max.(vs₀, vs₁); kw...)
end

# 1D oriented simplicial sets
#----------------------------

@present OrientedDeltaSchema1D <: DeltaCategory1D begin
  Orientation::AttrType
  edge_orientation::Attr(E,Orientation)
end

""" A one-dimensional oriented delta set.

Edges are oriented from source to target when `edge_orientation` is
true/positive and from target to source when it is false/negative.
"""
@acset_type OrientedDeltaSet1D(OrientedDeltaSchema1D,
                               index=[:∂v0,:∂v1]) <: AbstractDeltaSet1D

orientation(::Type{Val{1}}, s::HasDeltaSet1D, args...) =
  s[args..., :edge_orientation]
set_orientation!(::Type{Val{1}}, s::HasDeltaSet1D, e, orientation) =
  (s[e, :edge_orientation] = orientation)

function ∂_nz(::Type{Val{1}}, s::HasDeltaSet1D, e::Int)
  (edge_vertices(s, e), sign(1,s,e) * @SVector([1,-1]))
end

function d_nz(::Type{Val{0}}, s::HasDeltaSet1D, v::Int)
  e₀, e₁ = coface(1,0,s,v), coface(1,1,s,v)
  (lazy(vcat, e₀, e₁), lazy(vcat, sign(1,s,e₀), -sign(1,s,e₁)))
end

# 1D embedded simplicial sets
#----------------------------

@present EmbeddedDeltaSchema1D <: OrientedDeltaSchema1D begin
  Point::AttrType
  point::Attr(V, Point)
end

""" A one-dimensional, embedded, oriented delta set.
"""
@acset_type EmbeddedDeltaSet1D(EmbeddedDeltaSchema1D,
                               index=[:∂v0,:∂v1]) <: AbstractDeltaSet1D

""" Point associated with vertex of complex.
"""
point(s::HasDeltaSet, args...) = s[args..., :point]

struct CayleyMengerDet end

volume(::Type{Val{n}}, s::EmbeddedDeltaSet1D, x) where n =
  volume(Val{n}, s, x, CayleyMengerDet())
volume(::Type{Val{1}}, s::HasDeltaSet1D, e::Int, ::CayleyMengerDet) =
  volume(point(s, edge_vertices(s, e)))

# 2D simplicial sets
####################

@present DeltaCategory2D <: DeltaCategory1D begin
  Tri::Ob
  (∂e0, ∂e1, ∂e2)::Hom(Tri,E) # (∂₂(0), ∂₂(1), ∂₂(2))

  # Simplicial identities.
  ∂e1 ⋅ ∂v1 == ∂e2 ⋅ ∂v1 # ∂₂(1) ⋅ ∂₁(1) == ∂₂(2) ⋅ ∂₁(1) == v₀
  ∂e0 ⋅ ∂v1 == ∂e2 ⋅ ∂v0 # ∂₂(0) ⋅ ∂₁(1) == ∂₂(2) ⋅ ∂₁(0) == v₁
  ∂e0 ⋅ ∂v0 == ∂e1 ⋅ ∂v0 # ∂₂(0) ⋅ ∂₁(0) == ∂₂(1) ⋅ ∂₁(0) == v₂
end

""" Abstract type for C-sets containing a 2D delta set.
"""
@abstract_acset_type HasDeltaSet2D <: HasDeltaSet1D

""" Abstract type for 2D delta sets.
"""
@abstract_acset_type AbstractDeltaSet2D <: HasDeltaSet2D

""" A 2D delta set, aka semi-simplicial set.

The triangles in a semi-simpicial set can be interpreted in several ways.
Geometrically, they are triangles (2-simplices) whose three edges are directed
according to a specific pattern, determined by the ordering of the vertices or
equivalently by the simplicial identities. This geometric perspective is present
through the subpart names `∂e0`, `∂e1`, and `∂e2` and through the boundary map
[`∂`](@ref). Alternatively, the triangle can be interpreted as a
higher-dimensional link or morphism, going from two edges in sequence (which
might be called `src2_first` and `src2_last`) to a transitive edge (say `tgt2`).
This is the shape of the binary composition operation in a category.
"""
@acset_type DeltaSet2D(DeltaCategory2D,
                       index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2]) <: AbstractDeltaSet2D

triangles(s::HasDeltaSet2D) = parts(s, :Tri)
ntriangles(s::HasDeltaSet2D) = nparts(s, :Tri)
nsimplices(::Type{Val{2}}, s::HasDeltaSet2D) = ntriangles(s)

face(::Type{Val{(2,0)}}, s::HasDeltaSet2D, args...) = subpart(s, args..., :∂e0)
face(::Type{Val{(2,1)}}, s::HasDeltaSet2D, args...) = subpart(s, args..., :∂e1)
face(::Type{Val{(2,2)}}, s::HasDeltaSet2D, args...) = subpart(s, args..., :∂e2)

coface(::Type{Val{(2,0)}}, s::HasDeltaSet2D, args...) = incident(s, args..., :∂e0)
coface(::Type{Val{(2,1)}}, s::HasDeltaSet2D, args...) = incident(s, args..., :∂e1)
coface(::Type{Val{(2,2)}}, s::HasDeltaSet2D, args...) = incident(s, args..., :∂e2)

""" Boundary edges of a triangle.
"""
function triangle_edges(s::HasDeltaSet2D, t...)
  SVector(∂(2,0,s,t...), ∂(2,1,s,t...), ∂(2,2,s,t...))
end

""" Boundary vertices of a triangle.

This accessor assumes that the simplicial identities hold.
"""
function triangle_vertices(s::HasDeltaSet2D, t...)
  SVector(s[s[t..., :∂e1], :∂v1],
          s[s[t..., :∂e2], :∂v0],
          s[s[t..., :∂e1], :∂v0])
end

""" Add a triangle (2-simplex) to a simplicial set, given its boundary edges.

In the arguments to this function, the boundary edges have the order ``0 → 1``,
``1 → 2``, ``0 → 2``.

!!! warning

    This low-level function does not check the simplicial identities. It is your
    responsibility to ensure they are satisfied. By contrast, triangles added
    using the function [`glue_triangle!`](@ref) always satisfy the simplicial
    identities, by construction. Thus it is often easier to use this function.
"""
add_triangle!(s::HasDeltaSet2D, src2_first::Int, src2_last::Int, tgt2::Int; kw...) =
  add_part!(s, :Tri; ∂e0=src2_last, ∂e1=tgt2, ∂e2=src2_first, kw...)

""" Glue a triangle onto a simplicial set, given its boundary vertices.

If a needed edge between two vertices exists, it is reused (hence the "gluing");
otherwise, it is created.
"""
function glue_triangle!(s::HasDeltaSet2D, v₀::Int, v₁::Int, v₂::Int; kw...)
  add_triangle!(s, get_edge!(s, v₀, v₁), get_edge!(s, v₁, v₂),
                get_edge!(s, v₀, v₂); kw...)
end

function get_edge!(s::HasDeltaSet1D, src::Int, tgt::Int)
  es = edges(s, src, tgt)
  isempty(es) ? add_edge!(s, src, tgt) : first(es)
end

""" Glue a triangle onto a simplicial set, respecting the order of the vertices.
"""
function glue_sorted_triangle!(s::HasDeltaSet2D, v₀::Int, v₁::Int, v₂::Int; kw...)
  v₀, v₁, v₂ = sort(SVector(v₀, v₁, v₂))
  glue_triangle!(s, v₀, v₁, v₂; kw...)
end

# 2D oriented simplicial sets
#----------------------------

@present OrientedDeltaSchema2D <: DeltaCategory2D begin
  Orientation::AttrType
  edge_orientation::Attr(E,Orientation)
  tri_orientation::Attr(Tri,Orientation)
end

""" A two-dimensional oriented delta set.

Triangles are ordered in the cyclic order ``(0,1,2)`` when `tri_orientation` is
true/positive and in the reverse order when it is false/negative.
"""
@acset_type OrientedDeltaSet2D(OrientedDeltaSchema2D,
                               index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2]) <: AbstractDeltaSet2D

orientation(::Type{Val{2}}, s::HasDeltaSet2D, args...) =
  s[args..., :tri_orientation]
set_orientation!(::Type{Val{2}}, s::HasDeltaSet2D, t, orientation) =
  (s[t, :tri_orientation] = orientation)

function ∂_nz(::Type{Val{2}}, s::HasDeltaSet2D, t::Int)
  edges = triangle_edges(s,t)
  (edges, sign(2,s,t) * sign(1,s,edges) .* @SVector([1,-1,1]))
end

function d_nz(::Type{Val{1}}, s::HasDeltaSet2D, e::Int)
  sgn = sign(1, s, e)
  t₀, t₁, t₂ = coface(2,0,s,e), coface(2,1,s,e), coface(2,2,s,e)
  (lazy(vcat, t₀, t₁, t₂),
   lazy(vcat, sgn*sign(2,s,t₀), -sgn*sign(2,s,t₁), sgn*sign(2,s,t₂)))
end

# 2D embedded simplicial sets
#----------------------------

@present EmbeddedDeltaSchema2D <: OrientedDeltaSchema2D begin
  Point::AttrType
  point::Attr(V, Point)
end

""" A two-dimensional, embedded, oriented delta set.
"""
@acset_type EmbeddedDeltaSet2D(EmbeddedDeltaSchema2D,
                               index=[:∂v0,:∂v1,:∂e0,:∂e1,:∂e2]) <: AbstractDeltaSet2D

volume(::Type{Val{n}}, s::EmbeddedDeltaSet2D, x) where n =
  volume(Val{n}, s, x, CayleyMengerDet())
volume(::Type{Val{2}}, s::HasDeltaSet2D, t::Int, ::CayleyMengerDet) =
  volume(point(s, triangle_vertices(s,t)))

# General operators
###################

""" Wrapper for simplex or simplices of dimension `n`.

See also: [`V`](@ref), [`E`](@ref), [`Tri`](@ref).
"""
@parts_array_struct Simplex{n}

""" Vertex in simplicial set: alias for `Simplex{0}`.
"""
const V = Simplex{0}

""" Edge in simplicial set: alias for `Simplex{1}`.
"""
const E = Simplex{1}

""" Triangle in simplicial set: alias for `Simplex{2}`.
"""
const Tri = Simplex{2}

""" Wrapper for chain of oriented simplices of dimension `n`.
"""
@vector_struct SimplexChain{n}

const VChain = SimplexChain{0}
const EChain = SimplexChain{1}
const TriChain = SimplexChain{2}

""" Wrapper for discrete form, aka cochain, in simplicial set.
"""
@vector_struct SimplexForm{n}

const VForm = SimplexForm{0}
const EForm = SimplexForm{1}
const TriForm = SimplexForm{2}

""" Simplices of given dimension in a simplicial set.
"""
@inline simplices(n::Int, s::HasDeltaSet) = 1:nsimplices(Val{n}, s)

""" Number of simplices of given dimension in a simplicial set.
"""
@inline nsimplices(n::Int, s::HasDeltaSet) = nsimplices(Val{n}, s)

""" Face map and boundary operator on simplicial sets.

Given numbers `n` and `0 <= i <= n` and a simplicial set of dimension at least
`n`, the `i`th face map is implemented by the call

```julia
∂(n, i, s, ...)
```

The boundary operator on `n`-faces and `n`-chains is implemented by the call

```julia
∂(n, s, ...)
```

Note that the face map returns *simplices*, while the boundary operator returns
*chains* (vectors in the free vector space spanned by oriented simplices).
"""
@inline ∂(i::Int, s::HasDeltaSet, x::Simplex{n}) where n =
  Simplex{n-1}(face(Val{(n,i)}, s, x.data))
@inline ∂(n::Int, i::Int, s::HasDeltaSet, args...) =
  face(Val{(n,i)}, s, args...)

@inline coface(i::Int, s::HasDeltaSet, x::Simplex{n}) where n =
  Simplex{n+1}(coface(Val{(n+1,i)}, s, x.data))
@inline coface(n::Int, i::Int, s::HasDeltaSet, args...) =
  coface(Val{(n,i)}, s, args...)

∂(s::HasDeltaSet, x::SimplexChain{n}) where n =
  SimplexChain{n-1}(∂(Val{n}, s, x.data))
@inline ∂(n::Int, s::HasDeltaSet, args...) = ∂(Val{n}, s, args...)

function ∂(::Type{Val{n}}, s::HasDeltaSet, args...) where n
  operator_nz(Int, nsimplices(n-1,s), nsimplices(n,s), args...) do x
    ∂_nz(Val{n}, s, x)
  end
end

""" Alias for the face map and boundary operator [`∂`](@ref).
"""
const boundary = ∂

""" The discrete exterior derivative, aka the coboundary operator.
"""
d(s::HasDeltaSet, x::SimplexForm{n}) where n =
  SimplexForm{n+1}(d(Val{n}, s, x.data))
@inline d(n::Int, s::HasDeltaSet, args...) = d(Val{n}, s, args...)

function d(::Type{Val{n}}, s::HasDeltaSet, args...) where n
  operator_nz(Int, nsimplices(n+1,s), nsimplices(n,s), args...) do x
    d_nz(Val{n}, s, x)
  end
end

""" Alias for the coboundary operator [`d`](@ref).
"""
const coboundary = d

""" Alias for the discrete exterior derivative [`d`](@ref).
"""
const exterior_derivative = d

""" Orientation of simplex.
"""
orientation(s::HasDeltaSet, x::Simplex{n}) where n =
  orientation(Val{n}, s, x.data)
@inline orientation(n::Int, s::HasDeltaSet, args...) =
  orientation(Val{n}, s, args...)

@inline Base.sign(n::Int, s::HasDeltaSet, args...) = sign(Val{n}, s, args...)
Base.sign(::Type{Val{n}}, s::HasDeltaSet, args...) where n =
  numeric_sign.(orientation(Val{n}, s, args...))

numeric_sign(x) = sign(x)
numeric_sign(x::Bool) = x ? +1 : -1

""" Set orientation of simplex.
"""
@inline set_orientation!(n::Int, s::HasDeltaSet, args...) =
  set_orientation!(Val{n}, s, args...)

""" ``n``-dimensional volume of ``n``-simplex in an embedded simplicial set.
"""
volume(s::HasDeltaSet, x::Simplex{n}, args...) where n =
  volume(Val{n}, s, x.data, args...)
@inline volume(n::Int, s::HasDeltaSet, args...) = volume(Val{n}, s, args...)

""" Convenience function for linear operator based on structural nonzero values.
"""
operator_nz(f, ::Type{T}, m::Int, n::Int,
            x::Int, Vec::Type=SparseVector{T}) where T = fromnz(Vec, f(x)..., m)
operator_nz(f, ::Type{T}, m::Int, n::Int,
            vec::AbstractVector) where T = applynz(f, vec, m, n)
operator_nz(f, ::Type{T}, m::Int, n::Int,
            Mat::Type=SparseMatrixCSC{T}) where T = fromnz(f, Mat, m, n)

# Consistent orientation
########################

""" Consistently orient simplices in a simplicial set, if possible.

Two simplices with a common face are *consistently oriented* if they induce
opposite orientations on the shared face. This function attempts to consistently
orient all simplices of a given dimension and returns whether this has been
achieved. Each connected component is oriently independently using the helper
function [`orient_component!`](@ref).
"""
orient!(s::AbstractDeltaSet1D) = orient!(s, E)
orient!(s::AbstractDeltaSet2D) = orient!(s, Tri)

function orient!(s::HasDeltaSet, ::Type{Simplex{n}}) where n
  # Compute connected components as coequalizer of face maps.
  ndom, ncodom = nsimplices(n, s), nsimplices(n-1, s)
  face_maps = SVector{n+1}([ FinFunction(x -> ∂(n,i,s,x), ndom, ncodom)
                             for i in 0:n ])
  π = only(coequalizer(face_maps))

  # Choose an arbitrary representative of each component.
  reps = zeros(Int, length(codom(π)))
  for x in reverse(simplices(n, s))
    reps[π(∂(n,0,s,x))] = x
  end

  # Orient each component, starting at the chosen representative.
  init_orientation = one(eltype(orientation(n, s)))
  for x in reps
    orient_component!(s, Simplex{n}(x), init_orientation) || return false
  end
  true
end

""" Consistently orient simplices in the same connected component, if possible.

Given an ``n``-simplex and a choice of orientation for it, this function
attempts to consistently orient all ``n``-simplices that may be reached from it
by traversing ``(n-1)``-faces. The traversal is depth-first. If a consistent
orientation is possible, the function returns `true` and the orientations are
assigned; otherwise, it returns `false` and no orientations are changed.

If the simplicial set is not connected, the function [`orient!`](@ref) may be
more convenient.
"""
orient_component!(s::AbstractDeltaSet1D, e::Int, args...) =
  orient_component!(s, E(e), args...)
orient_component!(s::AbstractDeltaSet2D, t::Int, args...) =
  orient_component!(s, Tri(t), args...)

function orient_component!(s::HasDeltaSet, x::Simplex{n},
                           x_orientation::Orientation) where {n, Orientation}
  orientations = repeat(Union{Orientation,Nothing}[nothing], nsimplices(n, s))

  orient_stack = Vector{Pair{Int64, Orientation}}()

  push!(orient_stack, x[] => x_orientation)
  is_orientable = true
  while !isempty(orient_stack)
    x, target = pop!(orient_stack)
    current = orientations[x]
    if isnothing(current)
      # If not visited, set the orientation and add neighbors to stack.
      orientations[x] = target
      for i in 0:n, j in 0:n
        next = iseven(i+j) ? negate(target) : target
        for y in coface(n, j, s, ∂(n, i, s, x))
          y == x || push!(orient_stack, y=>next)
        end
      end
    elseif current != target
      is_orientable = false
      break
    end
  end

  if is_orientable
    component = findall(!isnothing, orientations)
    set_orientation!(n, s, component, orientations[component])
  end
  is_orientable
end

negate(x) = -x
negate(x::Bool) = !x

# Euclidean geometry
####################

""" ``n``-dimensional volume of ``n``-simplex spanned by given ``n+1`` points.
"""
function volume(points)
  CM = cayley_menger(points...)
  n = length(points) - 1
  sqrt(abs(det(CM)) / 2^n) / factorial(n)
end

""" Construct Cayley-Menger matrix for simplex spanned by given points.

For an ``n`-simplex, this is the ``(n+2)×(n+2)`` matrix that appears in the
[Cayley-Menger
determinant](https://en.wikipedia.org/wiki/Cayley-Menger_determinant).
"""
function cayley_menger(p0::V, p1::V) where V <: AbstractVector
  d01 = sqdistance(p0, p1)
  SMatrix{3,3}(0,  1,   1,
               1,  0,   d01,
               1,  d01, 0)
end
function cayley_menger(p0::V, p1::V, p2::V) where V <: AbstractVector
  d01, d12, d02 = sqdistance(p0, p1), sqdistance(p1, p2), sqdistance(p0, p2)
  SMatrix{4,4}(0,  1,   1,   1,
               1,  0,   d01, d02,
               1,  d01, 0,   d12,
               1,  d02, d12, 0)
end
function cayley_menger(p0::V, p1::V, p2::V, p3::V) where V <: AbstractVector
  d01, d12, d02 = sqdistance(p0, p1), sqdistance(p1, p2), sqdistance(p0, p2)
  d03, d13, d23 = sqdistance(p0, p3), sqdistance(p1, p3), sqdistance(p2, p3)
  SMatrix{5,5}(0,  1,   1,   1,   1,
               1,  0,   d01, d02, d03,
               1,  d01, 0,   d12, d13,
               1,  d02, d12, 0,   d23,
               1,  d03, d13, d23, 0)
end

""" Squared Euclidean distance between two points.
"""
sqdistance(x, y) = sum((x-y).^2)

end


text:
import Base.Order: lt;
using PolynomialRings
@ring! ℚ[p,n,d,q]
I = [p^5-n,p^10-d,p^25-q]
G = groebner_basis(I)
@show G
@show rem(p^317, G)

max_stars_repo_name: aabrandt/pennies

max_issues_repo_name: aabrandt/pennies

max_forks_repo_name: aabrandt/pennies


text:
# import Base: copyto!, isequal, iterate, getindex, length, ndims, size
#
# ##
# # TODO put it as a parameter?
# # constants
# learning_relations = [
#     IdentityIntervalRelation(BAlgebra)
#     ExistentialIntervalRelation{:G}(BAlgebra)
#     ExistentialIntervalRelation{:L}(BAlgebra)
#     ExistentialIntervalRelation{:A}(BAlgebra)
#     ExistentialIntervalRelation{:O}(BAlgebra)
#     ExistentialIntervalRelation{:E}(BAlgebra)
#     ExistentialIntervalRelation{:D}(BAlgebra)
#     ExistentialIntervalRelation{:E}(BAlgebra)
#     ExistentialIntervalRelation{:InvL}(BAlgebra)
#     ExistentialIntervalRelation{:InvA}(BAlgebra)
#     ExistentialIntervalRelation{:InvO}(BAlgebra)
#     ExistentialIntervalRelation{:InvE}(BAlgebra)
#     ExistentialIntervalRelation{:InvD}(BAlgebra)
#     ExistentialIntervalRelation{:InvE}(BAlgebra)
#     UniversalIntervalRelation{:G}(BAlgebra)
#     UniversalIntervalRelation{:L}(BAlgebra)
#     UniversalIntervalRelation{:A}(BAlgebra)
#     UniversalIntervalRelation{:O}(BAlgebra)
#     UniversalIntervalRelation{:E}(BAlgebra)
#     UniversalIntervalRelation{:D}(BAlgebra)
#     UniversalIntervalRelation{:E}(BAlgebra)
#     UniversalIntervalRelation{:InvL}(BAlgebra)
#     UniversalIntervalRelation{:InvA}(BAlgebra)
#     UniversalIntervalRelation{:InvO}(BAlgebra)
#     UniversalIntervalRelation{:InvE}(BAlgebra)
#     UniversalIntervalRelation{:InvD}(BAlgebra)
#     UniversalIntervalRelation{:InvE}(BAlgebra)
# ]
#
# order_relations = [≤, <, >, ≥]
#
# ##
# # TODO we need to see if there is a need for two algebras, i.e., perhaps AbstractHeytingValue -> AbstractAlgebraValue
# abstract type AbstractHeytingValue end
#
# struct HeytingValue{T<:Real} <: AbstractHeytingValue
#     value::T
# end
#
# function Base.show(io::IO, v::HeytingValue)
#     if v.value == 0.0
#         print(io, "⊥")
#     elseif v.value == 1.0
#         print(io, "⊤")
#     else
#         print(io, v.value)
#     end
# end
#
# ##
# # TODO <: HeytingAlgebra?
# mutable struct Node{T}
#     data::T
#     parent::Union{Node{<:AbstractProposition}, Node{<:AbstractRelation}, Node{HeytingValue}}
#     left::Union{Node{<:AbstractProposition}, Node{<:AbstractRelation}, Node{HeytingValue}}
#     right::Union{Node{<:AbstractProposition}, Node{<:AbstractRelation}, Node{HeytingValue}}
#     # Root constructor.
#     Node{T}(data) where T = new{T}(data)
#     Node(data) = new{typeof(data)}(data)
# end
#
# parentnode(currnode::Node) = getproperty(currnode.parent)
# leftnode(currnode::Node) = getproperty(currnode.left)
# rightnode(currnode::Node) = getproperty(currnode.right)
# parentnode!(currnode::Node, pnode::Node) = setproperty!(currnode.parent, pnode)
# leftnode!(currnode::Node, lnode::Node) = setproperty!(currnode.left, lnode)
# rightnode!(currnode::Node, rnode::Node) = setproperty!(currnode.right, rnode)
#
# ##
# abstract type AbstractRule end
#
# mutable struct ClassificationRule <: AbstractRule
#     antecedent::Node
#     consequent::String
# end
#
# antecedent(r::ClassificationRule) = r.antecedent
# consequent(r::ClassificationRule) = r.consequent
#
# function size(rule::ClassificationRule)
#     return _size(antecedent(rule))
# end
#
# function _size(node::Node)
#     s = 1
#     sl = isdefined(node, :left) ? _size(node.left) : 0
#     sr = isdefined(node, :right) ? _size(node.right) : 0
#     return s + sl + sr
# end
#
# function upwards_modaldepth(node::Node)
#     md = isdefined(node, :parent) ? upwards_modaldepth(node.parent) : 0
#     md += typeof(node.data) <:AbstractModalRelation ? 1 : 0
#     return md
# end
#
# function is_leaf(node::Node)
#     return !isdefined(node, :left) && !isdefined(node, :right) ? true : false
# end
#
# function is_internal_node(node::Node)
#     return !is_leaf(node) ? true : false
# end
#
# Base.show(io::IO, n::Node) = print(io, _repr_antecedent(n))
#
# function _repr_antecedent(n::Node)
#     if typeof(n.data) <: AbstractUniversalIntervalRelation
#         return repr(n.data) * "(" * _repr_antecedent(n.left) * " → " * _repr_antecedent(n.right) * ")"
#     elseif typeof(n.data) <: AbstractExistentialIntervalRelation
#         return repr(n.data) * "(" * _repr_antecedent(n.left) * " ∧ " * _repr_antecedent(n.right) * ")"
#     else # Propositions and Heyting values.
#         return repr(n.data)
#     end
# end
#
# function _repr_consequent(c::String)
#     return c
# end
#
# function _repr_classificationrule(r::ClassificationRule)
#     return _repr_antecedent(antecedent(r)) * " ⟹  " * _repr_consequent(consequent(r))
# end
#
# Base.show(io::IO, r::ClassificationRule) = print(io, _repr_classificationrule(r))
#
# ##
# mutable struct ClassificationRules
#     rules::Vector{ClassificationRule}
#     modalities::Vector{<:AbstractModalRelation}
#     minmd::Int64
#     maxnumrules::Int64
#
#     ClassificationRules(rules;
#         modalities::Vector{<:AbstractModalRelation}=learning_relations,
#         minmd::Int64=3,
#         maxnumrules::Int64=10) = new(rules, modalities, minmd, maxnumrules)
# end
#
# rules(Γ::ClassificationRules) = Γ.rules
# modalities(Γ::ClassificationRules) = Γ.modalities
# minmd(Γ::ClassificationRules) = Γ.minmd
# # maxmd(Γ::ClassificationRules) = Γ.maxmd
# maxnumrules(Γ::ClassificationRules) = Γ.maxnumrules
# classes(Γ::ClassificationRules) = collect(Set([consequent(ρ) for ρ ∈ Γ]))
# mean_size(Γ::ClassificationRules) = reduce(+,[size(ρ) for ρ ∈ rules(Γ)])/length(Γ)
# getindex(Γ::ClassificationRules, inds...) = getindex(rules(Γ), inds...)
#
#
# function _repr_rules(Γ::ClassificationRules)
#     res = ""
#     for ρ ∈ rules(Γ)
#         res *= _repr_classificationrule(ρ) * "\n"
#     end
#     return res[1:end-1] # remove the last "\n"
# end
#
# Base.show(io::IO, Γ::ClassificationRules) = print(io, _repr_rules(Γ))
# # p = LowLevelProposition(1, ≤, 0.7, 5)
#
# ##
# # Defined, otherwise it gives errors
# length(crules::ClassificationRules) = length(rules(crules))
# iterate(crules::ClassificationRules) = iterate(rules(crules))
# iterate(crules::ClassificationRules, i::Int64) = iterate(rules(crules), i)
# ndims(::Type{ClassificationRules}) = 1
# size(crules::ClassificationRules) = (length(crules),)
# copyto!(crules::ClassificationRules, b::Base.Broadcast.Broadcasted) = copyto!(rules(crules), b)
#
# ##
# function _rand_levelproposition(numattributes::Int64, relations::Array{Function, 1}, domains::Array{Array{Float64, N} where N, 1})
#     i = rand(1:numattributes)
#     fun = [min, max][rand(1:2)]
#     return Proposition(fun, i, relations[rand(1:length(relations))], domains[i][rand(1:length(domains[i]))])
# end
# FINO A QUA

##
# Abstract relation.
abstract type AbstractRelation end

##
# Abstract proposotional relation.
abstract type AbstractPropositionalRelation     <: AbstractRelation end
# Abstract modal relation.
abstract type AbstractModalRelation             <: AbstractRelation end

##
# Abstract temporal relation.
abstract type AbstractTemporalRelation          <: AbstractModalRelation end
# Abstract spatial relation.
abstract type AbstractSpatialRelation           <: AbstractModalRelation end
# Abstract spatio-temporal relation.
abstract type AbstractSpatioTemporalRelation    <: AbstractModalRelation end

##
# Abstract point (temporal) relation.
abstract type AbstractPointRelation             <: AbstractTemporalRelation end
# Abstract interval (temporal) relation.
abstract type AbstractIntervalRelation          <: AbstractTemporalRelation end
# Abstract rectangle (spatial) relation.
abstract type AbstractRectangleRelation         <: AbstractSpatialRelation end
# Abstract cube (spatio-temporal) relation.
abstract type AbstractCubeRelation              <: AbstractSpatioTemporalRelation end

##
# Abstract existential interval (temporal) relation.
abstract type AbstractExistentialIntervalRelation   <: AbstractIntervalRelation end
# Abstract universal interval (temporal) relation.
abstract type AbstractUniversalIntervalRelation     <: AbstractIntervalRelation end

##
#
struct IdentityIntervalRelation         <: AbstractIntervalRelation end
# Using the singleton design pattern; similar to Val{T}.
struct ExistentialIntervalRelation{T}   <: AbstractExistentialIntervalRelation end
struct UniversalIntervalRelation{T}     <: AbstractUniversalIntervalRelation end
ExistentialIntervalRelation(s::AbstractString)  = ExistentialIntervalRelation{Symbol(s)}()
ExistentialIntervalRelation(s::Symbol)          = ExistentialIntervalRelation{s}()
UniversalIntervalRelation(s::AbstractString)    = UniversalIntervalRelation{Symbol(s)}()
UniversalIntervalRelation(s::Symbol)            = UniversalIntervalRelation{s}()

##
Base.show(io::IO, ::IdentityIntervalRelation)           = print(io, "|=|")

Base.show(io::IO, ::ExistentialIntervalRelation{:G})    = print(io, "⟨G⟩")

Base.show(io::IO, ::ExistentialIntervalRelation{:L})    = print(io, "⟨L⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:A})    = print(io, "⟨A⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:O})    = print(io, "⟨O⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:E})    = print(io, "⟨E⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:D})    = print(io, "⟨D⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:B})    = print(io, "⟨B⟩")

Base.show(io::IO, ::ExistentialIntervalRelation{:InvL}) = print(io, "⟨InvL⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:InvA}) = print(io, "⟨InvA⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:InvO}) = print(io, "⟨InvO⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:InvE}) = print(io, "⟨InvE⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:InvD}) = print(io, "⟨InvD⟩")
Base.show(io::IO, ::ExistentialIntervalRelation{:InvB}) = print(io, "⟨InvB⟩")

Base.show(io::IO, ::UniversalIntervalRelation{:G})      = print(io, "[G]")

Base.show(io::IO, ::UniversalIntervalRelation{:L})      = print(io, "[L]")
Base.show(io::IO, ::UniversalIntervalRelation{:A})      = print(io, "[A]")
Base.show(io::IO, ::UniversalIntervalRelation{:O})      = print(io, "[O]")
Base.show(io::IO, ::UniversalIntervalRelation{:E})      = print(io, "[E]")
Base.show(io::IO, ::UniversalIntervalRelation{:D})      = print(io, "[D]")
Base.show(io::IO, ::UniversalIntervalRelation{:B})      = print(io, "[B]")

Base.show(io::IO, ::UniversalIntervalRelation{:InvL})   = print(io, "[InvL]")
Base.show(io::IO, ::UniversalIntervalRelation{:InvA})   = print(io, "[InvA]")
Base.show(io::IO, ::UniversalIntervalRelation{:InvO})   = print(io, "[InvO]")
Base.show(io::IO, ::UniversalIntervalRelation{:InvE})   = print(io, "[InvE]")
Base.show(io::IO, ::UniversalIntervalRelation{:InvD})   = print(io, "[InvD]")
Base.show(io::IO, ::UniversalIntervalRelation{:InvB})   = print(io, "[InvB]")

using ExprRules

G = @grammar begin
    S = ExistentialIntervalRelation{:G}(F) | UniversalIntervalRelation{:G}(F)
    F = ExistentialIntervalRelation{:L}(and(F,F))
    F = UniversalIntervalRelation{:L}(implies(F,F))
    F = implies(F,F)
    F = and(F,F)
    F = p(sum, A, ≤, a)
    A = "fever"
    # F = 1
    F = H
    H = 0.0 | 0.25 | 0.5 | 0.75 | 1.0
    a = |([i for i in 1:10])
end

p(f::Function, A, ≤, a) = Base.repr("f(A) ≤ a")
implies(a::Float64, b::Float64) = a ≤ b ? float(1) : b
and(a::Float64, b::Float64) = min(a,b)

using GBDTs

using MultivariateTimeSeries
X, y = read_data_labeled(joinpath(dirname(pathof(GBDTs)), "..", "data", "auslan_youtube8"));

grammar = @grammar begin
    b = G(bvec) | F(bvec) | G(implies(bvec,bvec))
    bvec = and(bvec, bvec)
    bvec = or(bvec, bvec)
    bvec = not(bvec)
    bvec = lt(rvec, rvec)
    bvec = lte(rvec, rvec)
    bvec = gt(rvec, rvec)
    bvec = gte(rvec, rvec)
    bvec = f_lt(x, xid, v, vid)
    bvec = f_lte(x, xid, v, vid)
    bvec = f_gt(x, xid, v, vid)
    bvec = f_gte(x, xid, v, vid)
    rvec = x[xid]
    xid = |(["x_1","y_1","z_1","roll_1","pitch_1","yaw_1","thumbbend_1","forebend_1","middlebend_1","ringbend_1","littlebend_1","x_2","y_2","z_2","roll_2","pitch_2","yaw_2","thumbbend_2","forebend_2","middlebend_2","ringbend_2","littlebend_2"])
    # xid = |(["Fce", "Fsfc", "G2", "MAD", "T3-2", "T13-12", "T16-15", "T13-5", "T14-6", "T15-7"])
    vid = |(1:10)
end

grammar = @grammar begin
    b = G(bvec) | F(bvec) | G(implies(bvec,bvec))
    bvec = and(bvec, bvec)
    bvec = or(bvec, bvec)
    bvec = not(bvec)
    bvec = lt(rvec, rvec)
    bvec = lte(rvec, rvec)
    bvec = gt(rvec, rvec)
    bvec = gte(rvec, rvec)
    bvec = f_lt(x, xid, v, vid)
    bvec = f_lte(x, xid, v, vid)
    bvec = f_gt(x, xid, v, vid)
    bvec = f_gte(x, xid, v, vid)
    rvec = x[xid]
    xid = |(["Fce", "Fsfc", "G2", "MAD", "T3-2", "T13-12", "T16-15", "T13-5", "T14-6", "T15-7"])
    vid = |(1:10)
end

G(v) = all(v)                                                #globally
F(v) = any(v)                                                #eventually
f_lt(x, xid, v, vid) = lt(x[xid], v[xid][vid])               #feature is less than a constant
f_lte(x, xid, v, vid) = lte(x[xid], v[xid][vid])             #feature is less than or equal to a constant
f_gt(x, xid, v, vid) = gt(x[xid], v[xid][vid])               #feature is greater than a constant
f_gte(x, xid, v, vid) = gte(x[xid], v[xid][vid])             #feature is greater than or equal to a constant

#workarounds for slow dot operators:
implies(v1, v2) = (a = similar(v1); a .= v2 .| .!v1)         #implies
not(v) = (a = similar(v); a .= .!v)                          #not
and(v1, v2) = (a = similar(v1); a .= v1 .& v2)               #and
or(v1, v2) = (a = similar(v1); a .= v1 .| v2)                #or
lt(x1, x2) = (a = Vector{Bool}(undef,length(x1)); a .= x1 .< x2)   #less than
lte(x1, x2) = (a = Vector{Bool}(undef,length(x1)); a .= x1 .≤ x2)  #less than or equal to
gt(x1, x2) = (a = Vector{Bool}(undef,length(x1)); a .= x1 .> x2)   #greater than
gte(x1, x2) = (a = Vector{Bool}(undef,length(x1)); a .= x1 .≥ x2)  #greater than or equal to

const v = Dict{String,Vector{Float64}}()
mins, maxes = minimum(X), maximum(X)
for (i,xid) in enumerate(names(X))
    v[xid] = collect(range(mins[i],stop=maxes[i],length=10))
end;

p = MonteCarlo(2, 5)

using Random; Random.seed!(1)
model = induce_tree(grammar, :b, p, X, y, 6);


text:
# This file is part of Kpax3. License is MIT.

#=
References:

Tajima F. and Nei M. (1984). Estimation of evolutionary distance between
nucleotide sequences. Mol Biol Evol 1 (3):269-285.
http://mbe.oxfordjournals.org/content/1/3/269

Tamura K. and Kumar S. (2002). Evolutionary Distance Estimation Under
Heterogeneous Substitution Pattern Among Lineages. Mol Biol Evol 19 (10):
1727-1736. http://mbe.oxfordjournals.org/content/19/10/1727
=#

#=
Description:

Compute Tajima Nei (1984) pairwise distances of protein sequences, with the
Tamura and Kumar (2002) correction for heterogeneous patterns.

Arguments:

  data::Matrix{UInt8}
    m-by-n data matrix, where m is the common sequence length and n is the
    sample size
  ref::Vector{UInt8}
    reference sequence, i.e. a vector of length M storing the values of
    homogeneous sites

Details:

Pyrrolysine (O) and Selenocysteine (U) are treated as valid amino acids if they
have not been given a missing value.

If a pairwise distance is equal to -1.0, it means that is wasn't possible to
compute it. This usually happens when the hypotheses of the underlying
evolutionary model are not satisfied.

Value:

  d::Vector{Float64}
    evolutionary distances. Vector length is equal to n * (n - 1) / 2. It
    contains the values of the lower triangular matrix of the full distance
    matrix, ordered by column.
    To access the distance between units i and j (i < j), use
    d[n * (i - 1) - div(i * (i - 1), 2) + j - i]
=#
function distaamtn84(data::Matrix{UInt8},
                     ref::Vector{UInt8})
  (m, n) = size(data)

  d = zeros(Float64, div(n * (n - 1), 2))

  tmpref = zeros(Float64, 127)
  tmpraw = zeros(Float64, 127)

  for i in 1:length(ref)
    tmpref[ref[i]] += 1
  end

  for col in 1:n, row in 1:m
    tmpraw[data[row, col]] += 1
  end

  aaset = Int['a', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'k', 'l', 'm', 'n', 'o',
              'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'y']

  gt = zeros(Float64, 127)
  gb = zeros(Float64, 127)
  st = 0.0
  sb = 0.0
  sc = 0.0

  for aa in aaset
    gt[aa] = tmpref[aa]
    gb[aa] = tmpraw[aa] + n * gt[aa]

    st += gt[aa]
    sb += gb[aa]
  end

  for aa in aaset
    gb[aa] /= sb
    sc += gb[aa]^2
  end

  v = 1 - sc

  idx = 1
  for j in 1:(n - 1), i in (j + 1):n
    d[idx] = aamtn84(i, j, data, gt, st, v)
    idx += 1
  end

  d
end

#=
Description:

Compute the Tajima Nei (1984) distance between two protein sequences, with the
Tamura and Kumar (2002) correction for heterogeneous patterns.

Arguments:

  i::Int
    index of first sequence
  j::Int
    index of second sequence
  data::Matrix{UInt8}
    m-by-n data matrix, where m is the common sequence length and n is the
    sample size
  gt::Vector{Float64}
    common absolute frequency for the 22 amino acids, i.e. the count of each
    amino acid at homogeneous sites
  h::Float64
    total number of homogeneous sites that are not missing
  v::Float64
    v = 1 - sum(pi^2), where pi is the proportion of amino acid i observed in
    the whole dataset

Value:

  d::Float64
    evolutionary distance between the two protein sequences
=#
function aamtn84(i::Int,
                 j::Int,
                 data::Matrix{UInt8},
                 gt::Vector{Float64},
                 h::Float64,
                 v::Float64)
  d = -1.0

  # effective length, i.e. total number of sites at which both sequences have
  # non-missing values
  n = fill(h, 3)

  # proportion of different elements
  p = 0.0

  # proportion of observed amino acids
  g1 = copy(gt)
  g2 = copy(gt)

  f = 0.0
  w = 0.0

  x1 = 0x00
  x2 = 0x00

  for b in 1:size(data, 1)
    x1 = data[b, i]
    x2 = data[b, j]

    if UInt8(96) < x1 < UInt8(123)
      g1[x1] += 1
      n[1] += 1
    end

    if UInt8(96) < x2 < UInt8(123)
      g2[x2] += 1
      n[2] += 1
    end

    if (UInt8(96) < x1 < UInt8(123)) && (UInt8(96) < x2 < UInt8(123))
      if x1 != x2
        p += 1
      end

      n[3] += 1
    end
  end

  if p > 0
    g1 /= n[1]
    g2 /= n[2]

    p /= n[3]

    for i in 1:length(g1)
      f += g1[i] * g2[i]
    end

    x = 1 - (p / (1 - f))

    if x > 0
      d = - v * log(x)
    end
  else
    # sequences are identical (or couldn't be compared because of missing
    # values. But in this case, by default, we consider them identical)
    d = 0.0
  end

  d
end


text:
"""
    SubscriberFlags
Some useful flags when dealing with subscribers. Describes the state of the system.
Particularly helpful when the subscriber is actively receiving messages in another
thread and you want to query the state of the subscriber.
"""
Base.@kwdef mutable struct SubscriberFlags
    "Is the subscriber currently waiting to receive data"
    isreceiving::Bool = false

    "Did the subscriber exit with an error"
    diderror::Bool = false

    "Has the subscriber received a message"
    hasreceived::Bool = false

    "# Of bytes of last message recieved"
    bytesrecieved::Int64 = 0
end

abstract type Subscriber end

Base.isopen(sub::Subscriber)::Nothing =
    error("The `isopen` method hasn't been implemented for your Subscriber yet!")
Base.close(sub::Subscriber)::Nothing =
    error("The `close` method hasn't been implemented for your Subscriber yet!")
# TODO: probably should delete this
forceclose(sub::Subscriber)::Nothing =
    error("The `forceclose` method hasn't been implemented for your Subscriber yet!")

getname(sub::Subscriber)::String = sub.name
getflags(sub::Subscriber)::SubscriberFlags = sub.flags
portstring(sub::Subscriber)::String = ""

# Keep track of newly recieved message
bytesreceived(sub::Subscriber)::Int64 = getflags(sub).bytesrecieved
has_new(sub::Subscriber)::Bool = getflags(sub).hasreceived

"""
Set the has received flag on subscriber to true
"""
function got_new!(sub::Subscriber)
    flags = getflags(sub)
    flags.hasreceived = true
    return flags.hasreceived
end

"""
Set the has received flag on subscriber to false
"""
function read_new!(sub::Subscriber)
    flags = getflags(sub)
    flags.hasreceived = false
    return flags.hasreceived
end

"""
Read in the byte data into the message container buf. Returns the number of bytes read
"""
function decode!(buf::ProtoBuf.ProtoType, bin_data::IOBuffer)
    bytes_read = length(bin_data.data)
    ProtoBuf.readproto(bin_data, buf)

    return bytes_read
end

function decode!(buf::AbstractVector{UInt8}, bin_data::IOBuffer)
    bytes_read = min(length(buf), bin_data.size)
    copyto!(buf, 1, bin_data.data, 1, bytes_read)

    return bytes_read
end

"""
Receive function, is expected to return true if a message was received, false otherwise
"""
function receive(sub::Subscriber, buf::MercuryMessage)::Bool
    throw(
        MercuryException(
            "The `receive` method hasn't been implemented for your Subscriber yet!",
        ),
    )
end

"""
Specifies a subcriber along with specific message type.
This is useful for tracking multiple messages at once
"""
struct SubscribedMessage
    # Handle case in which listening for byte array or for protobuf (see decode!)
    msg::MercuryMessage
    sub::Subscriber          # Note this is an abstract type
    name::String

    function SubscribedMessage(msg::MercuryMessage, sub::Subscriber; name = getname(sub))
        new(msg, sub, name)
    end
end

function receive(submsg::SubscribedMessage)::Bool
    # Check if we recieved a new message, if we do set has recieved flag to true
    if (receive(submsg.sub, submsg.msg))
        got_new!(submsg.sub)
    end
    return has_new(submsg.sub)
end

@inline getname(submsg::SubscribedMessage) = submsg.name
@inline getcomtype(submsg::SubscribedMessage) = getcomtype(submsg.sub)

function printstatus(sub::SubscribedMessage; indent = 0)
    prefix = " "^indent
    println(prefix, "Subscriber: ", getname(sub))
    println(prefix, "  Type: ", getcomtype(sub))
    println(prefix, "  Message Type: ", typeof(sub.msg))
    println(prefix, "  Has received? ", getflags(sub.sub).hasreceived)
    println(prefix, "  Is Open? ", isopen(sub.sub))
end

"""
    on_new(func::Function, submsg::SubscribedMessage)
Helpful function for executing code blocks when a SubscribedMessage type has recieved a
new message on its Subscriber's Socket. The function func is expected to have a signature
of `func(msg::ProtoBuf.ProtoType)` where msg is the message which `submsg` has recieved.

Example:
```
on_new(nodeio.subs[1]) do msg
    println(msg.pos_x)
end
```
"""
function on_new(func::Function, submsg::SubscribedMessage)
    if has_new(submsg.sub)
        func(submsg.msg)

        read_new!(submsg.sub)
    end
end


text:
include("@aml_parse/literals.jl")
include("@aml_parse/custom_code.jl")

"""
@aml parser function
"""
function aml_parse(expr::Expr)

    expr.head == :struct || error("Invalid usage of aml_parse")

    # reminder:
    # var is a symbol
    # var::T or anything more complex is an expression

    is_struct_mutable = expr.args[1]
    T = expr.args[2] # Type name +(curly braces)

    argsexpr = expr.args[3] # arguments of the type

    # TODO: optimize and fuse these
    areargs_inds = Base.findall(x->!(isa(x, LineNumberNode)), argsexpr.args)
    macronum = count(x-> isa(x, Tuple{Symbol, Expr}), argsexpr.args)

    data = argsexpr.args[areargs_inds]

    datanum = length(areargs_inds)
    argsnum = datanum - (1 + macronum) # 1 for struct name

    args_param = Vector{Union{Expr,Symbol}}(undef, argsnum) # Expr(:parameters)[]
    args_var = Vector{Union{Expr,Symbol}}(undef, argsnum)
    args_defaultvalue = Vector{Any}(missing, argsnum)
    args_type = Vector{Union{Type, Symbol, Expr}}(undef, argsnum)
    args_name =Vector{Union{Missing,String}}(missing, argsnum)
    args_function = Vector{Union{Missing, Symbol, Function}}(missing, argsnum)
    args_literaltype = Vector{Union{Missing, Type}}(missing, argsnum)
    args_custom_creator = Vector{Union{Nothing, Expr}}(nothing, argsnum+1) # +1 for the end
    args_custom_extractor = Vector{Union{Nothing, Expr}}(nothing, argsnum+1)
    args_custom_updater = Vector{Union{Nothing, Expr}}(nothing, argsnum+1)
    struct_name = "name"
    struct_nodetype = AbsDocOrNode
    struct_function = Array{Union{Missing, Symbol, Function},0}(missing)

    iMacro = 0
    for iData = 1:datanum # iterating over arguments of each type argument

        iArg = iData - (1+ iMacro)
        i = areargs_inds[iData] # used for argsexpr.args[i]

        ei = data[iData] # type argument element i

        ########################
        # Single struct name - "aml name"
        if isa(ei, String)

            # struct_function[1]=missing # function

            # Self-name checker
            if ei == "~"
                if T isa Symbol
                    struct_name = string(T)
                elseif T isa Expr && T.head == :curly
                    struct_name = string(T.args[1]) # S
                end
            else
                struct_name = ei  # Type aml name
            end

            argsexpr.args[i] =  nothing # removing "aml name" from expr args
            struct_nodetype = AbsNormal

        ################################################################
        # Literal Struct name - empty"aml name"
        elseif isa(ei, Tuple)
            ################################################################
            # Struct aml
            ########################
            # Literal only  empty"aml name"
            if isa(ei, Tuple{Type,String})

                # struct_function[1]=missing # function

                # Self-name checker
                if ei[2] == "~"
                    if T isa Symbol
                        struct_name = string(T)
                    elseif T isa Expr && T.head == :curly
                        struct_name = string(T.args[1]) # S
                    end
                else
                    struct_name = ei[2] # Type aml name
                end

                struct_nodetype = ei[1]
                struct_nodetype = aml_dispatch(struct_nodetype, struct_name)

                argsexpr.args[i] =  nothing # removing "aml name" from expr args

            # Custom Code
            elseif isa(ei, Tuple{Symbol, Expr})
                iMacro += 1
                # Row for code insertion (insert before iArg-th argument)
                if ei[1] == :creator
                    args_custom_creator[iArg] = ei[2]
                elseif ei[1]  == :extractor
                    args_custom_extractor[iArg] = ei[2]
                elseif ei[1]  == :updater
                    args_custom_updater[iArg] = ei[2]
                end

                argsexpr.args[i] =  nothing # removing custom code macro from expr args

            end

        elseif ei isa Expr && ei.head == :tuple
            ########################
            # Struct Function - "aml name", F
            if isa(ei.args[1], String) && isa(ei.args[2], Union{Symbol,Function}) # "aml name", F

                struct_function[1]=ei.args[2] # function

                # Self-name checker
                if ei.args[1] == "~"
                    if T isa Symbol
                        struct_name = string(T)
                    elseif T isa Expr && T.head == :curly
                        struct_name = string(T.args[1]) # S
                    end
                else
                    struct_name = ei.args[1] # Type aml name
                end

                struct_nodetype = AbsNormal
                argsexpr.args[i] =  nothing # removing "aml name" from expr args

            ########################
            # Literal and Struct Function - empty"aml name", F
            elseif isa(ei.args[1], Tuple)  && isa(ei.args[2], Union{Symbol,Function})

                struct_function[1]=ei.args[2] # function

                # Self-name checker
                if ei.args[1][2] == "~"
                    if T isa Symbol
                        struct_name = string(T)
                    elseif T isa Expr && T.head == :curly
                        struct_name = string(T.args[1]) # S
                    end
                else
                    struct_name = ei.args[1][2] # Type aml name
                end

                struct_nodetype = ei.args[1][1]
                struct_nodetype = aml_dispatch(struct_nodetype, struct_name)

                argsexpr.args[i] =  nothing # removing "aml name" from expr args
        ################################################################
        # Arguments
        ########################
            # No Def Value
            elseif ei.args[1] isa Union{Symbol,Expr} # var/var::T, "name"

                # Def Value
                # args_defaultvalue[iArg] = missing

                # Type Checker
                lhs = ei.args[1]
                if lhs isa Symbol #  var, "name"

                    var = ei.args[1]

                    args_type[iArg] = String # consider String as the type
                    args_param[iArg] = var
                    args_var[iArg] = var

                    argsexpr.args[i] = var  # removing "name",...

                elseif lhs isa Expr && lhs.head == :(::) && lhs.args[1] isa Symbol # var::T, "name"

                    var = lhs.args[1]
                    varType = lhs.args[2] # Type

                    args_type[iArg] = varType
                    args_param[iArg] =  var
                    args_var[iArg] =  var

                    argsexpr.args[i] = lhs  # removing "name",...

                end

                # Literal Checker
                if length(ei.args[2]) == 2 # literal

                    argAmlType = ei.args[2][1]
                    args_literaltype[iArg] = argAmlType # literal type

                    ni = ei.args[2][2]

                    # Self-name checker
                    if ni == "~"
                        args_name[iArg] = string(var)
                    else
                        args_name[iArg] = ni
                    end

                else
                    args_literaltype[iArg] = AbsNormal # non-literal

                    ni = ei.args[2]

                    # Self-name checker
                    if ni == "~"
                        args_name[iArg] = string(var)
                    else
                        args_name[iArg] = ni
                    end
                end

                # Function Checker
                if length(ei.args) == 3 && isa(ei.args[3], Union{Function, Symbol}) #  var/var::T, "name", f

                    fun = ei.args[3]   # function
                    args_function[iArg] = fun

                # else # function name isn't given
                    # args_function[iArg] =  missing
                end

            end  # end Tuple sub possibilities
        ################################################################
        # Def Value
        elseif ei isa Expr && ei.head == :(=) # def value provided

            # aml name Checker
            if ei.args[2] isa Expr && ei.args[2].head == :tuple # var/var::T = defVal, "name"

                # Def Value
                defVal = ei.args[2].args[1]

                args_defaultvalue[iArg] = defVal

                lhs = ei.args[1]

                argsexpr.args[i] = lhs # remove =defVal for type definition

                # Type Checker
                if lhs isa Symbol #  var = defVal, "name"

                    var = ei.args[1]

                    args_type[iArg] = String # consider String as the type
                    args_param[iArg] = Expr(:kw, var, defVal)
                    args_var[iArg] =  var

                    argsexpr.args[i] = var  # removing "name",...

                elseif lhs isa Expr && lhs.head == :(::) && lhs.args[1] isa Symbol # var::T = defVal, "name"

                    var = lhs.args[1]
                    varType = lhs.args[2] # Type

                    args_type[iArg] = varType
                    args_param[iArg] = Expr(:kw, var, defVal) # TODO also put type expression
                    args_var[iArg] =  var

                    argsexpr.args[i] = lhs  # removing "name",...

                end

                # Literal Checker
                if length(ei.args[2].args[2]) == 2 # literal

                    argAmlType = ei.args[2].args[2][1]
                    args_literaltype[iArg] = argAmlType # literal type

                    ni = ei.args[2].args[2][2]

                    # Self-name checker
                    if ni == "~"
                        args_name[iArg] = string(var)
                    else
                        args_name[iArg] = ni
                    end

                else
                    args_literaltype[iArg] = AbsNormal # non-literal

                    ni = ei.args[2].args[2]

                    # Self-name checker
                    if ni == "~"
                        args_name[iArg] = string(var)
                    else
                        args_name[iArg] = ni
                    end

                end

                # Function Checker
                if length(ei.args[2].args) == 3 && isa(ei.args[2].args[3], Union{Function, Symbol}) #  var/var::T  = defVal, "name", f

                    fun = ei.args[2].args[3]  # function
                    args_function[iArg] = fun

                # else # function name isn't given
                    # args_function[iArg] = missing
                end

            ########################
            #  No aml Name - But defVal
            else # var/var::T = defVal # ignored for creating aml

                # Type Checker
                lhs = ei.args[1]
                if lhs isa Symbol #  var = defVal

                    defVal = ei.args[2]

                    args_defaultvalue[iArg] = defVal
                    # args_name[iArg] = missing # ignored for creating aml
                    # args_function[iArg] =  missing # ignored for creating aml

                    var = ei.args[1]

                    args_type[iArg] = Any
                    args_param[iArg] = Expr(:kw, var, defVal)
                    args_var[iArg] =  var

                    argsexpr.args[i] = var # remove =defVal for type definition

                elseif lhs isa Expr && lhs.head == :(::) && lhs.args[1] isa Symbol # var::T = defVal

                    defVal = ei.args[2]

                    args_defaultvalue[iArg] = defVal
                    # args_name[iArg] = missing # ignored for creating aml
                    # args_function[iArg] = missing # ignored for creating aml

                    var = lhs.args[1]
                    varType = lhs.args[2] # Type

                    args_type[iArg] = varType
                    args_param[iArg] =  Expr(:kw, var, defVal) # TODO also put type expression
                    args_var[iArg] =  var

                    argsexpr.args[i] = lhs # remove =defVal for type definition

                else
                    # something else, e.g. inline inner constructor
                    #   F(...) = ...
                    continue
                end

            end

        ################################################################
        # No aml name - No defVal
        else  # var/var::T  # ignored for creating aml

            # Type Checker
            if ei isa Symbol #  var
                # args_name[iArg] = missing # argument ignored for aml
                # args_function[iArg] =  missing # ignored for creating aml

                args_type[iArg] = String

                var = ei

                args_param[iArg] = var
                args_var[iArg] =  var

            elseif ei.head == :(::) && ei.args[1] isa Symbol # var::T
                # args_name[iArg] = missing # argument ignored for aml
                # args_function[iArg] =  missing # ignored for creating aml

                var = ei.args[1]
                varType = ei.args[2] # Type

                args_type[iArg] = varType
                args_param[iArg] = var
                args_var[iArg] = var

            elseif ei.head == :block  # anything else should be evaluated again
                # can arise with use of @static inside type decl
                argsexpr, args_param, args_defaultvalue, args_type, args_var, args_name, args_function, args_literaltype, struct_name, struct_nodetype, struct_function, is_struct_mutable, args_custom_creator, args_custom_extractor, args_custom_updater, T = aml_parse(expr)
            else
                continue
            end

        end # end ifs
    end # endfor

    ########################
    # self closing tags checker
    if  struct_nodetype == AbsEmpty
        # add a field with nothing type
        push!(args_name, "content") # argument ignored for aml
        push!(args_type, Nothing)
        push!(args_function,missing)
        push!(args_literaltype, AbsIgnore)
        push!(args_param, Expr(:kw, :content, nothing))
        push!(args_var, :content)
        push!(args_defaultvalue, nothing)
        push!(argsexpr.args,:(content::Nothing))
        push!(args_custom_creator, nothing)
        push!(args_custom_extractor, nothing)
        push!(args_custom_updater, nothing)

        # argsexpr, args_param, args_defaultvalue, args_type, args_var, args_name, args_function, args_literaltype, struct_name, struct_nodetype, struct_function, is_struct_mutable, args_custom_creator, args_custom_extractor, args_custom_updater, T = aml_parse(expr)
    end

    ########################
    # aml::Node adder
    push!(argsexpr.args,:(aml::Union{Document,Node}))

    return argsexpr, args_param, args_defaultvalue, args_type, args_var, args_name, args_function, args_literaltype, struct_name, struct_nodetype, struct_function, is_struct_mutable, args_custom_creator, args_custom_extractor, args_custom_updater, T
end


text:
using .Tracker: @grad, TrackedArray, track, data

TransmutedDimsArray{T,N,P,Q,AT}(A::TrackedArray) where {T,N,P,Q,AT<:AbstractArray} = track(TransmutedDimsArray, A, P)

@grad function TransmutedDimsArray(A::AbstractArray, P)
    B = TransmutedDimsArray(data(A), P)
    Q = getinvperm(B)
    function transmute_back(Δ::AbstractArray)
        if unique_or_zero(P)
            (transmute(data(Δ), Q), nothing)
        else
            (_undiagonal(data(Δ), P, Q), nothing)
        end
    end
    B, transmute_back
end


text:
function build_ode_problem(
  model::Model;
  parameters::Vector{Pair{Symbol,Float64}} = Pair{Symbol,Float64}[],
  events_active::Union{Nothing, Vector{Pair{Symbol,Bool}}} = Pair{Symbol,Bool}[],
  events_save::Union{Tuple,Vector{Pair{Symbol, Tuple{Bool, Bool}}}} = (true,true),
  observables_::Union{Nothing,Vector{Symbol}} = nothing,
  saveat::Union{Nothing,AbstractVector} = nothing,
  tspan::Union{Nothing,Tuple} = nothing,
  save_scope::Bool = true,
  time_type::DataType = Float64
)
  # saveat and tspan
  _saveat, _tspan = saveat_tspan(saveat, tspan, time_type)

  # initial values and params
  init_func = model.init_func
  parameters_nt = NamedTuple(parameters)
  u0, params = init_values(init_func, merge(model.constants, parameters_nt))

  # check observables
  if !isnothing(observables_)
    records_ind = indexin(observables_, records(model))
    if any((x)-> x === nothing, records_ind)
      lost_observables_ind = findall((x)-> x===nothing, records_ind)
      lost_observables = observables_[lost_observables_ind]
      error("The following observables have not been found in the model: $lost_observables")
    end
  end

  # saving setup
  utype = eltype(u0)
  merged_observables = isnothing(observables_) ? observables(model) : observables_ # use default if not set
  saved_values = SavedValues(
    LArray{utype,1,Array{utype,1},Tuple(merged_observables)}[],
    time_type[],
    save_scope ? Symbol[] : nothing
    )
  saving_func = model.saving_generator(merged_observables)
  scb = saving_wrapper(saving_func, saved_values; saveat=_saveat, save_scope)

  # events
  ev_on_nt = !isnothing(events_active) ? NamedTuple(events_active) : NamedTuple()
  events = active_events(model.events, merge(model.events_active, ev_on_nt), events_save)
  cbs = CallbackSet(scb, events...)

  # problem setup
  return ODEProblem(
    model.ode_func, # ODE function
    u0, # u0
    _tspan, # tspan
    params; # constants and static
    callback = cbs # callback
  )
end

function saveat_tspan(saveat, tspan, time_type)

  if !isnothing(saveat) && !isempty(saveat)
    _saveat = collect_saveat(saveat)
    _tspan = (zero(time_type), time_type(maximum(_saveat)))
  elseif !isnothing(tspan)
    _saveat = time_type[]
    _tspan = (zero(time_type), time_type(last(tspan))) # tspan should start from zero?
  else
    error("Please, provide either `saveat` or `tspan` value.")
  end
  return (_saveat, _tspan)
end

collect_saveat(saveat::Tuple) = Float64[]
collect_saveat(saveat::Vector{S}) where S<:Real = Float64.(saveat)
collect_saveat(saveat::AbstractRange{S}) where S<:Real = Float64.(saveat)

function init_values(init_func, constants)
  u0, p0 = init_func(constants)
  lvcons = LVector(constants)
  return (u0,Params{typeof(lvcons),typeof(p0)}(lvcons, p0))
end

function  update_init_values(prob, init_func, x)
  constants = merge(NamedTuple(prob.p.constants),x)
  u0, p0 = init_values(init_func, constants)

  prob_upd = remake(prob; u0=u0, p=p0)
  #prob_upd.p.constants .= p0.constants
  #prob_upd.p.static .= p0.static

  return prob_upd
end

function remake_saveat(prob, measurements)
  saveat = unique([dp.t for dp in measurements])
  tspan = (prob.tspan[1], eltype(prob.tspan)(maximum(saveat)))
  scb_orig = prob.kwargs[:callback].discrete_callbacks[1].affect!
  utype = eltype(prob.u0)
  save_scope = scb_orig.save_scope

  saved_values = SavedValues(
    LArray{utype,1,Array{utype,1},observables(prob)}[],
    eltype(tspan)[],
    save_scope ? Symbol[] : nothing
    )
  save_func = scb_orig.save_func
  scb_new = saving_wrapper(save_func, saved_values; saveat, save_scope=save_scope)

  cbs = list_callbacks(prob)
  cb_set = CallbackSet(scb_new,cbs[1]...,cbs[2]...)
  remake(prob; callback=cb_set, tspan=tspan)
end

function list_callbacks(prob)
  discr_cbs = prob.kwargs[:callback].discrete_callbacks
  cont_cbs = prob.kwargs[:callback].continuous_callbacks
  (discr_cbs[1:end .!=1], cont_cbs)
end

observables(prob::SciMLBase.AbstractODEProblem) = observables(
  eltype(prob.kwargs[:callback].discrete_callbacks[1].affect!.saved_values.u)
)

observables(::Type{LArray{T,N,D,Syms}}) where {T,N,D,Syms} = Syms


text:
struct MatrixBingham{N,T,n,k} <: ParameterizedMeasure{N}
    par::NamedTuple{N,T}
    manifold::Stiefel{n,k,ℝ}
end
function MatrixBingham(n, k; kwargs...)
    par = NamedTuple(kwargs)
    return MatrixBingham(par, Stiefel(n, k))
end

Manifolds.base_manifold(d::MatrixBingham) = getfield(d, :manifold)

function MeasureTheory.basemeasure(μ::MatrixBingham{(:B,)})
    return normalize(Hausdorff(base_manifold(μ), μ.B))
end

function MeasureTheory.logdensity(d::MatrixBingham{(:B,)}, x::AbstractMatrix)
    B = d.B
    return dot(x, B, x) - log(pFq((k//2,), (n//2,), B))
end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

using Pkg.Artifacts

datadir = joinpath(artifact"acedata_v0.8.1", "ACEData-0.8.1")
testsdir = joinpath(datadir, "tests")
traindir = joinpath(datadir, "trainingsets")
potsdir = joinpath(datadir, "potentials")

# unzip the eam potential
if !isfile(joinpath(traindir, "w_eam4.fs"))
   eamzip = joinpath(traindir, "w_eam4.fs.zip")
   run(`unzip $eamzip -d $traindir`)
end


text:
"""
$(TYPEDSIGNATURES)

Plot the a 3D surface from the Beta mix together measured data as circles.
"""
function plotsurface3D(dobs,betamix ; yscal=1e2, zscal=2.0, markersize=3500,
                       displayfig=true, ymin=nothing,ymax=nothing,
                       savefig=false,outdir="")

    # markersize=350Makie.px

    println("\nPlotting 3D surface from Beta mix and measured data as circles.")
    println(" Scaling factors are $yscal for [$(dobs.protein)] and $zscal for enthalpy.")

    N = 300
    xmin = betamix.betpar.a
    xmax = betamix.betpar.b
    if ymin==nothing
        ymin = betamix.betpar.ymin
    end
    if ymax==nothing
        ymax = betamix.betpar.ymax
    end

    xv = collect(LinRange(xmin,xmax,N))
    yv = collect(LinRange(ymin,ymax,N))

    nx = length(xv)
    ny = length(yv)
    dcalc = zeros(nx,ny)
    for j=1:ny
        xy = [xv yv[j]*ones(nx)]
        dcalc[:,j] = forwmod2D(betamix.betpar,xy,betamix.modkonamp)
    end
    zmin = minimum(dcalc)
    zmax = maximum(dcalc)

    xobs = dobs.sdsprotcon[:,1]
    yobs = dobs.sdsprotcon[:,2]
    zobs = dobs.enthalpy

    ## Makie plot
    fig = GLMakie.Figure(resolution = (1000, 800))
    ax1 = GLMakie.Axis3(fig[1,1]) # Axis3 -> 3D

    GLMakie.surface!(ax1,xv,yscal*yv,zscal*dcalc,alpha=1.0,colormap=:rainbow1)

    GLMakie.scatter!(ax1,xobs,yscal*yobs,zscal*zobs,color=:black,markersize=markersize)

    ax1.xlabel = "[SDS]"
    ax1.ylabel = "[$(dobs.protein)]x$yscal"
    ax1.zlabel = "Enthalpy x $zscal"

    if savefig
        flname = joinpath(outdir,"surfaceplot.png")
        GLMakie.save(flname,fig)
    end

    GLMakie.center!(fig.scene)
    if displayfig
        display(fig)
    end

    return fig
end