text:
struct Cost{T}
    evaluate::Any
    gradient::Any
    hessian::Any
    num_gradient::Int
    num_hessian::Int
    sparsity::Vector{Vector{Int}}
    evaluate_cache::Vector{T}
    gradient_cache::Vector{T}
    hessian_cache::Vector{T}
end

function Cost(f::Function, num_state::Int, num_action::Int;
    num_parameter::Int=0,
    evaluate_hessian=false)

    #TODO: option to load/save methods
    @variables x[1:num_state], u[1:num_action], w[1:num_parameter]
    evaluate = f(x, u, w)
    gradient = Symbolics.gradient(evaluate, [x; u])
    num_gradient = num_state + num_action
    evaluate_func = eval(Symbolics.build_function([evaluate], x, u, w)[2])
    gradient_func = eval(Symbolics.build_function(gradient, x, u, w)[2])
    if evaluate_hessian
        hessian = Symbolics.sparsehessian(evaluate, [x; u])
        hessian_func = eval(Symbolics.build_function(hessian.nzval, x, u, w)[2])
        sparsity = [findnz(hessian)[1:2]...]
        num_hessian = length(hessian.nzval)
    else
        hessian_func = Expr(:null)
        sparsity = [Int[]]
        num_hessian = 0
    end

    return Cost(
        evaluate_func,
        gradient_func,
        hessian_func,
        num_gradient,
        num_hessian,
        sparsity,
        zeros(1),
        zeros(num_gradient),
        zeros(num_hessian))
end

Objective{T} = Vector{Cost{T}} where T

function cost(objective::Objective, states, actions, parameters)
    J = 0.0
    for (t, cost) in enumerate(objective)
        cost.evaluate(cost.evaluate_cache, states[t], actions[t], parameters[t])
        J += cost.evaluate_cache[1]
    end
    return J
end

function gradient!(gradient, indices, objective::Objective, states, actions, parameters)
    for (t, cost) in enumerate(objective)
        cost.gradient(cost.gradient_cache, states[t], actions[t], parameters[t])
        @views gradient[indices[t]] .+= cost.gradient_cache
        # fill!(cost.gradient_cache, 0.0) # TODO: confirm this is necessary
    end
end

function hessian!(hessian, indices, objective::Objective, states, actions, parameters, scaling)
    for (t, cost) in enumerate(objective)
        cost.hessian(cost.hessian_cache, states[t], actions[t], parameters[t])
        cost.hessian_cache .*= scaling
        @views hessian[indices[t]] .+= cost.hessian_cache
        # fill!(cost.hessian_cache, 0.0) # TODO: confirm this is necessary
    end
end

function sparsity_hessian(objective::Objective, num_state::Vector{Int}, num_action::Vector{Int})
    row = Int[]
    col = Int[]
    for (t, cost) in enumerate(objective)
        if !isempty(cost.sparsity[1])
            shift = (t > 1 ? (sum(num_state[1:t-1]) + sum(num_action[1:t-1])) : 0)
            push!(row, (cost.sparsity[1] .+ shift)...)
            push!(col, (cost.sparsity[2] .+ shift)...)
        end
    end
    return collect(zip(row, col))
end

function hessian_indices(objective::Objective, key::Vector{Tuple{Int,Int}}, num_state::Vector{Int}, num_action::Vector{Int})
    indices = Vector{Int}[]
    for (t, cost) in enumerate(objective)
        if !isempty(cost.sparsity[1])
            row = Int[]
            col = Int[]
            shift = (t > 1 ? (sum(num_state[1:t-1]) + sum(num_action[1:t-1])) : 0)
            push!(row, (cost.sparsity[1] .+ shift)...)
            push!(col, (cost.sparsity[2] .+ shift)...)
            rc = collect(zip(row, col))
            push!(indices, [findfirst(x -> x == i, key) for i in rc])
        else
            push!(indices, Int[])
        end
    end
    return indices
end

# num_gradient(objective::Objective{T}) where T = sum([obj.num_gradient  for obj in objective])
# num_hessian(objective::Objective{T}) where T  = sum([obj.num_hessian   for obj in objective])


text:
"""
Functions that help the automation of tests
"""

using Test
using Catlab.CategoricalAlgebra, Catlab.Present
using Catlab.CategoricalAlgebra.CSetDataStructures: struct_acset

# Functions to make schemas programmatically
############################################
xs(x::Int)::Symbol = Symbol("x$x")
xs(xx::AbstractVector{Int})::Vector{Symbol} = [Symbol("x$x") for x in xx]
es(x::Int)::Symbol = Symbol("e$x")
es(xx::AbstractVector{Int})::Vector{Symbol} = [Symbol("e$x") for x in xx]

"""
Create a CSet type specified by a graph
Vertices are x₁,x₂,..., edges are e₁, e₂,...
all edges are indexed
"""
function graph_to_cset(grph::StructACSet, name::Symbol)::Pair{Presentation,StructACSet}
  pres = Presentation(FreeSchema)
  xobs = [Ob(FreeSchema,xs(i)) for i in 1:nv(grph)]
  for x in xobs
    add_generator!(pres, x)
  end
  for (i,(src, tgt)) in enumerate(zip(grph[:src], grph[:tgt]))
    add_generator!(pres, Hom(es(i), xobs[src], xobs[tgt]))
  end

  expr = struct_acset(name, StructACSet, pres, index=es(1:ne(grph)))
  eval(expr)
  csettype = eval(name)
  return pres => Base.invokelatest(csettype)
end

  """Create n copies of a CSet based on a graph schema"""
function init_graphs(name::Symbol, schema::StructACSet, consts::Vector{Int},
             n::Int=2)::Pair{Presentation, Vector{StructACSet}}
  pres, cset = graph_to_cset(schema, name)
  for (i, con) in enumerate(consts)
    add_parts!(cset, Symbol("x$i"), con)
  end
  return pres => [deepcopy(cset) for _ in 1:n]
end

# Testing nauty.jl and CSetAutomorphisms on the same input
##########################################################

"""Confirm canonical hash tracks with whether two ACSets are iso"""
function test_iso(a::StructACSet,b::StructACSet, pres::Union{Nothing,Presentation}=nothing)::Test.Pass
  @test a != b  # confirm they're not literally equal
  eq = is_isomorphic(a,b)
  tst = (x,y) -> eq ? x==y : x!=y # hashes should be equal iff they're iso
  if !isnothing(pres)
    @test tst(canonical_hash(a,pres=pres), canonical_hash(b; pres=pres))
  end
  @test tst(canonical_hash(a), canonical_hash(b))
end

# Particular Schemas
####################
@present TheoryDecGraph(FreeSchema) begin
  E::Ob
  V::Ob
  src::Hom(E,V)
  tgt::Hom(E,V)

  X::AttrType
  dec::Attr(E,X)
end

@acset_type Labeled(TheoryDecGraph)

ext: jl


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    struct PartitionedSampling <: AbstractSamplingAlgorithm

*Experimental feature, not part of stable public API.*

A sampling algorithm that partitions parameter space into multiple subspaces and
samples/integrates them independently ([Caldwell et al.](https://arxiv.org/abs/2008.03098)).

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)
"""
@with_kw struct PartitionedSampling{
    TR<:AbstractDensityTransformTarget,
    S<:AbstractSamplingAlgorithm,
    E<:AbstractSamplingAlgorithm,
    I<:IntegrationAlgorithm,
    P<:SpacePartitioningAlgorithm
} <: AbstractSamplingAlgorithm
    trafo::TR = PriorToUniform()
    npartitions::Integer = 10
    sampler::S = MCMCSampling(strict=false)
    exploration_sampler::E = MCMCSampling(nchains=30, nsteps = 800, strict=false)
    partitioner::P = KDTreePartitioning()
    integrator::I = AHMIntegration()
    nmax_resampling::AbstractFloat = 5
end

export PartitionedSampling

function bat_sample_impl(rng::AbstractRNG, target::PosteriorDensity, algorithm::PartitionedSampling)

    density_notrafo = convert(AbstractDensity, target)
    shaped_density, trafo = bat_transform(algorithm.trafo, density_notrafo)
    density = unshaped(shaped_density)

    @info "Generating Exploration Samples"
    exploration_samples = bat_sample(density, algorithm.exploration_sampler).result

    @info "Constructing Partition Tree"
    partition_tree, cost_values = partition_space(exploration_samples, algorithm.npartitions, algorithm.partitioner)
    # Convert 'partition_tree' structure into a set of truncated posteriors:
    posteriors_array = convert_to_posterior(density, partition_tree, extend_bounds = algorithm.partitioner.extend_bounds)

    @info "Sampling Subspaces"
    iterator_subspaces = [
        [subspace_ind, posteriors_array[subspace_ind], algorithm.sampler] for subspace_ind in Base.OneTo(algorithm.npartitions)]
    samples_subspaces_run = pmap(inp -> sample_subspace(inp...), iterator_subspaces)

    unconv_mask = [samples_subspace.isvalid for samples_subspace in samples_subspaces_run] # returns "false" if subspace was not converged during tuning cycle
    unconv_ind = findall(x->x==false, unconv_mask)
    rep_sspace = !isempty(unconv_ind) # perform resampling if "true"

    if algorithm.nmax_resampling > 0
        samples_subspaces = samples_subspaces_run[unconv_mask]
    else
        samples_subspaces = samples_subspaces_run
    end

    rec_level = 1
    resampled_trees = []

    while rep_sspace && rec_level <= algorithm.nmax_resampling

        @info "Re-Sampling Subspaces: #$rec_level"

        posteriors_array_run = []

        for (ind, rep_ind) in enumerate(unconv_ind)

            exploration_samples_rep = bat_sample(
                samples_subspaces_run[rep_ind].samples,
                OrderedResampling(nsamples=algorithm.exploration_sampler.nsteps)
            ).result

            partition_tree_rep, _ = partition_space(exploration_samples_rep, 2, algorithm.partitioner)
            push!(resampled_trees, partition_tree_rep)
            posteriors_rep_array = convert_to_posterior(posteriors_array[rep_ind], partition_tree_rep, extend_bounds = algorithm.partitioner.extend_bounds)

            append!(posteriors_array_run, posteriors_rep_array)

        end

        iterator_subspaces = [
            [subspace_ind, posteriors_array_run[subspace_ind], algorithm.sampler] for subspace_ind in Base.OneTo(length(posteriors_array_run))]
        samples_subspaces_run = pmap(inp -> sample_subspace(inp...), iterator_subspaces)

        unconv_mask = [samples_subspace.isvalid for samples_subspace in samples_subspaces_run]
        unconv_ind = findall(x->x==false, unconv_mask)
        rep_sspace = !isempty(unconv_ind)

        rec_level += 1
        if rep_sspace && (rec_level > algorithm.nmax_resampling)
            @warn "Convergence is not reached. Try to increase the number of resampling cycles."
            append!(samples_subspaces, samples_subspaces_run)
        else
            append!(samples_subspaces, samples_subspaces_run[unconv_mask])
        end

        posteriors_array = posteriors_array_run

    end

    @info "Integrating Subspaces"
    samples_subspaces = pmap(inp -> integrate_subspace(inp, algorithm.integrator), samples_subspaces)

    @info "Combining Samples"
    samples = deepcopy(samples_subspaces[1].samples)
    info = deepcopy(samples_subspaces[1].info)
    # Save indices from different subspaces:
    info.samples_ind[1] = 1:length(samples)
    for subspace in samples_subspaces[2:end]
        start_ind, stop_ind = length(samples)+1, length(samples)+length(subspace.samples)
        subspace.info.samples_ind[1] = start_ind:stop_ind
        append!(samples, subspace.samples)
        append!(info, subspace.info)
    end

    samples_trafo = varshape(shaped_density).(samples)
    samples_notrafo = inv(trafo).(samples_trafo)

    return (
        result = samples_notrafo, result_trafo = samples_trafo, trafo = trafo,
        info = info,
        exp_samples = exploration_samples, part_tree = partition_tree,
        cost_values = cost_values,
        resampled_trees = resampled_trees
    )
end

function sample_subspace(
    space_id::Integer,
    posterior::PosteriorDensity,
    sampling_algorithm::A,
) where {N<:NamedTuple, A<:AbstractSamplingAlgorithm, I<:IntegrationAlgorithm}

    @info "Sampling subspace #$space_id"
    sampling_wc_start = Dates.Time(Dates.now())
    sampling_cpu_time = CPUTime.@CPUelapsed begin
        sampling_output = bat_sample(posterior, sampling_algorithm)
        samples_subspace = sampling_output.result
        isvalid = sampling_output.isvalid
    end
    sampling_wc_stop = Dates.Time(Dates.now())

    info_subspace = TypedTables.Table(
            sampling_cpu_time = [sampling_cpu_time],
            sampling_wc = [Dates.value(sampling_wc_start):Dates.value(sampling_wc_stop)],
            worker_id = [Distributed.myid()],
            n_threads = [Threads.nthreads()],
            samples_ind = [0:0],
            sum_weights = [sum(samples_subspace.weight)],
        )

    return (samples = samples_subspace, info = info_subspace, isvalid=isvalid)
end

function integrate_subspace(
    sampling_reuslt::N,
    integration_algorithm::I,
) where {N<:NamedTuple, A<:AbstractSamplingAlgorithm, I<:IntegrationAlgorithm}

    samples_subspace = sampling_reuslt.samples

    integration_wc_start = Dates.Time(Dates.now())
    integration_cpu_time = CPUTime.@CPUelapsed begin
        # ToDo: Use samples_subspace_trafo for integration instead of samples_subspace?
        integras_subspace = bat_integrate(samples_subspace, integration_algorithm).result
    end
    integration_wc_stop = Dates.Time(Dates.now())

    samples_subspace_reweighted = DensitySampleVector(
        (
            samples_subspace.v,
            samples_subspace.logd,
            integras_subspace.val .* samples_subspace.weight ./ sum(samples_subspace.weight),
            samples_subspace.info,
            samples_subspace.aux
        )
    )

    info_subspace = TypedTables.Table(
            density_integral = [integras_subspace],
            integration_cpu_time = [integration_cpu_time],
            integration_wc = [Dates.value(integration_wc_start):Dates.value(integration_wc_stop)],
        )

    return (samples = samples_subspace_reweighted, info = TypedTables.Table(info_subspace, sampling_reuslt.info))
end

function convert_to_posterior_resampled(
    posterior::PosteriorDensity,
    partition_tree::SpacePartTree,
    posterior_subspace::PosteriorDensity;
    extend_bounds::Bool=true
)

    # Note: this function is very inefficient
    # To improve it, truncation of truncated densities should be defined, or use space transformaton

    if extend_bounds
        # Exploration samples might not always cover properly tails of the distribution.
        # We will extend boudnaries of the partition tree with original bounds which are:
        lo_bounds = getprior(posterior_subspace).bounds.vol.lo
        hi_bounds = getprior(posterior_subspace).bounds.vol.hi
        extend_tree_bounds!(partition_tree, lo_bounds, hi_bounds)
    end

    #Get flattened rectangular parameter bounds from tree
    subspaces_rect_bounds = get_tree_par_bounds(partition_tree)

    n_params = size(subspaces_rect_bounds[1])[1]

    posterior_array = map(
        x -> begin
            bounds = HyperRectBounds(x[:,1], x[:,2])
            prior_dist = BAT.truncate_density(getprior(posterior), bounds)
            PosteriorDensity(getlikelihood(posterior), prior_dist)
        end,
    subspaces_rect_bounds)

    return posterior_array
end

function convert_to_posterior(posterior::PosteriorDensity, partition_tree::SpacePartTree; extend_bounds::Bool=true)

    if extend_bounds
        # Exploration samples might not always cover properly tails of the distribution.
        # We will extend boudnaries of the partition tree with original bounds which are:
        vol = spatialvolume(var_bounds(posterior))
        lo_bounds = vol.lo
        hi_bounds = vol.hi
        extend_tree_bounds!(partition_tree, lo_bounds, hi_bounds)
    end

    #Get flattened rectangular parameter bounds from tree
    subspaces_rect_bounds = get_tree_par_bounds(partition_tree)

    posterior_array = map(subspaces_rect_bounds) do x
        bounds = StructArray{Interval}((x[:,1], x[:,2]))
        BAT.truncate_density(posterior, bounds)
    end

    return posterior_array
end


text:
using Parameters

# defining a simple test search problem
struct GridPosition
    x_idx::Int
    y_idx::Int
end

struct GridAction
    dx::Int
    dy::Int
end

apply_action(s::GridPosition, a::GridAction) = GridPosition(s.x_idx + a.dx, s.y_idx + a.dy)

@with_kw mutable struct GridNavigationProblem <: SearchProblem{GridPosition, GridAction}
    grid_dimensions::Tuple{Int, Int} = (10, 10)
    obstacles::Vector{GridPosition} = []
    aspace::Vector{GridAction} = [GridAction(0, 0),
                                  GridAction(1, 0), GridAction(-1, 0),
                                  GridAction(0, 1), GridAction(0, -1)]
    _n_expanded::Int = 0
end

GraphSearchZero.start_state(p::GridNavigationProblem) = GridPosition(1, 1)
GraphSearchZero.is_goal_state(p::GridNavigationProblem, s::GridPosition) = s == GridPosition(p.grid_dimensions...)
on_grid(p::GridNavigationProblem, s::GridPosition) = (1 <= s.x_idx <= p.grid_dimensions[1]) && (1 <= s.y_idx <= p.grid_dimensions[2])

function GraphSearchZero.successors(p::GridNavigationProblem, s::GridPosition)
    p._n_expanded += 1
    successors::Vector{Tuple{GridPosition, GridAction, Int}} = []
    sizehint!(successors, length(p.aspace))
    for a in p.aspace
        sp = apply_action(s, a)
        if !on_grid(p, sp) || sp in p.obstacles
            continue
        end
        # we are solving a minimum time problem, the step cost is always 1
        push!(successors, (sp, a, 1))
    end
    return successors
end

@testset "SearchTest" begin
    # tuples of (dimensions, obstacles, optimal_nsteps, solvable, n_expanded)
    test_setups::Vector{Tuple{Tuple{Int, Int}, Vector, Union{Int, Nothing}, Bool, Union{Int, Nothing}}} =
    [
     # 1x1 grid, no action needed to reach the goal
     ((1, 1), [], 0, true, 0),
     # empty 10x10 grid, solvable in 18 steps
     ((10, 10), [], 18, true, 18),
     # 10x10 grid with wall in the middle that has a gap, solvable in 18 steps
     ((10, 10), [GridPosition(5, y) for y in 1:9], 18, true, 18),
     # 10x10 grid with wall all the way, not solvable
     ((10, 10), [GridPosition(5, y) for y in 1:10], nothing, false, 40),
     # 5x5 grid with two walls, gap at top and bottom, solvable in 16 steps
     ((5, 5), [[GridPosition(2, y) for y in 1:4]...,
               [GridPosition(4, y) for y in 2:5]...], 16, true, 16)
    ]

    for (test_dims, test_obstacles, optimal_nsteps, solvable, n_expanded) in test_setups

        p = GridNavigationProblem(grid_dimensions=test_dims,
                                  obstacles=test_obstacles)
        # using the manhattan distance as a heuristic
        h = (s::GridPosition) -> abs(s.x_idx - p.grid_dimensions[1]) + abs(s.y_idx - p.grid_dimensions[2])
        eps = 0.0001
        if solvable
            aseq, sseq = @inferred weighted_astar_search(p, h, eps)
            @test length(aseq) == optimal_nsteps
        else
            @test_throws InfeasibleSearchProblemError weighted_astar_search(p, h, eps)
        end
        @test p._n_expanded == n_expanded

        # also run type inference tests:
        if solvable
            @test @testblock quote
                pp = deepcopy(p)
                a = rand(pp.aspace)
                n = GraphSearchZero.root_node(pp)
                @inferred GraphSearchZero.expand(n, pp)
            end
        end
    end
end


text:
export  calculateMultiIndices,
        findUnivariateIndices

function calculateMultiIndices(d::Int, n::Int)
    # d denotes dimension of random variables/number of sources of uncertainty,
    # n the maximum degree of multivariate basis
    # function to calculate indices of multivariate basis following the algorithm
    # from "Spectral Methods for Uncertainty Quantiﬁcation; Le Maitre, Knio;2014" p.516-517
    # No::BigInt = factorial(BigInt(d+n))/(factorial(Bigint(d))*factorial(BigInt(n)));  #number of polynomials of multivariate basis
  n < 0 && throw(DomainError(n, "maximum degree must be non-negative"))
  d <= 0 && throw(DomainError(d, "number of uncertainties must be positive"))
  # catch case n == 0 --> No-d==0
  n == 0 && return zeros(Int64,1,d)
  # non-pathological cases begin here
  No = numberPolynomials(d,n)
  inds = vcat(zeros(Int64,1,d),Matrix(1I,d,d),zeros(Int64, No-d-1, d));  #initiate index matrix for basis
  pi=ones(Int64,No,d);

  for k=2:No
    g = 0;
    for l=1:d
      pi[k,l]=sum(pi[k-1,:])-g;
      g=g+pi[k-1,l];
    end
  end

  P=d+1;
  for k=2:n
    L = P;
    for j=1:d, m=L-pi[k,j]+1:L
        P=P+1;
        inds[P,:]=inds[m,:];
        inds[P,j]=inds[P,j]+1;
    end
  end

  return inds
end

"""
    computes the number of polynomials with a multivariate basis
    `(d+n)!/(d!+n!)`
"""

function numberPolynomials(d::Int64, n::Int64)
    x, y = max(d,n), min(d,n)
    return UInt128(prod(UInt128(x+1):UInt128(d+n)) ÷ factorial(UInt128(y)))
end

"""
    findUnivariateIndices(i::Int64,ind::Matrix{Int64})::Vector{Int64}
Given the multi-index `ind` this function returns all entries of the multivariate basis
that correspond to the `i`th univariate basis.
"""
function findUnivariateIndices(i::Int64,ind::Matrix{Int64})::Vector{Int64}
  l,p = size(ind)
  @assert i<=p "basis is $p-variate, you requested $i-variate"
  deg = ind[end,end]
  @assert deg>=0 "invalid degree"
  col = ind[:,i]
  myind = zeros(Int64,deg)
  for deg_ = 1:deg
      myind[deg_] = findfirst(x->x==deg_,col)
  end
  pushfirst!(myind,1)
end

#################################################################
# function calculateMultiIndices_interaction(nξ::Int,deg::Int,j::Int,p::Int)
#     inds = calculateMultiIndices(nξ,deg)
#     get_interaction(inds,j,p)
# end

# function get_interaction(inds::Matrix{Int},j::Int,p::Int)
#     j < 0 && throw(error("interaction order must be non-negative"))
#     j > size(inds,2) && throw(error("interaction order cannot be greater than number of uncertainties"))
#     Iterators.filter(x -> count(!iszero,x) == j && sum(x) == p, eachrow(inds))
#     # alternative -> collect
# end


text:
# CXXScopeSpec
CXXScopeSpec() = CXXScopeSpec(create_cxx_scope_spec())

"""
    create_cxx_scope_spec() -> CXCXXScopeSpec
Return a pointer to a `clang::CXXScopeSpec` object.
"""
function create_cxx_scope_spec()
    status = Ref{CXInit_Error}(CXInit_NoError)
    ss = clang_CXXScopeSpec_create(status)
    @assert status[] == CXInit_NoError
    return ss
end

dispose(x::CXXScopeSpec) = clang_CXXScopeSpec_dispose(x)

function clear(x::CXXScopeSpec)
    @check_ptrs x
    return clang_CXXScopeSpec_clear(x)
end

function getScopeRep(x::CXXScopeSpec)
    @check_ptrs x
    return NestedNameSpecifier(clang_CXXScopeSpec_getScopeRep(x))
end

function getBeginLoc(x::CXXScopeSpec)
    @check_ptrs x
    return SourceLocation(clang_CXXScopeSpec_getBeginLoc(x))
end

function getEndLoc(x::CXXScopeSpec)
    @check_ptrs x
    return SourceLocation(clang_CXXScopeSpec_getEndLoc(x))
end

function setBeginLoc(x::CXXScopeSpec, loc::SourceLocation)
    @check_ptrs x
    return clang_CXXScopeSpec_setBeginLoc(x, loc)
end

function setEndLoc(x::CXXScopeSpec, loc::SourceLocation)
    @check_ptrs x
    return clang_CXXScopeSpec_setEndLoc(x, loc)
end

function isEmpty(x::CXXScopeSpec)
    @check_ptrs x
    return clang_CXXScopeSpec_isEmpty(x)
end

function isNotEmpty(x::CXXScopeSpec)
    @check_ptrs x
    return clang_CXXScopeSpec_isNotEmpty(x)
end

function isInvalid(x::CXXScopeSpec)
    @check_ptrs x
    return clang_CXXScopeSpec_isInvalid(x)
end

function isValid(x::CXXScopeSpec)
    @check_ptrs x
    return clang_CXXScopeSpec_isValid(x)
end

function setRange(x::CXXScopeSpec, r::SourceRange)
    setBeginLoc(x, getBeginLoc(r))
    setEndLoc(x, getEndLoc(r))
    return nothing
end

getRange(x::CXXScopeSpec, r::SourceRange) = SourceRange(getBeginLoc(x), getEndLoc(x))


text:
function calc_drag(buttons, drag, indrag, tracked_mousebutton)
    # only track if still the same button is pressed
    if length(buttons) == 1 && (!indrag[] || tracked_mousebutton[] == first(buttons))
        if !indrag[]
            tracked_mousebutton[] = first(buttons); indrag[] = true
            drag[] = Mouse.down # just started, so dragging is still false
            return drag[]
        else
            drag[] = Mouse.pressed
            return drag[]
        end
    end
    # already on notpressed, no need for update
    if drag[] != Mouse.notpressed
        drag[] = indrag[] ? Mouse.up : Mouse.notpressed
    end
    indrag[] = false
    return drag[]
end

function mousedrag(scene::Scene, native_window)
    indrag = RefValue(false)
    tracked_mousebutton = RefValue(Mouse.left)
    drag = RefValue(Mouse.notpressed)
    events = scene.events
    onany(events.mouseposition, events.mousebuttons) do mp, buttons
        d = calc_drag(buttons, drag, indrag, tracked_mousebutton)
        if (d == Mouse.pressed) || (d != events.mousedrag[])
            events.mousedrag[] = d
        end
        return
    end
    return
end

function disconnect!(window::AbstractScreen, signal)
    disconnect!(to_native(window), signal)
end
window_area(scene, native_window) = not_implemented_for(native_window)
window_open(scene, native_window) = not_implemented_for(native_window)
mouse_buttons(scene, native_window) = not_implemented_for(native_window)
mouse_position(scene, native_window) = not_implemented_for(native_window)
mousedrag(scene, native_window) = not_implemented_for(native_window)
scroll(scene, native_window) = not_implemented_for(native_window)
keyboard_buttons(scene, native_window) = not_implemented_for(native_window)
unicode_input(scene, native_window) = not_implemented_for(native_window)
dropped_files(scene, native_window) = not_implemented_for(native_window)
hasfocus(scene, native_window) = not_implemented_for(native_window)
entered_window(scene, native_window) = not_implemented_for(native_window)

function register_callbacks(scene::Scene, native_window)

    window_area(scene, native_window)
    window_open(scene, native_window)
    mouse_buttons(scene, native_window)
    mouse_position(scene, native_window)
    mousedrag(scene, native_window)
    scroll(scene, native_window)
    keyboard_buttons(scene, native_window)
    unicode_input(scene, native_window)
    dropped_files(scene, native_window)
    hasfocus(scene, native_window)
    entered_window(scene, native_window)

end

button_key(x::Type{T}) where {T} = error("Must be a keyboard or mouse button. Found: $T")
button_key(x::Type{Keyboard.Button}) = :keyboardbuttons
button_key(x::Type{Mouse.Button}) = :mousebuttons
button_key(x::Set{T}) where {T} = button_key(T)
button_key(x::T) where {T} = button_key(T)

"""
    ispressed(scene, buttons)

Returns true if all `buttons` are pressed in the given `scene`. `buttons` can be
a `Vector` or `Tuple` of `Keyboard` buttons (e.g. `Keyboard.a`), `Mouse` buttons
(e.g. `Mouse.left`) and `nothing`.
"""
function ispressed(scene::SceneLike, button::Union{Vector, Tuple})
    all(x-> ispressed(scene, x), button)
end

# TODO this is a bit shady, but maybe a nice api!
# So you can use void whenever you don't care what is pressed
ispressed(scene::SceneLike, ::Nothing) = true

function ispressed(buttons::Set{T}, button::T) where T <: Union{Keyboard.Button, Mouse.Button}
    return button in buttons
end

function ispressed(buttons::Set{T}, button::Set{T}) where T <: Union{Keyboard.Button, Mouse.Button}
    return issubset(button, buttons)
end

"""
    ispressed(scene, button)

Returns true if `button` is pressed in the given `scene`. The `button` can be
a `Keyboard` button (e.g. `Keyboard.a`), a `Mouse` button (e.g. `Mouse.left`)
or `nothing`. In the latter case `true` is always returned.
"""
function ispressed(scene::SceneLike, button)
    buttons = getfield(events(scene), button_key(button))[]
    ispressed(buttons, button)
end

"""
Picks a mouse position.  Implemented by the backend.
"""
function pick end

"""
    onpick(func, plot)
Calls `func` if one clicks on `plot`.  Implemented by the backend.
"""
function onpick end


text:
#############################################################################
# SCS.jl
# Wrapper around the SCS solver https://github.com/cvxgrp/scs
# See http://github.com/JuliaOpt/SCS.jl
#############################################################################
# test/options.jl
# Tests the ability to pass options
#############################################################################

using Base.Test
using MathProgBase.SolverInterface
using SCS

# The normal test
A = [1.0 1.0 0.0 0.0 0.0;
     0.0 1.0 0.0 0.0 1.0;
     0.0 0.0 1.0 1.0 1.0]
collb = [0.0, 0.0, 0.0, 0.0, 0.0]
obj   = [3.0, 4.0, 4.0, 9.0, 5.0]
rowub = [ 5.0,  3.0,  9.0]
s = SCSSolver()
m = MathProgBase.ConicModel(s)
MathProgBase.loadproblem!(m, -obj, A, rowub, [(:NonNeg,1:3)],[(:NonNeg,1:5)])
MathProgBase.optimize!(m)
@test_approx_eq_eps MathProgBase.getobjval(m) -99.0 1e-3

# With eps = 1e-8, solution should be far more accurate
s = SCSSolver(eps=1e-8)
m = MathProgBase.ConicModel(s)
MathProgBase.loadproblem!(m, -obj, A, rowub, [(:NonNeg,1:3)],[(:NonNeg,1:5)])
MathProgBase.optimize!(m)
@test_approx_eq_eps MathProgBase.getobjval(m) -99.0 1e-5

# With a warmstart from the eps = 1e-8 solution, solution should be extremely accurate even after 1 iteration
push!(m.options, (:warm_start, true))
push!(m.options, (:max_iters, 1))
MathProgBase.optimize!(m)
@test_approx_eq_eps MathProgBase.getobjval(m) -99.0 1e-5

# Now let's do the same warmstart, but on a new instance of the same problem
primal_sol = m.primal_sol
dual_sol = m.dual_sol
slack = m.slack
s = SCSSolver(max_iters=1)
m = MathProgBase.ConicModel(s)
MathProgBase.loadproblem!(m, -obj, A, rowub, [(:NonNeg,1:3)],[(:NonNeg,1:5)])
MathProgBase.setwarmstart!(m, primal_sol; dual_sol = dual_sol, slack = slack)
MathProgBase.optimize!(m)
@test_approx_eq_eps MathProgBase.getobjval(m) -99.0 1e-5


text:
# -----------------------------------------------------------------------------------------
##############
## Basic (Behavioural tests)
##############

# ----------------------------------------------------------------------------------------------------------------------
###############
## Detailed Tests
###############

@testset "coerce arrays" begin
    A = rand(Int, 2, 3)
    z = rand(Char, 2, 3)
    y = Any[1.0 2; 3 4]
    @test scitype_union(coerce(A, Continuous)) == Continuous
    @test scitype_union(coerce(A, OrderedFactor)) <: OrderedFactor
    @test scitype_union(coerce(z, Multiclass)) <: Multiclass
    @test scitype_union(coerce(y, Count)) === Count

    # test fix for issue 39
    y = collect(Int64, 1:5)
    @test_throws ScientificTypes.CoercionError coerce(y, Float64)
    @test_throws ScientificTypes.CoercionError coerce(y, Textual)

    @testset "coersion of Real->OrderedFactor" begin
        v = [0.1, 0.2, 0.2, 0.3, missing, 0.1]
        w = [0.1, 0.2, 0.2, 0.3, 0.1]
        @test_logs((:info, r"Trying to coerce from `Union{Missing,"),
                    global cv = coerce(v, OrderedFactor))
        cw = coerce(w, OrderedFactor)
        @test all(skipmissing(unique(cv)) .== [0.1, 0.2, 0.3])
        @test all(unique(cw) .== [0.1, 0.2, 0.3])
    end

    @testset "coersion of Any->MultiClass" begin
        v2 = Symbol.(collect("asdfghj"))
        global v2c = @test_logs((:warn, r"Converting array"),
                                coerce(v2, Multiclass))
        @test scitype_union(v2c) == Multiclass{7}
        @test eltype(v2c) <: Union{Missing,CategoricalValue{String}}

        v3 = Any[1, 2, 3]
        coerce(v3, Multiclass) # should be no warning here

        # normal behaviour
        v1 = categorical([1,2,1,2,1,2,missing])
        v2 = collect("aksldjfalsdjkfslkjdfalksjdf")
        @test_logs((:info, r"Trying to coerce from `Union"),
                global v1c = coerce(v1, Multiclass))
        global v2c # otherwise julia complains...
        v2c = coerce(v2, Multiclass)
        @test scitype_union(v1c) == Union{Missing,Multiclass{2}}
        @test scitype_union(v2c) == Multiclass{7}
        @test eltype(v1c) <: Union{Missing,CategoricalValue{Int64}}
        @test eltype(v2c) <: CategoricalValue{Char}
    end

    @testset "coerce Categorical->Continuous" begin
        a = categorical(["a","b","a","b",missing])
        a1 = coerce(a, Union{Count,Missing})
        @test scitype_union(a1) == Union{Missing,Count}
        @test all(skipmissing(a1 .== [1, 2, 1, 2, missing]))
        a1 = coerce(a, Union{Continuous,Missing})
        @test scitype_union(a1) == Union{Missing,Continuous}
        @test all(skipmissing(a1 .== [1., 2., 1., 2., missing]))

        y = categorical(1:10, ordered=true)
        new_order = [4, 10, 9, 7, 6, 2, 8, 3, 1, 5]
        levels!(y, new_order)
        @test all(coerce(y, Count) .== sortperm(new_order))
        @test all(coerce(y, Count) .== [9, 6, 8, 1, 10, 5, 4, 7, 3, 2])

        y = categorical([1:10..., missing, 11], ordered=true)
        @test all(skipmissing(coerce(y, Union{Continuous, Missing}) .==
                            float([1:10...,missing,11])))
    end

    # see issue 9
    @testset "coerce Text=>Num" begin
        x = ["52", "15", "125"]
        cx = coerce(x, Count)
        @test cx == [52, 15, 125]
        cx = coerce(x, Continuous)
        @test cx == [52.0, 15.0, 125.0]
        x = ["52", "15", "125", missing]
        cx = coerce(x, Union{Count,Missing})
        @test all(cx .=== [52, 15, 125, missing])
        cx = coerce(x, Union{Missing,Continuous})
        @test all(cx .=== [52.0, 15.0, 125.0, missing])
        @test_logs (:info, "Trying to coerce from `Union{Missing, String}` to `Count`.\nCoerced to `Union{Missing,Count}` instead.") coerce(x, Count)
    end
end

@testset "coersion of Images" begin
    #1. Single GrayImage:
    img = rand(10, 10)
    @test coerce(img, GrayImage) |> size == size(img)
    @test scitype(coerce(img, GrayImage)) == GrayImage{10, 10}

    #2. Single colorimage
    img = rand(10, 10, 3)
    @test coerce(img, ColorImage) |> size == (10, 10)
    @test scitype(coerce(img, ColorImage)) == ColorImage{10, 10}

    #3. Collection: 3dcollection -> GrayImage
    imgs = rand(10, 10, 3)
    @test coerce(imgs, GrayImage) |> size == (3,)
    @test scitype(coerce(imgs, GrayImage)) == AbstractArray{GrayImage{10, 10},1}

    #4. Collection: 4dcollection -> GrayImage
    imgs = rand(10, 10, 1, 3)
    @test coerce(imgs, GrayImage) |> size == (3,)
    @test scitype(coerce(imgs, GrayImage)) == AbstractArray{GrayImage{10, 10},1}

    #5. Collection : 4dcollection -> ColorImage
    imgs = rand(10, 10, 3, 5)
    @test coerce(imgs, ColorImage) |> size == (5,)
    @test scitype(coerce(imgs, ColorImage)) == AbstractArray{ColorImage{10,10},1}
end

# TODO: CSV as testing dep is causing issues for julia 1.0. Perhaps
# the followinng can be re-instated when we have a new LTS release:

# if VERSION ≥ v"1.3.0-"
#    @testset "Coerce Col2" begin
#       X = Tables.table(ones(1_000, 2))
#       tmp = tempname()
#       CSV.write(tmp, X)
#       data = CSV.read(tmp, DataFrame)
#       # data.Column1 and data.Column2 are Column2 (as of CSV 5.19)
#       @test data.Column1 isa AbstractArray{<:AbstractFloat}
#       dc = coerce(data, autotype(data, :discrete_to_continuous))
#       @test scitype(dc) == Table{AbstractArray{Continuous,1}}
#       rm(tmp)
#    end
# end

@testset "misc" begin
    x = Any['a', 5]
    @test (@test_logs (:info, "Char value encountered, such value will be coerced according to the corresponding numeric value (e.g. 'A' to 65).") coerce(x, Count)) == [97, 5]
    x = categorical(['a','b','a','b'])
    @test coerce(x, Continuous) == [1.0,2.0,1.0,2.0]
    y = [missing, 1, 2]
    x = y[2:end]
    c = ScientificTypes._tighten_if_needed(x, eltype(x), true)
    @test c == [1, 2]
    @test eltype(c) == Int
    y = categorical([missing,1,2])
    x = y[2:end]
    @test eltype(x) >: Missing
    c = ScientificTypes._tighten_if_needed(x, eltype(x), true)
    @test c == [1,2]
    @test eltype(c) == Int
    c = coerce(x, OrderedFactor, tight=true)
    @test c == categorical([1,2])
    @test !(eltype(c) >: Missing)
end


text:
baremodule CityKind
   const Small = 0
   const Medium = 1
   const Large = 2
end

baremodule AreaKind
   const Open = 0
   const SubUrban =  1
   const Urban = 2
end

baremodule TerrainKind
  const A = 0  # hilly terrain with moderate-to-heavy tree densities
  const B = 1  # hilly environment but rare vegetation, or high vegetation but flat terrain
  const C = 2  # mostly flat terrain with light tree densities
end

baremodule LeeArea
  const FreeSpace = (2.0,45.0)
  const OpenArea = (4.35,49.0)
  const SubUrban= (3.84,61.7)
  const Philadelphia = (3.68, 70.0)
  const Newark = (4.31, 64.0)
  const Tokyo = (3.05, 84.0)
  const NewYorkCity = (4.8,77.0)
end


text:
function load_synonym_table(db::SQLite.DB; logger::Union{Nothing,SimpleLogger} = nothing)::VLResult

    # initialize -
    synonym_dictionary = Dict{String,String}()

    try

        # load the synonym table -
        sql_string = "SELECT * FROM LEXEME_SYNONYM_TABLE;"

        # execute the call -
        query_result = DBInterface.execute(db, sql_string)

        # turn the result into a table (and return for now so we can see what gets returned)
        df = DataFrame(query_result)

        # populate the synonym table -
        number_of_rows = size(df,1)
        for row_index = 1:number_of_rows
            key = df[row_index,:general_lexeme]
            value = df[row_index,:canonical_lexeme]
            synonym_dictionary[key] = value
        end

        # return -
        return VLResult(synonym_dictionary)
    catch error
        return VLResult(error)
    end
end

max_stars_repo_name: varnerlab/CellFreeModelGenerationKit

max_issues_repo_name: varnerlab/CellFreeModelGenerationKit

max_forks_repo_name: varnerlab/CellFreeModelGenerationKit


text:
#Define Simulation state that also contains the control variables using the DEDataVector
#interface from the DifferentialEquations package
import DifferentialEquations: ODEProblem

mutable struct SimState{Float64} <: DEDataVector{Float64}
    x::Vector{Float64} #State variables
    u::Vector{Float64} #Control Variables
    F::Type{<:AbstractFrame} #FrameType
end #SimState

SimState(x::Vector{Float64}, u::Vector{Float64}) = SimState(x, u, ICRF)

# Create New ODE interface that uses x_dot = f(x,u,t), instead of x_dot = f(x,t)
#Also, make sure that integration is done in ICRF reference frame
function ODEProblem(f::Function, x0::Vector{Float64}, u0::Vector{Float64}, t0::Float64, tend::Float64)
    y0 = SimState(x0,u0)
    return ODEProblem(f,y0,(t0,tend))
end


text:
module test_52_MonteCarloMeasurementsWithDistributions

using ModiaPlot

include("test_base/test_52_MonteCarloMeasurementsWithDistributions.jl")

end


text:
# This file is a part of ForwardDiffPullbacks.jl, licensed under the MIT License (MIT).

const ZeroLike = Union{ChainRulesCore.AbstractZero, Nothing}
const RealTangentLike = Union{Real, ZeroLike}
const RealDualLike = Union{Real, Nothing, Missing}


text:
# Find a "formula" for the `decide` function in `resolve!`
# Given the result of `matches` (`b` below) and the value of `static` (`a` below),
# the goal is to find the simplest formula involving few basic functions
# with the same "truth table" as `r` below, using three-valued logic

r(a, b) = # a: static, b: match
    if a === true
        coalesce(b, false)
    elseif a === false
        b === missing
    else
        coalesce(b, true)
    end

# give a name to some basic functions:
arg1(a, b) = a
arg2(a, b) = b
narg1(a, b) = !a
narg2(a, b) = !b
egal(a, b) = a === b # have to give a name, because of a Julia bug:
                     # https://github.com/JuliaLang/julia/issues/40612
negal(a, b) = a !== b
impl(a, b) = !a | b
equiv(a, b) = impl(a, b) & impl(b, a)
revcoalesce(a, b) = coalesce(b, a)

# composing directly these functions lead to too much compilation for brute force,
# so we use truth tables instead

struct TT
    t::Matrix{Int}
    n::Symbol # name
end

TT(n::Symbol = :?) = TT(zeros(Bool, 3, 3), n)

# init a truthtable from a binary function
function init!(tt::TT, f::Function)
    m = (true, false, missing)
    r(x) = x === missing ? 3 : x ? 1 : 2
    for a = 1:3, b = 1:3
        tt.t[a, b] = r(f(m[a], m[b]))
    end
    tt
end

# compose truthtables: tt(a, b) = o(i1(a, b), i2(a, b))
function init!(tt::TT, o::TT, i1::TT, i2::TT)
    for a=1:3, b=1:3
        tt.t[a, b] = o.t[i1.t[a, b], i2.t[a, b]]
    end
    tt
end

init!(tt::TT, src::TT) = (copy!(tt.t, src.t); tt)

TT(f::Function) = init!(TT(nameof(f)), f)
TT(o::TT, i1::TT, i2::TT) = init!(TT(), o, i1, i2)

Base.show(io::IO, tt::TT) = print(io, tt.n)

const R = TT(r)

# compute distance from R
function check(tt::TT)
    # TODO: add multi-args `count` to `Base` ?
    dst = 0
    for i in 1:9
        dst += tt.t[i] != R.t[i]
    end
    dst
end

function solve(simple=false)
    fns = TT.([arg1, arg2, narg1, narg2, coalesce, |, &, ==, !=,  egal, negal])
    !simple && push!(fns, TT(revcoalesce), TT(impl), TT(equiv))

    F = TT(); O2 = TT()

    c = 0
    iter = 0
    for fn in fns, o1 in fns, o2 in fns, i1 in fns, i2 in fns
        iter += 1
        init!(O2, o2, i1, i2)
        init!(F, fn, o1, O2)

        if check(F) == 0
            if simple
                issimple = false
                for arg in (:arg1, :arg2)
                    issimple |= any([fn, o1, i1, i2, o2]) do f
                        f.n == arg
                    end
                end
                issimple || continue
            end
            c += 1
            @show fn, o1, o2, i1, i2
        end
    end
    (niter = iter, nsolutions = c)
end

sol1(s, m) = s | m === coalesce(s == m, m)
sol5(s, m) = s | m === (coalesce(s, m) == m)
sol2(s, m) = m === coalesce(s, m) | !m
sol3(s, m) = m === coalesce(s, m) | (s == m)
sol4(s, m) = m === (coalesce(s, m) == s | m)

for sol in (sol1, sol2, sol3, sol4, sol5)
    t = TT(sol)
    @assert check(t) == 0
end

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
#creates a random partition of 1:S in n_folds folds
function create_folds(S::Int, n_folds::Int)
    samples_per_fold = S ÷ n_folds
    fold_samples = fill(samples_per_fold, n_folds-1) #temporary
    push!(fold_samples, S - samples_per_fold*(n_folds-1))
    idx_folds = randsplit(S, fold_samples...)
end

#creates train and test set where k is the fold to be used as test
function traintest_split(data, ind_folds::Vector, k::Int)
    ind_test = ind_folds[k]
    ind_train = vcat(ind_folds[1:end .!= k]...)
    data_test = subdata(data, ind_test)
    data_train = subdata(data, ind_train)
    return (data_train, ind_train), (data_test, ind_test)
end

#performs the entire crossvalidation and writes in the supplied variables the learning result
function crossvalidation!(data, model::Model, prior_var::Float64, nfolds::Int, training_set::Vector{Any}, test_set::Vector{Any}, model_parameters::Vector{Any},
     outcomes::Vector{Any}, L2reg::Vector{Any}, logP::Vector{Any}, obs::Vector{Any}, verbose::Bool=true; ftol_rel = 1e-7, Δλ = 1e-10, monitor = nothing)

     ind_folds = create_folds(data.S, nfolds)

    for k in 1:nfolds
        if verbose
            println("fold ", k, "/", nfolds, "    ************************")
        end

        #creating test set and training set
        (data_train, ind_train), (data_test, ind_test) = traintest_split(data, ind_folds, k)

        #setting starting point for learning
        starting_point!(model, data_train)

        #setting up regularization
        prior = set_prior(prior_var, model)

        #learning
        opt = set_optimizer(:LD_MMA, model, data_train, ftol_rel = ftol_rel, Δλ = Δλ)
        model, ret, loglike, prop = learn!(model, prior, opt, data_train; monitor = monitor)

        push!(training_set, (data_train, ind_train))
        push!(test_set, (data_test, ind_test))
        push!(model_parameters, model)
        push!(outcomes, ret)
        push!(L2reg, prior)
        push!(logP, loglike)
        push!(obs, prop)
    end
    return nothing
end

#computes the average curve and associate to it an error using the standard deviation of the mean
function avg_selectivity_curve(θcurves_train)
    avg_corr_train = Array{Float64, 2}(undef,size(θcurves_train[1]))

    for j in 1:size(avg_corr_train)[1]
        avg_corr_train[j,1] = θcurves_train[1][j,1]
        avg_corr_train[j,2] = mean([θcurves_train[k][j,2] for k in 1:length(θcurves_train)])
    end

    std_corr_train = Array{Float64, 1}(undef, size(avg_corr_train)[1])

    for j in 1:length(std_corr_train)
        std_corr_train[j] = std([θcurves_train[k][j,2] for k in 1:length(θcurves_train)])
    end
    return avg_corr_train, std_corr_train ./ sqrt(length(θcurves_train))
end

#computes energies and probabilities of a dataset using learnt parameters
function inference(model::Model, data)
    energies = [Energy(model.epistasis, model.x, s) for s in data.variants]
    p = [pb(model.epistasis, model.x, s) for s in data.variants]
    return (energies, p)
end

#computes energies and probabilities of a dataset using learnt parameters for all train-test instancies of the crossvalidation
function traintest_inference(model_parameters::Vector{Any}, data_train::Vector{Any}, data_test::Vector{Any})
    nfolds = length(model_parameters)
    etr = Vector{Vector{Float64}}(undef, nfolds)
    ptr = Vector{Vector{Float64}}(undef, nfolds)
    ete = Vector{Vector{Float64}}(undef, nfolds)
    pte = Vector{Vector{Float64}}(undef, nfolds)

    for k in 1:nfolds
    etr[k], ptr[k] = inference(model_parameters[k], data_train[k][1])
    ete[k], pte[k] = inference(model_parameters[k], data_test[k][1])
    end
    return (etr, ptr), (ete, pte)
end

#plots the probability distribution as function of energies for a dataset
function plot_distribution(energies::Vector{Float64}, ps::Vector{Float64})
    ind = sortperm(energies, rev=true)
    plot(energies[ind], ps[ind])
    xlabel("Es")
    ylabel("ps")
end

#plots Ns(t) vs Ns(t+1)/λ(t) to check for deterministic binding regime
function deterministic_binding(pb::Vector{Float64}, λ::Float64, data, t::Int; count_thr = 20)
    ind = findall(data.counts[:,1,t] .> count_thr);
    average = pb[ind] .* data.counts[ind,1,t]
    variance = pb[ind] .* (1.0 .- pb[ind]) .* data.counts[ind,1,t]
    line = LinRange(0.0, max(maximum(data.counts[ind,1,t+1] ./ λ), maximum(average)),10)
    errorbar(data.counts[ind,1,t+1] ./ λ, average, sqrt.(variance), ls="None", marker="s")
    corr11 = cor(data.counts[ind,1,t+1] ./ λ, average)

    plot(line,line, linewidth=1, color="red")
    title("(t=$t)   corr= $corr11")
    ylabel("N(s)ps")
    xlabel("Ns(t+1)/λ(t)")
end

#does the same of the previous function but plots the result in a logarithmic scale
function log_deterministic_binding(pb::Vector{Float64}, λ::Float64, data, t::Int; count_thr = 20)
    ind = findall(data.counts[:,1,t] .> count_thr);
    average = pb[ind] .* data.counts[ind,1,t]
    variance = pb[ind] .* (1.0 .- pb[ind]) .* data.counts[ind,1,t]
    line = LinRange(0.0, max(maximum(data.counts[ind,1,t+1] ./ λ), maximum(average)),10)
    xscale(:log)
    yscale(:log)
    errorbar(data.counts[ind,1,t+1] ./ λ, average, sqrt.(variance), ls="None", marker="s")
    corr11 = cor(data.counts[ind,1,t+1] ./ λ, average)

    plot(line,line, linewidth=1, color="red")
    title("(t=$t)   corr= $corr11")
    ylabel("N(s)ps")
    xlabel("Ns(t+1)/λ(t)")
end

export crossvalidation!, avg_selectivity_curve, inference, traintest_inference,
    plot_distribution, deterministic_binding, log_deterministic_deterministic


text:
module UnconstrainedProblems

### Sources
###
### [1] Ali, Khompatraporn, & Zabinsky: A Numerical Evaluation of Several Stochastic Algorithms on Selected Continuous Global Optimization Test
### Link: www.researchgate.net/profile/Montaz_Ali/publication/226654862_A_Numerical_Evaluation_of_Several_Stochastic_Algorithms_on_Selected_Continuous_Global_Optimization_Test_Problems/links/00b4952bef133a1a6b000000.pdf
###
### [2] Fletcher & Powell: A rapidly convergent descent method for minimization,

immutable OptimizationProblem
    name::AbstractString
    f::Function
    g!::Function
    h!::Function
    initial_x::Vector{Float64}
    solutions::Vector
    isdifferentiable::Bool
    istwicedifferentiable::Bool
end

examples = Dict{AbstractString, OptimizationProblem}()

##########################################################################
###
### Exponential Function
###
##########################################################################

function exponential(x::Vector)
    return exp((2.0 - x[1])^2) + exp((3.0 - x[2])^2)
end

function exponential_gradient!(x::Vector, storage::Vector)
    storage[1] = -2.0 * (2.0 - x[1]) * exp((2.0 - x[1])^2)
    storage[2] = -2.0 * (3.0 - x[2]) * exp((3.0 - x[2])^2)
end

function exponential_hessian!(x::Vector, storage::Matrix)
    storage[1, 1] = 2.0 * exp((2.0 - x[1])^2) * (2.0 * x[1]^2 - 8.0 * x[1] + 9)
    storage[1, 2] = 0.0
    storage[2, 1] = 0.0
    storage[2, 2] = 2.0 * exp((3.0 - x[1])^2) * (2.0 * x[2]^2 - 12.0 * x[2] + 19)
end

examples["Exponential"] = OptimizationProblem("Exponential",
                                              exponential,
                                              exponential_gradient!,
                                              exponential_hessian!,
                                              [0.0, 0.0],
                                              [2.0, 3.0],
                                              true,
                                              true)

##########################################################################
###
### Fletcher-Powell
###
### From [2]
### Source: A rapidly convergent descent method for minimization
###         Fletcher & Powell
##########################################################################

function fletcher_powell(x::Vector)
    function theta(x::Vector)
        if x[1] > 0
            return atan(x[2] / x[1]) / (2.0 * pi)
        else
            return (pi + atan(x[2] / x[1])) / (2.0 * pi)
        end
    end

    return 100.0 * (x[3] - 10.0 * theta(x))^2 +
            (sqrt(x[1]^2 + x[2]^2) - 1.0)^2 + x[3]^2
end

# TODO: Implement
function fletcher_powell_gradient!(x::Vector, storage::Vector)
    return
end

# TODO: Implement
function fletcher_powell_hessian!(x::Vector, storage::Matrix)
    return
end

examples["Fletcher-Powell"] = OptimizationProblem("Fletcher-Powell",
                                                  fletcher_powell,
                                                  fletcher_powell_gradient!,
                                                  fletcher_powell_hessian!,
                                                  [-1.0, 0.0, 0.0], # Same as in source
                                                  [1.0, 0.0, 0.0],
                                                  false,
                                                  false)

##########################################################################
###
### Himmelblau's Function
###
##########################################################################

function himmelblau(x::Vector)
    return (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
end

function himmelblau_gradient!(x::Vector, storage::Vector)
    storage[1] = 4.0 * x[1]^3 + 4.0 * x[1] * x[2] -
                  44.0 * x[1] + 2.0 * x[1] + 2.0 * x[2]^2 - 14.0
    storage[2] = 2.0 * x[1]^2 + 2.0 * x[2] - 22.0 +
                  4.0 * x[1] * x[2] + 4.0 * x[2]^3 - 28.0 * x[2]
end

function himmelblau_hessian!(x::Vector, storage::Matrix)
    storage[1, 1] = 12.0 * x[1]^2 + 4.0 * x[2] - 42.0
    storage[1, 2] = 4.0 * x[1] + 4.0 * x[2]
    storage[2, 1] = 4.0 * x[1] + 4.0 * x[2]
    storage[2, 2] = 12.0 * x[2]^2 + 4.0 * x[1] - 26.0
end

examples["Himmelblau"] = OptimizationProblem("Himmelblau",
                                             himmelblau,
                                             himmelblau_gradient!,
                                             himmelblau_hessian!,
                                             [2.0, 2.0],
                                             [3.0, 2.0],
                                             true,
                                             true)
##########################################################################
###
### Hosaki's Problem
###
### Problem 20 in [1]
##########################################################################

function hosaki(x::Vector)
    a = (1.0 - 8.0 * x[1] + 7.0 * x[1]^2 - (7.0 / 3.0) * x[1]^3 + (1.0 / 4.0) * x[1]^4)
    return a * x[2]^2 * exp(-x[2])
end

function hosaki_gradient!(x::Vector, storage::Vector)
    storage[1] = (x[1]^3 - 7.0 * x[1]^2 + 14.0 * x[1] - 8)* x[2]^2 * exp(-x[2])
    storage[2] = 2.0 * (1.0 - 8.0 * x[1] + 7.0 * x[1]^2 - (7.0 / 3.0) * x[1]^3 + (1.0 / 4.0) * x[1]^4) * x[2] * exp(-x[2]) - (1.0 - 8.0 * x[1] + 7.0 * x[1]^2 - (7.0 / 3.0) * x[1]^3 + (1.0 / 4.0) * x[1]^4) * x[2]^2 * exp(-x[2])
end

function hosaki_hessian!(x::Vector, storage::Matrix)
    storage[1, 1] = (3.0 * x[1]^2 - 14.0 * x[1] + 14.0) * x[2]^2 * exp(-x[2])
    storage[1, 2] = 2.0 * (x[1]^3 - 7.0 * x[1]^2 + 14.0 * x[1] - 8.0) * x[2] * exp(-x[2])  - (x[1]^3 - 7.0 * x[1]^2 + 14.0 * x[1] - 8.0) * x[2]^2 * exp(-x[2])
    storage[2, 1] =  2.0 * (x[1]^3 - 7.0 * x[1]^2 + 14.0 * x[1] - 8.0) * x[2] * exp(-x[2])  - (x[1]^3 - 7.0 * x[1]^2 + 14.0 * x[1] - 8.0) * x[2]^2 * exp(-x[2])
    storage[2, 2] = 2.0 * (1.0 - 8.0 * x[1] + 7.0 * x[1]^2 - (7.0 / 3.0) * x[1]^3 + (1.0 / 4.0) * x[1]^4) * exp(-x[2]) - 4.0 * ( 1.0 - 8.0 * x[1] + 7.0 *  x[1]^2 - (7.0 / 3.0) * x[1]^3 + (1.0 / 4.0) * x[1]^4) * x[2] * exp(-x[2]) + (1.0 - 8.0 * x[1] + 7.0 * x[1]^2 - (7.0 / 3.0) * x[1]^3 + (1.0 / 4.0) * x[1]^4) * x[2]^2 * exp(-x[2])
end

examples["Hosaki"] = OptimizationProblem("Hosaki",
                                         hosaki,
                                         hosaki_gradient!,
                                         hosaki_hessian!,
                                         [3.6, 1.9],
                                         [4.0, 2.0],
                                         true,
                                         true)

##########################################################################
###
### Large-Scale Quadratic
###
##########################################################################

function large_polynomial(x::Vector)
    res = zero(x[1])
    for i in 1:250
        res += (i - x[i])^2
    end
    return res
end

function large_polynomial_gradient!(x::Vector, storage::Vector)
    for i in 1:250
        storage[i] = -2.0 * (i - x[i])
    end
end

function large_polynomial_hessian!(x::Vector, storage::Matrix)
    for i in 1:250
        for j in i:250
            if i == j
                storage[i, j] = 2.0
            else
                storage[i, j] = 0.0
                storage[j, i] = 0.0
            end
        end
    end
end

examples["Large Polynomial"] = OptimizationProblem("Large Polynomial",
                                                   large_polynomial,
                                                   large_polynomial_gradient!,
                                                   large_polynomial_hessian!,
                                                   zeros(250),
                                                   float([1:250;]),
                                                   true,
                                                   true)

##########################################################################
###
### Parabola
###
##########################################################################

function parabola(x::Vector)
    return (1.0 - x[1])^2 + (2.0 - x[2])^2 + (3.0 - x[3])^2 +
            (5.0 - x[4])^2 + (8.0 - x[5])^2
end

function parabola_gradient!(x::Vector, storage::Vector)
    storage[1] = -2.0 * (1.0 - x[1])
    storage[2] = -2.0 * (2.0 - x[2])
    storage[3] = -2.0 * (3.0 - x[3])
    storage[4] = -2.0 * (5.0 - x[4])
    storage[5] = -2.0 * (8.0 - x[5])
end

function parabola_hessian!(x::Vector, storage::Matrix)
    for i in 1:5
        for j in 1:5
            if i == j
                storage[i, j] = 2.0
            else
                storage[i, j] = 0.0
            end
        end
    end
end

examples["Parabola"] = OptimizationProblem("Parabola",
                                           parabola,
                                           parabola_gradient!,
                                           parabola_hessian!,
                                           [0.0, 0.0, 0.0, 0.0, 0.0],
                                           [1.0, 2.0, 3.0, 5.0, 8.0],
                                           true,
                                           true)

##########################################################################
###
### Simple 4th-Degree Polynomial Example
###
##########################################################################

function polynomial(x::Vector)
    return (10.0 - x[1])^2 + (7.0 - x[2])^4 + (108.0 - x[3])^4
end

function polynomial_gradient!(x::Vector, storage::Vector)
    storage[1] = -2.0 * (10.0 - x[1])
    storage[2] = -4.0 * (7.0 - x[2])^3
    storage[3] = -4.0 * (108.0 - x[3])^3
end

function polynomial_hessian!(x::Vector, storage::Matrix)
    storage[1, 1] = 2.0
    storage[1, 2] = 0.0
    storage[1, 3] = 0.0
    storage[2, 1] = 0.0
    storage[2, 2] = 12.0 * (7.0 - x[2])^2
    storage[2, 3] = 0.0
    storage[3, 1] = 0.0
    storage[3, 2] = 0.0
    storage[3, 3] = 12.0 * (108.0 - x[3])^2
end

examples["Polynomial"] = OptimizationProblem("Polynomial",
                                             polynomial,
                                             polynomial_gradient!,
                                             polynomial_hessian!,
                                             [0.0, 0.0, 0.0],
                                             [10.0, 7.0, 108.0],
                                             true,
                                             true)

##########################################################################
###
### Powell (d=4)
###
### Problem 35 in [1]
### Difficult since the hessian is singular at the optimum
##########################################################################

function powell(x::Vector)
    return (x[1] + 10.0 * x[2])^2 + 5.0 * (x[3] - x[4])^2 +
            (x[2] - 2.0 * x[3])^4 + 10.0 * (x[1] - x[4])^4
end

function powell_gradient!(x::Vector, storage::Vector)
    storage[1] = 2.0 * (x[1] + 10.0 * x[2]) + 40.0 * (x[1] - x[4])^3
    storage[2] = 20.0 * (x[1] + 10.0 * x[2]) + 4.0 * (x[2] - 2.0 * x[3])^3
    storage[3] = 10.0 * (x[3] - x[4]) - 8.0 * (x[2] - 2.0 * x[3])^3
    storage[4] = -10.0 * (x[3] - x[4]) - 40.0 * (x[1] - x[4])^3
end

function powell_hessian!(x::Vector, storage::Matrix)
    storage[1, 1] = 2.0 + 120.0 * (x[1] - x[4])^2
    storage[1, 2] = 20.0
    storage[1, 3] = 0.0
    storage[1, 4] = -120.0 * (x[1] - x[4])^2
    storage[2, 1] = 20.0
    storage[2, 2] = 200.0 + 12.0 * (x[2] - 2.0 * x[3])^2
    storage[2, 3] = -24.0 * (x[2] - 2.0 * x[3])^2
    storage[2, 4] = 0.0
    storage[3, 1] = 0.0
    storage[3, 2] = -24.0 * (x[2] - 2.0 * x[3])^2
    storage[3, 3] = 10.0 + 48.0 * (x[2] - 2.0 * x[3])^2
    storage[3, 4] = -10.0
    storage[4, 1] = -120.0 * (x[1] - x[4])^2
    storage[4, 2] = 0.0
    storage[4, 3] = -10.0
    storage[4, 4] = 10.0 + 120.0 * (x[1] - x[4])^2
end

examples["Powell"] = OptimizationProblem("Powell",
                                         powell,
                                         powell_gradient!,
                                         powell_hessian!,
                                         [3.0, -1.0, 0.0, 1.0],
                                         [0.0, 0.0, 0.0, 0.0],
                                         true,
                                         true)

##########################################################################
###
### Rosenbrock (2D)
###
### Problem 38 in [1]
###
### Saddle point makes optimization difficult
##########################################################################

function rosenbrock(x::Vector)
    return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end

function rosenbrock_gradient!(x::Vector, storage::Vector)
    storage[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
    storage[2] = 200.0 * (x[2] - x[1]^2)
end

function rosenbrock_hessian!(x::Vector, storage::Matrix)
    storage[1, 1] = 2.0 - 400.0 * x[2] + 1200.0 * x[1]^2
    storage[1, 2] = -400.0 * x[1]
    storage[2, 1] = -400.0 * x[1]
    storage[2, 2] = 200.0
end

examples["Rosenbrock"] = OptimizationProblem("Rosenbrock",
                                             rosenbrock,
                                             rosenbrock_gradient!,
                                             rosenbrock_hessian!,
                                             [0.0, 0.0],
                                             [1.0, 1.0],
                                             true,
                                             true)

end # module


text:
# This file was generated, do not modify it. # hide
coerce!(df, :Class => OrderedFactor{2});

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using ITensors,
      Test

include("test_tagset.jl")
include("test_index.jl")
include("test_indexset.jl")
include("test_itensor.jl")
include("test_contract.jl")
include("test_trg.jl")
include("test_ctmrg.jl")
include("test_siteset.jl")
include("test_mps.jl")
include("test_mpo.jl")
include("test_autompo.jl")
include("test_svd.jl")


text:
module TestLeafArgument

using Comonicon.AST
using Comonicon.JuliaExpr
using Comonicon.JuliaExpr: emit, emit_body, emit_norm_body, emit_dash_body
using Test

function foo(a::Int)
    @test a == 3
end

cmd = Entry(;
    version = v"1.1.0",
    root = LeafCommand(; fn = :foo, name = "leaf", args = [Argument(; name = "a", type = Int)]),
)

eval(emit(cmd))

@testset "test leaf argument" begin
    @test command_main(["3"]) == 0
    @test command_main(["1.2"]) == 1
end

end

module Issue141

using Comonicon.AST
using Comonicon.JuliaExpr
using Comonicon.JuliaExpr: emit, emit_body, emit_norm_body, emit_dash_body
using Test

function main(xs...)
    @test xs == ("1", "2", "3")
end

cmd = Entry(
    LeafCommand(
        main,
        "main",
        Argument[],
        0,
        Argument("xs", Any, true, false, nothing, Description(nothing, nothing), LineNumberNode(0)),
        Dict{String,Flag}(),
        Dict{String,Option}(),
        Description("", ""),
        LineNumberNode(0),
    ),
    nothing,
    LineNumberNode(0),
)

eval(emit(cmd))

@testset "issue/#141" begin
    Issue141.command_main(["1", "2", "3"])
end

end


text:
function cell_connectivity(mesh::UniformMesh)
    ncells = number_of_elements(mesh)
    nfaces = faces_per_cell(mesh)
    connectivity = zeros(Int, nfaces, ncells)
    for cellid = 1:ncells
        connectivity[:, cellid] .= neighbors(mesh, cellid)
    end
    return connectivity
end

function dimension(mesh::UniformMesh)
    return CartesianMesh.dimension(mesh)
end

function construct_cell_maps(mesh::UniformMesh)
    ncells = number_of_elements(mesh)
    return [CellMap(element(mesh, i)...) for i = 1:ncells]
end

function reference_corner(mesh::UniformMesh)
    return mesh.x0
end

function reference_corner(mesh)
    return reference_corner(background_mesh(mesh))
end

function mesh_widths(mesh::UniformMesh)
    return mesh.widths
end

function elements_per_mesh_side(mesh::UniformMesh)
    return mesh.nelements
end

function mesh_widths(mesh)
    return mesh_widths(background_mesh(mesh))
end

function nodes_per_mesh_side(mesh)
    return nodes_per_mesh_side(background_mesh(mesh))
end

function cell_sign_to_row(s)
    (s == -1 || s == +1) ||
        error("Use ±1 to index into rows (i.e. phase), got index = $s")
    row = s == +1 ? 1 : 2
    return row
end

function levelset_normals(levelset, points, invjac)
    npts = size(points)[2]
    if npts == 0
        return []
    else
        g = hcat([gradient(levelset, points[:, i])' for i = 1:npts]...)
        normals = diagm(invjac) * g
        normalize_normals!(normals)
        return normals
    end
end

function normalize_normals!(normals)
    dim, npts = size(normals)
    for i = 1:npts
        n = normals[:, i]
        normals[:, i] .= n / norm(n)
    end
end

function points(quad::QuadratureRule)
    return quad.points
end

function weights(quad::QuadratureRule)
    return quad.weights
end

function opposite_face(faceid)
    if faceid == 1
        return 3
    elseif faceid == 2
        return 4
    elseif faceid == 3
        return 1
    elseif faceid == 4
        return 2
    else
        error("Expected faceid ∈ {1,2,3,4} got faceid = $faceid")
    end
end

function jacobian(mesh)
    return jacobian(background_mesh(mesh))
end

function inverse_jacobian(mesh)
    return inverse_jacobian(background_mesh(mesh))
end

function dimension(mesh)
    return dimension(background_mesh(mesh))
end

function face_determinant_jacobian(mesh)
    return face_determinant_jacobian(background_mesh(mesh))
end

function cell_connectivity(mesh, faceid, cellid)
    return cell_connectivity(background_mesh(mesh), faceid, cellid)
end

function reference_face_normals()
    n1 = [0.0, -1.0]
    n2 = [1.0, 0.0]
    n3 = [0.0, 1.0]
    n4 = [-1.0, 0.0]
    return [n1, n2, n3, n4]
end

function reference_bottom_face_midpoint()
    [0.0, -1.0]
end

function reference_right_face_midpoint()
    [1.0, 0.0]
end

function reference_top_face_midpoint()
    [0.0, 1.0]
end

function reference_left_face_midpoint()
    [-1.0, 0.0]
end

function reference_face_midpoints()
    [
        reference_bottom_face_midpoint(),
        reference_right_face_midpoint(),
        reference_top_face_midpoint(),
        reference_left_face_midpoint(),
    ]
end

function check_cellsign(cellsign)
    @assert cellsign == +1 || cellsign == -1 || cellsign == 0
end

function check_eta(eta)
    @assert eta == +1 || eta == -1 || eta == 0
end

function rotate_90(normals)
    rot = [
        0.0 -1.0
        1.0 0.0
    ]
    return rot * normals
end

function scale_area(tangents, invjac)
    den = sqrt.((tangents .^ 2)' * (invjac .^ 2))
    return 1.0 ./ den
end

function map_to_spatial(refpoints, refcellids, mesh)

    dim, numpts = size(refpoints)
    @assert length(refcellids) == numpts
    spatialpoints = zeros(dim, numpts)

    for idx = 1:numpts
        cellid = refcellids[idx]
        cellmap = cell_map(mesh, cellid)
        spatialpoints[:, idx] = cellmap(refpoints[:, idx])
    end
    return spatialpoints
end

function map_to_spatial_on_merged_mesh(
    refpoints,
    refcellids,
    levelsetsign,
    mesh,
)

    dim, numpts = size(refpoints)
    @assert length(refcellids) == numpts
    spatialpoints = zeros(dim, numpts)

    for idx = 1:numpts
        cellid = refcellids[idx]
        cellmap = cell_map(mesh, levelsetsign, cellid)
        spatialpoints[:, idx] = cellmap(refpoints[:, idx])
    end
    return spatialpoints
end

function map_to_reference(spatialpoints, cellids, mesh)
    referencepoints = similar(spatialpoints)
    for (idx, cellid) in enumerate(cellids)
        cellmap = cell_map(mesh, cellid)
        referencepoints[:, idx] = inverse(cellmap, spatialpoints[:, idx])
    end
    return referencepoints
end

function map_to_reference_on_merged_mesh(
    spatialpoints,
    cellids,
    levelsetsign,
    mesh,
)
    referencepoints = similar(spatialpoints)
    for (idx, cellid) in enumerate(cellids)
        cellmap = cell_map(mesh, levelsetsign, cellid)
        referencepoints[:, idx] = inverse(cellmap, spatialpoints[:, idx])
    end
    return referencepoints
end

function collect_normals_at_spatial_points(spatialpoints, cellids, levelset)
    dim, npts = size(spatialpoints)
    @assert length(cellids) == npts

    normals = zeros(dim, npts)
    mesh = background_mesh(levelset)
    invjac = inverse_jacobian(mesh)

    for (idx, cellid) in enumerate(cellids)
        load_coefficients!(levelset, cellid)
        cellmap = cell_map(mesh, cellid)
        refpoint = inverse(cellmap, spatialpoints[:, idx])

        normals[:, idx] = gradient(interpolater(levelset), refpoint)
    end
    normals = diagm(invjac) * normals
    normalize_normals!(normals)
    return normals
end

function hermite_dofs_per_node(dim)
    @assert dim == 2
    return 4
end

function hermite_interpolation_indices(dim)
    @assert dim == 2
    return [1,5,9,13]
end


text:
# unary
"""
    pyneg(x)

Equivalent to `-x` in Python.
"""
pyneg(x) = pynew(errcheck(@autopy x C.PyNumber_Negative(getptr(x_))))
"""
    pypos(x)

Equivalent to `+x` in Python.
"""
pypos(x) = pynew(errcheck(@autopy x C.PyNumber_Positive(getptr(x_))))
"""
    pyabs(x)

Equivalent to `abs(x)` in Python.
"""
pyabs(x) = pynew(errcheck(@autopy x C.PyNumber_Absolute(getptr(x_))))
"""
    pyinv(x)

Equivalent to `~x` in Python.
"""
pyinv(x) = pynew(errcheck(@autopy x C.PyNumber_Invert(getptr(x_))))
"""
    pyindex(x)

Convert `x` losslessly to an `int`.
"""
pyindex(x) = pynew(errcheck(@autopy x C.PyNumber_Index(getptr(x_))))
export pyneg, pypos, pyabs, pyinv, pyindex

# binary
"""
    pyadd(x, y)

Equivalent to `x + y` in Python.
"""
pyadd(x, y) = pynew(errcheck(@autopy x y C.PyNumber_Add(getptr(x_), getptr(y_))))
"""
    pysub(x, y)

Equivalent to `x - y` in Python.
"""
pysub(x, y) = pynew(errcheck(@autopy x y C.PyNumber_Subtract(getptr(x_), getptr(y_))))
"""
    pymul(x, y)

Equivalent to `x * y` in Python.
"""
pymul(x, y) = pynew(errcheck(@autopy x y C.PyNumber_Multiply(getptr(x_), getptr(y_))))
"""
    pymatmul(x, y)

Equivalent to `x @ y` in Python.
"""
pymatmul(x, y) = pynew(errcheck(@autopy x y C.PyNumber_MatrixMultiply(getptr(x_), getptr(y_))))
"""
    pyfloordiv(x, y)

Equivalent to `x // y` in Python.
"""
pyfloordiv(x, y) = pynew(errcheck(@autopy x y C.PyNumber_FloorDivide(getptr(x_), getptr(y_))))
"""
    pytruediv(x, y)

Equivalent to `x / y` in Python.
"""
pytruediv(x, y) = pynew(errcheck(@autopy x y C.PyNumber_TrueDivide(getptr(x_), getptr(y_))))
"""
    pymod(x, y)

Equivalent to `x % y` in Python.
"""
pymod(x, y) = pynew(errcheck(@autopy x y C.PyNumber_Remainder(getptr(x_), getptr(y_))))
"""
    pydivmod(x, y)

Equivalent to `divmod(x, y)` in Python.
"""
pydivmod(x, y) = pynew(errcheck(@autopy x y C.PyNumber_Divmod(getptr(x_), getptr(y_))))
"""
    pylshift(x, y)

Equivalent to `x << y` in Python.
"""
pylshift(x, y) = pynew(errcheck(@autopy x y C.PyNumber_Lshift(getptr(x_), getptr(y_))))
"""
    pyrshift(x, y)

Equivalent to `x >> y` in Python.
"""
pyrshift(x, y) = pynew(errcheck(@autopy x y C.PyNumber_Rshift(getptr(x_), getptr(y_))))
"""
    pyand(x, y)

Equivalent to `x & y` in Python.
"""
pyand(x, y) = pynew(errcheck(@autopy x y C.PyNumber_And(getptr(x_), getptr(y_))))
"""
    pyxor(x, y)

Equivalent to `x ^ y` in Python.
"""
pyxor(x, y) = pynew(errcheck(@autopy x y C.PyNumber_Xor(getptr(x_), getptr(y_))))
"""
    pyor(x, y)

Equivalent to `x | y` in Python.
"""
pyor(x, y) = pynew(errcheck(@autopy x y C.PyNumber_Or(getptr(x_), getptr(y_))))
export pyadd, pysub, pymul, pymatmul, pyfloordiv, pytruediv, pymod, pydivmod, pylshift, pyrshift, pyand, pyxor, pyor

# binary in-place
"""
    pyiadd(x, y)

In-place add. `x = pyiadd(x, y)` is equivalent to `x += y` in Python.
"""
pyiadd(x, y) = pynew(errcheck(@autopy x y C.PyNumber_InPlaceAdd(getptr(x_), getptr(y_))))
"""
    pyisub(x, y)

In-place subtract. `x = pyisub(x, y)` is equivalent to `x -= y` in Python.
"""
pyisub(x, y) = pynew(errcheck(@autopy x y C.PyNumber_InPlaceSubtract(getptr(x_), getptr(y_))))
"""
    pyimul(x, y)

In-place multiply. `x = pyimul(x, y)` is equivalent to `x *= y` in Python.
"""
pyimul(x, y) = pynew(errcheck(@autopy x y C.PyNumber_InPlaceMultiply(getptr(x_), getptr(y_))))
"""
    pyimatmul(x, y)

In-place matrix multiply. `x = pyimatmul(x, y)` is equivalent to `x @= y` in Python.
"""
pyimatmul(x, y) = pynew(errcheck(@autopy x y C.PyNumber_InPlaceMatrixMultiply(getptr(x_), getptr(y_))))
"""
    pyifloordiv(x, y)

In-place floor divide. `x = pyifloordiv(x, y)` is equivalent to `x //= y` in Python.
"""
pyifloordiv(x, y) = pynew(errcheck(@autopy x y C.PyNumber_InPlaceFloorDivide(getptr(x_), getptr(y_))))
"""
    pyitruediv(x, y)

In-place true division. `x = pyitruediv(x, y)` is equivalent to `x /= y` in Python.
"""
pyitruediv(x, y) = pynew(errcheck(@autopy x y C.PyNumber_InPlaceTrueDivide(getptr(x_), getptr(y_))))
"""
    pyimod(x, y)

In-place subtraction. `x = pyimod(x, y)` is equivalent to `x %= y` in Python.
"""
pyimod(x, y) = pynew(errcheck(@autopy x y C.PyNumber_InPlaceRemainder(getptr(x_), getptr(y_))))
"""
    pyilshift(x, y)

In-place left shift. `x = pyilshift(x, y)` is equivalent to `x <<= y` in Python.
"""
pyilshift(x, y) = pynew(errcheck(@autopy x y C.PyNumber_InPlaceLshift(getptr(x_), getptr(y_))))
"""
    pyirshift(x, y)

In-place right shift. `x = pyirshift(x, y)` is equivalent to `x >>= y` in Python.
"""
pyirshift(x, y) = pynew(errcheck(@autopy x y C.PyNumber_InPlaceRshift(getptr(x_), getptr(y_))))
"""
    pyiand(x, y)

In-place and. `x = pyiand(x, y)` is equivalent to `x &= y` in Python.
"""
pyiand(x, y) = pynew(errcheck(@autopy x y C.PyNumber_InPlaceAnd(getptr(x_), getptr(y_))))
"""
    pyixor(x, y)

In-place xor. `x = pyixor(x, y)` is equivalent to `x ^= y` in Python.
"""
pyixor(x, y) = pynew(errcheck(@autopy x y C.PyNumber_InPlaceXor(getptr(x_), getptr(y_))))
"""
    pyior(x, y)

In-place or. `x = pyior(x, y)` is equivalent to `x |= y` in Python.
"""
pyior(x, y) = pynew(errcheck(@autopy x y C.PyNumber_InPlaceOr(getptr(x_), getptr(y_))))
export pyiadd, pyisub, pyimul, pyimatmul, pyifloordiv, pyitruediv, pyimod, pyidivmod, pyilshift, pyirshift, pyiand, pyixor, pyior

# power
"""
    pypow(x, y, z=None)

Equivalent to `x ** y` or `pow(x, y, z)` in Python.
"""
pypow(x, y, z=pybuiltins.None) = pynew(errcheck(@autopy x y z C.PyNumber_Power(getptr(x_), getptr(y_), getptr(z_))))
"""
    pyipow(x, y, z=None)

In-place power. `x = pyipow(x, y)` is equivalent to `x **= y` in Python.
"""
pyipow(x, y, z=pybuiltins.None) = pynew(errcheck(@autopy x y z C.PyNumber_InPlacePower(getptr(x_), getptr(y_), getptr(z_))))
export pypow, pyipow


text:
#= = # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   fluxtable.txt
#
#   This file stores the F10.7 in different formats.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # = =#

struct _fluxtable_Structure
    F10obs::_space_indices_itp_constant{Float64}
    F10adj::_space_indices_itp_constant{Float64}
end

# Remote file: fluxtable.txt
_fluxtable = @RemoteFile(
    "ftp://ftp.seismo.nrcan.gc.ca/spaceweather/solar_flux/daily_flux_values/fluxtable.txt",
    file = "fluxtable.txt",
    dir = "$path_si_root",
    updates = :daily
)

# Optional variable that will store the `fluxtable.txt` data.
@OptionalData(_fluxtable_data, _fluxtable_Structure,
              "Run `init_space_indices()` with `:fluxtable` in `enabled_files` array to initialize required data.")

################################################################################
#                              Private Functions
################################################################################

"""
    _init_fluxtable(;force_download = false, local_path = nothing)

Initialize the data in the file `fluxtable.txt` by creating `_fluxtable_data`.
The initialization process is composed of:

1. Download the file, if it is necessary;
2. Parse the file;
3. Create the interpolations and the structures.

If the keyword `force_download` is `true`, then the file will always be
downloaded.

The user can also specify a location for the file using the keyword
`local_path`. If it is `nothing`, which is the default, then the file will be
downloaded.

"""
function _init_fluxtable(;force_download=false, local_path=nothing)
    # Update the remote files if no path is given.
    if local_path == nothing
        download(_fluxtable; force=force_download)
        local_path = path(_fluxtable)
    end

    push!(_fluxtable_data,   _parse_fluxtable(local_path))

    nothing
end

"""
    _parse_fluxtable(path::AbstractString)

Parse the `fluxtable.txt` file in `path` and return an instance of the structure
`_fluxtable_Structure` with the initialize interpolations.

"""
function _parse_fluxtable(path::AbstractString)
    # Allocate raw data.
    DT      = DateTime[]
    F10obs  = Float64[]
    F10adj  = Float64[]

    open(path) do file
        line_num = 0

        for ln in eachline(file)
            line_num += 1

            # Jump comments at the beginning of the file.
            (line_num in [1,2]) && continue

            # Get the tokens.
            tokens   = split(ln)
            num_toks = length(tokens)

            # Check if the line is valid.
            (num_toks != 7) && @error "Error parsing the line $line_num of fluxtable.txt"

            # We will save only the flux obtained at 20h. This is what
            # NRLMSISE-00 and JB2008 uses.
            fluxtime = parse(Int, tokens[2])
            (fluxtime != 200000) && continue

            # The DT of the data will be computed at noon. Hence, we will be
            # able to use the nearest-neighbor algorithm in the interpolations.
            year  = parse(Int, tokens[1][1:4])
            month = parse(Int, tokens[1][5:6])
            day   = parse(Int, tokens[1][7:8])
            DT_k  = DateTime(year, month, day, 12, 0, 0)

            # Get the raw data.

            if isempty(DT) || DT_k != DT[end]
                push!(DT,      DT_k)
                push!(F10obs,  parse(Float64, tokens[5]))
                push!(F10adj,  parse(Float64, tokens[6]))
            else

                DT[end] =  DT_k
                F10obs[end] = parse(Float64, tokens[5])
                F10adj[end] = parse(Float64, tokens[6])
            end
        end
    end

    # Create the interpolations for each parameter.
    knots = (datetime2julian.(DT),)
    # @show knots
    # @show F10obs
    # throw(error("test"))
    itp_F10obs  = interpolate(knots, F10obs,  Gridded(Constant()))
    itp_F10adj  = interpolate(knots, F10adj,  Gridded(Constant()))

    _fluxtable_Structure(itp_F10obs, itp_F10adj)
end


text:
module SNOW

using ForwardDiff
using ReverseDiff
using Zygote
using DiffResults
using FiniteDiff
using SparseArrays
using SparseDiffTools
using Requires
using Ipopt

include("derivatives.jl")

export ForwardAD, ReverseAD, RevZyg, ForwardFD, CentralFD, ComplexStep, UserDeriv
export DensePattern, SparsePattern

include("interface.jl")

export minimize
export Options

include("ipopt.jl")

export IPOPT

# conditionally load Snopt
function __init__()
    @require Snopt="0e9dc826-d618-11e8-1f57-c34e87fde2c0" include("snopt.jl")
end

end


text:
# Test indexes
@testset "Indexes" begin
    # Abstract types
    @test AbstractInfOptIndex isa DataType
    @test ObjectIndex isa DataType
    @test ObjectIndex <: AbstractInfOptIndex
    # IndependentParameterIndex
    @test IndependentParameterIndex <: ObjectIndex
    @test IndependentParameterIndex(1).value == 1
    # DependentParametersIndex
    @test DependentParametersIndex <: ObjectIndex
    @test DependentParametersIndex(1).value == 1
    # DependentParameterIndex
    @test DependentParameterIndex <: AbstractInfOptIndex
    idx = DependentParametersIndex(1)
    @test DependentParameterIndex(idx, 1).object_index == idx
    @test DependentParameterIndex(idx, 1).param_index == 1
    # InfiniteParameterIndex
    @test InfiniteParameterIndex <: AbstractInfOptIndex
    # FiniteParameterIndex
    @test FiniteParameterIndex <: ObjectIndex
    @test FiniteParameterIndex(1).value == 1
    # ParameterFunctionIndex
    @test ParameterFunctionIndex <: ObjectIndex
    @test ParameterFunctionIndex(1).value == 1
    # InfiniteVariableIndex
    @test InfiniteVariableIndex <: ObjectIndex
    @test InfiniteVariableIndex(1).value == 1
    # ReducedVariableIndex
    @test ReducedVariableIndex <: ObjectIndex
    @test ReducedVariableIndex(1).value == 1
    # PointVariableIndex
    @test PointVariableIndex <: ObjectIndex
    @test PointVariableIndex(1).value == 1
    # HoldVariableIndex
    @test HoldVariableIndex <: ObjectIndex
    @test HoldVariableIndex(1).value == 1
    # DerivativeIndex
    @test DerivativeIndex <: ObjectIndex
    @test DerivativeIndex(1).value == 1
    # MeasureIndex
    @test MeasureIndex <: ObjectIndex
    @test MeasureIndex(1).value == 1
    # ConstraintIndex
    @test ConstraintIndex <: ObjectIndex
    @test ConstraintIndex(1).value == 1
    # test CleverDict extensions
    @testset "CleverDict Extensions" begin
        # index_to_key
        @test MOIUC.index_to_key(IndependentParameterIndex, Int64(42)) == IndependentParameterIndex(42)
        @test MOIUC.index_to_key(DependentParametersIndex, Int64(42)) == DependentParametersIndex(42)
        @test MOIUC.index_to_key(FiniteParameterIndex, Int64(42)) == FiniteParameterIndex(42)
        @test MOIUC.index_to_key(ParameterFunctionIndex, Int64(42)) == ParameterFunctionIndex(42)
        @test MOIUC.index_to_key(InfiniteVariableIndex, Int64(42)) == InfiniteVariableIndex(42)
        @test MOIUC.index_to_key(ReducedVariableIndex, Int64(42)) == ReducedVariableIndex(42)
        @test MOIUC.index_to_key(PointVariableIndex, Int64(42)) == PointVariableIndex(42)
        @test MOIUC.index_to_key(HoldVariableIndex, Int64(42)) == HoldVariableIndex(42)
        @test MOIUC.index_to_key(DerivativeIndex, Int64(42)) == DerivativeIndex(42)
        @test MOIUC.index_to_key(MeasureIndex, Int64(42)) == MeasureIndex(42)
        @test MOIUC.index_to_key(ConstraintIndex, Int64(42)) == ConstraintIndex(42)
        # index_to_key
        @test MOIUC.key_to_index(MeasureIndex(42)) == 42
    end
    # test Base.length
    @testset "Base.length" begin
        @test length(MeasureIndex(1)) == 1
    end
    # test Base.broadcastable
    @testset "Base.broadcastable" begin
        @test Base.broadcastable(MeasureIndex(1)) isa Ref
    end
    # test Base.iterate
    @testset "Base.iterate" begin
        @test iterate(MeasureIndex(1)) == (MeasureIndex(1), true)
        @test iterate(MeasureIndex(1), true) isa Nothing
    end
    # test Base.hash
    @testset "Base.hash" begin
        @test Base.hash(MeasureIndex(1), UInt(1)) isa UInt
    end
end

# Test Infinite Sets
@testset "Infinite Sets" begin
    # Abstract types
    @test AbstractInfiniteSet isa DataType
    @test InfiniteScalarSet <: AbstractInfiniteSet
    @test InfiniteArraySet <: AbstractInfiniteSet
    # IntervalSet
    @test IntervalSet <: InfiniteScalarSet
    @test IntervalSet(0, 1) == IntervalSet(0.0, 1.0)
    @test_throws ErrorException IntervalSet(1, 0)
    # UniDistributionSet
    @test UniDistributionSet <: InfiniteScalarSet
    @test UniDistributionSet(Normal()) isa UniDistributionSet{<:Normal}
    # NonUnivariateDistribution
    @test InfiniteOpt.NonUnivariateDistribution <: Distribution
    # MultiDistributionSet
    @test MultiDistributionSet <: InfiniteArraySet
    @test MultiDistributionSet(Dirichlet([0.5, 0.5])) isa MultiDistributionSet{<:Dirichlet}
    # CollectionSet
    @test CollectionSet <: InfiniteArraySet
    @test CollectionSet([IntervalSet(0, 1), IntervalSet(2, 3)]) isa CollectionSet{IntervalSet}
end

# Test Generative Support Info
@testset "Generative Support Info" begin
    # Abstract types
    @test AbstractGenerativeInfo isa DataType
    # NoGenerativeSupports
    @test NoGenerativeSupports <: AbstractGenerativeInfo
    @test NoGenerativeSupports() isa NoGenerativeSupports
    @test NoGenerativeSupports() == NoGenerativeSupports()
    # UniformGenerativeInfo
    @test UniformGenerativeInfo <: AbstractGenerativeInfo
    @test UniformGenerativeInfo([0.0, 1.0], InternalLabel) isa UniformGenerativeInfo
    @test UniformGenerativeInfo([-2, 0, 1], InternalLabel, -2, 1).support_basis == [0, 2.0/3.0, 1]
    @test UniformGenerativeInfo([0.0, 1.0], InternalLabel) == UniformGenerativeInfo([2, 6], InternalLabel, 2, 6)
    @test UniformGenerativeInfo([0.0, 1.0], InternalLabel) != UniformGenerativeInfo([2, 5], InternalLabel, 2, 6)
    @test_throws ErrorException UniformGenerativeInfo([2, 6], InternalLabel)
end

# Test Derivative Evaluation Methods
@testset "Derivative Methods" begin
    # abstract types
    @test AbstractDerivativeMethod isa DataType
    @test GenerativeDerivativeMethod <: AbstractDerivativeMethod
    @test NonGenerativeDerivativeMethod <: AbstractDerivativeMethod
    # test orthogonal collocation (defined in derivative_evaluations.jl)
    @test OrthogonalCollocation <: GenerativeDerivativeMethod
    @test OrthogonalCollocation{GaussLobatto}(3, GaussLobatto()) isa OrthogonalCollocation
    @test OrthogonalCollocation(3, GaussLobatto()) == OrthogonalCollocation{GaussLobatto}(3, GaussLobatto())
    @test OrthogonalCollocation(3) == OrthogonalCollocation(3, GaussLobatto())
    @test_throws ErrorException OrthogonalCollocation(1)
    # test finite difference
    @test FDTechnique isa DataType
    @test Forward <: FDTechnique
    @test Central <: FDTechnique
    @test Backward <: FDTechnique
    @test FiniteDifference <: NonGenerativeDerivativeMethod
    @test FiniteDifference(Forward()) isa FiniteDifference
    @test FiniteDifference() == FiniteDifference(Backward())
    @test FiniteDifference(Forward(), false) isa FiniteDifference
end

# Test parameter datatypes
@testset "Parameters" begin
    # abstract types
    @test InfOptParameter <: AbstractVariable
    @test ScalarParameter <: InfOptParameter
    @test AbstractDataObject isa DataType
    # test IndependentParameter
    @test IndependentParameter <: ScalarParameter
    dict = SortedDict{Float64, Set{DataType}}(2 => Set([All]))
    method = FiniteDifference()
    info = NoGenerativeSupports()
    @test IndependentParameter(IntervalSet(0, 1), dict, 6, method, info).set isa IntervalSet
    # test FiniteParameter
    @test FiniteParameter <: ScalarParameter
    @test FiniteParameter(1) == FiniteParameter(1)
    # test DependentParameters
    @test DependentParameters <: InfOptParameter
    @test DependentParameters(CollectionSet([IntervalSet(0, 1)]),
                              Dict(zeros(1) => Set([All])), 6, [method]).set isa CollectionSet
    # test ScalarParameterData
    @test ScalarParameterData <: AbstractDataObject
    @test ScalarParameterData(FiniteParameter(42), 1, 1, "bob").name == "bob"
    # test MultiParameterData
    @test MultiParameterData <: AbstractDataObject
    params = DependentParameters(CollectionSet([IntervalSet(0, 1)]),
                                 Dict(zeros(1) => Set([All])), 6, [method])
    @test MultiParameterData(params, 1, 1:1, ["par[1]"]) isa MultiParameterData
end

# Test the InfiniteModel datatype
@testset "InfiniteModel" begin
    # test basic
    @test InfiniteModel <: JuMP.AbstractModel
    @test InfiniteModel().ext isa Dict
    # prepare optimizer constructor
    mockoptimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                             eval_objective_value=false)
    mockattributes = MOI.OptimizerWithAttributes(mockoptimizer, MOI.Silent() => true)
    # test optimizer constructors
    @test InfiniteModel(mockoptimizer).optimizer_constructor == mockoptimizer
    @test InfiniteModel(mockattributes).optimizer_constructor == mockoptimizer
    m = InfiniteModel();
    @test isa(Base.broadcastable(m), Base.RefValue{InfiniteModel})
    @test length(JuMP.object_dictionary(m)) == 0
    @test InfiniteModel(caching_mode = MOIU.MANUAL) isa JuMP.AbstractModel
    @test InfiniteModel(mockoptimizer,
                        caching_mode = MOIU.MANUAL) isa JuMP.AbstractModel
    # test accessors
    @test InfiniteOpt._last_param_num(m) == 0
    @test InfiniteOpt._param_object_indices(m) isa Vector{Union{IndependentParameterIndex, DependentParametersIndex}}
end

# Test reference variable datatypes
@testset "References" begin
    m = InfiniteModel();
    # Abstract types
    @test DispatchVariableRef <: AbstractVariableRef
    @test MeasureFiniteVariableRef <: DispatchVariableRef
    @test FiniteVariableRef <: MeasureFiniteVariableRef
    # test GeneralVariableRef
    @test GeneralVariableRef <: AbstractVariableRef
    @test GeneralVariableRef(m, 1, MeasureIndex) isa GeneralVariableRef
    # test IndependentParameterRef
    @test IndependentParameterRef <: DispatchVariableRef
    idx = IndependentParameterIndex(1)
    @test IndependentParameterRef(m, idx) isa IndependentParameterRef
    # test DependentParameterRef
    @test DependentParameterRef <: DispatchVariableRef
    idx = DependentParameterIndex(DependentParametersIndex(1), 1)
    @test DependentParameterRef(m, idx) isa DependentParameterRef
    # test FiniteParameterRef
    @test FiniteParameterRef <: FiniteVariableRef
    idx = FiniteParameterIndex(1)
    @test FiniteParameterRef(m, idx) isa FiniteParameterRef
    # InfiniteVariableRef
    @test InfiniteVariableRef <: DispatchVariableRef
    idx = InfiniteVariableIndex(1)
    @test InfiniteVariableRef(m, idx) isa InfiniteVariableRef
    # ReducedVariableRef
    @test ReducedVariableRef <: DispatchVariableRef
    idx = ReducedVariableIndex(1)
    @test ReducedVariableRef(m, idx) isa ReducedVariableRef
    # ParameterFunctionRef
    @test ParameterFunctionRef <: DispatchVariableRef
    idx = ParameterFunctionIndex(1)
    @test ParameterFunctionRef(m, idx) isa ParameterFunctionRef
    # PointVariableRef
    @test PointVariableRef <: FiniteVariableRef
    idx = PointVariableIndex(1)
    @test PointVariableRef(m, idx) isa PointVariableRef
    # HoldVariableRef
    @test HoldVariableRef <: FiniteVariableRef
    idx = HoldVariableIndex(1)
    @test HoldVariableRef(m, idx) isa HoldVariableRef
    # DerivativeRef
    @test DerivativeRef <: DispatchVariableRef
    idx = DerivativeIndex(1)
    @test DerivativeRef(m, idx) isa DerivativeRef
    # MeasureRef
    @test MeasureRef <: MeasureFiniteVariableRef
    idx = MeasureIndex(1)
    @test MeasureRef(m, idx) isa MeasureRef
    # InfOptConstraintRef
    @test InfOptConstraintRef isa DataType
    idx = ConstraintIndex(1)
    @test InfOptConstraintRef(m, idx) isa InfOptConstraintRef
end

# Test ParameterBounds
@testset "Parameter Bounds" begin
    # prepare the data
    m = InfiniteModel()
    idx = DependentParametersIndex(1)
    idx1 = DependentParameterIndex(idx, 1)
    idx2 = DependentParameterIndex(idx, 2)
    idx3 = IndependentParameterIndex(1)
    par1 = GeneralVariableRef(m, 1, DependentParameterIndex, 1)
    par2 = GeneralVariableRef(m, 1, DependentParameterIndex, 2)
    par3 = GeneralVariableRef(m, 1, IndependentParameterIndex, -1)
    pars = [par1, par2]
    # test datatype
    @testset "DataType" begin
        @test ParameterBounds isa UnionAll
        d = Dict(par3 => IntervalSet(0, 1))
        @test ParameterBounds(d) isa ParameterBounds{GeneralVariableRef}
        @test ParameterBounds() isa ParameterBounds{GeneralVariableRef}
    end
    # test _expand_parameter_tuple
    @testset "_expand_parameter_tuple" begin
        d = (par3 => IntervalSet(0, 1),)
        @test InfiniteOpt._expand_parameter_tuple(d) == Dict(d...)
    end
    # test _expand_parameter_dict(Dict{ParameterRef,IntervalSet}))
    @testset "_expand_parameter_dict (acceptable Form)" begin
        d = (par3 => IntervalSet(0, 1),)
        @test InfiniteOpt._expand_parameter_dict(d) == Dict(d...)
    end
    # test _expand_parameter_dict(Dict{Any,IntervalSet}))
    @testset "_expand_parameter_dict (Array Form)" begin
        d = (pars => IntervalSet(0, 1), par3 => IntervalSet(0, 1))
        @test isa(InfiniteOpt._expand_parameter_dict(d),
                  Dict{GeneralVariableRef, IntervalSet})
    end
    # test _expand_parameter_dict(Dict))
    @testset "_expand_parameter_dict (Fallback)" begin
        d = (pars => 1, par3 => 2)
        @test_throws ErrorException InfiniteOpt._expand_parameter_dict(d)
    end
    # test expansion definition
    @testset "ParameterBounds Expansion" begin
        d = (par3 => IntervalSet(0, 1),)
        @test ParameterBounds(d).intervals == Dict(d...)
        d = (pars => IntervalSet(0, 1), par3 => IntervalSet(0, 1))
        @test ParameterBounds(d).intervals isa Dict{GeneralVariableRef, IntervalSet}
    end
    pb = ParameterBounds((par3 => IntervalSet(0, 1),))
    # test intervals
    @testset "intervals" begin
        @test intervals(pb) == pb.intervals
    end
    # test Base.length
    @testset "Base.length" begin
        @test length(pb) == 1
    end
    # test Base.isempty
    @testset "Base.isempty" begin
        @test !isempty(pb)
        @test isempty(ParameterBounds())
    end
    # test Base.:(==)
    @testset "Base.:(==)" begin
        @test pb == ParameterBounds((par3 => IntervalSet(0, 1),))
        @test pb != ParameterBounds()
    end
    # test Base.copy
    @testset "Base.copy" begin
        @test copy(pb) == pb
    end
    # test Base.getindex
    @testset "Base.getindex" begin
        @test pb[par3] == IntervalSet(0, 1)
    end
    # test Base.setindex!
    @testset "Base.setindex!" begin
        @test (pb[par3] = IntervalSet(0, 2)) == IntervalSet(0, 2)
    end
    # test Base.haskey
    @testset "Base.haskey" begin
        @test haskey(pb, par3)
        @test !haskey(pb, par2)
    end
    # test Base.keys
    @testset "Base.keys" begin
        @test keys(pb) == keys(intervals(pb))
    end
    # test Base.iterate
    @testset "Base.iterate" begin
        @test [p for p in pb] == [par3 => IntervalSet(0, 2)]
    end
    # test Base.merge
    @testset "Base.merge" begin
        new_pb = ParameterBounds()
        @test merge(new_pb, pb) == pb
    end
    # test Base.merge!
    @testset "Base.merge!" begin
        new_pb = ParameterBounds()
        @test merge!(new_pb, pb) isa ParameterBounds
        @test new_pb == pb
    end
    # test Base.filter
    @testset "Base.filter" begin
        @test filter(e -> e[1] != par3, pb) == ParameterBounds()
    end
    # test Base.delete!
    @testset "Base.delete!" begin
        @test delete!(pb, par3) isa ParameterBounds
        @test length(pb) == 0
    end
end

# Test Infinite parameter functions
@testset "Infinite Parameter Functions" begin
    # useful data
    m = InfiniteModel()
    pref = GeneralVariableRef(m, 1, IndependentParameterIndex, -1)
    vt = IC.VectorTuple(pref)
    # test InfiniteParameterFunction
    @test InfiniteParameterFunction(sin, vt, [1], [1]) isa InfiniteParameterFunction
    # test ParameterFunctionData
    @test ParameterFunctionData(InfiniteParameterFunction(sin, vt, [1], [1])) isa ParameterFunctionData
end

# Test variable datatypes
@testset "Variables" begin
    # useful data
    m = InfiniteModel()
    num = Float64(0)
    sample_info = VariableInfo(true, num, true, num, true, num, true, num, true, true)
    func = (x) -> NaN
    inf_info = VariableInfo{Float64, Float64, Float64, Function}(true, num, true,
                 num, true, num, false, func, true, true)
    pref = GeneralVariableRef(m, 1, IndependentParameterIndex, -1)
    vref = GeneralVariableRef(m, 1, InfiniteVariableIndex)
    # Abstract variables
    @test InfOptVariable <: AbstractVariable
    # Infinite variable
    @test InfiniteVariable <: InfOptVariable
    @test InfiniteVariable(inf_info, IC.VectorTuple(pref), [1], [1], true) isa InfiniteVariable
    # Reduced variable
    @test ReducedVariable <: InfOptVariable
    @test ReducedVariable(vref, Dict(1 => Float64(2)), [1], [1]) isa ReducedVariable
    # Point variable
    @test PointVariable <: InfOptVariable
    @test PointVariable(sample_info, vref, Float64[1]) isa PointVariable
    # Hold variable
    @test HoldVariable <: InfOptVariable
    @test HoldVariable(sample_info, ParameterBounds()) isa HoldVariable
    # VariableData
    @test VariableData <: AbstractDataObject
    @test VariableData(HoldVariable(sample_info, ParameterBounds())) isa VariableData
end

# Test derivative datatypes
@testset "Derivatives" begin
    # useful data
    m = InfiniteModel()
    num = Float64(0)
    sample_info = VariableInfo(true, num, true, num, true, num, true, num, true, true)
    func = (x) -> NaN
    inf_info = VariableInfo{Float64, Float64, Float64, Function}(true, num, true,
                 num, true, num, false, func, true, true)
    pref = GeneralVariableRef(m, 1, IndependentParameterIndex, -1)
    vref = GeneralVariableRef(m, 1, InfiniteVariableIndex)
    dref = GeneralVariableRef(m, 1, DerivativeIndex)
    # derivative
    @test Derivative <: InfOptVariable
    @test Derivative(inf_info, true, vref, pref) isa Derivative
    # Reduced derivative
    @test ReducedVariable(dref, Dict(1 => Float64(2)), [1], [1]) isa ReducedVariable
    # Point derivative
    @test PointVariable(sample_info, dref, Float64[1]) isa PointVariable
    # VariableData
    @test VariableData(Derivative(inf_info, true, vref, pref)) isa VariableData
end

# Test the measure datatypes
@testset "Measures" begin
    # Helper data
    m = InfiniteModel()
    pref = GeneralVariableRef(m, 2, IndependentParameterIndex)
    w(t) = 1
    info = NoGenerativeSupports()
    # abstract types
    @test AbstractMeasureData isa DataType
    # DiscreteMeasureData
    @test DiscreteMeasureData <: AbstractMeasureData
    @test DiscreteMeasureData(pref, ones(2), ones(2), All, w, NaN, NaN, false) isa DiscreteMeasureData
    @test DiscreteMeasureData([pref], ones(2), ones(1, 2), All, w, [NaN], [NaN], false) isa DiscreteMeasureData
    # FunctionalDistributionSet
    @test FunctionalDiscreteMeasureData <: AbstractMeasureData
    @test FunctionalDiscreteMeasureData(pref, w, 2, All, info, w, NaN, NaN, false) isa FunctionalDiscreteMeasureData
    @test FunctionalDiscreteMeasureData([pref], w, 2, All, w, [NaN], [NaN], false) isa FunctionalDiscreteMeasureData
    @test FunctionalDiscreteMeasureData([pref], w, 2, All, info, w, [NaN], [NaN], false) isa FunctionalDiscreteMeasureData
    # Measure
    @test Measure isa UnionAll
    @test Measure(zero(AffExpr),
                  DiscreteMeasureData(pref, ones(2), ones(2), All, w, NaN, NaN, false),
                  [1], [1], false) isa Measure
    # MeasureData
    @test MeasureData <: AbstractDataObject
    @test MeasureData(Measure(zero(AffExpr), DiscreteMeasureData(pref,
                      ones(2), ones(2), All, w, NaN, NaN, false), [1], [1], true)) isa MeasureData
end

# Test the constraint datatypes
@testset "Constraints" begin
    # Setup
    m = InfiniteModel()
    pref = GeneralVariableRef(m, 2, IndependentParameterIndex, -1)
    bounds = ParameterBounds(Dict(pref => IntervalSet(0, 1)))
    # Bounded constraints
    @test BoundedScalarConstraint <: JuMP.AbstractConstraint
    @test BoundedScalarConstraint(zero(AffExpr), MOI.Integer(),
                bounds, bounds).bounds.intervals[pref].lower_bound == 0.0
    # ConstraintData
    @test ConstraintData <: AbstractDataObject
    @test ConstraintData(ScalarConstraint(zero(AffExpr), MOI.Integer()), [1], "",
                         MeasureIndex[], false) isa ConstraintData
end


text:
function dealias!(out::T,u::T,d::Real) where {T<:AbstractField}
    isrealspace(u) && fourier!(u)
    isrealspace(out) && setfourier!(out)
    kx = u.kx
    ky = u.ky
    kz = u.kz
    kmax2 = (d*kx[end])^2

    @inbounds for k in axes(u,3)
        kz2 = kz[k]*kz[k]
        for j in axes(u,2)
            ky2 = ky[j]*ky[j]
            @simd for i in axes(u,1)
                kx2 = kx[i]*kx[i]
                out[i,j,k] = ((kx2<kmax2) | (ky2<kmax2) | (kz2<kmax2))*u[i,j,k]
            end
        end
    end

    return out
end