text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
config = MOIT.TestConfig()

@testset "Scalar slack" begin
    MOI.empty!(mock)
    bridged_mock = MOIB.Constraint.ScalarSlack{Float64}(mock)

    x = MOI.add_variable(bridged_mock)
    y = MOI.add_variable(bridged_mock)
    f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}.([1.0, 2.0], [x, y]), 0.0)
    ci = MOI.add_constraint(bridged_mock, f, MOI.GreaterThan(0.0))

    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ f
    newf = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}.([2.0, 1.0], [x, y]), 0.0)
    MOI.set(bridged_mock, MOI.ConstraintFunction(), ci, newf)
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ newf
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == MOI.GreaterThan(0.0)
    MOI.set(bridged_mock, MOI.ConstraintSet(), ci, MOI.GreaterThan(1.0))
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == MOI.GreaterThan(1.0)
    MOI.modify(bridged_mock, ci, MOI.ScalarConstantChange{Float64}(1.0))
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}.([2.0, 1.0], [x, y]), 1.0)

    @testset "$attr" for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        MOI.set(bridged_mock, attr, ci, 2.0)
        @test MOI.get(bridged_mock, attr, ci) == 2.0
        bridge = MOIB.bridge(bridged_mock, ci)
        if attr isa MOI.ConstraintPrimalStart
            @test MOI.get(mock, MOI.VariablePrimalStart(), bridge.slack) == 2.0
            @test MOI.get(mock, attr, bridge.equality) == 0.0
        else
            @test MOI.get(mock, attr, bridge.equality) == 2.0
        end
    end

    test_delete_bridge(bridged_mock, ci, 2, (
        (MOI.SingleVariable, MOI.GreaterThan{Float64}, 0),
        (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}, 0)
    ))

    MOIT.basic_constraint_tests(bridged_mock, config,
                                include=[(F, S) for
                                F in [MOI.ScalarAffineFunction{Float64},
                                      MOI.ScalarQuadraticFunction{Float64}],
                                S in [MOI.GreaterThan{Float64},
                                      MOI.LessThan{Float64}]])

    # There are extra variables due to the bridge
    MOIU.set_mock_optimize!(mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1, 0, 1], con_basis =
            [(MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})  => [MOI.NONBASIC],
             (MOI.SingleVariable, MOI.GreaterThan{Float64})             => [MOI.BASIC, MOI.NONBASIC],
             (MOI.SingleVariable, MOI.LessThan{Float64})                => [MOI.NONBASIC]]))
    MOIT.linear2test(bridged_mock, MOIT.TestConfig(duals = false, basis = true))
    c1 = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}())
    @test length(c1) == 1
    @test MOI.get(bridged_mock, MOI.ConstraintBasisStatus(), c1[]) == MOI.NONBASIC

    MOIU.set_mock_optimize!(mock,
     (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 1.0, 2.0, 2.0]),
     (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [0.5, 0.5, 1.0, 1.0],
        (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) => [1, 0],
        (MOI.SingleVariable, MOI.GreaterThan{Float64})            => [1],
        (MOI.SingleVariable, MOI.LessThan{Float64})               => [0]))
    MOIT.linear11test(bridged_mock, MOIT.TestConfig(duals = false))

    c1 = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}())
    @test length(c1) == 1
    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), c1[]) ≈ 1.0
    @test MOI.get(bridged_mock, MOI.ConstraintDual(), c1[]) ≈ 1.0
    c2 = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}())
    @test length(c2) == 1
    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), c2[]) ≈ 1.0
    @test MOI.get(bridged_mock, MOI.ConstraintDual(), c2[]) ≈ 0.0

    loc = MOI.get(bridged_mock, MOI.ListOfConstraints())
    @test length(loc) == 2
    @test (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) in loc
    @test (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) in loc
    loc = MOI.get(mock, MOI.ListOfConstraints())
    @test length(loc) == 3
    @test (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) in loc
    @test (MOI.SingleVariable, MOI.LessThan{Float64}) in loc
    @test (MOI.SingleVariable, MOI.GreaterThan{Float64}) in loc

    for T in [Int, Float64], S in [MOI.GreaterThan{T}, MOI.GreaterThan{T}]
        for F in [MOI.ScalarAffineFunction{T}, MOI.ScalarQuadraticFunction{T}]
            @test MOIB.added_constraint_types(MOIB.Constraint.ScalarSlackBridge{T, F, S}) == [(F, MOI.EqualTo{T})]
        end
    end
end

@testset "Vector slack" begin
    MOI.empty!(mock)
    bridged_mock = MOIB.Constraint.VectorSlack{Float64}(mock)

    x = MOI.add_variable(bridged_mock)
    y = MOI.add_variable(bridged_mock)
    f = MOI.VectorAffineFunction(MOI.VectorAffineTerm.(1, MOI.ScalarAffineTerm.([1.0, 2.0], [x, y])), [0.0])
    ci = MOI.add_constraint(bridged_mock, f, MOI.Nonpositives(1))

    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ f
    newf = MOI.VectorAffineFunction(MOI.VectorAffineTerm.(1, MOI.ScalarAffineTerm.([2.0, 1.0], [x, y])), [0.0])
    MOI.set(bridged_mock, MOI.ConstraintFunction(), ci, newf)
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈ newf
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == MOI.Nonpositives(1)
    MOI.modify(bridged_mock, ci, MOI.VectorConstantChange([1.0]))
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), ci) ≈
        MOI.VectorAffineFunction(MOI.VectorAffineTerm.(1, MOI.ScalarAffineTerm.([2.0, 1.0], [x, y])), [1.0])

    @testset "$attr" for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        MOI.set(bridged_mock, attr, ci, [2.0])
        @test MOI.get(bridged_mock, attr, ci) == [2.0]
        bridge = MOIB.bridge(bridged_mock, ci)
        if attr isa MOI.ConstraintPrimalStart
            @test MOI.get(mock, MOI.VariablePrimalStart(), bridge.slack) == [2.0]
            @test MOI.get(mock, attr, bridge.equality) == [0.0]
        else
            @test MOI.get(mock, attr, bridge.equality) == [2.0]
        end
    end

    test_delete_bridge(bridged_mock, ci, 2, (
        (MOI.VectorOfVariables, MOI.Nonpositives, 0),
        (MOI.VectorAffineFunction{Float64}, MOI.Zeros, 0)
    ))

    fp = MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1,2,3], MOI.ScalarAffineTerm.([1.0, 2.0, 3.0], [x, y, y])), [0.0, 0.0, 0.0])
    cp = MOI.add_constraint(bridged_mock, fp, MOI.PowerCone(0.1))
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), cp) == MOI.PowerCone(0.1)
    MOI.set(bridged_mock, MOI.ConstraintSet(), cp, MOI.PowerCone(0.2))
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), cp) == MOI.PowerCone(0.2)

    MOIT.basic_constraint_tests(bridged_mock, config,
                                include=[(F, S) for
                                F in [MOI.VectorAffineFunction{Float64},
                                      MOI.VectorQuadraticFunction{Float64}],
                                S in [MOI.Nonnegatives,
                                      MOI.Nonpositives]
                                    ])

    # There are extra variables due to the bridge
    MOIU.set_mock_optimize!(mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [0, 0, 0, 0]),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [100, 0, 100, 0]),
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [100, -100, 100, -100],
        (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[1.], [1.]],
        (MOI.VectorOfVariables, MOI.Nonnegatives)      => [[1.]],
        (MOI.VectorOfVariables, MOI.Nonpositives)      => [[1.]]))
    MOIT.linear7test(bridged_mock, config)

    c1 = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives}())
    @test length(c1) == 1
    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), c1[]) ≈ [100.0]
    @test MOI.get(bridged_mock, MOI.ConstraintDual(), c1[]) ≈ [1.0]
    c2 = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Nonpositives}())
    @test length(c2) == 1
    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), c2[]) ≈ [-100.0]
    @test MOI.get(bridged_mock, MOI.ConstraintDual(), c2[]) ≈ [1.0]

    loc = MOI.get(bridged_mock, MOI.ListOfConstraints())
    @test length(loc) == 2
    @test (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) in loc
    @test (MOI.VectorAffineFunction{Float64}, MOI.Nonpositives) in loc
    loc = MOI.get(mock, MOI.ListOfConstraints())
    @test length(loc) == 3
    @test (MOI.VectorAffineFunction{Float64}, MOI.Zeros) in loc
    @test (MOI.VectorOfVariables, MOI.Nonnegatives) in loc
    @test (MOI.VectorOfVariables, MOI.Nonpositives) in loc

    for T in [Int, Float64], S in [MOI.Nonnegatives, MOI.Nonpositives]
        for F in [MOI.VectorAffineFunction{T}, MOI.VectorQuadraticFunction{T}]
            @test MOIB.added_constraint_types(MOIB.Constraint.VectorSlackBridge{T, F, S}) == [(F, MOI.Zeros)]
        end
    end
end


text:
# Curves

type NullCurve <: Curve end

type FittingCost <: CostFunction
  # value::Vector{Float64}
  # values::Vector{Float64}
  firstCashFlow::Vector{Int}
  curve::Curve
end

function FittingCost{C <: Curve}(size::Int, curve::C)
  # value = Vector{Float64}()
  # values = Vector{Float64}()
  firstCashFlow = zeros(Int, size)

  return FittingCost(firstCashFlow, curve)
end

# Interpolated Curve methods #
max_date{C <: InterpolatedCurve}(curve::C) = curve.dates[end]

function discount{C <: Curve}(curve::C, t::Float64)
  return discount_impl(curve, t)
end

function discount_impl{C <: InterpolatedCurve}(curve::C, t::Float64)
  calculate!(curve)
  if t <= curve.times[end]
    return QuantLib.Math.value(curve.interp, t)
  end

  # println("outside")
  # println(t)
  # println(curve.times)

  # do flat fwd extrapolation
end

function perform_calculations!{C <: InterpolatedCurve}(curve::C)
  _calculate!(curve.boot, curve)
  return curve
end

function value{C <: CostFunction, T}(cf::C, x::Vector{T})
  ref_date = cf.curve.referenceDate
  dc = cf.curve.dc
  squared_error = 0.0
  n = length(cf.curve.bonds)

  # for (i, bh) in enumerate(cf.curve.bonds)
  for i=1:length(cf.curve.bonds)
    bond = cf.curve.bonds[i].bond
    bond_settlement = get_settlement_date(bond)
    model_price = -accrued_amount(bond, bond_settlement)
    leg = bond.cashflows
    for k = cf.firstCashFlow[i]:length(leg.coupons)
      # @inbounds df = discount_function(cf.curve.fittingMethod, x, year_fraction(dc, ref_date, date(leg.coupons[k])))
      @inbounds model_price += amount(leg.coupons[k]) * discount_function(cf.curve.fittingMethod, x, year_fraction(dc, ref_date, date(leg.coupons[k])))
    end

    # redemption
    # @inbounds model_price += amount(leg.redemption) * discount_function(cf.curve.fittingMethod, x, year_fraction(dc, ref_date, date(leg.redemption)))

    # adjust NPV for forward settlement
    if bond_settlement != ref_date
      model_price /= discount_function(cf.curve.fittingMethod, x, year_fraction(dc, ref_date, bond_settlement))
    end

    market_price = bond.faceAmount
    price_error = model_price - market_price
    @inbounds weighted_error = cf.curve.fittingMethod.commons.weights[i] * price_error
    squared_error += weighted_error * weighted_error
  end

  return squared_error
end

ext: jl


text:
__precompile__()

module GSA_CME

# Write your package code here.
using PolyChaos
using LinearAlgebra
using CSV
using DelimitedFiles
using DataFrames
using Dates
using ColorSchemes
using Plots
using StatsPlots

include("gsaUtilities.jl")

end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/subsolver_config/config/cbc.jl
# Configuration adjustment subroutines for Cbc.
#############################################################################

function set_default_config!(::Val{:cbc}, ext::ExtensionType, m::T) where T <: MOI.AbstractOptimizer

    return nothing
end


text:
VARIABLES_NOT_DEFINED = "variables are not defined"
TIME_VARIABLE_INEXISTENT = "selected time variable does not exists"
PANEL_VARIABLE_INEXISTENT = "selected panel variable does not exists"
INVALID_METHOD = "invalid method, must be :fast (Float32) or :precise (Float64)"
SELECTED_VARIABLES_DOES_NOT_EXISTS = "selected variables does not exists"
SOME_VARIABLES_NOT_FOUND = "some variables were not found"
TIME_ERROR = "time variable has some gaps or repeated values or"
PANEL_ERROR = "panel variable has some gaps"
DATANAMES_REQUIRED = "required varnames in datanames option"


text:
using PyCall
pyimport_conda("geographiclib", "geographiclib", "conda-forge")

max_stars_repo_name: anowacki/TauPy


text:
using Optim: OnceDifferentiable

export solve_model

function OnceDifferentiable(ew::EstimationWrapper, theta::Vector)
    check_finite(theta)

    function f(x::Vector)
        check_finite(x)
        reset!(LocalEstObj(ew))
        return parallel_simloglik!(ew, x, false)
        # return serial_simloglik!(ew, x, false)
    end

    function fg!(g::Vector,x::Vector)
        check_finite(x)
        nll = parallel_simloglik!(ew, x, true)
        # nll = serial_simloglik!(ew, x, true)
        g .= grad(LocalEstObj(ew))
        check_finite(g)
        return nll
    end

    odfg = Optim.OnceDifferentiable(f, fg!, fg!, theta)
    return odfg
end

function update_invhess!(leo::LocalEstObj)
    rnk = rank(hess(leo))
    k = checksquare(hess(leo))
    if rnk == k
        invhess(leo) .= inv(hess(leo))
    else
        @warn "rank of hess(leo) = $rnk < $k"
    end
    return invhess(leo)
end

function invhessian!(ew::EstimationWrapper, theta)
    parallel_simloglik!(ew, theta, true)
    leo = LocalEstObj(ew)
    return update_invhess!(leo)
end

const OptimOpts = Optim.Options(allow_f_increases=true, show_trace=true, time_limit=30)

# see https://github.com/JuliaNLSolvers/LineSearches.jl/tree/master/src
bfgs(ew) =  Optim.BFGS(;
    # linesearch = Optim.BackTracking(order=3),
    linesearch = Optim.MoreThuente(),
    initial_invH = x -> invhessian!(ew, x)
)

nelder(ew) = Optim.NelderMead()

function solve_model(ew, theta; OptimOpts=OptimOpts, OptimMethod=bfgs(ew))
    check_finite(theta)
    leo = LocalEstObj(ew)
    theta0(leo) .= theta
    odfg  = OnceDifferentiable(ew, theta)
    resetcount!()
    res = optimize(odfg, theta, OptimMethod, OptimOpts)
    return res
end


text:
##Clustering DSM

function Clustering(DSM,label)

    powcc = 1;
    powbid = 0;
    powdep = 1;
    rand_accept = 30;
    rand_bid = 30;
    times = 5;
    stable_limit = 2;

    original_DSM = copy(DSM);
    original_label = copy(label);
    (DSMsize,) = size(original_DSM);

    Clustermatrix = zeros(Int64,DSMsize,DSMsize) + Diagonal(ones(Int64,DSMsize,DSMsize)) ;

    cDSM = copy(original_DSM);

    function ClusterBid(Clustermatrix,cDSM,DSMsize,powdep,powbid,element)

        bestCluster = 0;
        bestBid = -1;

        for i = 1 : DSMsize

            Clusterlist= findall(x-> x == 1 , Clustermatrix[i,:]);
            (Clustersize,) = size(Clusterlist);
            bid = 0;

            for j = 1 : Clustersize

                bid += cDSM[element,Clusterlist[j]];

            end

            bid = bid ^ powdep / Clustersize ^ powbid;

            if bid > bestBid

                 bestBid = bid;
                bestCluster = i;
            end

        end

        return bestCluster

    end

    function TCC(cDSM,DSMsize,Clustermatrix,powcc)

        intraCost=0;
        extraCost=0;

        for i = 1: DSMsize

            for j = 1 : DSMsize

                Clusterofi = findall(x-> x == 1 , Clustermatrix[:,i]);
                Clusterofj = findall(x-> x == 1 , Clustermatrix[:,j]);
                (Clustersizeofi,) = size(Clustermatrix[Clusterofi,:]);

                cost = cDSM[i,j] + cDSM[j,i];

                if Clusterofi == Clusterofj

                    intraCost += cost*Clustersizeofi^powcc;

                else

                    extraCost += cost*DSMsize^powcc;
                end

            end

        end

        totalCost = intraCost + extraCost;

        return totalCost
    end

    TCCost = TCC(cDSM,DSMsize,Clustermatrix,powcc)

    costlist = zeros(DSMsize*times,1);

    for i = 1 : DSMsize * times;

        preClustermatrix = Clustermatrix;
        element = rand(1:DSMsize);
        bestCluster = ClusterBid(Clustermatrix,cDSM,DSMsize,powdep,powbid,element);
        Clustermatrix[:,element] = zeros(DSMsize,1);
        Clustermatrix[bestCluster,element] = 1;
        newCost = TCC(cDSM,DSMsize,Clustermatrix,powcc);
        if newCost <= maximum([TCCost,rand_accept])
            TCCost = newCost;
        else
            Clustermatrix = preClustermatrix;
        end

        costlist[i] = newCost
    end

    orderCluster = [sum(Clustermatrix,dims=2) Clustermatrix]
    orderedCluster = sortslices(orderCluster,dims = 1, rev=true)
    reorderedCluster = orderedCluster[:,2:DSMsize+1]

    Order = [];

    for i = 1 : DSMsize
        Order = [Order ; findall(x-> x == 1 , reorderedCluster[i,:])]
    end

    ## reorder DSM
    I =zeros(DSMsize,DSMsize) + Diagonal(ones(DSMsize,DSMsize))
    P =zeros(DSMsize,DSMsize) + Diagonal(ones(DSMsize,DSMsize))
    Clusteredlabel = copy(original_label)
    copylabel= copy(original_label)

    ## Permutation matrix
    for i = 1 : DSMsize
        P[i,:] = I[Order[i],:]
        Clusteredlabel[i] = copylabel[Order[i]]
    end

    ## order DSM
    ClusteredDSM = P*original_DSM*P'

    return ClusteredDSM, Clusteredlabel;

end


text:
while true
    println("SPAM")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Distances
export pmf_unormalised, cardinality, eachinterseq, eachpath
export get_sample_space, get_normalising_const
export get_true_dist_dict, get_true_dist_vec
export get_entropy
export log_multinomial_ratio, myaddcounts!
export sample_frechet_mean, sample_frechet_mean_mem, sample_frechet_var

function sample_frechet_mean(
    data::Vector{T},
    d::Metric;
    show_progress::Bool=false,
    with_memory::Bool=false
    ) where T

    if with_memory
        return sample_frechet_mean_mem(
            data, d, show_progress=show_progress
        )
    elseif show_progress
        iter = Progress(length(data), 1)
    end
    z_best = Inf
    ind_best = 1
    n = length(data)
    for i in 1:n
        z = 0.0
        j = 1 + (i==1)
        while (z < z_best) & (j ≤ n)
            z += d(data[i], data[j])^2
            j += 1 + (j==(i-1))
        end
        if z < z_best
            z_best = copy(z)
            ind_best = i
        end
        if show_progress
            next!(iter)
        end
    end
    return data[ind_best], ind_best, z_best
end

function sample_frechet_mean_mem(
    data::InteractionSequenceSample{T},
    d::Metric;
    show_progress::Bool=false
    ) where T <:Union{Int, String}
    if show_progress
        iter = Progress(length(data), 1)
    end
    data_unq = unique(data)
    weights = countmap(data)
    z_best = Inf
    ind_best_unq = 1
    for i in eachindex(data_unq)
        z = 0.0
        j = 1 + (i==1)
        S1 = data_unq[i]
        while (z < z_best) & (j ≤ length(data_unq))
            S2 = data_unq[j]
            dist = d(S1,S2)
            z += weights[S2] * dist^2
            j += 1 + (j==(i-1))
        end
        if z < z_best
            z_best = copy(z)
            ind_best_unq = i
        end
        if show_progress
            next!(iter)
        end
    end
    ind_best = findfirst(x->x==data_unq[ind_best_unq], data)
    return data[ind_best], ind_best, z_best
end

function sample_frechet_var(
    data::InteractionSequenceSample{T},
    d::Metric;
    show_progress::Bool=false,
    with_memory::Bool=false
    ) where T <:Union{Int, String}

    S, i, out = sample_frechet_mean(data, d, show_progress=show_progress, with_memory=with_memory)
    return out/length(data)
end

function myaddcounts!(d::Dict{T}, s::T) where {T}
    d[s] = get(d, s, 0) + 1
end

# """
# For two vectors ``X = [x_1, \dots, x_N]`` and ``Y = [y_1, \dots, y_N]`` this calculates the following term
# ```math
#     log \left( \frac{N! w^\prime_1! \cdots w^\prime_{m}!}{M! w_1!\cdots w_{n}!}\right)
# ```
# where ``w_i`` and ``w_j^\prime`` represent the counts of (respectively) ``n`` and ``m`` unique entries of ``X`` and ``Y``.

# This comes into use when looking to sample multisets via marginalising out order information of sequences.
# """
function log_multinomial_ratio(x::AbstractVector, y::AbstractVector)
    if length(x) > length(y)
        x = log_multinomial_ratio(y,x)
        return -x
    end
    # Now we can assume length(x) ≤ length(y)
    z = 0.0
    dx = Dict{eltype(x), Int}()  # To store values seen in x and their counts
    dy = Dict{eltype(y), Int}()  # To store values seen in y and their counts

    # First sort the element counts terms
    @inbounds for i in eachindex(x)
        # @show x_val, y_val, typeof(x_val)
        myaddcounts!(dx, x[i])
        myaddcounts!(dy, y[i])
        z += log(dy[y[i]]) - log(dx[x[i]])
    end
    @inbounds for i=(length(x)+1):(length(y))
        myaddcounts!(dy, y[i])
        z += log(dy[y[i]]) - log(i)
    end
    return z
end

"""
Evaluate the unormalised probability of an interaction seq `x`
"""
function pmf_unormalised(
    model::Union{SIS, SIM},
    x::Vector{Path{Int}})

    return exp(- model.γ * model.dist(x, model.mode))
end

"""
Calculate the sample space cardinality.
"""
function cardinality(
    model::SIS
)::Int
    if (model.K_inner == Inf) | (model.K_outer == Inf)
        return Inf
    else
        V = length(model.V)
        num_paths = V * (V^model.K_inner - 1) / (V - 1)
        return num_paths * (num_paths^model.K_outer - 1) / (num_paths - 1)
    end
end

"""
Calculate the sample space cardinality.
"""
function cardinality(
    model::SPF
)
    if model.K == Inf
        return Inf
    else
        V = length(model.V)
        return Int(V * (V^model.K - 1) / (V - 1))
    end
end

"""
`eachinterseq(V, K, L)`

Returns an iterator over all interaction sequences over the vertex set `V`, with dimension bounds defined by `K` and `L`, specifically
* `V` = vertex set, must be a vector of unique strings or integers;
* `K` = max interaction length, must be an integer;
* `L` = max number of interactions, must be an integer.
"""
function eachinterseq(
    V::Vector{T},
    K_inner::Int,
    K_outer::Int) where T <:Union{Int, String}

    return Base.Iterators.flatten(
        [Base.Iterators.product([eachpath(V, K_inner) for j = 1:k]...) for k=1:K_outer ]
    )
end

function eachpath(V,K::Int)
    return Base.Iterators.flatten(
    [Base.Iterators.product([V for j=1:k]...) for k=1:K]
    )
end

"""
Returns vector with all elements in the sample space.
"""
function get_sample_space(model::SIS)
    z = Vector{Vector{Path{Int}}}()
    for I in eachinterseq(model.V, model.K_inner.u, model.K_outer.u)
        push!(z, [collect(p) for p in I])
    end
    return z
end

function get_sample_space(model::SIM)

    z = Vector{Vector{Path{Int}}}()
    for I in eachinterseq(model.V, model.K_inner.u, model.K_outer.l)
        push!(z, [collect(p) for p in I])
    end

    z = Multiset.(z)
    return unique(z)
end

"""
Returns vector with all elements in the sample space.
"""
function get_sample_space(
    model::SPF
)
    z = Vector{Path}()
    for P in eachpath(model.V, model.K)
        push!(z, collect(P))
    end
    return z

end

"""
Calculate the normalising constant of SIS
"""
function get_normalising_const(
    model::SIS
)

iter = Progress(cardinality(model), 1, "Evaluating normalising constant...")  # Loading bar. Minimum update interval: 1 second
    Z = 0.0
    for S in eachinterseq(model.V, model.K_inner, model.K_outer)
        Z += pmf_unormalised(model, [collect(p) for p in S])
        next!(iter)
    end
    return Z
end

"""
Calculate the true normalising constant.
"""
function get_normalising_const(
    model::SPF
    )

    @assert model.K < Inf "Model must be bounded (K<∞)"
    @assert typeof(model.K)==Int "K must be integer"

    Z = 0.0
    for i=1:model.K
        for P in eachpath(model.V, model.K)
            Z += exp( - model.γ * model.dist(collect(P), model.mode) )
        end
    end
    return Z
end

function get_true_dist_vec(model::SIS; show_progress=true)
    if show_progress
        x = Vector{Float64}()
        iter = Progress(cardinality(model), 1)  # minimum update interval: 1 second
        Z = 0.0 # Normalising constant
        for I in eachinterseq(model.V, model.K_inner.u, model.K_outer.u)
            val = pmf_unormalised(model, [collect(p) for p in I])  # evaluate unormalised probability
            Z += val
            push!(x, val)
            next!(iter)
        end
    else
        x = Vector{Float64}()
        Z = 0.0 # Normalising constant
        for I in eachinterseq(model.V, model.K_inner.u, model.K_outer.u)
            val = pmf_unormalised(model, [collect(p) for p in I])  # evaluate unormalised probability
            Z += val
            push!(x, val)
        end
    end
    return x / Z
end

function get_true_dist_dict(
    model::SIS;
    show_progress=true
)
    d = Dict{Vector{Path{Int}}, Float64}()
    Z = 0.0 # Normalising constant
    prob_val = 0.0
    if show_progress
        iter = Progress(cardinality(model), 1)  # minimum update interval: 1 second
        # val = Vector{Path}()
        for I in eachinterseq(model.V, model.K_inner, model.K_outer)
            val = [collect(p) for p in I]
            # @show val
            prob_val = pmf_unormalised(model, val)  # evaluate unormalised probability
            Z += prob_val
            d[val] = prob_val
            next!(iter)
        end
    else
        for I in eachinterseq(model.V, model.K_inner, model.K_outer)
            val = [collect(p) for p in I]
            prob_val = pmf_unormalised(model, val)  # evaluate unormalised probability
            Z += prob_val
            d[val] = prob_val
        end
    end
    map!(x -> x/Z, values(d)) # Normalise
    return d
end

function get_true_dist_dict(
    model::SIM;
    sample_space::Vector{Multiset{Path{Int}}}=get_sample_space(model)
)
    d = Dict{Multiset{Path{Int}}, Float64}()
    Z = 0.0 # Normalising constant
    prob_val = 0.0

    # sample_space = get_sample_space(model)

    for val in sample_space
        prob_val = pmf_unormalised(model, collect(val))  # collect() turns the multiset val into a vector for passing to pmf_unormalised()
        Z += prob_val
        d[val] = prob_val
    end
    map!(x -> x/Z, values(d)) # Normalise
    return d
end

function get_entropy(
    model::SPF;
    show_progress::Bool=true
    )

    if show_progress
        iter = Progress(
            cardinality(model), # How many iters
            1,  # At which granularity to update loading bar
            "Evaluating entropy....")  # Loading bar. Minimum update interval: 1 second
    end
    d, γ, V, K = (model.dist, model.γ, model.V, model.K) # Aliases
    Z, H = (0.0,0.0)
    for P in eachpath(V, K)
        d_tmp = d(collect(P), model.mode)
        Z += exp(-γ * d_tmp)
        H += - model.γ * d_tmp * exp(-γ * d_tmp)
        if show_progress
            next!(iter)
        end
    end
    return log(Z) - H/Z

end

max_stars_repo_path: src/Metric Models/utils.jl

max_issues_repo_path: src/Metric Models/utils.jl

max_forks_repo_path: src/Metric Models/utils.jl


text:
"""
spin12.jl - sampling robustness for the δf problem

This optimization uses the infidelity metric rather than
the standard diagonal LQR metric.
"""

WDIR = joinpath(@__DIR__, "../../")
include(joinpath(WDIR, "src", "spin", "spin.jl"))

using Altro
using HDF5
using LinearAlgebra
using RobotDynamics
using StaticArrays
using TrajectoryOptimization
const RD = RobotDynamics
const TO = TrajectoryOptimization

# paths
const EXPERIMENT_META = "spin"
const EXPERIMENT_NAME = "spin12"
const SAVE_PATH = abspath(joinpath(WDIR, "out", EXPERIMENT_META, EXPERIMENT_NAME))

# problem
const CONTROL_COUNT = 1
const STATE_COUNT = 2
const ASTATE_SIZE_BASE = STATE_COUNT * HDIM_ISO + 3 * CONTROL_COUNT
const SAMPLE_STATE_COUNT = 4
const SAMPLES_PER_STATE = 2
const SAMPLE_COUNT = SAMPLE_STATE_COUNT * SAMPLES_PER_STATE
const ASTATE_SIZE = ASTATE_SIZE_BASE + SAMPLE_COUNT * HDIM_ISO
const ACONTROL_SIZE = CONTROL_COUNT
# state indices
const STATE1_IDX = 1:HDIM_ISO
const STATE2_IDX = STATE1_IDX[end] + 1:STATE1_IDX[end] + HDIM_ISO
const INTCONTROLS_IDX = STATE2_IDX[end] + 1:STATE2_IDX[end] + CONTROL_COUNT
const CONTROLS_IDX = INTCONTROLS_IDX[end] + 1:INTCONTROLS_IDX[end] + CONTROL_COUNT
const DCONTROLS_IDX = CONTROLS_IDX[end] + 1:CONTROLS_IDX[end] + CONTROL_COUNT
const S1_IDX = DCONTROLS_IDX[end] + 1:DCONTROLS_IDX[end] + HDIM_ISO
const S2_IDX = S1_IDX[end] + 1:S1_IDX[end] + HDIM_ISO
const S3_IDX = S2_IDX[end] + 1:S2_IDX[end] + HDIM_ISO
const S4_IDX = S3_IDX[end] + 1:S3_IDX[end] + HDIM_ISO
const S5_IDX = S4_IDX[end] + 1:S4_IDX[end] + HDIM_ISO
const S6_IDX = S5_IDX[end] + 1:S5_IDX[end] + HDIM_ISO
const S7_IDX = S6_IDX[end] + 1:S6_IDX[end] + HDIM_ISO
const S8_IDX = S7_IDX[end] + 1:S7_IDX[end] + HDIM_ISO
# control indices
const D2CONTROLS_IDX = 1:CONTROL_COUNT

# model
struct Model <: AbstractModel
    h0_samples::Vector{SMatrix{HDIM_ISO, HDIM_ISO}}
end
@inline RD.state_dim(::Model) = ASTATE_SIZE
@inline RD.control_dim(::Model) = ACONTROL_SIZE

# This cost puts a gate error cost on
# the sample states and a LQR cost on the other terms.
# The hessian w.r.t the state and controls is constant.
struct Cost{N,M,T} <: TO.CostFunction
    Q::Diagonal{T, SVector{N,T}}
    R::Diagonal{T, SVector{M,T}}
    q::SVector{N, T}
    c::T
    hess_astate::Symmetric{T, SMatrix{N,N,T}}
    target_states::Array{SVector{HDIM_ISO, T}, 1}
    q_ss1::T
    q_ss2::T
    q_ss3::T
    q_ss4::T
end

function Cost(Q::Diagonal{T,SVector{N,T}}, R::Diagonal{T,SVector{M,T}},
              xf::SVector{N,T}, target_states::Array{SVector{HDIM_ISO}, 1},
              q_ss1::T, q_ss2::T, q_ss3::T, q_ss4::T) where {N,M,T}
    q = -Q * xf
    c = 0.5 * xf' * Q * xf
    hess_astate = zeros(N, N)
    # For reasons unknown to the author, throwing a -1 in front
    # of the gate error Hessian makes the cost function work.
    # This is strange, because the gate error Hessian has been
    # checked against autodiff.
    hess_state1 = -1 * q_ss1 * hessian_gate_error_iso2(target_states[1])
    hess_state2 = -1 * q_ss2 * hessian_gate_error_iso2(target_states[2])
    hess_state3 = -1 * q_ss3 * hessian_gate_error_iso2(target_states[3])
    hess_state4 = -1 * q_ss4 * hessian_gate_error_iso2(target_states[4])
    hess_astate[S1_IDX, S1_IDX] = hess_state1
    hess_astate[S2_IDX, S2_IDX] = hess_state2
    hess_astate[S3_IDX, S3_IDX] = hess_state3
    hess_astate[S4_IDX, S4_IDX] = hess_state4
    hess_astate[S5_IDX, S5_IDX] = hess_state1
    hess_astate[S6_IDX, S6_IDX] = hess_state2
    hess_astate[S7_IDX, S7_IDX] = hess_state3
    hess_astate[S8_IDX, S8_IDX] = hess_state4
    hess_astate += Q
    hess_astate = Symmetric(SMatrix{N, N}(hess_astate))
    return Cost{N,M,T}(Q, R, q, c, hess_astate, target_states, q_ss1, q_ss2, q_ss3, q_ss4)
end

@inline TO.state_dim(cost::Cost{N,M,T}) where {N,M,T} = N
@inline TO.control_dim(cost::Cost{N,M,T}) where {N,M,T} = M
@inline Base.copy(cost::Cost{N,M,T}) where {N,M,T} = Cost{N,M,T}(
    cost.Q, cost.R, cost.q, cost.c, cost.hess_astate,
    cost.target_states, cost.q_ss1, cost.q_ss2, cost.q_ss3, cost.q_ss4
)

@inline TO.stage_cost(cost::Cost{N,M,T}, astate::SVector{N}) where {N,M,T} = (
    0.5 * astate' * cost.Q * astate + cost.q'astate + cost.c
    + cost.q_ss1 * gate_error_iso2(astate, cost.target_states[1]; s1o=S1_IDX[1] - 1)
    + cost.q_ss2 * gate_error_iso2(astate, cost.target_states[2]; s1o=S2_IDX[1] - 1)
    + cost.q_ss3 * gate_error_iso2(astate, cost.target_states[3]; s1o=S3_IDX[1] - 1)
    + cost.q_ss4 * gate_error_iso2(astate, cost.target_states[4]; s1o=S4_IDX[1] - 1)
    + cost.q_ss1 * gate_error_iso2(astate, cost.target_states[1]; s1o=S5_IDX[1] - 1)
    + cost.q_ss2 * gate_error_iso2(astate, cost.target_states[2]; s1o=S6_IDX[1] - 1)
    + cost.q_ss3 * gate_error_iso2(astate, cost.target_states[3]; s1o=S7_IDX[1] - 1)
    + cost.q_ss4 * gate_error_iso2(astate, cost.target_states[4]; s1o=S8_IDX[1] - 1)
)

@inline TO.stage_cost(cost::Cost{N,M,T}, astate::SVector{N},
                      acontrol::SVector{M}) where {N,M,T} = (
    TO.stage_cost(cost, astate) + 0.5 * acontrol' * cost.R * acontrol
)

function TO.gradient!(E::TO.QuadraticCostFunction, cost::Cost{N,M,T},
                      astate::SVector{N,T}) where {N,M,T}
    E.q = (cost.Q * astate + cost.q + [
        @SVector zeros(ASTATE_SIZE_BASE);
        cost.q_ss1 * jacobian_gate_error_iso2(astate, cost.target_states[1]; s1o=S1_IDX[1] - 1);
        cost.q_ss2 * jacobian_gate_error_iso2(astate, cost.target_states[2]; s1o=S2_IDX[1] - 1);
        cost.q_ss3 * jacobian_gate_error_iso2(astate, cost.target_states[3]; s1o=S3_IDX[1] - 1);
        cost.q_ss4 * jacobian_gate_error_iso2(astate, cost.target_states[4]; s1o=S4_IDX[1] - 1);
        cost.q_ss1 * jacobian_gate_error_iso2(astate, cost.target_states[1]; s1o=S5_IDX[1] - 1);
        cost.q_ss2 * jacobian_gate_error_iso2(astate, cost.target_states[2]; s1o=S6_IDX[1] - 1);
        cost.q_ss3 * jacobian_gate_error_iso2(astate, cost.target_states[3]; s1o=S7_IDX[1] - 1);
        cost.q_ss4 * jacobian_gate_error_iso2(astate, cost.target_states[4]; s1o=S8_IDX[1] - 1);
    ])
    return false
end

function TO.gradient!(E::TO.QuadraticCostFunction, cost::Cost{N,M,T}, astate::SVector{N,T},
                      acontrol::SVector{M,T}) where {N,M,T}
    TO.gradient!(E, cost, astate)
    E.r = cost.R * acontrol
    E.c = 0
    return false
end

function TO.hessian!(E::TO.QuadraticCostFunction, cost::Cost{N,M,T}, astate::SVector{N,T}) where {N,M,T}
    E.Q = cost.hess_astate
    return true
end

function TO.hessian!(E::TO.QuadraticCostFunction, cost::Cost{N,M,T}, astate::SVector{N,T},
                     acontrol::SVector{M,T}) where {N,M,T}
    TO.hessian!(E, cost, astate)
    E.R = cost.R
    E.H .= 0
    return true
end

# dynamics
function RD.discrete_dynamics(::Type{RD.RK3}, model::Model, astate::StaticVector,
                              acontrols::StaticVector, time::Real, dt::Real) where {SC}
    negi_hc = astate[CONTROLS_IDX[1]] * NEGI_H1_ISO
    h_prop = exp((FQ_NEGI_H0_ISO + negi_hc) * dt)
    state1 = h_prop * astate[STATE1_IDX]
    state2 = h_prop * astate[STATE2_IDX]
    intcontrols = astate[INTCONTROLS_IDX[1]] + dt * astate[CONTROLS_IDX[1]]
    controls = astate[CONTROLS_IDX[1]] + dt * astate[DCONTROLS_IDX[1]]
    dcontrols = astate[DCONTROLS_IDX[1]] + dt * acontrols[D2CONTROLS_IDX[1]]

    hp_prop = exp((model.h0_samples[1] + negi_hc) * dt)
    hn_prop = exp((model.h0_samples[2] + negi_hc) * dt)
    s1 = hp_prop * astate[S1_IDX]
    s2 = hp_prop * astate[S2_IDX]
    s3 = hp_prop * astate[S3_IDX]
    s4 = hp_prop * astate[S4_IDX]
    s5 = hn_prop * astate[S5_IDX]
    s6 = hn_prop * astate[S6_IDX]
    s7 = hn_prop * astate[S7_IDX]
p    s8 = hn_prop * astate[S8_IDX]

    astate_ = [
        state1; state2; intcontrols; controls; dcontrols;
        s1; s2; s3; s4; s5; s6; s7; s8;
    ]

    return astate_
end

# main
function run_traj(;gate_type=zpiby2, evolution_time=18., solver_type=altro,
                  sqrtbp=false, integrator_type=rk3,
                  qs=[1e0, 1e0, 1e0, 1e-1, 1e0, 1e0, 1e0, 1e0, 1e-1],
                  dt_inv=Int64(1e1), smoke_test=false, constraint_tol=1e-8, al_tol=1e-4,
                  pn_steps=2, max_penalty=1e11, verbose=true, save=true,
                  max_iterations=Int64(2e5), fq_cov=FQ * 1e-2, benchmark=false)
    # model configuration
    h0_samples = Array{SMatrix{HDIM_ISO, HDIM_ISO}}(undef, SAMPLE_COUNT)
    h0_samples[1] = (FQ + fq_cov) * NEGI_H0_ISO
    h0_samples[2] = (FQ - fq_cov) * NEGI_H0_ISO
    model = Model(h0_samples)
    n = state_dim(model)
    m = control_dim(model)
    t0 = 0.
    tf = evolution_time

    # initial state
    x0 = SVector{n}([
        IS1_ISO_;
        IS2_ISO_;
        zeros(3 * CONTROL_COUNT);
        repeat([IS1_ISO_; IS2_ISO_; IS3_ISO_; IS4_ISO_], 2);
    ])

    # final state
    gate_unitary = GT_GATE_ISO[gate_type]
    target_states = Array{SVector{HDIM_ISO}, 1}(undef, 4)
    target_states[1] = gate_unitary * IS1_ISO_
    target_states[2] = gate_unitary * IS2_ISO_
    target_states[3] = gate_unitary * IS3_ISO_
    target_states[4] = gate_unitary * IS4_ISO_
    xf = SVector{n}([
        target_states[1];
        target_states[2];
        zeros(3 * CONTROL_COUNT);
        repeat([target_states[1]; target_states[2];
                target_states[3]; target_states[4]], 2);
    ])

    # control amplitude constraint
    x_max = fill(Inf, n)
    x_max[CONTROLS_IDX] .= MAX_CONTROL_NORM_0
    x_max = SVector{n}(x_max)
    x_min = fill(-Inf, n)
    x_min[CONTROLS_IDX] .= -MAX_CONTROL_NORM_0
    x_min = SVector{n}(x_min)

    # control amplitude constraint at boundary
    x_max_boundary = fill(Inf, n)
    x_max_boundary[CONTROLS_IDX] .= 0
    x_max_boundary = SVector{n}(x_max_boundary)
    x_min_boundary = fill(-Inf, n)
    x_min_boundary[CONTROLS_IDX] .= 0
    x_min_boundary = SVector{n}(x_min_boundary)

    # initial trajectory
    dt = dt_inv^(-1)
    N = Int(floor(evolution_time * dt_inv)) + 1
    U0 = [SVector{m}([
        fill(1e-4, CONTROL_COUNT);
    ]) for k = 1:N-1]
    X0 = [SVector{n}([
        fill(NaN, n);
    ]) for k = 1:N]
    Z = Traj(X0, U0, dt * ones(N))

    # cost function
    Q = Diagonal(SVector{n}([
        fill(qs[1], STATE_COUNT * HDIM_ISO); # ψ1, ψ2
        fill(qs[2], 1); # ∫a
        fill(qs[3], 1); # a
        fill(qs[4], 1); # ∂a
        fill(0, SAMPLE_COUNT * HDIM_ISO);
    ]))
    Qf = Q * N
    R = Diagonal(SVector{m}([
        fill(qs[9], CONTROL_COUNT); # ∂2a
    ]))
    # objective = LQRObjective(Q, R, Qf, xf, N)
    cost_k = Cost(Q, R, xf, target_states, qs[5], qs[6], qs[7], qs[8])
    cost_f = Cost(Qf, R, xf, target_states, N * qs[5], N * qs[6], N * qs[7], N * qs[8])
    objective = TO.Objective(cost_k, cost_f, N)

    # must satisfy control amplitude bound
    control_bnd = BoundConstraint(n, m, x_max=x_max, x_min=x_min)
    # must statisfy conrols start and end at 0
    control_bnd_boundary = BoundConstraint(n, m, x_max=x_max_boundary, x_min=x_min_boundary)
    # must reach target state, must have zero net flux
    target_astate_constraint = GoalConstraint(xf, [STATE1_IDX; STATE2_IDX; INTCONTROLS_IDX])
    # must obey unit norm.
    norm_constraints = [NormConstraint(n, m, 1, TO.Equality(), idxs) for idxs in (
        STATE1_IDX, STATE2_IDX, S1_IDX, S2_IDX, S3_IDX,
        S4_IDX, S5_IDX, S6_IDX, S7_IDX, S8_IDX,
    )]
    constraints = ConstraintList(n, m, N)
    add_constraint!(constraints, control_bnd, 2:N-2)
    add_constraint!(constraints, control_bnd_boundary, N-1:N-1)
    add_constraint!(constraints, target_astate_constraint, N:N);
    for norm_constraint in norm_constraints
        add_constraint!(constraints, norm_constraint, 2:N-1)
    end

    # solve problem
    prob = Problem{IT_RDI[integrator_type]}(model, objective, constraints,
                                            x0, xf, Z, N, t0, evolution_time)
    solver = ALTROSolver(prob)
    verbose_pn = verbose ? true : false
    verbose_ = verbose ? 2 : 0
    projected_newton = solver_type == altro ? true : false
    constraint_tolerance = solver_type == altro ? constraint_tol : al_tol
    iterations_inner = smoke_test ? 1 : 300
    iterations_outer = smoke_test ? 1 : 30
    n_steps = smoke_test ? 1 : pn_steps
    set_options!(solver, square_root=sqrtbp, constraint_tolerance=constraint_tolerance,
                 projected_newton_tolerance=al_tol, n_steps=n_steps,
                 penalty_max=max_penalty, verbose_pn=verbose_pn, verbose=verbose_,
                 projected_newton=projected_newton, iterations_inner=iterations_inner,
                 iterations_outer=iterations_outer, iterations=max_iterations)
    if benchmark
        benchmark_result = Altro.benchmark_solve!(solver)
    else
        benchmark_result = nothing
        Altro.solve!(solver)
    end

    # post-process
    acontrols_raw = TO.controls(solver)
    acontrols_arr = permutedims(reduce(hcat, map(Array, acontrols_raw)), [2, 1])
    astates_raw = TO.states(solver)
    astates_arr = permutedims(reduce(hcat, map(Array, astates_raw)), [2, 1])
    Q_raw = Array(Q)
    Q_arr = [Q_raw[i, i] for i in 1:size(Q_raw)[1]]
    Qf_raw = Array(Qf)
    Qf_arr = [Qf_raw[i, i] for i in 1:size(Qf_raw)[1]]
    R_raw = Array(R)
    R_arr = [R_raw[i, i] for i in 1:size(R_raw)[1]]
    cidx_arr = Array(CONTROLS_IDX)
    d2cidx_arr = Array(D2CONTROLS_IDX)
    cmax = TO.max_violation(solver)
    cmax_info = TO.findmax_violation(TO.get_constraints(solver))
    iterations_ = Altro.iterations(solver)

    result = Dict(
        "acontrols" => acontrols_arr,
        "controls_idx" => cidx_arr,
        "d2controls_dt2_idx" => d2cidx_arr,
        "evolution_time" => evolution_time,
        "astates" => astates_arr,
        "Q" => Q_arr,
        "Qf" => Qf_arr,
        "R" => R_arr,
        "cmax" => cmax,
        "cmax_info" => cmax_info,
        "dt" => dt,
        "sample_count" => SAMPLE_COUNT,
        "solver_type" => Integer(solver_type),
        "sqrtbp" => Integer(sqrtbp),
        "max_penalty" => max_penalty,
        "constraint_tol" => constraint_tol,
        "al_tol" => al_tol,
        "gate_type" => Integer(gate_type),
        "save_type" => Integer(jl),
        "integrator_type" => Integer(integrator_type),
        "iterations" => iterations_,
        "max_iterations" => max_iterations,
    )

    # save
    if save
        save_file_path = generate_file_path("h5", EXPERIMENT_NAME, SAVE_PATH)
        println("Saving this optimization to $(save_file_path)")
        h5open(save_file_path, "cw") do save_file
            for key in keys(result)
                write(save_file, key, result[key])
            end
        end
        result["save_file_path"] = save_file_path
    end

    result = benchmark ? benchmark_result : result

    return result
end

function forward_pass(save_file_path; integrator_type=rk6, gate_type=xpiby2)
    (evolution_time, d2controls, dt
     ) = h5open(save_file_path, "r+") do save_file
         save_type = SaveType(read(save_file, "save_type"))
         if save_type == jl
             d2controls_idx = read(save_file, "d2controls_dt2_idx")
             acontrols = read(save_file, "acontrols")
             d2controls = acontrols[:, d2controls_idx]
             dt = read(save_file, "dt")
             evolution_time = read(save_file, "evolution_time")
         elseif save_type == samplejl
             d2controls = read(save_file, "d2controls_dt2_sample")
             dt = DT_PREF
             ets = read(save_file, "evolution_time_sample")
             evolution_time = Integer(floor(ets / dt)) * dt
         end
         return (evolution_time, d2controls, dt)
     end
    rdi = IT_RDI[integrator_type]
    knot_count = Integer(floor(evolution_time / dt))

    if gate_type == xpiby2
        target_state1 = Array(XPIBY2_ISO_1)
        target_state2 = Array(XPIBY2_ISO_2)
    elseif gate_type == ypiby2
        target_state1 = Array(YPIBY2_ISO_1)
        target_state2 = Array(YPIBY2_ISO_2)
    elseif gate_type == zpiby2
        target_state1 = Array(ZPIBY2_ISO_1)
        target_state2 = Array(ZPIBY2_ISO_2)
    end

    model = Model(sample_count)
    n = state_dim(model)
    m = control_dim(model)
    time = 0.
    astate = SVector{n}([
        IS1;
        IS2;
        zeros(3 * CONTROL_COUNT);
        repeat([IS1; IS2], sample_count);
    ])
    acontrols = [SVector{m}([d2controls[i, 1],]) for i = 1:knot_count - 1]

    for i = 1:knot_count - 1
        astate = RD.discrete_dynamics(rdi, model, astate, acontrols[i], time, dt)
        time = time + dt
    end

    res = Dict(
        "astate" => astate,
        "target_state1" => target_state1,
        "target_state2" => target_state2,
    )

    return res
end

function state_diffs(save_file_path; gate_type=zpiby2)
    (astates,
     ) = h5open(save_file_path, "r") do save_file
        astates = read(save_file, "astates")
        return (astates,)
    end
    knot_count = size(astates, 1)
    fidelities = zeros(SAMPLE_COUNT)
    mse = zeros(SAMPLE_COUNT)
    gate_unitary = GT_GATE[gate_type]
    ts1 = gate_unitary * IS1
    ts2 = gate_unitary * IS2
    ts3 = gate_unitary * IS3
    ts4 = gate_unitary * IS4
    s1 = astates[end, S1_IDX]
    fidelities[1] = fidelity_vec_iso2(s1, ts1)
    d1 = s1 - ts1
    mse[1] = d1'd1
    s2 = astates[end, S2_IDX]
    fidelities[2] = fidelity_vec_iso2(s2, ts2)
    d2 = s2 - ts2
    mse[2] = d2'd2
    s3 = astates[end, S3_IDX]
    fidelities[3] = fidelity_vec_iso2(s3, ts3)
    d3 = s3 - ts3
    mse[3] = d3'd3
    s4 = astates[end, S4_IDX]
    fidelities[4] = fidelity_vec_iso2(s4, ts4)
    d4 = s4 - ts4
    mse[4] = d4'd4
    s5 = astates[end, S5_IDX]
    fidelities[5] = fidelity_vec_iso2(s5, ts1)
    d5 = s5 - ts1
    mse[5] = d5'd5
    s6 = astates[end, S6_IDX]
    fidelities[6] = fidelity_vec_iso2(s6, ts2)
    d6 = s6 - ts2
    mse[6] = d6'd6
    s7 = astates[end, S7_IDX]
    fidelities[7] = fidelity_vec_iso2(s7, ts3)
    d7 = s7 - ts3
    mse[7] = d7'd7
    s8 = astates[end, S8_IDX]
    fidelities[8] = fidelity_vec_iso2(s8, ts4)
    d8 = s8 - ts4
    mse[8] = d8'd8

    return (fidelities, mse)
end

max_stars_repo_name: SchusterLab/rbqoc

max_issues_repo_name: SchusterLab/rbqoc

max_forks_repo_name: SchusterLab/rbqoc


text:
import Molecules.PointGroup: findSEA, buildD, find_rotation_sets, checkSEA, find_rotations, is_there_ortho_c2, find_c2, find_point_group
import Molecules: center_of_mass, translate
using Molecules
using Molecules.PointGroup

r2 = 2.0^0.5
r3 = (3.0^0.5)/2.0
Square = [Atom(1, 1.0, [-0.5, 0.0, 0.0]), Atom(1, 1.0, [0.5, 0.0, 0.0]), Atom(1, 1.0, [0.0, 0.5, 0.0]), Atom(1, 1.0, [0.0, -0.5, 0.0])]
Triangle = [Atom(1, 1.0, [0.0, 0.0, 0.0]), Atom(1, 1.0, [1.0, 0.0, 0.0]), Atom(1, 1.0, [0.5, (3.0^0.5)/2.0, 0.0])]
Dodecagon = [Atom(1, 1.0, [r3, 0.5, 0.0]), Atom(1, 1.0, [0.5, r3, 0.0]), Atom(1, 1.0, [0.0, 1.0, 0.0]), Atom(1, 1.0, [-0.5, r3, 0.0]), Atom(1, 1.0, [-r3, 0.5, 0.0]), Atom(1, 1.0, [-1.0, 0.0, 0.0]), Atom(1, 1.0, [-r3, -0.5, 0.0]), Atom(1, 1.0, [-0.5, -r3, 0.0]), Atom(1, 1.0, [0.0, -1.0, 0.0]), Atom(1, 1.0, [0.5, -r3, 0.0]), Atom(1, 1.0, [r3, -0.5, 0.0]), Atom(1, 1.0, [1.0, 0.0, 0.0])]
Line = [Atom(1, 1.0, [-0.5, 0.0, 0.0]), Atom(1, 1.0, [0.5, 0.0, 0.0])]
RecPrism = [Atom(1, 1.0, [1.0, 1.0, 1.0]), Atom(1, 1.0, [-1.0, 1.0, 1.0]), Atom(1, 1.0, [1.0, -1.0, 1.0]), Atom(1, 1.0, [-1.0, -1.0, 1.0]), Atom(1, 1.0, [1.0, 1.0, -2.0]), Atom(1, 1.0, [-1.0, 1.0, -2.0]), Atom(1, 1.0, [1.0, -1.0, -2.0]), Atom(1, 1.0, [-1.0, -1.0, -2.0])]
RecAnti = [Atom(1, 1.0, [1.0, 1.0, 1.0]), Atom(1, 1.0, [1.0, -1.0, 1.0]), Atom(1, 1.0, [-1.0, 1.0, 1.0]), Atom(1, 1.0, [-1.0, -1.0, 1.0]), Atom(1, 1.0, [r2, 0.0, -2.0]), Atom(1, 1.0, [0.0, r2, -2.0]), Atom(1, 1.0, [-r2, 0.0, -2.0]), Atom(1, 1.0, [0.0, -r2, -2.0])]
RecPrismO = [Atom(1, 1.0, [1.0, 1.0, 1.0]), Atom(1, 1.0, [-1.0, 1.0, 1.0]), Atom(1, 1.0, [1.0, -1.0, 1.0]), Atom(1, 1.0, [-1.0, -1.0, 1.0]), Atom(8, 15.99, [1.0, 1.0, -2.0]), Atom(8, 15.99, [-1.0, 1.0, -2.0]), Atom(8, 15.99, [1.0, -1.0, -2.0]), Atom(8, 15.99, [-1.0, -1.0, -2.0])]

names = ["C1", "Ci", "Cs", "Cs_nonplanar", "C2", "C3", "S4", "C2v", "C3v", "C2h", "C3h", "D2", "D3", "D2d", "D4d", "D3h", "D4h", "D5h", "D6h", "D12h",
         "Cinfv", "Dinfh", "Th", "Td", "cube", "octahedron", "dodecahedron", "icosahedron"]
pgs = ["C1", "Ci", "Cs", "Cs", "C2", "C3", "S4", "C2v", "C3v", "C2h", "C3h", "D2", "D3", "D2d", "D4d", "D3h", "D4h", "D5h", "D6h", "D12h",
       "Cinfv", "Dinfh", "Th", "Td", "Oh", "Oh", "Ih", "Ih"]

@testset "Point Groups" begin
    for i = 1:size(names)[1]
        path = joinpath(@__DIR__, "sxyz/"*names[i]*".xyz")
        mol = Molecules.parse_file(path)
        s = find_point_group(mol)
        @test s == pgs[i]
    end
end
#out = []
#for i = eachindex(mols)
#    path = joinpath(@__DIR__, "xyz/"*mols[i]*".xyz")
#    mol = Molecules.parse_file(path)
#
#    D = buildD(mol)
#
#    sea = findSEA(D, 5)
#    push!(out, sea)
#
#end
#println(out)


text:
const REFERENCE_LOG_10_FREQUENCIES = 1.3:0.1:4.1

const EXPONENT_INTERPOLATOR = CubicSplineInterpolation(
    REFERENCE_LOG_10_FREQUENCIES,
    [
        0.532,
        0.506,
        0.480,
        0.455,
        0.432,
        0.409,
        0.387,
        0.367,
        0.349,
        0.330,
        0.315,
        0.301,
        0.288,
        0.276,
        0.267,
        0.259,
        0.253,
        0.250,
        0.246,
        0.244,
        0.243,
        0.243,
        0.243,
        0.242,
        0.242,
        0.245,
        0.254,
        0.271,
        0.301,
    ],
)

const LINEAR_INTERPOLATOR = CubicSplineInterpolation(
    REFERENCE_LOG_10_FREQUENCIES,
    [
        -31.6,
        -27.2,
        -23.0,
        -19.1,
        -15.9,
        -13.0,
        -10.3,
        -8.1,
        -6.2,
        -4.5,
        -3.1,
        -2.0,
        -1.1,
        -0.4,
        0.0,
        0.3,
        0.5,
        0.0,
        -2.7,
        -4.1,
        -1.0,
        1.7,
        2.5,
        1.2,
        -2.1,
        -7.1,
        -11.2,
        -10.7,
        -3.1,
    ],
)

const THRESHOLD_INTERPOLATOR = CubicSplineInterpolation(
    REFERENCE_LOG_10_FREQUENCIES,
    [
        78.5,
        68.7,
        59.5,
        51.1,
        44.0,
        37.5,
        31.5,
        26.5,
        22.1,
        17.9,
        14.4,
        11.4,
        8.6,
        6.2,
        4.4,
        3.0,
        2.2,
        2.4,
        3.5,
        1.7,
        -1.3,
        -4.2,
        -6.0,
        -5.4,
        -1.5,
        6.0,
        12.6,
        13.9,
        12.3,
    ],
)

function equivalent_sound_pressure(phon, frequency)
    @assert 20 <= phon <= 80
    log_10_frequency = log10(frequency / Hz)
    exponent = EXPONENT_INTERPOLATOR(log_10_frequency)
    linear = LINEAR_INTERPOLATOR(log_10_frequency)
    threshold = THRESHOLD_INTERPOLATOR(log_10_frequency)
    A_f =
        4.47e-3 * (10^(0.025 * phon) - 1.15) +
        (0.4 * (10^((threshold + linear) / 10 - 9)))^exponent
    ((10 / exponent * log10(A_f)) - linear + 94)dB * 20μPa
end

const MEDIUM_LOUDNESS = 40

const MAXIMUM_PRESSURE = equivalent_sound_pressure(MEDIUM_LOUDNESS, 20Hz)

"""
    equal_loudness(synthesizer::Map{<:Any, Tuple{Cycles}})

Change the volume of a synthesizer so that sounds played at different frequencies will have
the same perceived volume. Assumes that the map function has a period of 2π.

```jldoctest
julia> using AudioSchedules

julia> using Unitful: Hz

julia> soft = equal_loudness(Map(cos, Cycles(10000Hz)));

julia> first(make_series(soft, 44100Hz))
0.0007717680651318246
```

Technical details: uses the ISO 226:2003 curve for 40 phons. Scales output by a ratio
of the equivalent sound pressure at the current frequency to the equivalent sound pressure
at 20Hz (about as low as humans can hear).
"""
function equal_loudness(synthesizer::Map{<:Any,Tuple{Cycles}})
    Map(
        let multiplier =
                equivalent_sound_pressure(
                    MEDIUM_LOUDNESS,
                    (synthesizer.synthesizers[1].frequency),
                ) / MAXIMUM_PRESSURE
            function (amplitude)
                amplitude * multiplier
            end
        end,
        synthesizer,
    )
end
export equal_loudness

max_stars_repo_name: bramtayl/AudioSchedulers

max_forks_repo_name: bramtayl/AudioSchedulers


text:
"""
    find_all(A [,a...]; count=false)
A: string/array of elements of the same type

default   : Array containing the index (indices) of selected elements of A (default: all elements)

count=true: The number of indices found for selected elements of A (default: all elements)
#### Examples:
```
A = [:📑,:📌,:📢,:📌,:📞]
B = [1,2,3,2,5]
str = "aβcβd";
find_all(A) == find_all(B) == find_all(str)
true

find_all(A,:📌)
1-element Array{Array{Int64,1},1}:
 [2, 4]

find_all(str)
4-element Array{Array{Int64,1},1}:
 [1]
 [2, 4]
 [3]
 [5]

find_all(A; count=true)
4-element Array{Int64,1}:
 1
 2
 1
 1

str = "📑📌📢📌📞"
find_all(str,'📌')
1-element Array{Array{Int64,1},1}:
 [2, 4]
```
"""
function find_all(A::Union{String,AbstractArray{T,1}}, a::T...; count=false)  where T

    typeof(A) == String ? A = Base.collect(A) : nothing

    a == () ? a = Base.unique(A) : nothing

    o = [Base.findall(x -> x == a[i], A) for i in eachindex(a)]

    return count ? length.(o) : o

end

"""
    find_first(A [,a...]; dict=false)
The first index of selected Array element

A: string/array of elements of the same type

default  : Array containing the first index (indices) of selected elements of A (default: all elements)

dict=true: Dict for the first index (indices) of selected elements of A (default: all elements)
#### Examples:
```
A = [:📑,:📌,:📢,:📌,:📞]
B = [1,2,3,2,5]
str = "aβcβd";

find_first(A) == find_first(B) == find_first(str)
true

find_first(A,:📌)
1-element Array{Array{Int64,1},1}:
 2

find_last(A,:📌; dict=true)
1-element Array{Pair{Symbol,Int64},1}:
 :📌 => 2

find_last(A; dict=true)
4-element Array{Pair{Symbol,Int64},1}:
 :📑 => 1
 :📌 => 2
 :📢 => 3
 :📞 => 5

find_first(str)
4-element Array{Int64,1}:
 1
 2
 3
 5
```
"""
function find_first(A::Union{String,AbstractArray{T,1}}, a::T...; dict=false)  where T

    typeof(A) == String ? A = Base.collect(A) : nothing

    a == () ? a = Base.unique(A) : nothing

    o = [Base.findfirst(x -> x == a[i], A) for i in eachindex(a)]

    return dict ? [a[i] => o[i] for i in eachindex(a)] : o

end

"""
    find_last(A [,a...]; dict=false)

The last index of selected Array element

A: string/array of elements of the same type

default  : Array containing the lasst index (indices) of selected elements of A (default: all elements)

dict=true: Dict for the lasst index (indices) of selected elements of A (default: all elements)
#### Examples:
```
A = [:📑,:📌,:📢,:📌,:📞]
B = [1,2,3,2,5]
str = "aβcβd";
find_last(A) == find_first(B) == find_first(str)
true

find_last(A,:📌)
1-element Array{Array{Int64,1},1}:
 4

find_last(A,:📌; dict=true)
1-element Array{Pair{Symbol,Int64},1}:
 :📌 => 4

find_last(A; dict=true)
4-element Array{Pair{Symbol,Int64},1}:
 :📑 => 1
 :📌 => 4
 :📢 => 3
 :📞 => 5

find_last(str)
4-element Array{Int64,1}:
 1
 4
 3
 5
```
"""
function find_last(A::Union{String,AbstractArray{T,1}}, a::T...; dict=false)  where T

    typeof(A) == String ? A = Base.collect(A) : nothing

    a == () ? a = Base.unique(A) : nothing

    o = [Base.findlast(x -> x == a[i], A) for i ∈ eachindex(a)]

    return dict ? [a[i] => o[i] for i ∈ eachindex(a)] : o

end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: resource_groups_tagging_api
using AWS.Compat
using AWS.UUIDs

"""
    describe_report_creation()
    describe_report_creation(params::Dict{String,<:Any})

Describes the status of the StartReportCreation operation.  You can call this operation
only from the organization's management account and from the us-east-1 Region.

"""
function describe_report_creation(; aws_config::AbstractAWSConfig=global_aws_config())
    return resource_groups_tagging_api(
        "DescribeReportCreation"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function describe_report_creation(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return resource_groups_tagging_api(
        "DescribeReportCreation",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_compliance_summary()
    get_compliance_summary(params::Dict{String,<:Any})

Returns a table that shows counts of resources that are noncompliant with their tag
policies. For more information on tag policies, see Tag Policies in the Organizations User
Guide.  You can call this operation only from the organization's management account and
from the us-east-1 Region. This operation supports pagination, where the response can be
sent in multiple pages. You should check the PaginationToken response parameter to
determine if there are additional results available to return. Repeat the query, passing
the PaginationToken response parameter value as an input to the next request until you
recieve a null value. A null value for PaginationToken indicates that there are no more
results waiting to be returned.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"GroupBy"`: Specifies a list of attributes to group the counts of noncompliant resources
  by. If supplied, the counts are sorted by those attributes.
- `"MaxResults"`: Specifies the maximum number of results to be returned in each page. A
  query can return fewer than this maximum, even if there are more results still to return.
  You should always check the PaginationToken response value to see if there are more
  results. You can specify a minimum of 1 and a maximum value of 100.
- `"PaginationToken"`: Specifies a PaginationToken response value from a previous request
  to indicate that you want the next page of results. Leave this parameter empty in your
  initial request.
- `"RegionFilters"`: Specifies a list of Amazon Web Services Regions to limit the output
  to. If you use this parameter, the count of returned noncompliant resources includes only
  resources in the specified Regions.
- `"ResourceTypeFilters"`: Specifies that you want the response to include information for
  only resources of the specified types. The format of each resource type is
  service[:resourceType]. For example, specifying a resource type of ec2 returns all Amazon
  EC2 resources (which includes EC2 instances). Specifying a resource type of ec2:instance
  returns only EC2 instances. The string for each service name and resource type is the same
  as that embedded in a resource's Amazon Resource Name (ARN). Consult the  Amazon Web
  Services General Reference  for the following:   For a list of service name strings, see
  Amazon Web Services Service Namespaces.   For resource type strings, see Example ARNs.
  For more information about ARNs, see Amazon Resource Names (ARNs) and Amazon Web Services
  Service Namespaces.   You can specify multiple resource types by using a comma separated
  array. The array can include up to 100 items. Note that the length constraint requirement
  applies to each resource type filter.
- `"TagKeyFilters"`: Specifies that you want the response to include information for only
  resources that have tags with the specified tag keys. If you use this parameter, the count
  of returned noncompliant resources includes only resources that have the specified tag keys.
- `"TargetIdFilters"`: Specifies target identifiers (usually, specific account IDs) to
  limit the output by. If you use this parameter, the count of returned noncompliant
  resources includes only resources with the specified target IDs.
"""
function get_compliance_summary(; aws_config::AbstractAWSConfig=global_aws_config())
    return resource_groups_tagging_api(
        "GetComplianceSummary"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function get_compliance_summary(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return resource_groups_tagging_api(
        "GetComplianceSummary",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_resources()
    get_resources(params::Dict{String,<:Any})

Returns all the tagged or previously tagged resources that are located in the specified
Amazon Web Services Region for the account. Depending on what information you want
returned, you can also specify the following:    Filters that specify what tags and
resource types you want returned. The response includes all tags that are associated with
the requested resources.   Information about compliance with the account's effective tag
policy. For more information on tag policies, see Tag Policies in the Organizations User
Guide.    This operation supports pagination, where the response can be sent in multiple
pages. You should check the PaginationToken response parameter to determine if there are
additional results available to return. Repeat the query, passing the PaginationToken
response parameter value as an input to the next request until you recieve a null value. A
null value for PaginationToken indicates that there are no more results waiting to be
returned.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ExcludeCompliantResources"`: Specifies whether to exclude resources that are compliant
  with the tag policy. Set this to true if you are interested in retrieving information on
  noncompliant resources only. You can use this parameter only if the
  IncludeComplianceDetails parameter is also set to true.
- `"IncludeComplianceDetails"`: Specifies whether to include details regarding the
  compliance with the effective tag policy. Set this to true to determine whether resources
  are compliant with the tag policy and to get details.
- `"PaginationToken"`: Specifies a PaginationToken response value from a previous request
  to indicate that you want the next page of results. Leave this parameter empty in your
  initial request.
- `"ResourceARNList"`: Specifies a list of ARNs of resources for which you want to retrieve
  tag data. You can't specify both this parameter and any of the pagination parameters
  (ResourcesPerPage, TagsPerPage, PaginationToken) in the same request. If you specify both,
  you get an Invalid Parameter exception. If a resource specified by this parameter doesn't
  exist, it doesn't generate an error; it simply isn't included in the response. An ARN
  (Amazon Resource Name) uniquely identifies a resource. For more information, see Amazon
  Resource Names (ARNs) and Amazon Web Services Service Namespaces in the Amazon Web Services
  General Reference.
- `"ResourceTypeFilters"`: Specifies the resource types that you want included in the
  response. The format of each resource type is service[:resourceType]. For example,
  specifying a resource type of ec2 returns all Amazon EC2 resources (which includes EC2
  instances). Specifying a resource type of ec2:instance returns only EC2 instances.  The
  string for each service name and resource type is the same as that embedded in a resource's
  Amazon Resource Name (ARN). For the list of services whose resources you can use in this
  parameter, see Services that support the Resource Groups Tagging API. You can specify
  multiple resource types by using an array. The array can include up to 100 items. Note that
  the length constraint requirement applies to each resource type filter. For example, the
  following string would limit the response to only Amazon EC2 instances, Amazon S3 buckets,
  or any Audit Manager resource:  ec2:instance,s3:bucket,auditmanager
- `"ResourcesPerPage"`: Specifies the maximum number of results to be returned in each
  page. A query can return fewer than this maximum, even if there are more results still to
  return. You should always check the PaginationToken response value to see if there are more
  results. You can specify a minimum of 1 and a maximum value of 100.
- `"TagFilters"`: Specifies a list of TagFilters (keys and values) to restrict the output
  to only those resources that have tags with the specified keys and, if included, the
  specified values. Each TagFilter must contain a key with values optional. A request can
  include up to 50 keys, and each key can include up to 20 values.  Note the following when
  deciding how to use TagFilters:   If you don't specify a TagFilter, the response includes
  all resources that are currently tagged or ever had a tag. Resources that currently don't
  have tags are shown with an empty tag set, like this: \"Tags\": [].   If you specify more
  than one filter in a single request, the response returns only those resources that satisfy
  all filters.   If you specify a filter that contains more than one value for a key, the
  response returns resources that match any of the specified values for that key.   If you
  don't specify a value for a key, the response returns all resources that are tagged with
  that key, with any or no value. For example, for the following filters: filter1=
  {keyA,{value1}}, filter2={keyB,{value2,value3,value4}}, filter3= {keyC}:
  GetResources({filter1}) returns resources tagged with key1=value1
  GetResources({filter2}) returns resources tagged with key2=value2 or key2=value3 or
  key2=value4     GetResources({filter3}) returns resources tagged with any tag with the key
  key3, and with any or no value    GetResources({filter1,filter2,filter3}) returns resources
  tagged with (key1=value1) and (key2=value2 or key2=value3 or key2=value4) and (key3, any or
  no value)
- `"TagsPerPage"`: Amazon Web Services recommends using ResourcesPerPage instead of this
  parameter. A limit that restricts the number of tags (key and value pairs) returned by
  GetResources in paginated output. A resource with no tags is counted as having one tag (one
  key and value pair).  GetResources does not split a resource and its associated tags across
  pages. If the specified TagsPerPage would cause such a break, a PaginationToken is returned
  in place of the affected resource and its tags. Use that token in another request to get
  the remaining data. For example, if you specify a TagsPerPage of 100 and the account has 22
  resources with 10 tags each (meaning that each resource has 10 key and value pairs), the
  output will consist of three pages. The first page displays the first 10 resources, each
  with its 10 tags. The second page displays the next 10 resources, each with its 10 tags.
  The third page displays the remaining 2 resources, each with its 10 tags. You can set
  TagsPerPage to a minimum of 100 items up to a maximum of 500 items.
"""
function get_resources(; aws_config::AbstractAWSConfig=global_aws_config())
    return resource_groups_tagging_api(
        "GetResources"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function get_resources(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return resource_groups_tagging_api(
        "GetResources", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

"""
    get_tag_keys()
    get_tag_keys(params::Dict{String,<:Any})

Returns all tag keys currently in use in the specified Amazon Web Services Region for the
calling account. This operation supports pagination, where the response can be sent in
multiple pages. You should check the PaginationToken response parameter to determine if
there are additional results available to return. Repeat the query, passing the
PaginationToken response parameter value as an input to the next request until you recieve
a null value. A null value for PaginationToken indicates that there are no more results
waiting to be returned.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"PaginationToken"`: Specifies a PaginationToken response value from a previous request
  to indicate that you want the next page of results. Leave this parameter empty in your
  initial request.
"""
function get_tag_keys(; aws_config::AbstractAWSConfig=global_aws_config())
    return resource_groups_tagging_api(
        "GetTagKeys"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function get_tag_keys(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return resource_groups_tagging_api(
        "GetTagKeys", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

"""
    get_tag_values(key)
    get_tag_values(key, params::Dict{String,<:Any})

Returns all tag values for the specified key that are used in the specified Amazon Web
Services Region for the calling account. This operation supports pagination, where the
response can be sent in multiple pages. You should check the PaginationToken response
parameter to determine if there are additional results available to return. Repeat the
query, passing the PaginationToken response parameter value as an input to the next request
until you recieve a null value. A null value for PaginationToken indicates that there are
no more results waiting to be returned.

# Arguments
- `key`: Specifies the tag key for which you want to list all existing values that are
  currently used in the specified Amazon Web Services Region for the calling account.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"PaginationToken"`: Specifies a PaginationToken response value from a previous request
  to indicate that you want the next page of results. Leave this parameter empty in your
  initial request.
"""
function get_tag_values(Key; aws_config::AbstractAWSConfig=global_aws_config())
    return resource_groups_tagging_api(
        "GetTagValues",
        Dict{String,Any}("Key" => Key);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_tag_values(
    Key, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return resource_groups_tagging_api(
        "GetTagValues",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("Key" => Key), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    start_report_creation(s3_bucket)
    start_report_creation(s3_bucket, params::Dict{String,<:Any})

Generates a report that lists all tagged resources in the accounts across your organization
and tells whether each resource is compliant with the effective tag policy. Compliance data
is refreshed daily. The report is generated asynchronously. The generated report is saved
to the following location:
s3://example-bucket/AwsTagPolicies/o-exampleorgid/YYYY-MM-ddTHH:mm:ssZ/report.csv  You can
call this operation only from the organization's management account and from the us-east-1
Region.

# Arguments
- `s3_bucket`: The name of the Amazon S3 bucket where the report will be stored; for
  example:  awsexamplebucket  For more information on S3 bucket requirements, including an
  example bucket policy, see the example S3 bucket policy on this page.

"""
function start_report_creation(S3Bucket; aws_config::AbstractAWSConfig=global_aws_config())
    return resource_groups_tagging_api(
        "StartReportCreation",
        Dict{String,Any}("S3Bucket" => S3Bucket);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function start_report_creation(
    S3Bucket,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return resource_groups_tagging_api(
        "StartReportCreation",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("S3Bucket" => S3Bucket), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    tag_resources(resource_arnlist, tags)
    tag_resources(resource_arnlist, tags, params::Dict{String,<:Any})

Applies one or more tags to the specified resources. Note the following:   Not all
resources can have tags. For a list of services with resources that support tagging using
this operation, see Services that support the Resource Groups Tagging API. If the resource
doesn't yet support this operation, the resource's service might support tagging using its
own API operations. For more information, refer to the documentation for that service.
Each resource can have up to 50 tags. For other limits, see Tag Naming and Usage
Conventions in the Amazon Web Services General Reference.    You can only tag resources
that are located in the specified Amazon Web Services Region for the Amazon Web Services
account.   To add tags to a resource, you need the necessary permissions for the service
that the resource belongs to as well as permissions for adding tags. For more information,
see the documentation for each service.    Do not store personally identifiable information
(PII) or other confidential or sensitive information in tags. We use tags to provide you
with billing and administration services. Tags are not intended to be used for private or
sensitive data.   Minimum permissions  In addition to the tag:TagResources permission
required by this operation, you must also have the tagging permission defined by the
service that created the resource. For example, to tag an Amazon EC2 instance using the
TagResources operation, you must have both of the following permissions:    tag:TagResource
    ec2:CreateTags

# Arguments
- `resource_arnlist`: Specifies the list of ARNs of the resources that you want to apply
  tags to. An ARN (Amazon Resource Name) uniquely identifies a resource. For more
  information, see Amazon Resource Names (ARNs) and Amazon Web Services Service Namespaces in
  the Amazon Web Services General Reference.
- `tags`: Specifies a list of tags that you want to add to the specified resources. A tag
  consists of a key and a value that you define.

"""
function tag_resources(
    ResourceARNList, Tags; aws_config::AbstractAWSConfig=global_aws_config()
)
    return resource_groups_tagging_api(
        "TagResources",
        Dict{String,Any}("ResourceARNList" => ResourceARNList, "Tags" => Tags);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function tag_resources(
    ResourceARNList,
    Tags,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return resource_groups_tagging_api(
        "TagResources",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("ResourceARNList" => ResourceARNList, "Tags" => Tags),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    untag_resources(resource_arnlist, tag_keys)
    untag_resources(resource_arnlist, tag_keys, params::Dict{String,<:Any})

Removes the specified tags from the specified resources. When you specify a tag key, the
action removes both that key and its associated value. The operation succeeds even if you
attempt to remove tags from a resource that were already removed. Note the following:   To
remove tags from a resource, you need the necessary permissions for the service that the
resource belongs to as well as permissions for removing tags. For more information, see the
documentation for the service whose resource you want to untag.   You can only tag
resources that are located in the specified Amazon Web Services Region for the calling
Amazon Web Services account.    Minimum permissions  In addition to the tag:UntagResources
permission required by this operation, you must also have the remove tags permission
defined by the service that created the resource. For example, to remove the tags from an
Amazon EC2 instance using the UntagResources operation, you must have both of the following
permissions:    tag:UntagResource     ec2:DeleteTags

# Arguments
- `resource_arnlist`: Specifies a list of ARNs of the resources that you want to remove
  tags from. An ARN (Amazon Resource Name) uniquely identifies a resource. For more
  information, see Amazon Resource Names (ARNs) and Amazon Web Services Service Namespaces in
  the Amazon Web Services General Reference.
- `tag_keys`: Specifies a list of tag keys that you want to remove from the specified
  resources.

"""
function untag_resources(
    ResourceARNList, TagKeys; aws_config::AbstractAWSConfig=global_aws_config()
)
    return resource_groups_tagging_api(
        "UntagResources",
        Dict{String,Any}("ResourceARNList" => ResourceARNList, "TagKeys" => TagKeys);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function untag_resources(
    ResourceARNList,
    TagKeys,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return resource_groups_tagging_api(
        "UntagResources",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "ResourceARNList" => ResourceARNList, "TagKeys" => TagKeys
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
##### MEAN-VARIANCE OPTIMIZATION
#This is not supposed to produce great results, since mean and covariance
#is supposed to be very exact before good weights are found
using JuMP
using Gurobi
using StatsBase
using DataFrames
using CSV
include("SupportFunction.jl")
include("DataLoad.jl")
include("CEOSupportFunctions.jl")
println("Leeeeroooy Jenkins")
#Esben's path
#path = "$(homedir())/Documents/GitHub/Thesis/Data/IndexData"

#inputArg1 = 1   # should range from 1:5
#inputArg2 = 0   # should range from 0:24
inputArg1 = parse(Int64, ARGS[1])
inputArg2 = parse(Int64, ARGS[2])

path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/Data/IndexDataDiff/"
path = "/zhome/9f/d/88706/SpecialeCode/Thesis/Data/IndexDataDiff/"

trainingSize = 240
possibilities = 5
industries = ["NoDur", "Durbl", "Manuf", "Enrgy", "HiTec", "Telcm", "Shops", "Hlth", "Utils", "Other"]
industriesTotal = length(industries)

modelMatrix = zeros(industriesTotal, possibilities)
#noDurModel = [1 0 1 1 1]
#noDurModel = [1 0 1 0 1]
noDurModel = [1 0 0 1 1]
testModel = [1 0 0 1 1]
modelMatrix[1, :] = noDurModel
for i=2:industriesTotal
    modelMatrix[i, :] = noDurModel
end

XArrays = Array{Array{Float64, 2}}(industriesTotal)
YArrays = Array{Array{Float64, 2}}(industriesTotal)

XArrays, YArrays = generateXandYs(industries, modelMatrix)

nGammas = 5
standY = YArrays[1]
SSTO = sum((standY[i]-mean(standY[:]))^2 for i=1:length(standY))
lambdaValues = log.(logspace(100, SSTO/2, nGammas))

lambda4Values = [1e4,1e5,1e6,1e7,1e8]

nRows = size(standY)[1]
amountOfModels = nGammas^4

modelConfig = zeros(amountOfModels, 4)
counter = 1
for l1 = 1:nGammas
    for l2 = 1:nGammas
        for l3 = 1:nGammas
            for l4 = 1:nGammas
                #modelConfig[counter,:] = [lambdaValues[l1] lambdaValues[l2] lambdaValues[l3] lambdaValues[l4]]
                modelConfig[counter,:] = [lambdaValues[l1] lambdaValues[l2] lambdaValues[l3] lambda4Values[l4]]
                counter += 1
            end
        end
    end
end
modelConfig

#Initialization of parameters
w1N = repeat([1/11], outer = 11) #1/N weights
gamma = 2.4 #risk aversion
validationPeriod = 5
PMatrix = zeros(nRows-trainingSize, amountOfModels)
return1NMatrix = zeros(nRows-trainingSize, amountOfModels)
returnCEOMatrix = zeros(nRows-trainingSize, amountOfModels)
returnPerfectMatrix = zeros(nRows-trainingSize)

bestModelAmount = 5
bestModelConfigs = zeros(bestModelAmount, 4)
bestModelIndexes = zeros(bestModelAmount)

weightsPerfect = zeros(nRows-trainingSize, 11)
weightsCEO     = zeros(nRows-trainingSize, 11, amountOfModels)
#expectedReturnMatrix = zeros(nRows-trainingSize, 11)
#forecastErrors = zeros(nRows-trainingSize, 11)

rfRates = loadRiskFreeRate("NoDur", path)
rfRates = rfRates[:,1]

#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/CEO/"
path = "/zhome/9f/d/88706/SpecialeCode/Results/CEORFR/"
weightsPerfect = Array{Float64}(CSV.read(path*"weightsPerfect.csv",header=false, datarow=1, nullable=false))
returnPerfectMatrix = Array{Float64}(CSV.read(path*"returnPerfectMatrix.csv", header=false, datarow=1, nullable=false))

println("Starting CEO Validation loop")
t=inputArg1
trainingXArrays, trainingYArrays, validationXRows, validationY, OOSXArrays, OOSYArrays, OOSRow, OOSY = createDataSplits(XArrays, YArrays, t, trainingSize)

for m = 1+(inputArg2*25):25+(inputArg2*25)#:amountOfModels
    expectedReturns = zeros(11)
    println(m-(inputArg2*25), " out of ", 25)

    #CHANGES
    rfRatesVec = rfRates[t:(t+trainingSize-1)]
    betaArray, U = @time(runCEORFR(trainingXArrays, trainingYArrays, modelConfig[m, :], gamma, rfRatesVec))

    ##PREVIOUS
    #betaArray, U = @time(runCEO(trainingXArrays, trainingYArrays, modelConfig[m, :], gamma))

    expectedReturns[1:10] = generateExpectedReturns(betaArray, trainingXArrays, trainingYArrays, validationXRows)
    expectedReturns[11] = rfRates[t+trainingSize]

    #Need to send OOSRow to mean-variance optimization to get "perfect information" since validationY is the values in OOSRow[1:10]
    valY = zeros(11)
    for i = 1:10
        valY[i] = validationY[i][1]
    end
    valY[11] = rfRates[t+trainingSize]
    #return1N, returnCEO, wStar = performMVOptimization(expectedReturns, U, gamma, OOSRow[1][1:10], valY)
    rfRatesVec = rfRates[t:(t+trainingSize-1)]
    trainX = hcat(trainingXArrays[1][:,1:10], rfRatesVec)
    Sigma =  cov(trainX)
    F = lufact(Sigma)
    U = F[:U]  #Cholesky factorization of Sigma

    return1N, returnCEO, wStar = performMVOptimizationRISK(expectedReturns, U, gamma, valY, valY)
    weightsCEO[t, 1:11, m]        = wStar
    return1NMatrix[t, m]       = return1N
    returnCEOMatrix[t, m]      = returnCEO

    returnPerfect = returnPerfectMatrix[t]
    println("1N returns is $return1N, returnPerfect is $returnPerfect and returnCEO is $returnCEO")
    PMatrix[t, m] = calculatePvalue(return1N, returnPerfect, returnCEO)
    #trackReturn(returnCEOTotal, returnCEO)
end

# Writing files to be used in CEOParallelPart2
#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/CEO/"
path = "/zhome/9f/d/88706/SpecialeCode/Results/CEORFR/"
writedlm(path*string(inputArg1)*"_"*string(inputArg2)*"_weightsCEO.csv", weightsPerfect[t,:], ",") # do additional stuff here
writedlm(path*string(inputArg1)*"_"*string(inputArg2)*"_return1NMatrix.csv", return1NMatrix[t,:], ",")
writedlm(path*string(inputArg1)*"_"*string(inputArg2)*"_returnCEOMatrix.csv", returnCEOMatrix[t,:], ",")
writedlm(path*string(inputArg1)*"_"*string(inputArg2)*"_PMatrix.csv", PMatrix[t,:], ",")
#writedlm(path*string(inputArg1)*"_"*string(inputArg2)*"_ceoRFRErrors.csv", forecastErrors[t,:], ",")
#writedlm(path*string(inputArg1)*"_"*string(inputArg2)*"_ceoRFRForecasts.csv", expectedReturnMatrix[t,:], ",")

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_name: skipperafrosenborg/Thesis


text:
@testset "LocalField ..." begin
  include("LocalField/Conjugates.jl")
  include("LocalField/Poly.jl")
  include("LocalField/LocalField.jl")
  include("LocalField/neq.jl")
end


text:
#=

  The BUGS Book
  Example 3.3.2 Surgery (continued): beta-binomial analysis using BUGS

  Cf ~/jags/bugs_book_3_3_2.R
  """
  model {
    theta ~ dbeta(a, b)
    y ~ dbin(theta, n)
    Y.pred ~ dbin(theta, n.pred)
    P.crit <- step(Y.pred - n.crit + 0.5)
  }

  Output:
          Mean      SD  Naive SE Time-series SE
P.crit 0.41531 0.49278 0.0017422      0.0019429
Y.pred 1.49477 1.42409 0.0050349      0.0061002
theta  0.07453 0.04117 0.0001456      0.0002249
  """

  Summary Statistics
  parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

       theta    0.0748    0.0416     0.0004    0.0005   7624.8902    1.0000     3450.1765
      y_pred    1.4931    1.4276     0.0143    0.0144   9178.6177    1.0000     4153.2207
      p_crit    0.4177    0.4932     0.0049    0.0051   9070.0308    1.0001     4104.0863

  Cf ~/webppl/bugs_book_3_3_2.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function bugs_book_3_3_2(y=0)
    a = 3
    b = 27
    n = 10
    n_pred = 20
    n_crit = 2

    theta ~ Beta(a,b)
    y ~ Binomial(n,theta)
    y_pred ~ Binomial(n_pred,theta)
    p_crit ~ Dirac(y_pred + 0.5 > n_crit) #  toReal(y_pred - n_crit) + 0.5 > 0.0;

end

model = bugs_book_3_3_2()

num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 100_000)
# chns = sample(model, PG(15), 10_000)
chns = sample(model, SMC(10_000), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, NUTS(), 10_000)
# chns = sample(model, HMC(0.1,5), 10_000)

display(chns)
# display(plot(chns))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# SPDX-License-Identifier: MIT
# Author: Pierre Lairez

module Codegen

using IntervalArithmetic
using StaticArrays

import Base: +, *, -, one, convert, promote_rule

export coord, @genpredicate

@enum FilterEnum fastfp_flt zerotest_flt accuratefp_flt interval_flt exact_flt naive_flt nothing_flt

coord(x) = Tuple(x)
coord(x :: Complex) = reim(x)

function adddict(a, b)
    d = copy(a)
    for k in keys(b)
        if haskey(a, k)
            d[k] += b[k]
        else
            d[k] = b[k]
        end
    end
    return d
end

function mergedict(a, b)
    d = copy(a)
    for k in keys(b)
        if haskey(a, k)
            d[k] = union(d[k], b[k])
        else
            d[k] = b[k]
        end
    end
    return d
end

mutable struct Formula <: Number
    id :: Symbol
    head :: Symbol
    args :: Vector
    group :: Union{Nothing, Symbol}
end

function Formula(name :: Union{Expr, Symbol}, group = nothing)
    s = gensym()
    Formula(s, :sym, [name], group)
end

macro var(args...)
    vars = [Formula(v) for v in args]
    Expr(:block,
         [Expr(:(=), esc(args[i]), vars[i]) for i in 1:length(vars)]...)
end

macro point2(args...)
    vars = [SVector(Formula(:($v[1])), Formula(:($v[2]))) for v in args]
    Expr(:block,
         [Expr(:(=), esc(args[i]), vars[i]) for i in 1:length(vars)]...)
end

function +(f :: Formula, g :: Formula)
    Formula(gensym(), :+, [f, g], nothing)
end

function *(f :: Formula, g :: Formula)
    Formula(gensym(), :*, [f, g], nothing)
end

function -(f :: Formula, g :: Formula)
    Formula(gensym(), :-, [f, g], nothing)
end

function -(f :: Formula)
    Formula(gensym(), :-, [f], nothing)
end

function convert(::Type{Formula}, i :: T) where T <: Integer
    Formula(gensym(), :const, [i], nothing)
end

zero(::Type{Formula}) = convert(Formula, 0)

promote_rule(::Type{T}, ::Type{Formula}) where {T <: Integer} = Formula

function group!(args...)
    @gensym g
    for f in args
        f.group = g
    end
end

function accumulator(f :: Formula)
    if f.head == :sym
        ret = (deg = Dict(), groups = Dict(), bound = Inf, error = 0.0)
    elseif f.head == :+ || (f.head == :- && length(f.args) == 2)
        af = accumulator(f.args[1])
        ag = accumulator(f.args[2])
        @assert af.deg == ag.deg

        bound = nextfloat(af.bound + ag.bound)
        error = nextfloat(nextfloat(af.error + ag.error) + eps(bound)/2)
        ret = (deg = af.deg, groups = mergedict(af.groups, ag.groups), bound = bound, error = error)
    elseif f.head == :*
        af = accumulator(f.args[1])
        ag = accumulator(f.args[2])
        bound = nextfloat(af.bound * ag.bound)
        error = nextfloat(nextfloat(nextfloat(af.error*ag.bound) + nextfloat(af.bound*ag.error)) + eps(bound)/2)
        ret = (deg = adddict(af.deg, ag.deg), groups = mergedict(af.groups, ag.groups), bound = bound, error = error)
    elseif f.head == :-
        ret = accumulator(f.args[1])
    elseif f.head == :const
        ret = (deg = Dict(), groups = Dict(), bound = abs(float(f.args[1])), error = 0.0)
    else
        throw(DomainError())
    end

    if !isnothing(f.group)
        @assert f.head == :sym || f.head == :- && f.args[1].head == f.args[2].head == :sym
        ret = (deg = Dict(f.group => 1), groups = Dict(f.group => Set([f.id])), bound = 1.0, error = eps(1.0)/2)
    end

    return ret
end

function evalcode(f :: Formula, conv = nothing)
    code = []
    stack = [f]
    syms = Set{Symbol}()

    while !isempty(stack)
        e = pop!(stack)

        if e.id ∈ syms
            continue
        end

        s = (e.id)

        if e.head == :sym
            if isnothing(conv)
                push!(code, Expr(:(=), s, e.args[1]))
            else
                push!(code, Expr(:(=), s, conv(e.args[1])))
            end
        elseif e.head == :const
            push!(code, code, Expr(:(=), s, e.args[1]))
        else
            if all(a.id ∈ syms for a in e.args)
                push!(code, Expr(:(=), s, Expr(:call, e.head, ((a.id) for a in e.args)...)))
            else
                push!(stack, e, e.args...)
                continue
            end
        end

        push!(syms, e.id)
    end

    return quote
        $(code...)
        $(f.id)
    end
end

function fastfilter(f :: Formula ; withretcode :: Bool = false)
    code = []

    @gensym res

    let
        fpcode = quote
            $res = $(evalcode(f))
        end
        push!(code, fpcode)
    end

    # Now code contains the computation of the formula, nothing more. We can
    # reliably determine the sign of the result if its absolute value is larger
    # than acc.error*scaling. The scaling depends on the homogeneity structure
    # of the formula.

    # The format of acc.groups is Dict(g1=>[a,b,c,...], g2=>[d,e,f,...], ...),
    # where a, b, c, ... are input gates (or difference or inpur gates, which
    # are considered as new input gates). The formula is expected to be
    # multihomogeneous w.r.t. to each group. The degree of homogeneity w.r.t to
    # the group g is acc.deg[g]. Therefore, the appropriate scaling is the
    # product of all λ[g]^acc.deg[g], where λ[g] is the maximum
    # absolute value of an input in group g.

    # The accurate way of determining if abs(res) >
    # acc.error*Π λ[g]^acc.deg[g] is to do the computation with fp numbers.
    # But we can check quickly a sufficient condition just with integer
    # arithmetic on the exponents.

    # There are two fundamental things to check: that no overflow occured with computing res
    # and that no underflow occurs when computing acc.error*Π λ[g]^acc.deg[g].

    acc = accumulator(f)
    totdeg = UInt(sum(values(acc.deg)))

    @gensym signres

    let
        # we first try to determine if the absolute value of the result is
        # bigger than the maximal possible absolute error by looking only at the
        # exponents.

        emask = Base.exponent_mask(Float64)
        mantissabits = UInt(Base.Math.significand_bits(Float64))
        grouplogs = []
        for (idx, g) in acc.groups
            @gensym glog
            push!(grouplogs, glog)

            alllogs = [:(reinterpret(UInt64, $v) & $emask) for v in g]
            push!(code, :($glog = max($(alllogs...)) >> $mantissabits))
            # we could probably work without the shift, but it would make it
            # harder to think about overflows/underflows.

            if acc.deg[idx] != 1
                push!(code, :($glog *= $(acc.deg[idx])))
            end
        end

        errlog = ((reinterpret(UInt64, acc.error) & emask) >> mantissabits) -
            totdeg * ( Base.exponent_one(Float64) >> mantissabits ) + totdeg

        @gensym ε
        @gensym reslog
        @gensym rawres

        filter = quote
            $ε = $(Expr(:call, :+, errlog, grouplogs...))
            $rawres = reinterpret(UInt64, $res)
            $signres = Base.ifelse($rawres & $(Base.sign_mask(Float64)) == 0, 1, -1)
            $reslog = ($rawres & $emask) >> $mantissabits

            if $reslog != $(emask >> mantissabits) && # protect against Inf and Nan
                #$ε ≥ 0 &&       # protect against underflow in computation of ε
                $reslog > $ε    # main test
                # we should check that ε is not
                # negative. But since we are in unsigned arithmetic, ε < 0 means
                # actually that ε is big, so reslog > ε cannot hold.

                $(withretcode ? :(return ($signres, $fastfp_flt)) : :(return $signres))
            end
        end
        push!(code, filter)
    end

    let
        # Quick test is non conclusive. We check if there is a group in which
        # all the variables are zero, in which case the result is zero.

        filter = quote
            if $(Expr(:call, :|,
                      (Expr(:call, :&,
                            (Expr(:call, :(==), v, 0.0)
                             for v in g)...)
                       for g in values(acc.groups))...)
                 )

                $(withretcode ? :(return (0, $zerotest_flt)) : :(return 0))
            end
        end
        push!(code, filter)
    end

    let
        # We now go to the finer, more straight forward, test.

        groupabs = []
        for (idx, g) in acc.groups
            @gensym gabs
            push!(groupabs, gabs)
            allabs = [:(abs($v)) for v in g]
            push!(code, :($gabs = max($(allabs...))))

            if acc.deg[idx] != 1
                push!(code, :($gabs = $gabs^$(acc.deg[idx])))
            end
        end

        errabs =  acc.error * (1+eps(1.0))^totdeg

        @gensym ε

        filter = quote
            $ε = *($(groupabs...)) * $errabs
            if !issubnormal($ε) && $ε > 0 && isfinite($res) && abs($res) > $ε
                $(withretcode ? :(return ($signres, $accuratefp_flt)) : :(return $signres))
            end
        end
        push!(code, filter)
    end

    return quote
        $(code...)
    end

end

function ivfilter(f :: Formula ; withretcode :: Bool = false)

    @gensym ivres
    quote
        # Default rounding mode may have problems with underflows
        # See https://github.com/JuliaIntervals/IntervalArithmetic.jl/issues/302
        setrounding(Interval, :slow)
        # We now resort to interval arithmetic It is an interesting filter when
        # the data is made of exactly representable integers.
        $ivres = $(evalcode(f,  s -> :( interval($s) )))
        if $ivres < 0
            $(withretcode ? :(return (-1, $interval_flt)) : :(return -1))
        elseif $ivres > 0
            $(withretcode ? :(return (1, $interval_flt)) : :(return 1))
        elseif $ivres == 0
            $(withretcode ? :(return (0, $interval_flt)) : :(return 0))
        end
    end

end

function exfilter(f :: Formula ; withretcode :: Bool = false)

    @gensym exrec
    quote
        # Exact arithmetic. Always conclusive.
        $exrec = $(evalcode(f, s -> :( Rational{BigInt}($s) )))
        $(withretcode ? :(return (Int(sign($exrec)), $exact_flt)) : :(return Int(sign($exrec))))
    end

end

function naivefilter(f :: Formula  ; withretcode :: Bool = false)
    @gensym fpres
    quote
        # Exact arithmetic. Always conclusive.
        $fpres = $(evalcode(f))
        if $fpres < 0
            $(withretcode ? :(return (-1, $naive_flt)) : :(return -1))
        elseif $fpres > 0
            $(withretcode ? :(return (1, $naive_flt)) : :(return 1))
        else
            $(withretcode ? :(return (0, $naive_flt)) : :(return 0))
        end
    end
end

"""

Generate sign predicate for a function that compute a polynomial in the
coordinates of the arguments.

"""
macro genpredicate(args...)
    if first(args) == :nogeneric
        defgeneric = false
        fun = args[2]
    else
        defgeneric = true
        fun = args[1]
    end

    sig = fun.args[1]
    args = sig.args[2:length(sig.args)]

    @assert fun.head == :function

    nargs = []
    tupleconv = []
    input = []
    for a in args
        if isa(a, Symbol)
            push!(input, Formula(a))
            push!(nargs, :($a :: Float64))
        elseif isa(a, Expr) && a.head == :(::)
            v, dim = a.args
            push!(input, SVector((Formula(:($v[$i])) for i in 1:dim)...))
            push!(nargs, :($v :: NTuple{$dim, Float64}))
            push!(tupleconv, :($v = coord($v)))
        else
            throw(DomainError("Unknown argument $a"))
        end
    end

    base = string(sig.args[1])
    mainf = esc(Symbol(base))
    slowf = esc(Symbol(base, "_slow"))
    referencef = esc(Symbol(base, "_reference"))
    naivef = esc(Symbol(base, "_naive"))

    debug = s -> esc(Symbol(s.args[1], "_dbg"))

    nsig = (a for a in nargs)
    vars = (a.args[1] for a in nargs)

    # This is the fun part: we inject into the body of the function given in
    # argument the formal input. This gives an object of type Formula that
    # represents the polynomial.
    formula = Core.eval(__module__,
         Expr(:call,
              Expr(:->, Expr(:tuple, vars...), fun.args[2]),
              input...))

    if defgeneric
        genfun = quote
            function $(mainf)($(vars...))
                $(tupleconv...)
                return $(mainf)($(vars...))
            end
        end
    else
        genfun = quote end
    end

    quote
        function $(naivef)($(nsig...))
            $(naivefilter(formula))
        end

        function $(debug(naivef))($(nsig...))
            $(naivefilter(formula, withretcode=true))
        end

        function $(referencef)($(nsig...))
            $(exfilter(formula))
        end

        function $(debug(referencef))($(nsig...))
            $(exfilter(formula, withretcode=true))
        end

        function $(slowf)($(nsig...))
            $(ivfilter(formula))
            return $(referencef)($(vars...))
        end

        function $(debug(slowf))($(nsig...))
            $(ivfilter(formula, withretcode=true))
            return $(debug(referencef))($(vars...))
        end

        Core.@__doc__(@inline function $(mainf)($(nsig...))
            $(fastfilter(formula))
            return $(slowf)($(vars...))
        end)

        function $(debug(mainf))($(nsig...))
            $(fastfilter(formula, withretcode=true))
            return $(debug(slowf))($(vars...))
        end

        $(genfun)

    end

end

end