text:
using MCCGA
using Test

@testset "Bit operations" begin

    @testset "Bits to float" begin
        bitsofnumber = [
            0,
            0,
            1,
            1,
            1,
            1,
            1,
            1,
            1,
            0,
            0,
            1,
            0,
            0,
            0,
            1,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0
        ]

        floatlist = floats(bitsofnumber)

        @test floatlist isa Vector
        @test length(floatlist) == 1
        @test floatlist[1] == 1.1328125
    end

    @testset "Float to bits" begin
        fvalue = 76289.375
        fbits = bits(fvalue)

        @test length(fbits) == 32
        @test fbits == [0,1,0,0,0,1,1,1,1,0,0,1,0,1,0,1,0,0,0,0,0,0,0,0,1,0,1,1,0,0,0,0]
    end

    @testset "Bit-float conversations" begin
        tol = 0.001

        v = [2.71828, 3.14159265, 1.0, 5, -100]

        mybits = bits(v)

        myfloats = floats(mybits)

        @test length(mybits) == 32 * length(v)
        @test length(myfloats) == length(v)

        for i = 1:length(v)
            @test isapprox(v[i], myfloats[i], atol = tol)
        end
    end

    @testset "Sign bit test" begin
        mybits1 = bits(9.8)
        mybits2 = bits(-9.8)

        @test mybits1[1] == 0
        @test mybits2[1] == 1
    end

    @testset "Infinity bits" begin

        infbits = [
            0,
            1,
            1,
            1,
            1,
            1,
            1,
            1,
            1,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
        ]

        mybits1 = bits(Inf64)
        mybits2 = bits(Inf32)
        mybits3 = bits(Inf)

        @test mybits1 == infbits
        @test mybits2 == infbits
        @test mybits3 == infbits
    end

    @testset "NaN bits" begin
        nanbits = [
            0,
            1,
            1,
            1,
            1,
            1,
            1,
            1,
            1,
            1,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
            0,
        ]

        mybits1 = bits(NaN64)
        mybits2 = bits(NaN32)
        mybits3 = bits(NaN)

        @test mybits1 == nanbits
        @test mybits2 == nanbits
        @test mybits3 == nanbits
    end
end

@testset "Sampling" begin
    @testset "Convergency state" begin
        probvector = [1.0, 1.0, 1.0]
        mysample = MCCGA.sample(probvector)

        @test length(mysample) == 3
        @test mysample[1] == 1.0
        @test mysample[2] == 1.0
        @test mysample[3] == 1.0
    end

    @testset "Normal state" begin
        probvector = [0.5, 0.5, 0.5]
        mysample = MCCGA.sample(probvector)

        @test length(mysample) == 3
        @test mysample[1] >= 0.0
        @test mysample[2] >= 0.0
        @test mysample[3] >= 0.0
        @test mysample[1] <= 1.0
        @test mysample[2] <= 1.0
        @test mysample[3] <= 1.0
    end
end

@testset "Optimization of Pi and E" begin
    function f(x)
        return (x[1] - 3.14159265)^2 + (x[2] - exp(1.0))^2
    end

    lower = [-100.0, -100.0]
    upper = [100.0, 100.0]

    result = MCCGA.mccga(
        lower = lower,
        upper = upper,
        costfunction = f,
        popsize = 100,
        maxsamples = 10000,
    )

    @test result isa Dict

    finalsolution = result["final_solution"]

    @test finalsolution[1] >= 3
    @test finalsolution[1] <= 4

    @test finalsolution[2] >= 2
    @test finalsolution[2] <= 3

end

@testset "Chichinadze" begin

    tol = 0.001

    lower = [-30.0, -30.0]
    upper = [30.0, 30.0]

    result = MCCGA.mccga(
        lower = lower,
        upper = upper,
        costfunction = Functions.chichinadze,
        popsize = 200,
        maxsamples = 10000,
    )

    finalsolution = result["final_solution"]
    finalminimum = result["final_minimum"]

    @test isapprox(finalsolution[1], 5.90133, atol = tol)
    @test isapprox(finalsolution[2], 0.5, atol = tol)
    @test isapprox(finalminimum, -43.3159, atol = tol)
end

@testset "Levy" begin

    tol = 0.001

    lower = [-10.0, -10.0, -10]
    upper = [10.0, 10.0, 10]

    result = MCCGA.mccga(
        lower = lower,
        upper = upper,
        costfunction = Functions.levy,
        popsize = 100,
        maxsamples = 10000,
    )

    finalsolution = result["final_solution"]
    finalminimum = result["final_minimum"]

    @test isapprox(finalsolution[1], 1.0, atol = tol)
    @test isapprox(finalsolution[2], 1.0, atol = tol)
    @test isapprox(finalsolution[3], 1.0, atol = tol)
    @test isapprox(finalminimum, 0.0, atol = tol)
end

@testset "Ackley" begin

    tol = 0.001

    lower = [-32.768, -32.768, -32.768]
    upper = [32.768, 32.768, 32.768]

    result = MCCGA.mccga(
        lower = lower,
        upper = upper,
        costfunction = Functions.ackley,
        popsize = 100,
        maxsamples = 10000,
    )

    finalsolution = result["final_solution"]
    finalminimum = result["final_minimum"]

    @test isapprox(finalsolution[1], 0.0, atol = tol)
    @test isapprox(finalsolution[2], 0.0, atol = tol)
    @test isapprox(finalsolution[3], 0.0, atol = tol)
    @test isapprox(finalminimum, 0.0, atol = tol)
end

max_stars_repo_name: jbytecode/MCCGA

max_issues_repo_name: jbytecode/MCCGA

max_forks_repo_name: jbytecode/MCCGA


text:
module MediaCLI
    include("transfer.jl")
    include("deflate.jl")

    export
            MediaConfig,
            Episode,
            Season,
            TVShow,
            transfer,
            deflate,
            discover_episodes,
            discover_seasons,
            plex_denomination
end

ext: jl


text:
#!/bin/env julia

const filename = "latex2unicode_utf-8"
const include_emoji = false # set this to true if you want to keymap empjis as well

# We want to avoid situations in which the user types e.g. \delt and pauses,
# and the result is "∇t" because "\del" gets recognized and then there is some leftover "t".
# This allows us to get completions with <Tab> for example.
function fix_completions(completions::Dict{String,String})
    allinputs = sort!(collect(keys(completions)))
    new_completions = copy(completions)

    for input in allinputs
        chars = completions[input]
        l = length(input)
        longer = filter(x->startswith(x, input)&&length(x)>l, allinputs)
        n = length(longer)
        n == 0 && continue
        new_completions[input * "<Tab>"] = chars
        for other in longer
            add = other[1:(l+1)]
            get!(new_completions, add, add)
        end
    end
    return new_completions
end

function unicode_data()
    file = normpath(Sys.BINDIR, "..", "..", "doc", "UnicodeData.txt")
    names = Dict{UInt32, String}()
    open(file) do unidata
        for line in readlines(unidata)
            id, name, desc = split(line, ";")[[1, 2, 11]]
            codepoint = parse(UInt32, "0x$id")
            names[codepoint] = (name == "" ? desc : desc == "" ? name : "$name / $desc")
        end
    end
    return names
end

# Prepend a dotted circle ('◌' i.e. '\u25CC') to combining characters
function fix_combining_chars(char)
    cat = Base.Unicode.category_code(char)
    return string(cat == 6 || cat == 8 ? "◌" : "", char)
end

function table_entries(completions::Dict{String,String}, unicode_dict)
    latex = String[]
    code = String[]
    unicode = String[]
    desc = String[]

    for (input, chars) in sort!(collect(completions))
        code_points, unicode_names, characters = String[], String[], String[]
        if startswith(chars, "\\")
            push!(code_points, replace(chars, "\\" => "\\\\"))
            push!(unicode_names, "(Incomplete sequence)")
            push!(characters, "")
        else
            for char in chars
                push!(code_points, "<char-0x$(uppercase(string(UInt32(char), base = 16, pad = 5)))>")
                push!(unicode_names, get(unicode_dict, UInt32(char), "(No Unicode name)"))
                push!(characters, isempty(characters) ? fix_combining_chars(char) : "$char")
            end
        end
        push!(latex, replace(input, "\\"=>"\\\\"))
        push!(code, join(code_points))
        push!(unicode, join(characters))
        push!(desc, join(unicode_names, " + "))
    end
    return latex, code, unicode, desc
end

open("$filename.vim","w") do f
    print(f, """
        " This file is autogenerated from the script '$(basename(Base.source_path()))'
        " The symbols are based on Julia version $VERSION
        " The full generating script can be found in the comments at the bottom of this file,
        " and it can be extracted with:
        "
        "     \$ grep '^\">' $filename.vim | cut -c4- > $(basename(Base.source_path()))'
        "
        " To produce this keymap file you need to have Julia compilied from source, and
        " to run `make UnicodeData.txt` inside Julia's `doc` directory.
        " Then you can run:
        "
        "     \$ julia $(basename(Base.source_path())).jl
        "

        scriptencoding utf-8

        let b:keymap_name = "L2U"
        loadkeymap

        """)

    col_headers = ["\" Tab completion sequence", "Code point", "Character", "Unicode name"]

    orig_completions = include_emoji ? merge(
                REPL.REPLCompletions.latex_symbols,
                REPL.REPLCompletions.emoji_symbols
                ) :
                REPL.REPLCompletions.latex_symbols

    latex, code, unicode, desc =
        table_entries(
            fix_completions(orig_completions),
            unicode_data()
            )

    lw = max(length(col_headers[3]), maximum(map(length, latex)))
    cw = max(length(col_headers[1]), maximum(map(length, code)))
    uw = max(length(col_headers[2]), maximum(map(length, unicode)))
    dw = max(length(col_headers[4]), maximum(map(length, desc)))

    print_padded(l, c, u, d) = println(f, rpad(l, lw), " ", rpad(c, cw), " \" ", rpad(u, uw), " : ", d)

    print_padded(col_headers...)
    print_padded("\" " * "-"^(lw-2), "-"^cw, "-"^uw, "-"^dw)

    for (l, c, u, d) in zip(latex, code, unicode, desc)
        print_padded(l, c, u, d)
    end
    print_padded("\" " * "-"^(lw-2), "-"^cw, "-"^uw, "-"^dw)

    print(f, """

        " Below here is the script that was used to produce this file.

        """)
    for l in readlines(Base.source_path())
        println(f, "\"> ", l)
    end
    println(f)
end

max_stars_repo_name: vvbond/julia-vim

max_issues_repo_name: KevinGG/julia-vim

max_forks_repo_name: KevinGG/julia-vim


text:
using Test
using UnsafeArrays, BangBang
using LinearAlgebra

using MuJoCo
using MuJoCo: TESTMODELXML

using MuJoCo.MJCore
using MuJoCo.MJCore: mjtNum, mjNSOLVER, mjNREF

const TESTMODELMJB = joinpath(@__DIR__, "testmodel.mjb")

rm(TESTMODELMJB, force=true)
mj_saveModel(mj_loadXML(TESTMODELXML), TESTMODELMJB)

function makemd()
    pm = mj_loadXML(TESTMODELXML)
    pd = mj_makeData(pm)
    m = unsafe_load(pm)
    d = unsafe_load(pd)
    mj_resetData(m, d)
    mj_forward(m, d)
    m, d, pm, pd
end

function jlmakemd()
    pm = mj_loadXML(TESTMODELXML)
    pd = mj_makeData(pm)
    m = unsafe_load(pm)
    d = unsafe_load(pd)
    mj_resetData(m, d)
    mj_forward(m, d)
    m, d, jlModel(pm), jlData(pm, pd)
end

function test_uninit(T)
    zeroarg = T()
    posonly = T((getfield(zeroarg, name) for name in fieldnames(T))...)
    isequal(zeroarg, posonly) && hash(zeroarg) == hash(posonly)
end

@testset "MuJoCo.jl" begin
    include("MJCore/MJCore.jl")
    include("sugar.jl")
end


text:
"""
    solve_objbound_edge_cases(model::MOI.ModelLike, config::Config)

Test a variety of edge cases related to the ObjectiveBound attribute.
"""
function solve_objbound_edge_cases(model::MOI.ModelLike, config::Config)
    @testset "Min IP with constant" begin
        MOI.empty!(model)
        @test MOI.is_empty(model)
        MOIU.loadfromstring!(
            model,
            """
    variables: x
    minobjective: 2.0x + -1.0
    x >= 1.5
    x in Integer()
""",
        )
        x = MOI.get(model, MOI.VariableIndex, "x")
        test_model_solution(
            model,
            config;
            objective_value = 3.0,
            variable_primal = [(x, 2.0)],
        )
        if config.solve
            @test MOI.get(model, MOI.ObjectiveBound()) <= 3.0
        end
    end
    @testset "Max IP with constant" begin
        MOI.empty!(model)
        @test MOI.is_empty(model)
        MOIU.loadfromstring!(
            model,
            """
    variables: x
    maxobjective: 2.0x + 1.0
    x <= 1.5
    x in Integer()
""",
        )
        x = MOI.get(model, MOI.VariableIndex, "x")
        test_model_solution(
            model,
            config;
            objective_value = 3.0,
            variable_primal = [(x, 1.0)],
        )
        if config.solve
            @test MOI.get(model, MOI.ObjectiveBound()) >= 3.0
        end
    end
    @testset "Min LP with constant" begin
        MOI.empty!(model)
        @test MOI.is_empty(model)
        MOIU.loadfromstring!(
            model,
            """
    variables: x
    minobjective: 2.0x + -1.0
    x >= 1.5
""",
        )
        x = MOI.get(model, MOI.VariableIndex, "x")
        test_model_solution(
            model,
            config;
            objective_value = 2.0,
            variable_primal = [(x, 1.5)],
        )
        if config.solve
            @test MOI.get(model, MOI.ObjectiveBound()) <= 2.0
        end
    end
    @testset "Max LP with constant" begin
        MOI.empty!(model)
        @test MOI.is_empty(model)
        MOIU.loadfromstring!(
            model,
            """
    variables: x
    maxobjective: 2.0x + 1.0
    x <= 1.5
""",
        )
        x = MOI.get(model, MOI.VariableIndex, "x")
        test_model_solution(
            model,
            config;
            objective_value = 4.0,
            variable_primal = [(x, 1.5)],
        )
        if config.solve
            @test MOI.get(model, MOI.ObjectiveBound()) >= 4.0
        end
    end
end
unittests["solve_objbound_edge_cases"] = solve_objbound_edge_cases

function solve_unbounded_model(model::MOI.ModelLike, config::Config)
    MOI.empty!(model)
    x = MOI.add_variables(model, 5)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.DUAL_INFEASIBLE
    end
end
unittests["solve_unbounded_model"] = solve_unbounded_model

function solve_single_variable_dual_min(model::MOI.ModelLike, config::Config)
    MOI.empty!(model)
    x = MOI.add_variable(model)
    xl = MOI.add_constraint(model, x, MOI.GreaterThan(1.0))
    xu = MOI.add_constraint(model, x, MOI.LessThan(1.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    if config.solve && config.duals
        MOI.optimize!(model)
        @test isapprox(
            MOI.get(model, MOI.VariablePrimal(), x),
            1.0,
            atol = config.atol,
        )
        sl = MOI.get(model, MOI.ConstraintDual(), xl)
        su = MOI.get(model, MOI.ConstraintDual(), xu)
        @test isapprox(sl + su, 1.0, atol = config.atol)
        @test sl >= -config.atol
        @test su <= config.atol
    end
end
unittests["solve_single_variable_dual_min"] = solve_single_variable_dual_min

function solve_single_variable_dual_max(model::MOI.ModelLike, config::Config)
    MOI.empty!(model)
    x = MOI.add_variable(model)
    xl = MOI.add_constraint(model, x, MOI.GreaterThan(1.0))
    xu = MOI.add_constraint(model, x, MOI.LessThan(1.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    if config.solve && config.duals
        MOI.optimize!(model)
        @test isapprox(
            MOI.get(model, MOI.VariablePrimal(), x),
            1.0,
            atol = config.atol,
        )
        sl = MOI.get(model, MOI.ConstraintDual(), xl)
        su = MOI.get(model, MOI.ConstraintDual(), xu)
        @test isapprox(sl + su, -1.0, atol = config.atol)
        @test sl >= -config.atol
        @test su <= config.atol
    end
end
unittests["solve_single_variable_dual_max"] = solve_single_variable_dual_max

function solve_result_index(model::MOI.ModelLike, config::Config)
    atol = config.atol
    rtol = config.rtol
    MOI.empty!(model)
    x = MOI.add_variable(model)
    c = MOI.add_constraint(model, x, MOI.GreaterThan(1.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    if config.solve
        MOI.optimize!(model)
        result_count = MOI.get(model, MOI.ResultCount())
        function result_err(attr)
            return MOI.ResultIndexBoundsError{typeof(attr)}(attr, result_count)
        end
        result_index = result_count + 1
        @test MOI.get(model, MOI.ObjectiveValue(1)) ≈ 1.0 atol = atol rtol =
            rtol
        @test_throws result_err(MOI.ObjectiveValue(result_index)) MOI.get(
            model,
            MOI.ObjectiveValue(result_index),
        )
        if config.dual_objective_value
            @test MOI.get(model, MOI.DualObjectiveValue(1)) ≈ 1.0 atol = atol rtol =
                rtol
            @test_throws result_err(MOI.DualObjectiveValue(result_index)) MOI.get(
                model,
                MOI.DualObjectiveValue(result_index),
            )
        end
        @test MOI.get(model, MOI.PrimalStatus(1)) == MOI.FEASIBLE_POINT
        @test MOI.get(model, MOI.PrimalStatus(result_index)) == MOI.NO_SOLUTION
        @test MOI.get(model, MOI.VariablePrimal(1), x) ≈ 1.0 atol = atol rtol =
            rtol
        @test_throws result_err(MOI.VariablePrimal(result_index)) MOI.get(
            model,
            MOI.VariablePrimal(result_index),
            x,
        )
        @test MOI.get(model, MOI.ConstraintPrimal(1), c) ≈ 1.0 atol = atol rtol =
            rtol
        @test_throws result_err(MOI.ConstraintPrimal(result_index)) MOI.get(
            model,
            MOI.ConstraintPrimal(result_index),
            c,
        )
        if config.duals
            @test MOI.get(model, MOI.DualStatus(1)) == MOI.FEASIBLE_POINT
            @test MOI.get(model, MOI.DualStatus(result_index)) ==
                  MOI.NO_SOLUTION
            @test MOI.get(model, MOI.ConstraintDual(1), c) ≈ 1.0 atol = atol rtol =
                rtol
            @test_throws result_err(MOI.ConstraintDual(result_index)) MOI.get(
                model,
                MOI.ConstraintDual(result_index),
                c,
            )
        end
    end
end
unittests["solve_result_index"] = solve_result_index

function solve_farkas_equalto_upper(model::MOI.ModelLike, config::Config)
    MOI.empty!(model)
    x = MOI.add_variables(model, 2)
    clb = MOI.add_constraint.(model, x, MOI.GreaterThan(0.0))
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
        MOI.EqualTo(-1.0),
    )
    if config.solve && config.infeas_certificates
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
        @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
        clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
        c_dual = MOI.get(model, MOI.ConstraintDual(), c)
        @test clb_dual[1] > config.atol
        @test clb_dual[2] > config.atol
        @test c_dual[1] < -config.atol
        @test clb_dual ≈ [2, 1] .* -c_dual atol = config.atol rtol = config.rtol
    end
end
unittests["solve_farkas_equalto_upper"] = solve_farkas_equalto_upper

function solve_farkas_equalto_lower(model::MOI.ModelLike, config::Config)
    MOI.empty!(model)
    x = MOI.add_variables(model, 2)
    clb = MOI.add_constraint.(model, x, MOI.GreaterThan(0.0))
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-2.0, -1.0], x), 0.0),
        MOI.EqualTo(1.0),
    )
    if config.solve && config.infeas_certificates
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
        @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
        clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
        c_dual = MOI.get(model, MOI.ConstraintDual(), c)
        @test clb_dual[1] > config.atol
        @test clb_dual[2] > config.atol
        @test c_dual[1] > config.atol
        @test clb_dual ≈ [2, 1] .* c_dual atol = config.atol rtol = config.rtol
    end
end
unittests["solve_farkas_equalto_lower"] = solve_farkas_equalto_lower

function solve_farkas_lessthan(model::MOI.ModelLike, config::Config)
    MOI.empty!(model)
    x = MOI.add_variables(model, 2)
    clb = MOI.add_constraint.(model, x, MOI.GreaterThan(0.0))
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
        MOI.LessThan(-1.0),
    )
    if config.solve && config.infeas_certificates
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
        @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
        clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
        c_dual = MOI.get(model, MOI.ConstraintDual(), c)
        @test clb_dual[1] > config.atol
        @test clb_dual[2] > config.atol
        @test c_dual[1] < -config.atol
        @test clb_dual ≈ [2, 1] .* -c_dual atol = config.atol rtol = config.rtol
    end
end
unittests["solve_farkas_lessthan"] = solve_farkas_lessthan

function solve_farkas_greaterthan(model::MOI.ModelLike, config::Config)
    MOI.empty!(model)
    x = MOI.add_variables(model, 2)
    clb = MOI.add_constraint.(model, x, MOI.GreaterThan(0.0))
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-2.0, -1.0], x), 0.0),
        MOI.GreaterThan(1.0),
    )
    if config.solve && config.infeas_certificates
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
        @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
        clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
        c_dual = MOI.get(model, MOI.ConstraintDual(), c)
        @test clb_dual[1] > config.atol
        @test clb_dual[2] > config.atol
        @test c_dual[1] > config.atol
        @test clb_dual ≈ [2, 1] .* c_dual atol = config.atol rtol = config.rtol
    end
end
unittests["solve_farkas_greaterthan"] = solve_farkas_greaterthan

function solve_farkas_interval_upper(model::MOI.ModelLike, config::Config)
    MOI.empty!(model)
    x = MOI.add_variables(model, 2)
    clb = MOI.add_constraint.(model, x, MOI.GreaterThan(0.0))
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
        MOI.Interval(-2.0, -1.0),
    )
    if config.solve && config.infeas_certificates
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
        @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
        clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
        c_dual = MOI.get(model, MOI.ConstraintDual(), c)
        @test clb_dual[1] > config.atol
        @test clb_dual[2] > config.atol
        @test c_dual[1] < -config.atol
        @test clb_dual ≈ [2, 1] .* -c_dual atol = config.atol rtol = config.rtol
    end
end
unittests["solve_farkas_interval_upper"] = solve_farkas_interval_upper

function solve_farkas_interval_lower(model::MOI.ModelLike, config::Config)
    MOI.empty!(model)
    x = MOI.add_variables(model, 2)
    clb = MOI.add_constraint.(model, x, MOI.GreaterThan(0.0))
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-2.0, -1.0], x), 0.0),
        MOI.Interval(1.0, 2.0),
    )
    if config.solve && config.infeas_certificates
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
        @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
        clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
        c_dual = MOI.get(model, MOI.ConstraintDual(), c)
        @test clb_dual[1] > config.atol
        @test clb_dual[2] > config.atol
        @test c_dual[1] > config.atol
        @test clb_dual ≈ [2, 1] .* c_dual atol = config.atol rtol = config.rtol
    end
end
unittests["solve_farkas_interval_lower"] = solve_farkas_interval_lower

function solve_farkas_variable_lessthan(model::MOI.ModelLike, config::Config)
    MOI.empty!(model)
    x = MOI.add_variables(model, 2)
    clb = MOI.add_constraint.(model, x, MOI.LessThan(0.0))
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
        MOI.GreaterThan(1.0),
    )
    if config.solve && config.infeas_certificates
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
        @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
        clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
        c_dual = MOI.get(model, MOI.ConstraintDual(), c)
        @test clb_dual[1] < -config.atol
        @test clb_dual[2] < -config.atol
        @test c_dual[1] > config.atol
        @test clb_dual ≈ [2, 1] .* -c_dual atol = config.atol rtol = config.rtol
    end
end
unittests["solve_farkas_variable_lessthan"] = solve_farkas_variable_lessthan

function solve_farkas_variable_lessthan_max(
    model::MOI.ModelLike,
    config::Config,
)
    MOI.empty!(model)
    x = MOI.add_variables(model, 2)
    clb = MOI.add_constraint.(model, x, MOI.LessThan(0.0))
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], x), 0.0),
        MOI.GreaterThan(1.0),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x[1])
    if config.solve && config.infeas_certificates
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.INFEASIBLE
        @test MOI.get(model, MOI.DualStatus()) == MOI.INFEASIBILITY_CERTIFICATE
        clb_dual = MOI.get.(model, MOI.ConstraintDual(), clb)
        c_dual = MOI.get(model, MOI.ConstraintDual(), c)
        @test clb_dual[1] < -config.atol
        @test clb_dual[2] < -config.atol
        @test c_dual[1] > config.atol
        @test clb_dual ≈ [2, 1] .* -c_dual atol = config.atol rtol = config.rtol
    end
end
unittests["solve_farkas_variable_lessthan_max"] =
    solve_farkas_variable_lessthan_max

function solve_twice(model::MOI.ModelLike, config::Config)
    MOI.empty!(model)
    x = MOI.add_variable(model)
    c = MOI.add_constraint(model, x, MOI.GreaterThan(1.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    if config.solve
        MOI.optimize!(model)
        MOI.optimize!(model)
        MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        MOI.get(model, MOI.VariablePrimal(), x) == 1.0
    end
end
unittests["solve_twice"] = solve_twice


text:
export WLS
using Wavelets

"""
    WLS(surromethod::Surrogate = IAAFT(),
        rescale::Bool = true,
        wt::Wavelets.WT.OrthoWaveletClass = Wavelets.WT.Daubechies{16}())

A wavelet surrogate generated by taking the maximal overlap discrete
wavelet transform (MODWT) of the signal, shuffling detail
coefficients at each dyadic scale using the provided `surromethod`,
then taking the inverse transform to obtain a surrogate.

## Coefficient shuffling method

In contrast to the original
implementation where IAAFT is used, you may choose to use any surrogate
method from this package to perform the randomization of the detail
coefficients at each dyadic scale. Note: The iterative procedure after
the rank ordering step (step [v] in [^Keylock2006]) is not performed in
this implementation.

If `surromethod == IAAFT()`, the wavelet surrogates preserves the local
mean and variance structure of the signal, but randomises nonlinear
properties of the signal (i.e. Hurst exponents)[^Keylock2006]. These surrogates
can therefore be used to test for changes in nonlinear properties of the
original signal.

In contrast to IAAFT surrogates, the IAAFT-wavelet surrogates also
preserves nonstationarity. Using other `surromethod`s does not necessarily
preserve nonstationarity.

To deal with nonstationary signals, Keylock (2006) recommends using a
wavelet with a high number of vanishing moments. Thus, the default is to
use a Daubechies wavelet with 16 vanishing moments.

## Rescaling

If `rescale == true`, then surrogate values are mapped onto the
values of the original time series, as in the [`AAFT`](@ref) algorithm.
If `rescale == false`, surrogate values are not constrained to the
original time series values.

[^Keylock2006]: C.J. Keylock (2006). "Constrained surrogate time series with preservation of the mean and variance structure". Phys. Rev. E. 73: 036707. doi:10.1103/PhysRevE.73.036707.
"""
struct WLS{WT <: Wavelets.WT.OrthoWaveletClass, S <: Surrogate} <: Surrogate
    surromethod::Surrogate # should preserve values of the original series
    rescale::Bool
    wt::WT

    function WLS(method::S = IAAFT(), rescale::Bool = true, wt::WT = Wavelets.WT.Daubechies{16}()) where {S <: Surrogate, WT <: Wavelets.WT.OrthoWaveletClass}
        new{WT, S}(method, rescale, wt)
    end
end

function surrogenerator(x::AbstractVector{T}, method::WLS) where T
    wl = wavelet(method.wt)
    L = length(x)
    x_sorted = sort(x)

    # Wavelet coefficients (step [i] in Keylock)
    W = modwt(x, wl)
    Nscales = ndyadicscales(L)

    # Will contain surrogate realizations of the wavelet coefficients
    # at each scale (step [ii] in Keylock).
    sW = zeros(T, size(W))

    # We will also need a matrix to store the mirror images of the
    # surrogates (last part of step [ii])
    sWmirr = zeros(T, size(W))

    # Surrogate generators for each set of coefficients
    sgs = [surrogenerator(W[:, i], method.surromethod) for i = 1:Nscales]

    # Temporary array for the circular shift error minimizing step
    circshifted_s = zeros(T, size(W))
    circshifted_smirr = zeros(T, size(W))

    init = (wl = wl, W = W, Nscales = Nscales, L = L,
            sW = sW, sgs = sgs, sWmirr = sWmirr,
            circshifted_s = circshifted_s,
            circshifted_smirr = circshifted_smirr,
            x_sorted = x_sorted)

    return SurrogateGenerator(method, x, init)
end

function (sg::SurrogateGenerator{<:WLS})()
    fds = (:wl, :W, :Nscales, :L, :sW, :sgs, :sWmirr,
        :circshifted_s, :circshifted_smirr,
        :x_sorted)

    wl, W, Nscales, L, sW, sgs, sWmirr,
        circshifted_s, circshifted_smirr,
        x_sorted = getfield.(Ref(sg.init), fds)

    # Create surrogate versions of detail coefficients at each dyadic scale [first part of step (ii) in Keylock]
    for λ in 1:Nscales
        sW[:, λ] .= sgs[λ]()
    end
    # Mirror the surrogate coefficients [last part of step (ii) in Keylock]
    sWmirr .= reverse(sW, dims = 1)

    # In the original paper, surrogates and mirror images are matched to original
    # detail coefficients in a circular manner until some error criterion is
    # minimized. Then, the surrogate or its mirror image, depending on which provides
    # the best fit to the original coefficients, is chosen as the representative
    # for a particular dyadic scale. Here, we instead use maximal correlation as
    # the criterion for matching.
    optimal_shifts = zeros(Int, Nscales)
    optimal_shifts_mirr = zeros(Int, Nscales)
    maxcorrs = zeros(Nscales)
    maxcorrs_mirr = zeros(Nscales)

    for i in 0:L-1
        circshift!(circshifted_s, sW, (i, 0))
        circshift!(circshifted_smirr, sWmirr, (i, 0))

        for λ in 1:Nscales
            origW = W[:, λ]
            c = cor(origW, circshifted_s[:, λ])
            if c > maxcorrs[λ]
                maxcorrs[λ] = c
                optimal_shifts[λ] = i
            end

            c_mirr = cor(origW, circshifted_smirr[:, λ])
            if c_mirr > maxcorrs_mirr[λ]
                maxcorrs_mirr[λ] = c_mirr
                optimal_shifts_mirr[λ] = i
            end
        end
    end

    # Decide which coefficients are retained (either surrogate or mirror surrogate coefficients)
    R = zeros(size(W))
    for λ in 1:Nscales
        if maxcorrs[λ] >= maxcorrs_mirr[λ]
            R[:, λ] = circshift(sW[:, λ], optimal_shifts[λ])
        else
            R[:, λ] = circshift(sWmirr[:, λ], optimal_shifts_mirr[λ])
        end
    end

    s = imodwt(R, wl)

    if sg.method.rescale
        s[sortperm(s)] = x_sorted
    end

    return s

end


text:
# basic_maps.jl

"""
The identity map.
"""
struct IdentityMap{T} <: AbstractMap{T,T}
end

(m::IdentityMap)(x) = applymap(m, x)

applymap(map::IdentityMap{T}, x::T) where T = x

apply_inverse(map::IdentityMap{T}, y::T) where T = y

inv(m::IdentityMap) = m

islinear(::IdentityMap) = true

jacobian(m::IdentityMap{T}) where T = UnityMap{T,jac_type(T,T)}()

"The supertype of several maps that map to a constant value."
abstract type AbstractConstantMap{S,T} <: AbstractMap{S,T}
end

applymap(m::AbstractConstantMap, x) = constant(m)

islinear(::AbstractConstantMap) = true

jacobian(m::AbstractConstantMap{S,T}) where {S,T} = ZeroMap{S,jac_type(S,T)}()

"""
The zero map `f(x) = 0`.
"""
struct ZeroMap{S,T} <: AbstractConstantMap{S,T}
end

constant(m::ZeroMap{S,T}) where {S,T} = zero(T)

"""
The unity map `f(x) = 1`.
"""
struct UnityMap{S,T} <: AbstractConstantMap{S,T}
end

constant(m::UnityMap{S,T}) where {S,T} = one(T)

"""
The constant map `f(x) = c`.
"""
struct ConstantMap{S,T} <: AbstractConstantMap{S,T}
    c   ::  T
end

ConstantMap(c::T) where T = ConstantMap{T,T}(c)

constant(m::ConstantMap) = m.c