text:
"""
A generic ERGM for directed networks to be estimated using the pseudo likelihood
"""
struct  NetModeErgmPml <: Ergm
    ergmTermsString::String # needs to be compatible with R ergm package notation and names
    isDirected::Bool
    nErgmPar::Int
end
export NetModeErgmPml

NetModeErgmPml() = NetModeErgmPml("", false, 0)

NetModeErgmPml(ergmTermsString, isDirected)  = NetModeErgmPml(ergmTermsString, isDirected, 1 + count(i->(i=='+'), ergmTermsString))

name(x::NetModeErgmPml) = "NetModeErgmPml($(x.ergmTermsString))"

type_of_obs(model::NetModeErgmPml) =  Array{Float64, 2}

sample_ergm(model::NetModeErgmPml, N, parVec, nSample) = sample_ergm_RCall_sequence(model.ergmTermsString, N, parVec, nSample)[1][:,:,1,:]

sample_ergm_sequence(model::NetModeErgmPml, N, parVecSeq_T::Matrix, nSample) = sample_ergm_RCall_sequence(model.ergmTermsString, N, parVecSeq_T, nSample)

function change_stats(model::NetModeErgmPml, A::Matrix)
    ErgmRcall.get_change_stats(A,model.ergmTermsString)
end

function pseudo_loglikelihood_strauss_ikeda(par, changeStat, response, weights)
    # @debug "[pseudo_loglikelihood_strauss_ikeda][start][size par=$(size(par)), size changeStat = $(size(changeStat)), size response=$(size(response)), size weights=$(size(weights))]"

    H_vec = dropdims(sum(par.*changeStat', dims=1), dims=1)

    logPTot = sum(- weights .* ( ((1 .-response) .* H_vec) .+ log.(1 .+ exp.(-H_vec))) )

    #old version should be equal to current one
    # logit_P = sum(par.*changeStat', dims=1)

    # P = inv_logit.(logit_P)

    # Pbar = 1 .- P

    # logPVec = log.([response[i] == zero(response[i]) ? Pbar[i] : P[i] for i=1:length(response) ])

    # logPTot = sum(logPVec.*weights)
    # @debug "[pseudo_loglikelihood_strauss_ikeda][end]"

    return  logPTot
end

function grad_pseudo_loglikelihood_strauss_ikeda(par, changeStat, response, weights)

    H_vec = dropdims(sum(par.*changeStat', dims=1), dims=1)

    grad =  dropdims(sum((changeStat.*weights .* ( response .- (1 .+ exp.(-H_vec)).^(-1) ) ), dims=1), dims=1)

    return grad
end

function hessian_pseudo_loglikelihood_strauss_ikeda(par, changeStat, weights)

    H_vec = dropdims(sum(par.*changeStat', dims=1), dims=1)

    P_vec = (1 .+ exp.(-H_vec)).^(-1)

    nCSGroups, nErgmPar = size(changeStat) # number of groups of different change statistics
    prod = P_vec .* (1 .- P_vec).*weights
    hess = - [sum(prod.*changeStat[:,i].*changeStat[:,j]) for i = 1:nErgmPar, j=1:nErgmPar]

    return hess
end

function fisher_info_pseudo_loglikelihood_strauss_ikeda(par, changeStat,  weights)

    - hessian_pseudo_loglikelihood_strauss_ikeda(par, changeStat, weights)
end

obj_fun(par, changeStat, response, weights) =  pseudo_loglikelihood_strauss_ikeda(par, changeStat, response, weights)

grad_obj_fun(par, changeStat, response, weights) = grad_pseudo_loglikelihood_strauss_ikeda(par, changeStat, response, weights)

hessian_obj_fun(par, changeStat, weights) = hessian_pseudo_loglikelihood_strauss_ikeda(par, changeStat, weights)

fisher_info_obj_fun(par, changeStat,  weights) = fisher_info_pseudo_loglikelihood_strauss_ikeda(par, changeStat,  weights)

estimate_sequence(model::NetModeErgmPml, obsT::Array{Array{Float64, 2}}) = reduce(hcat, [ErgmRcall.get_one_mple(obst, model.ergmTermsString) for obst in obsT])

estimate_sequence(model::NetModeErgmPml, obsT::Array{T, 3} where T <: Real) = reduce(hcat, [ErgmRcall.get_one_mple(obsT[:,:,t], model.ergmTermsString) for t in 1:size(obsT)[3]])


text:
@testset "stepper.jl" begin
    @testset "Iterator sampling" begin
        Random.seed!(1234)
        as = []
        bs = []

        iter = AbstractMCMC.steps(MyModel(), MySampler())
        iter = AbstractMCMC.steps(MyModel(), MySampler(); a = 1.0) # `a` shouldn't do anything

        for (count, t) in enumerate(iter)
            if count >= 1000
                break
            end

            # don't save missing values
            t.a === missing && continue

            push!(as, t.a)
            push!(bs, t.b)
        end

        @test length(as) == length(bs) == 998

        @test mean(as) ≈ 0.5 atol=2e-2
        @test var(as) ≈ 1 / 12 atol=5e-3
        @test mean(bs) ≈ 0.0 atol=5e-2
        @test var(bs) ≈ 1 atol=5e-2

        @test Base.IteratorSize(iter) == Base.IsInfinite()
        @test Base.IteratorEltype(iter) == Base.EltypeUnknown()
    end
end

ext: jl


text:
# # Example 01: Relaxation experiment

# This script provides an example of fitting of a relaxation experiment.
# First import the packages.

using RHEOSafm
using RHEOS

using Plots
using Plots.PlotMeasures
gr();

# Next, the output file provided by the JPK software is imported. The file contains the time-force-displacement values and they are automatically converted into time-stress-strain.
# The function [`importJPK`](@ref) requires the file path, the radious of the indenter, and the segments of the curve to upload (i.g. extend, pause, retract).
input_file = joinpath(@__DIR__, "assets", "AFM_relaxation.txt")
R = 150e-9; # Radious of the indenter
interface = AFM(R);
data = importJPK(input_file, interface, sections = ["extend", "pause"]);

#md # !!! note "Beware!"
#md #     Currently RHEOSafm makes use of the Hertz contact model to convert force-displacement to stress-strain. Additional tip geometries will be added in next updates.
plot(data.t, data.σ, legend = false, xlabel = "Time", ylabel = "Stress", guidefont= 10, size = (300,300), lw=3, margin=16mm)

# To detect the point at which approximately contact occurs using RHEOSafm it is possible to:
# 1) define a force threshold
# 2) apply Hertz spherical contact model
# In this example the "threshold" method is used. An application of the Hertz method is availabel in example 02.
data_contact = contact_point(data, interface, "threshold", (threshold = 1e-8,));

#md # !!! compat "Note"
#md #     From this point, RHEOS functionalities are used. For more information refer to [RHEOS documentation](https://juliarheology.github.io/RHEOS.jl/stable/).

# To speed up the fitting procedure, the data points are downsampled.
d_downsample = resample(data_contact, -20);
# The relaxation curve is then fitted using a Fractional Standard Linear Solid model
SLS_model = modelstepfit(d_downsample, FractSLS_Zener, strain_imposed);
# Now we can extract the strain pattern
SLS_predict = extract(data_contact, strain_only);
# and calculate the stress based on the fitted model
SLS_predict = modelpredict(SLS_predict, SLS_model);
# Now we can plot data and model together for comparison
plot(data_contact.t, data_contact.σ, legend = true, xlabel = "Time", ylabel = "Stress", label = "Experimental data", guidefont= 10, size = (500,500), lw=3, margin=16mm)
plot!(SLS_predict.t, SLS_predict.σ, label = "Fitted",guidefont= 10, size = (400,400), lw=2)


text:
"""
    backest_sum!(MLM; isXSum, isZSum)

Back-estimate left-out sum contrast interactions for an MLM object.

# Arguments

- MLM = Mlm object

# Keyword arguments

- isXSum = boolean flag indicating whether a left-out `X` sum contrast
  interaction needs to be back-estimated. Defaults to `true`.
- isZSum = boolean flag indicating whether a left-out `Z` sum contrast
  interaction needs to be back-estimated. Defaults to `true`.

# Value

None; updates coefficient estimates in Mlm object in place.

# Some notes

This implementation assumes that the `X` predictor matrix
(when `isXSum` is set to `true`) and the `Z` predictor matrix
(when `isZSum` is set to `true`) consist of the sum contrasts for a single
categorical variable (and no other variables).

"""
function backest_sum!(MLM::Mlm; isXSum::Bool=true, isZSum::Bool=true)

    # Check to make sure X and/or Z intercepts were included for
    # interactions that need to be back-estimated
    if (isXSum==true) && (MLM.data.predictors.isXIntercept==false)
        error("No X intercept.")
    end
    if (isZSum==true) && (MLM.data.predictors.isZIntercept==false)
        error("No Z intercept.")
    end

    # Storing some larger computations to be used to calculate variances
    # Calculate and store transpose(Z)*Z
    ZTZ = transpose(get_Z(MLM.data))*get_Z(MLM.data)
    varLeft = inv(transpose(get_X(MLM.data))*get_X(MLM.data))
    varRight = transpose(ZTZ\(transpose((ZTZ\(transpose(get_Z(MLM.data))*
	                                     MLM.sigma))*get_Z(MLM.data))))

    # Back-estimate left-out X sum contrast interaction
    if isXSum==true
        C = transpose(vcat(0, -ones(MLM.data.p-1)))
        coeffX = C*MLM.B
        varX = transpose(kron_diag(C*varLeft*transpose(C), varRight))

        if isZSum==false
            MLM.B = vcat(MLM.B, coeffX)
            MLM.varB = vcat(MLM.varB, varX)
        end
    end

    # Back-estimate left-out Z sum contrast interaction
    if isZSum==true
        D = vcat(0, -ones(MLM.data.q-1))
        coeffZ = MLM.B*D
        varZ = transpose(kron_diag(varLeft, transpose(D)*varRight*D))

        if isXSum==false
            MLM.B = hcat(MLM.B, coeffZ)
            MLM.varB = hcat(MLM.varB, varZ)
        end
    end

    # Back-estimate left-out X and Z sum contrast interaction
    if (isXSum==true) && (isZSum==true)
        coeffXZ = C*MLM.B*D
        varXZ = transpose(kron_diag(C*varLeft*transpose(C),
                          transpose(D)*varRight*D))

        MLM.B = hcat(vcat(MLM.B, coeffX), vcat(coeffZ, coeffXZ))
        MLM.varB = hcat(vcat(MLM.varB, varX), vcat(varZ, varXZ))
    end
end

"""
    mlm_backest_sum(data; isXIntercept, isZIntercept, weights, targetType,
                    isXSum, isZSum)

Matrix linear model using least squares method that optionally back-estimates
left-out sum contrast interactions.

# Arguments

- data = RawData object

# Keyword arguments

- isXIntercept = boolean flag indicating whether or not to include an `X`
  intercept (row main effects). Defaults to `true`.
- isZIntercept = boolean flag indicating whether or not to include a `Z`
  intercept (column main effects). Defaults to `true`.
- weights = 1d array of floats to use as column weights for `Y`, or `nothing`.
  If the former, must be the same length as the number of columns of `Y`.
  Defaults to `nothing`.
- targetType = string indicating the target type toward which to shrink the
  error variance, or `nothing`. If the former, acceptable inputs are "A", "B",
  "C", and "D". Defaults to `nothing`.
    - "A": Target is identity matrix
    - "B": Target is diagonal matrix with constant diagonal
    - "C": Target is has same diagonal element, and same off-diagonal element
    - "D": Target is diagonal matrix with unequal entries
- isXSum = boolean flag indicating whether a left-out `X` sum contrast
  interaction needs to be back-estimated. Defaults to `true`.
- isZSum = boolean flag indicating whether a left-out `Z` sum contrast
  interaction needs to be back-estimated. Defaults to `true`.

# Value

None; updates coefficient estimates in Mlm object in place.

# Some notes

This implementation assumes that the `X` predictor matrix
(when `isXSum` is set to `true`) and the `Z` predictor matrix
(when `isZSum` is set to `true`) consisist of the sum contrasts for a single
categorical variable.

"""
function mlm_backest_sum(data::RawData;
                         isXIntercept::Bool=true, isZIntercept::Bool=true,
                         weights=nothing, targetType=nothing,
                         isXSum::Bool=true, isZSum::Bool=true)

    # Run matrix linear models
    MLM = mlm(data; isXIntercept=isXIntercept, isZIntercept=isZIntercept,
              weights=weights, targetType=targetType)

    # Back-estimate left-out sum contrast interactions
    backest_sum!(MLM; isXSum=isXSum, isZSum=isZSum)

    return MLM
end


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/05/2014

module WorldModel

export
    AbstractWorldModel,

    initialize,
    update,

    updateObjectState,
    getAllObjectStates,
    getAll,
    updateAll,

    AirSpace,
    ASWMState

import CommonInterfaces: initialize, update

using AbstractWorldModelImpl

using AirSpaceImpl

end


text:
#=
Author: Stankevich Andrey, MIPT, Russia, 2020

This package provides function templates for
exact solutions of the 2D transfer equation
Cauchy problem
See https://keldysh.ru/papers/2002/prep4/prep2002_4.html
=#

# helper function
inflate(f, xs, ys) = [f(x, y) for x in xs, y in ys]

####################
# radius vector norm

function r_template(x, y, x_00, y_00)
	return ((x .- x_00) .^ 2 + (y .- y_00) .^ 2) .^ 0.5
end

#################
# center location

function x0_constant_template(t, x_00)
	return x_00 * ones(size(t))
end

function dx0_dt_constant(t)
	return zeros(size(t))
end

function y0_constant_template(t, y_00)
	return y_00 * ones(size(t))
end

function dy0_dt_constant(t)
	return zeros(size(t))
end

function x0_linear_template(t, u_00, x_00)
	return x_00 .+ u_00 * t
end

function dx0_dt_linear_template(t, u_00)
	return u_00 * ones(size(t))
end

function y0_linear_template(t, u_00, y_00)
	return y_00 .+ u_00 * t
end

function dy0_dt_linear_template(t, u_00)
	return u_00 * ones(size(t))
end

function x0_circle_template(t, r_00, T_0)
	return r_00 * cos(2.0 * pi * t / T_0)
end

function dx0_dt_circle_template(t, r_00, T_0)
	return (2.0 * pi / T_0) * r_00 * sin(2.0 * pi * t / T_0) * (-1.0)
end

function y0_circle_template(t, r_00, T_0)
	return r_00 * sin(2.0 * pi * t / T_0)
end

function dy0_dt_circle_template(t, r_00, T_0)
	return (2.0 * pi / T_0) * r_00 * cos(2.0 * pi * t / T_0)
end

############################################
# parametric function a2 for transfer speeds

function a2_constant_template(t, T_1)
	return 2.0 * π * t / T_1
end

function da2_dt_constant_template(t, T_1)
	return 2.0 * π / T_1
end

function da2_increasing_template(t, T_1, c)
	return 2.0 * π * (t + c * t^2) / (T_1 + c * T_1^2)
end

function da2_dt_increasing_template(t, T_1, c)
	return 2.0 * π * (2 * c * t) /(T_1 + c * T_1^2)
end

function da2_decaying_template(t, T_1, T_2)
	return 2.0 * π * (1 - exp(-t/T_1)) / (1 - exp(-T_1/T_2))
end

function da2_dt_decaying_template(t, T_1, T_2)
	return 2.0 * π * exp(-t/T_1) / (1 - exp(-T_1/T_2))
end

############################################
# parametric function a1 for transfer speeds

function a1_linear(radius)
	return radius
end

function da1_dt_linear(radius)
	return radius .* 0
end

####################
# Initial conditions

function cone_template(r, phi, r_1, H_0)
	return H_0 * clamp.(1.0 .- r / r_1, 0.0, 1.0)
end

function cone_slice_template(r, phi, r_1, H_0, phi_0, phi_1)

	result = H_0 * clamp.(1.0 .- r / r_1, 0.0, 1.0)
	mask   = map(&, phi .<= phi_1, phi .>= phi_0)

	result[mask] .= 0.0

	return result
end

#############################
# transfer velocity functions
# should have a specific form

function v1_template(x, y, t, a1, da2_dt, y0, dx0_dt, r1)

	radius = r(x, y)

	lh = a1(radius) * (y .- y0(t)) / radius * da2_dt(t) + dx0_dt(t)

	return lh * (radius .<= r1)
end

function v2_template(x, y, t, a1, da1_dt, x0, dy0_dt, r1)

	radius = r(x, y)
	lh = -a1(radius) * (x .- x0(t)) / radius * da1_dt(t) + dy0_dt(t)

	return lh * (radius .<= r1)
end

#############################
# solution

import Base.Iterators.product

function solution_template(xs, ys, ts, R, a1, a2, x0, y0)

	grid = [(x, y) for x in xs, y in ys]

	result = []

	for t in ts

		x0_coord = x0(t)
		y0_coord = y0(t)

		radius = inflate(
			(x, y) -> ((x .- x0_coord) ^ 2 + (y .- y0_coord) ^ 2) ^ 0.5,
			xs, ys
		)

		rlh = [
			r != 0.0 ?
			acos((x .- x0_coord) / r) :
			0.5 * pi
			for (r, (x, y)) in zip(radius, grid)
		]

		rrh = [
			r != 0.0 ?
			(2.0 * pi .- acos((x .- x0_coord) / r)) :
			1.5 * pi
			for  (r, (x, y)) in zip(radius, grid)
		]

		rh   = zeros((length(xs), length(ys)))
		mask    = [y .<  y0_coord for (x, y) in grid]
		notmask	= [y .>= y0_coord for (x, y) in grid]
		rh[mask]    = rrh[mask]
		rh[notmask] = rlh[notmask]

		phi = [r != 0.0 ? a1(r) * a2(t) / r : 0.0 for r in radius]
		phi = phi + rh

		periods = floor.(phi / (2.0 * pi))
		phi     = phi - periods * 2.0 * pi

		push!(result, R(radius, phi))

	end

	return result

end

max_stars_repo_name: HSE-LAMBDA/RheologyReconstruction

max_issues_repo_name: HSE-LAMBDA/RheologyReconstruction

max_forks_repo_name: HSE-LAMBDA/RheologyReconstruction


text:
export BFFPS19

# ===============================================================
# Bogomolov, Forets, Frehse, Potomkin, Schilling.
# ===============================================================

function options_BFFPS19()
    return OptionSpec[
        # general options
        OptionSpec(:discretization, "forward", domain=String,
            info="model for bloating/continuous time analysis"),
        OptionSpec(:δ, 1e-2, domain=Float64, aliases=[:sampling_time],
            domain_check=(v  ->  v > 0.), info="time step"),
        OptionSpec(:vars, Int[], domain=AbstractVector{Int}, domain_check=(
            v  ->  length(v) > 0 && all(e -> e > 0, v)),
            info="variables of interest; default: all variables"),
        OptionSpec(:guards_proj, Vector(), domain=AbstractVector,
             info="internal storage of projected guard constraints"),
        OptionSpec(:partition, [Int[]],
            domain=AbstractVector{<:AbstractVector{Int}}, domain_check=
            ispartition,
            info="block partition; a block is represented by a vector " *
                 "containing its indices"),

        # discretization options
        OptionSpec(:sih_method, "concrete", domain=String,
            info="method to compute the symmetric interval hull in discretization"),
        OptionSpec(:exp_method, "base", domain=String,
            info="method to compute the matrix exponential"),
        OptionSpec(:assume_sparse, false, domain=Bool,
            info="use an analysis for sparse discretized matrices?"),

        # reachability options
        OptionSpec(:lazy_inputs_interval, lazy_inputs_interval_always,
            domain=Union{Int, Function},
            domain_check=(v  ->  !(v isa Int) || v >= -1),
            info="length of interval in which the inputs are handled as a " *
                 "lazy set (``-1`` for 'never'); may generally also be a " *
                 "predicate over indices; the default corresponds to ``-1``"),

        # approximation options
        OptionSpec(:block_options, nothing, domain=Any,
            info="short hand to set ':block_options_init' and " *
                 "':block_options_iter'"),
        OptionSpec(:block_options_init, nothing, domain=Any,
            info="option for the approximation of the initial states " *
                 "(during decomposition)"),
        OptionSpec(:block_options_iter, nothing, domain=Any,
            info="option for the approximation of the states ``X_k``, ``k>0``"),

        # convenience options
        OptionSpec(:assume_homogeneous, false, domain=Bool,
            info="ignore dynamic inputs during the analysis?"),
    ]
end

function normalization_BFFPS19!(𝑂::TwoLayerOptions)
    # :lazy_inputs_interval option: convert integers to functions
    if haskey_specified(𝑂, :lazy_inputs_interval)
        v = 𝑂[:lazy_inputs_interval]
        if v isa Int
            if v == -1
                𝑂.specified[:lazy_inputs_interval] = lazy_inputs_interval_never
            elseif v == 0
                𝑂.specified[:lazy_inputs_interval] = lazy_inputs_interval_always
            else
                𝑂.specified[:lazy_inputs_interval] = (k -> k % v == 0)
            end
        end
    end

    # :block_options options: use convenience option for '_init' and '_iter'
    if !haskey_specified(𝑂, :block_options_init) &&
            haskey_specified(𝑂, :block_options)
        𝑂.specified[:block_options_init] = 𝑂[:block_options]
    end
    if !haskey_specified(𝑂, :block_options_iter) &&
            haskey_specified(𝑂, :block_options)
        𝑂.specified[:block_options_iter] = 𝑂[:block_options]
    end

    nothing
end

function validation_BFFPS19(𝑂)
    # block_options
    for b_options in [:block_options, :block_options_init, :block_options_iter]
        if haskey_specified(𝑂, b_options)
            bo = 𝑂[b_options]
            if bo isa Type{<:LazySet} || bo isa Type{<:AbstractDirections}
                # uniform options
            elseif bo isa AbstractVector || bo isa Dict{Int, Any}
                # mapping
            elseif bo isa Real || bo isa Pair{<:UnionAll, <:Real}
                ε = bo isa Real ? bo : bo[2]
                if ε <= 0
                    throw(DomainError(ε, "the `$b_options` option must be " *
                                         "positive"))
                end
            elseif b_options == :block_options_iter && bo == nothing
                # no overapproximation
            else
                throw(DomainError(bo == nothing ? "nothing" : bo,
                    "the `$b_options` option does not accept the given input"))
            end
        end
    end

    nothing
end

"""
    BFFPS19 <: ContinuousPost

Implementation of the reachability algorithm for purely continuous linear
time-invariant systems using block decompositions by S. Bogomolov, M. Forets,
G. Frehse, A. Podelski, C. Schilling and F. Viry [1].

### Fields

- `options` -- an `Options` structure that holds the algorithm-specific options

### Notes

The following options are available:

```julia
$(print_option_spec(options_BFFPS19()))
```

### Algorithm

We refer to [1] for technical details.

[1] [Reachability analysis of linear hybrid systems via block
decomposition](https://arxiv.org/abs/1905.02458).
Sergiy Bogomolov, Marcelo Forets, Goran Frehse, Kostiantyn Potomkin, Christian Schilling.
"""
struct BFFPS19 <: ContinuousPost
    options::TwoLayerOptions

    function BFFPS19(𝑂::Options)
        normalized_𝑂 = validate_and_wrap_options(𝑂, options_BFFPS19();
            validation=validation_BFFPS19,
            normalization=normalization_BFFPS19!)
        return new(normalized_𝑂)
    end
end

# convenience constructor from pairs of symbols
BFFPS19(𝑂::Pair{Symbol,<:Any}...) = BFFPS19(Options(Dict{Symbol,Any}(𝑂)))

# default options
BFFPS19() = BFFPS19(Options())

include("reach.jl")
include("reach_blocks.jl")

init(𝒫::BFFPS19, 𝑆::AbstractSystem, 𝑂::Options) = init!(𝒫, 𝑆, copy(𝑂))

function init!(𝒫::BFFPS19, 𝑆::AbstractSystem, 𝑂::Options)
    # state dimension for (purely continuous or purely discrete systems)
    𝑂copy = copy(𝑂)
    𝑂copy[:n] = statedim(𝑆)

    # solver-specific options (adds default values for unspecified options)
    𝑂validated = validate_solver_options_and_add_default_values!(𝑂copy)

    # :vars option; default: all variables
    if haskey_specified(𝒫.options, :vars)
        𝑂validated[:vars] = 𝒫.options[:vars]
    else
        𝑂validated[:vars] = 1:𝑂validated[:n]
    end

    # :partition option: use 1D blocks
    if haskey_specified(𝒫.options, :partition)
        𝑂validated[:partition] = 𝒫.options[:partition]
    else
        𝑂validated[:partition] = [[i] for i in 1:𝑂validated[:n]]
    end

    opD = 𝒫.options[:opD]
    @assert opD isa DecomposedDiscretePost "this continuous post operator works " *
                                       "only with DecomposedDiscretePost"
    HS = 𝒫.options[:HS]
    constrained_dims = constrained_dimensions(HS)
    out_vars = opD.options[:out_vars]
    loc_id = 𝒫.options[:loc_id]
    temp_vars = unique([out_vars; constrained_dims[loc_id]])
    temp_vars = get_variables_from_relevant_blocks(𝑂validated[:partition], temp_vars)
    opD.options[:temp_vars] = temp_vars
    guards_constraints = [guard(HS, trans) for trans in out_transitions(HS, loc_id)]
    𝑂validated[:vars] = temp_vars
    𝑂validated[:guards_proj] = [project(c, temp_vars) for c in guards_constraints]
    𝑂validated[:blocks] = compute_blocks(𝑂validated[:vars], 𝑂validated[:partition])

    # :block_options_init & :block_options_iter options:
    # set default according to :partition
    if !haskey_specified(𝒫.options, :block_options_init)
        𝑂validated[:block_options_init] =
            compute_default_block_options(𝑂validated[:partition])
    end
    if !haskey_specified(𝒫.options, :block_options_iter)
        𝑂validated[:block_options_iter] =
            compute_default_block_options(𝑂validated[:partition])
    end

    if 𝑂validated[:project_reachset]
        𝑂validated[:output_function] = nothing
    else
        𝑂validated[:output_function] = 𝑂validated[:projection_matrix]
    end

    return 𝑂validated
end

"""
    post(𝒫::BFFPS19, 𝑆::AbstractSystem, 𝑂::Options)

Calculate the reachable states of the given initial value problem using `BFFPS19`.

### Input

- `𝒫` -- post operator of type `BFFPS19`
- `𝑆` -- sytem, initial value problem for a continuous ODE
- `𝑂` -- algorithm-specific options
"""
function post(𝒫::BFFPS19, 𝑆::AbstractSystem, 𝑂_input::Options)
    𝑂 = TwoLayerOptions(merge(𝑂_input, 𝒫.options.specified), 𝒫.options.defaults)

    @assert 𝑂[:mode] == "reach" "the mode $(𝑂[:mode]) is not supported"

    info("Reachable States Computation...")
    @timing begin
        Rsets = reach_mixed(𝑆, 𝑂)

        info("- Total")
    end

    # Projection
    if 𝑂[:project_reachset]
        info("Projection...")
        RsetsProj = @timing project(Rsets, 𝑂)
    else
        RsetsProj = Rsets
    end

    return ReachSolution(RsetsProj, 𝑂_input)
end

function get_variables_from_relevant_blocks(partition, vars)
    result = Vector{Int}()
    for ur in partition
        if any(v ∈ ur for v in vars)
            append!(result, collect(ur))
        end
    end
    return result
end


text:
export rule

@rule Bernoulli(:out, Marginalisation) (q_p::Beta, ) = begin
    rho_1 = clamp(exp(logmean(q_p)), tiny, huge)
    rho_2 = clamp(exp(mirroredlogmean(q_p)), tiny, huge)
    return Bernoulli(rho_1 / (rho_1 + rho_2))
end


text:
## Test implementing some math logic

# The other trig functions are handled by sympy.
# We find than Cos(Pi) is about 50 time faster than Sin(Pi), etc. in some tests
# So it would be good to implement more of this in Julia

# sympy does more at the moment.
@mkapprule Cos :nargs => 1

do_Cos(mx::Mxpr{:Cos},x) = cos_one_arg(mx,x)

# This is overwritten now anyway
# This does not work. But, If I eval it after starting SJulia, it does work
# Should dispatch on number of args, anyway.
# apprules(mx::Mxpr{:Sin}) = length(mx) == 1 ? sin_one_arg(mx,margs(mx)...) : mx

Cos_pi_coeff(mx::Mxpr{:Cos},c::T) where {T<:Integer} = iseven(c) ? 1 : -1

## Do multiples of pi/2, pi/3, and pi/4
function Cos_pi_coeff(mx::Mxpr{:Cos},c::Rational{T}) where T<:Integer
    n = c.num
    d = c.den
    d == 2 && return zero(d)
    if d == 3
        md = mod(n,6)
        if (md == 2 || md == 4)
            nmx = -1//2
        else
            nmx = 1//2
        end
        return nmx
    end
    if d == 4
        md = mod(n,8)
        if (md == 3 || md == 5)
            nmx = _moosq2
        else
            nmx = _oosq2
        end
        return nmx
    end
    return mx
end

Cos_pi_coeff(mx::Mxpr{:Cos},c::T) where {T<:AbstractFloat} = cospi(c)
Cos_pi_coeff(mx::Mxpr{:Cos},c) = mx

function cos_one_arg(mx::Mxpr{:Cos}, arg::Mxpr{:Times})
    if is_Complex(arg[1])  # find a good way to dispatch on type here.
        (r,i) = reim(arg[1])
        if r == 0  #  Cos(I x) -> Cosh(x)
            nargs = copy(margs(arg))
            nargs[1] = i
            return mxpr(:Cosh,mxpr(:Times,nargs))
        end
    end
    return length(arg) == 2 ? Cos_factor_arg(mx,arg.args...) :
    mx
end

# function cos_one_arg2{T<:Real}(mx::Mxpr{:Cos}, arg::Mxpr{:Times}, arg1::Complex{T})
# end
# function cos_one_arg2(mx::Mxpr{:Cos}, arg::Mxpr{:Times}, arg1)
# end

function Cos_factor_arg(mx::Mxpr{:Cos},f1::T,f2::SJSym) where T<:Number
    if f2 == :Pi
        return Cos_pi_coeff(mx,f1)
    else
        return mx
    end
end

Cos_factor_arg(mx::Mxpr{:Cos},f1,f2) = mx
@inline cos_one_arg(mx::Mxpr{:Cos},arg::Symbol) = arg == :Pi ? -1 : mx
@inline cos_one_arg(mx::Mxpr{:Cos},arg::T) where {T<:Integer} = arg == 0 ? 1 : mx
@inline cos_one_arg(mx::Mxpr{:Cos},x::T) where {T<:AbstractFloat} = cos(x)
@inline cos_one_arg(mx::Mxpr{:Cos},x::Complex{T}) where {T<:AbstractFloat} = cos(x)
@inline cos_one_arg(mx::Mxpr{:Cos},x::Mxpr{:ArcCos}) = length(x) == 1 ? x[1] : mx
function cos_one_arg(mx::Mxpr{:Cos},x::Mxpr{:ArcSin})
    res = mpow((1-mpow(x[1],2)),(1//2))
    return res
end
cos_one_arg(mx::Mxpr{:Cos},x) = mx

#### Sin

# This is overwritten somewhere. Does not work

sin_one_arg(mx::Mxpr{:Sin}, arg::T) where {T<:Integer} = arg == 0 ? 0 : mx

function sin_one_arg(mx::Mxpr{:Sin},arg::Mxpr{:Times})
    if is_Complex(arg[1])  # use dispatch!
        (r,i) = reim(arg[1])
        if r == 0  #  Sin(I x) -> I Sinh(x)
            nargs = copy(margs(arg))
            nargs[1] = i
            return mxpr(:Times, complex(0,1), mxpr(:Sinh,mxpr(:Times,nargs)))
        end
    end
#    todo
    mx
end

sin_one_arg(mx::Mxpr{:Sin},x) = mx

####


text:
########################################
## File Name: forward_backward_simulation_test.jl
## Author: Haruki Nishimura (hnishimura@stanford.edu)
## Date Created: 2020/05/18
## Description: Test code for src/forward_backward_simulation.jl
########################################

using LinearAlgebra
using Random

@testset "Forward-Backward Simulation Test" begin

# Localization Task
dtc = 0.01;
dto = 0.5;
dtexec = dto*0.8;
Q = Matrix(1.0I, 2, 2);
Cu = Matrix(1.0I, 2, 2);
transPos = TransModel_Pos();
ukfModel = UKFPosRange();
simPosUKF = SimulatePosRangeLocalization2D(transPos,ukfModel,dtc,dto,dtexec,Q,Cu);

rng = MersenneTwister(1234);
p = PhysPos(0.,zeros(2));
bVec = VecBelMvNormal([BelMvNormal(0.,rand(rng, 2),10*Matrix(1.0I, 2, 2)) for ii = 1:100]);
s = AugState(p,bVec);

UArray = [PosControl(round(dtc*t,digits=5),zeros(2)) for t = 0:Int64(2.0/dtc)-1];

Y,S,S_before_update = simulateForward(simPosUKF,s,UArray,MersenneTwister(1234));

@test length(S_before_update) == 4
@test isequal(S_before_update[1].t, 0.49)
@test isequal(S_before_update[1].p.pos, zeros(2))
@test isapprox(S_before_update[1].b[1].params, [0.590845, 0.766797, 10.0, 0.0, 0.0, 10.0], rtol=1e-5)
@test isequal(S_before_update[4].t, 1.99)
@test isequal(S_before_update[4].p.pos, zeros(2))
@test isapprox(S_before_update[4].b[1].params, [0.828148, 1.07501, 11.2299, -0.352906, -0.352906, 11.0388], rtol=1e-5)

CoSArray = simulateBackward(simPosUKF,S,S_before_update,UArray,Y)

@test length(CoSArray) == 201
@test CoSArray[1].t == 0.0
@test isapprox(CoSArray[1].p.pos, [752.687, 833.817], rtol=1e-5)
@test isapprox(CoSArray[1].b[1].params, [-11.2424, -15.1043, 8.97812, -0.435698, -0.782821, 8.91428], rtol=1e-5)

# Manipulation Task
dtc = 0.02;
dto = 0.2;
dtexec = dto*0.8;
QMP = [Matrix(1.0I, 6, 6) zeros(6,5); zeros(5,6) 0.00*Matrix(1.0I, 5, 5)];
RMP = Matrix(1.0I, 9, 9);
CsMP = QMP;
CuMP = Matrix(1.0I, 3, 3);
simMP = SimulateManipulate2D(dtc,dto,dtexec,QMP,RMP,CsMP,CuMP)

xMP = PhysManipulate2D(0.,[0.1,0.2],0.5,[0.3,0.4],0.6,3.,5.,[10.,11.],0.2);
uMP = MControl2D(0.,6.,7.,8.);
bMP = BelMvNormal(xMP.t,vec(xMP),diagm([10.,10.,1.,10.,10.,10.,1.,1.,10.,10.,0.1]));

UArray = [MControl2D(round(dtc*t,digits=5),0.,0.0,0.0) for t = 0:Int64(2.0/dtc)-1];

Y,S,S_before_update = simulateForward(simMP,bMP,UArray,MersenneTwister(13));

@test length(S_before_update) == 10
@test S_before_update[1].t == 0.18
@test isapprox(S_before_update[1].params[1:11], [0.153713, 0.271617, 0.608, 0.296419, 0.395226, 0.6, 3.0, 5.0, 10.0, 11.0, 0.2], rtol=1e-5)
@test S_before_update[10].t == 1.98
@test isapprox(S_before_update[10].params[1:11], [8.91097, -4.67647, 3.30309, 2.53044, -4.53245, 4.0016, 3.21738, 5.0, 10.3436, 14.5485, 0.0583465], rtol=1e-5)

CoSArray = simulateBackward(simMP,S,S_before_update,UArray,Y)

@test length(CoSArray) == 101
@test CoSArray[1].t == 0.0
@test isapprox(CoSArray[1].params[1:11], [-0.606155, 0.0748971, 170.315, 0.574117, -0.324966, -154.575, -0.469705, 0.0, -19.4836, 13.3016, -3.49951], rtol=1e-5)

end


text:
# make Toydata
"""
design_maxtrix(subspace,moduels::Group_Mod)

Generates the design matrix for the given subspace of constructs.
The Designmatrix has n x p dimension where n is the number of constructs and p the number of modules.
Depending on the application it needs to be transposed.
"""

function design_maxtrix(subspace,moduels::Group_Mod)
    # mod to index
    dict_mod = Dict(moduels.m .=> collect(1:length(moduels.m)))
    # setup design matrix
    dm = Array{Int,2}(undef,length(subspace),length(moduels))
    #make matrix
    for (row,construct) in enumerate(subspace)
        dm[row,:] = _design_vector(construct,dict_mod)
    end
    return dm
end
"""
_design_vector(construct,dict_mod)
internal function to genareted the design space
generates the degign vector one construct.
"""
function _design_vector(construct,dict_mod)
    word2vec = zeros(length(dict_mod))
    index = [dict_mod[mod] for mod in construct]
    word2vec[index] .+=1
    return word2vec
end

"""
go2lab(rng,subspace,Toy_data ; repeat = 3 ,sigma = 1)

Simulated the evaluation of constructs in the lab using the given Toy data set.
This function wraps the evaluation of multiple constructs. The mean and standard deviation are returned:

[`_eval_con`](@ref)
"""

function go2lab(rng,subspace,Toy_data ; repeat = 3 ,sigma = 1)
    μ = Vector{Float64}(undef,length(subspace))
    σ = Vector{Float64}(undef,length(subspace))
    for (i,construct) in enumerate(subspace)
         results =  _eval_con(rng,construct,Toy_data,sigma,repeat)
         μ[i] = Statistics.mean(results)
         σ[i] = Statistics.std(results)
    end
    return μ,σ
end
"""
ground_truth(construct,Toy_data)

Calculates the true value of a construct given the Toy Data set.
"""

ground_truth(construct,Toy_data) = sum([Toy_data[mod] for mod in construct])

"""
_eval_con(rng,construct,Toy_data,sigma,repeat)

Internal function to evaluated a single constructed.
Sigma sets the standard deviation of the normal distribution that sets that generates the activities.
Repeats is the number of time the same constructed is evaluated.

[`_eval_con`](@ref)
"""
function _eval_con(rng,construct,Toy_data,sigma,repeat)
  mean = ground_truth(construct,Toy_data)
  distr = Normal(mean,sigma)
  return  [rand(rng,distr) for rep in 1:repeat]
end

@doc raw"""

BO_Lin_wrapper(constructs,activities,sigma,moduels)
Wrapper fit a Bayesian linear regression model using the given constructs.

```math
β ∼ \mathcal{N}(0,1)
```
"""
function BO_Lin_wrapper(constructs,activities,sigma,moduels)
    #fall back to set prior
    n_moduels = length(moduels) |> Int
    mw, Λw = zeros(n_moduels), Diagonal(ones(n_moduels))
    Σ_noise = Diagonal(sigma)
    return BO_Lin_wrapper(constructs,activities,moduels,mw,Λw,Σ_noise)
end

"""
BO_Lin_wrapper(constructs,activities,moduels,mw,Λw,Σ_noise)
Wrapper fit a Bayesian linear regression model using the given constructs.
Custom prior can be used in this case, if n = number of modules.
Than mw = vector of mean values with length n
Λw is the precision matrix nxn
"""

function BO_Lin_wrapper(constructs,activities,moduels,mw,Λw,Σ_noise)
    #setup of the prior
    f = BayesianLinearRegressor(mw, Λw)
    #make design Matrix
    X_design = transpose(design_maxtrix(constructs,moduels))
    # setup the model with design matrix off set if requierd
    fX = f(X_design, Σ_noise)
    # posterior
    f_post = posterior(fX,activities)
    return f_post
end

"""
BO_Lin_wrapper(constructs,activities,moduels,mw,Λw,Σ_noise)
Wrapper fit a Bayesian linear regression model using the given constructs.
Allow the model to update the prior base on the posterior of the previous model.
"""

function BO_Lin_wrapper(f_post,constructs,activities,sigma,moduels)
    #fall back to set prior
    n_moduels = length(moduels) |> Int
    mw, Λw = f_post.mw, f_post.Λw
    Σ_noise = Diagonal(sigma)
    return BO_Lin_wrapper(constructs,activities,moduels,mw,Λw,Σ_noise)
end
"""
make_predictions(Unseen_space, mod)

Wrapper function to make predictions using the model

"""
make_predictions(Model_output,Unseen_space, mod) = mean.(marginals(Model_output(transpose(design_maxtrix(Unseen_space ,mod)),eps())))

"""
Linear_TS_sampling(f,moduels,constructs,n_samples,σ² = 10^-6)
Thomson sampler for linear models. Internally a design matrix is constructed.
For more information see documentation or [`thompson_sampling`](@ref)

"""

function Linear_TS_sampling(f,moduels,constructs,
                n_samples,σ² = 10^-6)
    @assert length(constructs) ≥ n_samples
    dm = transpose(design_maxtrix(constructs,moduels))
    # indices selected
    selected = Int[]
    all_test = copy(dm)
    Ndist = f(all_test, σ²)
    while length(selected) < n_samples
        f̂ = rand(Ndist)
        f̂[selected] .= - Inf
       i = argmax(f̂)
       push!(selected, i)
   end
   return constructs[selected]
end

"""
model_linear(design,n_first,n_samples,Toy_data,n_cycles)

Wrapper to execute multiple cycles of the Optimisation process efficiently.
A demonstration function only useful to multiple evaluated steps.
Function you probably don't need.

"""

function model_linear(design,n_first,n_samples,Toy_data,n_cycles)
    rng = MersenneTwister()
    # to save the output
    max_model = Vector{Float64}(undef,n_cycles)
    #first steps
    space = getspace(design)
    new_constructs = StatsBase.sample(rng,space,n_first)
    lab_μ,lab_σ = go2lab(rng,new_constructs,Toy_data)
    df = DataFrame(constructs = new_constructs,μ = lab_μ, σ = lab_σ);
    #get the full space explicitly
    Unseen_space = getspace(design,full = true).space
    for i in 1:n_cycles
        # make model
        Model_output = BO_Lin_wrapper(df.constructs,df.μ,df.σ,design.mod)
        # filter space
        Unseen_space = filter(x-> !(x in new_constructs) , Unseen_space)
        #sample with TS
        new_constructs = Linear_TS_sampling(Model_output,design.mod,Unseen_space,n_samples)
        #evaluated new points
        lab_μ,lab_σ = go2lab(rng,new_constructs,Toy_data);
        new_df = DataFrame(constructs = new_constructs, μ = lab_μ, σ = lab_σ)
        append!(df,new_df)
        max_model[i] = maximum(df.μ)
    end
    return max_model
end

"""
model_random(design,n_first,n_samples,Toy_data,n_cycles)
Simple random model to use as a baseline to compear with other models.
"""
function model_random(design,n_first,n_samples,Toy_data,n_cycles)
    rng = MersenneTwister()
    max_random = Vector{Float64}(undef,n_cycles)
    #first steps
    space = getspace(design)
    new_constructs = StatsBase.sample(rng,space,n_first)
    lab_μ,lab_σ = go2lab(rng,new_constructs,Toy_data)
    df = DataFrame(constructs = new_constructs,μ = lab_μ, σ = lab_σ);
    #get the full space explicitly
    Unseen_space = getspace(design,full = true).space
    for i in 1:n_cycles
        Unseen_space = filter(x-> !(x in new_constructs) , Unseen_space)
        new_constructs = StatsBase.sample(rng,Unseen_space,n_samples)
        lab_μ,lab_σ = go2lab(rng,new_constructs,Toy_data);
        new_df = DataFrame(constructs = new_constructs, μ = lab_μ, σ = lab_σ)
        append!(df,new_df)
        max_random[i] = maximum(df.μ)
    end
    return max_random
end

"""
update_prior_model_linear(design,n_first,n_samples,Toy_data,n_cycles)
Linear model that updates the posterior in every cycle
"""
function update_prior_model_linear(design,n_first,n_samples,Toy_data,n_cycles)
    rng = MersenneTwister()    # to save the output
    max_model = Vector{Float64}(undef,n_cycles)
    #first steps
    space = getspace(design)
    #save usemodel and new constructs for previous round (model = 1, new_df = 2)
    model_vector = Array{BayesianLinearRegressor}(undef,1)
    df_new_vector = Array{DataFrame}(undef,1)
    #start
    new_constructs = StatsBase.sample(rng,space,n_first)
    lab_μ,lab_σ = go2lab(rng,new_constructs,Toy_data)
    df = DataFrame(constructs = new_constructs,μ = lab_μ, σ = lab_σ);
    #get the full space explicitly
    Unseen_space = getspace(design,full = true).space
    for i in 1:n_cycles
        # make model, update prior and only use new data points
        if i == 1
            model_vector[1] = BO_Lin_wrapper(df.constructs,df.μ,df.σ,design.mod)
        else
            model_vector[1] = BO_Lin_wrapper(model_vector[1],df_new_vector[1].constructs,df_new_vector[1].μ,df_new_vector[1].σ,design.mod)
        end
        # filter space
        Unseen_space = filter(x-> !(x in new_constructs) , Unseen_space)
        #sample with TS
        new_constructs = Linear_TS_sampling(model_vector[1],design.mod,Unseen_space,n_samples)
        #evaluated new points
        lab_μ,lab_σ = go2lab(rng,new_constructs,Toy_data);
        df_new_vector[1] = DataFrame(constructs = new_constructs, μ = lab_μ, σ = lab_σ)
        append!(df,df_new_vector[1])
        max_model[i] = maximum(df.μ)
    end
    return max_model
end

"""
make_plot(plt,max_values::Matrix; label="label",color = :blue)

Generates the plots to compear different models.
"""
function make_plot(plt,max_values::Matrix; label="label",color = :blue)
    mean_m = [mean(i) for i in eachcol(max_values)]
    sigma_m = [std(i) for i in eachcol(max_values)]
    return plot!(plt,collect(1:size(max_values)[2]),mean_m,ribbon = 2*sigma_m./sqrt(size(max_values)[1]), label=label, color = color ,lw=1.5,fillalpha = 0.3)
end


text:
# Generate documentation with this command:
# (cd docs && julia --color=yes make.jl)

push!(LOAD_PATH, "..")

using Documenter
using DualNumbers

makedocs(; sitename="DualNumbers", format=Documenter.HTML(), modules=[DualNumbers])

deploydocs(; repo="github.com/eschnett/DualNumbers.jl.git", devbranch="main", push_preview=true)


text:
# diffusion model samples

struct DMBoundOutcome  <: ValueSupport end
struct DMBoundsOutcome <: ValueSupport end

const DMBoundSampleable = Sampleable{Univariate, DMBoundOutcome}
const DMBoundsSampleable = Sampleable{Multivariate, DMBoundsOutcome}

Base.eltype(::Type{DMBoundOutcome}) = Float64
Base.eltype(::Type{DMBoundsOutcome}) = Float64, Bool

Base.length(::DMBoundSampleable) = 1
Base.length(::DMBoundsSampleable) = 2

# general Euler-Maruyama sampler

sampler(d::AbstractDrift, b::AbstractBounds) = DMBoundsSampler(d, b)

struct DMBoundsSampler
    d::AbstractDrift
    b::AbstractBounds
    sqrtdt::Float64

    function DMBoundsSampler(d::AbstractDrift, b::AbstractBounds)
        dt = getdt(d)
        @assert dt == getdt(b)
        return new(d, b, √(dt))
    end
end

function rand(s::DMBoundsSampler)
    x, n, dt, maxn = 0.0, 1, getdt(s.d), min(getmaxn(s.d), getmaxn(s.b))
    while n < maxn
        x += getmu(s.d, n) * dt + s.sqrtdt * randn()
        n += 1
        if x >= getubound(s.b, n)
            return (n - 1) * dt, true
        elseif x <= getlbound(s.b, n)
            return (n - 1) * dt, false
        end
    end
    # no bound crossing until n - return random sample and Inf
    return Inf, rand() < 0.5
end

# faster sampler for constant bounds

# samples only fpt, assumes bounds at -1, 1, unit variance
abstract type DMConstSymBoundsFPTSampler end

struct DMConstSymBoundsNormExpSampler <: DMConstSymBoundsFPTSampler
    mu2::Float64
    a::Float64
    sqrtamu::Float64
    fourmu2π::Float64
    Cf1s::Float64
    CF1st::Float64
    Cf1l::Float64
    F1lt::Float64
    F1inf::Float64

    function DMConstSymBoundsNormExpSampler(mu::Real)
        mu2 = abs2(mu)
        tthresh = 0.12+0.5exp(-abs(mu)/3)
        a = (3 + √(9 + 4mu2)) / 6
        sqrtamu = √((a-1) * mu2 / a)
        fourmu2π = (4mu2 + abs2(π)) / 8
        Cf1s = √(a) * exp(-sqrtamu)
        Cf1l = π / 4fourmu2π
        CF1st = Cf1s * erfc(1 / √(2a * tthresh))
        F1lt = - expm1(-tthresh * fourmu2π)
        F1inf = CF1st + Cf1l * (1 - F1lt)
        return new(mu2, a, sqrtamu, fourmu2π, Cf1s, CF1st, Cf1l, F1lt, F1inf)
    end
end

function rand(s::DMConstSymBoundsNormExpSampler)
    while true
        P = s.F1inf * rand()
        if P <= s.CF1st
            # short-time series
            t = 1 / (2s.a * abs2(erfcinv(P / s.Cf1s)))
            !_acceptt(t, exp(- 1/(2s.a*t) - s.sqrtamu + s.mu2*t/2), 1 / 2t) || return t
        else
            # long-time series
            t = -log1p(- (P - s.CF1st)/s.Cf1l - s.F1lt) / s.fourmu2π
            π2t8 = abs2(π) * t / 8
            !_acceptt(t, exp(-π2t8), π2t8) || return t
        end
    end
end

struct DMConstSymBoundsInvNormSampler <: DMConstSymBoundsFPTSampler
    invabsmu::Float64
    invmu2::Float64

    function DMConstSymBoundsInvNormSampler(mu::Real)
        @assert mu > 0.0
        return new(1 / abs(mu), 1 / abs2(mu))
    end
end

function rand(s::DMConstSymBoundsInvNormSampler)
    while true
        t = _invgaussrand(s.invabsmu, s.invmu2)
        one2t = 1 / 2t
        if t < 2.5
            # short-time series, always accept for mu > 1500
            !_acceptt(t, exp(-one2t), one2t) && !(s.invabsmu < 0.000666) || return t
        else
            # long-time series
            Cl = -log(π/4) - 0.5log(twoπ)
            !_acceptt(t, exp(Cl - one2t - 3log(t)/2), abs2(π) * t / 8) || return t
        end
    end
end

function fastfptsampler(mu)
    absmu = abs(mu)
    return absmu < 1.0 ? DMConstSymBoundsNormExpSampler(absmu) :
                         DMConstSymBoundsInvNormSampler(absmu)
end

# sampler for inverse-Gamma distribution with lambda = 1, mean = mu, mu > 0
function _invgaussrand(mu, mu2)
    y = abs2(randn())
    x = mu + 0.5mu2 * y - 0.5mu * sqrt(4mu * y + mu2 * abs2(y))
    return rand() <= 1 / (1 + x / mu) ? x : mu2 / x
end

function _acceptt(t, f, c2)
    @assert c2 > 0.06385320297074884 # log(5/3) / 16, req. for convergence
    z = f * rand()
    b, twok = exp(-c2), 3
    while true
        z < b || return false
        b -= twok * exp(-c2 * abs2(twok))
        z > b || return true
        twok += 2
        b += twok * exp(-c2 * abs2(twok))
        twok += 2
    end
end

sampler(d::ConstDrift, b::ConstSymBounds) = DMConstSymBoundsSampler(d, b)

struct DMConstSymBoundsSampler
    b2::Float64
    pu::Float64
    fpts::DMConstSymBoundsFPTSampler

    function DMConstSymBoundsSampler(d::ConstDrift, b::ConstSymBounds)
        mu = getmu(d)
        theta = getbound(b)
        mutheta = theta * mu
        return new(abs2(theta), 1 / (1 + exp(-2mutheta)), fastfptsampler(mutheta))
    end
end
rand(s::DMConstSymBoundsSampler) = rand(s.fpts) * s.b2, rand() < s.pu

sampler(d::ConstDrift, b::ConstAsymBounds) = DMConstAsymBoundsSampler(d, b)

struct DMConstAsymBoundsSampler
    mu::Float64
    bup::Float64
    blo::Float64

    DMConstAsymBoundsSampler(d::ConstDrift, b::ConstAsymBounds) =
        new(getmu(d), getubound(b), getlbound(b))
end

function rand(s::DMConstAsymBoundsSampler)
    t, x = 0.0, 0.0
    while true
        xlo, xup = x - s.blo, s.bup - x
        if isapprox(xlo, xup)
            # symmetric bounds, diffusion model in [x - xup, x + xup]
            mutheta = xup * s.mu
            return t + abs2(xup) * rand(fastfptsampler(mutheta)),
                   rand() < 1 / (1 + exp(-2mutheta))
        elseif xlo > xup
            # x closer to upper bound, diffusion model in [x - xup, x + xup]
            mutheta = xup * s.mu
            t += abs2(xup) * rand(fastfptsampler(mutheta))
            rand() >= 1 / (1 + exp(-2mutheta)) || return t, true
            x -= xup
        else
            # x closer to lower bound, diffusion model in [x - xlo, x + xlo]
            mutheta = xlo * s.mu
            t += abs2(xlo) * rand(fastfptsampler(mutheta))
            rand() <= 1 / (1 + exp(-2mutheta)) || return t, false
            x += xlo
        end
    end
end


text:
function scenarios(n::Int; seed=nothing, Ni=1000)
  seed == nothing || Random.seed!(seed)

  if n == 0  # TEST CASE!
    return scenarios(1, seed=seed, Ni=100)
  elseif n == 1
    return CDE.Model.generate_samples(NC=Ni, NT=Ni, gammaC=0.1, gammaT=0.1,
                                      etaC=[.25, .75, 0, 0],
                                      etaT=[.05, .05, .9, 0],
                                      loc = [-1.5, 3.5, 5.1, 5],
                                      scale = [1.6, 1.76, 1.76, 1.6],
                                      df = [12, 10, 10, 15.],
                                      skew = [0, -10, -10, 0.])
  elseif n == 2
    return CDE.Model.generate_samples(NC=Ni, NT=Ni, gammaC=0.1, gammaT=0.1,
                                      etaC=[.05, .05, .5, .4],
                                      etaT=[.05, .05, .9, 0],
                                      loc = [-1.5, 3.5, 5.1, 4.3],
                                      scale = [1.6, 1.76, 1.76, 1.6],
                                      df = [12, 10, 10, 15.],
                                      skew = [0, -10, -10, -11.])
  elseif n == 3
    return CDE.Model.generate_samples(NC=Ni, NT=Ni, gammaC=0.1, gammaT=0.05,
                                      etaC=[.05, .05, .5, .4],
                                      etaT=[.05, .05, .5, .4],
                                      loc = [-1.5, 3.5, 5.1, 4.3],
                                      scale = [1.6, 1.76, 1.76, 1.6],
                                      df = [12, 10, 10, 15.],
                                      skew = [0, -10, -10, -11.])
  elseif n == 4
    return CDE.Model.generate_samples(NC=Ni, NT=Ni, gammaC=0.1, gammaT=0.1,
                                      etaC=[.05, .05, .5, .4],
                                      etaT=[.05, .05, .5, .4],
                                      loc = [-1.5, 3.5, 5.1, 4.3],
                                      scale = [1.6, 1.76, 1.76, 1.6],
                                      df = [12, 10, 10, 15.],
                                      skew = [0, -10, -10, -11.])
  elseif n == 5  # a test case
    return CDE.Model.generate_samples(NC=Ni, NT=Ni, gammaC=0.1, gammaT=0.1,
                                      etaC=[0, 1., 0, 0],
                                      etaT=[0, 0, 1., 0],
                                      loc = [-1.5, 3.5, 5.1, 4.3],
                                      scale = [1.6, 1.76, 1.76, 1.6],
                                      df = [12, 10, 10, 15.],
                                      skew = [0, -10, -10, -11.])
  else
    throw(ArgumentError("n=$(n) was provided. But n must be an integer between 0 and 5."))
  end
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
include("Include.jl")

data_dictionary, TXTL_dictionary, species_list, Mean, Upper, Lower = LoadDictionaries()

t0_exp = 0
tf_exp = 3
tstep_exp = 0.01
experimental_time_exp = [t0_exp:tstep_exp:tf_exp;]

metabolite_list = data_dictionary["list_of_metabolite_symbols"]
rxn_list = data_dictionary["list_of_reaction_strings"]

(n_species,n_rxn) = size(data_dictionary["stoichiometric_matrix"])
number_of_fluxes = find(rxn_list.=="tRNA_charging_M_val_L_c_CAT::16.0*M_val_L_c+16.0*M_atp_c+16.0*tRNA+16.0*M_h2o_c --> 16.0*M_val_L_c_tRNA+16.0*M_amp_c+16.0*M_ppi_c")[1]

t0 = 0
tf = 3
tstep = 0.01
experimental_time = t0:tstep:tf
length_time = length(experimental_time)
data_dictionary["tstep"] = tstep

# this determines which parameters are fitted to data
dataset_index = [15; 18; 22; 24; 27; 29; 31; 34; 38; 41; 42; 59; 60; 61; 66; 69; 70; 76; 79; 83; 84; 86; 88; 89; 106; 110; 114; 121; 123; 126; 128; 130; 132; 133; 134; 135; 140]
non_dataset_index = collect(1:n_species)
deleteat!(non_dataset_index,dataset_index)

species_list = data_dictionary["list_of_metabolite_symbols"]
species_list_dataset = species_list[dataset_index]
species_list_non_dataset = species_list[non_dataset_index]

species_constraint_index = copy(dataset_index)

plot_color = "black"

num_constraint_sets = 1#size(species_constraint_array,1)

# not currently used
parameter = ones(202,2)
parameter[:,1] = parameter[:,1]*110
parameter[:,2] = parameter[:,2]*0.03

initial_conditions = data_dictionary["initial_conditions"]

# set up initial conditions (use kinetic model data)
time_state_array = zeros(n_species,length_time)
time_state_array[:,1] = initial_conditions
time_flux_array = zeros(n_rxn,length_time)
FBA = data_dictionary["default_flux_bounds_array"]
SBA = data_dictionary["species_bounds_array"]

EXIT_FLAG = zeros(length_time-1,num_constraint_sets)

overall_state_array = zeros()
# loop for all constraint subsets
FVA_run = true

constraint_index_array = collect(1:num_constraint_sets)

Exit_flag = Int64[]
uptake_array = 0
flux_array = 0

if !isdir("Ens")
	mkdir("Ens")
end
if !isdir("Ens/Best")
	mkdir("Ens/Best")
end

syn_data = Dict{AbstractString,Any}()

fva_time_array = collect(0.0:0.25:3.0)

fluxes_to_optimize = [[2 3], [23 24], 26, 37, [73 74], 21, 75]

species_to_vary = readdlm("species_to_vary")
fixed_indices = Int64[]
for i in 1:length(species_list)
	if !in(species_list[i],species_to_vary)
		push!(fixed_indices,i)
	end
end

species_boolean = zeros(size(species_list))
for i in 1:length(species_list)
	species_boolean[i] = in(i,dataset_index)
end

switch_prob = .25

if ~isfile("Ens/num_dir")
    writedlm("Ens/num_dir",0)
end

alpha = 1
num_to_reset = 50

if isfile("Ens/Best/error")
	error_best = readdlm("Ens/Best/error")[1]
	species_boolean = readdlm("Ens/Best/species_boolean")
else
	error_best = readdlm("error_best")[1]
end
species_boolean_best = copy(species_boolean)
error = copy(error_best)

error_bc = 3608115.9638814344
flux_uncertainty_bc = 4000 # based on branch points only 20932.739199515625

num_dir = convert(Int64,readdlm("Ens/num_dir")[1]) # Number of pre-existing directories

num_iter = 10000
for i = num_dir+1:num_dir+num_iter
	# Reset periodically
	if rem(i,num_to_reset) == 0
		species_boolean_new = copy(species_boolean_best)
	else
		species_boolean_new = copy(species_boolean)
	end

	# Perturb constraint array
	for j in 1:length(species_boolean_new)
#		if !in(j,fixed_indices)
			if rand() < switch_prob
				species_boolean_new[j] = 1-species_boolean_new[j]
			end
#		end
	end

	species_boolean_new[fixed_indices] = 0

	num_diff = Int64(sum(abs(species_boolean_new-species_boolean)))

	# discretized dfba
	state_array = deepcopy(initial_conditions)
	index = 1

	species_constraint_index_new = find(species_boolean_new)

	syn_data_upper = zeros(n_species,length_time)
	syn_data_lower = zeros(n_species,length_time)
	syn_mean = zeros(n_species,length_time)
	for syn_index in species_constraint_index_new
		itp = interpolate((experimental_time_exp,),Mean[metabolite_list[syn_index]],Gridded(Linear()))
		Mean_itp = itp[experimental_time]
		itp = interpolate((experimental_time_exp,),Lower[metabolite_list[syn_index]],Gridded(Linear()))
		Lower_itp = itp[experimental_time]
		itp = interpolate((experimental_time_exp,),Upper[metabolite_list[syn_index]],Gridded(Linear()))
		Upper_itp = itp[experimental_time]
		syn_data_upper[syn_index,:] = Upper_itp
		syn_data_lower[syn_index,:] = Lower_itp
		syn_mean[syn_index,:] = Mean_itp
	end
	syn_data["mean"] = syn_mean
	syn_data["upper"] = syn_data_upper
	syn_data["lower"] = syn_data_lower

	Exit_flag = Int64[]
	uptake_array = 0
	flux_array = 0

	FVA_min = Dict()
	FVA_max = Dict()
	Percentage_failed_tps_min = Dict()
	Percentage_failed_tps_max = Dict()
	for rxn_idx in fluxes_to_optimize
		key = split(rxn_list[rxn_idx[1]],':')[1]
		FVA_min[key] = zeros(length(fva_time_array),191) # Only record 191 fluxes
		FVA_max[key] = zeros(length(fva_time_array),191) # Only record 191 fluxes
	end

	#dfba
	fva_time = experimental_time[1:end-1]
	for j in 1:length(fva_time)
		time_index = fva_time[j]
		# constraint
		data_dictionary["objective_coefficient_array"] = zeros(number_of_fluxes)
		data_dictionary["objective_coefficient_array"][171] = -1
		data_dictionary = Bounds(data_dictionary,TXTL_dictionary,state_array)
		species_constraints = calculate_constraints(state_array, parameter, species_constraint_index_new, syn_data, index, tstep)
		data_dictionary["species_bounds_array"] = species_constraints
		#  data_dictionary["default_flux_bounds_array"] = flux_constraints
		data_dictionary["state_array"] = state_array

		# solve the lp problem -
		(objective_value, flux_array, dual_array, uptake_array, exit_flag) = FluxDriver(data_dictionary)
		push!(Exit_flag,exit_flag)

		# mass balance for concentration
		state_array = state_array + (uptake_array * tstep) # cell free mass balanace
		state_array[101] = initial_conditions[101]
		Z0 = objective_value
		if FVA_run && in(time_index,fva_time_array)
			for rxn_idx in fluxes_to_optimize
				key = split(rxn_list[rxn_idx[1]],':')[1]
				if length(rxn_idx) == 1
					#constraint objective flux to value determined
					data_dictionary["default_flux_bounds_array"][171] = Z0

					# minimize flux
					data_dictionary["objective_coefficient_array"] = zeros(number_of_fluxes)
					data_dictionary["objective_coefficient_array"][rxn_idx] = 1
					(objective_value, flux_array, dual_array, uptake_array, exit_flag) = FluxDriver(data_dictionary)
					FVA_min[key][Int64(time_index*4+1),:] = flux_array[1:191]

					# maximize flux
					data_dictionary["objective_coefficient_array"] = zeros(number_of_fluxes)
					data_dictionary["objective_coefficient_array"][rxn_idx] = -1
					(objective_value, flux_array, dual_array, uptake_array, exit_flag) = FluxDriver(data_dictionary)
					FVA_max[key][Int64(time_index*4+1),:] = flux_array[1:191]
				else
					#constraint objective flux to value determined
					data_dictionary["default_flux_bounds_array"][171] = Z0

					# maximize reverse flux
					data_dictionary["objective_coefficient_array"] = zeros(number_of_fluxes)
					data_dictionary["objective_coefficient_array"][rxn_idx[2]] = -1
					(objective_value, flux_array, dual_array, uptake_array, exit_flag) = FluxDriver(data_dictionary)
					FVA_min[key][Int64(time_index*4+1),:] = flux_array[1:191]

					# maximize forward flux
					data_dictionary["objective_coefficient_array"] = zeros(number_of_fluxes)
					data_dictionary["objective_coefficient_array"][rxn_idx[1]] = -1
					(objective_value, flux_array, dual_array, uptake_array, exit_flag) = FluxDriver(data_dictionary)
					FVA_max[key][Int64(time_index*4+1),:] = flux_array[1:191]
				end
			end
		end # if FVA_run

		index = index + 1
		time_state_array[:,index] = state_array
		time_flux_array[:,index] = flux_array

	end # for time_index in [t0:tstep:tf-tstep;]

	if sum(5-Exit_flag) == 0
		error_vector = CalcError(Upper,Lower,Mean,experimental_time,time_state_array,data_dictionary)

		FVA_MIN = zeros(length(fva_time_array),length(fluxes_to_optimize))
		for j in 1:length(fluxes_to_optimize)
			rxn_idx = fluxes_to_optimize[j]
			key = split(rxn_list[rxn_idx[1]],':')[1]
			if length(rxn_idx) == 1
				FVA_MIN[:,j] = FVA_min[key][:,rxn_idx]
			else
				FVA_MIN[:,j] = FVA_min[key][:,rxn_idx[1]]-FVA_min[key][:,rxn_idx[2]]
			end
		end
		FVA_MAX = zeros(length(fva_time_array),length(fluxes_to_optimize))
		for j in 1:length(fluxes_to_optimize)
			rxn_idx = fluxes_to_optimize[j]
			key = split(rxn_list[rxn_idx[1]],':')[1]
			if length(rxn_idx) == 1
				FVA_MAX[:,j] = FVA_max[key][:,rxn_idx]
			else
				FVA_MAX[:,j] = FVA_max[key][:,rxn_idx[1]]-FVA_max[key][:,rxn_idx[2]]
			end
		end
		flux_uncertainty = norm(FVA_MAX-FVA_MIN)

		error_new = sum(error_vector)*flux_uncertainty
#		error_new = length(species_constraint_index_new)

		if sum(error_vector) > error_bc
			error_new *= 1e9
		end
		if flux_uncertainty > flux_uncertainty_bc
			error_new *= 1e9
		end
		if sum(species_boolean_new) == 0
			error_new *= 1e9
		end

		acc_prob = exp(alpha*(error-error_new)/error)

		error_round = round(error/10^floor(log(10,error)),3)*10^floor(log(10,error))
		error_new_round = round(error_new/10^floor(log(10,error_new)),3)*10^floor(log(10,error_new))
		error_best_round = round(error_best/10^floor(log(10,error_best)),3)*10^floor(log(10,error_best))
		acc_prob_round = round(acc_prob,2)

		# If a new best is achieved, overwrite parameters and cost and save the relevant information to Best directory
		if error_new < error_best
		    # Save to Best directory
		    writedlm("Ens/Best/species_boolean",species_boolean_new)
		    writedlm("Ens/Best/time_state_array",time_state_array)
		    writedlm("Ens/Best/exit_flag",Exit_flag)
		    writedlm("Ens/Best/state_error",sum(error_vector))
		    writedlm("Ens/Best/flux_uncertainty",flux_uncertainty)
		    writedlm("Ens/Best/error",error_new)
		    species_boolean_best = copy(species_boolean_new)
		    error_best = copy(error_new)
		    println("$i: error_new = $error_new_round, error = $error_round, best = $error_best_round, acc_prob = $acc_prob_round, num_diff = $num_diff, NEW BEST")
		else
		    println("$i: error_new = $error_new_round, error = $error_round, best = $error_best_round, acc_prob = $acc_prob_round, num_diff = $num_diff")
		end
		# If new cost is lower than previous cost, choose new params as reference point for perturbation
		if rand() < acc_prob # Accept all better sets and some worse sets
		    # Create directory for next sample and save the relevant information
		    num_dir += 1
			writedlm("Ens/num_dir",i) # Record new number of directories
		    mkdir("Ens/$num_dir")
		    writedlm("Ens/$num_dir/species_boolean",species_boolean_new)
		    writedlm("Ens/$num_dir/time_state_array",time_state_array)
		    writedlm("Ens/$num_dir/exit_flag",Exit_flag)
		    writedlm("Ens/$num_dir/state_error",sum(error_vector))
		    writedlm("Ens/$num_dir/flux_uncertainty",flux_uncertainty)
		    writedlm("Ens/$num_dir/error",error_new)
		    species_boolean = copy(species_boolean_new)
		    error = copy(error_new)
			for rxn_idx in fluxes_to_optimize
				key = split(rxn_list[rxn_idx[1]],':')[1]
				mkdir("Ens/$num_dir/$key")
				writedlm("Ens/$num_dir/$key/FVA_min",FVA_min[key])
				writedlm("Ens/$num_dir/$key/FVA_max",FVA_max[key])
			end
		end
	end
end


text:
struct MarginalDist{N,D<:Distribution,VS<:AbstractValueShape}
    dims::NTuple{N,Int}
    dist::D
    origvalshape::VS
end

function _get_edges(data::Tuple, nbins::Tuple{Vararg{<:Integer}}, closed::Symbol)
    return StatsBase.histrange(data, StatsBase._nbins_tuple(data, nbins), closed)
end

function _get_edges(data::Any, nbins::Integer, closed::Symbol)
    return StatsBase.histrange((data, ), StatsBase._nbins_tuple((data, ), (nbins,)), closed)[1]
end

function _get_edges(data::Any, nbins::Union{AbstractRange, Tuple{AbstractRange}}, closed::Symbol)
    return nbins
end

function get_marginal_dist(
    maybe_shaped_samples::DensitySampleVector,
    key::Union{Integer, Symbol, Expr};
    bins = 200,
    closed::Symbol = :left,
    filter::Bool = false
)
    samples = BAT.unshaped.(maybe_shaped_samples)

    if filter
        samples = BAT.drop_low_weight_samples(samples)
    end

    idx = asindex(maybe_shaped_samples, key)
    s = flatview(samples.v)[idx, :]

    edges = _get_edges(s, bins, closed)

    hist = fit(Histogram,
            s,
            FrequencyWeights(samples.weight),
            edges,
            closed = closed)

    uvbd = EmpiricalDistributions.UvBinnedDist(hist)
    marg = MarginalDist((idx,), uvbd, varshape(maybe_shaped_samples))

    return (result = marg, )
end

function get_marginal_dist(
    maybe_shaped_samples::DensitySampleVector,
    key::Union{NTuple{n,Integer}, NTuple{n,Union{Symbol, Expr}}} where n;
    bins = 200,
    closed::Symbol = :left,
    filter::Bool = false
)
    samples = unshaped.(maybe_shaped_samples)

    if filter
        samples = BAT.drop_low_weight_samples(samples)
    end

    idxs = asindex.(Ref(maybe_shaped_samples), key)
    s = Tuple(BAT.flatview(samples.v)[i, :] for i in idxs)

    edges = if isa(bins, Integer)
        _get_edges(s, (bins,), closed)
    else
        Tuple(_get_edges(s[i], bins[i], closed) for i in 1:length(bins))
    end

    hist = fit(Histogram,
            s,
            FrequencyWeights(samples.weight),
            edges,
            closed = closed)

    mvbd = EmpiricalDistributions.MvBinnedDist(hist)
    marg =  MarginalDist(idxs, mvbd, varshape(maybe_shaped_samples))

    return (result = marg, )
end

#for prior
function get_marginal_dist(
    prior::NamedTupleDist,
    key::Union{Integer, Symbol};
    bins = 200,
    edges = nothing,
    closed::Symbol = :left,
    nsamples::Integer = 10^6
)
    idx = asindex(prior, key)
    r = flatview(unshaped.(rand(prior, nsamples)))

    edges = _get_edges(r[idx, :], bins, closed)

    hist = fit(Histogram, r[idx, :], edges, closed = closed)

    uvbd = EmpiricalDistributions.UvBinnedDist(hist)
    marg = MarginalDist((idx,), uvbd, varshape(prior))

    return (result = marg, )
end

function get_marginal_dist(
    prior::NamedTupleDist,
    key::Union{NTuple{2, Symbol}, NTuple{2, Integer}};
    bins = 200,
    closed::Symbol = :left,
    nsamples::Integer = 10^6
)
    idxs = asindex.(Ref(prior), key)

    r = flatview(unshaped.(rand(prior, nsamples)))
    s = Tuple(r[i, :] for i in idxs)

    edges = if isa(bins, Integer)
        _get_edges(s, (bins,), closed)
    else
        Tuple(_get_edges(s[i], bins[i], closed) for i in 1:length(bins))
    end

    hist = fit(Histogram,
            s,
            edges,
            closed = closed)

    mvbd = EmpiricalDistributions.MvBinnedDist(hist)
    marg =  MarginalDist(idxs, mvbd, varshape(prior))

    return (result = marg, )
end

function get_marginal_dist(
    original::MarginalDist,
    vsel::NTuple{n, Int} where n
)
    original_hist = convert(Histogram, original.dist)
    dims = collect(1:ndims(original_hist.weights))
    vsel = Tuple(findfirst(x-> x == p, original.dims) for p in vsel)

    weights = sum(original_hist.weights, dims=setdiff(dims, vsel))
    weights = dropdims(weights, dims=Tuple(setdiff(dims, vsel)))

    edges = Tuple([original_hist.edges[p] for p in vsel])
    hist = StatsBase.Histogram(edges, weights, original_hist.closed)

    bd = if length(vsel) == 1
        EmpiricalDistributions.UvBinnedDist(hist)
    else
        EmpiricalDistributions.MvBinnedDist(hist)
    end

    marg = MarginalDist(vsel, bd, original.origvalshape)

    return (result = marg, )
end


text:
"""
    DAMMmat(Tₛ::Array{Float64, 1}, θ::Array{Float64, 1}, Rₛ::Array{Float64, 1}, r::Int64)

Generates a matrix of DAMM output for gridded inputs x and y
Inputs: soil temperature (Tₛ), soil moisture (θ), respiration (Rₛ), resolution (r)

# Examples:
```julia-repl
julia> Tₛ = collect(15.0:2.5:40.0)
julia> θ = collect(0.2:0.05:0.7)
julia> Rₛ = [1.0, 1.2, 1.5, 2.0, 2.7, 3.8, 4.9, 6.7, 4.1, 2.0, 0.4]
julia> r = 10
julia> out = DAMMmat(Tₛ, θ, Rₛ, r)
```

    DAMMmat(Tₛ::Array{Float64, 1}, θ::Array{Float64, 1}, Rₛ::Array{Float64, 1}, r::Int64, n::Int64)

Bin data by n quantiles

# Examples:
```julia-repl
julia> n = 4
julia> out = DAMMmat(Tₛ, θ, Rₛ, r, n)
```
"""
function DAMMmat(Tₛ::Array{Float64, 1}, θ::Array{Float64, 1}, Rₛ::Array{Float64, 1}, r::Int64)
  poro_val = maximum(θ) + 0.2 # will sometime fail to fit if equal maximum(θ)
  params = DAMMfit(hcat(Tₛ, θ), Rₛ, poro_val)
# create x y arrays, and z matrix, to plot DAMM surface
  x = collect(range(minimum(Tₛ), length=r, stop=maximum(Tₛ))) # T axis, °C from min to max
  y = collect(range(minimum(θ), length=r, stop=maximum(θ))) # M axis, % from min to max
  X = repeat(1:r, inner=r) # X for DAMM matrix
  Y = repeat(1:r, outer=r) # Y for DAMM matrix
  X2 = repeat(x, inner=r) # T values to fit DAMM on
  Y2 = repeat(y, outer=r) # M values to fit DAMM on
  xy = hcat(X2, Y2) # T and M matrix to create DAMM matrix
  DAMM_Matrix = Matrix(sparse(X, Y, DAMM(xy, params)))
  return sDAMMmat(poro_val, params, x, y, DAMM_Matrix)
end

function DAMMmat(Tₛ::Array{Float64, 1}, θ::Array{Float64, 1}, Rₛ::Array{Float64, 1}, r::Int64, n::Int64)
  Tmed, θmed, Rmed = qbins(Tₛ, θ, Rₛ, n)
  poro_val = maximum(θ) + 0.2 # will sometime fail to fit if equal maximum(θ)
  params = DAMMfit(hcat(Tmed, θmed), Rmed, poro_val)
# create x y arrays, and z matrix, to plot DAMM surface
  x = collect(range(minimum(Tmed), length=r, stop=maximum(Tmed))) # T axis, °C from min to max
  y = collect(range(minimum(θmed), length=r, stop=maximum(θmed))) # M axis, % from min to max
  X = repeat(1:r, inner=r) # X for DAMM matrix
  Y = repeat(1:r, outer=r) # Y for DAMM matrix
  X2 = repeat(x, inner=r) # T values to fit DAMM on
  Y2 = repeat(y, outer=r) # M values to fit DAMM on
  xy = hcat(X2, Y2) # T and M matrix to create DAMM matrix
  DAMM_Matrix = Matrix(sparse(X, Y, DAMM(xy, params)))
  return sDAMMmatq(poro_val, Tmed, θmed, Rmed, params, x, y, DAMM_Matrix)
end


text:
using TensorOperations
using LinearAlgebra

function RMP2(Alg::RMP2a)
    aoints = IntegralHelper{Float64}()
    rhf = RHF(aoints)
    moints = IntegralHelper(orbitals=rhf.orbitals)
    RMP2(moints, aoints, Alg)
end

function RMP2(aoints::IntegralHelper{Float64,E,AtomicOrbitals}, Alg::RMP2a) where E <: AbstractERI
    rhf = RHF(aoints)
    moints = IntegralHelper(orbitals=rhf.orbitals)
    RMP2(moints, aoints, Alg)
end

function RMP2(O::AbstractRestrictedOrbitals, Alg::RMP2a)
    aoints = IntegralHelper()
    moints = IntegralHelper(orbitals=O)
    RMP2(moints, aoints, Alg)
end

function RMP2(rhf::RHF, Alg::RMP2a)
    aoints = IntegralHelper()
    moints = IntegralHelper(orbitals=rhf.orbitals)
    RMP2(moints, aoints, Alg)
end

function RMP2(M::Molecule, Alg::RMP2a)
    aoints = IntegralHelper{Float64}(molecule=M, eri_type=SparseERI)
    rhf = RHF(M,aoints)
    moints = IntegralHelper(molecule = M, orbitals=rhf.orbitals)
    RMP2(moints, aoints, Alg)
end

function RMP2(moints::IntegralHelper{T1,Chonky,O}, aoints::IntegralHelper{T2,E,AtomicOrbitals}, Alg::RMP2a) where {T1<:AbstractFloat,
                                                                                T2<:AbstractFloat,O<:RHFOrbitals,E<:AbstractERI}
    mp_header()
    mo_from_ao!(moints, aoints, "OVOV")
    RMP2(moints, Alg)
end

function RMP2(moints::IntegralHelper{T1,E1,O}, aoints::IntegralHelper{T2,E2,AtomicOrbitals}, Alg::RMP2a) where {T1<:AbstractFloat,T2<:AbstractFloat,
                                                                                E1<:AbstractDFERI,E2<:AbstractDFERI,O<:RHFOrbitals}
    mp_header()
    mo_from_ao!(moints, aoints, "BOV")
    RMP2(moints, Alg)
end

function RMP2(moints::IntegralHelper{T1,Chonky,O}, aoints::IntegralHelper{T2,E,AtomicOrbitals}, Alg::RMP2a) where {T1<:AbstractFloat,
                                                                                T2<:AbstractFloat,O<:AbstractRestrictedOrbitals,E<:AbstractERI}
    mp_header()
    mo_from_ao!(moints, aoints, "Fia","OVOV")
    RMP2(moints, Alg)
end

function RMP2(moints::IntegralHelper{T1,E1,O}, aoints::IntegralHelper{T2,E2,AtomicOrbitals}, Alg::RMP2a) where {T1<:AbstractFloat,T2<:AbstractFloat,
                                                                                E1<:AbstractDFERI,E2<:AbstractDFERI,O<:AbstractRestrictedOrbitals}
    mp_header()
    mo_from_ao!(moints, aoints, "Fia","BOV")
    RMP2(moints, Alg)
end

function RMP2(ints::IntegralHelper{<:AbstractFloat,<:AbstractERI,<:AbstractRestrictedOrbitals}, Alg::RMP2a)

    # Check frozen core and inactive virtual
    ndocc = ints.molecule.Nα
    nvir = size(ints.orbitals.C,1) - ndocc
    core = Options.get("drop_occ")
    inac = Options.get("drop_vir")

    if core ≥ ndocc
        throw(FermiException("invalid number of frozen orbitals ($core) for $ndocc doubly occupied orbitals"))
    end
    if inac ≥ nvir
        throw(FermiException("invalid number of inactive virtual orbitals ($inac) for $nvir total virtual orbitals"))
    end

    # Collect RHF information

    output("  Starting MP2 computation")
    output(" Number of frozen orbitals:             {:d}" , core)
    output(" Number of inactive orbitals:           {:d}" , inac)
    output(" Number of correlated electron pairs:   {:d}", ndocc-core)
    output(" Number of correlated virtual orbitals: {:d}", nvir-inac)
    output(" ⇒ Total number of MP2 amplitudes:      {:d}", (ndocc-core)^2*(nvir-inac)^2)

    output(repeat("-",80))

    # Compute MP2 energy
    t = @elapsed Emp2 = RMP2_energy(ints, Alg)
    Eref = ints.orbitals.sd_energy

    output("   @Final RMP2 Correlation Energy {:>20.12f} Eₕ", Emp2)
    output("   Reference Energy               {:>20.12f} Eₕ", Eref)
    output("   @Final RMP2 Total Energy       {:>20.12f} Eₕ", Emp2+Eref)
    output(repeat("-",80))

    RMP2(Emp2, Emp2+Eref)
end

function RMP2_energy(ints::IntegralHelper{T,E,O}, Alg::RMP2a) where {T<:AbstractFloat, E<:AbstractERI,O<:AbstractRestrictedOrbitals}
    inac = Options.get("drop_vir")
    core = Options.get("drop_occ")
    ndocc = ints.molecule.Nα
    nbf = size(ints.orbitals.C,1)
    ϵo = ints["Fd"][(core+1):ndocc]
    ϵv = ints["Fd"][(ndocc+1):(nbf-inac)]

    Emp2 = RMP2_nonrhf_energy(ints, ϵo, ϵv)
    Emp2 += RMP2_rhf_energy(ints, ϵo, ϵv)

    return Emp2
end

function RMP2_energy(ints::IntegralHelper{T,E,RHFOrbitals}, Alg::RMP2a) where {T<:AbstractFloat, E<:AbstractERI}
    inac = Options.get("drop_vir")
    core = Options.get("drop_occ")
    ndocc = ints.molecule.Nα
    nbf = size(ints.orbitals.C,1)
    ϵo = ints["Fd"][(core+1):ndocc]
    ϵv = ints["Fd"][(ndocc+1):(nbf-inac)]

    Emp2 = RMP2_rhf_energy(ints, ϵo, ϵv)

    return Emp2
end

function RMP2_nonrhf_energy(ints::IntegralHelper{T,E,O}, ϵo::AbstractArray{T,1}, ϵv::AbstractArray{T,1}) where {T<:AbstractFloat,E<:AbstractERI, O<:AbstractRestrictedOrbitals}

    t = @elapsed begin
    output("Computing non-RHF contribution to the MP2 energy... ", ending="")

    Dia = [ϵo[i]-ϵv[a] for i = eachindex(ϵo), a = eachindex(ϵv)]
    Fia = ints["Fia"]
    Emp2 = sum(transpose(Fia)*(Fia ./ Dia))
    end
    output("Done in {:3.5f} seconds.", t)

    return Emp2
end

function RMP2_rhf_energy(ints::IntegralHelper{T,RIFIT,O}, ϵo::AbstractArray{T,1}, ϵv::AbstractArray{T,1}) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}
    Bov = ints["BOV"].data

    output(" Computing DF-MP2 Energy... ", ending="")
    v_size = length(ϵv)
    o_size = length(ϵo)

    # Pre-allocating arrays for threads
    BABs = [zeros(T, v_size, v_size) for i=1:Threads.nthreads()]
    BBAs = [zeros(T, v_size, v_size) for i=1:Threads.nthreads()]

    # Vector containing the energy contribution computed by each thread
    ΔMP2s = zeros(T,Threads.nthreads())
    TWO = T(2.0)
    ONE = one(T)
    t = @elapsed begin
        @sync for i in 1:o_size
            Threads.@spawn begin
            id = Threads.threadid()
            @views Bi = Bov[:,i,:]
            BAB = BABs[id]
            BBA = BBAs[id]
            for j in i:o_size

                @views Bj = Bov[:,j,:]
                @tensor BAB[a,b] = Bi[Q,a]*Bj[Q,b]
                transpose!(BBA,BAB)

                eij = ϵo[i] + ϵo[j]
                dmp2_1 = zero(T)
                dmp2_2 = zero(T)
                @fastmath for b in 1:v_size
                    @inbounds eij_b = eij - ϵv[b]
                    for a in 1:v_size
                        @inbounds begin
                            iajb = BAB[a,b]
                            d = iajb/(eij_b - ϵv[a])
                            dmp2_1 += d*iajb
                            dmp2_2 += d*BBA[a,b]
                        end
                    end
                end
                fac = i != j ? TWO : ONE
                dmp2 = TWO*dmp2_1 - dmp2_2
                ΔMP2s[id] += fac*dmp2
            end
            end
        end
        Emp2 = sum(ΔMP2s)
    end
    output("Done in {:5.5f} seconds.", t)
    return Emp2
end

function RMP2_rhf_energy(ints::IntegralHelper{T,Chonky,O}, ϵo::AbstractArray{T,1}, ϵv::AbstractArray{T,1}) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}
    output(" Computing MP2 Energy... ", ending="")
    ovov = ints["OVOV"]

    ΔMP2s = zeros(T,Threads.nthreads())
    TWO = T(2.0)
    t = @elapsed begin
        @sync for b in eachindex(ϵv)
            Threads.@spawn begin
            id = Threads.threadid()
            for a in eachindex(ϵv)
                @inbounds ϵ_ab = ϵv[a] + ϵv[b]
                for j in eachindex(ϵo)
                    @inbounds ϵ_abj = ϵo[j] - ϵ_ab
                    for i in eachindex(ϵo)
                        @fastmath @inbounds ΔMP2s[id] += ovov[i,a,j,b]*(TWO*ovov[i,a,j,b] - ovov[i,b,j,a]) / (ϵo[i]+ϵ_abj)
                    end
                end
            end
            end
        end
    end
    output("Done in {:5.5f} s\n", t)
    return sum(ΔMP2s)
end


text:
const NAME_FLAG_WIDTH = 12
const ARG_HELP_INDENT = NAME_FLAG_WIDTH + 2

function show_help(io::IO, parser::ArgumentParser; exit_when_done = !isinteractive())
    show_usage(io, parser)
    println(io)

    if parser.description !== nothing
        println_wrapped(io, parser.description)
        println(io)
    end

    if !isempty(parser.positional_args)
        print(io, "Positional arguments:\n")

        for arg in parser.positional_args
            help = something(arg.help, "$(arg.name) help")
            print_arg_help(io, arg.name, help)
        end

        println(io)
    end

    if !isempty(parser.optional_args)
        print(io, "Optional arguments:\n")

        for arg in parser.optional_args
            flags = format_flags(arg)
            help = something(arg.help, "$(arg.name) help")
            print_arg_help(io, flags, help)
        end

        println(io)
    end

    if parser.epilog !== nothing
        println_wrapped(io, parser.epilog)
    end

    exit_when_done && exit(0)
    nothing
end

function show_usage(io::IO, parser::ArgumentParser)
    if parser.usage !== nothing
        println_wrapped(io, "Usage: $(parser.usage)\n"; subsequent_indent = ARG_HELP_INDENT, break_long_words = false)
    else
        cmdline_name = length(PROGRAM_FILE) > 0 ? PROGRAM_FILE : "PROGRAM"
        prog = something(parser.prog, cmdline_name)
        options_str = join((format_usage_flag(arg) for arg in parser.optional_args), ' ')
        params_str = join((format_arg_name(arg, false) for arg in parser.positional_args), ' ')

        println_wrapped(io, "Usage: $prog $options_str $params_str"; subsequent_indent = ARG_HELP_INDENT, break_long_words = false)
    end
end

show_usage(parser::ArgumentParser) = show_usage(stdout::IO, parser)

function print_arg_help(io, name_or_flags::AbstractString, help::AbstractString)
    if length(name_or_flags) < NAME_FLAG_WIDTH
        pad = ARG_HELP_INDENT - length(name_or_flags) - 1  # the -1 at the end is for the initial space
        print(io, ' ', name_or_flags, ' '^pad)

        wrapped_help = wrap(help;
            initial_indent = ARG_HELP_INDENT,
            subsequent_indent = ARG_HELP_INDENT,) |> lstrip
        println(io, help)
    else
        println(io, ' ', name_or_flags)
        println_wrapped(io,
            help;
            initial_indent = ARG_HELP_INDENT,
            subsequent_indent = ARG_HELP_INDENT,)
    end
end

show_help(parser::ArgumentParser; exit_when_done::Bool = !isinteractive()) =
    show_help(stdout::IO, parser, exit_when_done = exit_when_done)

function format_flags(arg::Argument)
    arg.nargs === 0 && return join(arg.flags, ", ")
    return join([format_flag(flag, arg.nargs, arg_name(arg)) for flag in arg.flags], ", ")
end

format_usage_flag(arg::Argument) = format_usage_flag(arg.default_flag, arg.nargs, arg_name(arg), arg.required)

function format_usage_flag(flag, nargs, arg_name, required)
    required && return format_flag(flag, nargs, arg_name)
    return '[' * format_flag(flag, nargs, arg_name) * ']'
end

function format_flag(flag, nargs, arg_name)
    nargs === 0 && return flag
    flag_name_str = format_arg_name(arg_name, nargs)
    return "$flag $flag_name_str"
end

format_arg_name(arg::Argument, to_uppercase::Bool = true) =
    format_arg_name(arg_name(arg, to_uppercase), arg.nargs)

function format_arg_name(name, nargs::Integer)
    return join(fill(name, nargs), " ")
end

function format_arg_name(name, nargs)
    nargs === :? && return "[$name]"
    nargs === :* && return "[$name [$name ...]]"
    nargs === :+ && return "$name [$name ...]"

    throw(ArgumentError("Unexpected value for nargs: $(nargs)"))
end


text:
"""
    Codec

Codecs provide the encoding-decoding scheme in a flexible way.
Codecs are titled according to the general attributes that the codec scheme
would require, and NOT the actual encoding/decoding function itself.

This makes things customizable, and allows plug-and-play behaviour with
different variants of the same codec scheme, or for strikingly similar codec
schemes.

!!! note
    Every Codec object MUST have exactly one set of `encode_f::Function` and
    `decode_f::Function` fields.
"""
module Codec

include("simplecodec.jl")
export SimpleCodec

include("onehotonecold.jl")
export OneHotOneCold, build_onehot_vec

end


text:
println("**********Testing DataSet**********")

facts("Creating a DataSet") do
  ds = DataSet()
  @fact typeof(ds) => DataSet
  @fact isempty(ds) => true
end

facts("Indexing a DataSet") do
  cl1 = Cluster()
  cl2 = Cluster()
  cl3 = Cluster()

  fv1 = FeatureVector(["word1" => 4, "word2" => 3])
  fv2 = FeatureVector(["word1" => 1, "word2" => 6])
  fv3 = FeatureVector(["word1" => 3, "word2" => 2])

  cl1["fv1"] = fv1
  cl1["fv2"] = fv2
  cl1["fv3"] = fv3

  cl2["fv1"] = fv1*2
  cl2["fv2"] = fv2*2
  cl2["fv3"] = fv3*2

  ds = DataSet()
  ds["cl1"] = cl1
  ds["cl2"] = cl2
  ds["cl3"] = cl3

  @fact ds["cl1"] => cl1
  @fact ds["cl2"] => cl2
  @fact ds["cl3"] => cl3
end

facts("Getting Keys and Values") do
  ds = DataSet()

  @fact isempty(ds) => true

  cl1 = Cluster()
  cl2 = Cluster()
  cl3 = Cluster()
  cl4 = Cluster()

  fv1 = FeatureVector(["word1" => 4, "word2" => 3])
  fv2 = FeatureVector(["word1" => 1, "word2" => 6])
  fv3 = FeatureVector(["word1" => 3, "word2" => 2])

  cl1["fv1"] = fv1
  cl1["fv2"] = fv2
  cl1["fv3"] = fv3

  cl2["fv1"] = fv1*2
  cl2["fv2"] = fv2*2
  cl2["fv3"] = fv3*2

  ds["cl1"] = cl1
  ds["cl2"] = cl2
  ds["cl3"] = cl3

  k = collect(keys(ds))
  v = collect(values(ds))

  @fact isempty(ds) => false
  @fact "cl1" in k => true
  @fact "cl2" in k => true
  @fact "cl3" in k => true
  @fact "cl4" in k => false

  @fact cl1 in v => true
  @fact cl2 in v => true
  @fact cl3 in v => true
  @fact cl4 in v => false
end


text:
module _AttnEn
export AttnEn
using StatsBase: Histogram, fit
    """
        Av4, (Hxx,Hnn,Hxn,Hnx) = AttnEn(Sig)

    Returns the attention entropy (`Av4`) calculated as the average of the
    sub-entropies (`Hxx`,`Hxn`,`Hnn`,`Hnx`) estimated from the data sequence
    (`Sig`) using a base-2 logarithm.

        Av4, (Hxx, Hnn, Hxn, Hnx) = AttnEn(Sig::AbstractArray{T,1} where T<:Real; Logx::Real=2)

    Returns the attention entropy (`Av4`) and the sub-entropies (`Hxx`,`Hnn`,`Hxn`,`Hnx`)
    from the data sequence (`Sig`) where,
    Hxx:    entropy of local-maxima intervals
    Hnn:    entropy of local minima intervals
    Hxn:    entropy of intervals between local maxima and subsequent minima
    Hnx:    entropy of intervals between local minima and subsequent maxima

    # Arguments:
    `Logx`  - Logarithm base, a positive scalar
              (Enter 0 for natural logarithm)

    See also `EnofEn`, `SpecEn`, `XSpecEn`, `PermEn`, `MSEn`

    # References:
        [1] Jiawei Yang, et al.,
            "Classification of Interbeat Interval Time-series Using
            Attention Entropy."
            IEEE Transactions on Affective Computing
            (2020)

    """
    function AttnEn(Sig::AbstractArray{T,1} where T<:Real; Logx::Real=2)

    (Logx == 0) ? Logx = exp(1) : nothing
    N = size(Sig,1)
    (N > 10) ? nothing : error("Sig:   must be a numeric vector")
    (Logx>0) ? nothing : error("Logx:  must be a positive number > 0")

    Xmax = PkFind(Sig)
    Xmin = PkFind(-Sig)
    Txx = diff(Xmax)
    Tnn = diff(Xmin)
    Temp = diff(sort(vcat(Xmax, Xmin)))

    if isempty(Xmax)
        error("No local maxima found!")
    elseif isempty(Xmin)
        error("No local minima found!")
    end

    (Xmax[1]<Xmin[1]) ? (Txn = Temp[1:2:end]; Tnx = Temp[2:2:end]) :
            (Txn = Temp[2:2:end]; Tnx = Temp[1:2:end])

    Edges = -0.5:N
    Pnx = fit(Histogram,Tnx,Edges).weights
    Pnn = fit(Histogram,Tnn,Edges).weights
    Pxx = fit(Histogram,Txx,Edges).weights
    Pxn = fit(Histogram,Txn,Edges).weights

    Pnx = Pnx[Pnx.!=0]/size(Tnx,1)
    Pxn = Pxn[Pxn.!=0]/size(Txn,1)
    Pnn = Pnn[Pnn.!=0]/size(Tnn,1)
    Pxx = Pxx[Pxx.!=0]/size(Txx,1)

    Hxx = -sum(Pxx.*(log.(Logx,Pxx)))
    Hxn = -sum(Pxn.*(log.(Logx,Pxn)))
    Hnx = -sum(Pnx.*(log.(Logx,Pnx)))
    Hnn = -sum(Pnn.*(log.(Logx,Pnn)))
    Av4 = (Hnn + Hxx + Hxn + Hnx)/4

    return Av4, (Hxx,Hnn,Hxn,Hnx)
    end

    function PkFind(X)
        Nx = size(X,1)
        Indx = zeros(Int,Nx);
        for n = 2:Nx-1
            if X[n-1]< X[n] > X[n+1]
                Indx[n] = n

            elseif X[n-1] < X[n] == X[n+1]
                k = 1
                while (n+k)<Nx && X[n] == X[n+k]
                    k +=1
                end
                if X[n] > X[n+k]
                    Indx[n] = n + floor((k-1)/2)
                end
            end
        end
        Indx = Indx[Indx.!==0]
    return Indx
    end

end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

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

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

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

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
#Demo 4: Empty/undefined plots
#-------------------------------------------------------------------------------
using InspectDR
using Colors

#==Input
===============================================================================#
#Constants
#-------------------------------------------------------------------------------
black = RGB24(0, 0, 0)
white = RGB24(1, 1, 1)
red = RGB24(1, 0, 0)
green = RGB24(0, 1, 0)
blue = RGB24(0, 0, 1)

#Input data
#-------------------------------------------------------------------------------

#==Generate plot
===============================================================================#
mplot = InspectDR.Multiplot(title="Empty/Undefined Plots")
mplot.layout[:ncolumns] = 2
xlabel = "X-Axis (X-Unit)"
ylabel = "Y-Axis (Y-Unit)"
kwargs = (:xlabel=>xlabel, :ylabels=>[ylabel])

#Subplot 1
#-------------------------------------------------------------------------------
plot = add(mplot, InspectDR.Plot2D(:lin, :lin, title="xmax = ∞"; kwargs...))
wfrm = add(plot, Float64[0, Inf], Float64[0, 1])

#Subplot 2
#-------------------------------------------------------------------------------
plot = add(mplot, InspectDR.Plot2D(:lin, :lin, title="xmin = NaN"; kwargs...))
wfrm = add(plot, Float64[NaN, 1], Float64[0, 1])

#Subplot 3
#-------------------------------------------------------------------------------
plot = add(mplot, InspectDR.Plot2D(:lin, :lin, title="ymin = -∞"; kwargs...))
wfrm = add(plot, Float64[0, 1], Float64[-Inf, 0])

#Subplot 4
#-------------------------------------------------------------------------------
plot = add(mplot, InspectDR.Plot2D(:lin, :lin, title="ymax = NaN"; kwargs...))
wfrm = add(plot, Float64[0, 1], Float64[1, NaN])

#Subplot 5
#-------------------------------------------------------------------------------
plot = add(mplot, InspectDR.Plot2D(:lin, :lin, title="No Data"; kwargs...))

#Subplot 6
#-------------------------------------------------------------------------------
plot = add(mplot, InspectDR.Plot2D(:lin, :lin, title="xmin = xmax"; kwargs...))
wfrm = add(plot, Float64[2, 2], Float64[0, 1])

#Display
#-------------------------------------------------------------------------------
gplot = display(InspectDR.GtkDisplay(), mplot)

:DONE


text:
module Gravityfield
# ] activate .
export Unitful
export FreeUnits, 𝐋, ustrip, Quantity
export m ,  m² ,  m³ ,  s ,  s² ,  s³ ,  kg ,  kgm ,  kgm² ,  kgm³
export kgs ,  kgs² ,  kgs³, N

export γ,  m_earth ,  r_earth ,  r_sun ,  d_sun_earth ,  d_moon_earth, ρ_earth
export m_sun ,  m_moon ,  yr ,  ρ_sun
export V, ρ, Pos, Npos
export nest, slice_1, slice_2, slice_3, slice_4, slice_5, slice_6
export Body, Bodies, Sphere
export earthlikeSphere, earthlikeSphere_r_known, earthlikeSphere_m_known
export bodies_3d, bodies_2d, bodies_1d
export rvec, uvec, len, g, gvec, gfieldvals
export xygrids, sqgrids
import Base:show, show_circular, _truncate_at_width_or_chars

using Unitful
#  \bfL<tab> displays as an enclose question mark in Windows Julia,
# and in VScode terminal. On Notepad++, it shows as an Icelandic 'eth'
# Start from cygwin / mintty with a good font, or VSCode, in order
# to get the mathematical bold L for 'length'.
import Unitful: FreeUnits, isunitless, unit
# Time, length, mass
import Unitful: 𝐓, 𝐋, 𝐌
# Draft, going to separate the units to another package:
#@unit transportdistance "kgm" TransportDistance u"kg*m" false
#@unit turn          "τ"             Turn          2π*u"rad"     false
#@unit transportdistance "kgm" TransportDistance 1u"kg*m" false
#@unit transportdistance "kgm" TransportDistance 1u"kg*m" false
#@unit momo "kgm" TransportDistance 1u"kg*m" false
"Shorthands which may cause confusion in other contexts"
const m = 1.0u"m"
const m² = 1.0u"m^2"
const m³ = 1.0u"m^3"
const s = 1.0u"s"
const s² = 1.0u"s^2"
const s³ = 1.0u"s^3"
const kg = 1.0u"kg"
const kgm = kg*m
const kgm² = kg*m²
const kgm³ = kg*m³
const kgs = kg*s
const kgs² = kg*s²
const kgs³ = kg*s³
const N = kgm/s²

"Newton's gravity constant"
const γ=6.67428E-11m³/kgs²
const m_earth = 5.97E24kg
const r_earth = 6300000.0m
const r_sun = 695510000.0m
const d_sun_earth = 1.496e11m
const d_moon_earth = 384400000.0m
const m_sun = 1.989E30kg
const m_moon = 7.34767309e22kg
const yr = 31556926.0s
const ρ_sun = 1410.0kg/m³

"Volume of a sphere"
V(r) = 4π / 3 * r^3
"Density of a sphere"
ρ(;m=1., r=1.) = m / V(r)
const ρ_earth = ρ(m = m_earth, r = r_earth)

#=
We want to print quantities without a space between value and unit.
This enables copying output for easy redefinitions.
For easier reading, we'll also put a bit of colour on units.
julia> 1kg
1.0kg

julia> 1 kg
ERROR: syntax: extra token "kg" after end of expression
=#
function show(io::IO, x::Quantity)
    show(io,x.val)
    if !isunitless(unit(x))
        printstyled(io, unit(x); color=:cyan)
    end
    nothing
end

"Position vector"
const Pos = Vector{Quantity{T,𝐋,U}} where {T,U}
Pos(p::Pos) = Vector(p)

function show(io::IO, p::Pos) # short form
    typ = typeof(p)
    ioc = IOContext(io, :typeinfo => typ)
    show(ioc, ustrip(p))
    printstyled(ioc, unit(eltype(typ)); color=:cyan)
end
function show(io::IO, ::MIME"text/plain", p::Pos{T}) where T# long form
    typ = typeof(p)
    ioc = IOContext(io, :typeinfo => typ)
    print(ioc, "Pos{", T, "}(")
    show(ioc, ustrip(p))
    printstyled(ioc, unit(eltype(typ)); color=:cyan)
    print(ioc, ")")
end

"Nested positions. Each element is a position vector; the container can be a vector or matrix"
const Npos = Union{Array{S, 2}, Array{S,1}} where S<:Pos
Npos(p::Npos) = Array(p)
function show(io::IO, np::Npos) # short form
    np_unitless = map(p->ustrip(p), np)
    typ = typeof(np_unitless)
    ioc = IOContext(io, :typeinfo => typ)
    show(ioc, np_unitless)
    printstyled(ioc, unit(eltype(np[1])); color=:cyan)
    nothing
end
function show(io::IO, ::MIME"text/plain", np::Npos) # long form
    etyp = eltype(np[1])
    np_unitless = map(p->ustrip(p), np)
    typul = eltype(np_unitless[1])
    if ndims(np) == 1
        r= size(np)[1]
        println(io, r, "-element Npos{Pos{", typul, "}}:")
    else
        r, c = size(np)
        println(io, r, "x", c, " Npos{Pos{", typul, "}}:")
    end
    typ = typeof(np_unitless)
    ioc = IOContext(io, :typeinfo => typ)
    show(ioc, np_unitless)
    printstyled(io, unit(eltype(etyp)); color=:cyan)
end

const AA = AbstractArray
function nest(x::AA, y::AA, z::Number)
    Npos(map(x, y) do x, y
            [x, y, z]
        end)
end
function nest(x::AA, y::Number, z::Number)
    Npos(map(x) do x
            [x, y, z]
        end)
end
function nest(x::Number, y::AA, z::Number)
    Npos(map(y) do y
            [x, y, z]
        end)
end
function nest(x::Number, y::Number, z::AA)
    Npos(map(z) do z
            [x, y, z]
        end)
end

function nest(x::AA, y::Number)
    Npos(map(x) do x
            [x, y]
        end)
end
function nest(x::Number, y::AA)
    Npos(map(y) do y
            [x, y]
        end)
end
function nest(x::AA, y::AA, z::AA, α::AA, β::AA, γ::AA)
    Npos(map(x, y, z, α, β, γ) do x, y, z, α, β, γ
            [x, y, z, α, β, γ]
        end)
end
function nest(x::AA, y::AA, z::AA, α::AA, β::AA)
    Npos(map(x, y, z, α, β) do x, y, z, α, β
            [x, y, z, α, β]
        end)
end
function nest(x::AA, y::AA, z::AA, α::AA)
    Npos(map(x, y, z, α, β) do x, y, z, α
            [x, y, z, α]
        end)
end
function nest(x::AA, y::AA, z::AA)
    Npos(map(x, y, z) do x, y, z
            [x, y, z]
        end)
end
function nest(x::AA, y::AA)
    Npos(map(x, y) do x, y
            [x, y]
        end)
end
function nest(x::AA)
    Npos(map(x) do x
            [x]
        end)
end
"Slice out 1st component from a nested vector or matrix"
slice_1(np::Npos) = map(p-> p[1], np)
"Slice out 2nd component from a nested vector or matrix"
slice_2(np::Npos) = map(p-> p[2], np)
"Slice out 3rd component from a nested vector or matrix"
slice_3(np::Npos) = map(p-> p[3], np)
"Slice out 4th component from a nested vector or matrix"
slice_4(np::Npos) = map(p-> p[4], np)
"Slice out 5th component from a nested vector or matrix"
slice_5(np::Npos) = map(p-> p[5], np)
"Slice out sixth component from a nested vector or matrix"
slice_6(np::Npos) = map(p-> p[6], np)

abstract type Body end
"""
Constructor defaults to a spherical earth. Alternatively, use
```
    earthlikeSphere(p = my_p, m = my_p)
    earthlikeSphere(p = my_p, r = my_r)
```
"""
Base.@kwdef struct Sphere <: Body
    p::Pos = [0., 0., 0.]m
    m::Number = m_earth
    r::Number = r_earth
    Sphere(p, m, r) = m > zero(m) && r > zero(r) ? new(p, m, r) : error("no negative m or r")
end

show(io::IO, b::Sphere) = print(io, "Sphere(p = ", b.p, ", m = ", b.m, ", r = ", b.r, ")")

"Alternative constructors of Sphere for when mass or radius is unknown"
earthlikeSphere_m_known(;p = [0, 0, 0]m, m = m_earth) = Sphere(p, m , (3m / (4π * ρ_earth))^(1/3))
earthlikeSphere_r_known(;p = [0, 0, 0]m, r = r_earth) = Sphere(p, ρ_earth * V(r), r)

const Bodies = Vector{T} where T <:Body
Bodies(bds::Bodies) = Vector(bds)
function bodies_3d()
    earth = Sphere()
    moon = earthlikeSphere_m_known(p = [d_moon_earth, 0m, 0m], m = m_moon)
    sun = Sphere([-d_sun_earth, 0m, 0m], m_sun, (3m_sun / (4π * ρ_sun))^(1/3))
    Bodies([sun, earth, moon])
end
function bodies_2d()
    earth = Sphere(p = [0,0]m)
    moon = earthlikeSphere_m_known(p = [d_moon_earth, 0m], m = m_moon)
    sun = Sphere([-d_sun_earth, 0m], m_sun, (3m_sun / (4π * ρ_sun))^(1/3))
    Bodies([sun, earth, moon])
end
function bodies_1d()
    earth = Sphere(p = [0]m)
    moon = earthlikeSphere_m_known(p = [d_moon_earth], m = m_moon)
    sun = Sphere([-d_sun_earth], m_sun, (3m_sun / (4π * ρ_sun))^(1/3))
    Bodies([sun, earth, moon])
end

function show(io::IO, ::MIME"text/plain", bds::Bodies) # long form.
    # show more descriptively, with one line per body. Based on Dict.
    recur_io = IOContext(io, :SHOWN_SET => bds)
    limit::Bool = get(io, :limit, false)
    if !haskey(io, :compact)
        recur_io = IOContext(recur_io, :compact => true)
    end
    summary(io, bds)
    isempty(bds) && return
    print(io, ":")
    show_circular(io, bds) && return
    if limit
        sz = displaysize(io)
        rows, cols = sz[1] - 3, sz[2]
        rows < 2   && (print(io, " …"); return)
        cols < 12  && (cols = 12) # Minimum widths of 2 for key, 4 for value
        cols -= 6 # Subtract the widths of prefix "  " separator " => "
        rows -= 1 # Subtract the summary

        # determine max key width to align the output, caching the strings
        ps = Vector{AbstractString}(undef, min(rows, length(bds)))
        ms = Vector{AbstractString}(undef, min(rows, length(bds)))
        rs = Vector{AbstractString}(undef, min(rows, length(bds)))
        plen = 0
        mlen = 0
        rlen = 0
        for (i, b) in enumerate(bds)
            i > rows && break
            st = sprint(show, b.p, context=recur_io, sizehint=0)
            ps[i] = "p = " * st
            st = sprint(show,  b.m, context=recur_io, sizehint=0)
            ms[i] = " m = " * st
            st = sprint(show, b.r, context=recur_io, sizehint=0)
            rs[i] = " r = " * st
            plen = clamp(length(ps[i]), plen, cols)
            mlen = clamp(length(ms[i]), mlen, cols)
            rlen = clamp(length(rs[i]), rlen, cols)
        end
        if plen > max(div(cols, 2), cols - mlen - rlen)
            plen = max(cld(cols, 3), cols - mlen - rlen)
        end
    else
        rows = cols = typemax(Int)
    end

    for (i, bd) in enumerate(bds)
        print(io, "\n  ")
        i == rows < length(bds) && (print(io, rpad("⋮", plen), ", ⋮"); break)
        if limit
            pos = rpad(_truncate_at_width_or_chars(ps[i], plen, "\r\n"), plen)
        else
            st = sprint(show, bd.p, context=recur_io, sizehint=0)
            pos = "p = " * st
        end
        print(recur_io, pos)
        if limit
            val = rpad(_truncate_at_width_or_chars(ms[i], cols - plen, "\r\n"), mlen)
            print(io, val)
            val = _truncate_at_width_or_chars(rs[i], cols - plen - mlen, "\r\n")
            print(io, val)
        else
            st = " m = " * sprint(show, bd.m, context=recur_io, sizehint=0) *
                " r = " * sprint(show, bd.r, context=recur_io, sizehint=0)
            print(recur_io, st)
        end
    end
end

"Vector from position (of) a to position (of) b."
rvec(a::T, b::T) where T<: Body = b.p - a.p
rvec(a::Pos, b::Body) = b.p - a
rvec(a::Body, b::Pos) = b - a.p
rvec(a::T, b::T) where T<: Pos = b - a
"Unity vector of vector"
uvec(p::Pos) = len(p) > zero(eltype(p)) ? p / len(p) : ustrip(zero(p))
"Unity vector from position (of) a to position (of) b."
uvec(a::T, b::T) where T<: Body = uvec(b.p - a.p)
uvec(a::Pos, b::Body) = uvec(b.p - a)
uvec(a::Body, b::Pos) = uvec(b - a.p)
uvec(a::T, b::T) where T<: Pos = uvec(b - a)

"Euclidean vector length"
len(p) =  sqrt(sum(p.*p))
"Euclidean length between positions (of) a and b"
len(a, b) = len(rvec(a, b))

"Gravitational acceleration magnitude"
g(m::Number, r::Number) = r > zero(r) ? γ * m / r^2 : zero(γ) / one(r)^2  # Test thoroughly for overflow

"Gravity vector on thisbody or position due to otherbody or bodies"
gvec(p::Pos, body::Body) = uvec(p, body) * g(body.m, len(p, body))
gvec(thisbody::Body, otherbody::Body) = gvec(thisbody.p, otherbody)
function gvec(p::Pos, sphere::Sphere)
    l = len(p, sphere)
    if l > sphere.r || l == 0.
        return uvec(p, sphere) * g(sphere.m, l)
    else
        # Inside sphere
        return uvec(p, sphere) * g(sphere.m, sphere.r) * l / sphere.r
    end
end
function gvec(p::Pos, bodies::Bodies)
    gvec_one = b -> gvec(p, b)
    sum(gvec_one.(bodies))
end
gvec(body::Body, bodies::Bodies) = gvec(body.p, bodies)
gvec(vp::Npos, bodies::Bodies) = map(p -> gvec(p, bodies), vp)

"Gravity field strength [length / time^2], nested position vectors"
gfieldvals(vp::Npos, bodies::Bodies) = map(p -> len(gvec(p, bodies)), vp)

"Return matrices mat_x, mat_y, size n * n"
function xygrids(;xmin = -3., xmax = 3., ymin, ymax, n = 100)
    mat_x = repeat(range(xmin, xmax, length = n)', n, 1)
    mat_y = repeat(range(ymin, ymax, length = n), 1, n)
    mat_x, mat_y
end
"100x100 grids around origo"
sqgrids(mi, ma) = xygrids( xmin = mi,xmax = ma, ymin = -(ma-mi)/2, ymax = (ma-mi)/2, n = 100)

nothing
end

max_stars_repo_name: hustf/Gravityfield

max_issues_repo_name: hustf/Gravityfield

max_forks_repo_name: hustf/Gravityfield