text:
@testset "deterministic_distribution_model" begin
    table = Array{
        Vector{NamedTuple{(:nextstate, :reward, :prob),Tuple{Int,Float64,Float64}}},
        2,
    }(
        undef,
        2,
        2,
    )

    table[1] = [
        (nextstate = 1, reward = 1.0, prob = 0.5),
        (nextstate = 2, reward = 0.0, prob = 0.5),
    ]
    table[2] = [
        (nextstate = 2, reward = -1.0, prob = 0.25),
        (nextstate = 1, reward = 0.0, prob = 0.75),
    ]
    table[3] = [
        (nextstate = 2, reward = -1.0, prob = 0.25),
        (nextstate = 1, reward = 0.5, prob = 0.75),
    ]
    table[4] = [
        (nextstate = 1, reward = 1.0, prob = 0.75),
        (nextstate = 2, reward = -0.5, prob = 0.25),
    ]

    m = DeterministicDistributionModel(table)

    @test length(observation_space(m)) == 2
    @test length(action_space(m)) == 2

    @test m(1, 1) == table[1, 1]
    @test m(1, 2) == table[1, 2]
    @test m(2, 1) == table[2, 1]
    @test m(2, 2) == table[2, 2]
end


text:
using RNGTest
using RNG

function bigcrush{T<:Number}(rng::RNG.AbstractRNG{T})
    p = RNGTest.wrap(r, T)
    RNGTest.bigcrushTestU01(p)
end

function speed_test{T<:Number}(rng::RNG.AbstractRNG{T}, n=100_000_000)
    A = Array{T}(n)
    rand!(rng, A)
    elapsed = @elapsed rand!(rng, A)
    elapsed * 1e9 / n * 8 / sizeof(T)
end

function bigcrush(rng::MersenneTwister)
    p = RNGTest.wrap(r, UInt64)
    RNGTest.bigcrushTestU01(p)
end

function speed_test(rng::MersenneTwister, n=100_000_000)
    T = UInt64
    A = Array{T}(n)
    rand!(rng, A)
    elapsed = @elapsed rand!(rng, A)
    elapsed * 1e9 / n * 8 / sizeof(T)
end

function test_all(rng::Base.AbstractRNG, n=100_000_000)
    fo = open("$TEST_NAME.log", "w")
    redirect_stdout(fo)
    println(TEST_NAME)
    speed = speed_test(rng, n)
    @printf "Speed Test: %.3f ns/64 bits\n" speed
    flush(fo)
    bigcrush(rng)
    close(fo)
end

ext: jl


text:
using Fun1d
using QuadGK
using Random
using Test

# Random rationals
function Base.rand(rng::AbstractRNG,
                   ::Random.SamplerType{Rational{T}}) where {T}
    return Rational{T}(T(rand(rng, -1000:1000)) // 1000)
end

const Rat128 = Rational{Int128}

@testset "Domain" begin
    S = Rat128
    domain = Domain{S}(0, 1)
    @test Fun1d.invariant(domain)
end

Random.seed!(0)
@testset "Grid" begin
    S = Rat128
    domain = Domain{S}(0, 1)
    grid = Grid(domain, 10)
    @test Fun1d.invariant(grid)

    # Test some specific basis functions
    dx = spacing(grid)
    @test dx > 0

    @test basis(grid, 1, 0 * dx) == 1
    @test basis(grid, 1, 1 // 2 * dx) == 1 // 2
    @test basis(grid, 1, dx) == 0

    @test basis(grid, 2, 0 * dx) == 0
    @test basis(grid, 2, 1 // 2 * dx) == 1 // 2
    @test basis(grid, 2, dx) == 1
    @test basis(grid, 2, 3 // 2 * dx) == 1 // 2
    @test basis(grid, 2, 2 * dx) == 0

    # Test that basis functions form a position basis: At grid point
    # `i`, the `i`-th basis function must be `1`, all other basis
    # function must be `0`.
    for i in 1:(grid.ncells + 1)
        x = location(grid, i)
        for j in 1:(grid.ncells + 1)
            y = basis(grid, j, x)
            @test y == (i == j)
        end
    end

    # Test that basis functions form a partition of unity: At each
    # point in the domain, the basis functions must sum up to `1`
    for i in 1:10
        x = lincom(1, grid.domain.xmin, 1000000, grid.domain.xmax,
                   rand(1:1000000))
        sumy = zero(S)
        for j in 1:(grid.ncells + 1)
            y = basis(grid, j, x)
            @test 0 ≤ y ≤ 1
            sumy += y
        end
        @test sumy == 1
    end

    # TODO: test compact support
end

Random.seed!(0)
@testset "GridFuns are vectors" begin
    S = Rat128
    T = Rat128
    domain = Domain{S}(0, 1)
    grid = Grid(domain, 10)

    for n in 1:10
        f = rand(GridFun{S,T}, grid)
        g = rand(GridFun{S,T}, grid)
        h = rand(GridFun{S,T}, grid)
        z = zero(GridFun{S,T}, grid)
        a = rand(T)
        b = rand(T)

        @test Fun1d.invariant(f)
        @test Fun1d.invariant(g)
        @test Fun1d.invariant(h)
        @test Fun1d.invariant(z)

        @test map(identity, f) == f
        @test map(x -> x + b, map(x -> a * x, f)) == map(x -> a * x + b, f)
        @test map(+, f, g) == f + g

        @test iszero(z)
        @test f == f
        @test +f == f
        @test f + z == f
        @test z + f == f
        @test -(-f) == f
        @test f + (-f) == z
        @test f - g == f + (-g)
        @test f + g == g + f
        @test (f + g) + h == f + (g + h)

        @test zero(T) * f == z
        @test one(T) * f == f
        @test (-one(T)) * f == -f
        @test a * f == f * a
        @test a * (f + g) == a * f + a * g
        @test (a + b) * f == a * f + b * f
        @test (a * b) * f == a * (b * f)
        if !iszero(a)
            @test f / a == f * inv(a)
            @test a \ f == inv(a) * f
        end
    end
end

Random.seed!(0)
@testset "GridFun sample/evaluate" begin
    S = Rat128
    T = Rat128
    domain = Domain{S}(0, 1)
    grid = Grid(domain, 10)

    m = rand(T)
    b = rand(T)
    fun(x) = m * x + b
    gf = sample(T, grid, fun)
    @test Fun1d.invariant(gf)

    for n in 1:10
        x = lincom(1, grid.domain.xmin, 1000000, grid.domain.xmax,
                   rand(1:1000000))
        @test evaluate(gf, x) == fun(x)
    end
end

Random.seed!(0)
@testset "GridFun project/evaluate" begin
    S = Float64
    T = Float64
    domain = Domain{S}(0, 1)
    grid = Grid(domain, 10)

    m = rand(T)
    b = rand(T)
    fun(x) = m * x + b
    gf = project(T, grid, fun)
    @test Fun1d.invariant(gf)

    atol = eps(T)^(T(3) / 4)
    for n in 1:10
        x = lincom(1, grid.domain.xmin, 1000000, grid.domain.xmax,
                   rand(1:1000000))
        @test isapprox(evaluate(gf, x), fun(x); atol=atol)
    end
end

Random.seed!(0)
@testset "GridFun convergence" begin
    S = Float64
    T = Float64
    domain = Domain{S}(0, 1)
    grid1 = Grid(domain, 20)
    grid2 = Grid(domain, 40)

    atol = eps(T)^(T(3) / 4)
    function int(f)
        I, E = quadgk(f, domain.xmin, domain.xmax; atol=atol)
        return I
    end
    norm2(f) = sqrt(int(x -> abs2(f(x))))

    fun(x) = sinpi(x)
    gf1 = sample(T, grid1, fun)
    gf2 = sample(T, grid2, fun)

    e1 = norm2(x -> evaluate(gf1, x) - fun(x))
    e2 = norm2(x -> evaluate(gf2, x) - fun(x))
    ratio = e1 / e2
    @test isapprox(ratio, 4; rtol=0.01)
end

@testset "GridFun derivative linearity" begin
    S = Rat128
    T = Rat128
    domain = Domain{S}(0, 1)
    grid = Grid(domain, 10)

    for n in 1:10
        f = rand(GridFun{S,T}, grid)
        g = rand(GridFun{S,T}, grid)
        z = zero(GridFun{S,T}, grid)
        a = rand(T)

        @test deriv(z) == z
        @test deriv(f + g) == deriv(f) + deriv(g)
        @test deriv(a * f) == a * deriv(f)

        @test deriv2(z) == z
        @test deriv2(f + g) == deriv2(f) + deriv2(g)
        @test deriv2(a * f) == a * deriv2(f)

        @test deriv2_sphere(z) == z
        @test deriv2_sphere(f + g) == deriv2_sphere(f) + deriv2_sphere(g)
        @test deriv2_sphere(a * f) == a * deriv2_sphere(f)
    end
end

@testset "GridFun integral accuracy" begin
    S = Float64
    T = Float64
    domain = Domain{S}(0, 1)
    grid1 = Grid(domain, 20)
    grid2 = Grid(domain, 40)

    atol = eps(T)^(T(3) / 4)

    fun(x) = sinpi(x)
    gf1 = sample(T, grid1, fun)
    gf2 = sample(T, grid2, fun)

    i0 = 2 / T(π)
    i1 = integrate(gf1)
    i2 = integrate(gf2)
    iratio = (i1 - i0) / (i2 - i0)
    @test isapprox(iratio, 4; atol=0.01)

    n0 = 1 / sqrt(T(2))
    n1 = norm2(gf1)
    n2 = norm2(gf2)
    # nratio = (n1 - n0) / (n2 - n0)
    # @test isapprox(nratio, 4; atol=0.01)
    # Strangely, this result is much more accurate than expected. I
    # won't complain.
    @test isapprox(n1, n0; atol=atol)
    @test isapprox(n2, n0; atol=atol)
end

@testset "GridFun derivative accuracy" begin
    S = Float64
    T = Float64
    domain = Domain{S}(0, 1)
    grid1 = Grid(domain, 20)
    grid2 = Grid(domain, 40)

    atol = eps(T)^(T(3) / 4)
    dx = spacing(grid1)
    function int(f)
        # Omit outer boundary to get good convergence
        I, E = quadgk(f, domain.xmin + dx, domain.xmax - dx; atol=atol)
        return I
    end
    norm2(f) = sqrt(int(x -> abs2(f(x))))

    fun(x) = sinpi(x)
    gf1 = sample(T, grid1, fun)
    gf2 = sample(T, grid2, fun)

    dfun(x) = π * cospi(x)
    dgf1 = deriv(gf1)
    dgf2 = deriv(gf2)

    d2fun(x) = -T(π)^2 * sinpi(x)
    d2gf1 = deriv2(gf1)
    d2gf2 = deriv2(gf2)

    # ef1 = norm2(x -> evaluate(gf1, x) - fun(x))
    # ef2 = norm2(x -> evaluate(gf2, x) - fun(x))
    edf1 = norm2(x -> evaluate(dgf1, x) - dfun(x))
    edf2 = norm2(x -> evaluate(dgf2, x) - dfun(x))
    ed2f1 = norm2(x -> evaluate(d2gf1, x) - d2fun(x))
    ed2f2 = norm2(x -> evaluate(d2gf2, x) - d2fun(x))

    ratio = edf1 / edf2
    @test isapprox(ratio, 4; rtol=0.01)
    ratio2 = ed2f1 / ed2f2
    @test isapprox(ratio2, 4; rtol=0.01)
end

@testset "GridFun spherical derivative accuracy" begin
    S = Float64
    T = Float64
    domain = Domain{S}(0, 1)
    grid1 = Grid(domain, 20)
    grid2 = Grid(domain, 40)

    atol = eps(T)^(T(3) / 4)
    dx = spacing(grid1)
    function int(f)
        # Omit outer boundary to get good convergence
        I, E = quadgk(f, domain.xmin, domain.xmax - dx; atol=atol)
        return I
    end
    # norm2(f) = sqrt(int(r -> 4 * T(π) * r^2 * abs2(f(r))))
    norm2(f) = sqrt(int(r -> abs2(f(r))))

    fun(r) = cospi(r)
    gf1 = sample(T, grid1, fun)
    gf2 = sample(T, grid2, fun)

    d2fun(r) = -T(π)^2 * (cospi(r) + 2 * sinc(r))
    d2gf1 = deriv2_sphere(gf1)
    d2gf2 = deriv2_sphere(gf2)

    ed2f1 = norm2(r -> evaluate(d2gf1, r) - d2fun(r))
    ed2f2 = norm2(r -> evaluate(d2gf2, r) - d2fun(r))

    ratio2 = ed2f1 / ed2f2
    @test isapprox(ratio2, 4; rtol=0.02)
end


text:
"""
    genPopulation!(genj::GeneticJuliaStructure)

Generates the initial population if it has not been previously generated. If it
has already been generated from a previous experiment, it does nothing. This is
to continue evolving a previously evolved population.

# Arguments
- `genj::GeneticJuliaStructure`: the main structure, which specified the generation method
    and its arguments, if any.

See also: [`setGenerator`](@ref)
"""
function genPopulation!(genj::GeneticJuliaStructure)
    if !isdefined(genj, :_population)
        genj._population = genPopulation_(genj._generator, genj._experimentInfo)
    end
end # function

"""
    evaluate!(genj::GeneticJuliaStructure, population::Array{Individual})

Evaluates the fitness of the individuals of a population and sets their values.

# Arguments
- `genj::GeneticJuliaStructure`: the main structure, which contains the Evaluator that
    specifies the fitness functions to evaluate for each individual.
- `population::Array{Individual}`: the population that is going to be evaluated.

See also: [`setEvaluator`](@ref)
"""
function evaluate!(genj::GeneticJuliaStructure, population::Array{Individual})
    evaluate!_(genj._evaluator, genj._experimentInfo, population, genj._stopCondition)
end # function

"""
    selectParents(genj::GeneticJuliaStructure)

Selects a set of individuals of the population.

# Arguments
- `genj::GeneticJuliaStructure`: the main structure, which specifies the selection method,
    its arguments, if any, and the number of individuals to select.

# Returns
The set of selected individuals.

See also: [`setSelector`](@ref)
"""
function selectParents(genj::GeneticJuliaStructure)
    selectParents_(genj._selector, genj._population, genj._evaluator, genj._experimentInfo._rng)
end # function

"""
    selectParents(genj::GeneticJuliaStructure, population::Array{Individual})

Selects a set of individuals of a population.

# Arguments
- `genj::GeneticJuliaStructure`: the main structure, which specifies the selection method,
    its arguments, if any, and the number of individuals to select.
- `population::Array{Individual}`: the population from which individuals are
    going to be selected.

# Returns
The set of selected individuals.

See also: [`setSelector`](@ref)
"""
function selectParents(genj::GeneticJuliaStructure, population::Array{Individual})
    selectParents_(genj._selector, population, genj._evaluator, genj._experimentInfo._rng)
end # function

"""
    cross(genj::GeneticJuliaStructure, selected::Array{Individual})

Performs a crossover method between a set of selected individuals. If crossover
probability is set to 0, it returns the selected individuals unchanged.

# Arguments
- `genj::GeneticJuliaStructure`: the main structure, which specifies the crossover method,
    its arguments, if any, the number of parents it takes and the number of children
    it returns.
- `selected::Array{Individual}`: the selected individuals that are going to
    be crossed one another.

# Returns
The offspring after crossing the selected individuals. The selected individuals
unchanged if crossover probability equals 0.

See also: [`setCrossoverOperator`](@ref)
"""
function cross(genj::GeneticJuliaStructure, selected::Array{Individual})
    if genj._crossoverOp._probability > 0.0
        return cross_(genj._crossoverOp, selected, genj._experimentInfo)
    else
        return selected
    end
end # function

"""
    cross(genj::GeneticJuliaStructure, selected::Array{Individual})

Performs a mutation method over a set of selected individuals. If mutation
probability is set to 0, it returns the selected individuals unchanged.

# Arguments
- `genj::GeneticJuliaStructure`: the main structure, which specifies the mutation method
    and its arguments, if any.
- `selected::Array{Individual}`: the selected individuals that are going to
    be mutated.

# Returns
The mutated offspring. The selected individuals unchanged if mutation
probability equals 0.

See also: [`setMutationOperator`](@ref)
"""
function mutate(genj::GeneticJuliaStructure, offspring::Array{Individual})
    if genj._mutationOp._probability > 0.0
        mutate_(genj._mutationOp, offspring, genj._experimentInfo)
    else
        return offspring
    end
end # function

"""
    replacePopulation!(genj::GeneticJuliaStructure, offspring::Array{Individual})

Replaces the population selecting individuals from the current population and its
offspring according to a criterium.

# Arguments
- `genj::GeneticJuliaStructure`: the main structure, which specifies the replacement method
    and its arguments, if any.
- `offspring::Array{Individual}`: the individuals that resulted from crossover and
    mutation operations.

See also: [`setReplacementOperator`](@ref)
"""
function replacePopulation!(genj::GeneticJuliaStructure, offspring::Array{Individual})
    genj._population = replacePopulation_(genj._replacementOp, genj._population,
        offspring, genj._evaluator, genj._experimentInfo._rng)
end # function

"""
    reached(genj::GeneticJuliaStructure)

Checks if any of the stop conditions has been fulfilled.

# Arguments
- `genj::GeneticJuliaStructure`: the main structure, which contains the stop conditions.

# Returns
`true` if any of the stop conditions has been fulfilled, `false` otherwise.

See also: [`setStopCondition`](@ref)
"""
function reached(genj::GeneticJuliaStructure)

    if genj._stopCondition._maxIterations != -1
        notifyIteration(genj._stopCondition)
    end

    if genj._stopCondition._maxIterNotImproving != -1
        currentBestFitness = genj._population[getBestIndividual(genj._population, genj._evaluator._compareFunction)]

        if !genj._evaluator._compareFunction(currentBestFitness, GenJ._stopCondition._bestFitness)
            notifyIterNotImproving(genj._stopCondition)
        else
            genj._stopCondition._bestFitness = currentBestFitness
            resetIterNotImproving(genj._stopCondition)
        end
    end

    reached_(genj._stopCondition)
end # function

"""
    initTime(genj::GeneticJuliaStructure)

Initiates the timer.

# Arguments
- `genj::GeneticJuliaStructure`: the main structure, which contains the time as a stop condition.

See also: [`setStopCondition`](@ref), [`reached`](@ref)
"""
function initTime(genj::GeneticJuliaStructure)
    if genj._stopCondition._maxTime != Inf
        setTime(genj._stopCondition, time())
    end
end # function

"""
    initBestFitness(genj::GeneticJuliaStructure)

Sets the best fitness of the current population.

# Arguments
- `genj::GeneticJuliaStructure`: the main structure.

See also: [`setStopCondition`](@ref), [`reached`](@ref)
"""
function initBestFitness(genj::GeneticJuliaStructure)
    if genj._stopCondition._maxIterNotImproving != -1
        fitnesses = getFitness(genj._population)
        genj._stopCondition._bestFitness = fitnesses[getBestIndividual(genj._population, genj._evaluator._compareFunction)]
    end
end # function

"""
    getBestIndividual(genj::GeneticJuliaStructure, n::Integer = 1)

Returns the indexes of the best individuals of the population.

# Arguments
- `genj::GeneticJuliaStructure`: the main structure.
- `n::Integer = 1`: number of best individuals to be returned.
"""
function getBestIndividual(genj::GeneticJuliaStructure, n::Integer = 1)
    return getBestIndividual(genj._population, genj._evaluator._compareFunction, n)
end # function

"""
    printInformation(genj::GeneticJuliaStructure; printFitness::Bool = true,
                     printBestFitness::Bool = true, printMeanFitness::Bool = true,
                     printSTDFitness::Bool = true)

Prints all the information of the experiment.
"""
function printInformation(genj::GeneticJuliaStructure; printFitness::Bool = true,
                          printBestFitness::Bool = true, printMeanFitness::Bool = true,
                          printSTDFitness::Bool = true, outputFile::String = "")
    if isdefined(genj._experimentInfo, :_experimentSummary)
        if hasGlobalFitnessFunction(genj._evaluator)
            printInformationWithGlobal_(genj._experimentInfo._experimentSummary, printFitness=printFitness,
                                        printBestFitness=printBestFitness, printMeanFitness=printMeanFitness,
                                        printSTDFitness=printSTDFitness, outputFile=outputFile)
        else
            printInformation_(genj._experimentInfo._experimentSummary, printFitness=printFitness,
                              printBestFitness=printBestFitness, printMeanFitness=printMeanFitness,
                              printSTDFitness=printSTDFitness, outputFile=outputFile)
        end
    else
        error("The experiment has not been set with an information summary (See \"setExperimentSummary\")")
    end

    return nothing
end # function

"""
    printLastInformation(genj::GeneticJuliaStructure)

Prints the last bit of information collected of the experiment.
"""
function printLastInformation(genj::GeneticJuliaStructure, currGen::Integer = 0)
    if currGen == 0
        currGen = getCurrentIteration(genj._stopCondition)
    end
    if hasGlobalFitnessFunction(genj._evaluator)
        printLastInformationWithGlobal_(genj._experimentInfo._experimentSummary, currGen)
    else
        printLastInformation_(genj._experimentInfo._experimentSummary, currGen)
    end

    return nothing
end # function

"""
    saveResults(genj::GeneticJuliaStructure)

documentation
"""
function saveResults(genj::GeneticJuliaStructure)
    if isdefined(genj._experimentInfo, :_experimentSummary)
        batchSize = getBatchSize(genj._experimentInfo._experimentSummary)
        currIter = getCurrentIteration(genj._stopCondition)

        if batchSize == 0 && reached_(genj._stopCondition) || currIter%batchSize == 0
            if displayFitness(genj._experimentInfo._experimentSummary)
                if hasGlobalFitnessFunction(genj._evaluator)
                    saveFitness(genj._experimentInfo._experimentSummary,
                                getFitness(genj._population),
                                currIter, getGlobalFitness(genj._population))
                else
                    saveFitness(genj._experimentInfo._experimentSummary,
                                getFitness(genj._population),
                                currIter)
                end
            end

            if displayBestFitness(genj._experimentInfo._experimentSummary)
                bestIndividual = genj._population[getBestIndividual(genj)[1]]
                bestFitness = getFitness(bestIndividual, all=true)
                bestGenotype = getGenotype(bestIndividual)

                if typeof(bestGenotype) <: GAGenotype
                    bestRep = getRepresentation(bestGenotype)
                elseif typeof(bestGenotype) <: GEGenotype
                    bestRep = getPhenotype(genj._experimentInfo._gpExperimentInfo._grammar, bestGenotype)
                else
                    bestRep = getPhenotype(bestGenotype)
                end

                if hasGlobalFitnessFunction(genj._evaluator)
                    saveBestFitness(genj._experimentInfo._experimentSummary,
                                    bestFitness, bestRep,
                                    currIter, getGlobalFitness(bestIndividual))
                else
                    saveBestFitness(genj._experimentInfo._experimentSummary,
                                    bestFitness, bestRep,
                                    currIter)
                end
            end

            if displayMeanFitness(genj._experimentInfo._experimentSummary)
                if hasGlobalFitnessFunction(genj._evaluator)
                    saveMeanFitness(genj._experimentInfo._experimentSummary,
                                    getFitness(genj._population),
                                    currIter, getGlobalFitness(genj._population))
                else
                    saveMeanFitness(genj._experimentInfo._experimentSummary,
                                    getFitness(genj._population),
                                    currIter)
                end
            end

            if displaySTDFitness(genj._experimentInfo._experimentSummary)
                if hasGlobalFitnessFunction(genj._evaluator)
                    saveSTDFitness(genj._experimentInfo._experimentSummary,
                                    getFitness(genj._population),
                                    currIter, getGlobalFitness(genj._population))
                else
                    saveSTDFitness(genj._experimentInfo._experimentSummary,
                                    getFitness(genj._population),
                                    currIter)
                end
            end

            if printDuringExperiment(genj._experimentInfo._experimentSummary)
                outputFile = getOutputFile(genj._experimentInfo._experimentSummary)
                io = outputFile != "" ? open(outputFile, "a") : Base.stdout
                flush(io)

                println(io, "GENERATION ", currIter, ":")
                println(io)
                io == Base.stdout || close(io)

                printLastInformation(genj, currIter)

                println(io, "============================")
                println(io)
                println(io)

                io == Base.stdout || close(io)
            end
        end
    end

    return nothing
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text: largez02() = largez(2)


text:
include("utils.jl")

@testset "Euclidean" begin
    E = Euclidean(3)
    Ec = Euclidean(3; field=ℂ)
    EM = Manifolds.MetricManifold(E, Manifolds.EuclideanMetric())
    @test repr(E) == "Euclidean(3; field = ℝ)"
    @test repr(Ec) == "Euclidean(3; field = ℂ)"
    @test repr(Euclidean(2, 3; field = ℍ)) == "Euclidean(2, 3; field = ℍ)"
    @test is_default_metric(EM)
    @test is_default_metric(E, Manifolds.EuclideanMetric())
    @test Manifolds.default_metric_dispatch(E,Manifolds.EuclideanMetric()) === Val{true}()
    p = zeros(3)
    @test det_local_metric(EM, p) == one(eltype(p))
    @test log_local_metric_density(EM, p) == zero(eltype(p))
    @test project!(E, p, p) == p
    @test manifold_dimension(Ec) == 2*manifold_dimension(E)
    X = zeros(3)
    X[1] = 1.0
    Y = similar(X)
    project!(E, Y, p, X)
    @test Y==X

    # real manifold does not allow complex values
    @test_throws DomainError is_manifold_point(Ec,[:a,:b,:b], true)
    @test_throws DomainError is_manifold_point(E,[1.0, 1.0im, 0.0], true)
    @test_throws DomainError is_manifold_point(E,[1], true)
    @test_throws DomainError is_tangent_vector(Ec,[:a,:b,:b], [1.0, 1.0, 0.0], true)
    @test_throws DomainError is_tangent_vector(E,[1.0, 1.0im, 0.0], [1.0, 1.0, 0.0], true) # real manifold does not allow complex values
    @test_throws DomainError is_tangent_vector(E,[1], [1.0, 1.0, 0.0], true)
    @test_throws DomainError is_tangent_vector(E, [0.0, 0.0, 0.0], [1.0], true)
    @test_throws DomainError is_tangent_vector(E, [0.0, 0.0, 0.0], [1.0, 0.0, 1.0im], true)
    @test_throws DomainError is_tangent_vector(Ec, [0.0, 0.0, 0.0], [:a, :b, :c], true)

    @test E^2 === Euclidean(3, 2)
    @test ^(E,2) === Euclidean(3, 2)
    @test E^(2,) === Euclidean(3, 2)
    @test Ec^(4,5) === Euclidean(3, 4, 5; field = ℂ)

    manifolds = [ E, EM, Ec ]
    types = [ Vector{Float64}, ]
    TEST_FLOAT32 && push!(types, Vector{Float32})
    TEST_DOUBLE64 && push!(types, Vector{Double64})
    TEST_STATIC_SIZED && push!(types, MVector{3, Float64})

    types_complex = [ Vector{ComplexF64}, ]
    TEST_FLOAT32 && push!(types_complex, Vector{ComplexF32})
    TEST_DOUBLE64 && push!(types_complex, Vector{ComplexDF64})
    TEST_STATIC_SIZED && push!(types_complex, MVector{3, ComplexF64})

    for M in manifolds
        basis_types = if M == E
            (DefaultOrthonormalBasis(), ProjectedOrthonormalBasis(:svd), DiagonalizingOrthonormalBasis([1.0, 2.0, 3.0]))
        elseif M == Ec
            (DefaultOrthonormalBasis(), DiagonalizingOrthonormalBasis([1.0, 2.0, 3.0]))
        else
            ()
        end
        for T in types
            @testset "$M Type $T" begin
                pts = [convert(T, [1.0, 0.0, 0.0]),
                       convert(T, [0.0, 1.0, 0.0]),
                       convert(T, [0.0, 0.0, 1.0])]
                test_manifold(
                    M,
                    pts,
                    test_reverse_diff = isa(T, Vector),
                    test_project_point = true,
                    test_project_tangent = true,
                    test_musical_isomorphisms = true,
                    test_vector_transport = true,
                    test_mutating_rand = isa(T, Vector),
                    point_distributions = [Manifolds.projected_distribution(M, Distributions.MvNormal(zero(pts[1]), 1.0))],
                    tvector_distributions = [Manifolds.normal_tvector_distribution(M, pts[1], 1.0)],
                    basis_types_vecs = basis_types,
                    basis_types_to_from = basis_types,
                    basis_has_specialized_diagonalizing_get = true,
                    test_vee_hat = isa(M, Euclidean),
                )
            end
        end
    end
    for T in types_complex
        @testset "Complex Euclidean, type $T" begin
            pts = [convert(T, [1.0im, -1.0im, 1.0]),
                   convert(T, [0.0, 1.0, 1.0im]),
                   convert(T, [0.0, 0.0, 1.0])]
            test_manifold(
                Ec,
                pts,
                test_reverse_diff = isa(T, Vector),
                test_project_tangent = true,
                test_musical_isomorphisms = true,
                test_vector_transport = true,
                test_vee_hat = false,
            )
        end
    end

    @testset "hat/vee" begin
        E = Euclidean(3, 2)
        p = collect(reshape(1.0:6.0, (3, 2)))
        X = collect(reshape(7.0:12.0, (3, 2)))
        @test hat(E, p, vec(X)) ≈ X
        Y = allocate(X)
        @test hat!(E, Y, p, vec(X)) === Y
        @test Y ≈ X
        @test vee(E, p, X) ≈ vec(X)
        Y = allocate(vec(X))
        @test vee!(E, Y, p, X) === Y
        @test Y ≈ vec(X)
    end

    @testset "Number systems power" begin
        @test ℝ^2 === Euclidean(2)
        @test ℝ^(2,3) === Euclidean(2, 3)

        @test ℂ^2 === Euclidean(2; field = ℂ)
        @test ℂ^(2,3) === Euclidean(2, 3; field = ℂ)

        @test ℍ^2 === Euclidean(2; field = ℍ)
        @test ℍ^(2,3) === Euclidean(2, 3; field = ℍ)
    end
end


text:
export Webhook

"""
A Webhook.
More details [here](https://discordapp.com/developers/docs/resources/webhook#webhook-object).
"""
struct Webhook
    id::Snowflake
    guild_id::Union{Snowflake, Missing}
    channel_id::Snowflake
    user::Union{User, Missing}
    name::Union{String, Nothing}
    avatar::Union{String, Nothing}
    token::String
end
@boilerplate Webhook :dict :docs :lower :merge


text:
# Test all the reduced functions

@testset "Test reduced_grad_x_grad_xd and reduced_hess_x_grad_xd with no active dimensions" begin
    Nx_tab = [1;2;4;8]
    Ne = 50
    m = 10

    for Nx in Nx_tab

        X = randn(Nx, Ne) .* randn(Nx, Ne) .+ randn(Nx)
        B = MultiBasis(CstProHermiteBasis(m-2), Nx)
        idx = zeros(Int64, 1, Nx)
        coeff = randn(size(idx,1))
        f = ExpandedFunction(B, idx, coeff)

        # Test reduced_grad_x_grad_xd
        dψ = grad_x_grad_xd(f, X)
        dψr = reduced_grad_x_grad_xd(f, X)

        @test norm(dψ[:,f.dim] - dψr)<1e-8

        # Test reduced_hess_x_grad_xd
        d2ψ = hess_x_grad_xd(f, X)
        d2ψr = reduced_hess_x_grad_xd(f, X)

        @test norm(d2ψ[:,f.dim, f.dim] - d2ψr)<1e-8

    end
end

@testset "Test reduced_grad_x_grad_xd and reduced_hess_x_grad_xd with dimensions" begin
    Nx_tab = [1;2;2;2;4;6;8]
    Ne = 50
    m = 10
    Nψ = 5

    idx_tab = (reshape([0; 1; 2; 3; 4], (5, 1)),
               [0 0; 0 1; 1 0; 1 1; 1 2],
               [0 0; 0 1; 0 2; 0 3; 0 4],
               [0 0; 1 0; 2 0; 3 0; 4 0],
                [1  0  0  0;
                 2  0  0  0;
                 0  0  1  0;
                 0  0  2  0;
                 3  0  0  0],
                [0  0  0  0  0  1;
                 1  0  2  0  3  2;
                 1  0  0  0  0  0;
                 2  0  0  0  0  0;
                 1  0  2  2  1  1],
               [0  0  0  0  0  1  0  0;
                0  0  0  0  0  0  0  1;
                0  0  0  0  0  0  0  2;
                0  0  0  0  0  2  0  0;
                0  0  1  0  0  0  0  0])

    for (i, Nx) in enumerate(Nx_tab)
        X = randn(Nx, Ne) .* randn(Nx, Ne) .+ randn(Nx)
        B = MultiBasis(CstProHermiteBasis(m-2), Nx)
        idx = idx_tab[i]

        coeff = randn(size(idx,1))
        f = ExpandedFunction(B, idx, coeff)

        # Test reduced_grad_x_grad_xd
        dψ = grad_x_grad_xd(f, X)
        dψr = reduced_grad_x_grad_xd(f, X)

        @test norm(dψ[:,f.dim] - dψr)<1e-8

        # Test reduced_hess_x_grad_xd
        d2ψ = hess_x_grad_xd(f, X)
        d2ψr = reduced_hess_x_grad_xd(f, X)

        @test norm(d2ψ[:,f.dim, f.dim] - d2ψr)<1e-8

    end
end


text:
function Base.show(io::IO, ::MIME"text/plain", flux_res::FluxVariabilitySummary)

    longest_biomass_len =
        maximum(length(k) for k in keys(flux_res.biomass_fluxes); init = 0)
    longest_exchange_len =
        maximum(length(k) for k in keys(flux_res.exchange_fluxes); init = 0)
    word_pad_len = max(longest_biomass_len, longest_exchange_len)

    longest_biomass_flux_len = maximum(
        length(string(round(v[1], digits = 4))) for v in values(flux_res.biomass_fluxes);
        init = 0,
    )
    longest_exchange_flux_len = max(
        maximum(
            length(string(round(v[1], digits = 4))) for
            v in values(flux_res.exchange_fluxes);
            init = 0,
        ),
        length("Lower bound  "),
    )
    number_pad_len = max(longest_biomass_flux_len, longest_exchange_flux_len)

    println(
        io,
        "Biomass",
        _pad_spaces(length("Biomass"), word_pad_len + 3),
        "Lower bound  ",
        _pad_spaces("Lower bound  ", number_pad_len),
        "Upper bound",
    )
    for (k, v) in flux_res.biomass_fluxes
        lb = isnothing(v[1]) ? " " : string(round(v[1], digits = 4))
        ub = isnothing(v[2]) ? " " : string(round(v[1], digits = 4))
        println(
            io,
            "  ",
            k,
            ":",
            _pad_spaces(k, word_pad_len),
            lb,
            _pad_spaces(lb, number_pad_len),
            ub,
        )
    end

    println(io, "Exchange")
    ex_ids = collect(keys(flux_res.exchange_fluxes))
    vs = [
        abs(flux_res.exchange_fluxes[k][1]) + abs(flux_res.exchange_fluxes[k][2]) for
        k in ex_ids
    ]
    idxs = sortperm(vs, rev = true)
    for k in ex_ids[idxs]
        v = flux_res.exchange_fluxes[k]
        lb = isnothing(v[1]) ? " " : string(round(v[1], digits = 4))
        ub = isnothing(v[2]) ? " " : string(round(v[2], digits = 4))
        println(
            io,
            "  ",
            k,
            ":",
            _pad_spaces(k, word_pad_len),
            lb,
            _pad_spaces(lb, number_pad_len),
            ub,
        )
    end
end


text:
using SpecialMatrices
using Compat.Test
using Compat.LinearAlgebra

include("companion.jl") #Companion matrix
include("frobenius.jl") #Frobenius matrix
include("strang.jl") #Strang matrix
include("hankel.jl") #Hankel matrix
include("hilbert.jl") #Hilbert matrix
include("toeplitz.jl") #Toeplitz matrix
include("vandermonde.jl") #Vandermonde matrix


text:
abstract type AbstractLShapedSolver{F, T <: Real, A <: AbstractVector, M <: LQSolver, S <: LQSolver} <: AbstractStructuredModel end

nscenarios(lshaped::AbstractLShapedSolver) = lshaped.nscenarios

# Initialization #
# ======================================================================== #
function init!(lshaped::AbstractLShapedSolver, subsolver::SubSolver)
    # Initialize progress meter
    lshaped.progress.thresh = lshaped.parameters.τ
    # Cap bundle size
    lshaped.parameters.bundle = min(lshaped.parameters.bundle, lshaped.nscenarios)
    # Prepare the master optimization problem
    prepare_master!(lshaped)
    # Finish initialization based on solver traits
    init_subproblems!(lshaped, subsolver)
    init_solver!(lshaped)
end
# ======================================================================== #

# Functions #
# ======================================================================== #
function set_params!(lshaped::AbstractLShapedSolver; kwargs...)
    for (k,v) in kwargs
        setfield!(lshaped.parameters,k,v)
    end
end

function update_solution!(lshaped::AbstractLShapedSolver)
    ncols = decision_length(lshaped.stochasticprogram)
    nb = nbundles(lshaped)
    x = getsolution(lshaped.mastersolver)
    lshaped.mastervector[:] = x[1:ncols+nb]
    lshaped.x[1:ncols] = x[1:ncols]
    lshaped.θs[:] = x[ncols+1:ncols+nb]
    nothing
end

function calculate_estimate(lshaped::AbstractLShapedSolver)
    return lshaped.c⋅lshaped.x + sum(lshaped.θs)
end

function current_objective_value(lshaped::AbstractLShapedSolver,Qs::AbstractVector)
    return lshaped.c⋅lshaped.x + sum(Qs)
end
current_objective_value(lshaped) = current_objective_value(lshaped,lshaped.subobjectives)

function get_objective_value(lshaped::AbstractLShapedSolver)
    if !isempty(lshaped.Q_history)
        return lshaped.Q_history[end]
    else
        return calculate_objective_value(lshaped)
    end
end

function prepare_master!(lshaped::AbstractLShapedSolver{true})
    # θs
    for i = 1:nbundles(lshaped)
        MPB.addvar!(lshaped.mastersolver.lqmodel,-Inf,Inf,0.0)
        if typeof(lshaped.mastersolver.optimsolver) == GurobiSolver
            updatemodel!(lshaped.mastersolver.lqmodel)
        end
        push!(lshaped.mastervector,-1e10)
        push!(lshaped.θs,-1e10)
    end
end

function prepare_master!(lshaped::AbstractLShapedSolver{false})
    # θs
    for i = 1:nbundles(lshaped)
        MPB.addvar!(lshaped.mastersolver.lqmodel,-Inf,Inf,1.0)
        if typeof(lshaped.mastersolver.optimsolver) == GurobiSolver
            updatemodel!(lshaped.mastersolver.lqmodel)
        end
        push!(lshaped.mastervector,-1e10)
        push!(lshaped.θs,-1e10)
    end
end

function resolve_subproblems!(lshaped::AbstractLShapedSolver{F,T}) where {F, T <: Real}
    # Update subproblems
    update_subproblems!(lshaped.subproblems,lshaped.x)
    # Solve sub problems
    cut_bundle = CutBundle(T)
    for subproblem ∈ lshaped.subproblems
        cut::SparseHyperPlane{T} = subproblem()
        if lshaped.parameters.bundle == 1
            add_cut!(lshaped,cut)
            update_objective!(lshaped,cut)
        else
            add_to_bundle!(lshaped,cut_bundle,cut)
            if length(cut_bundle) == lshaped.parameters.bundle
                aggregated_cut = aggregate!(cut_bundle)
                add_cut!(lshaped,aggregated_cut)
                update_objective!(lshaped,aggregated_cut)
                lshaped.subobjectives[aggregated_cut.id] = cut_bundle.q
                cut_bundle.q = zero(T)
            end
        end
    end
    if lshaped.parameters.bundle > 1 && length(cut_bundle) > 0 && cut_bundle.q < Inf
        # Add remaining bundle
        aggregated_cut = aggregate!(cut_bundle)
        add_cut!(lshaped,aggregated_cut)
        update_objective!(lshaped,aggregated_cut)
        lshaped.subobjectives[aggregated_cut.id] = cut_bundle.q
    end
    # Return current objective value
    return current_objective_value(lshaped)
end

function iterate_nominal!(lshaped::AbstractLShapedSolver{F}) where F
    # Resolve all subproblems at the current optimal solution
    Q = resolve_subproblems!(lshaped)
    if Q == Inf && !F
        return :Infeasible
    end
    if Q == -Inf
        return :Unbounded
    end
    lshaped.solverdata.Q = Q
    # Update the optimization vector
    take_step!(lshaped)
    # Resolve master
    try
        solve_problem!(lshaped,lshaped.mastersolver)
    catch
        # Master problem could not be solved for some reason.
        @unpack Q,θ = lshaped.solverdata
        gap = abs(θ-Q)/(abs(Q)+1e-10)
        @warn "Master problem could not be solved, solver returned status $(status(lshaped.mastersolver)). The following relative tolerance was reached: $(@sprintf("%.1e",gap)). Aborting procedure."
        return :StoppedPrematurely
    end
    if status(lshaped.mastersolver) == :Infeasible
        @warn "Master is infeasible. Aborting procedure."
        return :Infeasible
    end
    # Update master solution
    update_solution!(lshaped)
    lshaped.solverdata.θ = calculate_estimate(lshaped)
    # Log progress
    log!(lshaped)
    # Check optimality
    if check_optimality(lshaped)
        # Optimal
        lshaped.solverdata.Q = calculate_objective_value(lshaped,lshaped.x)
        push!(lshaped.Q_history,lshaped.solverdata.Q)
        return :Optimal
    end
    # Project (if applicable)
    project!(lshaped)
    # Check optimality if level sets are used
    if hastrait(lshaped,LV)
        lshaped.solverdata.θ = calculate_estimate(lshaped)
        if check_optimality(lshaped)
            # Optimal
            lshaped.solverdata.Q = calculate_objective_value(lshaped,lshaped.x)
            push!(lshaped.Q_history,lshaped.solverdata.Q)
            # Final log
            if lshaped.parameters.log
                current_gap = gap(lshaped)
                ProgressMeter.update!(lshaped.progress,current_gap,
                                      showvalues = [
                                          ("Objective",Q),
                                          ("Gap",current_gap),
                                          ("Number of cuts",length(lshaped.cuts))
                                      ])
            end
            return :Optimal
        end
    end
    # Just return a valid status for this iteration
    return :Valid
end

function log!(lshaped::AbstractLShapedSolver)
    @unpack Q,θ = lshaped.solverdata
    push!(lshaped.Q_history,Q)
    push!(lshaped.θ_history,θ)
    lshaped.solverdata.iterations += 1

    log_regularization!(lshaped)

    if lshaped.parameters.log
        current_gap = gap(lshaped)
        ProgressMeter.update!(lshaped.progress,current_gap,
                              showvalues = [
                                  ("Objective",Q),
                                  ("Gap",current_gap),
                                  ("Number of cuts",length(lshaped.cuts))
                              ])
    end
end

function log!(lshaped::AbstractLShapedSolver, t::Integer)
    @unpack Q,θ = lshaped.solverdata
    lshaped.Q_history[t] = Q
    lshaped.θ_history[t] = θ

    log_regularization!(lshaped,t)

    if lshaped.parameters.log
        current_gap = gap(lshaped)
        ProgressMeter.update!(lshaped.progress,current_gap,
                              showvalues = [
                                  ("Objective",Q),
                                  ("Gap",current_gap),
                                  ("Number of cuts",length(lshaped.cuts))
                              ])
    end
end

function check_optimality(lshaped::AbstractLShapedSolver)
    @unpack τ = lshaped.parameters
    @unpack θ = lshaped.solverdata
    return θ > -Inf && gap(lshaped) <= τ
end
# ======================================================================== #

# Cut functions #
# ======================================================================== #
active(lshaped::AbstractLShapedSolver, hyperplane::HyperPlane) = active(hyperplane, lshaped.x, lshaped.parameters.τ)
active(lshaped::AbstractLShapedSolver, cut::HyperPlane{OptimalityCut}) = optimal(cut, lshaped.x, lshaped.θs[cut.id], lshaped.parameters.τ)
satisfied(lshaped::AbstractLShapedSolver, hyperplane::HyperPlane) = satisfied(hyperplane, lshaped.x, lshaped.parameters.τ)
satisfied(lshaped::AbstractLShapedSolver, cut::HyperPlane{OptimalityCut}) = satisfied(cut, lshaped.x, lshaped.θs[cut.id], lshaped.parameters.τ)
violated(lshaped::AbstractLShapedSolver, hyperplane::HyperPlane) = !satisfied(lshaped, hyperplane)
gap(lshaped::AbstractLShapedSolver, hyperplane::HyperPlane) = gap(hyperplane, lshaped.x)
gap(lshaped::AbstractLShapedSolver, cut::HyperPlane{OptimalityCut}) = gap(cut, lshaped.x, lshaped.θs[cut.id])

add_cut!(lshaped::AbstractLShapedSolver, cut::HyperPlane) = add_cut!(lshaped, cut, lshaped.subobjectives)

function add_cut!(lshaped::AbstractLShapedSolver, cut::HyperPlane{OptimalityCut}, subobjectives::AbstractVector, Q::Real)
    θ = lshaped.θs[cut.id]
    @unpack τ = lshaped.parameters
    # Update objective
    subobjectives[cut.id] = Q
    # Check if cut gives new information
    if θ > -Inf && abs(θ-Q) <= τ*(1+abs(Q))
        # Optimal with respect to this subproblem
        return false
    end
    # Add optimality cut
    process_cut!(lshaped, cut)
    MPB.addconstr!(lshaped.mastersolver.lqmodel, lowlevel(cut)...)
    if typeof(lshaped.mastersolver.optimsolver) == GurobiSolver
        updatemodel!(lshaped.mastersolver.lqmodel)
    end
    push!(lshaped.cuts, cut)
    return true
end
add_cut!(lshaped::AbstractLShapedSolver, cut::HyperPlane{OptimalityCut}, subobjectives::AbstractVector, x::AbstractVector) = add_cut!(lshaped, cut,subobjectives, cut(x))
add_cut!(lshaped::AbstractLShapedSolver, cut::HyperPlane{OptimalityCut}, subobjectives::AbstractVector) = add_cut!(lshaped, cut, subobjectives, lshaped.x)

function add_cut!(lshaped::AbstractLShapedSolver, cut::HyperPlane{FeasibilityCut}, subobjectives::AbstractVector, Q::Real)
    # Ensure that there is no false convergence
    subobjectives[cut.id] = Q
    # Add feasibility cut
    process_cut!(lshaped, cut)
    MPB.addconstr!(lshaped.mastersolver.lqmodel, lowlevel(cut)...)
    if typeof(lshaped.mastersolver.optimsolver) == GurobiSolver
        updatemodel!(lshaped.mastersolver.lqmodel)
    end
    push!(lshaped.cuts, cut)
    return true
end
add_cut!(lshaped::AbstractLShapedSolver, cut::HyperPlane{FeasibilityCut}, subobjectives::AbstractVector) = add_cut!(lshaped, cut, subobjectives, Inf)

function add_cut!(lshaped::AbstractLShapedSolver, cut::HyperPlane{Infeasible}, subobjectives::AbstractVector)
    @warn "Subproblem $cut.id is infeasible, procedure will abort."
    subobjectives[cut.id] = Inf
    return true
end

function add_cut!(lshaped::AbstractLShapedSolver, cut::HyperPlane{Unbounded}, subobjectives::AbstractVector)
    @warn "Subproblem $cut.id is unbounded, procedure will abort."
    subobjectives[cut.id] = -Inf
    return true
end

function add_to_bundle!(lshaped::AbstractLShapedSolver, bundle::CutBundle, cut::HyperPlane)
    add_cut!(lshaped, cut)
    bundle.q += cut(lshaped.x)
end

function add_to_bundle!(lshaped::AbstractLShapedSolver, bundle::CutBundle, cut::HyperPlane{OptimalityCut})
    push!(bundle.cuts, cut)
    bundle.q += cut(lshaped.x)
end

update_objective!(lshaped::AbstractLShapedSolver, cut::HyperPlane) = nothing
function update_objective!(lshaped::AbstractLShapedSolver{true}, cut::HyperPlane{OptimalityCut})
    # Ensure that θi is included in minimization if feasibility cuts are used
    c = MPB.getobj(lshaped.mastersolver.lqmodel)
    if c[length(lshaped.x)+cut.id] == 0.0
        c[length(lshaped.x)+cut.id] = 1.0
        MPB.setobj!(lshaped.mastersolver.lqmodel,c)
    end
end

function show(io::IO, lshaped::AbstractLShapedSolver)
    println(io, typeof(lshaped).name.name)
    println(io, "State:")
    show(io, lshaped.solverdata)
    println(io, "Parameters:")
    show(io, lshaped.parameters)
end

function show(io::IO, ::MIME"text/plain", lshaped::AbstractLShapedSolver)
    show(io, lshaped)
end
# ======================================================================== #

# Plot recipe #
# ======================================================================== #
@recipe f(lshaped::AbstractLShapedSolver) = lshaped, -1
@recipe function f(lshaped::AbstractLShapedSolver, time::Real; showθ = false)
    length(lshaped.Q_history) > 0 || error("No solution data. Has solver been run?")
    showθ && (length(lshaped.θ_history) > 0 || error("No solution data. Has solver been run?"))
    Qmin = showθ ? minimum(lshaped.θ_history) : minimum(lshaped.Q_history)
    Qmax = maximum(lshaped.Q_history)
    increment = std(lshaped.Q_history)

    linewidth --> 4
    linecolor --> :black
    tickfontsize := 14
    tickfontfamily := "sans-serif"
    guidefontsize := 16
    guidefontfamily := "sans-serif"
    titlefontsize := 22
    titlefontfamily := "sans-serif"
    xlabel := time == -1 ? "Iteration" : "Time [s]"
    ylabel := "Q"
    ylims --> (Qmin-increment,Qmax+increment)
    if time == -1
        xlims --> (1, length(lshaped.Q_history)+1)
        xticks --> 1:5:length(lshaped.Q_history)
    else
        xlims --> (0, time)
        xticks --> linspace(0, time, ceil(Int,length(lshaped.Q_history)/5))
    end
    yticks --> Qmin:increment:Qmax
    xformatter := (d) -> @sprintf("%.1f",d)
    yformatter := (d) -> begin
        if abs(d) <= sqrt(eps())
            "0.0"
        elseif (log10(abs(d)) < -2.0 || log10(abs(d)) > 3.0)
            @sprintf("%.4e",d)
        elseif log10(abs(d)) > 2.0
            @sprintf("%.1f",d)
        else
            @sprintf("%.2f",d)
        end
    end

    @series begin
        label --> "Q"
        seriescolor --> :black
        if time == -1
            1:1:length(lshaped.Q_history), lshaped.Q_history
        else
            linspace(0,time,length(lshaped.Q_history)), lshaped.Q_history
        end
    end

    if showθ
        @series begin
            label --> "θ"
            linestyle --> :dash
            seriescolor --> :red
            linecolor --> :red
            linewidth --> 2
            if time == -1
                1:1:length(lshaped.θ_history), lshaped.θ_history
            else
                linspace(0,time,length(lshaped.θ_history)), lshaped.θ_history
            end
        end
    end
end
# ======================================================================== #


text:
println("===============================================")
println("Testing definition.jl")
println("===============================================")

function _RheoTimeData_explicit_nolog()
    a = Vector{RheoFloat}(1.0:1.0:3.0)
    b = Vector{RheoFloat}(4.0:1.0:6.0)
    c = Vector{RheoFloat}(7.0:1.0:9.0)

    data = RheoTimeData(a, b, c, nothing)

    data.σ==a && data.ϵ==b && data.t==c && isnothing(data.log)
end
@test _RheoTimeData_explicit_nolog()

function _RheoTimeData_const_nolog()
    eps = Vector{RheoFloat}(1.0:1.0:3.0)
    sig = Vector{RheoFloat}(4.0:1.0:6.0)
    t = Vector{RheoFloat}(7.0:1.0:9.0)

    data = RheoTimeData(strain=eps, stress=sig, t=t, savelog=false)
    showlog(data)

    getstress(data)==sig && getstrain(data)==eps && gettime(data)==t && isnothing(data.log)
end
@test _RheoTimeData_const_nolog()

function _RheoTimeData_const_nostrain()
    sig = Vector{RheoFloat}(4.0:1.0:6.0)
    t = Vector{RheoFloat}(7.0:1.0:9.0)

    data = RheoTimeData(σ=sig, t=t)
    showlog(data)

    data.σ==sig && data.ϵ==RheoFloat[] && data.t==t && data.log[1].info.type == stress_only::TimeDataType
end
@test _RheoTimeData_const_nostrain()

function _operators_logs()
    d1=strainfunction(timeline(0:0.5:10),t->exp(-t))
    d2=strainfunction(timeline(0:0.5:10),t->1-exp(-t))
    d=2*d1 - (-d2) + d2

    d3=rheologrun(d.log)

    (d3.ϵ == d.ϵ) && all([ abs(e-2.)<=eps(RheoFloat) for e in d.ϵ ])
end
@test _operators_logs()

function _RheoFreqData_explicit_nolog()
    a = Vector{RheoFloat}(1.0:1.0:3.0)
    b = Vector{RheoFloat}(4.0:1.0:6.0)
    c = Vector{RheoFloat}(7.0:1.0:9.0)

    data = RheoFreqData(a, b, c, nothing)

    data.Gp==a && data.Gpp==b && data.ω==c && isnothing(data.log)
end
@test _RheoFreqData_explicit_nolog()

function _RheoFreqData_const_nolog()
    a = Vector{RheoFloat}(1.0:1.0:3.0)
    b = Vector{RheoFloat}(4.0:1.0:6.0)
    c = Vector{RheoFloat}(7.0:1.0:9.0)

    data = RheoFreqData(omega=a, Gp=b, Gpp=c , savelog=false)
    showlog(data)

    getfreq(data)==a && getstorage(data)==b && getloss(data)==c && isnothing(data.log)
end
@test _RheoFreqData_const_nolog()

function _rheoconvert()
    vi64=Int64(1)
    ai64=[vi64,vi64]
    arf=Vector{RheoFloat}([1,2,3])
    rheoconvert(RheoFloat(1.0))===rheoconvert(vi64) && typeof(rheoconvert(ai64)) == Vector{RheoFloat} &&
        arf === rheoconvert(arf)  &&  ai64 !== rheoconvert(ai64)
end
@test _rheoconvert()

function _union_strain1stress2()
    time_instance = timeline(t_start=0.0, t_end=15.0, step=0.2)
    imposed_strain = strainfunction(time_instance, sawtooth(offset=5.0, amp=2, period=5))
    imposed_stress = stressfunction(time_instance, stairs(offset=3.0, amp=0.5, width=4.0))

    combined = imposed_strain|imposed_stress

    (combined.σ==imposed_stress.σ) && (combined.ϵ==imposed_strain.ϵ) && (combined.t==imposed_stress.t)
end
@test _union_strain1stress2()

function _union_stress1strain2()
    time_instance = timeline(t_start=0.0, t_end=15.0, step=0.2)
    imposed_strain = strainfunction(time_instance, sawtooth(offset=5.0, amp=2, period=5))
    imposed_stress = stressfunction(time_instance, stairs(offset=3.0, amp=0.5, width=4.0))

    combined = imposed_stress|imposed_strain

    (combined.σ==imposed_stress.σ) && (combined.ϵ==imposed_strain.ϵ) && (combined.t==imposed_stress.t)
end
@test _union_stress1strain2()

function _freeze_params()
    SLS2_mod = freeze_params( SLS2, G₀=2, η₂=3.5)

    relaxmod(SLS2, 1, [2,1,2,3,3.5]) == relaxmod(SLS2_mod, 1, [1,2,3])
end
@test _freeze_params()

function _getparams()
    m = RheoModel(Maxwell, k=1, eta=2)
    p1 = getparams(m)
    p2 = getparams(m,unicode=false)
    (p1.k==1.0) && (p1.η==2.0) &&  (p2.k==1.0) && (p2.eta==2.0)
end
@test _getparams()

function _getparams_dict()
    m = RheoModel(Maxwell, k=1, eta=2)
    p = dict(getparams(m,unicode=false))
    (p[:k]==1.0) && (p[:eta]==2.0)
end
@test _getparams_dict()

#
#   The moduli functions on arrays are extensively tested as part of the model tests
#

function _scalar_moduli()
    m=RheoModel(Spring, k=2)
    (relaxmod(m))(1) == relaxmod(m, 1) == relaxmod(Spring, k=2, 1) == 2. &&
        (creepcomp(m))(1) == creepcomp(m, 1) == creepcomp(Spring, k=2, 1) == 0.5 &&
        (storagemod(m))(1) == storagemod(m, 1) == storagemod(Spring, k=2, 1) == 2. &&
        (lossmod(m))(1) == lossmod(m, 1) == lossmod(Spring, k=2, 1) == 0.0 &&
        (dynamicmod(m))(1) == dynamicmod(m, 1) == dynamicmod(Spring, k=2, 1) == 2.0 + 0.0*im
end
@test _scalar_moduli()


text:
using Images

println(pwd())
fn_img = "./examples/rakuten.png"
img = load(fn_img)
@show img

max_stars_repo_name: Soyukke/FindDifferences

max_issues_repo_name: Soyukke/FindDifferences

max_forks_repo_name: Soyukke/FindDifferences


text:
using Documenter
using LatticeRules

makedocs(;
    modules=[LatticeRules],
    authors="Pieterjan Robbe <pieterjan.robbe@kuleuven.be> and contributors",
    repo="https://github.com/PieterjanRobbe/LatticeRules.jl/blob/{commit}{path}#L{line}",
    sitename="LatticeRules.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://PieterjanRobbe.github.io/LatticeRules.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/PieterjanRobbe/LatticeRules.jl",
)


text:
using Manopt, Manifolds, ManifoldsBase, LinearAlgebra, Test

@testset "Conjugate Gradient coefficient rules" begin
    F(x) = norm(x)^2
    ∇F(x) = 2*x
    M = Euclidean(2)
    P = GradientProblem(M,F,∇F)
    x0 = [0.0, 1.0]
    sC = StopAfterIteration(1)
    s = ConstantStepsize(1.0)
    retr = ExponentialRetraction()
    vtm = ParallelTransport()

    ∇1 = [1.0, 1.0]
    δ1 = [0.0, 2.0]
    ∇2 = [1.0, 1.5]
    δ2 = [0.5, 2.0]
    diff = ∇2-∇1

    dU = SteepestDirectionUpdateRule()
    O = ConjugateGradientDescentOptions(x0, sC, s, dU, retr, vtm)
    @test O.coefficient(P,O,1) == 0

    dU = ConjugateDescentCoefficient()
    O = ConjugateGradientDescentOptions(x0, sC, s, dU, retr, vtm)
    O.∇ = ∇1
    O.δ = δ1
    # for the first case we get zero
    @test O.coefficient(P,O,1) == 0.0
    O.∇ = ∇2
    O.δ = δ2
    @test O.coefficient(P,O,2) == dot(∇2,∇2) / dot(-δ2, ∇1)

    dU = DaiYuanCoefficient()
    O = ConjugateGradientDescentOptions(x0, sC, s, dU, retr, vtm)
    O.∇ = ∇1
    O.δ = δ1
    # for the first case we get zero
    @test O.coefficient(P,O,1) == 0.0
    O.∇ = ∇2
    O.δ = δ2
    @test O.coefficient(P,O,2) == dot(∇2,∇2) / dot(δ2, ∇2-∇1)

    dU = FletcherReevesCoefficient()
    O = ConjugateGradientDescentOptions(x0, sC, s, dU, retr, vtm)
    O.∇ = ∇1
    O.δ = δ1
    # for the first case we get zero
    @test O.coefficient(P,O,1) == 1.0
    O.∇ = ∇2
    O.δ = δ2
    @test O.coefficient(P,O,2) == dot(∇2,∇2) / dot(∇1,∇1)

    dU = HagerZhangCoefficient()
    O = ConjugateGradientDescentOptions(x0, sC, s, dU, retr, vtm)
    O.∇ = ∇1
    O.δ = δ1
    # for the first case we get zero
    @test O.coefficient(P,O,1) == 0.0
    O.∇ = ∇2
    O.δ = δ2
    denom = dot(δ1, diff)
    ndiffsq = dot(diff,diff)
    @test O.coefficient(P,O,2) == dot(diff,∇2)/denom - 2*ndiffsq*dot(δ1,∇2)/denom^2

    dU = HeestenesStiefelCoefficient()
    O = ConjugateGradientDescentOptions(x0, sC, s, dU, retr, vtm)
    O.∇ = ∇1
    O.δ = δ1
    @test O.coefficient(P,O,1) == 0.0
    O.∇ = ∇2
    O.δ = δ2
    @test O.coefficient(P,O,2) == dot(diff,∇2)/dot(δ1,diff)

    dU = LiuStoreyCoefficient()
    O = ConjugateGradientDescentOptions(x0, sC, s, dU, retr, vtm)
    O.∇ = ∇1
    O.δ = δ1
    @test O.coefficient(P,O,1) == 0.0
    O.∇ = ∇2
    O.δ = δ2
    @test O.coefficient(P,O,2) == -dot(∇2,diff)/dot(δ1,∇1)

    dU = PolakRibiereCoefficient()
    O = ConjugateGradientDescentOptions(x0, sC, s, dU, retr, vtm)
    O.∇ = ∇1
    O.δ = δ1
    @test O.coefficient(P,O,1) == 0.0
    O.∇ = ∇2
    O.δ = δ2
    @test O.coefficient(P,O,2) == dot(∇2,diff)/dot(∇1,∇1)
end
@testset "Conjugate Gradient runs – Low Rank matrix approx" begin
    A = Diagonal([2.0,1.1,1.0])
    M = Sphere(size(A,1)-1)
    F(x) = x'*A*x
    euclidean_∇F(x) = 2*A*x
    ∇F(x) = project(M,x,euclidean_∇F(x))

    x0 = [2.0,0.0,2.0]/sqrt(8.)
    xOpt = conjugate_gradient_descent(
        M,
        F,
        ∇F,
        x0;
        stepsize=ArmijoLinesearch(),
        coefficient=FletcherReevesCoefficient(),
        stopping_criterion=StopAfterIteration(15),
        #debug=[:Iteration," | ",:Cost," | ",:Iterate," |\n"],
    )
    @test isapprox(F(xOpt),minimum(eigvals(A)), atol=1.2*10^-4)
    @test isapprox(xOpt, eigvecs(A)[:,size(A,1)], atol=1.7*10^-2)
end


text:
"""
    YASStyle()

Formatting style based on [YASGuide](https://github.com/jrevels/YASGuide)
and [JuliaFormatter#198](https://github.com/domluna/JuliaFormatter.jl/issues/198).

Configurable options with different defaults to [`DefaultStyle`](@ref) are:
- `always_for_in` = true
- `whitespace_ops_in_indices` = true
- `remove_extra_newlines` = true
- `import_to_using` = true
- `pipe_to_function_call` = true
- `short_to_long_function_def` = true
- `always_use_return` = true
- `whitespace_in_kwargs` = false
"""
struct YASStyle <: AbstractStyle
    innerstyle::Union{Nothing,AbstractStyle}
end
YASStyle() = YASStyle(nothing)
@inline getstyle(s::YASStyle) = s.innerstyle === nothing ? s : s.innerstyle

function options(style::YASStyle)
    return (;
        always_for_in = true,
        whitespace_ops_in_indices = true,
        remove_extra_newlines = true,
        import_to_using = true,
        pipe_to_function_call = true,
        short_to_long_function_def = true,
        always_use_return = true,
        whitespace_in_kwargs = false,
    )
end

function nestable(::YASStyle, cst::CSTParser.EXPR)
    (CSTParser.defines_function(cst) || is_assignment(cst)) && return false
    (cst[2].val == "=>" || cst[2].val == "->") && return false
    return true
end

function p_import(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = p_import(DefaultStyle(ys), cst, s)
    idx = findfirst(n -> n.typ === PLACEHOLDER, t.nodes)
    if idx !== nothing && is_colon(t[idx-1])
        t[idx] = Whitespace(1)
    end
    t
end

function p_using(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = p_import(ys, cst, s)
    t.typ = Using
    t
end

function p_export(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = p_import(ys, cst, s)
    t.typ = Export
    t
end

function p_curly(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = FST(Curly, cst, nspaces(s))
    for (i, a) in enumerate(cst)
        n = pretty(style, a, s)
        if CSTParser.is_comma(a) && i == length(cst) - 1
            continue
        elseif CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Placeholder(0), s)
        else
            add_node!(t, n, s, join_lines = true)
        end
    end
    t
end

function p_braces(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = p_curly(ys, cst, s)
    t.typ = Braces
    t
end

function p_tuple(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = FST(TupleN, cst, nspaces(s))
    for (i, a) in enumerate(cst)
        if CSTParser.is_comma(a) && i + 1 == length(cst)
            n = pretty(style, a, s)
            if length(cst.args) == 1
                add_node!(t, n, s, join_lines = true)
            elseif !is_closer(cst[i+1])
                add_node!(t, n, s, join_lines = true)
                add_node!(t, Placeholder(1), s)
            end
        elseif CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, pretty(style, a, s), s, join_lines = true)
            add_node!(t, Placeholder(1), s)
        elseif is_binary(a) && a[2].val == "="
            add_node!(t, p_kw(style, a, s), s, join_lines = true)
        else
            add_node!(t, pretty(style, a, s), s, join_lines = true)
        end
    end
    t
end

function p_tuple(ys::YASStyle, nodes::Vector{CSTParser.EXPR}, s::State)
    style = getstyle(ys)
    t = FST(TupleN, nspaces(s))
    for (i, a) in enumerate(nodes)
        n = pretty(style, a, s)
        if CSTParser.is_comma(a) && i + 1 == length(nodes)
            if length(nodes) == 1
                add_node!(t, n, s, join_lines = true)
            elseif !is_closer(nodes[i+1])
                add_node!(t, n, s, join_lines = true)
                add_node!(t, Placeholder(1), s)
            end
        elseif CSTParser.is_comma(a) && i < length(nodes) && !is_punc(nodes[i+1])
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Placeholder(1), s)
        else
            add_node!(t, n, s, join_lines = true)
        end
    end
    t
end

# Brackets
function p_invisbrackets(
    ys::YASStyle,
    cst::CSTParser.EXPR,
    s::State;
    nonest = false,
    nospace = false,
)
    style = getstyle(ys)
    t = p_invisbrackets(DefaultStyle(style), cst, s, nonest = nonest, nospace = nospace)
    for (i, n) in enumerate(t.nodes)
        if n.typ === PLACEHOLDER
            t[i] = Whitespace(length(n))
        end
    end
    t
end

function p_call(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = FST(Call, cst, nspaces(s))
    for (i, a) in enumerate(cst)
        n = pretty(style, a, s)
        if CSTParser.is_comma(a) && i + 1 == length(cst)
            continue
        elseif CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Placeholder(1), s)
        else
            add_node!(t, n, s, join_lines = true)
        end
    end
    if !parent_is(cst, is_function_or_macro_def)
        separate_kwargs_with_semicolon!(t)
    end
    t
end
function p_vect(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = p_call(ys, cst, s)
    t.typ = Vect
    t
end

function p_vcat(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = FST(Vcat, nspaces(s))
    st = cst.head === :vcat ? 1 : 2

    for (i, a) in enumerate(cst)
        n = pretty(style, a, s)
        if !is_closer(a) && i > st
            join_lines = i == st + 1 ? true : t.endline == n.startline
            if join_lines && i != st + 1
                add_node!(t, Placeholder(1), s)
            end

            add_node!(t, n, s; join_lines = join_lines)
            if i != length(cst) - 1
                has_semicolon(s.doc, n.startline) &&
                    add_node!(t, InverseTrailingSemicolon(), s)
            end
        else
            add_node!(t, n, s, join_lines = true)
        end
    end
    t
end
function p_typedvcat(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = p_vcat(ys, cst, s)
    t.typ = TypedVcat
    t
end

function p_ref(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = FST(RefN, cst, nspaces(s))
    nospace = !s.opts.whitespace_ops_in_indices

    for (i, a) in enumerate(cst)
        if CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, pretty(style, a, s), s, join_lines = true)
            add_node!(t, Placeholder(1), s)
        elseif is_binary(a) || is_chain(a) || a.head === :brackets || a.head === :comparison
            n = pretty(style, a, s, nonest = true, nospace = nospace)
            add_node!(t, n, s, join_lines = true)
        else
            add_node!(t, pretty(style, a, s), s, join_lines = true)
        end
    end
    t
end

function p_comprehension(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = FST(Comprehension, cst, nspaces(s))

    if is_block(cst[2]) || (cst[2].head === :Generator && is_block(cst[2][1]))
        t.nest_behavior = AlwaysNest
    end

    add_node!(t, pretty(style, cst[1], s), s, join_lines = true)
    add_node!(t, pretty(style, cst[2], s), s, join_lines = true)
    add_node!(t, pretty(style, cst[3], s), s, join_lines = true)
    t
end

function p_typedcomprehension(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = FST(TypedComprehension, cst, nspaces(s))

    if is_block(cst[3]) || (cst[3].head === :generator && is_block(cst[3][1]))
        t.nest_behavior = AlwaysNest
    end

    add_node!(t, pretty(style, cst[1], s), s, join_lines = true)
    add_node!(t, pretty(style, cst[2], s), s, join_lines = true)
    add_node!(t, pretty(style, cst[3], s), s, join_lines = true)
    add_node!(t, pretty(style, cst[4], s), s, join_lines = true)
    t
end

function p_macrocall(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    fst = p_macrocall(DefaultStyle(style), cst, s)
    fst.typ == MacroCall || return fst
    is_closer(cst[end]) || return fst

    # remove initial and last placeholders
    # @call(PLACEHOLDER, args..., PLACEHOLDER) -> @call(args...)
    idx = findfirst(n -> n.typ === PLACEHOLDER, fst.nodes)
    idx !== nothing && (fst[idx] = Whitespace(0))
    idx = findlast(n -> n.typ === PLACEHOLDER, fst.nodes)
    idx !== nothing && (fst[idx] = Whitespace(0))

    return fst
end

function p_whereopcall(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = FST(Where, cst, nspaces(s))

    add_node!(t, pretty(style, cst[1], s), s)
    add_node!(t, Whitespace(1), s)
    add_node!(t, pretty(style, cst[2], s), s, join_lines = true)
    add_node!(t, Whitespace(1), s)

    add_braces =
        !CSTParser.is_lbrace(cst[3]) &&
        cst.parent.head !== :curly &&
        cst[3].head !== :curly &&
        cst[3].head !== :bracescat

    if add_braces
        brace = FST(PUNCTUATION, -1, t.endline, t.endline, "{")
        add_node!(t, brace, s, join_lines = true)
    end

    for i = 3:length(cst)
        a = cst[i]
        n = is_binary(a) ? pretty(style, a, s, nospace = true) : pretty(style, a, s)

        if CSTParser.is_comma(a) && i == length(cst) - 1
            continue
        elseif CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Placeholder(0), s)
        else
            add_node!(t, n, s, join_lines = true)
        end
    end

    if add_braces
        brace = FST(PUNCTUATION, -1, t.endline, t.endline, "}")
        add_node!(t, brace, s, join_lines = true)
    end
    t
end

function p_generator(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = FST(Generator, cst, nspaces(s))
    has_for_kw = false

    for (i, a) in enumerate(cst)
        n = pretty(style, a, s)
        if CSTParser.iskeyword(a)
            if !has_for_kw && a.head === :FOR
                has_for_kw = true
            end

            incomp =
                parent_is(a, is_iterable, ignore = n -> is_gen(n) || n.head === :brackets)

            if is_block(cst[i-1])
                add_node!(t, Newline(), s)
            elseif incomp
                add_node!(t, Placeholder(1), s)
            else
                add_node!(t, Whitespace(1), s)
            end
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Whitespace(1), s)

            if !is_gen(cst[i+1])
                tupargs = CSTParser.EXPR[]
                for j = i+1:length(cst)
                    push!(tupargs, cst[j])
                end

                # verify this is not another for loop
                any(b -> b.head === :FOR, tupargs) && continue

                tup = p_tuple(style, tupargs, s)
                add_node!(t, tup, s, join_lines = true)

                if has_for_kw
                    for nn in tup.nodes
                        eq_to_in_normalization!(nn, s.opts.always_for_in)
                    end
                end
                break
            end
        elseif CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Placeholder(1), s)
        else
            add_node!(t, n, s, join_lines = true)
        end

        has_for_kw && eq_to_in_normalization!(n, s.opts.always_for_in)
    end

    t
end

function p_filter(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = p_generator(ys, cst, s)
    t.typ = Filter
    t
end

function p_flatten(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = p_generator(ys, cst, s)
    t.typ = Flatten
    t
end


text:
module Preconditioners

using AlgebraicMultigrid, Compat, LimitedLDLFactorizations
const AMG = AlgebraicMultigrid
const LLDL = LimitedLDLFactorizations

using LinearAlgebra, SparseArrays
import LinearAlgebra: ldiv!, \, *, mul!

abstract type AbstractPreconditioner end

function get_data(A)
    if A isa Symmetric || A isa Hermitian
        @warn("Using the data field of the symmetric/Hermitian matrix input.")
        return A.data
    else
        return A
    end
end

include("incompletecholesky.jl")
include("diagonal.jl")
include("amg.jl")

export  CholeskyPreconditioner,
        EmptyCholeskyPreconditioner,
        DiagonalPreconditioner,
        EmptyDiagonalPreconditioner,
        AMGPreconditioner,
        SmoothedAggregation,
        RugeStuben,
        UpdatePreconditioner!

end # module


text:
function check_data_file_hash(path, expected_checksum)
    actual_checksum = bytes2hex(sha256(open(path)))
    if actual_checksum != expected_checksum
        @error "The downloaded data file `$path' seems to be different from the expected one. Tests will likely fail." actual_checksum expected_checksum
    end
end

function download_data_file(url, path, hash)
    if isfile(path)
        check_data_file_hash(path, hash)
        @info "using cached `$path'"
        return path
    end

    Downloads.download(url, path)
    check_data_file_hash(path, hash)
    return path
end

isdir("downloaded") || mkdir("downloaded")
df(s) = joinpath("downloaded", s)

model_paths = Dict{String,String}(
    "iJO1366.json" => download_data_file(
        "http://bigg.ucsd.edu/static/models/iJO1366.json",
        df("iJO1366.json"),
        "9376a93f62ad430719f23e612154dd94c67e0d7c9545ed9d17a4d0c347672313",
    ),
    "iJO1366.mat" => download_data_file(
        "http://bigg.ucsd.edu/static/models/iJO1366.mat",
        df("iJO1366.mat"),
        "b5cfe21b6369a00e45d600b783f89521f5cc953e25ee52c5f1d0a3f83743be30",
    ),
    "iJO1366.xml" => download_data_file(
        "http://bigg.ucsd.edu/static/models/iJO1366.xml",
        df("iJO1366.xml"),
        "d6d9ec61ef6f155db5bb2f49549119dc13b96f6098b403ef82ea4240b27232eb",
    ),
    "ecoli_core_model.xml" => download_data_file(
        "http://systemsbiology.ucsd.edu/sites/systemsbiology.ucsd.edu/files/Attachments/Images/downloads/Ecoli_core/ecoli_core_model.xml",
        df("ecoli_core_model.xml"),
        "78692f8509fb36534f4f9b6ade23b23552044f3ecd8b48d84d484636922ae907",
    ),
    "e_coli_core.json" => download_data_file(
        "http://bigg.ucsd.edu/static/models/e_coli_core.json",
        df("e_coli_core.json"),
        "7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
    ),
    "e_coli_core.xml" => download_data_file(
        "http://bigg.ucsd.edu/static/models/e_coli_core.xml",
        df("e_coli_core.xml"),
        "b4db506aeed0e434c1f5f1fdd35feda0dfe5d82badcfda0e9d1342335ab31116",
    ),
    "e_coli_core.mat" => download_data_file(
        "http://bigg.ucsd.edu/static/models/e_coli_core.mat",
        df("e_coli_core.mat"),
        "478e6fa047ede1b248975d7565208ac9363a44dd64aad1900b63127394f4175b",
    ),
    "iJR904.mat" => download_data_file(
        "http://bigg.ucsd.edu/static/models/iJR904.mat",
        df("iJR904.mat"),
        "d17be86293d4caafc32b829da4e2d0d76eb45e1bb837e0138327043a83e20c6e",
    ),
)


text:
module TestMacroMLJ

using MLJBase, Test, Distances

# No type, no default
@mlj_model mutable struct A1
    a
end
a = A1()
@test ismissing(a.a)
a.a = 5
@test a.a == 5

# No type, with default
@mlj_model mutable struct A1b
    a = 5
end
a = A1b()
@test a.a == 5
a.a = "hello"
@test a.a == "hello"

# If a type is given but no default value is given, then the macro tries to fill
# a default value; either 0 if it's a Number type, or an empty string and otherwise fails.
@mlj_model mutable struct A1c
    a::Int
end
a = A1c()
@test a.a == 0
a = A1c(a=7)
@test a.a == 7
@test_throws InexactError A1c(a=5.3)
@test_throws MethodError A1c(a="hello")

# Type is given and default is given
@mlj_model mutable struct A1d
    a::Int = 5
end
a = A1d()
@test a.a == 5
a = A1d(a=7)
@test a.a == 7

# No type is given but a default and constraint
@mlj_model mutable struct A1e
    a = 5::(_ > 0)
end
a = A1e()
@test a.a == 5
a = A1e(a=7)
@test a.a == 7
@test @test_logs (:warn, "Constraint `model.a > 0` failed; using default: a=5.") A1e(a=-1).a==5
a = A1e(a=7.5)
@test a.a == 7.5

# Type is given with default and constraint
@mlj_model mutable struct A1f
    a::Int = 5::(_ > 0)
end
a = A1f()
@test a.a == 5
a = A1f(a=7)
@test a.a == 7
@test_throws InexactError A1f(a=7.5)
@test @test_logs (:warn, "Constraint `model.a > 0` failed; using default: a=5.") A1f(a=-1).a==5

abstract type FooBar end
@mlj_model mutable struct B1a <: FooBar
    a::Symbol = :auto::(_ in (:auto, :semi))
end
b = B1a()
@test b.a == :auto
b = B1a(a=:semi)
@test b.a == :semi
@test @test_logs (:warn, "Constraint `model.a in (:auto, :semi)` failed; using default: a=:auto.") B1a(a=:autos).a == :auto
@test_throws MethodError B1a(b="blah")

# == dependence on other types
@mlj_model mutable struct B1b
    a::SemiMetric = Euclidean()::(_ isa Metric)
end
@test B1b().a isa Euclidean
@test @test_logs (:warn, "Constraint `model.a isa Metric` failed; using default: a=Euclidean().") B1b(a=BhattacharyyaDist()).a isa Euclidean

@mlj_model mutable struct B1c
    a::SemiMetric = Euclidean()
end
@test B1c().a isa Euclidean

# Implicit defaults
@mlj_model mutable struct Ca
    a::String
end
@test Ca().a == ""
@mlj_model mutable struct Cb
    a::Any
end
@test Cb().a === missing
@mlj_model mutable struct Cc
    a::Union{Nothing,Int}
end
@test Cc().a === nothing
@mlj_model mutable struct Cd
    a::Union{Missing,Int}
end
@test Cd().a === missing

end
true


text:
# Test basic methods
@testset "Basics" begin
    # initialize model and semi_infinite variable
    m = InfiniteModel()
    @infinite_parameter(m, a in [0, 1])
    @infinite_parameter(m, b[1:2] in [0, 1], independent = true)
    @infinite_parameter(m, c[1:2] in [0, 1])
    @variable(m, ivref, Infinite(a, b, c))
    eval_supps = Dict{Int, Float64}(1 => 0.5, 3 => 1, 4 => 0, 5 => 0)
    var = SemiInfiniteVariable(ivref, eval_supps, [2], [2])
    object = VariableData(var, "var")
    idx = SemiInfiniteVariableIndex(1)
    vref = SemiInfiniteVariableRef(m, idx)
    gvref = GeneralVariableRef(m, 1, SemiInfiniteVariableIndex)
    bad_vref = SemiInfiniteVariableRef(m, SemiInfiniteVariableIndex(-1))
    # JuMP.owner_model
    @testset "JuMP.owner_model" begin
        @test owner_model(vref) === m
        @test owner_model(gvref) === m
    end
    # JuMP.index
    @testset "JuMP.index" begin
        @test index(vref) == idx
        @test index(gvref) == idx
    end
    # dispatch_variable_ref
    @testset "dispatch_variable_ref" begin
        @test isequal(dispatch_variable_ref(m, idx), vref)
        @test isequal(dispatch_variable_ref(gvref), vref)
    end
    # _add_data_object
    @testset "_add_data_object" begin
        @test InfiniteOpt._add_data_object(m, object) == idx
    end
    # _data_dictionary
    @testset "_data_dictionary" begin
        @test InfiniteOpt._data_dictionary(m, SemiInfiniteVariable) === m.semi_infinite_vars
        @test InfiniteOpt._data_dictionary(vref) === m.semi_infinite_vars
        @test InfiniteOpt._data_dictionary(gvref) === m.semi_infinite_vars
    end
    # JuMP.is_valid
    @testset "JuMP.is_valid" begin
        @test is_valid(m, vref)
        @test is_valid(m, gvref)
    end
    # _data_object
    @testset "_data_object" begin
        @test InfiniteOpt._data_object(vref) === object
        @test InfiniteOpt._data_object(gvref) === object
        @test_throws ErrorException InfiniteOpt._data_object(bad_vref)
    end
    # _core_variable_object
    @testset "_core_variable_object" begin
        @test InfiniteOpt._core_variable_object(vref) === var
        @test InfiniteOpt._core_variable_object(gvref) === var
    end
    # _set_core_variable_object
    @testset "_set_core_variable_object" begin
        @test InfiniteOpt._set_core_variable_object(vref, var) isa Nothing
    end
    # _measure_dependencies
    @testset "_measure_dependencies" begin
        @test InfiniteOpt._measure_dependencies(vref) == MeasureIndex[]
        @test InfiniteOpt._measure_dependencies(gvref) == MeasureIndex[]
    end
    # _constraint_dependencies
    @testset "_constraint_dependencies" begin
        @test InfiniteOpt._constraint_dependencies(vref) == InfOptConstraintIndex[]
        @test InfiniteOpt._constraint_dependencies(gvref) == InfOptConstraintIndex[]
    end
    # _derivative_dependencies
    @testset "_derivative_dependencies" begin
        @test InfiniteOpt._derivative_dependencies(vref) == DerivativeIndex[]
        @test InfiniteOpt._derivative_dependencies(gvref) == DerivativeIndex[]
    end
    # _object_numbers
    @testset "_object_numbers" begin
        @test InfiniteOpt._object_numbers(vref) == [2]
    end
    # _parameter_numbers
    @testset "_parameter_numbers" begin
        @test InfiniteOpt._parameter_numbers(vref) == [2]
    end
    # infinite_variable_ref
    @testset "infinite_variable_ref" begin
        @test isequal(infinite_variable_ref(vref), ivref)
        @test isequal(infinite_variable_ref(gvref), ivref)
    end
    # eval_supports
    @testset "eval_supports" begin
        @test eval_supports(vref) === eval_supps
        @test eval_supports(gvref) === eval_supps
    end
    # raw_parameter_refs
    @testset "raw_parameter_refs" begin
        @test isequal(raw_parameter_refs(vref), IC.VectorTuple(b[1]))
        @test isequal(raw_parameter_refs(gvref), IC.VectorTuple(b[1]))
    end
    # parameter_refs
    @testset "parameter_refs" begin
        @test isequal(parameter_refs(vref), (b[1],))
        @test isequal(parameter_refs(gvref), (b[1],))
    end
    # parameter_list
    @testset "parameter_list" begin
        @test isequal(parameter_list(vref), [b[1]])
        @test isequal(parameter_list(gvref), [b[1]])
    end
    # JuMP.set_name
    @testset "JuMP.set_name" begin
        # test normal
        @test isa(set_name(vref, "new"), Nothing)
        @test name(vref) == "new"
        @test_throws ErrorException set_name(bad_vref, "")
        # test default
        @test isa(set_name(gvref, "a"), Nothing)
        @test name(vref) == "a"
    end
    # JuMP.name
    @testset "JuMP.name" begin
        @test name(vref) == "a"
        @test name(bad_vref) == ""
    end
    # _make_variable_ref
    @testset "_make_variable_ref" begin
        @test isequal(InfiniteOpt._make_variable_ref(m, idx), gvref)
    end
    # _var_name_dict
    @testset "_var_name_dict" begin
        @test InfiniteOpt._var_name_dict(m) isa Nothing
    end
    # _update_var_name_dict
    @testset "_update_var_name_dict" begin
        m.name_to_var = Dict{String, ObjectIndex}()
        dict = InfiniteOpt._data_dictionary(vref)
        @test InfiniteOpt._update_var_name_dict(m, dict) isa Nothing
        @test InfiniteOpt._var_name_dict(m) == Dict(name(vref) => idx)
        m.name_to_var = nothing
    end
    # parameter_by_name
    @testset "JuMP.variable_by_name" begin
        # test normal
        @test isequal(variable_by_name(m, "a"), gvref)
        @test variable_by_name(m, "test") isa Nothing
        # prepare variable with same name
        idx2 = SemiInfiniteVariableIndex(2)
        @test InfiniteOpt._add_data_object(m, object) == idx2
        vref2 = SemiInfiniteVariableRef(m, idx2)
        @test set_name(vref2, "a") isa Nothing
        # test multiple name error
        @test_throws ErrorException variable_by_name(m, "a")
    end
    # _delete_data_object
    @testset "_delete_data_object" begin
        @test InfiniteOpt._delete_data_object(vref) isa Nothing
        @test length(InfiniteOpt._data_dictionary(vref)) == 1
        @test !is_valid(m, vref)
    end
end

# Test definition methods
@testset "Definition" begin
    # initialize model and semi_infinite variable
    m = InfiniteModel()
    @infinite_parameter(m, a in [0, 1])
    @infinite_parameter(m, b[1:2] in [0, 1], independent = true)
    @infinite_parameter(m, c[1:2] in [0, 1])
    @variable(m, ivref, Infinite(a, b, c))
    eval_supps = Dict{Int, Float64}(1 => 0.5, 3 => 1, 4 => 0, 5 => 0)
    info = VariableInfo(false, NaN, false, NaN, false, NaN, false, NaN, false, false)
    info2 = VariableInfo(false, NaN, false, NaN, false, NaN, false, NaN, false, true)
    # test _process_value
    @testset "_process_value" begin
        @test isequal(InfiniteOpt._process_value(0), 0)
        @test isequal(InfiniteOpt._process_value(0a + 1), 1)
        @test isequal(InfiniteOpt._process_value(1a), a)
        @test isequal(InfiniteOpt._process_value(a + 2), a + 2)
    end
    # test SemiInfinite{V, T}
    @testset "SemiInfinite{T, V}" begin
        @test isequal(SemiInfinite(ivref, 0, b, c).infinite_variable_ref, ivref)
        @test SemiInfinite(ivref, 0, [b[1], 0], c).parameter_values isa IC.VectorTuple
    end
    # test JuMP.build_variable for macros
    @testset "JuMP.build_variable macro wrapper" begin
        # test errors
        @test_throws ErrorException build_variable(error, info,
                                                   SemiInfinite(ivref, 0, b, a),
                                                   bad = 42)
        @test_throws ErrorException build_variable(error, info2,
                                                   SemiInfinite(ivref, 0, b, a))
        @test_throws ErrorException build_variable(error, info, SemiInfinite(0))
        @test_throws ErrorException build_variable(error, info,
                                                   SemiInfinite(ivref, 0, b,
                                                   dispatch_variable_ref.(c)))
        @test_throws ErrorException build_variable(error, info,
                                                   SemiInfinite(ivref, 0, b))
        # test normal
        @test isequal(build_variable(error, info,
                             SemiInfinite(ivref, 0, b, c)).infinite_variable_ref, ivref)
        @test build_variable(error, info,
                             SemiInfinite(ivref, 0, b, c)).eval_supports == Dict(1 => 0.0)
    end
    # test JuMP.build_variable
    @testset "JuMP.build_variable" begin
        # test errors
        dvref = InfiniteOpt._make_variable_ref(m, FiniteVariableIndex(1))
        @test_throws ErrorException build_variable(error, dvref, eval_supps)
        eval_supps[6] = 1
        @test_throws ErrorException build_variable(error, ivref, eval_supps)
        delete!(eval_supps, 6)
        eval_supps[4] = 2
        @test_throws ErrorException build_variable(error, ivref, eval_supps)
        eval_supps[4] = 0
        # test normal
        @test isequal(build_variable(error, ivref, eval_supps).infinite_variable_ref, ivref)
        @test build_variable(error, ivref, eval_supps).eval_supports === eval_supps
        @test build_variable(error, ivref, eval_supps).object_nums == [2]
        @test build_variable(error, ivref, eval_supps, check = false).object_nums == [2]
        @test build_variable(error, ivref, eval_supps, check = false).parameter_nums == [2]
    end
    # test JuMP.add_variable
    @testset "JuMP.add_variable" begin
        # test errors
        m2 = InfiniteModel()
        var = build_variable(error, ivref, eval_supps, check = false)
        @test_throws VariableNotOwned{InfiniteVariableRef} add_variable(m2, var)
        # test normal
        idx = SemiInfiniteVariableIndex(1)
        vref = SemiInfiniteVariableRef(m, idx)
        gvref = GeneralVariableRef(m, 1, SemiInfiniteVariableIndex)
        @test isequal(add_variable(m, var), gvref)
        @test name(vref) == ""
        @test eval_supports(vref) === eval_supps
        @test supports(a) == [0.5]
        @test supports(b[2]) == [1]
        @test supports(c) == zeros(2, 1)
        @test supports(b[1]) == []
        @test InfiniteOpt._object_numbers(vref) == [2]
        @test InfiniteOpt._semi_infinite_variable_dependencies(ivref) == [idx]
        # test with set name (redundant add)
        @test isequal(add_variable(m, var, "cat"), gvref)
        @test name(vref) == "cat"
        # test with partial supports
        eval_supps2 = Dict{Int, Float64}(1 => 0.5, 3 => 0, 5 => 1)
        var = build_variable(error, ivref, eval_supps2, check = false)
        idx = SemiInfiniteVariableIndex(2)
        vref = SemiInfiniteVariableRef(m, idx)
        gvref = GeneralVariableRef(m, 2, SemiInfiniteVariableIndex)
        @test isequal(add_variable(m, var), gvref)
        @test InfiniteOpt._data_object(vref).name == ""
        @test supports(a) == [0.5]
        @test supports(b[2]) == [0, 1]
        @test supports(b[1]) == []
        @test supports(c) == zeros(2, 1)
    end
    # test macro definition
    @testset "Macro Definition" begin
        # anonymous definition
        vref = GeneralVariableRef(m, 3, SemiInfiniteVariableIndex)
        @test isequal(@variable(m, variable_type = SemiInfinite(ivref, 0, [b[1], 0], c)), vref)
        @test isequal(parameter_refs(vref), (b[1], c))
        @test eval_supports(vref) == Dict(1 => 0.0, 3 => 0.0)
        @test supports(a) == [0, 0.5]
        @test supports(b[2]) == [0, 1]
        @test supports(b[1]) == []
        # explicit definition (redundant)
        vref = GeneralVariableRef(m, 3, SemiInfiniteVariableIndex)
        @test isequal(@variable(m, test, SemiInfinite(ivref, 0, [b[1], 0], c)), vref)
        @test isequal(parameter_refs(vref), (b[1], c))
        @test eval_supports(vref) == Dict(1 => 0.0, 3 => 0.0)
        @test name(vref) == "test"
        # array definition
        vrefs = [GeneralVariableRef(m, i, SemiInfiniteVariableIndex) for i in 4:5]
        @test isequal(@variable(m, [i = 1:2], SemiInfinite(ivref, i - 1, b, c)), vrefs)
        @test isequal(parameter_refs(vrefs[1]), (b, c))
        @test eval_supports.(vrefs) == [Dict(1 => 0.0), Dict(1 => 1.0)]
    end
    # test restriciton definition
    @testset "Restriction" begin
        # test errors
        @test_throws ErrorException restrict(ivref, 0, b)
        @test_throws ErrorException ivref(0, b)
        # test normal wth restrict
        vref = GeneralVariableRef(m, 6, SemiInfiniteVariableIndex)
        @test isequal(restrict(ivref, 0.5, [b[1], 0], c), vref)
        @test isequal(parameter_refs(vref), (b[1], c))
        @test eval_supports(vref) == Dict(1 => 0.5, 3 => 0.0)
        # test normal functionally
        vref = GeneralVariableRef(m, 6, SemiInfiniteVariableIndex)
        @test isequal(ivref(0.5, [b[1], 0], c), vref)
        @test isequal(parameter_refs(vref), (b[1], c))
        @test eval_supports(vref) == Dict(1 => 0.5, 3 => 0.0)
    end
end

# Test the info methods
@testset "Info" begin
    # initialize model and semi_infinite variable
    m = InfiniteModel()
    @infinite_parameter(m, a in [0, 1])
    @infinite_parameter(m, b[1:2] in [0, 1], independent = true)
    @infinite_parameter(m, c[1:2] in [0, 1])
    @variable(m, 0 <= ivref1 <= 1, Infinite(a, b, c), Int)
    @variable(m, ivref2 == 1, Infinite(a, b, c), Bin, start = 0)
    eval_supps = Dict{Int, Float64}(1 => 0.5, 3 => 1, 4 => 0, 5 => 0)
    var1 = build_variable(error, ivref1, eval_supps, check = false)
    gvref1 = add_variable(m, var1)
    rvref1 = dispatch_variable_ref(gvref1)
    var2 = build_variable(error, ivref2, eval_supps, check = false)
    gvref2 = add_variable(m, var2)
    rvref2 = dispatch_variable_ref(gvref2)
    dvref1 = dispatch_variable_ref(ivref1)
    dvref2 = dispatch_variable_ref(ivref2)
    # has_lower_bound
    @testset "JuMP.has_lower_bound" begin
        @test has_lower_bound(rvref1)
        @test has_lower_bound(gvref1)
        @test !has_lower_bound(rvref2)
    end
    # lower_bound
    @testset "JuMP.lower_bound" begin
        @test lower_bound(rvref1) == 0
        @test lower_bound(gvref1) == 0
        @test_throws ErrorException lower_bound(rvref2)
    end
    # _lower_bound_index
    @testset "InfiniteOpt._lower_bound_index" begin
        @test InfiniteOpt._lower_bound_index(rvref1) == InfiniteOpt._lower_bound_index(dvref1)
        @test_throws ErrorException InfiniteOpt._lower_bound_index(rvref2)
    end
    # LowerBoundRef
    @testset "JuMP.LowerBoundRef" begin
        @test LowerBoundRef(rvref1) == LowerBoundRef(dvref1)
        @test_throws ErrorException LowerBoundRef(rvref2)
        @test_throws ErrorException LowerBoundRef(gvref2)
    end
    # has_upper_bound
    @testset "JuMP.has_upper_bound" begin
        @test has_upper_bound(rvref1)
        @test has_upper_bound(gvref1)
        @test !has_upper_bound(rvref2)
    end
    # upper_bound
    @testset "JuMP.upper_bound" begin
        @test upper_bound(rvref1) == 1
        @test upper_bound(gvref1) == 1
        @test_throws ErrorException upper_bound(rvref2)
    end
    # _upper_bound_index
    @testset "InfiniteOpt._upper_bound_index" begin
        @test InfiniteOpt._upper_bound_index(rvref1) == InfiniteOpt._upper_bound_index(dvref1)
        @test_throws ErrorException InfiniteOpt._upper_bound_index(rvref2)
    end
    # UpperBoundRef
    @testset "JuMP.UpperBoundRef" begin
        @test UpperBoundRef(rvref1) == UpperBoundRef(dvref1)
        @test_throws ErrorException UpperBoundRef(rvref2)
        @test_throws ErrorException UpperBoundRef(gvref2)
    end
    # is_fixed
    @testset "JuMP.is_fixed" begin
        @test is_fixed(rvref2)
        @test !is_fixed(rvref1)
        @test !is_fixed(gvref1)
    end
    # fix_value
    @testset "JuMP.fix_value" begin
        @test fix_value(rvref2) == 1
        @test_throws ErrorException fix_value(rvref1)
        @test_throws ErrorException fix_value(gvref1)
    end
    # _fix_index
    @testset "InfiniteOpt._fix_index" begin
        @test InfiniteOpt._fix_index(rvref2) == InfiniteOpt._fix_index(dvref2)
        @test_throws ErrorException InfiniteOpt._fix_index(rvref1)
    end
    # FixRef
    @testset "JuMP.FixRef" begin
        @test FixRef(rvref2) == FixRef(dvref2)
        @test FixRef(rvref2) == FixRef(dvref2)
        @test_throws ErrorException FixRef(rvref1)
    end
    # start_value
    @testset "JuMP.start_value" begin
        @test_throws ErrorException start_value(rvref2)
        @test_throws ErrorException start_value(gvref2)
    end
    # start_value_function
    @testset "start_value_function" begin
        @test start_value_function(rvref2)([0]) == 0
        @test start_value_function(gvref2)([0]) == 0
    end
    # is_binary
    @testset "JuMP.is_binary" begin
        @test is_binary(rvref2)
        @test is_binary(gvref2)
        @test !is_binary(rvref1)
    end
    # _binary_index
    @testset "InfiniteOpt._binary_index" begin
        @test InfiniteOpt._binary_index(rvref2) == InfiniteOpt._binary_index(dvref2)
        @test_throws ErrorException InfiniteOpt._binary_index(rvref1)
    end
    # BinaryRef
    @testset "JuMP.BinaryRef" begin
        @test BinaryRef(rvref2) == BinaryRef(dvref2)
        @test BinaryRef(gvref2) == BinaryRef(dvref2)
        @test_throws ErrorException BinaryRef(rvref1)
    end
    # is_integer
    @testset "JuMP.is_integer" begin
        @test is_integer(rvref1)
        @test !is_integer(rvref2)
        @test !is_integer(gvref2)
    end
    # _integer_index
    @testset "InfiniteOpt._integer_index" begin
        @test InfiniteOpt._integer_index(rvref1) == InfiniteOpt._integer_index(dvref1)
        @test_throws ErrorException InfiniteOpt._integer_index(rvref2)
    end
    # IntegerRef
    @testset "JuMP.IntegerRef" begin
        @test IntegerRef(rvref1) == IntegerRef(dvref1)
        @test IntegerRef(gvref1) == IntegerRef(dvref1)
        @test_throws ErrorException IntegerRef(rvref2)
    end
end

# test usage methods
@testset "Usage" begin
    # initialize model and stuff
    m = InfiniteModel()
    @infinite_parameter(m, t in [0, 1])
    @infinite_parameter(m, x[1:2] in [-1, 1])
    @variable(m, y, Infinite(t, x))
    eval_supps = Dict{Int, Float64}(1 => 0.5, 3 => 1)
    var = build_variable(error, y, eval_supps, check = false)
    gvref = add_variable(m, var)
    vref = dispatch_variable_ref(gvref)
    # test used_by_measure
    @testset "used_by_measure" begin
        @test !used_by_measure(vref)
        push!(InfiniteOpt._measure_dependencies(vref), MeasureIndex(1))
        @test used_by_measure(gvref)
        @test used_by_measure(vref)
        empty!(InfiniteOpt._measure_dependencies(vref))
    end
    # test used_by_constraint
    @testset "used_by_constraint" begin
        @test !used_by_constraint(vref)
        push!(InfiniteOpt._constraint_dependencies(vref), InfOptConstraintIndex(1))
        @test used_by_constraint(gvref)
        @test used_by_constraint(vref)
        empty!(InfiniteOpt._constraint_dependencies(vref))
    end
    # test used_by_derivative
    @testset "used_by_derivative" begin
        @test !used_by_derivative(vref)
        push!(InfiniteOpt._derivative_dependencies(vref), DerivativeIndex(1))
        @test used_by_derivative(gvref)
        @test used_by_derivative(vref)
        empty!(InfiniteOpt._derivative_dependencies(vref))
    end
    # test used_by_objective
    @testset "used_by_objective" begin
        @test !used_by_objective(gvref)
        @test !used_by_objective(vref)
    end
    # test is_used
    @testset "is_used" begin
        # test not used
        @test !is_used(gvref)
        # test used by constraint and/or measure
        push!(InfiniteOpt._constraint_dependencies(vref), InfOptConstraintIndex(1))
        @test is_used(y)
        empty!(InfiniteOpt._constraint_dependencies(vref))
        # test used by derivative
        func = (x) -> NaN
        num = 0.
        info = VariableInfo{Float64, Float64, Float64, Function}(true, num, true,
                                         num, true, num, false, func, true, true)
        deriv = Derivative(info, true, y, t)
        object = VariableData(deriv)
        idx = DerivativeIndex(1)
        dref = DerivativeRef(m, idx)
        @test InfiniteOpt._add_data_object(m, object) == idx
        push!(InfiniteOpt._derivative_dependencies(vref), idx)
        @test !is_used(vref)
        push!(InfiniteOpt._constraint_dependencies(dref), InfOptConstraintIndex(2))
        @test is_used(vref)
        empty!(InfiniteOpt._derivative_dependencies(vref))
    end
end

# Test variable queries for models
@testset "Model Queries" begin
    # initialize model, parameter, and variables
    m = InfiniteModel()
    @infinite_parameter(m, a in [0, 1])
    @infinite_parameter(m, b[1:2] in [0, 1], independent = true)
    @infinite_parameter(m, c[1:2] in [0, 1])
    @variable(m, ivref, Infinite(a, b, c))
    @variable(m, pvref, Point(ivref, 0, [0, 0], [0, 0]))
    @variable(m, hvref)
    eval_supps = Dict{Int, Float64}(1 => 0.5, 3 => 1, 4 => 0, 5 => 0)
    var = build_variable(error, ivref, eval_supps, check = false)
    rvref = add_variable(m, var)
    # num_variables
    @testset "JuMP.num_variables" begin
        @test num_variables(m, InfiniteVariable) == 1
        @test num_variables(m, PointVariable) == 1
        @test num_variables(m, FiniteVariable) == 1
        @test num_variables(m, SemiInfiniteVariable) == 1
        @test num_variables(m) == 4
    end
    # all_variables
    @testset "JuMP.all_variables" begin
        @test all_variables(m, InfiniteVariable) == [ivref]
        @test all_variables(m, PointVariable) == [pvref]
        @test all_variables(m, FiniteVariable) == [hvref]
        @test all_variables(m, SemiInfiniteVariable) == [rvref]
        @test all_variables(m) == [ivref, rvref, pvref, hvref]
    end
end


text:
module BilinearModels

include("fit.jl"); using .Fit
include("infer.jl"); using .Infer
include("common.jl"); using .Common
include("outcome_nb.jl"); using .Outcome_NB
include("sim.jl"); using .Sim

end # module BilinearModels


text:
using HypothesisTests, Test

@testset "Permutation" begin
@test isapprox(pvalue(ExactPermutationTest([1,2,3],[4,5,6],mean), tail=:both), 0.1)
@test isapprox(pvalue(ExactPermutationTest([4,5,6],[1,2,3],mean), tail=:both), 0.1)
@test isapprox(pvalue(ExactPermutationTest([1,2,3],[4,5,6],mean), tail=:left), 0.05)
@test isapprox(pvalue(ExactPermutationTest([4,5,6],[1,2,3],mean), tail=:right), 0.05)

Random.seed!(12345)
@test isapprox(pvalue(ApproximatePermutationTest([1,2,3],[4,5,6],mean,100), tail=:both), 0.09)
@test isapprox(pvalue(ApproximatePermutationTest([4,5,6],[1,2,3],mean,100), tail=:both), 0.04)
@test isapprox(pvalue(ApproximatePermutationTest([1,2,3],[4,5,6],mean,100), tail=:left), 0.08)
@test isapprox(pvalue(ApproximatePermutationTest([4,5,6],[1,2,3],mean,100), tail=:right), 0.07)
end


text:
using Documenter, ApproxInferenceBase

makedocs(
    modules = [ApproxInferenceBase],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Julia Approx Inference members",
    sitename = "ApproxInferenceBase.jl",
    pages = Any["index.md"]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)
deploydocs(
    repo = "github.com/JuliaApproxInference/ApproxInferenceBase.jl.git",
)