text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestConstraintGeomean

using Test

using MathOptInterface
const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_GeoMean_basic()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged_mock = MOI.Bridges.Constraint.GeoMean{Float64}(mock)
    MOI.Test.runtests(
        bridged_mock,
        MOI.Test.Config(),
        include = [
            "test_basic_VectorOfVariables_GeometricMeanCone",
            "test_basic_VectorAffineFunction_GeometricMeanCone",
            "test_basic_VectorQuadraticFunction_GeometricMeanCone",
        ],
    )
    return
end

function test_conic_GeometricMeanCone_VectorAffineFunction()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.GeoMean{Float64}(mock)
    var_primal = vcat(ones(4), 2, √2, √2)
    rsoc_duals = [[√2, √2, -2] / 3, [1, 1, -√2] / 3, [1, 1, -√2] / 3]
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            var_primal,
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1, -inv(3)],
            (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone) => rsoc_duals,
        )
    MOI.Test.test_conic_GeometricMeanCone_VectorOfVariables(
        bridged_mock,
        config,
    )
    MOI.empty!(bridged_mock)
    MOI.Test.test_conic_GeometricMeanCone_VectorAffineFunction(
        bridged_mock,
        config,
    )
    var_names = ["t", "x", "y", "z", "x21", "x11", "x12"]
    MOI.set(
        mock,
        MOI.VariableName(),
        MOI.get(mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    lessthan = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64},
            MOI.LessThan{Float64},
        }(),
    )
    @test length(lessthan) == 2
    MOI.set.(mock, MOI.ConstraintName(), lessthan, ["lessthan1", "lessthan2"])
    rsoc = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.RotatedSecondOrderCone,
        }(),
    )
    @test length(rsoc) == 3
    MOI.set.(mock, MOI.ConstraintName(), rsoc, ["rsoc21", "rsoc11", "rsoc12"])
    s = """
    variables: t, x, y, z, x21, x11, x12
    lessthan1: t + -0.5 * x21 in MathOptInterface.LessThan(0.0)
    lessthan2: x + y + z in MathOptInterface.LessThan(3.0)
    rsoc11: [1.0x, y, x11] in MathOptInterface.RotatedSecondOrderCone(3)
    rsoc12: [z, 0.5 * x21, x12] in MathOptInterface.RotatedSecondOrderCone(3)
    rsoc21: [1.0x11, x12, x21] in MathOptInterface.RotatedSecondOrderCone(3)
    maxobjective: t
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        mock,
        model,
        var_names,
        ["lessthan1", "lessthan2", "rsoc11", "rsoc12", "rsoc21"],
    )
    var_names = ["t", "x", "y", "z"]
    MOI.set(
        bridged_mock,
        MOI.VariableName(),
        MOI.get(bridged_mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    lessthan = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64},
            MOI.LessThan{Float64},
        }(),
    )
    @test length(lessthan) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), lessthan[1], "lessthan")
    geomean = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.GeometricMeanCone,
        }(),
    )
    @test length(geomean) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), geomean[1], "geomean")
    s = """
    variables: t, x, y, z
    lessthan: x + y + z in MathOptInterface.LessThan(3.0)
    geomean: [1.0t, x, y, z] in MathOptInterface.GeometricMeanCone(4)
    maxobjective: t
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        bridged_mock,
        model,
        var_names,
        ["lessthan", "geomean"],
    )
    # set primal/dual start is not yet implemented for GeoMean bridge
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.GeometricMeanCone,
            }(),
        ),
    )
    _test_delete_bridge(
        bridged_mock,
        ci,
        4,
        (
            (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 1),
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 0),
        ),
    )
    return
end

function test_conic_GeometricMeanCone_VectorAffineFunction_2()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.GeoMean{Float64}(mock)
    var_primal = vcat(ones(10), 4.0, fill(2 * √2, 2), fill(2, 4), fill(√2, 8))
    rsoc_duals = vcat(
        [[1 / √2, 1 / √2, -1] * 4 / 9],
        fill([1, 1, -√2] * 2 / 9, 2),
        fill([1 / √2, 1 / √2, -1] * 2 / 9, 4),
        fill([1, 1, -√2] * 1 / 9, 8),
    )
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            var_primal,
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1],
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                fill(-inv(9), 9),
            (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone) => rsoc_duals,
        )
    MOI.Test.test_conic_GeometricMeanCone_VectorOfVariables_2(
        bridged_mock,
        config,
    )
    MOI.empty!(bridged_mock)
    MOI.Test.test_conic_GeometricMeanCone_VectorAffineFunction_2(
        bridged_mock,
        config,
    )
    x_names = ["x$(i)" for i in 1:9]
    l1_names = ["x1$(i)" for i in 1:8]
    l2_names = ["x2$(i)" for i in 1:4]
    l3_names = ["x3$(i)" for i in 1:2]
    var_names = vcat("t", x_names, "x41", l3_names, l2_names, l1_names)
    MOI.set(
        mock,
        MOI.VariableName(),
        MOI.get(mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    lessthan = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64},
            MOI.LessThan{Float64},
        }(),
    )
    @test length(lessthan) == 1
    MOI.set(mock, MOI.ConstraintName(), lessthan[1], "lessthan")
    equalto = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64},
            MOI.EqualTo{Float64},
        }(),
    )
    @test length(equalto) == 9
    equalto_names = ["equalto$(i)" for i in 1:9]
    MOI.set.(mock, MOI.ConstraintName(), equalto, equalto_names)
    rsoc = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.RotatedSecondOrderCone,
        }(),
    )
    @test length(rsoc) == 15
    rsoc_names = vcat(
        "rsoc41",
        ["rsoc3$(i)" for i in 1:2],
        ["rsoc2$(i)" for i in 1:4],
        ["rsoc1$(i)" for i in 1:8],
    )
    MOI.set.(mock, MOI.ConstraintName(), rsoc, rsoc_names)
    varnames_str = vcat("t", [", " * v for v in var_names[2:end]])
    s = """
    variables: $(varnames_str...)
    lessthan: t + -0.25 * x41 in MathOptInterface.LessThan(0.0)
    equalto1: 1.0x1 in MathOptInterface.EqualTo(1.0)
    equalto2: 1.0x2 in MathOptInterface.EqualTo(1.0)
    equalto3: 1.0x3 in MathOptInterface.EqualTo(1.0)
    equalto4: 1.0x4 in MathOptInterface.EqualTo(1.0)
    equalto5: 1.0x5 in MathOptInterface.EqualTo(1.0)
    equalto6: 1.0x6 in MathOptInterface.EqualTo(1.0)
    equalto7: 1.0x7 in MathOptInterface.EqualTo(1.0)
    equalto8: 1.0x8 in MathOptInterface.EqualTo(1.0)
    equalto9: 1.0x9 in MathOptInterface.EqualTo(1.0)
    rsoc41: [1.0x31, x32, x41] in MathOptInterface.RotatedSecondOrderCone(3)
    rsoc31: [1.0x21, x22, x31] in MathOptInterface.RotatedSecondOrderCone(3)
    rsoc32: [1.0x23, x24, x32] in MathOptInterface.RotatedSecondOrderCone(3)
    rsoc21: [1.0x11, x12, x21] in MathOptInterface.RotatedSecondOrderCone(3)
    rsoc22: [1.0x13, x14, x22] in MathOptInterface.RotatedSecondOrderCone(3)
    rsoc23: [1.0x15, x16, x23] in MathOptInterface.RotatedSecondOrderCone(3)
    rsoc24: [1.0x17, x18, x24] in MathOptInterface.RotatedSecondOrderCone(3)
    rsoc11: [1.0x1, x2, x11] in MathOptInterface.RotatedSecondOrderCone(3)
    rsoc12: [1.0x3, x4, x12] in MathOptInterface.RotatedSecondOrderCone(3)
    rsoc13: [1.0x5, x6, x13] in MathOptInterface.RotatedSecondOrderCone(3)
    rsoc14: [1.0x7, x8, x14] in MathOptInterface.RotatedSecondOrderCone(3)
    rsoc15: [1.0x9, 0.25 * x41, x15] in MathOptInterface.RotatedSecondOrderCone(3)
    rsoc16: [0.25 * x41, 0.25 * x41, x16] in MathOptInterface.RotatedSecondOrderCone(3)
    rsoc17: [0.25 * x41, 0.25 * x41, x17] in MathOptInterface.RotatedSecondOrderCone(3)
    rsoc18: [0.25 * x41, 0.25 * x41, x18] in MathOptInterface.RotatedSecondOrderCone(3)
    maxobjective: t
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        mock,
        model,
        var_names,
        vcat("lessthan", equalto_names, rsoc_names),
    )
    var_names = vcat("t", ["x$(i)" for i in 1:9])
    MOI.set(
        bridged_mock,
        MOI.VariableName(),
        MOI.get(bridged_mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    equalto = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64},
            MOI.EqualTo{Float64},
        }(),
    )
    @test length(equalto) == 9
    equalto_names = ["equalto$(i)" for i in 1:9]
    MOI.set.(bridged_mock, MOI.ConstraintName(), equalto, equalto_names)
    geomean = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.GeometricMeanCone,
        }(),
    )
    @test length(geomean) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), geomean[1], "geomean")

    s = """
    variables: t, x1, x2, x3, x4, x5, x6, x7, x8, x9
    equalto1: 1.0x1 in MathOptInterface.EqualTo(1.0)
    equalto2: 1.0x2 in MathOptInterface.EqualTo(1.0)
    equalto3: 1.0x3 in MathOptInterface.EqualTo(1.0)
    equalto4: 1.0x4 in MathOptInterface.EqualTo(1.0)
    equalto5: 1.0x5 in MathOptInterface.EqualTo(1.0)
    equalto6: 1.0x6 in MathOptInterface.EqualTo(1.0)
    equalto7: 1.0x7 in MathOptInterface.EqualTo(1.0)
    equalto8: 1.0x8 in MathOptInterface.EqualTo(1.0)
    equalto9: 1.0x9 in MathOptInterface.EqualTo(1.0)
    geomean: [1.0t, x1, x2, x3, x4, x5, x6, x7, x8, x9] in MathOptInterface.GeometricMeanCone(10)
    maxobjective: t
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        bridged_mock,
        model,
        var_names,
        vcat(equalto_names, "geomean"),
    )
    # set primal/dual start is not yet implemented for GeoMean bridge
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.GeometricMeanCone,
            }(),
        ),
    )
    _test_delete_bridge(
        bridged_mock,
        ci,
        10,
        (
            (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 0),
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 0),
        ),
    )
    return
end

function test_conic_GeometricMeanCone_VectorAffineFunction_3()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.GeoMean{Float64}(mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [2, 2],
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-2, -2],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[0]],
        )
    MOI.Test.test_conic_GeometricMeanCone_VectorOfVariables_3(
        bridged_mock,
        config,
    )
    MOI.empty!(bridged_mock)
    MOI.Test.test_conic_GeometricMeanCone_VectorAffineFunction_3(
        bridged_mock,
        config,
    )
    MOI.set(
        mock,
        MOI.VariableName(),
        MOI.get(mock, MOI.ListOfVariableIndices()),
        ["t", "x"],
    )
    lessthan = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64},
            MOI.LessThan{Float64},
        }(),
    )
    @test length(lessthan) == 2
    MOI.set.(mock, MOI.ConstraintName(), lessthan, ["lessthan1", "lessthan2"])
    rsoc = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.RotatedSecondOrderCone,
        }(),
    )
    @test isempty(rsoc)
    nonneg = MOI.get(
        mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.Nonnegatives,
        }(),
    )
    @test length(nonneg) == 1
    MOI.set(mock, MOI.ConstraintName(), nonneg[1], "nonneg")

    s = """
    variables: t, x
    lessthan1: t + -1.0x in MathOptInterface.LessThan(0.0)
    lessthan2: 1.0x in MathOptInterface.LessThan(2.0)
    nonneg: [1.0x] in MathOptInterface.Nonnegatives(1)
    maxobjective: 2.0 * t
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        mock,
        model,
        ["t", "x"],
        ["lessthan1", "lessthan2", "nonneg"],
    )
    MOI.set(
        bridged_mock,
        MOI.VariableName(),
        MOI.get(bridged_mock, MOI.ListOfVariableIndices()),
        ["t", "x"],
    )
    lessthan = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64},
            MOI.LessThan{Float64},
        }(),
    )
    @test length(lessthan) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), lessthan[1], "lessthan")
    geomean = MOI.get(
        bridged_mock,
        MOI.ListOfConstraintIndices{
            MOI.VectorAffineFunction{Float64},
            MOI.GeometricMeanCone,
        }(),
    )
    @test length(geomean) == 1
    MOI.set(bridged_mock, MOI.ConstraintName(), geomean[1], "geomean")

    s = """
    variables: t, x
    lessthan: 1.0x in MathOptInterface.LessThan(2.0)
    geomean: [1.0t, x] in MathOptInterface.GeometricMeanCone(2)
    maxobjective: 2.0 * t
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        bridged_mock,
        model,
        ["t", "x"],
        ["lessthan", "geomean"],
    )
    # set primal/dual start is not yet implemented for GeoMean bridge
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.GeometricMeanCone,
            }(),
        ),
    )
    _test_delete_bridge(
        bridged_mock,
        ci,
        2,
        (
            (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 1),
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 0),
        ),
    )
    return
end

end  # module

TestConstraintGeomean.runtests()

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
#!/usr/bin/env julia

using Test, Printf, EAGO, MathOptInterface, Cbc #, MINLPTests

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

include(joinpath(@__DIR__, "moit_tests.jl"))
#include(joinpath(@__DIR__, "minlp_tests.jl"))

#include(joinpath(@__DIR__, "branch_bound.jl"))
#include(joinpath(@__DIR__, "domain_reduction.jl"))
#include(joinpath(@__DIR__, "semiinfinite_tests.jl"))

#, JuMP, Ipopt, ForwardDiff
#using IntervalArithmetic, SpecialFunctions
#using ForwardDiff: Dual, Partials
#include("script_optimizer.jl")

ext: jl


text:
using NormalizeQuantiles

using Base.Test

# write your own tests here
@test 1 == 1

if VERSION >= v"0.4.0-"

testfloat = [ 3.0 2.0 8.0 1.0 ; 4.0 5.0 6.0 2.0 ; 9.0 7.0 8.0 3.0 ; 5.0 2.0 8.0 4.0 ]
r=normalizeQuantiles(testfloat)
@test mean(r[:,1]) >= 4.8124 && mean(r[:,1]) <= 4.8126
@test mean(r[:,2]) >= 4.8124 && mean(r[:,2]) <= 4.8126
@test mean(r[:,3]) >= 4.8124 && mean(r[:,3]) <= 4.8126
@test mean(r[:,4]) >= 4.8124 && mean(r[:,4]) <= 4.8126

sa=SharedArray(Nullable{Float64},(size(testfloat,1),size(testfloat,2)));
sa[:]=testfloat[:]
r=normalizeQuantiles(sa)
r=convert(Array{Float64},reshape([get(r[i]) for i=1:length(r)],size(r)))
@test mean(r[:,1]) >= 4.8124 && mean(r[:,1]) <= 4.8126
@test mean(r[:,2]) >= 4.8124 && mean(r[:,2]) <= 4.8126
@test mean(r[:,3]) >= 4.8124 && mean(r[:,3]) <= 4.8126
@test mean(r[:,4]) >= 4.8124 && mean(r[:,4]) <= 4.8126

testfloat = [ 3.5 2.0 8.1 1.0 ; 4.5 5.0 6.0 2.0 ; 9.0 7.6 8.2 3.0 ; 5.0 2.0 8.0 4.0 ]
r=normalizeQuantiles(testfloat)
@test mean(r[:,1]) >= 4.93124 && mean(r[:,1]) <= 4.93125
@test mean(r[:,2]) >= 4.93124 && mean(r[:,2]) <= 4.93125
@test mean(r[:,3]) >= 4.93124 && mean(r[:,3]) <= 4.93125
@test mean(r[:,4]) >= 4.93124 && mean(r[:,4]) <= 4.93125

sa=SharedArray(Nullable{Float64},(size(testfloat,1),size(testfloat,2)));
sa[:]=testfloat[:]
r=normalizeQuantiles(sa)
r=convert(Array{Float64},reshape([get(r[i]) for i=1:length(r)],size(r)))
@test mean(r[:,1]) >= 4.93124 && mean(r[:,1]) <= 4.93125
@test mean(r[:,2]) >= 4.93124 && mean(r[:,2]) <= 4.93125
@test mean(r[:,3]) >= 4.93124 && mean(r[:,3]) <= 4.93125
@test mean(r[:,4]) >= 4.93124 && mean(r[:,4]) <= 4.93125

testfloat=[ 3.0 2.0 1.0 ; 4.0 5.0 6.0 ; 9.0 7.0 8.0 ; 5.0 2.0 8.0 ]
check=[ 2.0 3.0 2.0 ; 4.0 6.0 4.0 ; 8.0 8.0 7.0 ; 6.0 3.0 7.0 ]
qn=normalizeQuantiles(testfloat)
@test qn == check
sa=SharedArray(Nullable{Float64},(size(testfloat,1),size(testfloat,2)));
sa[:]=testfloat[:]
qn=normalizeQuantiles(sa)
qn=convert(Array{Float64},reshape([get(qn[i]) for i=1:length(qn)],size(qn)))
@test qn == check

testint = [ 1 1 1 ; 1 1 1 ; 1 1 1 ]
qn=normalizeQuantiles(testint)
@test qn == testint
sa=SharedArray(Nullable{Int},(size(testint,1),size(testint,2)));
sa[:]=testint[:]
qn=normalizeQuantiles(sa)
qn=convert(Array{Float64},reshape([get(qn[i]) for i=1:length(qn)],size(qn)))
@test qn == testint
sa=SharedArray(Int,(size(testint,1),size(testint,2)));
sa[:]=testint[:]
qn=normalizeQuantiles(sa)

dafloat=Array{Nullable{Float64}}(testfloat)
dafloat[2,2]=Nullable{Float64}()
qn=normalizeQuantiles(dafloat)
@test isnull(qn[2,2])
@test get(qn[1,2])==3.5
@test get(qn[2,1])==5.0
sa=SharedArray(Nullable{Float64},(size(dafloat,1),size(dafloat,2)));
sa[:]=dafloat[:]
qn=normalizeQuantiles(sa)
@test isnull(qn[2,2])
@test get(qn[1,2])==3.5
@test get(qn[2,1])==5.0

dafloat[2,:]=Nullable{Float64}()
qn=normalizeQuantiles(dafloat)
@test isnull(qn[2,1])
@test isnull(qn[2,2])
@test isnull(qn[2,3])
sa=SharedArray(Nullable{Float64},(size(dafloat,1),size(dafloat,2)));
sa[:]=dafloat[:]
qn=normalizeQuantiles(sa)
@test isnull(qn[2,1])
@test isnull(qn[2,2])
@test isnull(qn[2,3])

dafloat[3,1:2]=Nullable{Float64}()
qn=normalizeQuantiles(dafloat)
@test isnull(qn[3,1])
@test isnull(qn[3,2])
sa=SharedArray(Nullable{Float64},(size(dafloat,1),size(dafloat,2)));
sa[:]=dafloat[:]
qn=normalizeQuantiles(sa)
@test isnull(qn[3,1])
@test isnull(qn[3,2])

dafloat[1,:]=Nullable{Float64}()
dafloat[2,:]=Nullable{Float64}()
dafloat[3,:]=Nullable{Float64}()
dafloat[4,:]=Nullable{Float64}()
qn = normalizeQuantiles(dafloat)
@test isnull(qn[1,1])
@test isnull(qn[1,2])
@test isnull(qn[1,3])
@test isnull(qn[2,1])
@test isnull(qn[2,2])
@test isnull(qn[2,3])
@test isnull(qn[3,1])
@test isnull(qn[3,2])
@test isnull(qn[3,3])
@test isnull(qn[4,1])
@test isnull(qn[4,2])
@test isnull(qn[4,3])
sa=SharedArray(Nullable{Float64},(size(dafloat,1),size(dafloat,2)));
sa[:]=dafloat[:]
qn=normalizeQuantiles(sa)
@test isnull(qn[1,1])
@test isnull(qn[1,2])
@test isnull(qn[1,3])
@test isnull(qn[2,1])
@test isnull(qn[2,2])
@test isnull(qn[2,3])
@test isnull(qn[3,1])
@test isnull(qn[3,2])
@test isnull(qn[3,3])
@test isnull(qn[4,1])
@test isnull(qn[4,2])
@test isnull(qn[4,3])

testfloat = [ 2.0 2.0 8.0 0.0 7.0 ]
a=Array(Nullable{Float64},(size(testfloat,1),size(testfloat,2)));
a[:]=testfloat[:]
a[4]=Nullable{Float64}()
(r,m)=sampleRanks(a,tiesMethod=tmMin,naIncreasesRank=true,resultMatrix=true)
r[4]=Nullable{Int}(0)
r=[ Int(get(x)) for x in r ]
@test r==Array{Int}([1,1,4,0,2])

testfloat = [ 2.0 2.0 8.0 0.0 7.0 ]
a=Array(Nullable{Float64},(size(testfloat,1),size(testfloat,2)));
a[:]=testfloat[:]
a[4]=Nullable{Float64}()
(r,m)=sampleRanks(a,tiesMethod=tmOrder,naIncreasesRank=true,resultMatrix=true)
r[4]=Nullable{Int}(0)
r=[ Int(get(x)) for x in r ]
@test r==Array{Int}([1,2,5,0,3])

testfloat = [ 2.0 2.0 8.0 0.0 7.0 ]
a=Array(Nullable{Float64},(size(testfloat,1),size(testfloat,2)));
a[:]=testfloat[:]
a[4]=Nullable{Float64}()
(r,m)=sampleRanks(a,tiesMethod=tmMin,naIncreasesRank=false,resultMatrix=true)
r[4]=Nullable{Int}(0)
r=[ Int(get(x)) for x in r ]
@test r==Array{Int}([1,1,3,0,2])

testfloat = [ 5.0 2.0 4.0 3.0 1.0 ]
a=Array(Nullable{Float64},(size(testfloat,1),size(testfloat,2)));
a[:]=testfloat[:]
(r,m)=sampleRanks(a,tiesMethod=tmMin,naIncreasesRank=true,resultMatrix=true)
r=[ Int(get(x)) for x in r ]
@test r==Array{Int}([5,2,4,3,1])

testfloat = [ 2.0 2.0 0.0 2.0 2.0 ]
a=Array(Nullable{Float64},(size(testfloat,1),size(testfloat,2)));
a[:]=testfloat[:]
a[3]=Nullable{Float64}()
(r,m)=sampleRanks(a,tiesMethod=tmMin,naIncreasesRank=true,resultMatrix=true)
r[3]=Nullable{Int}(0)
r=[ Int(get(x)) for x in r ]
@test r==Array{Int}([1,1,0,1,1])

testfloat = [ 2.0 2.0 0.0 2.0 4.0 ]
a=Array(Nullable{Float64},(size(testfloat,1),size(testfloat,2)));
a[:]=testfloat[:]
a[3]=Nullable{Float64}()
(r,m)=sampleRanks(a,tiesMethod=tmMin,naIncreasesRank=true,resultMatrix=true)
r[3]=Nullable{Int}(0)
r=[ Int(get(x)) for x in r ]
@test r==Array{Int}([1,1,0,1,3])

testfloat = [ 2.0 2.0 0.0 2.0 4.0 ]
a=Array(Nullable{Float64},(size(testfloat,1),size(testfloat,2)));
a[:]=testfloat[:]
a[3]=Nullable{Float64}()
(r,m)=sampleRanks(a,tiesMethod=tmMin,naIncreasesRank=false,resultMatrix=true)
r[3]=Nullable{Int}(0)
r=[ Int(get(x)) for x in r ]
@test r==Array{Int}([1,1,0,1,2])

testfloat = [ 2.0 2.0 0.0 3.0 4.0 ]
a=Array(Nullable{Float64},(size(testfloat,1),size(testfloat,2)));
a[:]=testfloat[:]
a[3]=Nullable{Float64}()
(r,m)=sampleRanks(a,tiesMethod=tmMin,naIncreasesRank=true,resultMatrix=true)
r[3]=Nullable{Int}(0)
r=[ Int(get(x)) for x in r ]
@test r==Array{Int}([1,1,0,3,4])

testfloat = [ 2.0 2.0 0.0 3.0 4.0 ]
a=Array(Nullable{Float64},(size(testfloat,1),size(testfloat,2)));
a[:]=testfloat[:]
a[3]=Nullable{Float64}()
(r,m)=sampleRanks(a,tiesMethod=tmMin,naIncreasesRank=false,resultMatrix=true)
r[3]=Nullable{Int}(0)
r=[ Int(get(x)) for x in r ]
@test r==Array{Int}([1,1,0,2,3])

testfloat = [ 0.0 2.0 5.0 3.0 4.0 ]
a=Array(Nullable{Float64},(size(testfloat,1),size(testfloat,2)));
a[:]=testfloat[:]
a[1]=Nullable{Float64}()
(r,m)=sampleRanks(a,tiesMethod=tmMin,naIncreasesRank=true,resultMatrix=true)
r[1]=Nullable{Int}(0)
r=[ Int(get(x)) for x in r ]
@test r==Array{Int}([0,2,5,3,4])

testfloat = [ 0.0 2.0 5.0 3.0 4.0 ]
a=Array(Nullable{Float64},(size(testfloat,1),size(testfloat,2)));
a[:]=testfloat[:]
a[1]=Nullable{Float64}()
(r,m)=sampleRanks(a,tiesMethod=tmMin,naIncreasesRank=false,resultMatrix=true)
r[1]=Nullable{Int}(0)
r=[ Int(get(x)) for x in r ]
@test r==Array{Int}([0,1,4,2,3])

testfloat = [ 2.0 2.0 2.0 2.0 2.0 ]
a=Array(Nullable{Float64},(size(testfloat,1),size(testfloat,2)));
a[:]=testfloat[:]
(r,m)=sampleRanks(a,tiesMethod=tmMin,naIncreasesRank=true,resultMatrix=true)
r=[ Int(get(x)) for x in r ]
@test r==Array{Int}([1,1,1,1,1])

testfloat = [ 2.0 2.0 2.0 2.0 2.0 ]
a=Array(Nullable{Float64},(size(testfloat,1),size(testfloat,2)));
a[:]=testfloat[:]
(r,m)=sampleRanks(a,tiesMethod=tmReverse,naIncreasesRank=true,resultMatrix=true)
r=[ Int(get(x)) for x in r ]
@test r==Array{Int}([5,4,3,2,1])

end # if VERSION >= v"0.4.0-"

if VERSION < v"0.4.0-"

using DataArrays

testfloat = [ 3.0 2.0 8.0 1.0 ; 4.0 5.0 6.0 2.0 ; 9.0 7.0 8.0 3.0 ; 5.0 2.0 8.0 4.0 ]
dafloat=DataArray(testfloat)
r=normalizeQuantiles(dafloat)
@test mean(r[:,1]) >= 4.8124 && mean(r[:,1]) <= 4.8126
@test mean(r[:,2]) >= 4.8124 && mean(r[:,2]) <= 4.8126
@test mean(r[:,3]) >= 4.8124 && mean(r[:,3]) <= 4.8126
@test mean(r[:,4]) >= 4.8124 && mean(r[:,4]) <= 4.8126
sa=SharedArray(Float64,(size(dafloat,1),size(dafloat,2)));
sa[:]=dafloat[:]
r=normalizeQuantiles(sa)
@test mean(r[:,1]) >= 4.8124 && mean(r[:,1]) <= 4.8126
@test mean(r[:,2]) >= 4.8124 && mean(r[:,2]) <= 4.8126
@test mean(r[:,3]) >= 4.8124 && mean(r[:,3]) <= 4.8126
@test mean(r[:,4]) >= 4.8124 && mean(r[:,4]) <= 4.8126

testfloat = [ 3.5 2.0 8.1 1.0 ; 4.5 5.0 6.0 2.0 ; 9.0 7.6 8.2 3.0 ; 5.0 2.0 8.0 4.0 ]
dafloat=DataArray(testfloat)
r=normalizeQuantiles(dafloat)
@test mean(r[:,1]) >= 4.93124 && mean(r[:,1]) <= 4.93125
@test mean(r[:,2]) >= 4.93124 && mean(r[:,2]) <= 4.93125
@test mean(r[:,3]) >= 4.93124 && mean(r[:,3]) <= 4.93125
@test mean(r[:,4]) >= 4.93124 && mean(r[:,4]) <= 4.93125
sa=SharedArray(Float64,(size(dafloat,1),size(dafloat,2)));
sa[:]=dafloat[:]
r=normalizeQuantiles(sa)
@test mean(r[:,1]) >= 4.93124 && mean(r[:,1]) <= 4.93125
@test mean(r[:,2]) >= 4.93124 && mean(r[:,2]) <= 4.93125
@test mean(r[:,3]) >= 4.93124 && mean(r[:,3]) <= 4.93125
@test mean(r[:,4]) >= 4.93124 && mean(r[:,4]) <= 4.93125

testfloat = [ 3.0 2.0 1.0 ; 4.0 5.0 6.0 ; 9.0 7.0 8.0 ; 5.0 2.0 8.0 ]
check = [ 2.0 3.0 2.0 ; 4.0 6.0 4.0 ; 8.0 8.0 7.0 ; 6.0 3.0 7.0 ]
qn = normalizeQuantiles(testfloat)
@test qn == check
sa=SharedArray(Float64,(size(testfloat,1),size(testfloat,2)));
sa[:]=testfloat[:]
sacheck=SharedArray(Float64,(size(check,1),size(check,2)));
sacheck[:]=check[:]
qn = normalizeQuantiles(sa)
@test qn == sacheck

testint = [ 1 1 1 ; 1 1 1 ; 1 1 1 ]
qn = normalizeQuantiles(testint)
@test qn == testint

dafloat = DataArray(testfloat)
dacheck = DataArray(check)
qn = normalizeQuantiles(dafloat)
@test qn == check
sa=SharedArray(Float64,(size(dafloat,1),size(dafloat,2)));
sa[:]=dafloat[:]
sacheck=SharedArray(Float64,(size(dacheck,1),size(dacheck,2)));
sacheck[:]=dacheck[:]
qn = normalizeQuantiles(sa)
@test qn == sacheck

daint = DataArray(testint)
qn = normalizeQuantiles(daint)
@test qn == daint

dafloat[2,2]=NA
qn = normalizeQuantiles(dafloat)
@test isa(qn[2,2],NAtype)
@test qn[1,2]==3.5
@test qn[2,1]==5.0

dafloat[2,:]=NA
qn = normalizeQuantiles(dafloat)
@test isa(qn[2,1],NAtype)
@test isa(qn[2,2],NAtype)
@test isa(qn[2,3],NAtype)

dafloat[3,1:2]=NA
qn = normalizeQuantiles(dafloat)
@test isa(qn[3,1],NAtype)
@test isa(qn[3,2],NAtype)

dafloat[1,:]=NA
dafloat[2,:]=NA
dafloat[3,:]=NA
dafloat[4,:]=NA
qn = normalizeQuantiles(dafloat)
@test isa(qn[1,1],NAtype)
@test isa(qn[1,2],NAtype)
@test isa(qn[1,3],NAtype)
@test isa(qn[2,1],NAtype)
@test isa(qn[2,2],NAtype)
@test isa(qn[2,3],NAtype)
@test isa(qn[3,1],NAtype)
@test isa(qn[3,2],NAtype)
@test isa(qn[3,3],NAtype)
@test isa(qn[4,1],NAtype)
@test isa(qn[4,2],NAtype)
@test isa(qn[4,3],NAtype)

end # if VERSION < v"0.4.0-"


text:
abstract type AbstractTikzPlot end
abstract type AbstractTikzAxis end
abstract type AbstractTikzFigure end


text:
#==============================================================================#
# AWSSNS.jl
#
# SNS API. See http://aws.amazon.com/documentation/sns/
#
# Copyright OC Technology Pty Ltd 2014 - All rights reserved
#==============================================================================#

__precompile__()

module AWSSNS

export sns_list_topics, sns_delete_topic, sns_create_topic, sns_subscribe_sqs,
       sns_subscribe_email, sns_subscribe_lambda, sns_publish,
       sns_list_subscriptions, sns_unsubscribe, send_sms

using AWSCore
using SymDict
using Retry
using AWSSQS

sns_arn(aws::AWSConfig, topic_name) = arn(aws, "sns", topic_name)

sns(aws::AWSConfig, action, args) = AWSCore.Services.sns(aws, action, args)

function sns(aws::AWSConfig, action::String, topic::String; args...)

    sns(aws, action, merge(stringdict(args),
                           Dict("Name" => topic,
                                "TopicArn" => sns_arn(aws, topic))))
end

"""
    sns_list_topics([::AWSConfig])

Returns a list of topic names.
"""
function sns_list_topics(aws::AWSConfig=default_aws_config())
    l = String[]
    while true
        if isempty(l)
            r = sns(aws, "ListTopics", [])
        else
            r = sns(aws, "ListTopics", ["NextToken" => r["NextToken"]])
        end
        for t in r["Topics"]
            push!(l, String(split(t["TopicArn"],":")[6]))
        end
        if !haskey(r, "NextToken") || r["NextToken"] == nothing
            break
        end
    end
    return l
end

"""
   sns_create_topic([::AWSConfig], topic_name)

Create a named topic.
"""
function sns_create_topic(aws::AWSConfig, topic_name)
    sns(aws, "CreateTopic", topic_name)
end

sns_create_topic(name) = sns_create_topic(default_aws_config(), name)

"""
   sns_delete_topic([::AWSConfig], topic_name)

Delete a named topic.
"""
function sns_delete_topic(aws::AWSConfig, topic_name)
    sns(aws, "DeleteTopic", topic_name)
end

sns_delete_topic(name) = sns_delete_topic(default_aws_config(), name)

"""
    sns_publish([::AWSConfig], topic_name, message, subject="No Subject")

Send a `message` to a named topic (with optional `subject`).
"""
function sns_publish(aws::AWSConfig, topic_name, message, subject="No Subject")

    if length(subject) > 100
        subject = subject[1:100]
    end
    sns(aws, "Publish", topic_name, Message = message, Subject = subject)
end

sns_publish(a...) = sns_publish(default_aws_config(), a...)

"""
    send_sms([::AWSConfig], number, message)

Send SMS `message` to `number`.
"""
function send_sms(aws::AWSConfig, number, message)

    sns(aws, "Publish", ["PhoneNumber" => number,
                         "Message" => message])
end

send_sms(number, message) = send_sms(default_aws_config(), number, message)

"""
    sns_subscribe_sqs([::AWSConfig], topic_name, queue; raw=false)

Connect SQS `queue` to `topic_name`.
"""
function sns_subscribe_sqs(aws::AWSConfig, topic_name, queue; raw=false)

    r = sns(aws, "Subscribe", ["Name" => topic_name,
                               "TopicArn" => sns_arn(aws, topic_name),
                               "Endpoint" => sqs_arn(queue),
                               "Protocol" => "sqs"])

    if raw
        sns(aws, "SetSubscriptionAttributes", topic_name,
                  SubscriptionArn = r["SubscriptionArn"],
                  AttributeName = "RawMessageDelivery",
                  AttributeValue = "true")
    end

#=
    import AWSSQS: sqs, sqs_get_queue

    This is probably a bad idea.
    It might wipe out a policy already attached to the queue

    q = sqs_get_queue(aws, queue)

    sqs(q, Dict(
        "Action" => "SetQueueAttributes",
        "Attribute.Name" => "Policy",
        "Attribute.Value" => """{
          "Version": "2008-10-17",
          "Statement": [
            {
              "Effect": "Allow",
              "Principal": {
                "AWS": "*"
              },
              "Action": "SQS:SendMessage",
              "Resource": "$(arn(aws, "sqs", queue))",
              "Condition": {
                "ArnEquals": {
                  "aws:SourceArn": "$(sns_arn(aws, topic_name))"
                }
              }
            }
          ]
        }"""
    ))
=#
end

function sns_subscribe_sqs(a...; b...)
    sns_subscribe_sqs(default_aws_config(), a...; b...)
end

"""
    sns_subscribe_email([::AWSConfig], topic_name, email)

Connect `email` to `topic_name`.
"""
function sns_subscribe_email(aws::AWSConfig, topic_name, email)

    sns(aws, "Subscribe", topic_name, Endpoint = email, Protocol = "email")
end

sns_subscribe_email(a...) = sns_subscribe_email(default_aws_config(), a...)

import AWSLambda: lambda

"""
    sns_subscribe_lambda([::AWSConfig], topic_name, lambda_name)

Connect `lambda_name` to `topic_name`.
"""
function sns_subscribe_lambda(aws::AWSConfig, topic_name, lambda_name)

    if occursin(r"^arn", lambda_name)
        lambda_arn = lambda_name
    else
        lambda_arn = arn(aws, "lambda", "function:$lambda_name")
    end

    sns(aws, "Subscribe", topic_name, Endpoint = lambda_arn, Protocol = "lambda")

# FIXME does not scale. End up with:
#   ERROR: LoadError: 400 -- The final policy size is bigger than the limit.
#=
    @protected try

        lambda(laws, "POST"; path="$lambda_name/policy", query=Dict(
               "Action" => "lambda:InvokeFunction",
               "Principal" => "sns.amazonaws.com",
               "SourceArn" => sns_arn(aws, topic_name),
               "StatementId" => "sns_$topic_name"))
    catch e
        @ignore if e.code == "409" end
    end
=#
end

sns_subscribe_lambda(a...) = sns_subscribe_lambda(default_aws_config(), a...)

"""
    sns_list_subscriptions([::AWSConfig], topic_name)

List endpoints that are subscribed to `topic_name`.
"""
function sns_list_subscriptions(aws::AWSConfig, topic_name)

    r = sns(aws, "ListSubscriptionsByTopic", topic_name)
    r["Subscriptions"]
end

sns_list_subscriptions(a) = sns_list_subscriptions(default_aws_config(), a)

"""
    sns_unsubscribe([::AWSConfig], topic_name, SubscriptionArn)
    sns_unsubscribe([::AWSConfig], topic_name, ::Regex)

Disconnect `SubscriptionArn` from `topic_name`.
"""
function sns_unsubscribe(aws::AWSConfig, topic_name, SubscriptionArn)

    sns(aws, "Unsubscribe", topic_name, SubscriptionArn = SubscriptionArn)
end

function sns_unsubscribe(aws::AWSConfig, topic_name, pattern::Regex)

    for s in sns_list_subscriptions(aws, topic_name)
        if occursin(pattern, s["Endpoint"])
            sns_unsubscribe(aws, topic_name, s["SubscriptionArn"])
        end
    end
end

sns_unsubscribe(a...) = sns_unsubscribe(default_aws_config(), a...)

end # module AWSSNS

#==============================================================================#
# End of file.
#==============================================================================#


text:
struct SumMean{T<:NTuple{N,Mean} where N} <: CompositeMean
    means::T
end

SumMean(means::Mean...) = SumMean(means)

mean(sm::SumMean, x::AbstractVector) = sum(mean(m, x) for m in components(sm))

get_param_names(sm::SumMean) = composite_param_names(components(sm), :sm)

function grad_mean(sm::SumMean, x::AbstractVector)
    dm = Array{Float64}(undef, num_params(sm))
    v = 1
    for m in components(sm)
        w = v + num_params(m)
        dm[v:(w - 1)] = grad_mean(m, x)
        v = w
    end
    dm
end

Base.:+(m1::SumMean, m2::Mean) = SumMean(m1.means..., m2)
Base.:+(m1::SumMean, m2::SumMean) = SumMean(m1.means..., m2.means...)
Base.:+(m1::Mean, m2::Mean) = SumMean(m1, m2)
Base.:+(m1::Mean, m2::SumMean) = SumMean(m1, m2.means...)


text:
"""
    Barrier(num_threads)
Create a thread-safe barrier that blocks all tasks until `num_threads`
tasks call `sync`.
"""
mutable struct Barrier
    notify::Threads.Condition
    tar_cnt::Int
    cur_cnt::Int
    Barrier(n_threads) = n_threads > 0 ? new(Threads.Condition(), n_threads, 0) : throw(ArgumentError("Barrier size must be > 0"))
end

"""
    sync(b::Barrier)

Stops all threads until they all reached the barrier and call `sync`.

See @Barrier
"""
function sync(b::Barrier)
    lock(b.notify)
    try
        b.cur_cnt += 1

        if b.cur_cnt == b.tar_cnt
            notify(b.notify, true)
            b.cur_cnt = 0
        else
            wait(b.notify)
        end
    finally
        unlock(b.notify)
    end
    return nothing
end


text:
@doc doc"""
Résolution des problèmes de minimisation sous contraintes d'égalités

# Syntaxe
```julia
Lagrangien_Augmente(algo,fonc,contrainte,gradfonc,hessfonc,grad_contrainte,
			hess_contrainte,x0,options)
```

# Entrées
  * **algo** 		   : (String) l'algorithme sans contraintes à utiliser:
    - **"newton"**  : pour l'algorithme de Newton
    - **"cauchy"**  : pour le pas de Cauchy
    - **"gct"**     : pour le gradient conjugué tronqué
  * **fonc** 		   : (Function) la fonction à minimiser
  * **contrainte**	   : (Function) la contrainte [x est dans le domaine des contraintes ssi ``c(x)=0``]
  * **gradfonc**       : (Function) le gradient de la fonction
  * **hessfonc** 	   : (Function) la hessienne de la fonction
  * **grad_contrainte** : (Function) le gradient de la contrainte
  * **hess_contrainte** : (Function) la hessienne de la contrainte
  * **x0** 			   : (Array{Float,1}) la première composante du point de départ du Lagrangien
  * **options**		   : (Array{Float,1})
    1. **epsilon** 	   : utilisé dans les critères d'arrêt
    2. **tol**         : la tolérance utilisée dans les critères d'arrêt
    3. **itermax** 	   : nombre maximal d'itération dans la boucle principale
    4. **lambda0**	   : la deuxième composante du point de départ du Lagrangien
    5. **mu0,tho** 	   : valeurs initiales des variables de l'algorithme

# Sorties
* **xmin**		   : (Array{Float,1}) une approximation de la solution du problème avec contraintes
* **fxmin** 	   : (Float) ``f(x_{min})``
* **flag**		   : (Integer) indicateur du déroulement de l'algorithme
   - **0**    : convergence
   - **1**    : nombre maximal d'itération atteint
   - **(-1)** : une erreur s'est produite
* **niters** 	   : (Integer) nombre d'itérations réalisées

# Exemple d'appel
```julia
using LinearAlgebra
f(x)=100*(x[2]-x[1]^2)^2+(1-x[1])^2
gradf(x)=[-400*x[1]*(x[2]-x[1]^2)-2*(1-x[1]) ; 200*(x[2]-x[1]^2)]
hessf(x)=[-400*(x[2]-3*x[1]^2)+2  -400*x[1];-400*x[1]  200]
algo = "gct" # ou newton|gct
x0 = [1; 0]
options = []
contrainte(x) =  (x[1]^2) + (x[2]^2) -1.5
grad_contrainte(x) = [2*x[1] ;2*x[2]]
hess_contrainte(x) = [2 0;0 2]
output = Lagrangien_Augmente(algo,f,contrainte,gradf,hessf,grad_contrainte,hess_contrainte,x0,options)
```
"""
function Lagrangien_Augmente(algo,fonc::Function,contrainte::Function,gradfonc::Function,
	hessfonc::Function,grad_contrainte::Function,hess_contrainte::Function,x0,options)

	if options == []
		epsilon = 1e-8
		tol = 1e-5
		itermax = 1000
		lambda0 = 2
		mu0 = 100
		tho = 2
	else
		epsilon = options[1]
		tol = options[2]
		itermax = options[3]
		lambda0 = options[4]
		mu0 = options[5]
		tho = options[6]
    end

    # Initialisation des parametres
	mu = mu0
	lambda = lambda0
	betha = 0.9
	alpha = 0.1
	etha_c = 0.1258925
	epsilon0 = 1
    etha = etha_c/(mu^alpha)

    xk = x0
    xk_1 = x0
    iter = 0
    flag = -1

    while (true)
        iter +=1
        # Definition de la fonction Lagrangien
        Lx(x) = fonc(x) + lambda'*contrainte(x) + (mu/2)*(norm(contrainte(x)))^2
        grad_x_Lx(x) = gradfonc(x) + lambda'*grad_contrainte(x) + mu*grad_contrainte(x)*contrainte(x)
        Hess_x_Lx(x) = hessfonc(x) + lambda'*hess_contrainte(x) + mu*hess_contrainte(x)*contrainte(x) +
            mu*grad_contrainte(x)*grad_contrainte(x)'

        # Calcul de xk
        if algo == "newton"
          option = [itermax,epsilon,tol]
          xk, ~ = Algorithme_De_Newton(Lx,grad_x_Lx,Hess_x_Lx,xk,option)
        elseif algo == "cauchy" || algo == "gct"
          option = [10,0.5,2,0.25,0.75,2,1000,epsilon,tol]
          xk, ~ = Regions_De_Confiance(algo,Lx,grad_x_Lx,Hess_x_Lx,xk,[])
        end

        " Test convergence globale"
        gradL(x,l) = gradfonc(x) + l'*grad_contrainte(xk)
        if norm(gradL(xk,lambda)) <= max(epsilon, tol*norm(gradL(x0,lambda0)))
          flag = 0
          break
        elseif norm(contrainte(xk)) <= max(tol, tol*norm(contrainte(x0)))
          flag = 0
          break
        elseif iter == itermax
          flag = 1
          break
        else # pas de convergence globale
          if norm(contrainte(xk))<=etha
            lambda = lambda + mu*contrainte(xk)
            epsilon = epsilon/mu
            etha = etha/(mu^betha)
          else
            mu = tho*mu
            epsilon = epsilon0/mu
            etha = etha_c/(mu^alpha)
          end # if

        end # if

    end # while
	xmin = xk
	fxmin = fonc(xmin)

	return xmin,fxmin,flag,iter
end

max_stars_repo_path: 2A/S7/UE Optimisation et RO/Optimisation/Optinum/src/Lagrangien_Augmente.jl

max_issues_repo_name: Sajid-Badr/ENSEEIHT

max_forks_repo_name: Sajid-Badr/ENSEEIHT


text:
using GeoPhyInv
using LinearAlgebra
using Random

# # Intro

#md # ```@docs
#md # SrcWav
#md # ```

# # Examples
# Define an acquisition geometry.
mgrid = [range(0.0, stop=10.,step=0.1), range(0.0, stop=30.,step=0.2)];
ageom=AGeom(mgrid, SSrcs(10), Srcs(10), Recs(10));

# Need a time grid.
tgrid=range(0, stop=1.0, step=0.1);

# Lets initialize records for `:P` field.
srcwav=SrcWav(tgrid, ageom, [:P]);

# Fill the `:P` field of 3rd supersource with random numbers.
Random.randn!(srcwav[3][:P]);

# Often we want to populate the same source wavelet to all
# the supersources and sources.
x=randn(length(tgrid));
update!(srcwav, [:P,], x);

# Populate two different wavelets for first and second supersources.
x1=randn(length(tgrid));
x2=randn(length(tgrid));
update!(srcwav[1], [:P,], x1);
update!(srcwav[2], [:P,], x2);

# Scale `srcwav` by a scalar overwriting it in-place.
rmul!(srcwav, 2.0);

# # Methods
# Most of the methods listed below are also applicable to individual elements of `srcwav`.

#md # ```@docs
#md # update!(GeoPhyInv.VNamedD, ::Vector{Symbol}, ::AbstractArray)
#md # Base.reverse!(::GeoPhyInv.VNamedD)
#md # Base.iszero(::GeoPhyInv.VNamedD)
#md # Base.isequal(::GeoPhyInv.VNamedD, ::GeoPhyInv.VNamedD)
#md # GeoPhyInv.issimilar(::GeoPhyInv.VNamedD, ::GeoPhyInv.VNamedD)
#md # Base.vec(::GeoPhyInv.VNamedD)
#md # Random.randn!(::GeoPhyInv.VNamedD)
#md # Base.fill!(::GeoPhyInv.VNamedD, ::Float64)
#md # LinearAlgebra.dot(::GeoPhyInv.VNamedD, ::GeoPhyInv.VNamedD)
#md # LinearAlgebra.rmul!(::GeoPhyInv.VNamedD, ::Number)
#md # Base.copyto!(::GeoPhyInv.VNamedD, ::GeoPhyInv.VNamedD)
#md # Base.copyto!(::AbstractVector{Float64}, ::GeoPhyInv.VNamedD)
#md # Base.copyto!(::GeoPhyInv.VNamedD, ::AbstractVector{Float64})
#md # ```


text:
## Julia based

using FlashWeave

# Load  metadata files
meta_data_path = "files/FlashWeave_data/metadata_flashweave_friendly.csv"

# Load SPECIES data files
species_data_path = "files/FlashWeave_data/flashweave_friendly_species_relab.csv"
####### Species Model #######
species_results = learn_network(species_data_path, meta_data_path, sensitive=true, heterogeneous=false, verbose=true, feed_forward= true, alpha=0.05, FDR=true)
save_network("data/Co_Abundance_Analysis/species_network_output.edgelist", species_results)
save_network("data/Co_Abundance_Analysis/species_network_output.gml", species_results)

# Load GENUS data files
genus_data_path = "files/FlashWeave_data/flashweave_friendly_genus_relab.csv"
####### Species Model #######
genus_results = learn_network(genus_data_path, meta_data_path, sensitive=true, heterogeneous=false, verbose=true, feed_forward= true, alpha=0.05, FDR=true)
save_network("data/Co_Abundance_Analysis/genus_network_output.edgelist", genus_results)
save_network("data/Co_Abundance_Analysis/genus_network_output.gml", genus_results)

# Load PHYLUM data files
phylum_data_path = "files/FlashWeave_data/flashweave_friendly_phylum_relab.csv"
####### Species Model #######
phylum_results = learn_network(phylum_data_path, meta_data_path, sensitive=true, heterogeneous=false, verbose=true, feed_forward= true, alpha=0.05, FDR=true)
save_network("data/Co_Abundance_Analysis/phylum_network_output.edgelist", phylum_results)
save_network("data/Co_Abundance_Analysis/phylum_network_output.gml", phylum_results)

# Load PATHWAY data files
pathways_data_path = "files/FlashWeave_data/flashweave_friendly_pathways_relab.csv"
####### PATHWAY Model #######
pathways_results = learn_network(pathways_data_path, meta_data_path, sensitive=true, heterogeneous=false, verbose=true, feed_forward= true, alpha=0.05, FDR=true)
save_network("data/Co_Abundance_Analysis/pathways_network_output.edgelist", pathways_results)
save_network("data/Co_Abundance_Analysis/pathways_network_output.gml", pathways_results)

# Load ENZYME data files
enzymes_data_path = "files/FlashWeave_data/flashweave_friendly_enzymes_relab.csv"
####### ENZYME Model #######
enzymes_results = learn_network(enzymes_data_path, meta_data_path, sensitive=true, heterogeneous=true, verbose=true, feed_forward= true, alpha=0.05, FDR=true)
save_network("data/Co_Abundance_Analysis/enzymes_network_output.edgelist", enzymes_results)
save_network("data/Co_Abundance_Analysis/enzymes_network_output.gml", enzymes_results)

# Load Kegg Ortholog data files
kos_data_path = "files/FlashWeave_data/flashweave_friendly_kos_relab.csv"
####### KO Model #######
kos_results = learn_network(kos_data_path, meta_data_path, sensitive=true, heterogeneous=true, verbose=true, feed_forward= true, alpha=0.05, FDR=true)
save_network("data/Co_Abundance_Analysis/kos_network_output.edgelist", kos_results)
save_network("data/Co_Abundance_Analysis/kos_network_output.gml", kos_results)


text:
import LinearAlgebra.I

idmatrix(n::Int) = Matrix{Float64}(I, n, n)

seq(start, stop; length = nothing) = range(start, stop = stop, length = length)


text:
Random.seed!(209)

@testset "Deprecated spaces" begin
    # GraphSpace
    @test_deprecated (@test_throws TypeError Space(Agents.Graph(1)))
    # GridSpace
    @test_deprecated (@test_throws TypeError Space((5, 1)))
    # ContinuousSpace
    @test_deprecated Space(5, (a, m) -> nothing)
end

@testset "0D grids" begin
    @test GridSpace((1,)).graph == Agents.Graph(1)
end

@testset "1D grids" begin
    a = GridSpace((5, 1))
    ae = collect(Agents.LightGraphs.edges(a.graph))

    b = GridSpace((5, 1), periodic = true)
    be = collect(Agents.LightGraphs.edges(b.graph))

    @test ae == [
        Agents.LightGraphs.Edge(1, 2),
        Agents.LightGraphs.Edge(2, 3),
        Agents.LightGraphs.Edge(3, 4),
        Agents.LightGraphs.Edge(4, 5),
    ]
    @test be == [
        Agents.LightGraphs.Edge(1, 2),
        Agents.LightGraphs.Edge(1, 5),
        Agents.LightGraphs.Edge(2, 3),
        Agents.LightGraphs.Edge(3, 4),
        Agents.LightGraphs.Edge(4, 5),
    ]
end

@testset "2D grids" begin
    @test GridSpace((2, 1)).graph == Agents.Graph(2, 1)

    a = GridSpace((2, 3))
    b = GridSpace((2, 3), periodic = true) # 2D grid

    @test a.dimensions == (2, 3)
    @test b.dimensions == (2, 3)
    @test length(a.agent_positions) == 6
    @test length(b.agent_positions) == 6

    @test nv(a) == 6
    @test ne(a) == 7
    @test nv(b) == 6
    @test ne(b) == 9

    ae = collect(Agents.LightGraphs.edges(a.graph))
    be = collect(Agents.LightGraphs.edges(b.graph))

    @test ae == [
        Agents.LightGraphs.Edge(1, 2),
        Agents.LightGraphs.Edge(1, 3),
        Agents.LightGraphs.Edge(2, 4),
        Agents.LightGraphs.Edge(3, 4),
        Agents.LightGraphs.Edge(3, 5),
        Agents.LightGraphs.Edge(4, 6),
        Agents.LightGraphs.Edge(5, 6),
    ]
    @test be == [
        Agents.LightGraphs.Edge(1, 2),
        Agents.LightGraphs.Edge(1, 3),
        Agents.LightGraphs.Edge(1, 5),
        Agents.LightGraphs.Edge(2, 4),
        Agents.LightGraphs.Edge(2, 6),
        Agents.LightGraphs.Edge(3, 4),
        Agents.LightGraphs.Edge(3, 5),
        Agents.LightGraphs.Edge(4, 6),
        Agents.LightGraphs.Edge(5, 6),
    ]
end

@testset "2D Moore" begin
    a = GridSpace((3, 3), moore = true)
    @test Agents.nv(a) == 9
    @test Agents.ne(a) == 20

    a = GridSpace((3, 4), moore = true)
    @test Agents.nv(a) == 12
    @test Agents.ne(a) == 29

    b = GridSpace((3, 2), moore = true, periodic = true)
    @test Agents.nv(b) == 6
    @test Agents.ne(b) == 15

    b = GridSpace((3, 3), moore = true, periodic = true)
    @test Agents.nv(b) == 9
    @test Agents.ne(b) == 36
end

@testset "3D grid" begin
    g1 = GridSpace((2, 3, 2))
    g2 = GridSpace((2, 3, 3))
    g3 = GridSpace((2, 3, 2), periodic = true)
    g4 = GridSpace((2, 3, 3), periodic = true)
    @test Agents.ne(g1) == 20
    @test Agents.ne(g2) == 33
    @test Agents.ne(g3) == 24
    @test Agents.ne(g4) == 39
end

@testset "grid coord/vertex conversions" begin
    @test Agents.coord2vertex((2, 2, 1), (3, 4, 1)) == 5
    @test Agents.coord2vertex((2, 2), (3, 4)) == 5
    @test Agents.coord2vertex((2, 2, 1), (3, 4, 1)) == 5
    @test Agents.coord2vertex((2, 2), (3, 4)) == 5
    @test Agents.vertex2coord(5, (3, 4, 1)) == (2, 2, 1)
    @test Agents.vertex2coord(5, (3, 4)) == (2, 2)

    @test Agents.coord2vertex((2, 2, 1), (2, 3, 3)) == 4
    @test Agents.coord2vertex((2, 2, 2), (2, 3, 3)) == 10
    @test Agents.coord2vertex((2, 2, 3), (2, 3, 3)) == 16
    @test Agents.coord2vertex(((1, 3, 1)), (2, 3, 3)) == 5
    @test Agents.coord2vertex((1, 3, 2), (2, 3, 3)) == 11
    @test Agents.coord2vertex((1, 3, 3), (2, 3, 3)) == 17
    @test Agents.coord2vertex((2,), (2, 3, 3)) == 2
    @test Agents.coord2vertex(2, (2, 3, 3)) == 2

    @test Agents.vertex2coord(5, (2, 3, 3)) == (1, 3, 1)
    @test Agents.vertex2coord(7, (2, 3, 3)) == (1, 1, 2)
    @test Agents.vertex2coord(13, (2, 3, 3)) == (1, 1, 3)
    @test Agents.vertex2coord(13, (2, 3, 3)) == (1, 1, 3)
    @test Agents.vertex2coord(15, (2, 3, 3)) == (1, 2, 3)
    @test Agents.vertex2coord(18, (2, 3, 3)) == (2, 3, 3)
    @test Agents.vertex2coord(18, (2, 3, 3)) == (2, 3, 3)

    Random.seed!(648)
    model = ABM(Agent3, GridSpace((5,5)))
    agent = Agent3(1, (1,1), 5.5)
    add_agent!(agent, model)
    @test Agents.coord2vertex(agent, model) == 20
    @test Agents.coord2vertex((1, 3), model) == 11
    @test Agents.coord2vertex(15, model) == 15
    @test Agents.vertex2coord((2,3), model) == (2,3)
    @test_throws MethodError Agents.coord2vertex((1, 3, 7), model) == 11
    @test_throws ErrorException Agents.vertex2coord(3, GraphSpace(complete_digraph(5)))
end

@testset "nodes" begin
    Random.seed!(782)
    space = GridSpace((3, 3))
    model = ABM(Agent1, space)
    @test has_empty_nodes(model)
    for node in nodes(model)
        if rand() > 0.7
            add_agent!(node, model)
        end
    end
    @test pick_empty(model) == 4
    @test nodes(model) == 1:9
    @test nodes(model, by = :random) == [9, 1, 8, 6, 3, 5, 7, 2, 4]
    @test nodes(model, by = :population) == [1, 5, 6, 9, 2, 3, 4, 7, 8]
    @test length(get_node_contents(5, model)) > length(get_node_contents(7, model))
    @test_throws ErrorException nodes(model, by = :notreal)

    iter = NodeIterator(model)
    @test length(iter) == 9
    @test first(NodeIterator(model))[1] == model[1].pos
end

@testset "Neighbors" begin
    undirected = ABM(Agent5, GraphSpace(path_graph(5)))
    @test node_neighbors(3, undirected) == [2, 4]
    @test node_neighbors(1, undirected) == [2]
    @test node_neighbors(5, undirected) == [4]
    @test node_neighbors(3, undirected; neighbor_type = :out) ==
    node_neighbors(3, undirected; neighbor_type = :in) ==
    node_neighbors(3, undirected; neighbor_type = :all) ==
    node_neighbors(3, undirected; neighbor_type = :default)
    add_agent!(1, undirected, rand())
    add_agent!(2, undirected, rand())
    add_agent!(3, undirected, rand())
    # We expect id 2 to be included for a grid based search
    @test sort!(space_neighbors(2, undirected)) == [1, 2, 3]
    # But to be excluded if we are looking around it.
    @test sort!(space_neighbors(undirected[2], undirected)) == [1, 3]

    directed = ABM(Agent5, GraphSpace(path_digraph(5)))
    @test node_neighbors(3, directed) == [4]
    @test node_neighbors(1, directed) == [2]
    @test node_neighbors(5, directed) == []
    @test node_neighbors(3, directed; neighbor_type = :default) ==
          node_neighbors(3, directed)
    @test node_neighbors(3, directed; neighbor_type = :in) == [2]
    @test node_neighbors(3, directed; neighbor_type = :out) == [4]
    @test sort!(node_neighbors(3, directed; neighbor_type = :all)) == [2, 4]
    add_agent!(1, directed, rand())
    add_agent!(2, directed, rand())
    add_agent!(3, directed, rand())
    @test sort!(space_neighbors(2, directed)) == [2, 3]
    @test sort!(space_neighbors(2, directed; neighbor_type=:in)) == [1, 2]
    @test sort!(space_neighbors(2, directed; neighbor_type=:all)) == [1, 2, 3]
    @test space_neighbors(directed[2], directed) == [3]
    @test space_neighbors(directed[2], directed; neighbor_type=:in) == [1]
    @test sort!(space_neighbors(directed[2], directed; neighbor_type=:all)) == [1, 3]

    gridspace = ABM(Agent3, GridSpace((3, 3)))
    @test node_neighbors((2, 2), gridspace) == [(2, 1), (1, 2), (3, 2), (2, 3)]
    @test node_neighbors((1, 1), gridspace) == [(2, 1), (1, 2)]
    a = add_agent!((2, 2), gridspace, rand())
    add_agent!((3, 2), gridspace, rand())
    @test space_neighbors((1, 2), gridspace) == [1]
    @test sort!(space_neighbors((1, 2), gridspace, 2)) == [1, 2]
    @test_throws MethodError space_neighbors((1, 2), gridspace, 1.5)
    @test sort!(space_neighbors((2, 2), gridspace)) == [1, 2]
    @test space_neighbors(a, gridspace) == [2]

    continuousspace = ABM(Agent6, ContinuousSpace(2; extend = (1, 1)))
    a = add_agent!((0.5, 0.5), continuousspace, (0.2, 0.1), 0.01)
    b = add_agent!((0.6, 0.5), continuousspace, (0.1, -0.1), 0.01)
    @test_throws MethodError node_neighbors(1, continuousspace)
    @test space_neighbors(a, continuousspace, 0.05) == []
    @test space_neighbors(a, continuousspace, 0.1) == [2]
    @test sort!(space_neighbors((0.55, 0.5), continuousspace, 0.05)) == [1, 2]
    move_agent!(a, continuousspace)
    move_agent!(b, continuousspace)
    @test space_neighbors(a, continuousspace, 0.1) == []
    # Checks for type instability #208
    @test typeof(space_neighbors(a, continuousspace, 0.1)) <: Vector{Int}
    @test typeof(space_neighbors((0.55, 0.5), continuousspace, 0.05)) <: Vector{Int}
end


text:
using ScatterMatrixPlots,DataFrames
using Base.Test

# write your own tests here
@test 1 == 1

#データセット作成
x=1:10
y1=2x
y2=2x-1
label1=["a","a","a","b","b","b","b","c","c","c"]
label2=["1","2","3","1","2","3","1","2","3","1"]
df1=DataFrame(x=x,y1=y1,y2=y2,label1=label1)
df2=DataFrame(x=x,y1=y1,y2=y2,label1=label1,label2=label2)

#ファイル表示まで一括
ScatterMatrixPlot(df1)
ScatterMatrixPlot(df1, colorido="label1")
ScatterMatrixPlot(df2)
ScatterMatrixPlot(df2, colorido="label1")

#詳細設定
ScatterMatrixPlot(df1, colorido="label1",xwidth=10cm,ywidth=15cm)

max_stars_repo_name: wakakusa/ScatterMatrixPlots

max_issues_repo_name: wakakusa/ScatterMatrixPlots

max_forks_repo_name: wakakusa/ScatterMatrixPlots


text:
using ACE, JuLIP, ACEbonds, LinearAlgebra
using ACEbonds: bonds
using ACE: State

##

at = bulk(:Si, cubic=true) * 3
set_pbc!(at, false )
rattle!(at, 0.3)

rbcut = 3.0
rcut = 5.0
zcut = 7.0
env_filter = (r, z) -> (r < rcut) && (abs(z) < zcut )

it = bonds(at, rbcut, rbcut/2+zcut, env_filter)

(i, j, rr0, Js, Rs, Zs), state = iterate(it)

for (i1, j1, rr01, Js1, Rs1, Zs1) in bonds(at, rbcut, rbcut/2+zcut, env_filter)
   if i1 == length(at) ÷ 2
      i = i1; j = j1; rr0 = rr01; Js = Js1; Rs = Rs1; Zs = Zs1
      break
   end
end

# now we can transform it to an ACE.jl compatible environment
env = ACEbonds.eucl2cyl(rr0, at.Z[i], at.Z[j], Rs, Zs)

# ## Manual debugging of the bond environments

# using GLMakie, Makie

# function plot_bondenv(i, j, Js, at)
#    X = at.X
#    xb = [ X[i][1], X[j][1] ]
#    yb = [ X[i][2], X[j][2] ]
#    zb = [ X[i][3], X[j][3] ]
#    colb = [:red, :red]

#    Xenv = X[Js[2:end]]
#    xe = [ x[1] for x in Xenv ]
#    ye = [ x[2] for x in Xenv ]
#    ze = [ x[3] for x in Xenv ]
#    cole = [:blue for _ in Xenv]

#    Xrest = X[setdiff(1:length(X), union([i,j], Js))]
#    xr = [ x[1] for x in Xrest ]
#    yr = [ x[2] for x in Xrest ]
#    zr = [ x[3] for x in Xrest ]
#    colr = [(:green, 0.3) for _ in Xrest]

#    scene = meshscatter([xb; xe; xr], [yb; ye; yr], [zb; ze; zr],
#                        color = [colb; cole; colr],
#                        markersize=0.7)
#    lines!(xb, yb, zb, color=:red, linewidth=10 )
#    return scene
# end

# scene = plot_bondenv(i, j, Js, at);
# window = display(scene)
# # close(window) to close it.


text:
# Export
export mdp_fs
export satellite_plan_mdp_fs

function mdp_reachable_states(problem::SatPlanningProblem, state::SatMDPState, action::Opportunity)
    return SatMDPState[baseline_step(problem, state, action)]
end

function mdp_depth_first_search(problem::SatPlanningProblem, state::SatMDPState, depth::Integer)

    # Early exit
    if depth == 0
        return Noop(t_start=state.time), 0.0
    end

    astar, vstar = Done(t_start=state.time), -Inf

    for a in problem.lt_feasible_actions[(state.last_cdo_action.id, state.last_action.id)]
        v = POMDPs.reward(problem, state, a)

        for sp in mdp_reachable_states(problem, state, a)
            ap, vp = mdp_depth_first_search(problem, sp, depth-1)

            # Strict Markov Update
            # v = v + problem.solve_discount*vp

            # Get distance action is in the future
            t_diff = abs(ap.t_start - state.time)

            # Semi-markov update
            v = v + vp*problem.solve_discount^(t_diff)
        end

        if v > vstar
            astar, vstar = a, v
        end
    end

    return astar, vstar
end

function mdp_fs(problem::SatPlanningProblem, state::SatMDPState)

    # println("Foward Search Step")
    action, value = mdp_depth_first_search(problem, state, problem.solve_depth)

    if typeof(action) == Collect
        if action.location.id in state.request_ids
            # println("Planned duplicate collect: $(action) - $(action.location)")
        end
    end

    # if typeof(action) == Contact
    #     println("Taking Contact: $action")
    # end

    if typeof(action) != Done && action.t_start == state.time
        # throw(ErrorException("Took action that isn't advancing time...\n $state - $action"))
        @debug "Took action that isn't advancing time...\n $state - $action"
        return state, Done(t_start=state.time), 0
    end

    # Step to next state with selected action
    state = baseline_step(problem, state, action)

    return state, action, value
end

function satellite_plan_mdp_fs(problem::SatPlanningProblem)

    # Sort opportunities
    sort!(problem.opportunities, by = x -> x.t_start)

    # Set initial state
    init_opp = problem.opportunities[1]
    state = SatMDPState(time=init_opp.t_start, last_action=init_opp, last_cdo_action=init_opp)

    states = SatMDPState[state]
    plan = Opportunity[state.last_action]
    total_reward = 0.0

    while true
        state, action, reward = mdp_fs(problem, state)

        # Break from search if terminal state
        if typeof(action) == Done
            break
        end

        # Add state and action to plan
        push!(states, state)
        push!(plan, action)
        total_reward += reward
    end

    return states, plan, total_reward
end


text: potential_energy(s::MosiSystem, model::MosiModel) = potential_energy_function(model, positions(s))


text:
# This file contains helper functions for testing the various
# GraphIO formats

using LightGraphs

graphs = Dict{String,Graph}(
    "graph1"    => CompleteGraph(5),
    "graph2"    => PathGraph(6),
    "graph3"    =>WheelGraph(4)
    )
digraphs = Dict{String,DiGraph}(
    "digraph1"   => CompleteDiGraph(5),
    "digraph2"   => PathDiGraph(6),
    "digraph3"   => WheelDiGraph(4)
)
allgraphs = merge(graphs, digraphs)

function gettempname()
    (f, fio) = mktemp()
    close(fio)
    return f
end

function read_test(format::LightGraphs.AbstractGraphFormat, g::LightGraphs.AbstractGraph, gname::String="g",
    fname::AbstractString=""; testfail=false)
    @test loadgraph(fname, gname, format) == g
    if testfail
        @test_throws Union{ArgumentError, ErrorException} loadgraph(fname, "badgraphXXX", format)
    end
    @test loadgraphs(fname, format)[gname] == g
end

function read_test_mult(format::LightGraphs.AbstractGraphFormat, d::Dict{String,G}, fname::AbstractString="") where G<: AbstractGraph
    rd = loadgraphs(fname, format)
    @test rd == d

end

function write_test(format::LightGraphs.AbstractGraphFormat, g::LightGraphs.AbstractGraph, gname::String="g",
    fname::AbstractString=gettempname(); remove=true, silent=false)
    @test savegraph(fname, g, gname, format) == 1
    if remove
        rm(fname)
    elseif !silent
        info("graphio/write_test: Left temporary file at: $fname")
    end
end

function write_test(format::LightGraphs.AbstractGraphFormat, d::Dict{String,G},
    fname::AbstractString=gettempname(); remove=true, silent=false) where G <: LightGraphs.AbstractGraph
    @test savegraph(fname, d, format) == length(d)
    if remove
        rm(fname)
    elseif !silent
        info("graphio/write_test: Left temporary file at: $fname")
    end
end

function readback_test(format::LightGraphs.AbstractGraphFormat, g::LightGraphs.AbstractGraph, gname="graph",
                        fname=gettempname(); remove=true, testfail=false)
    @test savegraph(fname, g, format) == 1
    @test loadgraphs(fname, format)[gname] == g
    @test loadgraph(fname, gname, format) == g
    if testfail
        @test_throws Union{ArgumentError, ErrorException} loadgraph(fname, "badgraphXXX", format)
    end
    if remove
        rm(fname)
    else
        info("graphio/readback_test: Left temporary file at: $fname")
    end
end


text:
"""
Stores information associated with a quadrature rule. zeta contains the positions of quadrature points.
"""
struct TriTetQuadRule{N_zeta, N_point}
    zeta::NTuple{N_point, SVector{N_zeta, Float64}}
    w::NTuple{N_point, Float64}
    the_degree::Int64
    function TriTetQuadRule{N_zeta, N_point}(a, b, c)  where {N_zeta, N_point}
        return new{N_zeta, N_point}(a, b, c)
    end
end

"""
Outputs the XiaoGimbutas triangle quadrature rules from the Python package quadpy.
"""
function getTriQuadRule(n_rule::Int64)
    if n_rule == 1
        n_quad_points = 1
        p = Vector{SVector{3, Float64}}(undef, n_quad_points)
        p[1] = SVector{3, Float64}(0.33333333333333331483, 0.33333333333333331483, 0.33333333333333331483)
        w = zeros(n_quad_points)
        w[1] = 1.00000000000000000000
        the_degree = 1
        return TriTetQuadRule{3, n_quad_points}(Tuple(p), Tuple(w), the_degree)
    elseif n_rule == 2
        n_quad_points = 3
        p = Vector{SVector{3, Float64}}(undef, n_quad_points)
        p[1] = SVector{3, Float64}(0.16666666666666674068, 0.66666666666666651864, 0.16666666666666674068)
        p[2] = SVector{3, Float64}(0.66666666666666651864, 0.16666666666666674068, 0.16666666666666674068)
        p[3] = SVector{3, Float64}(0.16666666666666674068, 0.16666666666666674068, 0.66666666666666651864)
        w = zeros(n_quad_points)
        w[1] = 0.33333333333333331483
        w[2] = 0.33333333333333331483
        w[3] = 0.33333333333333331483
        the_degree = 2
        return TriTetQuadRule{3, n_quad_points}(Tuple(p), Tuple(w), the_degree)
    elseif n_rule == 3
        n_quad_points = 6
        p = Vector{SVector{3, Float64}}(undef, n_quad_points)
        p[1] = SVector{3, Float64}(0.44594849091596489021, 0.10810301816807021957, 0.44594849091596489021)
        p[2] = SVector{3, Float64}(0.09157621350977085406, 0.81684757298045829188, 0.09157621350977085406)
        p[3] = SVector{3, Float64}(0.10810301816807021957, 0.44594849091596489021, 0.44594849091596489021)
        p[4] = SVector{3, Float64}(0.81684757298045829188, 0.09157621350977085406, 0.09157621350977085406)
        p[5] = SVector{3, Float64}(0.44594849091596489021, 0.44594849091596489021, 0.10810301816807021957)
        p[6] = SVector{3, Float64}(0.09157621350977085406, 0.09157621350977085406, 0.81684757298045829188)
        w = zeros(n_quad_points)
        w[1] = 0.22338158967801147181
        w[2] = 0.10995174365532188465
        w[3] = 0.22338158967801147181
        w[4] = 0.10995174365532188465
        w[5] = 0.22338158967801147181
        w[6] = 0.10995174365532188465
        the_degree = 3
        return TriTetQuadRule{3, n_quad_points}(Tuple(p), Tuple(w), the_degree)
    elseif n_rule == 4
        n_quad_points = 6
        p = Vector{SVector{3, Float64}}(undef, n_quad_points)
        p[1] = SVector{3, Float64}(0.44594849091596489021, 0.10810301816807021957, 0.44594849091596489021)
        p[2] = SVector{3, Float64}(0.09157621350977085406, 0.81684757298045829188, 0.09157621350977085406)
        p[3] = SVector{3, Float64}(0.10810301816807021957, 0.44594849091596489021, 0.44594849091596489021)
        p[4] = SVector{3, Float64}(0.81684757298045829188, 0.09157621350977085406, 0.09157621350977085406)
        p[5] = SVector{3, Float64}(0.44594849091596489021, 0.44594849091596489021, 0.10810301816807021957)
        p[6] = SVector{3, Float64}(0.09157621350977085406, 0.09157621350977085406, 0.81684757298045829188)
        w = zeros(n_quad_points)
        w[1] = 0.22338158967801147181
        w[2] = 0.10995174365532188465
        w[3] = 0.22338158967801147181
        w[4] = 0.10995174365532188465
        w[5] = 0.22338158967801147181
        w[6] = 0.10995174365532188465
        the_degree = 4
        return TriTetQuadRule{3, n_quad_points}(Tuple(p), Tuple(w), the_degree)
    elseif n_rule == 5
        n_quad_points = 7
        p = Vector{SVector{3, Float64}}(undef, n_quad_points)
        p[1] = SVector{3, Float64}(0.33333333333333331483, 0.33333333333333331483, 0.33333333333333337034)
        p[2] = SVector{3, Float64}(0.10128650732345639840, 0.79742698535308720320, 0.10128650732345634289)
        p[3] = SVector{3, Float64}(0.47014206410511505396, 0.05971587178976989208, 0.47014206410511505396)
        p[4] = SVector{3, Float64}(0.79742698535308720320, 0.10128650732345639840, 0.10128650732345634289)
        p[5] = SVector{3, Float64}(0.05971587178976989208, 0.47014206410511505396, 0.47014206410511505396)
        p[6] = SVector{3, Float64}(0.10128650732345639840, 0.10128650732345639840, 0.79742698535308720320)
        p[7] = SVector{3, Float64}(0.47014206410511505396, 0.47014206410511505396, 0.05971587178976989208)
        w = zeros(n_quad_points)
        w[1] = 0.22500000000000000555
        w[2] = 0.12593918054482713953
        w[3] = 0.13239415278850619195
        w[4] = 0.12593918054482713953
        w[5] = 0.13239415278850619195
        w[6] = 0.12593918054482713953
        w[7] = 0.13239415278850619195
        the_degree = 5
        return TriTetQuadRule{3, n_quad_points}(Tuple(p), Tuple(w), the_degree)
    else
        error()
    end
end

"""
Outputs the Yu tetrahedron quadrature rules from the Python package quadpy.
"""
function getTetQuadRule(n_rule::Int64)
    if n_rule == 1
        n_quad_points = 4
        p = Vector{SVector{4, Float64}}(undef, n_quad_points)
        p[1] = SVector{4, Float64}(0.13819660112501500060, 0.13819660112501500060, 0.58541019662495497045, 0.13819660112501508387)
        p[2] = SVector{4, Float64}(0.13819660112501500060, 0.58541019662495497045, 0.13819660112501500060, 0.13819660112501508387)
        p[3] = SVector{4, Float64}(0.58541019662495497045, 0.13819660112501500060, 0.13819660112501500060, 0.13819660112501508387)
        p[4] = SVector{4, Float64}(0.13819660112501500060, 0.13819660112501500060, 0.13819660112501500060, 0.58541019662495497045)
        w = zeros(n_quad_points)
        w[1] = 0.25000000000000000000
        w[2] = 0.25000000000000000000
        w[3] = 0.25000000000000000000
        w[4] = 0.25000000000000000000
        the_degree = 2
        return TriTetQuadRule{4, n_quad_points}(Tuple(p), Tuple(w), the_degree)
    elseif n_rule == 2
        n_quad_points = 5
        p = Vector{SVector{4, Float64}}(undef, n_quad_points)
        p[1] = SVector{4, Float64}(0.25000000000000000000, 0.25000000000000000000, 0.25000000000000000000, 0.25000000000000000000)
        p[2] = SVector{4, Float64}(0.16666666666666665741, 0.16666666666666665741, 0.50000000000000000000, 0.16666666666666674068)
        p[3] = SVector{4, Float64}(0.16666666666666665741, 0.50000000000000000000, 0.16666666666666665741, 0.16666666666666674068)
        p[4] = SVector{4, Float64}(0.50000000000000000000, 0.16666666666666665741, 0.16666666666666665741, 0.16666666666666674068)
        p[5] = SVector{4, Float64}(0.16666666666666665741, 0.16666666666666665741, 0.16666666666666665741, 0.50000000000000000000)
        w = zeros(n_quad_points)
        w[1] = -0.80000000000000004441
        w[2] = 0.45000000000000001110
        w[3] = 0.45000000000000001110
        w[4] = 0.45000000000000001110
        w[5] = 0.45000000000000001110
        the_degree = 3
        return TriTetQuadRule{4, n_quad_points}(Tuple(p), Tuple(w), the_degree)
    elseif n_rule == 3
        n_quad_points = 16
        p = Vector{SVector{4, Float64}}(undef, n_quad_points)
        p[1] = SVector{4, Float64}(0.07611903264425430582, 0.07611903264425430582, 0.77164290206723706866, 0.07611903264425434745)
        p[2] = SVector{4, Float64}(0.07611903264425430582, 0.77164290206723706866, 0.07611903264425430582, 0.07611903264425423643)
        p[3] = SVector{4, Float64}(0.77164290206723706866, 0.07611903264425430582, 0.07611903264425430582, 0.07611903264425423643)
        p[4] = SVector{4, Float64}(0.07611903264425430582, 0.07611903264425430582, 0.07611903264425430582, 0.77164290206723706866)
        p[5] = SVector{4, Float64}(0.40423391346726439410, 0.11970052779780189500, 0.07183164526766931679, 0.40423391346726444961)
        p[6] = SVector{4, Float64}(0.11970052779780189500, 0.40423391346726439410, 0.07183164526766931679, 0.40423391346726444961)
        p[7] = SVector{4, Float64}(0.40423391346726439410, 0.40423391346726439410, 0.07183164526766931679, 0.11970052779780193664)
        p[8] = SVector{4, Float64}(0.11970052779780189500, 0.07183164526766931679, 0.40423391346726439410, 0.40423391346726433859)
        p[9] = SVector{4, Float64}(0.40423391346726439410, 0.07183164526766931679, 0.40423391346726439410, 0.11970052779780182561)
        p[10] = SVector{4, Float64}(0.07183164526766931679, 0.40423391346726439410, 0.40423391346726439410, 0.11970052779780182561)
        p[11] = SVector{4, Float64}(0.40423391346726439410, 0.07183164526766931679, 0.11970052779780189500, 0.40423391346726433859)
        p[12] = SVector{4, Float64}(0.07183164526766931679, 0.40423391346726439410, 0.11970052779780189500, 0.40423391346726433859)
        p[13] = SVector{4, Float64}(0.40423391346726439410, 0.40423391346726439410, 0.11970052779780189500, 0.07183164526766927516)
        p[14] = SVector{4, Float64}(0.07183164526766931679, 0.11970052779780189500, 0.40423391346726439410, 0.40423391346726433859)
        p[15] = SVector{4, Float64}(0.40423391346726439410, 0.11970052779780189500, 0.40423391346726439410, 0.07183164526766927516)
        p[16] = SVector{4, Float64}(0.11970052779780189500, 0.40423391346726439410, 0.40423391346726439410, 0.07183164526766927516)
        w = zeros(n_quad_points)
        w[1] = 0.05037379410012281927
        w[2] = 0.05037379410012281927
        w[3] = 0.05037379410012281927
        w[4] = 0.05037379410012281927
        w[5] = 0.06654206863329238664
        w[6] = 0.06654206863329238664
        w[7] = 0.06654206863329238664
        w[8] = 0.06654206863329238664
        w[9] = 0.06654206863329238664
        w[10] = 0.06654206863329238664
        w[11] = 0.06654206863329238664
        w[12] = 0.06654206863329238664
        w[13] = 0.06654206863329238664
        w[14] = 0.06654206863329238664
        w[15] = 0.06654206863329238664
        w[16] = 0.06654206863329238664
        the_degree = 4
        return TriTetQuadRule{4, n_quad_points}(Tuple(p), Tuple(w), the_degree)
    elseif n_rule == 4
        n_quad_points = 17
        p = Vector{SVector{4, Float64}}(undef, n_quad_points)
        p[1] = SVector{4, Float64}(0.25000000000000000000, 0.25000000000000000000, 0.25000000000000000000, 0.25000000000000000000)
        p[2] = SVector{4, Float64}(0.08945436401412733274, 0.08945436401412733274, 0.73163690795761793240, 0.08945436401412742988)
        p[3] = SVector{4, Float64}(0.08945436401412733274, 0.73163690795761793240, 0.08945436401412733274, 0.08945436401412742988)
        p[4] = SVector{4, Float64}(0.73163690795761793240, 0.08945436401412733274, 0.08945436401412733274, 0.08945436401412742988)
        p[5] = SVector{4, Float64}(0.08945436401412733274, 0.08945436401412733274, 0.08945436401412733274, 0.73163690795761793240)
        p[6] = SVector{4, Float64}(0.42143943106625220141, 0.13258109993846570185, 0.02454003792902989534, 0.42143943106625214590)
        p[7] = SVector{4, Float64}(0.13258109993846570185, 0.42143943106625220141, 0.02454003792902989534, 0.42143943106625214590)
        p[8] = SVector{4, Float64}(0.42143943106625220141, 0.42143943106625220141, 0.02454003792902989534, 0.13258109993846567409)
        p[9] = SVector{4, Float64}(0.13258109993846570185, 0.02454003792902989534, 0.42143943106625220141, 0.42143943106625214590)
        p[10] = SVector{4, Float64}(0.42143943106625220141, 0.02454003792902989534, 0.42143943106625220141, 0.13258109993846578512)
        p[11] = SVector{4, Float64}(0.02454003792902989534, 0.42143943106625220141, 0.42143943106625220141, 0.13258109993846578512)
        p[12] = SVector{4, Float64}(0.42143943106625220141, 0.02454003792902989534, 0.13258109993846570185, 0.42143943106625225692)
        p[13] = SVector{4, Float64}(0.02454003792902989534, 0.42143943106625220141, 0.13258109993846570185, 0.42143943106625225692)
        p[14] = SVector{4, Float64}(0.42143943106625220141, 0.42143943106625220141, 0.13258109993846570185, 0.02454003792902992309)
        p[15] = SVector{4, Float64}(0.02454003792902989534, 0.13258109993846570185, 0.42143943106625220141, 0.42143943106625214590)
        p[16] = SVector{4, Float64}(0.42143943106625220141, 0.13258109993846570185, 0.42143943106625220141, 0.02454003792902992309)
        p[17] = SVector{4, Float64}(0.13258109993846570185, 0.42143943106625220141, 0.42143943106625220141, 0.02454003792902992309)
        w = zeros(n_quad_points)
        w[1] = 0.18841855673654109293
        w[2] = 0.06703858372604275251
        w[3] = 0.06703858372604275251
        w[4] = 0.06703858372604275251
        w[5] = 0.06703858372604275251
        w[6] = 0.04528559236327398679
        w[7] = 0.04528559236327398679
        w[8] = 0.04528559236327398679
        w[9] = 0.04528559236327398679
        w[10] = 0.04528559236327398679
        w[11] = 0.04528559236327398679
        w[12] = 0.04528559236327398679
        w[13] = 0.04528559236327398679
        w[14] = 0.04528559236327398679
        w[15] = 0.04528559236327398679
        w[16] = 0.04528559236327398679
        w[17] = 0.04528559236327398679
        the_degree = 5
        return TriTetQuadRule{4, n_quad_points}(Tuple(p), Tuple(w), the_degree)
    elseif n_rule == 5
        n_quad_points = 29
        p = Vector{SVector{4, Float64}}(undef, n_quad_points)
        p[1] = SVector{4, Float64}(0.25000000000000000000, 0.25000000000000000000, 0.25000000000000000000, 0.25000000000000000000)
        p[2] = SVector{4, Float64}(0.05742691731735682231, 0.05742691731735682231, 0.82771924804792951225, 0.05742691731735682925)
        p[3] = SVector{4, Float64}(0.05742691731735682231, 0.82771924804792951225, 0.05742691731735682231, 0.05742691731735682925)
        p[4] = SVector{4, Float64}(0.82771924804792951225, 0.05742691731735682231, 0.05742691731735682231, 0.05742691731735682925)
        p[5] = SVector{4, Float64}(0.05742691731735682231, 0.05742691731735682231, 0.05742691731735682231, 0.82771924804792951225)
        p[6] = SVector{4, Float64}(0.23129854365191471199, 0.05135188412556340731, 0.48605102857060716870, 0.23129854365191471199)
        p[7] = SVector{4, Float64}(0.04756909881472289742, 0.29675381296902597628, 0.60810798940152821501, 0.04756909881472293211)
        p[8] = SVector{4, Float64}(0.05135188412556340731, 0.23129854365191471199, 0.48605102857060716870, 0.23129854365191471199)
        p[9] = SVector{4, Float64}(0.29675381296902597628, 0.04756909881472289742, 0.60810798940152821501, 0.04756909881472293211)
        p[10] = SVector{4, Float64}(0.23129854365191471199, 0.23129854365191471199, 0.48605102857060716870, 0.05135188412556335180)
        p[11] = SVector{4, Float64}(0.04756909881472289742, 0.04756909881472289742, 0.60810798940152821501, 0.29675381296902603179)
        p[12] = SVector{4, Float64}(0.05135188412556340731, 0.48605102857060716870, 0.23129854365191471199, 0.23129854365191471199)
        p[13] = SVector{4, Float64}(0.29675381296902597628, 0.60810798940152821501, 0.04756909881472289742, 0.04756909881472293211)
        p[14] = SVector{4, Float64}(0.23129854365191471199, 0.48605102857060716870, 0.23129854365191471199, 0.05135188412556346282)
        p[15] = SVector{4, Float64}(0.04756909881472289742, 0.60810798940152821501, 0.04756909881472289742, 0.29675381296902592076)
        p[16] = SVector{4, Float64}(0.48605102857060716870, 0.23129854365191471199, 0.23129854365191471199, 0.05135188412556346282)
        p[17] = SVector{4, Float64}(0.60810798940152821501, 0.04756909881472289742, 0.04756909881472289742, 0.29675381296902592076)
        p[18] = SVector{4, Float64}(0.23129854365191471199, 0.48605102857060716870, 0.05135188412556340731, 0.23129854365191482302)
        p[19] = SVector{4, Float64}(0.04756909881472289742, 0.60810798940152821501, 0.29675381296902597628, 0.04756909881472282109)
        p[20] = SVector{4, Float64}(0.48605102857060716870, 0.23129854365191471199, 0.05135188412556340731, 0.23129854365191482302)
        p[21] = SVector{4, Float64}(0.60810798940152821501, 0.04756909881472289742, 0.29675381296902597628, 0.04756909881472282109)
        p[22] = SVector{4, Float64}(0.23129854365191471199, 0.23129854365191471199, 0.05135188412556340731, 0.48605102857060722421)
        p[23] = SVector{4, Float64}(0.04756909881472289742, 0.04756909881472289742, 0.29675381296902597628, 0.60810798940152821501)
        p[24] = SVector{4, Float64}(0.48605102857060716870, 0.05135188412556340731, 0.23129854365191471199, 0.23129854365191471199)
        p[25] = SVector{4, Float64}(0.60810798940152821501, 0.29675381296902597628, 0.04756909881472289742, 0.04756909881472293211)
        p[26] = SVector{4, Float64}(0.23129854365191471199, 0.05135188412556340731, 0.23129854365191471199, 0.48605102857060722421)
        p[27] = SVector{4, Float64}(0.04756909881472289742, 0.29675381296902597628, 0.04756909881472289742, 0.60810798940152821501)
        p[28] = SVector{4, Float64}(0.05135188412556340731, 0.23129854365191471199, 0.23129854365191471199, 0.48605102857060722421)
        p[29] = SVector{4, Float64}(0.29675381296902597628, 0.04756909881472289742, 0.04756909881472289742, 0.60810798940152821501)
        w = zeros(n_quad_points)
        w[1] = 0.09040129046014749870
        w[2] = 0.01911983427899123858
        w[3] = 0.01911983427899123858
        w[4] = 0.01911983427899123858
        w[5] = 0.01911983427899123858
        w[6] = 0.04361493840666567923
        w[7] = 0.02581167596199161057
        w[8] = 0.04361493840666567923
        w[9] = 0.02581167596199161057
        w[10] = 0.04361493840666567923
        w[11] = 0.02581167596199161057
        w[12] = 0.04361493840666567923
        w[13] = 0.02581167596199161057
        w[14] = 0.04361493840666567923
        w[15] = 0.02581167596199161057
        w[16] = 0.04361493840666567923
        w[17] = 0.02581167596199161057
        w[18] = 0.04361493840666567923
        w[19] = 0.02581167596199161057
        w[20] = 0.04361493840666567923
        w[21] = 0.02581167596199161057
        w[22] = 0.04361493840666567923
        w[23] = 0.02581167596199161057
        w[24] = 0.04361493840666567923
        w[25] = 0.02581167596199161057
        w[26] = 0.04361493840666567923
        w[27] = 0.02581167596199161057
        w[28] = 0.04361493840666567923
        w[29] = 0.02581167596199161057
        the_degree = 6
        return TriTetQuadRule{4, n_quad_points}(Tuple(p), Tuple(w), the_degree)
    else
        error()
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

# TODO PR: Remove this line; this is a dummy change to trigger CI builds.

name = "Fontconfig"
version = v"2.13.93"

# Collection of sources required to build FriBidi
sources = [
    ArchiveSource("https://www.freedesktop.org/software/fontconfig/release/fontconfig-$(version).tar.xz",
                  "ea968631eadc5739bc7c8856cef5c77da812d1f67b763f5e51b57b8026c1a0a0"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/fontconfig-*/

# Ensure that `${includedir}` is..... included
export CPPFLAGS="-I${includedir}"

FLAGS=()
if [[ "${target}" == *-linux-* ]] || [[ "${target}" == *-freebsd* ]]; then
    FLAGS+=(--with-add-fonts="/usr/local/share/fonts")
    if [[ "${target}" == *-linux-* ]]; then
        FLAGS+=(--with-cache-dir="/var/cache/fontconfig")
    elif [[ "${target}" == *-freebsd* ]]; then
        FLAGS+=(--with-cache-dir="/var/db/fontconfig")
    fi
elif [[ "${target}" == *-apple-* ]]; then
    FLAGS+=(--with-add-fonts="/System/Library/Fonts,/Library/Fonts,~/Library/Fonts,/System/Library/Assets/com_apple_MobileAsset_Font4,/System/Library/Assets/com_apple_MobileAsset_Font5")
fi

# Apply MinGW patches: https://github.com/msys2/MINGW-packages/tree/33f847297fe429d145cd9d72cb1fbbc574431cc5/mingw-w64-fontconfig
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/0001-fix-config-linking.all.patch"
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/0002-fix-mkdir.mingw.patch"
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/0004-fix-mkdtemp.mingw.patch"
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/0005-fix-setenv.mingw.patch"
autoreconf
./configure --prefix=$prefix --build=${MACHTYPE} --host=$target --disable-docs "${FLAGS[@]}"

# Disable tests
sed -i 's,all-am: Makefile $(PROGRAMS),all-am:,' test/Makefile

make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libfontconfig", :libfontconfig),
    ExecutableProduct("fc-cache", :fc_cache),
    ExecutableProduct("fc-cat", :fc_cat),
    ExecutableProduct("fc-conflist", :fc_conflist),
    ExecutableProduct("fc-list", :fc_list),
    ExecutableProduct("fc-match", :fc_match),
    ExecutableProduct("fc-pattern", :fc_pattern),
    ExecutableProduct("fc-query", :fc_query),
    ExecutableProduct("fc-scan", :fc_scan),
    ExecutableProduct("fc-validate", :fc_validate),
    FileProduct("etc/fonts/fonts.conf", :fonts_conf),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    HostBuildDependency("gperf_jll"),
    Dependency("FreeType2_jll"),
    Dependency("Bzip2_jll", v"1.0.8"; compat="1.0.8"),
    Dependency("Zlib_jll"),
    Dependency("Libuuid_jll"),
    Dependency("Expat_jll", v"2.2.10"; compat="~2.2.7"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat="1.6",
               init_block = """
ENV["FONTCONFIG_FILE"] = get(ENV, "FONTCONFIG_FILE", fonts_conf)
    ENV["FONTCONFIG_PATH"] = get(ENV, "FONTCONFIG_PATH", dirname(ENV["FONTCONFIG_FILE"]))
""")

max_stars_repo_name: luoq/Yggdrasil

max_issues_repo_name: luoq/Yggdrasil

max_forks_repo_name: luoq/Yggdrasil


text:
module QuantumESPRESSO

include("Inputs.jl")
include("Outputs.jl")
include("Commands.jl")

end


text:
__precompile__()
module ATA

using Printf
using CSV
using DelimitedFiles
using DataFrames
using Distributions
using Distributed
using Random
using LinearAlgebra
using JLD2
using StatsBase
using Requires
using Psychometrics
using FileIO

include("structs/structs.jl")
include("log/log.jl")
include("utils.jl")
include("build/build.jl")
include("opt.jl")
include("print/print.jl")
include("compact_ata.jl")

function __init__()
    @require JuMP = "4076af6c-e467-56ae-b986-b466b2749572" include("jump/jump.jl")

    @require Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76" begin
        function add_Cbc!(m)
            JuMP.set_optimizer(m, Cbc.Optimizer)
        end
    end
    @require Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b" begin
        function add_Gurobi!(m)
            JuMP.set_optimizer(m, Gurobi.Optimizer)
        end
    end
    @require CPLEX = "a076750e-1247-5638-91d2-ce28b192dca0" begin
        function add_CPLEX!(m)
            JuMP.set_optimizer(m, CPLEX.Optimizer)
        end
    end
    @require GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6" begin
        function add_GLPK!(m)
            JuMP.set_optimizer(m, GLPK.Optimizer)
        end
    end
    @require Xpress = "9e70acf3-d6c9-5be6-b5bd-4e2c73e3e054" begin
        function add_Xpress!(m)
            JuMP.set_optimizer(m, Gurobi.Optimizer)
        end
    end
    @require KNITRO = "67920dd8-b58e-52a8-8622-53c4cffbe346" begin
        function add_KNITRO!(m)
            JuMP.set_optimizer(m, KNITRO.Optimizer)
        end
    end
    @require Juniper = "2ddba703-00a4-53a7-87a5-e8b9971dde84" begin
        function add_Juniper!(m)
            JuMP.set_optimizer(m, Juniper.Optimizer)
        end
    end
    @require SCIP = "bdc5fa1a-e2a0-4ec3-bc77-c59bd234be9e" begin
        function add_SCIP!(m)
            JuMP.set_optimizer(m, SCIP.Optimizer)
        end
    end
    @require ATAPlot = "372623d9-1dd6-4d20-a513-3f20705132c0" begin
        include("plot.jl")
    end
    @require ATADash = "236b7dbe-4167-40cb-a459-bcf8ce4b2cbd" begin
        function ata_app()
            ATADash.ata_app()
        end
    end
end

export mean,
    start_ata,
    add_bank!,
    add_constraints!,
    add_overlap!,
    add_obj_fun!,
    assemble!,
    print_results,
    plot_results,
    print_infos,
    print_last_info,
    load_design!,
    item_info,
    item_char,
    resp_gen,
    student_likelihood,
    fs_to_items,#_mycopy,  fast_sort!, _myqle#, optimize
    ata_app,
    compact_ata

ATA

end


text:
struct GitRepo
    path::String
end

max_stars_repo_name: simonbyrne/GitX

max_issues_repo_name: simonbyrne/GitX

max_forks_repo_name: simonbyrne/GitX


text:
# Used "Conjugate Bayesian analysis of the Gaussian distribution" by Murphy as
# a reference.  Note that there were some typos in that document so the code
# here may not correspond exactly.

immutable NormalInverseWishart <: Distribution
    dim::Int
    zeromean::Bool
    mu::Vector{Float64}
    kappa::Float64              # This scales precision (inverse covariance)
    Lamchol::Cholesky{Float64}  # Covariance matrix (well, sqrt of one)
    nu::Float64

    function NormalInverseWishart(mu::Vector{Float64}, kappa::Real,
                                  Lamchol::Cholesky{Float64}, nu::Real)
        # Probably should put some error checking in here
        d = length(mu)
        zmean::Bool = true
        for i = 1:d
            if mu[i] != 0.
                zmean = false
                break
            end
        end
        new(d, zmean, mu, float64(kappa), Lamchol, float64(nu))
    end
end

function NormalInverseWishart(mu::Vector{Float64}, kappa::Real,
                              Lambda::Matrix{Float64}, nu::Real)
    NormalInverseWishart(mu, kappa, cholfact(Lambda), nu)

end

function insupport(::Type{NormalInverseWishart}, x::Vector{Float64}, Sig::Matrix{Float64})
    return (all(isfinite(x)) &&
           size(Sig, 1) == size(Sig, 2) &&
           isApproxSymmmetric(Sig) &&
           size(Sig, 1) == length(x) &&
           hasCholesky(Sig))
end

pdf(niw::NormalInverseWishart, x::Vector{Float64}, Sig::Matrix{Float64}) =
        exp(logpdf(niw, x, Sig))

function logpdf(niw::NormalInverseWishart, x::Vector{Float64}, Sig::Matrix{Float64})
    if !insupport(NormalInverseWishart, x, Sig)
        return -Inf
    else
        p = size(x, 1)

        nu = niw.nu
        kappa = niw.kappa
        mu = niw.mu
        Lamchol = niw.Lamchol
        hnu = 0.5 * nu
        hp = 0.5 * p

        # Normalization
        logp::Float64 = hnu * logdet(Lamchol)
        logp -= hnu * p * log(2.)
        logp -= lpgamma(p, hnu)
        logp -= hp * (log(2.*pi) - log(kappa))

        # Inverse-Wishart
        logp -= (hnu + hp + 1.) * logdet(Sig)
        logp -= 0.5 * trace(Sig \ (Lamchol[:U]' * Lamchol[:U]))

        # Normal
        z = niw.zeromean ? x : x - mu
        logp -= 0.5 * kappa * invquad(PDMat(Sig), z)

        return logp

    end
end

function rand(niw::NormalInverseWishart)
    Sig = rand(InverseWishart(niw.nu, niw.Lamchol))
    mu = rand(MvNormal(niw.mu, Sig ./ niw.kappa))
    return (mu, Sig)
end


text:
"""
    GAN{P<:CMD,G<:ACD,D<:ACD}([zlength::Int,p::P], g::G, d::D)

The Generative Adversarial Network.

# Arguments
* `zlength`: Length of latent vector
* `p`: Prior
* `g`: Generator
* `d`: Discriminator

# Example
Create a GAN with standard normal prior with:
```julia-repl
julia> gen = CMeanGaussian{DiagVar}(Dense(2,4),NoGradArray(ones(Float32,4)))
CMeanGaussian{DiagVar}(mapping=Dense(2, 4), σ2=4-element Array{Float32,1}

julia> disc = CMeanGaussian{DiagVar}(Dense(4,1,σ),NoGradArray(ones(Float32,1)))
CMeanGaussian{DiagVar}(mapping=Dense(4, 1, σ), σ2=1-element Array{Float32,1}

julia> gan = GAN(4, gen, disc)
GAN:
 prior   = (Gaussian(μ=4-element Array{Float32,1}, σ2=4-element Array{Float32...)
 generator = (CMeanGaussian{DiagVar}(mapping=Dense(2, 4), σ2=4-element Array{Flo...)
 discriminator = (CMeanGaussian{DiagVar}(mapping=Dense(4, 1, σ), σ2=1-element Array...)
```
"""
struct GAN{P<:ContinuousMultivariateDistribution,G<:ACD,D<:ACD} <: AbstractGAN
	prior::P
	generator::G
	discriminator::D
end

Flux.@functor GAN

function Flux.trainable(m::GAN)
    (generator=m.generator, discriminator=m.discriminator)
end

function GAN(zlength::Int, g::ACD, d::ACD)
    T = eltype(first(Flux.params(g)))
    prior = DistributionsAD.TuringMvNormal(zeros(T,zlength), ones(T,zlength))
    GAN(prior, g, d)
end

"""
	generator_loss(m::GAN, z::AbstractArray)
	generator_loss(m::GAN, batchsize::Int)

Loss of the GAN generator. The input is either the random code `z` or
`batchsize` of samples to generate from the model prior and compute the loss
from.
"""
function generator_loss(m::GAN, z::AbstractArray)
    T = eltype(m.prior)
    x = mean(m.generator,z)
    y = mean(m.discriminator, x)
    generator_loss(mean(m.discriminator, mean(m.generator,z)))
end

generator_loss(m::GAN, batchsize::Int) =
    generator_loss(m, rand(m.prior, batchsize))

"""
	discriminator_loss(m::GAN, x::AbstractArray[, z::AbstractArray])

Loss of the GAN discriminator given a batch of training samples `x` and latent
prior samples `z`.  If z is not given, it is automatically generated from the
model prior.
"""
function discriminator_loss(m::GAN, x::AbstractArray, z::AbstractArray)
    T = eltype(m.prior)
    st = mean(m.discriminator,x)
    sg = mean(m.discriminator, mean(m.generator,z))
    discriminator_loss(st, sg)
end

discriminator_loss(m::GAN, x::AbstractArray) =
    discriminator_loss(m, x, rand(m.prior, size(x,2)))

function Base.show(io::IO, m::AbstractGAN)
    p = repr(m.prior)
    p = sizeof(p)>70 ? "($(p[1:70-3])...)" : p
    g = repr(m.generator)
    g = sizeof(g)>70 ? "($(g[1:70-3])...)" : g
    d = repr(m.discriminator)
    d = sizeof(d)>70 ? "($(d[1:70-3])...)" : d
    msg = """$(nameof(typeof(m))):
     prior   = $(p)
     generator = $(g)
     discriminator = $(d)
    """
    print(io, msg)
end


text:
#===============================================================================

    keystring_constants.jl

    Author: Tom Price
    Date: June 2019

    header file for QpcrAnalysis.jl module
    defines alphanumeric Dict keys

===============================================================================#

const WATER_KEY = "water"
const FLUORESCENCE_VALUE_KEY = "fluorescence_value"
const CYCLE_NUM_KEY = "cycle_num"
const WELL_NUM_KEY = "well_num"
const CHANNEL_KEY = "channel"
const RAW_DATA_KEY = "raw_data"
const CATEG_WELL_VEC_KEY = "categ_well_vec"
const CTRL_WELL_DICT_KEY = "ctrl_well_dict"
const CQ_METHOD_KEY = "cq_method"
const CALIBRATION_INFO_KEY = "calibration_info"
const BASELINE_METHOD_KEY = "baseline_method"
const SIGMOID_KEY = "sigmoid"
const LINEAR_KEY = "linear"
const MEDIAN_KEY = "median"
const TEMPERATURE_KEY = "temperature"
const WELL_KEY = "well"
const TARGET_KEY = "target"
const CQ_KEY = "cq"
const QUANTITY_KEY = "quantity"
const SAMPLE_KEY = "sample"
const BASELINE_KEY = "baseline"
const EXCITATION_KEY = "excitation"
const LID_TEMP_KEY = "lid_temp"
const ELAPSED_TIME_KEY = "elapsed_time"
const HEAT_BLOCK_ZONE_1_TEMP_KEY = "heat_block_zone_1_temp"
const HEAT_BLOCK_ZONE_2_TEMP_KEY = "heat_block_zone_2_temp"

# const STEP_ID_KEY = "step_id"
# const RAMP_ID_KEY = "ramp_id"
# const STAGE_ID_KEY = "stage_id"

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
@testset "Test Line-Search Stopping with LineModel" begin

  include("rosenbrock.jl")
  nlp = ADNLPModel(rosenbrock, 0.5ones(6))
  #We create a LineModel from nlp.meta.x0
  #The pb can be a LineModel defined in [SolverTools.jl](https://github.com/JuliaSmoothOptimizers/SolverTools.jl)
  d = - grad(nlp, nlp.meta.x0)
  pb = LineModel(nlp, nlp.meta.x0, d)
  state = OneDAtX(0.0, f₀ = obj(pb, 0.), g₀ = grad(pb, 0.))
  @test state.x == 0.0
  @test !isnan(state.f₀)
  @test !isnan(state.g₀)
  @test isnan(state.fx)
  @test isnan(state.gx)

  τ₀ = 1e-4 #Armijo parameter
  meta = StoppingMeta(optimality_check = (pb, lsatt) -> lsatt.fx - lsatt.f₀ - lsatt.g₀ * lsatt.x * τ₀,
                      tol_check = (a,b,c) -> 0.0,
                      tol_check_neg = (atol, rtol, opt0) -> -Inf)
  src  = StopRemoteControl(resources_check = false)

  stp = NLPStopping(pb, meta, src, state)

  #Independent of x
  @test !Stopping._infeasibility_check!(stp, 1.)
  @test !Stopping._resources_check!(stp, 1.)
  #Depend on x
  #Stopping._unbounded_problem_check!(stp, 1.) #doesn't work if there is no state.fx

  @test start!(stp) #NaN is optimal
  @test isnan(stp.current_state.current_score)
  @test stop!(stp)
  @test !stp.meta.infeasible

  #Simulate a backtracking
  reinit!(stp)
  reset!(pb)
  t = 1.0
  i = 0
  OK = update_and_start!(stp, x = t, fx = obj(stp.pb, t), gx = grad(stp.pb, t))
  while !OK
    t /= 2
    i += 1
    OK = update_and_stop!(stp, x = t, fx = obj(stp.pb, t), gx = grad(stp.pb, t))
  end
  @test i == stp.meta.nb_of_stop
  @test neval_obj(stp.pb) == i+1
  @test neval_grad(stp.pb) == i+1
  @test stp.meta.optimal
  @test :Optimal in status(stp, list = true)

  #With the nlp as a main_stp
  #The issue is that pb.meta.x0 is a vector of size 1!, so NLPStopping(pb) is not possible
  stp2 = NLPStopping(pb, OneDAtX(0., f₀ = obj(pb, 0.), g₀ = grad(pb, 0.)), main_stp = NLPStopping(nlp), optimality_check = armijo)

  @test update_and_start!(stp2, fx = 0.0, gx = 0.0)
  @test update_and_stop!(stp2, fx = 0.0, gx = 0.0)

  #Check an OR condition
  lsatt = OneDAtX(0.0, f₀ = obj(pb, 0.), g₀ = grad(pb, 0.), (NaN, NaN))
  @test isnan.(lsatt.current_score) == (true, true)
  @test lsatt.x == 0.0

  #Check an AND condition
  lsatt = OneDAtX(0.0, f₀ = obj(pb, 0.), g₀ = grad(pb, 0.), [Inf, Inf])
  τ₀ = 1e-4 #Armijo parameter
  τ₁ = 0.99 #Wolfe parameter
  meta = StoppingMeta(optimality_check = (pb, lsatt) -> vcat(lsatt.fx - lsatt.f₀ - lsatt.g₀ * lsatt.x * τ₀,
                                                             abs(lsatt.gx) - τ₁ * abs(lsatt.g₀)),
                      tol_check = (a,b,c) -> 0.0,
                      tol_check_neg = (atol, rtol, opt0) -> -Inf)
  src  = StopRemoteControl(resources_check = false)

  stp = NLPStopping(pb, meta, src, lsatt)

  #Simulate a backtracking
  reinit!(stp)
  reset!(pb)
  t = 1.0
  i = 0
  OK = update_and_start!(stp, x = t, fx = obj(stp.pb, t), gx = grad(stp.pb, t))
  while !OK
    t /= 2
    i += 1
    OK = update_and_stop!(stp, x = t, fx = obj(stp.pb, t), gx = grad(stp.pb, t))
  end
  @test i == stp.meta.nb_of_stop
  @test neval_obj(stp.pb) == i+1
  @test neval_grad(stp.pb) == i+1
  @test stp.meta.optimal
  @test :Optimal in status(stp, list = true)

  reinit!(stp, rcounters = true, rstate = true)
  fill_in!(stp, 1.)
  @test !isnan(stp.current_state.fx)
  @test !isnan(stp.current_state.gx)
  @test !isnan(stp.current_state.f₀)
  @test !isnan(stp.current_state.g₀)

end


text:
using FactCheck

facts("Math module tests") do
  @fact 1 --> 1 "Primary numbers less than 100 search"
  @pending 1-->isempty([]) ""
  context("Trigonometric functions") do
    @fact 1 --> 1 "Sine"
    @fact 1 --> 2 "Cosine"
	@fact 1 --> 1 "Tangent"
  end
  @fact 2 --> 2 "Units of measurment"
end

facts("Words processing tests") do
  @fact error("test2") --> 1 "Make palindrome"
  @fact 1--> 2 "Make anagrame"
  context("Translate") do
	facts("From English") do
		@fact 1--> 1 "From English to Spanish"
		@fact 1--> 2 "From English to German"
	end
    facts("To English") do
      @fact 1--> 1 "From Spanish to English"
      @pending 1--> 2 "From German to English"
    end
  end
end


text:
using Base.Filesystem
const fs = Base.Filesystem

path_base = dirname(pwd())
path_sufix = "/Datasets/Dataset_Bovespa/"
path_destiny = "Unzipped/"

location = path_base * path_sufix * path_destiny

file_prefix = "COTAHIST.A"
file_start_year = 1995
file_end_year = 2000

file_prefix_desired = "COTAHIST_A"
file_extension = ".TXT"

for year in range(file_start_year, stop=file_end_year)
    path_original = location * file_prefix * string(year)
    path_renamed = location * file_prefix_desired * string(year) * file_extension
    fs.mv(path_original, path_renamed)
end

path_original = location * file_prefix_desired * "2001"
path_renamed = location * file_prefix_desired * "2001" * file_extension
fs.mv(path_original, path_renamed)

println("Files from 1995 to 2001 were successfully converted!!")