text:
################################################################################
#                          EVALUATION METHOD DATATYPES
################################################################################
"""
    AbstractIntegralMethod

An abstract type for integral evaluation methods use in combination with
`integral` and `generate_integral_data`.
"""
abstract type AbstractIntegralMethod end

"""
    Automatic <: AbstractIntegralMethod

An integral evaluation type for automically selecting an appropriate integral
evaluation method. Contains no fields.
"""
struct Automatic <: AbstractIntegralMethod end

"""
    AbstractUnivariateMethod <: AbstractIntegralMethod

An abstract type for integral evaluation methods for 1-dimensional integrals.
"""
abstract type AbstractUnivariateMethod <: AbstractIntegralMethod end

"""
    FiniteGuassQuad <: AbstractUnivariateMethod

Abstract type for Finite 1-D integral evaluation using a variety of Gaussian Quadrature methods over the entire domain.
"""
abstract type FiniteGaussQuad <: AbstractUnivariateMethod end

"""
    UniTrapezoid <: AbstractUnivariateMethod

An integral evalution method that uses the trapezoid rule to in combination
with all parameter supports available when the integral is expanded and/or when
the infinite model is optimized, whichever comes first. Note this method will
ignore the `num_supports` keyword argument. The upper and lower bounds of the integral
will automatically be added as supports. Note this is valid only for finite
integral domains. Contains no fields.
"""
struct UniTrapezoid <: AbstractUnivariateMethod end

"""
    UniMCSampling <: AbstractUnivariateMethod

An integral evaluation method that uses uniform Monte Carlo sampling to
approximate the integral. This variant will add more supports to the model as
needed to satisfy `num_supports` and it will include all supports with the
`MCSample` label up till the integral is expanded and/or when
the infinite model is optimized, whichever comes first. Note this is valid only
for finite integral domains. Contains no fields.
"""
struct UniMCSampling <: AbstractUnivariateMethod end

"""
    UniIndepMCSampling <: AbstractUnivariateMethod

An integral evaluation method that uses uniform Monte Carlo sampling to
approximate the integral similar to [`UniMCSampling`](@ref MeasureToolbox.UniMCSampling).
However, this variant will generate its own set of supports and ignore all other
supports with the `MCSample` label. Note this is valid only for finite integral
domains. This is not compatible with individual dependent parameters. Contains no fields.
"""
struct UniIndepMCSampling <: AbstractUnivariateMethod end

"""
    Quadrature <: AbstractUnivariateMethod

A general integral evaluation method that will automatically select the
appropriate quadrature method to approximate the integral. Please note that this
will generate a unique set of parameter supports and will ignore existing supports
when the integral is evaluated and thus should be used with caution. However,
this method is able to handle infinite and semi-infinite integral domains.
This is not compatible with individual dependent parameters. Contains no fields.
"""
struct Quadrature <: AbstractUnivariateMethod end

"""
    FEGaussLobatto <: AbstractUnivariateMethod

Integral evaluation method that allows for the user to specify supports to be included in
quadrature evaluation. The upper and lower bounds of the integral will automatically
be added as supports. This method uses Gauss Lobatto quadrature to decompose the overall Integral
into smaller integrals that span the user defined supports as follows:

``\\int_{x_1}^{x_3} f(x) dx = \\int_{x_1}^{x_2} f(x) dx + \\int_{x_2}^{x_3} f(x) dx``

where the integrals are evaluated using Gauss Lobatto quadrature:

``\\int f(x) dx \\approx \\sum_{i=1}^{n} \\alpha_i f(x_i)``

"""
struct FEGaussLobatto <: AbstractUnivariateMethod end

"""
    GaussLegendre <: FiniteGaussQuad

An integral evaulation method that uses Gauss-Legendre quadrature to
evaluate integrals. This is valid for finite integral domains. Note this will
generate its own set of supports and will ignore other parameter supports.
This is not compatible with individual dependent parameters. Contains no fields.
"""
struct GaussLegendre <: FiniteGaussQuad end

"""
    GaussRadau <: FiniteGaussQuad

An integral evaulation method that uses Gauss-Radau quadrature to
evaluate integrals. This is valid for finite integral domains. Note this will
generate its own set of supports and will ignore other parameter supports.
This is not compatible with individual dependent parameters. Contains no fields.
"""
struct GaussRadau <: FiniteGaussQuad end

"""
    GaussLobatto <: FiniteGaussQuad

An integral evaulation method that uses Gauss-Lobatto quadrature to
evaluate integrals. This is valid for finite integral domains. Note this will
generate its own set of supports and will ignore other parameter supports.
This is not compatible with individual dependent parameters. Contains no fields.
"""
struct GaussLobatto <: FiniteGaussQuad end

"""
    GaussJacobi <: FiniteGaussQuad

An integral evaulation method that uses Gauss-Jacobi quadrature to
evaluate integrals. It will take this form:

``\\int_{-1}^{1} f(x) (1-x)^\\alpha (1+x)^\\beta dx \\approx \\sum_{i=1}^{n} \\alpha_i f(x_i)``

Where,

`` (1-x)^\\alpha (1+x)^\\beta``

is the weight function.
This is valid for finite integral domains. This requires the user
to input the alpha and beta shape parameters for their function. This will then
generate its own set of supports and will ignore other parameter supports.
This is not compatible with individual dependent parameters.
If α or β is < -1, an error will be returned.

**Fields**
- `α::Float64`: Shape parameter that must be > -1
- `β::Float64`: Shape parameter that must be > -1
"""
struct GaussJacobi <: FiniteGaussQuad
    α::Float64
    β::Float64
    function GaussJacobi(α::Real, β::Real)
        if α < -1 || β < -1
            error("α and β must be greater than -1 for GaussJacobi")
        end
        return new(Float64(α), Float64(β))
    end
end

"""
    GaussChebyshev <: FiniteGaussQuad

An integral evaulation method that uses Gauss-Chebyshev quadrature to
evaluate integrals. This is valid for finite integral domains. This requires the user
to input the order of Guass-Chebyshev Quadrature they want to use.
If the order is not between 1 and 4 an error will be returned.
The integral evaluated is as follows:

``\\int_{-1}^{1} f(x) w(x) \\approx \\sum_{i=1}^{n} \\alpha_i f(x_i)``

The weight functions are as follows:

- 1st order: ``w(x)  =  \\frac{1}{\\sqrt{1-x^2}}``
- 2nd order: ``w(x) = {\\sqrt{1-x^2}}``
- 3rd order: ``w(x) = \\sqrt{(1+x)/(1-x)}``
- 4th order: ``w(x) = \\sqrt{(1-x)/(1+x)}``

This will then generate its own set of supports and will ignore other parameter supports.
This is not compatible with individual dependent parameters.

**Fields**
- `order::Int`: Specifies the order of Gauss-Chebyshev Quadrature. Must be between 1 and 4.
"""
struct GaussChebyshev <: FiniteGaussQuad
    order::Int
    function GaussChebyshev(order::Int)
        if order < 1 || order > 4
            error("Order must be between 1 and 4 for GaussChebyshev")
        end
        return new(order)
    end
end

"""
    GaussLaguerre <: AbstractUnivariateMethod

An integral evaulation method that uses Gauss-Laguerre quadrature to
evaluate integrals. This is valid for semi-infinite integral domains.

This method evaluates the following integral:

``\\int_{0}^{+∞} f(x) e^{-x} \\approx \\sum_{i=1}^{n} \\alpha_i f(x_i)``

Using the weight function:

``w(x) = e^{-x}``

Note this will generate its own set of supports and will ignore other parameter supports.
This is not compatible with individual dependent parameters.
"""
struct GaussLaguerre <: AbstractUnivariateMethod end

"""
    GaussHermite <: AbstractUnivariateMethod

An integral evaulation method that uses Gauss-Hermite quadrature to
evaluate integrals. This is valid for infinite integral domains.
It will take this form:

``\\int_{-∞}^{∞} f(x) e^{-x^2} \\approx \\sum_{i=1}^{n} \\alpha_i f(x_i)``

Using the weight function:
``w(x)`` = ``e^{-x^2}``

Note this will generate its own set of supports and will ignore other parameter supports.
This is not compatible with individual dependent parameters.
"""
struct GaussHermite <: AbstractUnivariateMethod end

"""
    AbstractMultivariateMethod <: AbstractIntegralMethod

An abstract type for integral evaluation methods for multi-dimensional integrals.
"""
abstract type AbstractMultivariateMethod <: AbstractIntegralMethod end

"""
    MultiMCSampling <: AbstractMultivariateMethod

An integral evaluation method that uses uniform Monte Carlo sampling to
approximate the integral. This variant will add more supports to the model as
needed to satisfy `num_supports` and it will include all supports with the
`MCSample` label up till the integral is expanded and/or when
the infinite model is optimized, whichever comes first. Note this is valid only
for finite integral domains. If an array of independent infinite parameters
is specified, they must use the same amount of supports. Contains no fields.
"""
struct MultiMCSampling <: AbstractMultivariateMethod end

"""
    MultiIndepMCSampling <: AbstractMultivariateMethod

An integral evaluation method that uses uniform Monte Carlo sampling to
approximate the integral similar to [`MultiMCSampling`](@ref MeasureToolbox.MultiMCSampling).
However, this variant will generate its own set of supports and ignore all other
supports with the `MCSample` label. Note this is valid only for finite integral
domains. Contains no fields.
"""
struct MultiIndepMCSampling <: AbstractMultivariateMethod end

################################################################################
#                             INTERNAL SUPPORT LABELS
################################################################################
"""
    InternalGaussLobatto <: InfiniteOpt.InternalLabel

A support label Gauss Lobatto points that are used as generative supports.
"""
struct InternalGaussLobatto <: InfiniteOpt.InternalLabel end

################################################################################
#                       DATA GENERATION METHODS (UNIVARIATE)
################################################################################
"""
    generate_integral_data(
        prefs::Union{InfiniteOpt.GeneralVariableRef, Vector{InfiniteOpt.GeneralVariableRef}},
        lower_bounds::Union{Real, Vector{<:Real}},
        upper_bounds::Union{Real, Vector{<:Real}},
        method::V; [num_supports::Int = InfiniteOpt.DefaultNumSupports,
        weight_func::Function = InfiniteOpt.default_weight,
        extra_kwargs...]
        )::InfiniteOpt.AbstractMeasureData where {V <: AbstractIntegralMethod}

Generate the appropriate concrete realization of `AbstractMeasureData` using
`method`. Here `prefs`, `lower_bounds`, and `upper_bounds` will always have a
1 to 1 correspondence when this is called from `integral`. Please refer
to the method docstrings for an explanation of each one.

User-defined method extensions should first
define a concrete `method` type inheriting from `AbstractUnivariateMethod` or
`AbstractMultivariateMethod` as appropriate and then implement extend this
method using that type for `method`.
"""
function generate_integral_data end

# General fallback
function generate_integral_data(prefs, lb, ub, method; kwargs...)
    error("`generate_integral_data` is not defined for method type `$(method)` " *
          "in accordance with the arguments given.")
end

# Single pref with Automatic
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
    lower_bound::Real,
    upper_bound::Real,
    method::Automatic;
    num_supports::Int = InfiniteOpt.DefaultNumSupports,
    weight_func::Function = InfiniteOpt.default_weight
    )::InfiniteOpt.AbstractMeasureData
    is_depend = InfiniteOpt._index_type(pref) == InfiniteOpt.DependentParameterIndex
    inf_bound_num = (lower_bound == -Inf) + (upper_bound == Inf)
    if inf_bound_num == 0 # finite interval
        return generate_integral_data(pref, lower_bound, upper_bound, UniTrapezoid(),
                                     num_supports = num_supports,
                                     weight_func = weight_func)
    elseif inf_bound_num == 1 && !is_depend # semi-infinite interval
        return generate_integral_data(pref, lower_bound, upper_bound,
                                     GaussLaguerre(), num_nodes = num_supports,
                                     weight_func = weight_func)
    elseif inf_bound_num == 2 && !is_depend # infinite interval
        return generate_integral_data(pref, lower_bound, upper_bound, GaussHermite(),
                                  num_nodes = num_supports,
                                  weight_func = weight_func)
    else
        error("Cannot generate measure data for individual dependent parameters ",
              "with infinite or semi-infinite domains.")
    end
end

# Univariate trapezoid coefficient function
function _trapezoid_coeff(supps::Vector{<:Real})::Vector{Float64}
    len = length(supps)
    len >= 2 || error("Cannot invoke `UniTrapezoid` on integral if there are ",
                      "less than 2 supports added to its infinite parameter. Ensure ",
                      "all infinite parameters have supports.")
    coeffs = Vector{Float64}(undef, len)
    coeffs[1] = (supps[2] - supps[1]) / 2
    for i in 2:len-1
        @inbounds coeffs[i] = (supps[i + 1] - supps[i - 1]) / 2
    end
    coeffs[end] = (supps[end] - supps[end - 1]) / 2
    return coeffs
end

# Single pref trapezoid rule
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
                                lower_bound::Real,
                                upper_bound::Real,
                                method::UniTrapezoid;
                                num_supports::Int = InfiniteOpt.DefaultNumSupports,
                                weight_func::Function = InfiniteOpt.default_weight
                                )
    return InfiniteOpt.FunctionalDiscreteMeasureData(pref, _trapezoid_coeff, 0,
                                                     InfiniteOpt.All,
                                                     InfiniteOpt.NoGenerativeSupports(),
                                                     weight_func, lower_bound,
                                                     upper_bound, false)
end

# Function for generating coefficients for FE Gauss Lobatto Quadrature
function _lobatto_coeff(supps::Vector{<:Real}, num_nodes::Int)::Vector{Float64}
    len = length(supps)
    len >= 2 || error("Cannot invoke `FEGaussLobatto` on integral if there are ",
                      "less than 2 supports added to its infinite parameter. Ensure ",
                      "all infinite parameters have supports.")
    sort!(supps)
    (_, coeffs) = FastGaussQuadrature.gausslobatto(num_nodes)
    coeff_vector = zeros((len - 1) * (num_nodes - 1) + 1)
    idx = 1:num_nodes
    for i in 1:len-1
        @inbounds coeff_vector[idx] += ((supps[i+1] - supps[i]) / 2 * coeffs)
        idx = idx .+ (num_nodes - 1)
    end
    return coeff_vector
end

# Single Pref Finite Gauss Lobatto Quadrature
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
    lower_bound::Real,
    upper_bound::Real,
    method::FEGaussLobatto;
    num_nodes::Int = 3,
    weight_func::Function = InfiniteOpt.default_weight
    )
    is_depend = InfiniteOpt._index_type(pref) == InfiniteOpt.DependentParameterIndex
    if lower_bound == -Inf || upper_bound == Inf
        error("The `FEGaussLobatto` can only be applied on finite intervals.")
    elseif num_nodes < 2
        error("For `FEGaussLobatto`, `num_nodes` must be at least 2.")
    elseif num_nodes > 2 && is_depend
        error("Cannot use more than 2 nodes with `FEGaussLobatto` with respect ",
              "to a dependent infinite parameter.")
    end
    (nodes, _) = FastGaussQuadrature.gausslobatto(num_nodes)
    info = UniformGenerativeInfo(nodes[2:end-1], InternalGaussLobatto, -1, 1)
    coeff_func = (supps) -> _lobatto_coeff(supps, num_nodes)
    return InfiniteOpt.FunctionalDiscreteMeasureData(pref, coeff_func, 0,
                                                     InfiniteOpt.All, info,
                                                     weight_func, lower_bound,
                                                     upper_bound, false)
end

# Useful error function for dependent parameters
function _ensure_independent_param(pref::InfiniteOpt.GeneralVariableRef,
                                   method)::Nothing
    if InfiniteOpt._index_type(pref) == InfiniteOpt.DependentParameterIndex
        error("Cannot generate measure data for individual dependent parameters ",
              "via `$(method)`. Try using `UniTrapezoid` for finite domains.")
    end
    return
end

# Single pref general quadrature dispatch
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
                                lower_bound::Real,
                                upper_bound::Real,
                                method::Quadrature;
                                num_nodes::Int = 3,
                                weight_func::Function = InfiniteOpt.default_weight
                                )
    _ensure_independent_param(pref, method)
    inf_bound_num = (lower_bound == -Inf) + (upper_bound == Inf)
    if inf_bound_num == 0 # finite interval
        method = FEGaussLobatto()
    elseif inf_bound_num == 1 # semi-infinite interval
        method = GaussLaguerre()
    else # infinite interval
        method = GaussHermite()
    end
    return generate_integral_data(pref, lower_bound, upper_bound, method,
                                  num_nodes = num_nodes, weight_func = weight_func)
end

## Helper methods for generating coefficients and nodes of Gauss quadrature methods
# Gauss-Legendre
function _make_nodes_weights(method::GaussLegendre, num_nodes::Int)
    return FastGaussQuadrature.gausslegendre(num_nodes)
end

# Gauss-Lobatto
function _make_nodes_weights(method::GaussLobatto, num_nodes::Int)
    return FastGaussQuadrature.gausslobatto(num_nodes)
end

# Gauss-Radau
function _make_nodes_weights(method::GaussRadau, num_nodes::Int)
    return FastGaussQuadrature.gaussradau(num_nodes)
end

# Gauss-Chebyshev
function _make_nodes_weights(method::GaussChebyshev, num_nodes::Int)
    return FastGaussQuadrature.gausschebyshev(num_nodes, method.order)
end

# Gauss-Jacobi
function _make_nodes_weights(method::GaussJacobi, num_nodes::Int)
    return FastGaussQuadrature.gaussjacobi(num_nodes, method.α, method.β)
end

# Single pref Finite Gauss Quadrature
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
                                lower_bound::Real,
                                upper_bound::Real,
                                method::FiniteGaussQuad;
                                num_nodes::Int = 3,
                                weight_func::Function = InfiniteOpt.default_weight
                                )
    _ensure_independent_param(pref, method)
    if lower_bound == -Inf || upper_bound == Inf
        @warn("The `$(typeof(method))` method can only be applied on finite intervals, " *
              "switching to an appropriate method.")
        return generate_integral_data(pref, lower_bound, upper_bound, Quadrature(),
                                      num_nodes = num_nodes,
                                      weight_func = weight_func)
    end
    (supports, coeffs) = _make_nodes_weights(method, num_nodes)
    supports = (upper_bound - lower_bound) / 2 * supports .+ (upper_bound + lower_bound) / 2
    coeffs = (upper_bound - lower_bound) / 2 * coeffs
    return InfiniteOpt.DiscreteMeasureData(pref, coeffs, supports,
                                           InfiniteOpt.generate_unique_label(),
                                           weight_func, lower_bound, upper_bound,
                                           false)
end

# Single pref Gauss-Laguerre
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
                                lower_bound::Real,
                                upper_bound::Real,
                                method::GaussLaguerre;
                                num_nodes::Int = 3,
                                weight_func::Function = InfiniteOpt.default_weight
                                )
    _ensure_independent_param(pref, method)
    if upper_bound == Inf
        if lower_bound == -Inf
            @warn("`GaussLaguerre` quadrature can only be applied on semi-infinite intervals, " *
                  "switching to an appropriate method.")
            return generate_integral_data(pref, lower_bound, upper_bound, GaussHermite(),
                                          num_nodes = num_nodes,
                                          weight_func = weight_func)
        end
        (supports, coeffs) = FastGaussQuadrature.gausslaguerre(num_nodes)
        supports .+= lower_bound
    elseif lower_bound == -Inf
        (supports, coeffs) = FastGaussQuadrature.gausslaguerre(num_nodes)
        supports = -supports .+ upper_bound
    else
        @warn("`GaussLaguerre` quadrature can only be applied on semi-infinite intervals, " *
              "switching to an appropriate method.")
        return generate_integral_data(pref, lower_bound, upper_bound, GaussLegendre(),
                                      num_nodes = num_nodes,
                                      weight_func = weight_func)
    end

    return InfiniteOpt.DiscreteMeasureData(pref, coeffs, supports,
                                           InfiniteOpt.generate_unique_label(),
                                           weight_func, lower_bound, upper_bound,
                                           false)
end

# Single pref Gauss-Hermite
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
                                lower_bound::Real,
                                upper_bound::Real,
                                method::GaussHermite;
                                num_nodes::Int = 3,
                                weight_func::Function = InfiniteOpt.default_weight
                                )
    _ensure_independent_param(pref, method)
    if lower_bound != -Inf || upper_bound != Inf
        @warn("`GaussHermite` quadrature can only be applied on infinite intervals, " *
              "switching to an appropriate method.")
        return generate_integral_data(pref, lower_bound, upper_bound, Quadrature(),
                                      num_nodes = num_nodes,
                                      weight_func = weight_func)
    end
    (supports, coeffs) = FastGaussQuadrature.gausshermite(num_nodes)
    return InfiniteOpt.DiscreteMeasureData(pref, coeffs, supports,
                                           InfiniteOpt.generate_unique_label(),
                                           weight_func, lower_bound, upper_bound,
                                           false)
end

# Single pref uniform Monte Carlo sampling
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
                                lower_bound::Real,
                                upper_bound::Real,
                                method::UniMCSampling;
                                num_supports::Int = InfiniteOpt.DefaultNumSupports,
                                weight_func::Function = InfiniteOpt.default_weight
                                )
    # check and process the arguments
    if lower_bound == -Inf || upper_bound == Inf
        error("`UniMCSampling` is not supported for (semi-)infinite intervals.")
    elseif InfiniteOpt._index_type(pref) == InfiniteOpt.DependentParameterIndex
        if !(num_supports in [0, InfiniteOpt.DefaultNumSupports])
            @warn("Cannot specify a nonzero `num_supports` when making an " *
                  "integral with repsect to an individual dependent parameter. " *
                  "Setting `num_supports = 0`.")
        end
        num_supports = 0
    end
    # make the coefficient function
    function _coeffs(supps::Vector{<:Real})::Vector{Float64}
        len = length(supps)
        return (upper_bound - lower_bound) * ones(len) / len
    end
    # prepare the data
    return InfiniteOpt.FunctionalDiscreteMeasureData(pref, _coeffs, num_supports,
                                                     InfiniteOpt.MCSample,
                                                     InfiniteOpt.NoGenerativeSupports(),
                                                     weight_func,
                                                     lower_bound, upper_bound,
                                                     false)
end

# Single pref independent Monte Carlo sampling
function generate_integral_data(pref::InfiniteOpt.GeneralVariableRef,
                                lower_bound::Real,
                                upper_bound::Real,
                                method::UniIndepMCSampling;
                                num_supports::Int = InfiniteOpt.DefaultNumSupports,
                                weight_func::Function = InfiniteOpt.default_weight
                                )
    _ensure_independent_param(pref, method)
    if lower_bound == -Inf || upper_bound == Inf
        error("`UniMCSampling` is not applicable to (semi-)infinite intervals.")
    end
    domain = InfiniteOpt.IntervalDomain(lower_bound, upper_bound)
    supports, _ = InfiniteOpt.generate_support_values(domain, InfiniteOpt.MCSample,
                                                      num_supports = num_supports)
    coeffs = (upper_bound - lower_bound) * ones(num_supports) / num_supports
    return InfiniteOpt.DiscreteMeasureData(pref, coeffs, supports,
                                           InfiniteOpt.MCSample, weight_func,
                                           lower_bound, upper_bound, false)
end

################################################################################
#                       DATA GENERATION METHODS (MULTIVARIATE)
################################################################################
# Vector of prefs with Automatic
function generate_integral_data(prefs::Vector{InfiniteOpt.GeneralVariableRef},
    lower_bounds::Vector{<:Real},
    upper_bounds::Vector{<:Real},
    method::Automatic;
    num_supports::Int = InfiniteOpt.DefaultNumSupports,
    weight_func::Function = InfiniteOpt.default_weight
    )::InfiniteOpt.AbstractMeasureData
    return generate_integral_data(prefs, lower_bounds, upper_bounds,
                                  MultiMCSampling(), num_supports = num_supports,
                                  weight_func = weight_func)
end

# Vector of prefs using Monte Carlo samples
function generate_integral_data(prefs::Vector{InfiniteOpt.GeneralVariableRef},
    lower_bounds::Vector{<:Real},
    upper_bounds::Vector{<:Real},
    method::MultiMCSampling;
    num_supports::Int = InfiniteOpt.DefaultNumSupports,
    weight_func::Function = InfiniteOpt.default_weight
    )::InfiniteOpt.AbstractMeasureData
    # check bounds
    if any(lb == -Inf for lb in lower_bounds) || any(ub == Inf for ub in upper_bounds)
        error("`MultiMCSampling` is not applicable to (semi-)infinite intervals.")
    end
    # make the coefficients function
    vol = prod(upper_bounds .- lower_bounds)
    function _coeffs(supps::Array{<:Real, 2})::Vector{Float64}
        num_samples = size(supps, 2)
        return vol * ones(num_samples) / num_samples
    end
    # make the data
    return InfiniteOpt.FunctionalDiscreteMeasureData(prefs, _coeffs,
                                                     num_supports,
                                                     InfiniteOpt.MCSample,
                                                     weight_func,
                                                     lower_bounds, upper_bounds,
                                                     false)
end

# Vector of prefs with independent Monte Carlo supports
function generate_integral_data(prefs::Vector{InfiniteOpt.GeneralVariableRef},
    lower_bounds::Vector{<:Real},
    upper_bounds::Vector{<:Real},
    method::MultiIndepMCSampling;
    num_supports::Int = InfiniteOpt.DefaultNumSupports,
    weight_func::Function = InfiniteOpt.default_weight
    )::InfiniteOpt.AbstractMeasureData
    # check bounds
    if any(lb == -Inf for lb in lower_bounds) || any(ub == Inf for ub in upper_bounds)
        error("`MultiIndepMCSampling` is not applicable to (semi-)infinite intervals.")
    end
    # prepare the supports
    domains = [InfiniteOpt.IntervalDomain(lower_bounds[i], upper_bounds[i])
            for i in eachindex(lower_bounds)]
    domain = InfiniteOpt.CollectionDomain(domains)
    supports, _ = InfiniteOpt.generate_support_values(domain, InfiniteOpt.MCSample,
                                                      num_supports = num_supports)
    # prepare the coefficients
    coeffs = prod(upper_bounds .- lower_bounds) * ones(num_supports) / num_supports
    # prepare the data
    return InfiniteOpt.DiscreteMeasureData(prefs, coeffs, supports,
                                           InfiniteOpt.MCSample, weight_func,
                                           lower_bounds, upper_bounds, false)
end

################################################################################
#                           UNIVARIATE INTEGRALS
################################################################################
# Define default keyword arguments for 1-D integrals
const UniIntegralDefaults = Dict{Symbol, Any}(:eval_method => Automatic())

"""
    uni_integral_defaults()::Dict{Symbol, Any}

Get the default keyword argument values for defining one-dimensional integrals.

```julia-repl
julia> uni_integral_defaults()
Dict{Symbol,Any} with 1 entry:
  :eval_method => Automatic()
```
"""
uni_integral_defaults()::Dict{Symbol, Any} = UniIntegralDefaults

"""
    set_uni_integral_defaults(; kwargs...)::Nothing

Set the default keyword argument settings for one-dimensional integrals.
The keyword arguments of this function will be recorded in the default keyword
argument dictionary. These will determine the default keyword argument values
when calling [`integral`](@ref MeasureToolbox.integral(::JuMP.AbstractJuMPScalar,::InfiniteOpt.GeneralVariableRef,::Real, ::Real))
with a single infinite parameter.

**Example**
```julia-repl
julia> uni_integral_defaults()
Dict{Symbol,Any} with 1 entry:
  :eval_method => Automatic()

julia> set_uni_integral_defaults(num_supports = 5, eval_method = Quadrature(),
                                 new_kwarg = true)

julia> uni_integral_defaults()
Dict{Symbol,Any} with 3 entries:
  :new_kwarg             => true
  :num_supports          => 5
  :eval_method           => Quadrature()
```
"""
function set_uni_integral_defaults(; kwargs...)::Nothing
    merge!(UniIntegralDefaults, kwargs)
    return
end

"""
    clear_uni_integral_defaults()::Nothing

Clears and resets the keyword argument defaults for univariate integrals to their
default state.

**Example**
```julia-repl
julia> uni_integral_defaults()
Dict{Symbol,Any} with 3 entries:
  :new_kwarg             => true
  :num_supports          => 5
  :eval_method           => Quadrature()

julia> clear_uni_integral_defaults()

julia> uni_integral_defaults()
Dict{Symbol,Any} with 1 entry:
  :eval_method => Automatic()
```
"""
function clear_uni_integral_defaults()::Nothing
    empty!(UniIntegralDefaults)
    UniIntegralDefaults[:eval_method] = Automatic()
    return
end
"""
    integral(expr::JuMP.AbstractJuMPScalar,
             pref::GeneralVariableRef,
             [lower_bound::Real = lower_bound(pref),
             upper_bound::Real = upper_bound(pref);
             kwargs...])::GeneralVariableRef

Returns a measure reference that evaluates the integral of `expr` with respect
to infinite parameter `pref` from `lower_bound` to `upper_bound`. This thus
considers integrals of the form: ``\\int_{p \\in P} expr(p) w(p) dp`` where ``p``
is an infinite parameter and ``w`` is the weight function is 1 by default. This
function provides a high-level interface that ultimately constructs an appropriate
concrete form of [`AbstractMeasureData`](@ref) via [`generate_integral_data`](@ref)
in accordance with the keyword arugment `eval_method` that is then used with
[`measure`](@ref). Note that it is preferred to call [`@integral`](@ref) when
`expr` is not just a single variable reference. Errors for bad bound input.

The keyword arguments are as follows:
- `eval_method::AbstractUnivariateMethod`: Used to determine the
    numerical evaluation scheme. Possible choices include:
    - [`Automatic()`](@ref MeasureToolbox.Automatic)
    - [`UniTrapezoid()`](@ref MeasureToolbox.UniTrapezoid)
    - [`UniMCSampling()`](@ref MeasureToolbox.UniMCSampling)
    - [`UniIndepMCSampling()`](@ref MeasureToolbox.UniIndepMCSampling)
    - [`Quadrature()`](@ref MeasureToolbox.Quadrature)
    - [`GaussHermite()`](@ref MeasureToolbox.GaussHermite)
    - [`GaussLegendre()`](@ref MeasureToolbox.GaussLegendre)
    - [`FEGaussLobatto()`](@ref MeasureToolbox.FEGaussLobatto)
    - [`GaussLageurre()`](@ref MeasureToolbox.GaussLaguerre)
    - [`GaussLobatto()`](@ref MeasureToolbox.GaussLobatto)
    - [`GaussChebyshev(order)`](@ref MeasureToolbox.GaussChebyshev)
    - [`GaussRadau()`](@ref MeasureToolbox.GaussRadau)
    - [`GaussJacobi(α, β)`](@ref MeasureToolbox.GaussJacobi)
- `num_supports`: The minimum number of supports to be generated (if used by
    `eval_method`)
- `weight_func`: ``w(p)`` above with parameter value inputs and scalar output

See [`set_uni_integral_defaults`](@ref) to update the default keyword argument
values for all one-dimensional integral calls.

**Example**
```julia-repl
julia> @infinite_parameter(model, x in [0, 1])
x

julia> @variable(model, f, Infinite(x))
f(x)

julia> int = integral(f, x)
∫{x ∈ [0, 1]}[f(x)]

julia> expand(int)
0.2 f(0.8236475079774124) + 0.2 f(0.9103565379264364) + 0.2 f(0.16456579813368521) + 0.2 f(0.17732884646626457) + 0.2 f(0.278880109331201)
```
"""
function integral(
    expr::JuMP.AbstractJuMPScalar,
    pref::InfiniteOpt.GeneralVariableRef,
    lower_bound::Real = NaN,
    upper_bound::Real = NaN;
    kwargs...
    )::InfiniteOpt.GeneralVariableRef
    # check parameter formatting
    InfiniteOpt._check_params(pref)
    # fill in bounds if needed
    domain = infinite_domain(pref)
    if isnan(lower_bound) && JuMP.has_lower_bound(domain)
        lower_bound = JuMP.lower_bound(pref)
    end
    if isnan(upper_bound) && JuMP.has_upper_bound(domain)
        upper_bound = JuMP.upper_bound(pref)
    end
    # ensure valid bounds
    if lower_bound >= upper_bound
        error("Invalid integral bounds, ensure that lower_bound < upper_bound.")
    elseif !InfiniteOpt.supports_in_domain([lower_bound, upper_bound], domain)
        error("Integral bounds violate the infinite domain.")
    end
    # prepare the keyword arguments and make the measure data
    processed_kwargs = merge(uni_integral_defaults(), kwargs)
    eval_method = pop!(processed_kwargs, :eval_method)
    data = generate_integral_data(pref, lower_bound, upper_bound,
                                  eval_method; processed_kwargs...)
    # make the measure
    return InfiniteOpt.measure(expr, data, name = "integral")
end

"""
    ∫(expr::JuMP.AbstractJuMPScalar,
      pref::GeneralVariableRef,
      [lower_bound::Real = NaN,
      upper_bound::Real = NaN;
      kwargs...])::GeneralVariableRef

A convenient wrapper for [`integral`](@ref). The `∫` unicode symbol is produced
via `\\int`.
"""
function ∫(
    expr::JuMP.AbstractJuMPScalar,
    pref::InfiniteOpt.GeneralVariableRef,
    lower_bound::Real = NaN,
    upper_bound::Real = NaN;
    kwargs...
    )::InfiniteOpt.GeneralVariableRef
    return integral(expr, pref, lower_bound, upper_bound; kwargs...)
end

################################################################################
#                            MULTIVARIATE INTEGRALS
################################################################################
# Define default keyword arguments for multi-D integrals
const MultiIntegralDefaults = Dict{Symbol, Any}(:eval_method => Automatic())

"""
    multi_integral_defaults()::Dict{Symbol, Any}

Get the default keyword argument values for defining multi-dimensional integrals.

```julia-repl
julia> multi_integral_defaults()
Dict{Symbol,Any} with 1 entry:
  :eval_method => Automatic()
```
"""
multi_integral_defaults()::Dict{Symbol, Any} = MultiIntegralDefaults

"""
    set_multi_integral_defaults(; kwargs...)::Nothing

Set the default keyword argument settings for multi-dimesnional integrals.
The keyword arguments of this function will be recorded in the default keyword
argument dictionary. These will determine the default keyword argument values
when calling [`integral`](@ref MeasureToolbox.integral(::JuMP.AbstractJuMPScalar,::AbstractArray{GeneralVariableRef},::Union{Real, AbstractArray{<:Real}}, ::Union{Real, AbstractArray{<:Real}}))
with an array of infinite parameters.

**Example**
```julia-repl
julia> multi_integral_defaults()
Dict{Symbol,Any} with 1 entry:
  :eval_method => Automatic()

julia> set_multi_integral_defaults(num_supports = 5, new_kwarg = true)

julia> multi_integral_defaults()
Dict{Symbol,Any} with 3 entries:
  :new_kwarg             => true
  :num_supports          => 5
  :eval_method           => Automatic()
```
"""
function set_multi_integral_defaults(; kwargs...)::Nothing
    merge!(MultiIntegralDefaults, kwargs)
    return
end

"""
    clear_multi_integral_defaults()::Nothing

Clears and resets the keyword argument defaults for multivariate integrals to their
default state.

**Example**
```julia-repl
julia> multi_integral_defaults()
Dict{Symbol,Any} with 3 entries:
  :new_kwarg             => true
  :num_supports          => 5
  :eval_method           => Automatic()

julia> clear_multi_integral_defaults()

julia> multi_integral_defaults()
Dict{Symbol,Any} with 1 entry:
  :eval_method => Automatic()
```
"""
function clear_multi_integral_defaults(; kwargs...)::Nothing
    empty!(MultiIntegralDefaults)
    MultiIntegralDefaults[:eval_method] = Automatic()
    return
end

Base.isnan(arr::AbstractArray{<:Real})::Bool = all(Base.isnan.(arr))
"""
    integral(expr::JuMP.AbstractJuMPScalar,
             prefs::AbstractArray{GeneralVariableRef},
             [lower_bounds::Union{Real, AbstractArray{<:Real}} = [lower_bound(pref)...],
             upper_bounds::Union{Real, AbstractArray{<:Real}} = [upper_bound(pref)...];
             kwargs...])::GeneralVariableRef

Returns a measure reference that evaluates the integral of `expr` with respect
to infinite parameters `prefs` from `lower_bounds` to `upper_bounds`. This thus
considers integrals of the form: ``\\int_{p \\in P} expr(p) w(p) dp`` where ``p``
is an infinite parameter and ``w`` is the weight function is 1 by default. This
function provides a high-level interface that ultimately constructs an appropriate
concrete form of [`AbstractMeasureData`](@ref) via [`generate_integral_data`](@ref)
in accordance with the keyword arugment `eval_method` that is then used with
[`measure`](@ref). Note that it is preferred to call [`@integral`](@ref) when
`expr` is not just a single variable reference. Errors when the container types
and dimensions do not match or the bounds are invalid.

The keyword arguments are as follows:
- `eval_method::AbstractMultivariateMethod`: Used to determine the
    numerical evaluation scheme. Possible choices include:
    - [`Automatic()`](@ref MeasureToolbox.Automatic)
    - [`MultiMCSampling()`](@ref MeasureToolbox.MultiMCSampling)
    - [`MultiIndepMCSampling()`](@ref MeasureToolbox.MultiIndepMCSampling)
- `num_supports`: The minimum number of supports to be generated (if used by
    `eval_method`)
- `weight_func`: ``w(p)`` above with parameter value inputs and scalar output

See [`set_multi_integral_defaults`](@ref) to update the default keyword argument
values for all multi-dimensional integral calls.

**Example**
```julia-repl
julia> @infinite_parameter(model, x[1:2] in [0, 1], independent = true);

julia> @variable(model, f, Infinite(x));

julia> int = integral(f, x)
∫{x ∈ [0, 1]^2}[f(x)]
```
"""
function integral(
    expr::JuMP.AbstractJuMPScalar,
    prefs::AbstractArray{InfiniteOpt.GeneralVariableRef},
    lower_bounds::Union{Real, AbstractArray{<:Real}} = NaN,
    upper_bounds::Union{Real, AbstractArray{<:Real}} = NaN;
    kwargs...
    )::InfiniteOpt.GeneralVariableRef
    # check parameter formatting
    InfiniteOpt._check_params(prefs)
    # fill in the lower bounds if needed
    if isnan(lower_bounds) && JuMP.has_lower_bound(first(prefs))
        lbs = map(p -> JuMP.lower_bound(p), prefs)
    elseif !isnan(lower_bounds) && lower_bounds isa Real
        lbs = map(p -> lower_bounds, prefs)
    else
        lbs = lower_bounds
    end
    # fill in the upper bounds if needed
    if isnan(upper_bounds) && JuMP.has_upper_bound(first(prefs))
        ubs = map(p -> JuMP.upper_bound(p), prefs)
    elseif !isnan(upper_bounds) && upper_bounds isa Real
        ubs = map(p -> upper_bounds, prefs)
    else
        ubs = upper_bounds
    end
    # do initial bound checks
    if InfiniteOpt._keys(prefs) != InfiniteOpt._keys(lbs) ||
       InfiniteOpt._keys(prefs) != InfiniteOpt._keys(ubs)
        error("Array keys of infinite parameters and bounds do not match.")
    end
    # process the arrays
    inds = InfiniteOpt.Collections.indices(prefs)
    vect_prefs = InfiniteOpt.Collections.vectorize(prefs, inds)
    vect_lbs = InfiniteOpt.Collections.vectorize(lbs, inds)
    vect_ubs = InfiniteOpt.Collections.vectorize(ubs, inds)
    # ensure valid bounds
    if any(vect_lbs[i] >= vect_ubs[i] for i in eachindex(vect_lbs))
        error("Invalid integral bounds, ensure that lower_bounds < upper_bounds.")
    end
    dprefs = map(p -> dispatch_variable_ref(p), vect_prefs)
    InfiniteOpt._check_bounds_in_domain(dprefs, vect_lbs, vect_ubs)
    # prepare the keyword arguments and make the measure data
    processed_kwargs = merge(multi_integral_defaults(), kwargs)
    eval_method = pop!(processed_kwargs, :eval_method)
    data = generate_integral_data(vect_prefs, vect_lbs, vect_ubs,
                                  eval_method; processed_kwargs...)
    # make the measure
    return InfiniteOpt.measure(expr, data, name = "integral")
end

"""
    ∫(expr::JuMP.AbstractJuMPScalar,
      prefs::AbstractArray{GeneralVariableRef},
      [lower_bounds::Union{Real, AbstractArray{<:Real}} = NaN,
      upper_bounds::Union{Real, AbstractArray{<:Real}} = NaN;
      kwargs...])::GeneralVariableRef

A convenient wrapper for [`integral`](@ref). The unicode symbol `∫` is produced
via `\\int`.
"""
function ∫(expr::JuMP.AbstractJuMPScalar,
           prefs::AbstractArray{InfiniteOpt.GeneralVariableRef},
           lower_bounds::Union{Real, AbstractArray{<:Real}} = NaN,
           upper_bounds::Union{Real, AbstractArray{<:Real}} = NaN;
           kwargs...)::InfiniteOpt.GeneralVariableRef
    return integral(expr, prefs, lower_bounds, upper_bounds; kwargs...)
end

"""
    @integral(expr::JuMP.AbstractJuMPScalar,
              prefs::Union{GeneralVariableRef, AbstractArray{GeneralVariableRef}},
              [lower_bounds::Union{Real, AbstractArray{<:Real}} = default_bounds,
              upper_bounds::Union{Real, AbstractArray{<:Real}} = default_bounds;
              kwargs...])::GeneralVariableRef

An efficient wrapper for [`integral`](@ref integral(::JuMP.AbstractJuMPScalar, ::InfiniteOpt.GeneralVariableRef, ::Real, ::Real))
and [`integral`](@ref integral(::JuMP.AbstractJuMPScalar, ::AbstractArray{InfiniteOpt.GeneralVariableRef}, ::Union{Real, AbstractArray{<:Real}}, ::Union{Real, AbstractArray{<:Real}})).
Please see the above doc strings for more information.
"""
macro integral(expr, prefs, args...)
    _error(str...) = InfiniteOpt._macro_error(:integral, (expr, prefs, args...),
                                              str...)
    extra, kwargs, _, _ = InfiniteOpt._extract_kwargs(args)
    if length(extra) != 0 && length(extra) != 2
        _error("Incorrect number of positional arguments for @integral. " *
               "Must provide both bounds or no bounds.")
    end
    expression = :( JuMP.@expression(InfiniteOpt._Model, $expr) )
    mref = :( integral($expression, $prefs, $(extra...); ($(kwargs...))) )
    return esc(mref)
end

"""
    @∫(expr::JuMP.AbstractJuMPScalar,
       prefs::Union{GeneralVariableRef, AbstractArray{GeneralVariableRef}},
       [lower_bounds::Union{Real, AbstractArray{<:Real}} = default_bounds,
       upper_bounds::Union{Real, AbstractArray{<:Real}} = default_bounds;
       kwargs...])::GeneralVariableRef

A convenient wrapper for [`@integral`](@ref). The unicode symbol `∫` is produced
via `\\int`.
"""
macro ∫(expr, prefs, args...)
    _error(str...) = InfiniteOpt._macro_error(:∫, (expr, prefs, args...), str...)
    extra, kwargs, _, _ = InfiniteOpt._extract_kwargs(args)
    if length(extra) != 0 && length(extra) != 2
        _error("Incorrect number of positional arguments for @integral. " *
               "Must provide both bounds or no bounds.")
    end
    expression = :( JuMP.@expression(InfiniteOpt._Model, $expr) )
    mref = :( integral($expression, $prefs, $(extra...); ($(kwargs...))) )
    return esc(mref)
end


text:
export BlocksPrism
function BlocksPrism(PrismData,c2::Float64,d::Float64,δ::Float64)
    PrA = PrismData[1]
    sIndex = PrismData[2]

    TempTetra = Array{Array}(undef,0)
    TempPyra = Array{Array}(undef,0)
    TempAmp = Array{Real}(undef,0)
    PosAmp = 1
    for i in sIndex
        (ja1,jb1,jd1t,je1t,jf1,ja2,jb2,jd2t,je2,jf2,jc,jg) = InverseSuperIndex(i,12)
        if jd1t == jd1 && je1t == je1 && jd2t == jd2
            push!(TempAmp,PrA[PosAmp])
            push!(TempTetra,[ja1 jb1 jg])
            push!(TempPyra,[jf1 ja2 jb2 je2 jf2 jc])
        end
        PosAmp += 1
    end
    PosTetra = unique(TempTetra)
    PosPyra = unique(TempPyra)
    SizeA = length(PosTetra)
    SizeB = length(PosPyra)
    mat = zeros(SizeA,SizeB)
    for i in 1:length(TempTetra)
        qa = findall(x -> x == TempTetra[i], PosTetra)
        qb = findall(x -> x == TempPyra[i], PosPyra)
        mat[qa[1],qb[1]] = TempAmp[i]
    end

    if SizeA == 0
        mat = 0
    end
    return mat, PosTetra, PosPyra
end

#### SVD embedding map, by extracting the edges we don't want to keep

export EmbeddingMap1
function EmbeddingMap1(PrismData,trunc,truncEmbedding1)
	sIndex,Amp = PrismEff(PrismData,trunc)

	#first embedding map, SVD on j1
	TempLine = Array{Real}(undef,0)
	TempRest = Array{Real}(undef,0)
	TempAmp = Array{Real}(undef,0)
	pos = 1
	for i in sIndex
		(β,α1,α2,α3,d,α1t,α2t,α3t,dt,j1,q,j2,u,ut,v,vt,q1,q1t) = InverseSuperIndex(i,18)
		push!(TempAmp,Amp[pos])
		push!(TempLine,j1)
		push!(TempRest,SuperIndex([β,α1,α2,α3,d,α1t,α2t,α3t,dt,q,j2,u,ut,v,vt,q1,q1t]))
		pos +=1
	end

	PosLine = unique(TempLine)
	PosRest = unique(TempRest)
	SizeL = length(PosLine)
	SizeR = length(PosRest)
	mat = zeros(SizeR,SizeL)
	for i in 1:length(TempLine)
        qb = findall(x -> x == TempLine[i], PosLine)
        qa = findall(x -> x == TempRest[i], PosRest)
        mat[qa[1],qb[1]] = TempAmp[i]
    end
    if SizeL == 0
        mat = 0
    end

	U, s, V = svd(mat)
	V = adjoint(V)
	sEmb1 = numchop(s)
	AmpEmb1 = U[:,1:truncEmbedding1]*sqrt.(sEmb1[1:truncEmbedding1])
	TruncVEmb1 = transpose(sqrt.(sEmb1[1:truncEmbedding1]))*V[1:truncEmbedding1,:]

	# second embedding, on j2,  should be equivalent to the first one. multiplying by the V
	TempLine = Array{Real}(undef,0)
	TempRest = Array{Real}(undef,0)
	TempAmp = Array{Real}(undef,0)
	pos = 1
	for i in PosRest
		(β,α1,α2,α3,d,α1t,α2t,α3t,dt,q,j2,u,ut,v,vt,q1,q1t) = InverseSuperIndex(i,17)
		push!(TempAmp,Amp[pos])
		push!(TempLine,j2)
		push!(TempRest,SuperIndex([β,α1,α2,α3,d,α1t,α2t,α3t,dt,q,u,ut,v,vt,q1,q1t]))
		pos +=1
	end

	PosLine = unique(TempLine)
	PosRest = unique(TempRest)
	SizeL = length(PosLine)
	SizeR = length(PosRest)
	mat = zeros(SizeR,SizeL)
	for i in 1:length(TempLine)
		qb = findall(x -> x == TempLine[i], PosLine)
		qa = findall(x -> x == TempRest[i], PosRest)
		mat[qa[1],qb[1]] = TempAmp[i]
	end
	if SizeL == 0
		mat = 0
	end

	U, s, V = svd(mat)
	V = adjoint(V)
	sEmb2 = numchop(s)
	AmpEmb2 = U[:,1:truncEmbedding1]*sqrt.(sEmb2[1:truncEmbedding1])
	TruncVEmb2 = transpose(sqrt.(sEmb2[1:truncEmbedding1]))*V[1:truncEmbedding1,:]

	return AmpEmb2, PosRest, TruncVEmb1,TruncVEmb2, sEmb1,sEmb2

end

export EmbeddingMap1Commutator
function EmbeddingMap1Commutator(PrismData,trunc,truncEmbedding1)
	sIndex,Amp = PrismEff(PrismData,trunc)

	#first embedding map, SVD on j1
	TempLine = Array{Real}(undef,0)
	TempRest = Array{Real}(undef,0)
	TempAmp = Array{Real}(undef,0)
	pos = 1
	for i in sIndex
		(β,α1,α2,α3,d,α1t,α2t,α3t,dt,j1,q,j2,u,ut,v,vt,q1,q1t) = InverseSuperIndex(i,18)
		push!(TempAmp,Amp[pos])
		push!(TempLine,j2)
		push!(TempRest,SuperIndex([β,α1,α2,α3,d,α1t,α2t,α3t,dt,j1,q,u,ut,v,vt,q1,q1t]))
		pos +=1
	end

	PosLine = unique(TempLine)
	PosRest = unique(TempRest)
	SizeL = length(PosLine)
	SizeR = length(PosRest)
	mat = zeros(SizeR,SizeL)
	for i in 1:length(TempLine)
        qb = findall(x -> x == TempLine[i], PosLine)
        qa = findall(x -> x == TempRest[i], PosRest)
        mat[qa[1],qb[1]] = TempAmp[i]
    end
    if SizeL == 0
        mat = 0
    end

	U, s, V = svd(mat)
	V = adjoint(V)
	sEmb1 = numchop(s)
	AmpEmb1 = U[:,1:truncEmbedding1]*sqrt.(sEmb1[1:truncEmbedding1])
	TruncVEmb1 = transpose(sqrt.(sEmb1[1:truncEmbedding1]))*V[1:truncEmbedding1,:]

	# second embedding, on j2,  should be equivalent to the first one. multiplying by the V
	TempLine = Array{Real}(undef,0)
	TempRest = Array{Real}(undef,0)
	TempAmp = Array{Real}(undef,0)
	pos = 1
	for i in PosRest
		(β,α1,α2,α3,d,α1t,α2t,α3t,dt,j1,q,u,ut,v,vt,q1,q1t) = InverseSuperIndex(i,17)
		push!(TempAmp,Amp[pos])
		push!(TempLine,j1)
		push!(TempRest,SuperIndex([β,α1,α2,α3,d,α1t,α2t,α3t,dt,q,u,ut,v,vt,q1,q1t]))
		pos +=1
	end

	PosLine = unique(TempLine)
	PosRest = unique(TempRest)
	SizeL = length(PosLine)
	SizeR = length(PosRest)
	mat = zeros(SizeR,SizeL)
	for i in 1:length(TempLine)
		qb = findall(x -> x == TempLine[i], PosLine)
		qa = findall(x -> x == TempRest[i], PosRest)
		mat[qa[1],qb[1]] = TempAmp[i]
	end
	if SizeL == 0
		mat = 0
	end

	U, s, V = svd(mat)
	V = adjoint(V)
	sEmb2 = numchop(s)
	AmpEmb2 = U[:,1:truncEmbedding1]*sqrt.(sEmb2[1:truncEmbedding1])
	TruncVEmb2 = transpose(sqrt.(sEmb2[1:truncEmbedding1]))*V[1:truncEmbedding1,:]

	return AmpEmb2, PosRest, TruncVEmb1,TruncVEmb2, sEmb1,sEmb2

end

function Move2_2diag(PrismData,trunc)
	sIndex,Amp = PrismEff(PrismData,trunc)

	NewData = Array{Real}(undef,0)
	NewsIndex = Array{Real}(undef,0)

	TempIndexWithoutqWithqMove = Array{Real}(undef,0)
	TempAmp = Array{Real}(undef,0)

	for i in 1:length(sIndex)
		(β,α1,α2,α3,d,α1t,α2t,α3t,dt,j1,q,j2,u,ut,v,vt,q1,q1t) = InverseSuperIndex(sIndex[i],18)
		for qMove in 0:0.5:y
			if delta(qMove,q1,ut)==1 && delta(qMove,q1t,v)==1 && delta(qMove,α3,dt) ==1
				push!(TempIndexWithoutqWithqMove,SuperIndex([β,α1,α2,α3,d,α1t,α2t,α3t,dt,j1,j2,u,ut,v,vt,q1,q1t,qMove]))
				push!(TempAmp,Fsymb(ut,q1t,q,v,q1,qMove)*Amp[i])
			end
		end
	end

	UniqueTempIndexWithoutqWithqMove = unique(TempIndexWithoutqWithqMove)
	for i in 1:length(UniqueTempIndexWithoutqWithqMove)
		pos = findall(j-> j == UniqueTempIndexWithoutqWithqMove[i],TempIndexWithoutqWithqMove)
		push!(NewData,sum(TempAmp[pos]))
	end

#	p = sortperm(vec(TempIndexWithoutqWithqMove))
#	TempIndexWithoutqWithqMove = TempIndexWithoutqWithqMove[p]
#	TempPos = TempPos[p]
#
#	UniqueTempIndexWithoutqWithqMove = unique(TempIndexWithoutqWithqMove)
#	CountUniqueTempIndexWithoutqWithqMove = Array{Real}(undef,length(UniqueTempIndexWithoutqWithqMove))
#	for i in 1:length(UniqueTempIndexWithoutqWithqMove)
#		CountUniqueTempIndexWithoutqWithqMove[i] = count(j->j == UniqueTempIndexWithoutqWithqMove[i],TempIndexWithoutqWithqMove)
#	end

	# Comparaison = 0
	# for i in TempIndexWithoutqWithqMove
	# 	if Comparaison == i
	# 		(β,α1,α2,α3,d,α1t,α2t,α3t,dt,j1,j2,u,ut,v,vt,q1,q1t,qMove) = InverseSuperIndex(i,18)
	# 		sol += Fsymb(ut,q1t,qa,v,q1,qMove)*Amp[j]
	# 	end
	# 	Comparaison = i
	# end
	#
	# for i in UniqueTempIndexWithoutqWithqMove
	# 	sol = 0
	# 	cnt = 1
	# 	if TempIndexWithoutqWithqMove[cnt] == i
	# 		sol += Fsymb(ut,q1t,qa,v,q1,qMove)*Amp[j]

	# for qMove in 0:0.5:y, i in sIndex
	# 	(β,α1,α2,α3,d,α1t,α2t,α3t,dt,j1,q,j2,u,ut,v,vt,q1,q1t) = InverseSuperIndex(i,18)
	# 	if delta(qMove,q1,ut)==1 && delta(qMove,q1t,v)==1 && delta(qMove,α3,dt) ==1
	# 		sol = 0
	# 		for j in 1:length(sIndex)
	# 			(βa,α1a,α2a,α3a,da,α1ta,α2ta,α3ta,dta,j1a,qa,j2a,ua,uta,va,vta,q1a,q1ta) = InverseSuperIndex(sIndex[j],18)
	# 			if β==βa && α1a==α1 && α2a==α2 && α3==α3a && da==d && α1ta==α1t && α2ta==α2t && α3ta==α3ta && dta==dt && j1a==j1 && j2a==j2 && ua==u && uta==ut && va==v && vta==vt && q1a==q1 && q1ta==q1t
	# 				sol += Fsymb(ut,q1t,qa,v,q1,qMove)*Amp[j]
	# 			end
	# 		end
	# 		push!(NewData,sol)
	# 		push!(NewsIndex,SuperIndex([β,α1,α2,α3,d,α1t,α2t,α3t,dt,j1,qMove,j2,u,ut,v,vt,q1,q1t]))
	# 	end
	# end

	return NewData,UniqueTempIndexWithoutqWithqMove
end

export PrismUVTruncated
function PrismUVTruncated(PrismData,trunc)
	FullU, FullV, FullTetraindex, FullPyraindex, FullCindex,s = FullsvdPrismTruncated(PrismData,trunc)

	AmpUV = Array{Float64}(undef,0)
	PosUVFaceExt = Array{Int64}(undef,0)
	PosUVFaceTrue = Array{Int64}(undef,0)
	PosUVFaceFalse = Array{Int64}(undef,0)
	SizeC = length(FullCindex)
	for i in 1:SizeC, j in 1:SizeC
		(dA,bA,aA) = FullCindex[i]
		(dB,bB,aB) = FullCindex[j]
		#glue along dA/dB
		if dA == dB
			TempA = FullTetraindex[i]
			TempB = FullPyraindex[j]
			for iA in 1:length(TempA), jB in 1:length(TempB)
				(j1A,uA,vA) = TempA[iA]
				(j2B,uB,j3B,j4B,vB,j5B) = TempB[jB]
				#and glue along uA/uB and vA/vB
				if uA == uB && vA == vB
					push!(AmpUV,FullU[i][iA]*FullV[j][jB]/(visqrt(dA)*visqrt(uA)*visqrt(vB)))
					push!(PosUV12,SuperIndex([vA,j2B,j5B]))
					push!(PosUVFace1,SuperIndex([bA,j1A,j4B,j3B,uA]))
					push!(PosUVFace2,SuperIndex([bA,aA,bB,aB,dA]))
				end
			end
		end
	end
	return AmpUV,PosUV12,PosUVFace1,PosUVFace2, s
end

ext: jl


text:
import Dates
function save(fname,L::Solution)
  date=string(Dates.now(Dates.UTC))
  open(fname,"w") do f
    write(f,"# Solution version 0\n")
    write(f,"#"*date*"\n")
    write(f,"params=[")
    for (key,value) in L.params
      write(f,"(:$(string(key)),$value),\n")
    end
    write(f,"]\n")
    write(f,"eigval=:$(L.eigval)\n")

    write(f,"v=")
    vector_write(f,L.v)
    write(f,"]\n")

    write(f,"v_adj=")
    vector_write(f,L.v_adj)

    write(f,"[eigval_pert]\n")
    for (key,value) in L.eigval_pert
      write(f,"\t[eigval_pert.$key]\n")
      if typeof(value)<:Tuple
        write(f,"\t\tnum=")
        vector_write(f,value[1])
        write(f,"\t\tden=")
        vector_write(f,value[2])
      else
        write(f,"\t\tnum=")
        vector_write(f,value)
      end
    end

    write(f,"[v_pert]\n")
    for (key,value) in L.v_pert
      write(f,"\t[v_pert.$key]\n")
      if typeof(value)<:Tuple
        write(f,"\t\t[v_pert.$key.num]\n")
        for (idx,val) in enumerate(value[1])
          write(f,"\t\t\t[v_pert.$key.num.$idx]\n")
          write(f,"\t\t\tv=")
          vector_write(f,val)
        end

        write(f,"\t\t[v_pert.$key.den]\n")
        for (idx,val) in enumerate(value[2])
          write(f,"\t\t\t[v_pert.$key.den.$idx]\n")
          write(f,"\t\t\tv=")
          vector_write(f,val)
        end
      else
        write(f,"\t\t[v_pert.$key.num]\n")
        for (idx,val) in enumerate(value)
          write(f,"\t\t\t[v_pert.$key.num.$idx]\n")
          write(f,"\t\t\tv=")
          vector_write(f,val)
        end
      end
    end
    # for (idx,val) in enumerate(L.v_pert)
    #   write(f,"\t[v_pert.$idx]\n")
    #   write(f,"\t\tv=")
    #   vector_write(f,val)
    # end
  end
end
#TODO: save v_pert
function vector_write(f,V)
  write(f,"[")
  for v in V #TODO: mehr Spalten
    if imag(v)>=0
      write(f,"$(real(v))+$(imag(v))im,")
    else
      write(f,"$(real(v))$(imag(v))im,")
    end
  end
  write(f,"]\n")
  return
end

#include("toml.jl")

function read_sol(f)
    D=read_toml(f)
    eigval_pert=Dict{Symbol,Any}()
    for (key,value) in D["/eigval_pert"]
      if haskey(value,"den")
        eigval_pert[Symbol(key[2:end])]=(value["num"],value["den"])
      else
        eigval_pert[Symbol(key[2:end])]=value["num"]
      end
    end
    params=Dict{Symbol,Complex{Float64}}()
    for (sym,val) in D["params"]
      params[sym]=val
    end
    eigval=D["eigval"]
    v=D["v"]
    v_adj=D["v_adj"]

    v_pert=Dict{Symbol,Any}()
    for (key,value) in D["/v_pert"]
      if haskey(value,"/den")
        N = length(value["/den"])
        B=Array{Array{Complex{Float64}},1}(undef,N)
        for idx=1:N
          B[idx] = value["/den"]["/$idx"]["v"]
        end
        N = length(value["/num"])
        A=Array{Array{Complex{Float64}},1}(undef,N)
        for idx=1:N
          A[idx] = value["/num"]["/$idx"]["v"]
        end
        v_pert[Symbol(key[2:end])]=A,B
      else
        N=length(value["/num"])
        V=Array{Array{Complex{Float64}},1}(undef,N)
        for idx=1:N
          V[idx] = value["/num"]["/$idx"]["v"]
        end
        v_pert[Symbol(key[2:end])]=V
      end
    end
    #TODO: proper type initialization
    # N=length(D["/v_pert"])
    # v_pert=Array{Array{Complex{Float64}},1}(undef,N)
    # for idx=1:N
    #   v_pert[idx]=D["/v_pert"]["/$idx"]["v"]
    # end
    return Solution(params,v,v_adj,eigval,eigval_pert,v_pert)
end


text:
function update_state!(
    status::State{xfgh_indiv},
    parameters::ECA,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
)
    K = parameters.K
    I = randperm(parameters.N)
    N = parameters.N
    population = status.population
    D = size(problem.bounds, 2)

    a = problem.bounds[1, :]
    b = problem.bounds[2, :]

    ε = 0.0

    p = status.f_calls / options.f_calls_limit

    if parameters.ε > 0.0
        ε = parameters.ε * (p - 1.0)^4
    end

    feasible_solutions = findall( s->s.is_feasible, status.population )

    fs = fvals(status)
    fx_max = maximum( abs.(fs) )

    vios = map(sum_violations, population)
    fs += 2.0*fx_max*vios
    γ = maximum( abs.(fs) )
    weights = 2.0γ .- fs

    # For each elements in Population
    for i = 1:parameters.N

        # current
        x = status.population[i].x

        # generate U masses
        U_ids = getU_ids(parameters.K, I, i, parameters.N, feasible_solutions)
        U = population[U_ids]

        # generate center of mass
        mass = weights[U_ids]
        c = center(U, mass)
        u_worst = argmin(mass)
        u_best = argmax(mass)

        # stepsize
        η = parameters.η_max * rand()

        # u: worst element in U
        u = U[u_worst].x
        v = U[u_best].x

        # current-to-center/bin
        y = x .+ η .* (c .- u)

        mask = rand(length(y)) .< 1.0 / length(y)
        y[mask] = v[mask]

        evo_boundary_repairer!(y, c, problem.bounds)

        sol = create_solution(y, problem, ε=ε)
        status.f_calls += 1

        if is_better(sol, status.best_sol)
            status.best_sol = sol
        end

        # stored but not used until replacement step
        push!(status.population, sol)

        status.stop = stop_check(status, information, options)
        status.stop && break
    end

    if length(status.population) == N
        # no resize population
        return
    end

    # replacement step
    sort!(status.population, lt = is_better)
    deleteat!(status.population, N+1:length(status.population))

end

function final_stage!(
    status::State{xfgh_indiv},
    parameters::ECA,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
)
    status.final_time = time()

end


text:
import DrWatson: quickactivate
quickactivate(@__DIR__, "Chemostat_Heerden2013")

@time begin
    import SparseArrays
    using Base.Threads
    using Serialization

    # -------------------------------------------------------------------
    # Run add https://github.com/josePereiro/Chemostat_Heerden2013.jl in
    # the Julia Pkg REPL to install the package, then you must activate
    # the package enviroment (see README)
    import Chemostat_Heerden2013
    const ChHd = Chemostat_Heerden2013
    const Hd  = ChHd.HeerdenData;
    const BD  = ChHd.BegData;
    const iJR = ChHd.iJR904

    # -------------------------------------------------------------------
    # run add "https://github.com/josePereiro/Chemostat" in the
    # julia Pkg REPL for installing the package
    import Chemostat
    import Chemostat.LP.MathProgBase
    const Ch = Chemostat
    const ChU = Ch.Utils
    const ChSS = Ch.SteadyState
    const ChLP = Ch.LP
    const ChEP = Ch.MaxEntEP
    const ChSU = Ch.SimulationUtils

    import UtilsJL
    const UJL = UtilsJL
    UJL.set_cache_dir(iJR.cachedir())
end

## ----------------------------------------------------------------------------
# globals
const WLOCK = ReentrantLock()
const DAT_FILE_PREFFIX = "maxent_ep_dat"

## ----------------------------------------------------------------------------
const ME_Z_OPEN_G_OPEN          = :ME_Z_OPEN_G_OPEN           # Do not use extra constraints
const ME_MAX_POL                = :ME_MAX_POL                 #
const ME_MAX_POL_B0             = :ME_MAX_POL_B0                 #
const ME_Z_EXPECTED_G_EXPECTED  = :ME_Z_EXPECTED_G_EXPECTED    # Match ME and Dy biom average and constraint av_ug
const ME_Z_EXPECTED_G_BOUNDED   = :ME_Z_EXPECTED_G_BOUNDED    # Match ME and Dy biom average and constraint av_ug
const ME_Z_EXPECTED_G_MOVING    = :ME_Z_EXPECTED_G_MOVING     #
const ME_Z_FIXXED_G_BOUNDED     = :ME_Z_FIXXED_G_BOUNDED      # Fix biom around observed

## -------------------------------------------------------------------
function dat_file(;kwargs...)
    fname = UJL.mysavename(DAT_FILE_PREFFIX, "jls"; kwargs...)
    iJR.procdir(fname)
end

function check_cache(;kwargs...)
    thid = threadid()
    datfile = dat_file(;kwargs...)
    isfile(datfile) && lock(WLOCK) do
        @info("Cached loaded (skipping)",
            datfile, thid
        ); println()
    end; isfile(datfile)
end

# ## -------------------------------------------------------------------
# # ME_Z_EXPECTED_G_MOVING
# include("2.0.1_ME_Z_EXPECTED_G_MOVING.jl")

# ## -------------------------------------------------------------------
# # ME_Z_EXPECTED_G_BOUNDED
# include("2.0.2_ME_Z_EXPECTED_G_BOUNDED.jl")

# ## -------------------------------------------------------------------
# # ME_Z_FIXXED_G_BOUNDED
# include("2.0.3_ME_Z_FIXXED_G_BOUNDED.jl")

# ## -------------------------------------------------------------------
# # ME_Z_OPEN_G_OPEN
# # It was computed in ME_Z_EXPECTED_G_BOUNDED
# include("2.0.4_ME_Z_OPEN_G_OPEN.jl")

## -------------------------------------------------------------------
# ME_MAX_POL
include("2.0.5_ME_MAX_POL.jl")

## ----------------------------------------------------------------------------
# ME_MAX_POL_B0
include("2.0.7_ME_MAX_POL_B0.jl")


text:
# Imaginary Time Evolution
# I don't think this works!
# #+HTML: Source
# #+HTML:

# [[file:~/.julia/dev/MatrixProductStates/README.org::*Imaginary%20Time%20Evolution][Imaginary Time Evolution:1]]
# Fixme! this does not appear to find ground states!

function _MPO_handed_time_evolver_odd(hs::Vector{Matrix{T}}, τ, L, d, reim) where {T}
    tensors = Array{T, 4}[]

    if iseven(L)
        I=reshape([1 0; 0 1],(1,1,2,2))
        push!(tensors,I)
    end

    for h in hs
        if reim =="imag"
            O = exp(-τ*h)
        elseif reim=="real"
            O = exp(-1.0im*τ*h)
        end
        @cast P[(σⁱ, σⁱ′), (σⁱ⁺¹, σⁱ⁺¹′)] |= O[(σⁱ, σⁱ⁺¹), (σⁱ′, σⁱ⁺¹′)] (σⁱ:d, σⁱ′:d)
        U, S, V = svd(P)

        @cast U[1, k, σⁱ, σⁱ′]     := U[(σⁱ, σⁱ′), k] * √(S[k])      (σⁱ:d)
        @cast Ū[k, 1, σⁱ⁺¹, σⁱ⁺¹′] := √(S[k]) * V'[k, (σⁱ⁺¹, σⁱ⁺¹′)] (σⁱ⁺¹:d)
        push!(tensors, U, Ū)
    end

    if isodd(L)
        I=reshape([1 0; 0 1],(1,1,2,2))
        push!(tensors,I)
    end

    MPO{L, T}(tensors)
end

function _MPO_handed_time_evolver_even(hs::Vector{Matrix{T}}, τ, L, d, reim) where {T}
    tensors = Array{T, 4}[]

    if isodd(L)
        I=reshape([1 0; 0 1],(1,1,2,2))
        push!(tensors,I)
    end

    for h in hs
        if reim =="imag"
            O = exp(-τ*h)
        elseif reim=="real"
            O = exp(-1.0im*τ*h)
        end
        @cast P[(σⁱ, σⁱ′), (σⁱ⁺¹, σⁱ⁺¹′)] |= O[(σⁱ, σⁱ⁺¹), (σⁱ′, σⁱ⁺¹′)] (σⁱ:d, σⁱ′:d)
        U, S, V = svd(P)

        @cast U[1, k, σⁱ, σⁱ′]     := U[(σⁱ, σⁱ′), k] * √(S[k])      (σⁱ:d)
        @cast Ū[k, 1, σⁱ⁺¹, σⁱ⁺¹′] := √(S[k]) * V'[k, (σⁱ⁺¹, σⁱ⁺¹′)] (σⁱ⁺¹:d)
        push!(tensors, U, Ū)
    end

    if iseven(L)
        I=reshape([1 0; 0 1],(1,1,2,2))
        push!(tensors,I)
    end

    MPO{L, T}(tensors)
end

function MPO_time_evolvers(h1::Matrix, hi::Matrix, hL::Matrix, τ, L, d, reim)

    if iseven(L)
        odd_hs  = [h1, [hi for _ in 3:2:(L-2)]...,hL]
        even_hs = [hi for i in 2:2:(L-2)]
    else
        odd_hs  = [h1, [hi for _ in 3:2:(L-2)]...] #size=3
        even_hs = [[hi for i in 2:2:(L-2)]...,hL] #size=3
    end

    Uodd  = _MPO_handed_time_evolver_odd(odd_hs, τ, L, d, reim)
    Ueven = _MPO_handed_time_evolver_even(even_hs, τ, L, d, reim)
    Uodd, Ueven
end

function imag_time_evolution(ψ::MPS{L, T}, h1::Matrix{T}, hi::Matrix{T}, hL::Matrix{T}, H::MPO{L,T}, β, N, Dcut) where {L, T}
    @warn "This probably still doesn't work!"
    τ = β/N
    d = length(ψ[1][1, 1, :])
    ϕ = ψ  # Ground state guess
    dir = left
    Uodd, Ueven = MPO_time_evolvers(h1, hi, hL, τ, L, d, "imag")
    for i in 1:N
        println("inverse temperature = ",i*τ)
        println((ϕ'*(H*ϕ))/ϕ'ϕ)

        ϕ1, dir = compress(Uodd  * ϕ,  dir, Dcut=Dcut)
        ϕ,  dir = compress(Ueven * ϕ1, dir, Dcut=Dcut)
        #ϕ,  dir = compress(Uodd  * ϕ2, dir, Dcut=Dcut)
    end
    ϕ
end

function real_time_evolution(ψ::MPS{L, T}, h1::Matrix{T}, hi::Matrix{T}, hL::Matrix{T}, H::MPO{L,T}, Time, N, Dcut) where {L, T}
    @warn "This probably still doesn't work!"
    τ = Time/N
    d = length(ψ[1][1, 1, :])
    ϕb = ψ
    dir = left
    Uodd, Ueven = MPO_time_evolvers(h1, hi, hL, τ, L, d, "real")
    time_list=Any[]
    expect_list=Any[]
    for i in 1:N
        t = i * τ
        expect = (ϕb' * (H * ϕb)) / ϕb'ϕb
        println("time = ",t)
        println(expect)
        push!(time_list,t)
        push!(expect_list,expect)

        Uoddϕ = Uodd * ϕb
        ϕa, dir = compress(Uoddϕ,  dir, Dcut=Dcut)
        #ϕb, dir = variational_compress(ϕa, dir, Uoddϕ)
        #ϕa, dir = variational_compress(ϕb, dir, Uoddϕ)
        Uevenϕ = Ueven * ϕa
        ϕb,  dir = compress(Uevenϕ, dir, Dcut=Dcut)
        #ϕa, dir = variational_compress(ϕb, dir, Uoddϕ)
        #ϕb, dir = variational_compress(ϕa, dir, Uoddϕ)
    end
    ϕb, time_list, expect_list
end

function real_timedep_time_evolution(ψ::MPS{L, T}, h1::Vector{Matrix{T}}, hi::Vector{Matrix{T}}, hL::Vector{Matrix{T}}, H::MPO{L,T}, Time, N, Dcut) where {L, T}
    @warn "This probably still doesn't work!"
    τ = Time/N
    d = length(ψ[1][1, 1, :])
    ϕb = ψ
    dir = left
    time_list=Any[]
    expect_list=Any[]
    for i in 1:Int(N)
        println(i)
        Uodd, Ueven = MPO_time_evolvers(h1[i], hi[i], hL[i], τ, L, d, "real")
        t = i * τ
        expect = (ϕb' * (H * ϕb)) / ϕb'ϕb
        println("time = ",t)
        println(expect)
        push!(time_list,t)
        push!(expect_list,expect)

        Uoddϕ = Uodd * ϕb
        ϕa, dir = compress(Uoddϕ,  dir, Dcut=Dcut)
        ϕb, dir = variational_compress(ϕa, dir, Uoddϕ)
        ϕa, dir = variational_compress(ϕb, dir, Uoddϕ)
        Uevenϕ = Ueven * ϕa
        ϕb, dir = compress(Uevenϕ, dir, Dcut=Dcut)
        ϕa, dir = variational_compress(ϕb, dir, Uevenϕ)
        ϕb, dir = variational_compress(ϕa, dir, Uevenϕ)
    end
    ϕb, time_list, expect_list
end

# Time Evolution:1 ends here


text:
module BinnedPool

# binned memory pool allocator
#
# the core design is a pretty simple:
# - bin allocations into multiple pools according to their size (see `poolidx`)
# - when requested memory, check the pool for unused memory, or allocate dynamically
# - conversely, when released memory, put it in the appropriate pool for future use
#
# to avoid memory hogging and/or trashing the Julia GC:
# - keep track of used and available memory, in order to determine the usage of each pool
# - keep track of each pool's usage, as well as a window of previous usages
# - regularly release memory from underused pools (see `reclaim(false)`)
#
# possible improvements:
# - context management: either switch contexts when performing memory operations,
#                       or just use unified memory for all allocations.
# - per-device pools

# TODO: move the management thread one level up, to be shared by all allocators

using ..CuArrays
using ..CuArrays: @pool_timeit

using CUDAdrv

using Base: @lock
using Base.Threads: SpinLock

## tunables

const MAX_POOL = 2^27 # 128 MiB

const USAGE_WINDOW = 5

# min and max time between successive background task iterations.
# when the pool usages don't change, scan less regularly.
#
# together with USAGE_WINDOW, this determines how long it takes for objects to get reclaimed
const MIN_DELAY = 1.0
const MAX_DELAY = 5.0

## block of memory

struct Block
    ptr::CuPtr{Nothing}
    sz::Int
end

Base.pointer(block::Block) = block.ptr
Base.sizeof(block::Block) = block.sz

@inline function actual_alloc(sz)
    ptr = CuArrays.actual_alloc(sz)
    block = ptr === nothing ? nothing : Block(ptr, sz)
end

function actual_free(block::Block)
    CuArrays.actual_free(pointer(block))
    return
end

## infrastructure

const pool_lock = ReentrantLock()

const pools_used = Vector{Set{Block}}()
const pools_avail = Vector{Vector{Block}}()

poolidx(n) = ceil(Int, log2(n))+1
poolsize(idx) = 2^(idx-1)

@assert poolsize(poolidx(MAX_POOL)) <= MAX_POOL "MAX_POOL cutoff should close a pool"

function create_pools(idx)
  if length(pool_usage) >= idx
    # fast-path without taking a lock
    return
  end

  @lock pool_lock begin
    while length(pool_usage) < idx
      push!(pool_usage, 1)
      push!(pool_history, initial_usage)
      push!(pools_used, Set{Block}())
      push!(pools_avail, Vector{Block}())
    end
  end
end

## pooling

const initial_usage = Tuple(1 for _ in 1:USAGE_WINDOW)

const pool_usage = Vector{Float64}()
const pool_history = Vector{NTuple{USAGE_WINDOW,Float64}}()

const freed = Vector{Block}()
const freed_lock = SpinLock()

# scan every pool and manage the usage history
#
# returns a boolean indicating whether any pool is active (this can be a false negative)
function scan()
  GC.gc(false) # quick, incremental collection

  active = false

  @lock pool_lock begin
    @inbounds for pid in 1:length(pool_history)
      nused = length(pools_used[pid])
      navail = length(pools_avail[pid])
      history = pool_history[pid]

      if nused+navail > 0
        usage = pool_usage[pid]
        current_usage = nused / (nused + navail)

        # shift the history window with the recorded usage
        history = pool_history[pid]
        pool_history[pid] = (Base.tail(pool_history[pid])..., usage)

        # reset the usage with the current one
        pool_usage[pid] = current_usage

        if usage != current_usage
          active = true
        end
      else
        pool_usage[pid] = 1
        pool_history[pid] = initial_usage
      end
    end
  end

  active
end

# reclaim unused buffers
function reclaim(target_bytes::Int=typemax(Int); full::Bool=true)
  repopulate()

  @lock pool_lock begin
    # find inactive buffers
    @pool_timeit "scan" begin
      pools_inactive = Vector{Int}(undef, length(pools_avail)) # pid => buffers that can be freed
      if full
        # consider all currently unused buffers
        for (pid, avail) in enumerate(pools_avail)
          pools_inactive[pid] = length(avail)
        end
      else
        # only consider inactive buffers
        @inbounds for pid in 1:length(pool_usage)
          nused = length(pools_used[pid])
          navail = length(pools_avail[pid])
          recent_usage = (pool_history[pid]..., pool_usage[pid])

          if navail > 0
            # reclaim as much as the usage allows
            reclaimable = floor(Int, (1-maximum(recent_usage))*(nused+navail))
            pools_inactive[pid] = reclaimable
          else
            pools_inactive[pid] = 0
          end
        end
      end
    end

    # reclaim buffers (in reverse, to discard largest buffers first)
    @pool_timeit "reclaim" begin
      freed_bytes = 0
      for pid in reverse(eachindex(pools_inactive))
        bytes = poolsize(pid)
        avail = pools_avail[pid]

        bufcount = pools_inactive[pid]
        @assert bufcount <= length(avail)
        for i in 1:bufcount
          block = pop!(avail)

          actual_free(block)

          freed_bytes += bytes
          if freed_bytes >= target_bytes
            return freed_bytes
          end
        end
      end
      return freed_bytes
    end
  end
end

# repopulate the "available" pools from the list of freed blocks
function repopulate()
  blocks = @lock freed_lock begin
    isempty(freed) && return
    blocks = Set(freed)
    empty!(freed)
    blocks
  end

  @lock pool_lock begin
    for block in blocks
      pid = poolidx(sizeof(block))

      @inbounds used = pools_used[pid]
      @inbounds avail = pools_avail[pid]

      # mark the buffer as available
      delete!(used, block)
      push!(avail, block)

      # update pool usage
      current_usage = length(used) / (length(used) + length(avail))
      pool_usage[pid] = max(pool_usage[pid], current_usage)
    end
  end

  return
end

function pool_alloc(bytes, pid=-1)
  block = nothing

  # NOTE: checking the pool is really fast, and not included in the timings
  @lock pool_lock begin
    if pid != -1 && !isempty(pools_avail[pid])
      block = pop!(pools_avail[pid])
    end
  end

  if block === nothing
    @pool_timeit "1. try alloc" begin
      block = actual_alloc(bytes)
    end
  end

  if block === nothing
    @pool_timeit "2a. gc (incremental)" begin
      GC.gc(false)
    end

    @pool_timeit "2b. repopulate" begin
      repopulate()
    end

    @lock pool_lock begin
      if pid != -1 && !isempty(pools_avail[pid])
        block = pop!(pools_avail[pid])
      end
    end
  end

  # TODO: we could return a larger allocation here, but that increases memory pressure and
  #       would require proper block splitting + compaction to be any efficient.

  if block === nothing
    @pool_timeit "3. reclaim unused" begin
      reclaim(bytes)
    end

    @pool_timeit "4. try alloc" begin
      block = actual_alloc(bytes)
    end
  end

  if block === nothing
    @pool_timeit "5a. gc (full)" begin
      GC.gc(true)
    end

    @pool_timeit "5b. repopulate" begin
      repopulate()
    end

    @lock pool_lock begin
      if pid != -1 && !isempty(pools_avail[pid])
        block = pop!(pools_avail[pid])
      end
    end
  end

  if block === nothing
    @pool_timeit "6. reclaim unused" begin
      reclaim(bytes)
    end

    @pool_timeit "7. try alloc" begin
      block = actual_alloc(bytes)
    end
  end

  if block === nothing
    @pool_timeit "8. reclaim everything" begin
      reclaim()
    end

    @pool_timeit "9. try alloc" begin
      block = actual_alloc(bytes)
    end
  end

  if block !== nothing && pid != -1
    @inbounds used = pools_used[pid]
    @inbounds avail = pools_avail[pid]

    # mark the buffer as used
    push!(used, block)

    # update pool usage
    current_usage = length(used) / (length(avail) + length(used))
    pool_usage[pid] = max(pool_usage[pid], current_usage)
  end

  return block
end

function pool_free(block)
    # we don't do any work here to reduce pressure on the GC (spending time in finalizers)
    # and to simplify locking (preventing concurrent access during GC interventions)
  @lock freed_lock begin
    push!(freed, block)
  end
end

## interface

const allocated = Dict{CuPtr{Nothing},Block}()

function init()
  create_pools(30) # up to 512 MiB

  if haskey(ENV, "CUARRAYS_MANAGED_POOL")
    Base.depwarn("The CUARRAYS_MANAGED_POOL environment flag is deprecated, please use JULIA_CUDA_MEMORY_POOL instead.", :__init_memory__)
    ENV["JULIA_CUDA_MEMORY_POOL_MANAGED"] = ENV["CUARRAYS_MANAGED_POOL"]
  end

  managed = parse(Bool, get(ENV, "JULIA_CUDA_MEMORY_POOL_MANAGED", "true"))
  if managed
    delay = MIN_DELAY
    @async begin
      while true
        @pool_timeit "background task" begin
          if scan()
            delay = MIN_DELAY
          else
            delay = min(delay*2, MAX_DELAY)
          end

          reclaim(full=false)
        end

        sleep(delay)
      end
    end
  end
end

function alloc(bytes)
  # only manage small allocations in the pool
  block = if bytes <= MAX_POOL
    pid = poolidx(bytes)
    create_pools(pid)
    alloc_bytes = poolsize(pid)
    pool_alloc(alloc_bytes, pid)
  else
    pool_alloc(bytes)
  end

  if block !== nothing
    ptr = pointer(block)
    @lock pool_lock begin
      allocated[ptr] = block
    end
    return ptr
  else
    return nothing
  end
end

function free(ptr)
  block = @lock pool_lock begin
    block = allocated[ptr]
    delete!(allocated, ptr)
    block
  end
  bytes = sizeof(block)

  # was this a pooled buffer?
  if bytes <= MAX_POOL
    pid = poolidx(bytes)
    @assert pid <= length(pools_used)
    @assert pid == poolidx(sizeof(block))
    pool_free(block)
  else
    actual_free(block)
  end

  return
end

function used_memory()
  sz = 0
  @lock pool_lock for (pid, pl) in enumerate(pools_used)
    bytes = poolsize(pid)
    sz += bytes * length(pl)
  end

  return sz
end

function cached_memory()
  sz = 0
  @lock pool_lock for (pid, pl) in enumerate(pools_avail)
    bytes = poolsize(pid)
    sz += bytes * length(pl)
  end

  return sz
end

dump() = return

end


text:
using FileTools
using Test
using DataFrames
using Dates
using DelimitedFiles: readdlm

# List of test files. Run the test from FileTools.jl folder
tests = ["asciifile","atmacsfile","baytapfile","campbell","dwdfile","dygraphs",
		 "eopfile","et0calc","ggpfile","gpcpdfile","gravityeffectfile",
		 "hydrusfile","otherfile","tidefile","tsffile"];
# Run all tests in the list
for i in tests
	include("$(i).jl")
end


text:
@testset verbose=true "Random" begin
    using Random
    Random.seed!(1234)
    rng = MersenneTwister(1234)

    @testset "$BT" for BT in FloatTypes
        for T in [Quaternion{BT}, Rotor{BT}, QuatVec{BT}]
            @test randn(T) isa T
            @test randn(T, 4) isa Vector{T}
            @test randn(T, 2, 3) isa Matrix{T}
            @test randn(T, 2, 3, 4) isa Array{T}
            # @test randn(T, ()) isa T
            @test randn(T, (4,)) isa Vector{T}
            @test randn(T, (2, 3)) isa Matrix{T}
            @test randn(T, (2, 3, 4)) isa Array{T}

            @test randn(rng, T) isa T
            @test randn(rng, T, 4) isa Vector{T}
            @test randn(rng, T, 2, 3) isa Matrix{T}
            @test randn(rng, T, 2, 3, 4) isa Array{T}
            # @test randn(rng, T, ()) isa T
            @test randn(rng, T, (4,)) isa Vector{T}
            @test randn(rng, T, (2, 3)) isa Matrix{T}
            @test randn(rng, T, (2, 3, 4)) isa Array{T}
        end

        N = (BT === BigFloat ? 10_000 : 1_000_000)
        error = maximum(abs, 1.0 .- abs2.(randn(Rotor{BT}, N)))
        @test error ≈ zero(BT) atol=4eps(BT)
    end
end


text:
import Pkg; Pkg.activate(@__DIR__); Pkg.instantiate();
#Write Altro folder instead of @__DIR__

using TrajectoryOptimization
using Plots, LinearAlgebra

#number individual systems
agent_number=6;

#defining product system dynamics
function dynamics!(ẋ::AbstractVector{T},x::AbstractVector{T},u::AbstractVector{T}) where T
	for k=1:agent_number
		ẋ[3*k-2] = u[2*k-1]*cos(x[3*k]);
		ẋ[3*k-1] = u[2*k-1]*sin(x[3*k]);
	    ẋ[3*k] = u[2*k];
	end
end
#Total number of states
n=agent_number*3
#Total number of control inputs
m=agent_number*2

#defining descrete model
model = Model(dynamics!,n,m)
model_d = rk4(model)
n = model.n;
m = model.m;

T = Float64;

#inital point for product system
x0 = [
6.0;1.2;0.0;
3.0;1.2;0.0;
0.0;1.2;0.0;
8;3.9;-0.64;
6;5.4;-0.64;
4.0;6.9;-0.64
]
#goal point for product system
xf = [
17.0;1.2;0.0;
17.0;1.2;0.0;
17.0;1.2;0.0;
17.0;1.2;0.0;
17.0;1.2;0.0;
17.0;1.2;0.0
]

goal = goal_constraint(xf);

#these 2 parameters are very important
#number of points in trajectory(number of time samples)
N = 110
#sampling time
dt = 0.1
#total time
tf = (N-1)*dt

#LQR Cost Matrices (same as default)
Q = 0.01*Diagonal(I,n)
Qf = 1000.0*Diagonal(I,n)
R = 0.01*Diagonal(I,m)
obj = LQRObjective(Q,R,Qf,xf,N)

#defining Constraint for each point in time
constraints = Constraints(N)

#creating upperbound and lower bound for control input with bound constraint
u_min = Array{Float64}(undef,m)
u_max = Array{Float64}(undef,m)
u_c=2;
for k=1:m
	u_min[k]=-u_c;
	u_max[k]=u_c;
end
bnd = BoundConstraint(n,m,u_min=u_min,u_max=u_max)
for k = 1:N-1
    constraints[k] += bnd
 end

#collision constraint

function collision_constraint(x)
	ϵ=1.22;
	M=1000;
	temp=ϵ+0.1
	for i=1:agent_number
		for j=i+1:agent_number
			if((x[3*i-2]<15) &&  (x[3*j-2]<15))
				temp=max((x[3*i-2]-x[3*j-2])^2 , (x[3*i-1]-x[3*j-1])^2);
			else
				continue;
			end
			M=min(temp,M);
		end
	end
	return -(M - ϵ^2)
end

#These 3 constraints are for avoiding the road sides

function constraint1(x)
	M=-10000
	temp=100
	for i=1:agent_number
		if( (3*x[3*i-2]+4*x[3*i-1] <=9.15*4) && (x[3*i-1]>=1.75) )
			temp=min(((3*x[3*i-2]+4*x[3*i-1]-9.15*4)/5)^2,(x[3*i-1]-1.75)^2);
		else
			 temp=-min(((3*x[3*i-2]+4*x[3*i-1]-9.15*4)/5)^2,(x[3*i-1]-1.75)^2) ;
		 end
		 M=max(M,temp);
    end
    return M;
end

function constraint2(x)
	M=-10000
	temp=100
	for i=1:agent_number
		if( (3*x[3*i-2]+4*x[3*i-1] >=10.65*4) && (x[3*i-1]>=1.75))
			temp=min(((3*x[3*i-2]+4*x[3*i-1]-10.65*4)/5)^2 , (x[3*i-1]-1.75)^2)
		else
			temp= -min(((3*x[3*i-2]+4*x[3*i-1]-10.65*4)/5)^2,(x[3*i-1]-1.75)^2);
		 end
		 M=max(M,temp);
    end
    return M
end

function constraint3(x)
	M=-10000
	temp=100
	for i=1:agent_number
		if( (x[3*i-1]<=0.65) )
			temp=(x[3*i-1]-0.65)^2
		else
			 temp= -(x[3*i-1]-0.65)^2;
		 end
		 M=max(M,temp);
    end
    return M;
end

function AllConstraints(c,x,u)
    c[1]=collision_constraint(x);
    c[2]=constraint1(x);
	c[3]=constraint2(x);
	c[4]=constraint3(x);
	return nothing
end
obs= Constraint{Inequality}(AllConstraints,n,m,4,:obs)

for k = 1:N-1
     constraints[k] += obs
end

prob = Problem(model_d,obj,constraints=constraints,x0=x0,N=N,dt=dt,xf=xf)

# options for ALTRO Solver
max_con_viol = 1.0e-8
verbose=false

opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,
    iterations=300)

opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,
    opts_uncon=opts_ilqr,
    iterations=40,
    cost_tolerance=1.0e-5,
    cost_tolerance_intermediate=1.0e-4,
    constraint_tolerance=max_con_viol,
    penalty_scaling=10.,
    penalty_initial=1.)

opts_pn = ProjectedNewtonSolverOptions{T}(verbose=verbose,
    feasibility_tolerance=max_con_viol,
    solve_type=:feasible)

opts_altro = ALTROSolverOptions{T}(verbose=verbose,
    opts_al=opts_al,
    R_inf=1.0e-8,
    resolve_feasible_problem=false,
    opts_pn=opts_pn,
    projected_newton=true,
    projected_newton_tolerance=1.0e-4);

	solver = AugmentedLagrangianSolver(prob)

    @time begin
    solution=solve!(prob, opts_altro); # solve with ALTRO
    end

#writing the solutions

	io = open("nom_tr.txt", "w")
    for i in prob.U
		for k=1:m-1
        	print(io,i[k])
        	print(io," ")
		end
		println(io,i[m])
    end
    close(io)

#Parsing the nominal trajectory to different arrays

	x=Array{Array{Float64,1}}(undef,agent_number)
	z=Array{Array{Float64,1}}(undef,agent_number)
	u1=Array{Array{Float64,1}}(undef,agent_number)
	u2=Array{Array{Float64,1}}(undef,agent_number)
	for i=1:agent_number
    	x[i] = [prob.X[k][3*i-2] for k = 1:N]
    	z[i] = [prob.X[k][3*i-1] for k = 1:N]
		u1[i]=[prob.U[k][2*i-1] for k = 1:N-1]
		u2[i]=[prob.U[k][2*i] for k = 1:N-1]
	end

#Plot and Animate

	a=plot(xlim=(-1,20),ylim=(-1,20))

	plot!(x[1],z[1],xlabel="x1",ylabel="x2",label="1",legend=:topleft,color=:blue,width=2,ratio=:equal,title="Multi-Agent Example Trajectories",alpha=0.5)
	plot!(x[2],z[2],label="2",legend=:topleft,width=2,ratio=:equal,alpha=0.5)
	plot!(x[3],z[3],label="3",legend=:topleft,width=2,ratio=:equal,alpha=0.5)
	plot!(x[4],z[4],label="4",legend=:topleft,width=2,ratio=:equal,alpha=0.5)
	plot!(x[5],z[5],label="5",legend=:topleft,width=2,ratio=:equal,alpha=0.5)
	plot!(x[6],z[6],label="6",legend=:topleft,width=2,ratio=:equal,alpha=0.5)

	plot!((x[1][1],z[1][1]),marker=:circle,color=:red,label="")
    plot!((x[1][end],z[1][end]),marker=:utriangle,color=:green,label="")

	plot!((x[2][1],z[2][1]),marker=:circle,color=:red,label="")
    plot!((x[2][end],z[2][end]),marker=:utriangle,color=:green,label="")

	plot!((x[3][1],z[3][1]),marker=:circle,color=:red,label="")
    plot!((x[3][end],z[3][end]),marker=:utriangle,color=:green,label="")

	plot!((x[4][1],z[4][1]),marker=:circle,color=:red,label="")
    plot!((x[4][end],z[4][end]),marker=:utriangle,color=:green,label="")

	plot!((x[5][1],z[5][1]),marker=:circle,color=:red,label="")
    plot!((x[5][end],z[5][end]),marker=:utriangle,color=:green,label="")

	plot!((x[6][1],z[6][1]),marker=:circle,color=:red,label="")
    plot!((x[6][end],z[6][end]),marker=:utriangle,color=:green,label="")

	plot!([-1.0,18.0],[0.0,0.0],color=:black,label="")
	plot!([-1.0,8.0],[2.4,2.4],color=:black,label="")
	plot!([12.0,18.0],[2.4,2.4],color=:black,label="")
	plot!([-1.0,8.0],[9.15,2.4],color=:black,label="")
	plot!([-1.0,12.0],[12.15,2.4],color=:black,label="")

	plot!([-1.0,18.0],[0.6,0.6],color=:red,label="")
	plot!([-1.0,9.8],[1.8,1.8],color=:red,label="")
	plot!([11.8,18.0],[1.8,1.8],color=:red,label="")
	plot!([-1.0,9.8],[9.9,1.8],color=:red,label="")
	plot!([-1.0,11.8],[11.4,1.8],color=:red,label="")

	anim = @animate for i in 1:N
		plot((x[1][i],z[1][i]),marker=:circle,xlim=(0,15),ylim=(0,15))
		plot!((x[2][i],z[2][i]),marker=:circle)
		plot!((x[3][i],z[3][i]),marker=:circle)
		plot!((x[4][i],z[4][i]),marker=:circle)
		plot!((x[5][i],z[5][i]),marker=:circle)
		plot!((x[6][i],z[6][i]),marker=:circle)

		plot!([0.0,18.0],[0.0,0.0],color=:black,label="")
		plot!([0.0,8.0],[2.4,2.4],color=:black,label="")
		plot!([12.0,18.0],[2.4,2.4],color=:black,label="")
		plot!([0.0,8.0],[8.4,2.4],color=:black,label="")
		plot!([0.0,12.0],[11.4,2.4],color=:black,label="")

		plot!([0.0,18.0],[0.6,0.6],color=:red,label="")
		plot!([0.0,9.8],[1.8,1.8],color=:red,label="")
		plot!([11.8,18.0],[1.8,1.8],color=:red,label="")
		plot!([0.0,9.8],[9.15,1.8],color=:red,label="")
		plot!([0.0,11.8],[10.65,1.8],color=:red,label="")
	end

max_stars_repo_name: MehrdadZareian/GMARA

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.
mutates(::ANode{:store!}) = true

arg(x::ANode{:axpy!}) = arg(x.input[1])
arg(x::ANode{:store!}) = arg(x.input[1])

function store!(x::ANode, y::ANode)
  @assert size(x) == size(y)
  return ANode(:store!, (x, y), size(x))
end

function do_forward!(n::ANode{:store!}, y, x)
  @assert length(y) == length(x)
  copy!(y, x)
end

@cuda function do_forward!(n::ANode{:store!}, y::CudaArray, x)
  @assert length(y) == length(x)
  copy!(y, x; stream=STREAM[1].handle)
end

mutates(::ANode{:axpy!}) = true

function axpy!(y::ANode, x::ANode, alpha::Real=1)
  @assert size(x) == size(y)
  return store!(y, y + alpha*x)
end

function do_forward!(n::ANode{:axpy!}, y, x)
  @assert length(y) == length(x)
  BLAS.axpy!(length(y), convert(eltype(y), n.arg), x, 1, y, 1)
end


text:
function puzzle06(path::String=joinpath(@__DIR__, "..", "data/06.txt"))
    data = readdlm(path, Int)
    n = size(data, 2)
    configurations = Dict()
    notfound = true
    iter = 1

    while notfound
        val, pos = findmax(data)
        data[pos] = 0
        for offset in 1:val
            i = offset + pos[2]
            r = rem(i, n)
            i = r==0 ? n : r
            data[i] += 1
        end

        dhash = hash(data)

        if haskey(configurations, dhash)
            notfound = false
            partone = iter
            parttwo = iter - configurations[dhash]
            return [partone, parttwo]
        end

        configurations[dhash] = iter
        iter += 1

    end

end


text:
using Test

using MARY_fRG

include("test_basics/test_points.jl")

include("test_basics/test_lines.jl")

include("test_basics/test_triangles.jl")

include("test_basics/test_rectangle.jl")

include("test_basics/test_hexagon.jl")

#绘图的测试
#include("test_drawers/test_polygon.jl")

include("test_triangulated/test_square.jl")

include("test_triangulated/test_hexagon.jl")

#refine
#include("test_triangulated/test_refine.jl")

include("test_fermi/test_brlu.jl")

include("test_fermi/test_latt.jl")

include("test_flowequ/test_bubble.jl")

include("test_flowequ/test_dl_ec.jl")


text:
using Documenter
using ResumableFunctions

makedocs(
  format = :html,
  sitename = "ResumableFunctions",
  authors = "Ben Lauwens",
  pages = [
    "Home" => "index.md",
    "Manual" => "manual.md",
    "Internals" => "internals.md",
    "Library" => "library.md"
  ]
)

deploydocs(
  julia = "1.0",
  repo = "github.com/BenLauwens/ResumableFunctions.jl",
  target = "build",
  deps   = nothing,
  make   = nothing
)


text:
""" Parse relation expressions in Julia syntax into undirected wiring diagrams.
"""
module RelationalPrograms
export RelationDiagram, UntypedRelationDiagram, TypedRelationDiagram,
  @relation, parse_relation_diagram

using MLStyle: @match

using ...CategoricalAlgebra.CSets, ...Present
using ...WiringDiagrams.UndirectedWiringDiagrams
using ...WiringDiagrams.UndirectedWiringDiagrams: TheoryUWD, TheoryTypedUWD

# Data types
############

""" Abstract type for UWDs created by [`@relation`](@ref) macro.
"""
@abstract_acset_type RelationDiagram <: AbstractUWD
@abstract_acset_type _UntypedRelationDiagram <: RelationDiagram
@abstract_acset_type _TypedRelationDiagram <: RelationDiagram

""" Untyped UWD created by [`@relation`](@ref) macro.
"""
const UntypedRelationDiagram{Name,VarName} =
  _UntypedRelationDiagram{S, Tuple{Name,VarName}} where S

""" Typed UWD created by [`@relation`](@ref) macro.
"""
const TypedRelationDiagram{T,Name,VarName} =
  _TypedRelationDiagram{S, Tuple{T,Name,VarName}} where S

@present TheoryRelationDiagram <: TheoryUWD begin
  (Name, VarName)::AttrType
  name::Attr(Box, Name)
  variable::Attr(Junction, VarName)
end

@present TheoryTypedRelationDiagram <: TheoryTypedUWD begin
  (Name, VarName)::AttrType
  name::Attr(Box, Name)
  variable::Attr(Junction, VarName)
end

@acset_type UntypedUnnamedRelationDiagram(TheoryRelationDiagram,
  index=[:box, :junction, :outer_junction], unique_index=[:variable]) <: _UntypedRelationDiagram
@acset_type TypedUnnamedRelationDiagram(TheoryTypedRelationDiagram,
  index=[:box, :junction, :outer_junction], unique_index=[:variable]) <: _TypedRelationDiagram

@present TheoryNamedRelationDiagram <: TheoryRelationDiagram begin
  port_name::Attr(Port, Name)
  outer_port_name::Attr(OuterPort, Name)
end

@present TheoryTypedNamedRelationDiagram <: TheoryTypedRelationDiagram begin
  port_name::Attr(Port, Name)
  outer_port_name::Attr(OuterPort, Name)
end

@acset_type UntypedNamedRelationDiagram(TheoryNamedRelationDiagram,
  index=[:box, :junction, :outer_junction], unique_index=[:variable]) <: _UntypedRelationDiagram
@acset_type TypedNamedRelationDiagram(TheoryTypedNamedRelationDiagram,
  index=[:box, :junction, :outer_junction], unique_index=[:variable]) <: _TypedRelationDiagram

function UntypedRelationDiagram{Name,VarName}(
    nports::Int; port_names=nothing) where {Name,VarName}
  if isnothing(port_names)
    return UntypedUnnamedRelationDiagram{Name,VarName}(nports)
  end
  d = UntypedNamedRelationDiagram{Name,VarName}(nports)
  set_subpart!(d, :outer_port_name, port_names)
  return d
end

function TypedRelationDiagram{T,Name,VarName}(
    ports::AbstractVector; port_names=nothing) where {T,Name,VarName}
  if isnothing(port_names)
    return TypedUnnamedRelationDiagram{T,Name,VarName}(ports)
  end
  d = TypedNamedRelationDiagram{T,Name,VarName}(ports)
  set_subpart!(d, :outer_port_name, port_names)
  return d
end

RelationDiagram(nports::Int; kw...) =
  UntypedRelationDiagram{Symbol,Symbol}(nports; kw...)
RelationDiagram(ports::AbstractVector{T}; kw...) where T =
  TypedRelationDiagram{T,Symbol,Symbol}(ports; kw...)

# Relation macro
################

""" Construct an undirected wiring diagram using relation notation.

Unlike the `@program` macro for directed wiring diagrams, this macro departs
significantly from the usual semantics of the Julia programming language.
Function calls with ``n`` arguments are now interpreted as assertions that an
``n``-ary relation holds at a particular point. For example, the composite of
binary relations ``R ⊆ X × Y`` and ``S ⊆ Y × Z`` can be represented as an
undirected wiring diagram by the macro call

```julia
@relation (x,z) where (x::X, y::Y, z::Z) begin
  R(x,y)
  S(y,z)
end
```

In general, the context in the `where` clause defines the set of junctions in
the diagram and variable sharing defines the wiring of ports to junctions. If
the `where` clause is omitted, the set of junctions is inferred from the
variables used in the macro call.

The ports and junctions of the diagram may be typed or untyped, and the ports
may be named or unnamed. Thus, four possible types of undirected wiring diagrams
may be returned, with the type determined by the form of relation header:

1. Untyped, unnamed: `@relation (x,z) where (x,y,z) ...`
2. Typed, unnamed: `@relation (x,z) where (x::X, y::Y, z::Z) ...`
3. Untyped, named: `@relation (out1=x, out2=z) where (x,y,z) ...`
4. Typed, named: `@relation (out=1, out2=z) where (x::X, y::Y, z::Z) ...`

All four types of diagram are subtypes of [`RelationDiagram`](@ref).
"""
macro relation(exprs...)
  :(parse_relation_diagram($((QuoteNode(expr) for expr in exprs)...)))
end

""" Parse an undirected wiring diagram from a relation expression.

For more information, see the corresponding macro [`@relation`](@ref).
"""
function parse_relation_diagram(expr::Expr)
  @match expr begin
    Expr(:function, head, body) => parse_relation_diagram(head, body)
    Expr(:->, head, body) => parse_relation_diagram(head, body)
    _ => error("Not a function or lambda expression")
  end
end

function parse_relation_diagram(head::Expr, body::Expr)
  # Parse variables and their types from context.
  outer_expr, all_vars, all_types = @match head begin
    Expr(:where, expr, context) => (expr, parse_relation_context(context)...)
    _ => (head, nothing, nothing)
  end
  var_types = if isnothing(all_types) # Untyped case.
    vars -> length(vars)
  else # Typed case.
    var_type_map = Dict{Symbol,Symbol}(zip(all_vars, all_types))
    vars -> getindex.(Ref(var_type_map), vars)
  end

  # Create wiring diagram and add outer ports and junctions.
  _, outer_port_names, outer_vars = parse_relation_call(outer_expr)
  isnothing(all_vars) || outer_vars ⊆ all_vars ||
    error("One of variables $outer_vars is not declared in context $all_vars")
  d = RelationDiagram(var_types(outer_vars), port_names=outer_port_names)
  if isnothing(all_vars)
    new_vars = unique(outer_vars)
    add_junctions!(d, var_types(new_vars), variable=new_vars)
  else
    add_junctions!(d, var_types(all_vars), variable=all_vars)
  end
  set_junction!(d, ports(d, outer=true),
                incident(d, outer_vars, :variable), outer=true)

  # Add box to diagram for each relation call.
  body = Base.remove_linenums!(body)
  for expr in body.args
    name, port_names, vars = parse_relation_call(expr)
    isnothing(all_vars) || vars ⊆ all_vars ||
      error("One of variables $vars is not declared in context $all_vars")
    box = add_box!(d, var_types(vars), name=name)
    if !isnothing(port_names)
      set_subpart!(d, ports(d, box), :port_name, port_names)
    end
    if isnothing(all_vars)
      new_vars = setdiff(unique(vars), d[:variable])
      add_junctions!(d, var_types(new_vars), variable=new_vars)
    end
    set_junction!(d, ports(d, box), incident(d, vars, :variable))
  end
  return d
end

function parse_relation_context(context)
  terms = @match context begin
    Expr(:tuple) => return (Symbol[], nothing)
    Expr(:tuple, terms...) => terms
    _ => error("Invalid syntax in relation context $context")
  end
  vars = map(terms) do term
    @match term begin
      Expr(:(::), var::Symbol, type::Symbol) => (var => type)
      var::Symbol => var
      _ => error("Invalid syntax in term $expr of context")
    end
  end
  if vars isa AbstractVector{Symbol}
    (vars, nothing)
  elseif vars isa AbstractVector{Pair{Symbol,Symbol}}
    (first.(vars), last.(vars))
  else
    error("Context $context mixes typed and untyped variables")
  end
end

function parse_relation_call(call)
  @match call begin
    Expr(:call, name::Symbol, Expr(:parameters, args)) =>
      (name, parse_relation_kw_args(args)...)
    Expr(:call, name::Symbol) => (name, nothing, Symbol[])
    Expr(:call, name::Symbol, args...) =>
      (name, parse_relation_inferred_args(args)...)

    Expr(:tuple, Expr(:parameters, args...)) =>
      (nothing, parse_relation_kw_args(args)...)
    Expr(:tuple) => (nothing, nothing, Symbol[])
    Expr(:tuple, args...) => (nothing, parse_relation_inferred_args(args)...)
    Expr(:(=), args...) => (nothing, parse_relation_inferred_args([call])...)

    _ => error("Invalid syntax in relation $call")
  end
end

function parse_relation_kw_args(args)
  args = map(args) do arg
    @match arg begin
      Expr(:kw, name::Symbol, var::Symbol) => (name => var)
      _ => error("Expected name as keyword argument")
    end
  end
  (first.(args), last.(args))
end

function parse_relation_inferred_args(args)
  @assert !isempty(args) # Need at least one argument to infer named/unnamed.
  args = map(args) do arg
    @match arg begin
      Expr(:kw, name::Symbol, var::Symbol) => (name => var)
      Expr(:(=), name::Symbol, var::Symbol) => (name => var)
      var::Symbol => var
      _ => error("Expected name as positional or keyword argument")
    end
  end
  if args isa AbstractVector{Symbol}
    (nothing, args)
  elseif args isa AbstractVector{Pair{Symbol,Symbol}}
    (first.(args), last.(args))
  else
    error("Relation mixes named and unnamed arguments $args")
  end
end

end


text:
function get_circle(obj::Edge)
    return jcall(obj, "getCircle", Circle, ())
end

function get_end(obj::Edge)
    return jcall(obj, "getEnd", Vertex, ())
end

function get_length(obj::Edge)
    return jcall(obj, "getLength", jdouble, ())
end

function get_point_at(obj::Edge, arg0::jdouble)
    return jcall(obj, "getPointAt", Vector3D, (jdouble,), arg0)
end

function get_start(obj::Edge)
    return jcall(obj, "getStart", Vertex, ())
end