text:
# libpurex:
# sampling rules

# We organise them in two levels
# - sampling rule; a factory for sampling rule states
# - sampling rule state; keeps track of i.e. tracking information etc.

include("regret.jl");

# Tracking helper types

struct CTracking
    sumw;
    CTracking(N) = new(
        Float64.(N)
    ); # makes a copy of the starting situation
end

abbrev(_::Type{CTracking}) = "C";

# add a weight vector and track it
function track(t::CTracking, N, w)
    @assert all(N .≤ t.sumw.+1) "N $N  sumw $(t.sumw)";
    @assert sum(N) ≈ sum(t.sumw);

    t.sumw .+= w;
    argmin(N ./ t.sumw);
end

struct DTracking
    DTracking(N) = new();
end

abbrev(_::Type{DTracking}) = "D";

function track(t::DTracking, N, w)
    argmin(N .- sum(N).*w);
end

# Wrapper to add forced exploration to a tracking rule
struct ForcedExploration
    t;
end

function track(fe::ForcedExploration, N, w)
    t = sum(N);
    K = length(N);
    undersampled = N .≤ sqrt(t) .- K/2;
    if any(undersampled)
        track(fe.t, N, undersampled/sum(undersampled));
    else
        track(fe.t, N, w);
    end
end

# optimistic gradients
function optimistic_gradient(pep, hμ, t, N, λs)
    # TODO: there may need to be log(log(t)) terms here
    # TODO: is it correct to put uppers and lowers also through ξs? (consensus: hell NO)

    [let dist = getexpfam(pep, k),
     ↑ = dup(dist, hμ[k], log(t)/N[k]),
     ↓ = ddn(dist, hμ[k], log(t)/N[k])
     max(d(dist, ↑, λs[k]), d(dist, ↓, λs[k]), log(t)/N[k])
     end
     for k in eachindex(hμ)];
end

struct RoundRobin # used as factory and state
end

long(sr::RoundRobin) = "Uniform";
abbrev(sr::RoundRobin) = "RR";

function start(sr::RoundRobin, N)
    return sr;
end

function nextsample(sr::RoundRobin, pep, istar, ξ, N, S)
    return 1+(sum(N) % length(N));
end

struct FixedWeights # used as factory and state
    w;
    function FixedWeights(w)
        @assert all(w .≥ 0) && sum(w) ≈ 1 "$w not in simplex";
        new(w)
    end
end

long(sr::FixedWeights) = "Oracle Weigths";
abbrev(sr::FixedWeights) = "opt";

function start(sr::FixedWeights, N)
    return sr;
end

function nextsample(sr::FixedWeights, pep, istar, ξ, N, S)
    argmin(N .- sum(N).*sr.w);
end

struct FictitiousPlay
end

long(sr::FictitiousPlay) = "Fictitious Play";
abbrev(sr::FictitiousPlay) = "FP";

struct FictitiousPlayState
    sum∇;
    FictitiousPlayState(K) = new(zeros(K));
end

function start(sr::FictitiousPlay, N)
    FictitiousPlayState(length(N));
end

function nextsample(sr::FictitiousPlayState, pep, istar, ξ, N, S)
    K = length(N);
    t = sum(N);
    hμ = S./N; # emp. estimates

    # best response λ-player to all past
    _, (_, λs), (_, ξs) = glrt(pep, N, hμ);

    #println("N =$N");
    #println("hμ=$hμ");
    #println("λs=$λs");
    #println("ξs=$ξs");

    ∇ = optimistic_gradient(pep, hμ, t, N, λs);
    sr.sum∇ .+= ∇;

    #println("$t, $N because $∇ summed to $(sr.sum∇)");

    argmax(sr.sum∇); # best response for k player to all past
end

struct TrackAndStop
    TrackingRule;
end

long(sr::TrackAndStop) = "Track-and-Stop " * abbrev(sr.TrackingRule);
abbrev(sr::TrackAndStop) = "T-" * abbrev(sr.TrackingRule);

struct TrackAndStopState
    t;
    TrackAndStopState(TrackingRule, N) = new(
        ForcedExploration(TrackingRule(N))
    );
end

function start(sr::TrackAndStop, N)
    TrackAndStopState(sr.TrackingRule, N);
end

function nextsample(sr::TrackAndStopState, pep, istar, ξ, N, S)
    K = length(N);
    t = sum(N);

    # oracle at ξ (the closest feasible bandit model)
    _, w = oracle(pep, ξ);

    # tracking
    return track(sr.t, N, w);
end

struct OptimisticTrackAndStop
    TrackingRule;
end

long(sr::OptimisticTrackAndStop) = "Optimistic TaS " * abbrev(sr.TrackingRule);
abbrev(sr::OptimisticTrackAndStop) = "O-" * abbrev(sr.TrackingRule);

struct OptimisticTrackAndStopState
    t;
    OptimisticTrackAndStopState(TrackingRule, N) = new(
        TrackingRule(N)
    );
end

function start(sr::OptimisticTrackAndStop, N)
    OptimisticTrackAndStopState(sr.TrackingRule, N);
end

function nextsample(sr::OptimisticTrackAndStopState, pep, istar, ξ, N, S)
    K = length(N);
    t = sum(N);

    # optimistic oracle at ξ (the closest feasible bandit model)
    # TODO: think about this ξ, should it be \hat μ? Or what?
    _, w = optimistic_oracle(pep, ξ, N);

    # tracking
    return track(sr.t, N, w);
end

struct Menard
    TrackingRule;
    scale;
end

long(sr::Menard) = "Menard " * abbrev(sr.TrackingRule);
abbrev(sr::Menard) = "M-" * abbrev(sr.TrackingRule);

function start(sr::Menard, N)
    MenardState(sr.TrackingRule, sr.scale, N);
end

struct MenardState
    h;
    t;
    MenardState(TrackingRule, scale, N) = new(
        FixedShare(length(N), S=scale),
        TrackingRule(N)
    );
end

function nextsample(sr::MenardState, pep, istar, ξ, N, S)
    K = length(N);
    t = sum(N);
    hμ = S./N; # emp. estimates

    # query the learner
    w = act(sr.h);

    # best response λ-player
    _, (_, λs), (_, _) = glrt(pep, w, hμ);

    # gradient
    ∇ = [d(getexpfam(pep, k), hμ[k], λs[k]) for k in eachindex(hμ)];

    # update learner
    incur!(sr.h, -∇);

    # tracking
    return track(sr.t, N, w);
end

struct UnaBomb
    TrackingRule;
end

long(sr::UnaBomb) = "UnaBomb " * abbrev(sr.TrackingRule);
abbrev(sr::UnaBomb) = "U-" * abbrev(sr.TrackingRule);

struct UnaBombState
    h; # one online learner in total
    t;
    UnaBombState(TrackingRule, N) = new(
        AdaHedge(length(N)),
        TrackingRule(N)
    );
end

function start(sr::UnaBomb, N)
    UnaBombState(sr.TrackingRule, N);
end

function nextsample(sr::UnaBombState, pep, istar, ξ, N, S)
    K = length(N);

    t = sum(N);

    hμ = S./N; # emp. estimates
    #println("hμ=$hμ");

    # query the learner
    w = act(sr.h);

    # best response λ-player
    _, (_, λs), (_, ξs) = glrt(pep, w, hμ);

    #println("λs=$λs");
    #println("ξs=$ξs");

    #println("↑s=$([dup(dist, hμ[k], log(t)/N[k]) for k in 1:K])");
    #println("↓s=$([ddn(dist, hμ[k], log(t)/N[k]) for k in 1:K])");

    # feed linear loss plus optimism back to learner
    # println("d=$([d(dist, hμ[k], λs[k]) for k in 1:K])");

    # TODO: there may need to be log(log(t)) terms here
    # TODO: is it correct to put uppers and lowers also through ξs?
    ∇ = optimistic_gradient(pep, hμ, t, N, λs);
    #println("∇=$∇");
    incur!(sr.h, -∇);

    # tracking
    return track(sr.t, N, w);
end

struct DaBomb
    TrackingRule;
    M;
end

long(sr::DaBomb) = "DaBomb " * abbrev(sr.TrackingRule);
abbrev(sr::DaBomb) = "D-" * abbrev(sr.TrackingRule);

struct DaBombState
    hs; # one online learner per answer
    t;
    DaBombState(TrackingRule, N, M) = new(
        map(x -> AdaHedge(length(N)), 1:M),
        TrackingRule(N)
    );
end

function start(sr::DaBomb, N)
    DaBombState(sr.TrackingRule, N, sr.M);
end

function nextsample(sr::DaBombState, pep, istar, ξ, N, S)
    K = length(N);

    t = sum(N);

    hμ = S./N; # emp. estimates
    #println("hμ=$hμ");

    # query the learner
    w = act(sr.hs[istar]);

    # best response λ-player
    _, (_, λs), (_, ξs) = glrt(pep, w, hμ);

    #println("λs=$λs");
    #println("ξs=$ξs");

    #println("↑s=$([dup(dist, hμ[k], log(t)/N[k]) for k in 1:K])");
    #println("↓s=$([ddn(dist, hμ[k], log(t)/N[k]) for k in 1:K])");

    # feed linear loss plus optimism back to learner
    # println("d=$([d(dist, hμ[k], λs[k]) for k in 1:K])");

    ∇ = optimistic_gradient(pep, hμ, t, N, λs)

    #println("∇=$∇");
    incur!(sr.hs[istar], -∇);

    # tracking
    return track(sr.t, N, w);
end


text:
#
# Input types
#
abstract type ABCInput end

abstract type ABCRejectionInput <: ABCInput end

RepetitiveTraining(; rt_iterations::Int64=0, rt_extra_training_points::Int64=1, rt_sample_size::Int64=1000) =
    RepetitiveTraining(rt_iterations, rt_extra_training_points, rt_sample_size)

struct DistanceSimulationInput
    reference_summary_statistic::AbstractArray{Float64,1}
    simulator_function::Function
    summary_statistic::Function
    distance_metric::Function
end

"""
	AbstractEmulatorRetraining

Subtypes of this abstract type control the additional retraining procedure that may or may not be carried out before each iteration of emulation-based ABC. Custom strategies may be implemented by creating new subtypes of this type and new [`abc_retrain_emulator`](@ref) methods for them.

The following implementations are shipped:
- [`IncrementalRetraining`](@ref)
- [`PreviousPopulationRetraining`](@ref)
- [`PreviousPopulationThresholdRetraining`](@ref)
- [`NoopRetraining`](@ref)
"""
abstract type AbstractEmulatorRetraining end

"""
	IncrementalRetraining <: AbstractEmulatorRetraining

This emulator retraining strategy samples extra particles from the previous population, and adds them to the set of design points that were used on the previous iteration. The new design points are filtered according to the threshold. The combined set of design points is used to train the new emulator.

# Fields
- `design_points`: number of design points to add on each iteration
- `max_simulations`: maximum number of simulations to perform during re-trainging on each iteration
"""
struct IncrementalRetraining <: AbstractEmulatorRetraining
    design_points::Int
    max_simulations::Int
end

"""
	PreviousPopulationRetraining <: AbstractEmulatorRetraining

This emulator retraining strategy samples extra particles from the previous population, and uses them to re-train the emulator from scratch. No filtering of the new design points is performed. Design points from the previous iteration are discarded.
"""
struct PreviousPopulationRetraining <: AbstractEmulatorRetraining end

"""
	PreviousPopulationThresholdRetraining <: AbstractEmulatorRetraining

This emulator retraining strategy samples extra particles from the previous population, and uses them to re-train the emulator from scratch. Design points from the previous iteration are discarded. This strategy allows to control how many design points are sampled with distance to the reference data below the threshold.

# Fields:
- `n_design_points`: number of design points
- `n_below_threshold`: number of design points below the threshold
- `max_iter`: maximum number of simulations to perform on each re-training iteration
"""
struct PreviousPopulationThresholdRetraining <: AbstractEmulatorRetraining
    n_design_points::Int
    n_below_threshold::Int
    max_iter::Int
end

"""
	NoopRetraining <: AbstractEmulatorRetraining

A sentinel retraining strategy that does not do anything. When used, the emulator is trained only once at the start of the process.
"""
struct NoopRetraining <: AbstractEmulatorRetraining end

"""
    AbstractEmulatorTraining

Subtypes of this abstract type control how the emulator is trained for emulation-based ABC algorithms (rejection and SMC).
At the moment, only [`DefaultEmulatorTraining`](@ref) is shipped. Custom emulator training procedure
can be implemented by creating new subtypes of this type and overriding [`train_emulator`](@ref) for them.

A typical use case would be trying to control the behaviour of [`gp_train`](@ref) more tightly,
or not using it altogeather (e.g. using another optimisation package).
"""
abstract type AbstractEmulatorTraining end

"""
	DefaultEmulatorTraining <: AbstractEmulatorTraining

# Fields
- `kernel::AbstractGPKernel`: the kernel ([`AbstractGPKernel`](@ref)) that will be used with the Gaussian Process ([`GPModel`](@ref)). Defaults to [`SquaredExponentialArdKernel`](@ref).

[`train_emulator`](@ref) method with this argument type calls [`gp_train`](@ref) with default arguments.
"""
struct DefaultEmulatorTraining{K<:AbstractGPKernel} <: AbstractEmulatorTraining
    kernel::K
end
DefaultEmulatorTraining() = DefaultEmulatorTraining(SquaredExponentialArdKernel())

struct EmulatorTrainingInput{ET<:AbstractEmulatorTraining}
    distance_simulation_input::DistanceSimulationInput
    design_points::Int64
    emulator_training::ET
end
EmulatorTrainingInput(dsi::DistanceSimulationInput) = EmulatorTrainingInput(dsi, DefaultEmulatorTraining())
EmulatorTrainingInput(n_design_points, reference_summary_statistic, simulator_function, summary_statistic, distance_metric, et=DefaultEmulatorTraining()) =
    EmulatorTrainingInput(DistanceSimulationInput(
        reference_summary_statistic, simulator_function,
        build_summary_statistic(summary_statistic), distance_metric),
        n_design_points, et)

"""
    AbstractEmulatedParticleSelection

Subtypes of this type control the criteria that determine what particles are included in the posterior for emulation-based ABC. Custom strategies
can be implemented by creating new subtypes of this type and overriding [`abc_select_emulated_particles`](@ref) for them.

Three implementations are shipped:
- [`MeanEmulatedParticleSelection`](@ref)
- [`MeanVarEmulatedParticleSelection`](@ref)
- [`PosteriorSampledEmulatedParticleSelection`](@ref)
"""
abstract type AbstractEmulatedParticleSelection end

"""
    MeanEmulatedParticleSelection <: AbstractEmulatedParticleSelection

When this strategy is used, the particles for which only the *mean* value returned by
[`gp_regression`](@ref) is below the ABC threshold are included in the posterior.
Variance is not checked.
"""
struct MeanEmulatedParticleSelection <: AbstractEmulatedParticleSelection end

"""
    MeanVarEmulatedParticleSelection <: AbstractEmulatedParticleSelection

When this strategy is used, the particles for which both *mean and standard deviation* returned by
[`gp_regression`](@ref) is below the ABC threshold are included in the posterior. In pseudocode:
```julia
means, vars = gp_regression(parameters, gpm)
accepted_indices = findall((means .<= threshold) .& (sqrt.(vars) .<= threshold))
```

The rationale behind using this selection strategy is to take into account the "level of uncertainty"
about the regression prediction that is provided by the Gaussian Process in form of standard deviation.
So, even if the mean of the GP is below the threshold, but the GP is "uncertain" about it
(i.e. the variance is high), this particle will not be included in the posterior distribution of ABC.
It is a more stringent acceptance criteria than `MeanEmulatedParticleSelection`.

# Fields
- `variance_threshold_factor`: scaling factor, by which the ABC threshold is multiplied
before checking the standard deviation. Defaults to 1.0.
"""
struct MeanVarEmulatedParticleSelection <: AbstractEmulatedParticleSelection
    variance_threshold_factor::Float64
end
MeanVarEmulatedParticleSelection() = MeanVarEmulatedParticleSelection(1.0)

"""
    PosteriorSampledEmulatedParticleSelection <: AbstractEmulatedParticleSelection

When this strategy is used, the distance is sampled from the GP posterior of the [`gp_regression`](@ref)
object. If the sampled distance is below the threshold the particle is accepted.

# Fields
- `use_diagonal_covariance`: if `true`, the GP posterior covariance will be approximated by its
diagonal elements only. Defaults to `false`.
"""
struct PosteriorSampledEmulatedParticleSelection <: AbstractEmulatedParticleSelection
    use_diagonal_covariance::Bool
end
PosteriorSampledEmulatedParticleSelection() = PosteriorSampledEmulatedParticleSelection(false)

struct SimulatedABCRejectionInput <: ABCRejectionInput
    n_params::Int64
    n_particles::Int64
    threshold::Float64
    priors::AbstractArray{ContinuousUnivariateDistribution,1}
    distance_simulation_input::DistanceSimulationInput
    max_iter::Int
end

struct EmulatedABCRejectionInput{CUD<:ContinuousUnivariateDistribution, EPS<:AbstractEmulatedParticleSelection} <: ABCRejectionInput
	n_params::Int64
	n_particles::Int64
	threshold::Float64
	priors::AbstractArray{CUD,1}
	batch_size::Int64
    max_iter::Int64
    emulator_training_input::EmulatorTrainingInput
    selection::EPS
end

abstract type ABCSMCInput <: ABCInput end

struct SimulatedABCSMCInput <: ABCSMCInput
    n_params::Int64
    n_particles::Int64
    threshold_schedule::AbstractArray{Float64,1}
    priors::AbstractArray{ContinuousUnivariateDistribution,1}
    distance_simulation_input::DistanceSimulationInput
    max_iter::Int
end

SimulatedABCRejectionInput(smc_input::SimulatedABCSMCInput) =
    SimulatedABCRejectionInput(smc_input.n_params,
        smc_input.n_particles,
        smc_input.threshold_schedule[1],
        smc_input.priors,
        smc_input.distance_simulation_input,
        smc_input.max_iter)

struct EmulatedABCSMCInput{CUD<:ContinuousUnivariateDistribution,
        ER<:AbstractEmulatorRetraining, EPS<:AbstractEmulatedParticleSelection} <: ABCSMCInput
    n_params::Int64
    n_particles::Int64
    threshold_schedule::AbstractArray{Float64,1}
    priors::AbstractArray{CUD,1}
    batch_size::Int64
    max_iter::Int64
    emulator_training_input::EmulatorTrainingInput
    emulator_retraining::ER
    selection::EPS
end

#
# Tracker types
#
abstract type ABCSMCTracker end

mutable struct SimulatedABCSMCTracker <: ABCSMCTracker
    n_params::Int64
    n_accepted::AbstractArray{Int64,1}
    n_tries::AbstractArray{Int64,1}
    threshold_schedule::AbstractArray{Float64,1}
	can_continue::Bool # we used to return this flag from iterateABCSMC!, but it got broken in julia 1.0 - see https://github.com/JuliaLang/julia/issues/29805
    population::AbstractArray{AbstractArray{Float64,2},1}
    distances::AbstractArray{AbstractArray{Float64,1},1}
    weights::AbstractArray{StatsBase.Weights,1}
    priors::AbstractArray{ContinuousUnivariateDistribution,1}
    distance_simulation_input::DistanceSimulationInput
    max_iter::Int64
end

mutable struct EmulatedABCSMCTracker{CUD<:ContinuousUnivariateDistribution, ET,
        ER<:AbstractEmulatorRetraining, EPS<:AbstractEmulatedParticleSelection} <: ABCSMCTracker
    n_params::Int64
    n_accepted::AbstractArray{Int64,1}
    n_tries::AbstractArray{Int64,1}
    threshold_schedule::AbstractArray{Float64,1}
	can_continue::Bool # we used to return this flag from iterateABCSMC!, but it got broken in julia 1.0 - see https://github.com/JuliaLang/julia/issues/29805
    population::AbstractArray{AbstractArray{Float64,2},1}
    distances::AbstractArray{AbstractArray{Float64,1},1}
    weights::AbstractArray{StatsBase.Weights,1}
    priors::AbstractArray{CUD,1}
    emulator_training_input::EmulatorTrainingInput
    emulator_retraining::ER
    selection::EPS
    batch_size::Int64
    max_iter::Int64
    emulators::AbstractArray{ET,1}
end

#
# Output types
#
abstract type ABCOutput end

"""
    ABCRejectionOutput

A container for the output of a rejection-ABC computation.

# Fields
- `n_params::Int64`: The number of parameters to be estimated.
- `n_accepted::Int64`: The number of accepted parameter vectors (particles) in the posterior.
- `n_tries::Int64`: The total number of parameter vectors (particles) that were tried.
- `threshold::Float64`: The maximum distance from the summarised model output to summarised observed data for a parameter vector to be included in the posterior.
- `population::AbstractArray{Float64,2}`: The parameter vectors (particles) in the posterior. Size: (`n_accepted`, `n_params`).
- `distances::AbstractArray{Float64,1}`: The distances for each parameter vector (particle) in the posterior to the observed data in summary statistic space. Size: (`n_accepted`).
- `weights::StatsBase.Weights`: The weight of each parameter vector (particle) in the posterior.
"""
abstract type ABCRejectionOutput <: ABCOutput end

struct EmulatedABCRejectionOutput{ET} <: ABCRejectionOutput
    n_params::Int64
    n_accepted::Int64
    n_tries::Int64
    threshold::Float64
    population::AbstractArray{Float64,2}
    distances::AbstractArray{Float64,1}
    weights::StatsBase.Weights
    emulator::ET
end

struct SimulatedABCRejectionOutput <: ABCRejectionOutput
    n_params::Int64
    n_accepted::Int64
    n_tries::Int64
    threshold::Float64
    population::AbstractArray{Float64,2}
    distances::AbstractArray{Float64,1}
    weights::StatsBase.Weights
end

"""
    ABCSMCOutput

A container for the output of a rejection-ABC computation.

# Fields
- `n_params::Int64`: The number of parameters to be estimated.
- `n_accepted::Int64`: The number of accepted parameter vectors (particles) in the posterior.
- `n_tries::Int64`: The total number of parameter vectors (particles) that were tried.
- `threshold_schedule::AbstractArray{Float64}`: A set of maximum distances from the summarised model output to summarised observed data for a parameter vector to be included in the posterior. Each distance will be used in a single run of the ABC-SMC algorithm.
- `population::AbstractArray{Float64,2}`: The parameter vectors (particles) in the posterior. Size: (`n_accepted`, `n_params`).
- `distances::AbstractArray{Float64,1}`: The distances for each parameter vector (particle) in the posterior to the observed data in summary statistic space. Size: (`n_accepted`).
- `weights::StatsBase.Weights`: The weight of each parameter vector (particle) in the posterior.
"""
abstract type ABCSMCOutput <: ABCOutput end

struct SimulatedABCSMCOutput <: ABCSMCOutput
    n_params::Int64
    n_accepted::AbstractArray{Int64,1}
    n_tries::AbstractArray{Int64,1}
    threshold_schedule::AbstractArray{Float64,1}
    population::AbstractArray{AbstractArray{Float64,2},1}
    distances::AbstractArray{AbstractArray{Float64,1},1}
    weights::AbstractArray{StatsBase.Weights,1}
end

struct EmulatedABCSMCOutput <: ABCSMCOutput
    n_params::Int64
    n_accepted::AbstractArray{Int64,1}
    n_tries::AbstractArray{Int64,1}
    threshold_schedule::AbstractArray{Float64,1}
    population::AbstractArray{AbstractArray{Float64,2},1}
    distances::AbstractArray{AbstractArray{Float64,1},1}
    weights::AbstractArray{StatsBase.Weights,1}
    emulators::AbstractArray{GPModel,1}
end

function checkABCInput(input::ABCInput)
	if input.n_params != length(input.priors)
		throw(ArgumentError("There are $(input.n_params) unknown parameters but $(length(input.priors)) priors were provided"))
	end
end

ext: jl


text:
function parse_n_walkers(n)
    nums = map(x -> parse(Int, x), split(n, ':'))
    if length(nums) == 1
        return nums[1]:nums[1]
    else
        return Colon()(nums...)
    end
end

function parse_style(s)
    return if startswith(s, "int")
        IsStochasticInteger()
    elseif startswith(s, "semi")
        IsDynamicSemistochastic()
    elseif startswith(s, "explosive")
        IsExplosive()
    elseif startswith(s, "double")
        IsDynamicSemistochastic(compression=Rimu.StochasticStyles.DoubleOrNothingWithThreshold())
    elseif startswith(s, "without_r")
        IsDynamicSemistochastic(spawning=Rimu.StochasticStyles.WithoutReplacement())
    elseif startswith(s, "bernoulli")
        IsDynamicSemistochastic(spawning=Rimu.StochasticStyles.Bernoulli())
    else
        error("Unknown style '", s, "'")
    end
end

function unpack_commandline_args(args)
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--num_walkers"
        help="Number of walkers to use."
        required=true

        "--id"
        help="The id of the computation."
        required=true

        "--style"
        help="Stochastic style to use."
        default="semi"

        "--dir"
        help="Output directory."
        default="."

        "--initiator"
        help="Use initiator?"
        arg_type=Bool
        default=false

        "--steps"
        help="Record this many steps."
        arg_type=Int
        default=60_000

        "--dt"
        help="Timestep size."
        arg_type=Float64
        default=1e-3

        "--continuation"
        help="Do continuation runs between walker numbers?"
        arg_type=Bool
        default=true

        "--warmup"
        help="Warm up on the first walker number for this many steps. Only applicable if \"--continuation=true\""
        arg_type=Int
        default=10_000
    end
    res = Dict{Symbol,Any}()
    for (k, v) in parse_args(args, s)
        res[Symbol(k)] = v
    end
    return res
end

"""
    plateau(ham; kwargs...)
    plateau(ham; ARGS)

Perform FCIQMC for various target walkers and save all results to a directory of files.
Uses continuations to reduce equilibration times.

# Arguments:

* `num_walkers` (required): range of walkers to run the computation with.
* `id` (required): the id. A directory with that name will be created.
* `style`: the stochastic style ("int", "semi", or "early").
* `dir="."`: Output directory.
* `initiator=false`: Use initiators?
* `steps=60_000`: Record this many steps.
* `dt=1e-3`: Timestep size.
* `continuation`: Do continuation runs between walker numbers?
* `warmup=10_000`: Do this many steps on the first walker number before measuring.

Arguments can also be passed in an array of strings. Example usage:

```julia
# script.jl
using Sekira
H = HubbardMom1D(BoseFS((0,0,0,7,0,0,0)); u=5)
Sekira.plateau(H, ARGS)
```

```
\$ mpirun julia script.jl --num_walkers=100:100:500 --style="int" --id="test"
```

This will create a directory `test` where output files are saved.
"""
plateau(ham, args; kwargs...) = plateau(ham; unpack_commandline_args(args)..., kwargs...)
function plateau(
    ham; num_walkers, style, id,
    dir=".",
    initiator=false,
    steps=60_000,
    dt=1e-3,
    continuation=true,
    warmup=10_000,
    kwargs...,
)
    dvec_type = initiator ? InitiatorDVec : DVec
    style = parse_style(style)
    num_walkers = parse_n_walkers(num_walkers)
    return plateau(
        ham, num_walkers, style, id, dir, dvec_type, steps, dt, continuation, warmup;
        kwargs...
    )
end

function plateau(
    ham, num_walkers, style, id, dir, dvec_type, steps, dτ, continuation, warmup; kwargs...
)
    @mpi_root begin
        @info "Run started"
        root = initialize_run(dir, id)
    end

    dv = MPIData(dvec_type(starting_address(ham) => 10; style))

    # Set up.
    s_strat = DoubleLogUpdate(targetwalkers=num_walkers[1])
    maxlength = 2 * maximum(num_walkers)

    if continuation
        @mpi_root @info "Warming up."
        params = RunTillLastStep(; laststep=warmup, dτ)
        lomc!(ham, dv; s_strat, params, maxlength, kwargs...)
    end
    prev_file = "__warmup__"

    for n_target in num_walkers
        @mpi_root @info "Computing targetwalkers=$n_target"
        s_strat = DoubleLogUpdate(targetwalkers=n_target)
        if continuation && params.step == params.laststep
            # Only do continuation if previous run did not overflow maxlength.
            params.step = 0
            params.laststep = steps
        else
            empty!(localpart(dv))
            localpart(dv)[starting_address(ham)] = 10
            params = RunTillLastStep(; laststep=steps, dτ)
        end
        time = @elapsed df = lomc!(ham, dv; s_strat, params, maxlength, kwargs...).df

        @mpi_root begin
            schema = Tables.schema(df)
            style = StochasticStyle(localpart(dv))
            dvec_type = typeof(localpart(dv))
            continued = continuation ? prev_file : nothing
            metadata = (
                ; ham, n_target, steps, time, schema, style, dvec_type, params, continued,
            )
            @info "Done in $time seconds."
            show_metadata(stdout, metadata)
            save_run(root, "$n_target", df, metadata)
            prev_file = "$n_target"
        end
    end
end

function show_metadata(io::IO, metadata)
    for (k, v) in pairs(metadata)
        println(io, rpad(k, 10), " : ", v)
    end
end


text:
#==============================================================================
    Transpose 5 rows of DataFrame for easy viewing

    TODO:       kwarg to randomly select rows

    author:     mahiki@users.noreply.github.com
=============================================================================#

using DataFrames, Dates

"""
    df_transpose(df::DataFrame; n=5)

Take the first n rows of dataframe and transpose to read columns sideways, useful for viewing wide dataframes

# Examples
```jldoctest
julia> df = DataFrame(a = 1:2, b = [1.0, π], c = [Date("20191002", "yyyymmdd"), Date("21120101", "yyyymmdd")], d = [1//2, missing]);
julia> df_transpose(df)
4×3 DataFrame
 Row │ variable  1           2
     │ String    Any         Any
─────┼──────────────────────────────────
   1 │ a         1           2
   2 │ b         1.0         3.14159
   3 │ c         2019-10-02  2112-01-01
   4 │ d         1//2        missing
```
"""
function df_transpose(df::DataFrame; n=5)
    dfn = df[1:n,:]
    colnames = names(dfn)
    dfn[!, :id] = 1:size(dfn, 1)
    dfi = stack(dfn, colnames)
    unstack(dfi, :variable, :id, :value)
end

max_stars_repo_name: mahiki/DesertIslandDisk

max_issues_repo_name: mahiki/DesertIslandDisk

max_forks_repo_name: mahiki/DesertIslandDisk


text:
@testset "Mapping" begin
  data1D = readgeotable(joinpath(datadir,"data1D.tsv"), delim='\t', coordnames=(:x,))
  data2D = readgeotable(joinpath(datadir,"data2D.tsv"), delim='\t', coordnames=(:x,:y))

  @testset "NearestMapping" begin
    grid1D = RegularGrid(100)
    mappings = map(data1D, grid1D, (:value,), NearestMapping())
    @test mappings[:value] == Dict(100=>11,81=>9,11=>2,21=>3,91=>10,51=>6,61=>7,71=>8,31=>4,41=>5,1=>1)

    grid2D = RegularGrid(100,100)
    mappings = map(data2D, grid2D, (:value,), NearestMapping())
    @test mappings[:value] == Dict(5076=>3,2526=>1,7551=>2)

    ps2D = PointSet([25. 50. 75.; 25. 75. 50.])
    mappings = map(data2D, ps2D, (:value,), NearestMapping())
    @test mappings[:value] == Dict(2=>2,3=>3,1=>1)
  end

  @testset "CopyMapping" begin
    d = georef((z=rand(10),), rand(2,10))
    g = RegularGrid(10,10)

    # copy data to first locations in domain
    mappings = map(d, g, (:z,), CopyMapping())
    @test mappings[:z] == Dict(i=>i for i in 1:10)

    # copy data to last locations in domain
    mappings = map(d, g, (:z,), CopyMapping(91:100))
    @test mappings[:z] == Dict(j=>i for (i,j) in enumerate(91:100))

    # copy first 3 data points to last 3 domain locations
    mappings = map(d, g, (:z,), CopyMapping(1:3, 98:100))
    @test mappings[:z] == Dict(j=>i for (i,j) in enumerate(98:100))
  end
end


text:
# Functions for checking the validity of quantum elements such as states, unitarties,...

export isQuantumState, isHermitian, isUnitary, isPPT;

""" `b = isQuantumState(rho)` or `b = isQuantumState(rho,prec)`

Checks whether the input matrix rho is a valid quantum state, that is, a matrix that is positive semidefinite and has trace 1. *prec* is the precision that small eigenvalues are dealt with: If the input matrix has a negative eigenvalue that is no smaller than *-prec*, then this eigenvalue is considered zero.

Returns true/false.

"""
function isQuantumState(rho::AbstractMatrix, prec::Number=0.00001)

	# Check whether it's hermitian
	if !isHermitian(rho, prec)
		print("isQuantumState: Matrix is not Hermitian.\n")
		return false;
	end

	# check for positivity, real is to deal with small numerical errors
	l = real(eig(rho)[1]);
	if (minimum(l) < - prec)
		print("isQuantumState: Matrix is not positive semidefinite.\n")
		return false;
	end

	# Check the normalization
	if (real(trace(rho)) < 1 - prec) || (real(trace(rho)) > 1 + prec)
		print("isQuantumState: Matrix does not have trace one.\n");
		return false;
	end

	return true;
end

""" `b = isHermitian(rho)` or `b = is Hermitian(rho, prec)`

Checks whether the input matrix rho is a Hermitian matrix. Entries are compares up to precision *prec* on average.

Returns true/false.
"""

function isHermitian(rho::AbstractMatrix, prec::Number=0.00001)

	# check whether the matrix is equal to its conjugate transpose
	(d1,d2) = size(rho);
	diff = sum(rho - rho');
	if(real(diff) + imag(diff) > prec * d1)
		return false;
	end
	return true;
end

""" `b = isUnitary(U)`

Checks whether the input matrix is unitary.

Returns true/false.

"""

function isUnitary(U::AbstractMatrix)

	# Check dimensions
	(d,da) = size(U);
	if (d != da)
		print("isUnitary: Input is not a square matrix.\n");
		return false;
	end

	# Check unitarity
	if (round(U*U',6) != eye(d))
		return false;
	end
	if (round(U'*U,6) != eye(d))
		return false;
	end

	return true;
end

""" `b = isPPT(rhoAB, nA, nB)`

Checks whether the input *rhoAB* is PPT across dimension *nA* and *nB*.

Returns true/false
"""

function isPPT(rho::AbstractMatrix,nA::Number,nB::Number)

	# defined the cutoff when we consider something to be positive
	epsilon = 10.0^(-15);

	# Check whether it's hermitian
	@assert isQuantumState(rho) "Input is not a quantum state."

	# Check whether dimensions match
	(da, db) = size(rho)
	@assert (da ==  nA * nB) "Input does not match given dimensions."

	rhoPT = partialtranspose(rho, 2, [nA, nB]);

	# check for positivity
	l = eig(rhoPT)[1]
	if (minimum(l) < - epsilon)
		return false
	end

	return true
end


text:
country() = executor(data["en"]["faker"]["address"]["country"])
city_prefix() = executor(data["en"]["faker"]["address"]["city_prefix"])
city_suffix()= executor(data["en"]["faker"]["address"]["city_suffix"])
street_prefix()= executor(data["en"]["faker"]["address"]["street_prefix"])
street_suffix()= executor(data["en"]["faker"]["address"]["street_suffix"])
secondary_address() = bothify(executor(data["en"]["faker"]["address"]["secondary_address"]))
state() = executor(data["en"]["faker"]["address"]["state"])
state_abbr() = executor(data["en"]["faker"]["address"]["state_abbr"])
municipality() = executor(data["en"]["faker"]["address"]["municipality"])
street_name() = executor(data["en"]["faker"]["address"]["street_name"])
street_address() = executor(data["en"]["faker"]["address"]["street_address"])
address() = ("$(street_address()) $(city()), $(state_abbr()) $(postcode())")
postcode() = bothify(executor(data["en"]["faker"]["address"]["postcode"]))
building_number() = bothify(executor(data["en"]["faker"]["address"]["building_number"]))
city() = executor(data["en"]["faker"]["address"]["city"])

function geo_coordinate(;center="None", radius=0.001)

        if center=="None"
            center=0
        end

        if center == 0
            return Float16(rand(-17.99 : 17.99))
        else
            return Float16(rand( center - radius : center + radius))
        end
end

latitude() = geo_coordinate() / 2.0
longitude() = geo_coordinate()/ 1.0


text:
import GLFW
using ModernGL

include("util.jl")

GLFW.Init()
window = GLFW.CreateWindow(500, 500, "Simple Example")
GLFW.MakeContextCurrent(window)

data = GLfloat[
    0.0, 0.5,
    0.5, -0.5,
    -0.5,-0.5
]

vao = glGenVertexArray()
glBindVertexArray(vao)
vbo = glGenBuffer()
glBindBuffer(GL_ARRAY_BUFFER, vbo)
glBufferData(GL_ARRAY_BUFFER, sizeof(data), data, GL_STATIC_DRAW)

const vsh = """
    #version 300 es
    precision highp float;
    in vec2 position;

    void main() {
        gl_Position = vec4(position, 0.0, 1.0);
    }
"""

const fsh = """
    #version 300 es
    precision highp float;
    out vec4 outColor;

    void main() {
        outColor = vec4(1.0, 1.0, 1.0, 1.0);
    }
"""

vertexShader = createShader(vsh, GL_VERTEX_SHADER)
fragmentShader = createShader(fsh, GL_FRAGMENT_SHADER)
program = createShaderProgram(vertexShader, fragmentShader)
glUseProgram(program)
positionAttribute = glGetAttribLocation(program, "position");
glEnableVertexAttribArray(positionAttribute)
glVertexAttribPointer(positionAttribute, 2, GL_FLOAT, GL_FALSE, 0, C_NULL)

while !GLFW.WindowShouldClose(window)
    glClearColor(0.0, 0.0, 0.0, 1.0)
    glClear(GL_COLOR_BUFFER_BIT)
    glDrawArrays(GL_TRIANGLES, 0, 3)
    GLFW.SwapBuffers(window)
    GLFW.PollEvents()
end

GLFW.Terminate()

max_stars_repo_name: richardanaya/julia-opengl

max_issues_repo_name: richardanaya/julia-opengl

max_forks_repo_name: richardanaya/julia-opengl


text:
function add_default_eop1980_history_loaders(::Type{FramesFactory}, arg0::JString, arg1::JString, arg2::JString, arg3::JString, arg4::JString)
    return jcall(FramesFactory, "addDefaultEOP1980HistoryLoaders", void, (JString, JString, JString, JString, JString), arg0, arg1, arg2, arg3, arg4)
end

function add_default_eop2000_history_loaders(::Type{FramesFactory}, arg0::JString, arg1::JString, arg2::JString, arg3::JString, arg4::JString)
    return jcall(FramesFactory, "addDefaultEOP2000HistoryLoaders", void, (JString, JString, JString, JString, JString), arg0, arg1, arg2, arg3, arg4)
end

function add_eop_history_loader(::Type{FramesFactory}, arg0::IERSConventions, arg1::EOPHistoryLoader)
    return jcall(FramesFactory, "addEOPHistoryLoader", void, (IERSConventions, EOPHistoryLoader), arg0, arg1)
end

function clear_eop_history_loaders(::Type{FramesFactory})
    return jcall(FramesFactory, "clearEOPHistoryLoaders", void, ())
end

function find_eop(::Type{FramesFactory}, arg0::Frame)
    return jcall(FramesFactory, "findEOP", EOPHistory, (Frame,), arg0)
end

function get_cirf(::Type{FramesFactory}, arg0::IERSConventions, arg1::jboolean)
    return jcall(FramesFactory, "getCIRF", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_ecliptic(::Type{FramesFactory}, arg0::IERSConventions)
    return jcall(FramesFactory, "getEcliptic", Frame, (IERSConventions,), arg0)
end

function get_eme2000(::Type{FramesFactory})
    return jcall(FramesFactory, "getEME2000", FactoryManagedFrame, ())
end

function get_eop_history(::Type{FramesFactory}, arg0::IERSConventions, arg1::jboolean)
    return jcall(FramesFactory, "getEOPHistory", EOPHistory, (IERSConventions, jboolean), arg0, arg1)
end

function get_frame(::Type{FramesFactory}, arg0::Predefined)
    return jcall(FramesFactory, "getFrame", Frame, (Predefined,), arg0)
end

function get_frames(::Type{FramesFactory})
    return jcall(FramesFactory, "getFrames", LazyLoadedFrames, ())
end

function get_gcrf(::Type{FramesFactory})
    return jcall(FramesFactory, "getGCRF", Frame, ())
end

function get_gtod(::Type{FramesFactory}, arg0::IERSConventions, arg1::jboolean)
    return jcall(FramesFactory, "getGTOD", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_gtod(::Type{FramesFactory}, arg0::jboolean)
    return jcall(FramesFactory, "getGTOD", FactoryManagedFrame, (jboolean,), arg0)
end

function get_icrf(::Type{FramesFactory})
    return jcall(FramesFactory, "getICRF", Frame, ())
end

function get_itrf(::Type{FramesFactory}, arg0::IERSConventions, arg1::jboolean)
    return jcall(FramesFactory, "getITRF", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_itrf(::Type{FramesFactory}, arg0::ITRFVersion, arg1::IERSConventions, arg2::jboolean)
    return jcall(FramesFactory, "getITRF", VersionedITRF, (ITRFVersion, IERSConventions, jboolean), arg0, arg1, arg2)
end

function get_itrf_equinox(::Type{FramesFactory}, arg0::IERSConventions, arg1::jboolean)
    return jcall(FramesFactory, "getITRFEquinox", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_mod(::Type{FramesFactory}, arg0::IERSConventions)
    return jcall(FramesFactory, "getMOD", FactoryManagedFrame, (IERSConventions,), arg0)
end

function get_mod(::Type{FramesFactory}, arg0::jboolean)
    return jcall(FramesFactory, "getMOD", FactoryManagedFrame, (jboolean,), arg0)
end

function get_non_interpolating_transform(::Type{FramesFactory}, arg0::Frame, arg1::Frame, arg2::AbsoluteDate)
    return jcall(FramesFactory, "getNonInterpolatingTransform", Transform, (Frame, Frame, AbsoluteDate), arg0, arg1, arg2)
end

function get_non_interpolating_transform(::Type{FramesFactory}, arg0::Frame, arg1::Frame, arg2::FieldAbsoluteDate)
    return jcall(FramesFactory, "getNonInterpolatingTransform", FieldTransform, (Frame, Frame, FieldAbsoluteDate), arg0, arg1, arg2)
end

function get_pz9011(::Type{FramesFactory}, arg0::IERSConventions, arg1::jboolean)
    return jcall(FramesFactory, "getPZ9011", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_teme(::Type{FramesFactory})
    return jcall(FramesFactory, "getTEME", FactoryManagedFrame, ())
end

function get_tirf(::Type{FramesFactory}, arg0::IERSConventions)
    return jcall(FramesFactory, "getTIRF", FactoryManagedFrame, (IERSConventions,), arg0)
end

function get_tirf(::Type{FramesFactory}, arg0::IERSConventions, arg1::jboolean)
    return jcall(FramesFactory, "getTIRF", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_tod(::Type{FramesFactory}, arg0::IERSConventions, arg1::jboolean)
    return jcall(FramesFactory, "getTOD", FactoryManagedFrame, (IERSConventions, jboolean), arg0, arg1)
end

function get_tod(::Type{FramesFactory}, arg0::jboolean)
    return jcall(FramesFactory, "getTOD", FactoryManagedFrame, (jboolean,), arg0)
end

function get_veis1950(::Type{FramesFactory})
    return jcall(FramesFactory, "getVeis1950", FactoryManagedFrame, ())
end

function set_eop_continuity_threshold(::Type{FramesFactory}, arg0::jdouble)
    return jcall(FramesFactory, "setEOPContinuityThreshold", void, (jdouble,), arg0)
end


text:
module AccessMacro

using MacroTools
export @access

get_fieldname(expr::Expr)   = expr.args[1]
get_fieldname(symb::Symbol) = symb

function field_to_function_expr(f, T)
    s = gensym()
    :(function $f($s :: $T)
        $s.$f
    end)
    #:( ($f)($s :: $T) = $s.$f )
end

macro access(expr)
    return_expr = Expr(:block)
    @capture(expr, struct (T_{P__}|T_) fields__ end | type (T_{P__}|T_) fields__ end | mutable struct (T_{P__}|T_) fields__ end | immutable (T_{P__}|T_) fields__ end)
    push!(return_expr.args, expr)
    push!(return_expr.args, field_to_function_expr.(get_fieldname.(fields), T)...)
    push!(return_expr.args, T)
    esc(return_expr)
end

end # module


text:
const TIMER_OUTPUTS = true
const being_timed = Ref{Bool}(false)

if TIMER_OUTPUTS
    using TimerOutputs

    macro timer(name, expr)
        :(if being_timed[]
              @timeit $(esc(name)) $(esc(expr))
          else
              $(esc(expr))
          end)
    end

    macro iftimer(expr)
        esc(expr)
    end

else
    macro timer(name, expr)
        esc(expr)
    end

    macro iftimer(expr)
    end
end

export @timer
export @iftimer


text:
@testset "geometricflux.jl" begin

    @testset "GraphConv" begin
        gnn = SeaPearl.GraphConv(3 => 3)

        graphs = Matrix.(adjacency_matrix.([random_regular_graph(10, 4) for i = 1:4]))
        nodeFeatures = [rand(Float32, 3, 10) for i = 1:4]
        featuredGraphs = [SeaPearl.FeaturedGraph(g; nf=nf) for (g, nf) in zip(graphs, nodeFeatures)]
        batchedFG = SeaPearl.BatchedFeaturedGraph{Float32}(featuredGraphs)

        classicOutput = (cat(map(fg->fg.nf,gnn.(featuredGraphs))...; dims=3))
        batchedOutput = gnn(batchedFG).nf

        @test all(abs.(classicOutput .- batchedOutput) .< 1e-5)

        grad = gradient(batchedFG) do inp
            gnn(inp).nf[1, 3, 2]
        end
        grad = grad[1]

        @test all(grad.nf[:, :, [1, 3, 4]] .== 0)
        @test !all(grad.nf[:, :, 2] .== 0)
        @test all(grad.graph[:, :, [1, 3, 4]] .== 0)
        @test !all(grad.graph[:, :, 2] .== 0)
    end
end


text:
# Tests for distances
function t_mds()

N = 10	# variables
M = 25  # samples

vectordata = rand(M)
matrixdata = rand(N,M)
labels = round.(2*rand(M))
multilabels = round.(2*rand(2,M))

A01 = j4pr.datacell(vectordata)
A02 = j4pr.datacell(vectordata, labels)
A03 = j4pr.datacell(vectordata, multilabels)
A04 = j4pr.datacell(matrixdata)
A05 = j4pr.datacell(matrixdata, labels)
A06 = j4pr.datacell(matrixdata, multilabels)

for tr in [A01, A02, A03, vectordata]
	Test.@test try
		tr |> j4pr.mds(1);
		true
	catch
		false
	end
end

for tr in [A04, A05, A06, matrixdata]
	Test.@test try
		tr |> j4pr.ica(2);
		true
	catch
		false
	end

end

end


text:
"""
```
visualize_flow(flow, ColorBased(), RasterConvention())
visualize_flow(flow, ColorBased(), CartesianConvention())
```

Returns an image that matches the dimensions of the input matrix and that depicts the
orientation and magnitude of the optical flow vectors using the HSV color space.

# Details
Hue encodes angle between optical flow vector and the x-axis in the image plane;
saturation encodes the ratio between the individual vector magnitudes and the maximum magnitude
among the whole motion field, and the values always equal one.

# Arguments
The flow parameter needs to be a two-dimensional arrays of length-2 vectors (of type SVector)
which represent the displacement of each pixel. The displacement can be expressed in a
coordinate system based on a RasterConvention() (i.e. rows and columns of a matrix)
or a CartesianConvention().

# Example

Compute the HSV encoded visualization of flow vectors in (row, column) convention.
```julia

using ImageTracking

hsv = visualize_flow(flow, ColorBased(), RasterConvention())

imshow(RGB.(hsv))
```

Compute the HSV encoded visualization of flow vectors in (x, y) convention.
```julia

using ImageTracking

hsv = visualize_flow(flow, ColorBased(), CartesianConvention())

imshow(RGB.(hsv))
```

# References
[1] S. Baker, D. Scharstein, JP Lewis, S. Roth, M.J. Black, and R. Szeliski. A database and
evaluation methodology for optical flow. International Journal of Computer Vision, 92(1):1–31, 2011.
"""
function visualize_flow(flow::Array{SVector{2, Float64}, 2}, method::ColorBased, convention::RasterConvention)
	# Convert from (row,column) to (x,y) convention.
	map!(x-> SVector(last(x),first(x)), flow, flow)

    # Display optical flow as an image, with hue encoding the orientation and
    # saturation encoding the relative magnitude.
    max_norm = maximum(map(norm, flow))
    normalised_flow = map(PolarFromCartesian(), flow / max_norm)
    hsv = zeros(HSV{Float32},size(flow))

    for i in eachindex(flow)
        hsv[i] = HSV((normalised_flow[i].θ + pi) * 180 / pi, normalised_flow[i].r, 1)
    end

    return hsv
end

function visualize_flow(flow::Array{SVector{2, Float64}, 2}, method::ColorBased, convention::CartesianConvention)
    # Display optical flow as an image, with hue encoding the orientation and
    # saturation encoding the relative magnitude.
    max_norm = maximum(map(norm, flow))
    normalised_flow = map(PolarFromCartesian(), flow / max_norm)
    hsv = zeros(HSV{Float32},size(flow))

    for i in eachindex(flow)
        hsv[i] = HSV((normalised_flow[i].θ + pi) * 180 / pi, normalised_flow[i].r, 1)
    end

    return hsv
end

function is_flow_known(p::SVector{2, Float64})
    return !isnan(first(p)) && !isnan(last(p)) && abs(first(p)) < 1e9 && abs(last(p)) < 1e9
end

"""
```
read_flow_file(file_path)
```

Returns a 2-Dimensional array that depicts the optical flow.

# Details
Reads simple .flo file format which is used for optical flow evaluation.
The .flo file should store 2-band float image for horizontal (u) and vertical (v) flow components.
Floats are stored in little endian format.

bytes contain -
  0-3     tag: "PIEH" in ASCII, which in little endian happens to be the float 202021.25
          (just a sanity check that floats are represented correctly)
  4-7     width as an integer
  8-11    height as an integer
  12-end  data (width*height*2*4 bytes total)
          the float values for u and v, interleaved, in row order, i.e.,
          u[row0,col0], v[row0,col0], u[row0,col1], v[row0,col1], ...

# Arguements
The file_path parameter needs to be a string, storing the address of .flo file.

# Example

Read flow file from .flo file and visualize the flow in HSV
```julia

using ImageTracking

file_path = "./Example_flow.flo"
flow = read_flow_file(file_path)
result = visualize_flow(flow, ColorBased(), RasterConvention())
imshow(RGB.(result))
```

# References
[1] http://vision.middlebury.edu/flow/data/
"""
# first four bytes to be same in little endian
# check this when reading a .flo file
FLOAT_TAG = 202021.25
function read_flow_file(file_path::String)

    io = open(file_path, "r")
    seekstart(io)

    tag = read(io, Float32) # reading the initial tag

    if tag != FLOAT_TAG
        print("read_from_file: Problem reading file ", filepath)
        return -1
    end

    width  = read(io, Int32)
    height = read(io, Int32)

    flow = Array{SVector{2, Float64}, 2}(undef,height,width)
    for i in 1:height
        for j in 1:width
            if eof(io) == false
                fx = read(io, Float32)
                fy = read(io, Float32)
                flow[i, j] = SVector{2, Float64}(fx,fy)
            end
        end
    end

    close(io)
    return flow
end

"""
```
evaluate_flow_error(ground_truth_flow, estimated_flow, EndpointError())
```

Returns a 2-Dimensional array that matches the dimensions of the input flow vector and
that depicts the end point error between the estimated flow and the ground truth flow.

# Details
If the estimated flow at a point is (u0, v0) and ground truth flow is (u1, v1), then
error will be sqrt[(u0 - u1)^2 + (v0 - v1)^2] at that point.

# Arguments
The flow parameters needs to be two-dimensional arrays of length-2 vectors (of type SVector)
which represent the displacement of each pixel.

# Example

Compute the end point error between two flows.
```julia

using ImageTracking

result = evaluate_flow_error(ground_truth_flow, estimated_flow, EndpointError())

imshow(result)
```

# References
[1] S. Baker, D. Scharstein, JP Lewis, S. Roth, M.J. Black, and R. Szeliski. A database and
evaluation methodology for optical flow. International Journal of Computer Vision, 92(1):1–31, 2011.
"""
function evaluate_flow_error(ground_truth_flow::Array{SVector{2, Float64}, 2}, estimated_flow::Array{SVector{2, Float64}, 2}, error_type::EndpointError)
    result = similar(Array{Float64, 2}, axes(estimated_flow))
    for i in CartesianIndices(estimated_flow)
        p1 = estimated_flow[i]
        p2 = ground_truth_flow[i]

        if is_flow_known(p1) && is_flow_known(p2)
            δ = p1 - p2
            result[i] = sqrt(sum(δ.^2))
        else
            result[i] = NaN
        end
    end

    return result
end

"""
```
evaluate_flow_error(ground_truth_flow, estimated_flow, AngularError())
```

Returns a 2-Dimensional array that matches the dimensions of the input flow vector and
that depicts the angle error between the estimated flow and the ground truth flow.

# Details
If the estimated flow at a point is (u0, v0) and ground truth flow is (u1, v1), it
calculates the angle between (u0, v0, 1) and (u1, v1, 1) vectors as measure for error.

# Arguments
The flow parameters needs to be two-dimensional arrays of length-2 vectors (of type SVector)
which represent the displacement of each pixel.

# Example

Compute the angle error between two flows.
```julia

using ImageTracking

result = evaluate_flow_error(ground_truth_flow, estimated_flow, AngularError())

imshow(result)
```

# References
[1] S. Baker, D. Scharstein, JP Lewis, S. Roth, M.J. Black, and R. Szeliski. A database and
evaluation methodology for optical flow. International Journal of Computer Vision, 92(1):1–31, 2011.
"""
function evaluate_flow_error(ground_truth_flow::Array{SVector{2, Float64}, 2}, estimated_flow::Array{SVector{2, Float64}, 2}, error_type::AngularError)
    result = similar(Array{Float64, 2}, axes(estimated_flow))
    for i in CartesianIndices(estimated_flow)
        p1 = ground_truth_flow[i]
        p2 = estimated_flow[i]
        if is_flow_known(p2) && is_flow_known(p1)
            p1 = push(p1, 1.0)
            p2 = push(p2, 1.0)
            cosine = dot(p1, p2)/(norm(p1)*norm(p2))
            if cosine > 1
                result[i] = 0
            else
                result[i] = acos(cosine)
            end
            if result[i] <= 1e-5
                result[i] = 0
            end
        else
            result[i] = NaN
        end
    end

    return result
end

"""
```
calculate_statistics(error)
```

Quantifies the flow error in terms of a mean and standard deviation as well as two additional
statistics that measure "robustness".

# Details
RX and AX are the robustness statistics. RX denotes the percentage of pixels that have
an error measure over X. AX denotes the accuracy of the error measure at the Xth percentile.

# Arguements
The error parameters needs to be two-dimensional arrays of length-2 vectors (of type SVector)
which represent the error between two flows.

# Example

Compute the Mean, SD, RX, AX stats of the flow error calculated using endpoint error method.
```julia

using ImageTracking

error = evaluate_flow_error(ground_truth_flow, estimated_flow, EndpointError())
mean, sd, rx, ax = calculate_statistics(error)
```

Compute the mean, SD, RX, AX stats of the flow error calculated using angula error method.
```julia

using ImageTracking

error = evaluate_flow_error(ground_truth_flow, estimated_flow, AngularError())
mean, sd, rx, ax = calculate_statistics(error)
```

# References
[1] S. Baker, D. Scharstein, JP Lewis, S. Roth, M.J. Black, and R. Szeliski. A database and
evaluation methodology for optical flow. International Journal of Computer Vision, 92(1):1–31, 2011.
"""
function calculate_statistics(error::Array{Float64, 2})
    R_thresholds = [0.5, 1.0, 2.0, 3.0, 5.0]
    A_thresholds = [0.5, 0.75, 0.95]

    # mean and standard deviations
    mean_value = mean(skipmissing(error))
    std_value = std(skipmissing(error); mean=mean_value)

    # RX stats
    RX_count = Dict()
    for rx in R_thresholds
        count_above_thresh = count(x->(x > rx), error)
        RX_count[rx] = count_above_thresh / (size(error)[1]*size(error)[2])
    end

    # AX stats
    AX_count = Dict()
    edges, counts = imhist(error)
    max_val = maximum(counts)
    total_pix = length(error)
    pixel_cumulative_count = cumsum(counts)
    for ax in A_thresholds
        cutoff = floor(ax*3.0 + 0.5)
        count_below_thresh = count(x->(x < cutoff), pixel_cumulative_count)
        AX_count[ax] = count_below_thresh / (length(edges)*max_val)
    end

    return mean_value, std_value, RX_count, AX_count
end


text:
using ElectrochemicalKinetics
using Plots
using NLsolve
using QuadGK

# define some constants and default parameters
kB = 8.617e-5
T = 298
muoA = 0.02
muoB = 0.03
Ω = 0.1

# our familiar thermodynamic functions
hs(x;Ω=Ω) = x*(1-x)*Ω # enthalpy of mixing
s(x) = -kB*(x*log(x) + (1-x)*log(1-x)) # entropy per particle
g_thermo(x; Ω=Ω, muoA=muoA, muoB=muoB, T=T) = hs(x;Ω=Ω) - T*s(x)+ muoA*(1-x) + muoB*x # Gibbs free energy per particle
μ_thermo(x; Ω=Ω, muoA=muoA, muoB=muoB, T=T) = (1-2*x)*Ω + kB*T*log(x/(1-x)) + muoB-muoA # chemical potential

"""
µ and g with kinetic constributions, can be modeled using any <:KineticModel object

These functions return single-argument functions (to easily use common-tangent function below while
still being able to swap out model parameters by calling "function-builders" with different arguments).
"""
µ_kinetic(I, km::KineticModel; Ω=Ω, muoA=muoA, muoB=muoB, T=T) =
    x -> µ_thermo(x; Ω=Ω, muoA=muoA, muoB=muoB, T=T) + fit_overpotential((1-x)*km, I)[1]
g_kinetic(I, km::KineticModel; Ω=Ω, muoA=muoA, muoB=muoB, T=T) =
    x -> g_thermo(x; Ω=Ω, muoA=muoA, muoB=muoB, T=T) + quadgk(x->fit_overpotential((1-x)*km, I), 0, x)[1][1]

# zeros of this function correspond to pairs of x's satisfying the common tangent condition for a given µ function
function common_tangent(x, I, km::KineticModel; Ω=Ω, muoA=muoA, muoB=muoB, T=T)
    g = g_kinetic(I, km; Ω=Ω, muoA=muoA, muoB=muoB, T=T)
    µ = µ_kinetic(I, km; Ω=Ω, muoA=muoA, muoB=muoB, T=T)
    [(g(x[2])-g(x[1]))/(x[2]-x[1]) - μ(x[1]), μ(x[2])-μ(x[1])]
end

function find_phase_boundaries(I, km::KineticModel; Ω=Ω, muoA=muoA, muoB=muoB, T=T)
    function myct!(storage, x)
        res = common_tangent(x, I, km; Ω=Ω, muoA=muoA, muoB=muoB, T=T)
        storage[1] = res[1]
        storage[2] = res[2]
    end
    # TODO: grad! here from AD
    x1 = nlsolve(myct!, [0.05 0.95])
    x1.zero
end

bv = ButlerVolmer(300)
I_vals = 10 .^ (1.1:0.025:3.1)

# # this line takes a few seconds with B-V but aaages with anything else...
# pbs = [find_phase_boundaries(I, bv, T=330) for I in I_vals]

# plot(vcat(pbs...), hcat(I_vals, I_vals), label="phase boundary")
# xlabel!("x")
# ylabel!("I")


text:
using Documenter, SaferIntegers

makedocs(
    modules = [SaferIntegers],
    sitename = "SaferIntegers.jl",
    authors = "Jeffrey Sarnoff",
    pages = Any[
        "Overview" => "index.md",
        "Highlights" => "highlights.md",
        "Basic Guide" => "basicguide.md",
        "Supported Operations" => "supports.md",
        "How To Use" => "howtouse.md",
        "Safe Shifts" => "safeshifts.md",
        "Safer Rationals" => "rationals.md",
        "Benchmarks" => "benchmarks.md",
        "Refs" => "references.md"
    ]
)

deploydocs(
    repo = "github.com/JeffreySarnoff/SaferIntegers.jl.git",
    target = "build"
)


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/evaluator/passes.jl
# Functions used to compute forward pass of nonlinear functions which include:
# set_value_post, overwrite_or_intersect, forward_pass_kernel, associated blocks
#############################################################################

xnum_yset(b, x, y) = is_num(b, x) && !is_num(b, y)
xset_ynum(b, x, y) = !is_num(b, x) && is_num(b, y)
xy_num(b, x, y) = is_num(b, x) && is_num(b, y)
xyset(b, x, y) = !(is_num(b, x) || is_num(b, y))

function varset(::Type{MC{N,T}}, i, x_cv, x_cc, l, u) where {V,N,T<:RelaxTag}
    v = seed_gradient(i, Val(N))
    return MC{N,T}(x_cv, x_cc, Interval{Float64}(l, u), v, v, false)
end

function fprop!(t::Relax, vt::Variable, g::DAT, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
    i = first_index(g, k)
    x = val(b, i)
    l = lbd(b, i)
    u = ubd(b, i)
    if l == u
        b[k] = x
    else
        z = varset(MC{N,T}, rev_sparsity(g, i, k), x, x, l, u)
        if !first_eval(t, b)
            z = z ∩ interval(b, k)
        end
        b[k] = z
    end
    nothing
end

function fprop!(t::Relax, ex::Subexpression, g::DAT, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
    x =  first_index(g, k)
    if subexpression_is_num(b, x)
        b[k] = subexpression_num(b, x)
    else
        b[k] = subexpression_set(b, x)
    end
end

for (F, f) in ((PLUS, :+), (MIN, :min), (MAX, :max), (DIV, :/), (ARH, :arh))
    @eval function fprop_2!(t::Relax, v::Val{$F}, g::DAT, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
        x = child(g, 1, k)
        y = child(g, 2, k)
        if !xy_num(b, x, y)
            if xyset(b, x, y)
                z = ($f)(set(b, x), set(b, y))
            elseif xset_ynum(b, x, y)
                z = ($f)(set(b, x), num(b, y))
            else
                z = ($f)(num(b, x), set(b, y))
            end
            b[k] = cut(z, set(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
        else
            b[k] = ($f)(num(b, x), num(b, y))
        end
    end
end

function fprop!(t::Relax, v::Val{MINUS}, g::DAT, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
    x = child(g, 1, k)
    if is_binary(g, k)
        y = child(g, 2, k)
        if !xy_num(b, x, y)
            if xyset(b, x, y)
                z = set(b, x) - set(b, y)
            elseif xset_ynum(b, x, y)
                z = set(b, x) - num(b, y)
            else
                z = num(b, x) - set(b, y)
            end
            b[k] = cut(z, set(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
        else
            b[k] = num(b, x) - num(b, y)
        end
    else
        if is_num(b, x)
            b[k] = -num(b, x)
        else
            z = -set(b, x)
            b[k] = cut(-set(b, x), set(b, k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
        end
    end
end

function fprop_2!(t::Relax, v::Val{MULT}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}

    x = child(g, 1, k)
    y = child(g, 2, k)

    if !xy_num(b, x, y)
        if xyset(b, x, y)
            #@show "2 arity multiplication..."
            xv = set(b, x)
            yv = set(b, y)
            if b.use_apriori_mul
                dp = b.dp
                dP = b.dP
                s = sparsity(g, 1)

                xr = info(b, x)
                yr = info(b, y)

                xr_cv = cv(xr)                   # GOOD
                xr_cc = cc(xr)                   # GOOD
                yr_cv = cv(yr)                   # GOOD
                yr_cc = cc(yr)                   # GOOD
                xr_cv_grad = cv_grad(xr)         # GOOD
                xr_cc_grad = cc_grad(xr)         # GOOD
                yr_cv_grad = cv_grad(yr)         # GOOD
                yr_cc_grad = cc_grad(yr)         # GOOD

                xrn_cc = -xr_cc
                yrn_cc = -yr_cc
                xrn_cc_grad = -xr_cc_grad
                yrn_cc_grad = -yr_cc_grad

                t3 = affine_expand_del(dP, xr_cv, xr_cv_grad, s)
                t4 = affine_expand_del(dP, yr_cv, yr_cv_grad, s)
                s3 = affine_expand_del(dP, xrn_cc, xrn_cc_grad, s)
                s4 = affine_expand_del(dP, yrn_cc, yrn_cc_grad, s)

                z = xv*yv
                #@show z
                wIntv = z.Intv
                if (t3 < xv.Intv.hi) || (t4 < yv.Intv.hi)
                    t1 = affine_expand_del(dp, xr_cv, xr_cv_grad, s)
                    t2 = affine_expand_del(dp, yr_cv, yr_cv_grad, s)
                    za_l = McCormick.mult_apriori_kernel(xv, yv, wIntv, t1, t2, t3, t4, xr_cv_grad, yr_cv_grad)
                    #@show za_l
                    z = z ∩ za_l
                end
                #=
                if (s3 > -xv.Intv.lo) || (s4 > -yv.Intv.lo)
                    s1 = affine_expand_del(dp, xrn_cc, xrn_cc_grad, s)
                    s2 = affine_expand_del(dp, yrn_cc, yrn_cc_grad, s)
                    za_u = McCormick.mult_apriori_kernel(-xv, -yv, wIntv, s1, s2, s3, s4, xrn_cc_grad, yrn_cc_grad)
                    @show za_u
                    z = z ∩ za_u
                end
                =#
                #@show z
            else
                z = xv*yv
            end
        elseif xset_ynum(b, x, y)
            z = set(b, x)*num(b, y)
        else
            z = num(b, x)*set(b, y)
        end
        b[k] = cut(z, set(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
    else
        b[k] = num(b, x)*num(b, y)
    end
end

function fprop_n!(t::Relax, v::Val{PLUS}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    z = zero(MC{N,T})
    znum = 0.0
    numval = true
    for i in children(g, k)
        if is_num(b, i)
            znum += num(b, i)
        else
            numval = false
            z += set(b, i)
        end
    end
    if numval
        b[k] = znum
    else
        z += znum
        b[k] = cut(z, set(b, k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
    end
end

function fprop_n!(t::Relax, v::Val{MIN}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    z = Inf*one(MC{N,T})
    znum = Inf
    numval = true
    for i in children(g, k)
        if is_num(b, i)
            znum = min(znum, num(b, i))
        else
            numval = false
            z = min(z, set(b, i))
        end
    end
    if numval
        b[k] = znum
    else
        z = min(z, znum)
        b[k] = cut(z, set(b, k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
    end
end

function fprop_n!(t::Relax, v::Val{MAX}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    z = -Inf*one(MC{N,T})
    znum = -Inf
    numval = true
    for i in children(g, k)
        if is_num(b, i)
            znum = max(znum, num(b, i))
        else
            numval = false
            z = max(z, set(b, i))
        end
    end
    if numval
        b[k] = znum
    else
        z = max(z, znum)
        b[k] = cut(z, set(b, k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
    end
end

function fprop_n!(t::Relax, ::Val{MULT}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    #println("start MULT N")
    #println(" ")

    dp = b.dp
    dP = b.dP
    s = sparsity(g, 1)

    #@show s
    #@show dp
    #@show dP

    z = one(MC{N,T})
    zt = one(MC{N,T})
    zr = one(MC{N,T})
    znum = one(Float64)
    firstset = true
    numval = true
    for (q,i) in enumerate(children(g, k))
        #println(" ")
        #@show q
        if is_num(b, i)
            znum = znum*num(b, i)
        else
            numval = false
            x = set(b, i)
            xr = info(b, i)
            #@show x
            #@show xr
            if b.use_apriori_mul
                if !firstset
                    xr_cv = cv(xr)
                    zr_cv = cv(zr)
                    xr_cv_grad = cv_grad(xr)
                    zr_cv_grad = cv_grad(zr)

                    xrn_cc = -cc(xr)
                    zrn_cc = -cc(zr)
                    xrn_cc_grad = -cc_grad(xr)
                    zrn_cc_grad = -cc_grad(zr)

                    t3 = affine_expand_del(dP, xr_cv, xr_cv_grad, s)
                    t4 = affine_expand_del(dP, zr_cv, zr_cv_grad, s)
                    s3 = affine_expand_del(dP, xrn_cc, xrn_cc_grad, s)
                    s4 = affine_expand_del(dP, zrn_cc, zrn_cc_grad, s)

                    wIntv = x.Intv*z.Intv
                    zm = x*z
                    #@show zm
                    if (t3 < x.Intv.hi) || (t4 < z.Intv.hi)
                        # Relaxations are incorrect for za_l
                        t1 = affine_expand_del(dp, xr_cv, xr_cv_grad, s)
                        t2 = affine_expand_del(dp, zr_cv, zr_cv_grad, s)

                        #=
                        println(" ")
                        @show x
                        @show xr
                        @show z
                        @show zr

                        println(" ")
                        @show wIntv
                        @show xr_cv_grad
                        @show zr_cv_grad

                        println(" ")
                        @show t1
                        @show t2
                        @show t3
                        @show t4

                        println(" ")
                        @show s
                        @show dp
                        @show dP
                        =#

                        za_l = McCormick.mult_apriori_kernel(x, z, wIntv, t1, t2, t3, t4, xr_cv_grad, zr_cv_grad)
                        #@show za_l
                        #@show zm
                        #@show za_l
                        zm = zm ∩ za_l
                    end

                    if (s3 > -x.Intv.lo) || (s4 > -z.Intv.lo)
                        #@show "improved overestimators"
                        s1 = affine_expand_del(dp, xrn_cc, xrn_cc_grad, s)
                        s2 = affine_expand_del(dp, zrn_cc, zrn_cc_grad, s)
                        #=
                        println(" ")
                        @show x
                        @show xr
                        @show z
                        @show zr

                        println(" ")
                        @show -x
                        @show -z

                        println(" ")
                        @show wIntv
                        @show xrn_cc
                        @show zrn_cc
                        @show xrn_cc_grad
                        @show zrn_cc_grad

                        println(" ")
                        @show s1
                        @show s2
                        @show s3
                        @show s4

                        println(" ")
                        @show s
                        @show dp
                        @show dP
                        =#

                        za_u = McCormick.mult_apriori_kernel(-x, -z, wIntv, s1, s2, s3, s4, xrn_cc_grad, zrn_cc_grad)
                        #@show za_u
                        #@show zm
                        zm = zm ∩ za_u
                    end

                    z = zm
                    zr = xr*zr
                else
                    firstset = false
                    zr = xr
                    z = x
                end
            else
                z = z*x
            end
            zt = zt*x
            #@show z
            #@show zt
        end
    end
    if numval
        b[k] = znum
    else
        z = z*znum
        b[k] = cut(z, set(b, k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
    end
end

for F in (PLUS, MULT, MIN, MAX)
    @eval function fprop!(t::Relax, v::Val{$F}, g::DAT, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
        is_binary(g, k) ? fprop_2!(Relax(), Val($F), g, b, k) : fprop_n!(Relax(), Val($F), g, b, k)
    end
end
function fprop!(t::Relax, v::Val{DIV}, g::DAT, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
    fprop_2!(Relax(), v, g, b, k)
end

function fprop!(t::Relax, v::Val{POW}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    x = child(g, 1, k)
    y = child(g, 2, k)
    if is_num(b, y) && isone(num(b, y))
        b[k] = set(b, x)
    elseif is_num(b,y) && iszero(num(b, y))
        b[k] = one(V)
    elseif !xy_num(b, x, y)
        if xyset(b, x, y)
            z = set(b, x)^set(b, y)
        elseif xset_ynum(b, x, y)
            z = set(b, x)^num(b, y)
        else
            z = num(b, x)^set(b, y)
        end
        b[k] = cut(z, set(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
    else
        b[k] = num(b, x)^num(b, y)
    end
end

function fprop!(t::Relax, v::Val{USER}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    f = user_univariate_operator(g, first_index(g, k))
    x = child(g, 1, k)
    if is_num(b, x)
        b[k] = f(num(b, x))
    else
        z = f(set(b, x))
        b[k] = cut(z, set(b, k), b.ic.v, zero(Float64), sparsity(g, k), b.cut, b.post)
    end
end

function fprop!(t::Relax, v::Val{USERN}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    mv = user_multivariate_operator(g, first_index(g, k))
    n = arity(g, k)
    set_input = _set_input(b, n)
    num_input = _num_input(b, n)
    anysets = false
    i = 1
    for c in children(g, k)
        if is_num(b, c)
            x = num(b, c)
            if !isinf(x)
                set_input[i] = MC{N,T}(x)
                num_input[i] = x
            end
        else
            set_input[i] = set(b, c)
            anysets = true
        end
        i += 1
    end
    if anysets
        z = MOI.eval_objective(mv, set_input)::MC{N,T}
        b[k] = cut(z, set(b, k), b.ic.v, zero(Float64), sparsity(g,k), b.cut, b.post)
    else
        b[k] = MOI.eval_objective(mv, num_input)
    end
end

for ft in UNIVARIATE_ATOM_TYPES
    f = UNIVARIATE_ATOM_DICT[ft]
    (f == :user || f == :+ || f == :-) && continue
    @eval function fprop!(t::Relax, v::Val{$ft}, g::DAT, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
        x = child(g, 1, k)
        if is_num(b, x)
            return b[k] = ($f)(num(b, x))
        else
            z = ($f)(set(b, x))
            b[k] = cut(z, set(b, k), b.ic.v, zero(Float64), sparsity(g,k), b.cut, b.post)
        end
    end
end

for (F, f) in ((LOWER_BND, :lower_bnd), (UPPER_BND, :upper_bnd))
    @eval function fprop!(t::Relax, v::Val{$F}, g::DAT, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
        y = child(g, 2, k)
        if is_num(b, y)
            z = set(b, child(g, 1, k))
            z = ($f)(z, num(b, y))
            b[k] = cut(z, set(b, k), b.ic.v, zero(Float64), sparsity(g, k), b.cut, b.post)
        end
    end
end

function fprop!(t::Relax, v::Val{BND}, g::DAT, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag}
    z = set(b, child(g, 1, k))
    y = child(g, 2, k)
    r = child(g, 3, k)
    if is_num(b, y) && is_num(b, r)
        z = bnd(z, num(b, y),num(b, r))
    end
    b[k] = cut(z, set(b, k), b.ic.v, zero(Float64), sparsity(g,k), b.cut, b.post)
end

function f_init!(t::Relax, g::DAT, b::RelaxCache)
    for k = node_count(g):-1:1
        c = node_class(g, k)
        (c == EXPRESSION)    && fprop!(t, Expression(), g, b, k)
        (c == VARIABLE)      && fprop!(t, Variable(), g, b, k)
        (c == SUBEXPRESSION) && fprop!(t, Subexpression(), g, b, k)
        b._info[k] = set(b, k)
    end
    for k = node_count(g):-1:1
        if is_binary(g, k)
            x = child(g, 1, k)
            y = child(g, 2, k)
        end
    end
    nothing
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions