text:
export BrickletUVLightUVLightCallbackThreshold
struct BrickletUVLightUVLightCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletUVLightIdentity
struct BrickletUVLightIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletUVLight
"""
Measures UV light
"""
mutable struct BrickletUVLight <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletUVLight(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_uv_light")
        deviceInternal = package.BrickletUVLight(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_uv_light
"""
    $(SIGNATURES)

Returns the UV light intensity of the sensor.
The sensor has already weighted the intensity with the erythemal
action spectrum to get the skin-affecting irradiation.

To get UV index you just have to divide the value by 250. For example, a UV
light intensity of 500 is equivalent to an UV index of 2.

If you want to get the intensity periodically, it is recommended to use the
:cb:`UV Light` callback and set the period with
:func:`Set UV Light Callback Period`.
"""
function get_uv_light(device::BrickletUVLight)
    return device.deviceInternal.get_uv_light()
end

export set_uv_light_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`UV Light` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`UV Light` callback is only triggered if the intensity has changed since
the last triggering.
"""
function set_uv_light_callback_period(device::BrickletUVLight, period)
    device.deviceInternal.set_uv_light_callback_period(period)
end

export get_uv_light_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set UV Light Callback Period`.
"""
function get_uv_light_callback_period(device::BrickletUVLight)
    return device.deviceInternal.get_uv_light_callback_period()
end

export set_uv_light_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`UV Light Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the intensity is *outside* the min and max values"
 "'i'",    "Callback is triggered when the intensity is *inside* the min and max values"
 "'<'",    "Callback is triggered when the intensity is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the intensity is greater than the min value (max is ignored)"
"""
function set_uv_light_callback_threshold(device::BrickletUVLight, option, min, max)
    device.deviceInternal.set_uv_light_callback_threshold(option, min, max)
end

export get_uv_light_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set UV Light Callback Threshold`.
"""
function get_uv_light_callback_threshold(device::BrickletUVLight)
    return device.deviceInternal.get_uv_light_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callbacks

* :cb:`UV Light Reached`,

are triggered, if the thresholds

* :func:`Set UV Light Callback Threshold`,

keep being reached.
"""
function set_debounce_period(device::BrickletUVLight, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletUVLight)
    return device.deviceInternal.get_debounce_period()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletUVLight)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletUVLight, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
LATEXIFY(x) = md"$(latexify(x))"

LATEXIFY(s::AbstractString, x) = md"$(s) $(latexify(x))"

L = LATEXIFY

export LATEXIFY, L

ext: jl


text:
@testset "442.find-all-duplicates-in-an-array.jl" begin
    @test find_duplicates_442([4, 3, 2, 7, 8, 2, 3, 1]) == [2, 3]
end


text:
"""
    VerticalStructure

Vertical memory structure. Decomposes stochastic program by stages.

"""
struct VerticalStructure{N, M, SP <: NTuple{M, AbstractScenarioProblems}} <: AbstractBlockStructure{N}
    decisions::NTuple{N, Decisions}
    first_stage::JuMP.Model
    scenarioproblems::SP

    function VerticalStructure(decisions::NTuple{N, Decisions}, scenarioproblems::NTuple{M,AbstractScenarioProblems}) where {N, M}
        M == N - 1 || error("Inconsistent number of stages $N and number of scenario types $M")
        SP = typeof(scenarioproblems)
        return new{N,M,SP}(decisions, Model(), scenarioproblems)
    end
end

function StochasticStructure(decisions::NTuple{N, Decisions}, scenario_types::ScenarioTypes{M}, instantiation::Union{Vertical, DistributedVertical}) where {N, M}
    scenarioproblems = ntuple(Val(M)) do i
        ScenarioProblems(scenario_types[i], instantiation)
    end
    return VerticalStructure(decisions, scenarioproblems)
end

function StochasticStructure(decisions::NTuple{N, Decisions}, scenarios::NTuple{M, Vector{<:AbstractScenario}}, instantiation::Union{Vertical, DistributedVertical}) where {N, M}
    scenarioproblems = ntuple(Val(M)) do i
        ScenarioProblems(scenarios[i], instantiation)
    end
    return VerticalStructure(decisions, scenarioproblems)
end

# Base overloads #
# ========================== #
function Base.print(io::IO, structure::VerticalStructure{N}) where N
    print(io, "Stage 1\n")
    print(io, "============== \n")
    print(io, structure.first_stage)
    for stage in 2:N
        print(io, "\nStage $stage\n")
        print(io, "============== \n")
        for (id, subproblem) in enumerate(subproblems(structure, stage))
            @printf(io, "Subproblem %d (p = %.2f):\n", id, probability(structure, stage, id))
            print(io, subproblem)
            print(io, "\n")
        end
    end
end
function Base.print(io::IO, structure::VerticalStructure{2})
    print(io, "First-stage \n")
    print(io, "============== \n")
    print(io, structure.first_stage)
    print(io, "\nSecond-stage \n")
    print(io, "============== \n")
    for (id, subproblem) in enumerate(subproblems(structure, 2))
        @printf(io, "Subproblem %d (p = %.2f):\n", id, probability(structure, 2, id))
        print(io, subproblem)
        print(io, "\n")
    end
end
# ========================== #

# MOI #
# ========================== #
function MOI.get(structure::VerticalStructure, attr::MOI.AbstractModelAttribute)
    return MOI.get(backend(structure.first_stage), attr)
end
function MOI.get(structure::VerticalStructure, attr::MOI.AbstractVariableAttribute, index::MOI.VariableIndex)
    return MOI.get(backend(structure.first_stage), attr, index)
end
function MOI.get(structure::VerticalStructure, attr::Type{MOI.VariableIndex}, name::String)
    return MOI.get(backend(structure.first_stage), attr, name)
end
function MOI.get(structure::VerticalStructure, attr::MOI.AbstractConstraintAttribute, cindex::MOI.ConstraintIndex)
    return MOI.get(backend(structure.first_stage), attr, cindex)
end
function MOI.get(structure::VerticalStructure, attr::ScenarioDependentModelAttribute)
    n = num_scenarios(structure, attr.stage)
    1 <= attr.scenario_index <= n || error("Scenario index $attr.scenario_index not in range 1 to $n.")
    return MOI.get(scenarioproblems(structure, attr.stage), attr)
end
function MOI.get(structure::VerticalStructure, attr::ScenarioDependentVariableAttribute, index::MOI.VariableIndex)
    n = num_scenarios(structure, attr.stage)
    1 <= attr.scenario_index <= n || error("Scenario index $attr.scenario_index not in range 1 to $n.")
    return MOI.get(scenarioproblems(structure, attr.stage), attr, index)
end
function MOI.get(structure::VerticalStructure, attr::ScenarioDependentConstraintAttribute, ci::MOI.ConstraintIndex)
    n = num_scenarios(structure, attr.stage)
    1 <= attr.scenario_index <= n || error("Scenario index $attr.scenario_index not in range 1 to $n.")
    return MOI.get(scenarioproblems(structure, attr.stage), attr, ci)
end

function MOI.set(structure::VerticalStructure, attr::MOI.AbstractModelAttribute, value)
    MOI.set(backend(structure.first_stage), attr, value)
    return nothing
end
function MOI.set(structure::VerticalStructure, attr::MOI.Silent, flag)
    # Silence master
    MOI.set(backend(structure.first_stage), attr, flag)
    # Silence subproblems
    MOI.set(scenarioproblems(structure), attr, flag)
    return nothing
end
function MOI.set(structure::VerticalStructure, attr::MOI.AbstractVariableAttribute,
                 index::MOI.VariableIndex, value)
    MOI.set(backend(structure.first_stage), attr, index, value)
    return nothing
end
function MOI.set(structure::VerticalStructure, attr::MOI.AbstractConstraintAttribute,
                 ci::MOI.ConstraintIndex, value)
    MOI.set(backend(structure.first_stage), attr, ci, value)
    return nothing
end
function MOI.set(structure::VerticalStructure, attr::ScenarioDependentModelAttribute, value)
    n = num_scenarios(structure, attr.stage)
    1 <= attr.scenario_index <= n || error("Scenario index $attr.scenario_index not in range 1 to $n.")
    MOI.set(scenarioproblems(subproblem, attr.stage), attr, value)
    return nothing
end
function MOI.set(structure::VerticalStructure, attr::ScenarioDependentVariableAttribute,
                 index::MOI.VariableIndex, value)
    n = num_scenarios(structure, attr.stage)
    1 <= attr.scenario_index <= n || error("Scenario index $attr.scenario_index not in range 1 to $n.")
    MOI.set(scenarioproblems(structure, attr.stage), attr, index, value)
    return nothing
end
function MOI.set(structure::VerticalStructure, attr::ScenarioDependentConstraintAttribute,
                 ci::MOI.ConstraintIndex, value)
    n = num_scenarios(structure, attr.stage)
    1 <= attr.scenario_index <= n || error("Scenario index $attr.scenario_index not in range 1 to $n.")
    MOI.set(scenarioproblems(structure, attr.stage), attr, ci, value)
    return nothing
end

function MOI.is_valid(structure::VerticalStructure, index::MOI.VariableIndex, stage::Integer)
    stage == 1 || error("No scenario index specified.")
    return MOI.is_valid(backend(structure.first_stage), index)
end
function MOI.is_valid(structure::VerticalStructure, ci::MOI.ConstraintIndex, stage::Integer)
    stage == 1 || error("No scenario index specified.")
    return MOI.is_valid(backend(structure.first_stage), ci)
end
function MOI.is_valid(structure::VerticalStructure{N}, ci::MOI.ConstraintIndex, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    return MOI.is_valid(scenarioproblems(structure, stage), ci, scenario_index)
end

function MOI.add_constraint(structure::VerticalStructure, f::SingleDecision, s::MOI.AbstractSet)
    return MOI.add_constraint(backend(structure.first_stage), f, s)
end
function MOI.add_constraint(structure::VerticalStructure, f::SingleDecision, s::MOI.AbstractSet, stage::Integer, scenario_index::Integer)
    stage == 1 && error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    return MOI.add_constraint(scenarioproblems(structure, stage), f, s, scenario_index)
end

function MOI.delete(structure::VerticalStructure{N}, index::MOI.VariableIndex, stage::Integer) where N
    stage == 1 || error("No scenario index specified.")
    JuMP.delete(structure.first_stage, DecisionRef(structure.first_stage, index))
    # Remove known decisions
    for s in 2:N
        JuMP.delete(scenarioproblems(structure, s), index)
    end
    return nothing
end
function MOI.delete(structure::VerticalStructure{N}, indices::Vector{MOI.VariableIndex}, stage::Integer) where N
    stage == 1 || error("No scenario index specified.")
    JuMP.delete(structure.first_stage, DecisionRef.(structure.first_stage, indices))
    # Remove known decisions
    for s in 2:N
        JuMP.delete(scenarioproblems(structure, s), indices)
    end
    return nothing
end
function MOI.delete(structure::VerticalStructure, ci::MOI.ConstraintIndex, stage::Integer)
    stage == 1 || error("No scenario index specified.")
    MOI.delete(backend(structure.first_stage), ci)
    return nothing
end
function MOI.delete(structure::VerticalStructure, cis::Vector{<:MOI.ConstraintIndex}, stage::Integer)
    stage == 1 || error("No scenario index specified.")
    MOI.delete(backend(structure.first_stage), cis)
    return nothing
end
function MOI.delete(structure::VerticalStructure{N}, ci::MOI.ConstraintIndex, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    MOI.delete(scenarioproblems(structure, stage), ci, scenario_index)
    return nothing
end
function MOI.delete(structure::VerticalStructure{N}, cis::Vector{<:MOI.ConstraintIndex}, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    MOI.delete(scenarioproblems(structure, stage), cis, scenario_index)
    return nothing
end

# JuMP #
# ========================== #
function JuMP.fix(structure::VerticalStructure{N}, index::MOI.VariableIndex, stage::Integer, val::Number) where N
    dref = DecisionRef(structure.first_stage, index)
    fix(dref, val)
    update_decision_state!(scenarioproblems(structure), index, Known)
    return nothing
end
function JuMP.fix(structure::VerticalStructure{N}, index::MOI.VariableIndex, stage::Integer, scenario_index::Integer, val::Number) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    fix(scenarioproblems(structure, stage), index, scenario_index, val)
    return nothing
end
function JuMP.unfix(structure::VerticalStructure{N}, index::MOI.VariableIndex, stage::Integer) where N
    dref = DecisionRef(structure.first_stage, index)
    unfix(dref)
    return nothing
end
function JuMP.unfix(structure::VerticalStructure{N}, index::MOI.VariableIndex, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    unfix(scenarioproblems(structure, stage), index, scenario_index)
    return nothing
end
function JuMP.objective_function_type(structure::VerticalStructure)
    return jump_function_type(structure.first_stage,
                              MOI.get(structure, MOI.ObjectiveFunctionType()))
end

function JuMP.objective_function(structure::VerticalStructure, FunType::Type{<:AbstractJuMPScalar})
    MOIFunType = moi_function_type(FunType)
    func = MOI.get(structure.first_stage,
                   MOI.ObjectiveFunction{MOIFunType}())::MOIFunType
    return jump_function(structure.first_stage, func)
end

function JuMP.objective_function(structure::AbstractBlockStructure, proxy::JuMP.Model, stage::Integer, FunType::Type{<:AbstractJuMPScalar})
    if stage == 1
        return objective_function(structure, FunType)
    else
        return objective_function(proxy, FunType)
    end
end

function JuMP._moi_optimizer_index(structure::VerticalStructure, index::VI)
    return decision_index(backend(structure.first_stage), index)
end
function JuMP._moi_optimizer_index(structure::VerticalStructure, ci::CI)
    return decision_index(backend(structure.first_stage), ci)
end

function JuMP.set_objective_coefficient(structure::VerticalStructure{N}, index::VI, var_stage::Integer, stage::Integer, coeff::Real) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    var_stage <= stage || error("Can only modify coefficient in current stage of decision or subsequent stages from where decision is taken.")
    if var_stage == 1
        if stage == 1
            dref = DecisionRef(structure.first_stage, index)
            set_objective_coefficient(structure.first_stage, dref, coeff)
        else
            set_objective_coefficient(scenarioproblems(structure, stage), index, coeff)
        end
    else
        for scenario_index in 1:num_scenarios(structure, stage)
            set_objective_coefficient(scenarioproblems(structure, stage), index, scenario_index, coeff)
        end
    end
    return nothing
end
function JuMP.set_objective_coefficient(structure::VerticalStructure{N}, index::VI, var_stage::Integer, stage::Integer, scenario_index::Integer, coeff::Real) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    if var_stage == 1
        set_objective_coefficient(scenarioproblems(structure, stage), index, scenario_index, coeff)
    else
        set_objective_coefficient(scenarioproblems(structure, stage), index, scenario_index, coeff)
    end
    return nothing
end

function JuMP.set_normalized_coefficient(structure::VerticalStructure,
                                         ci::CI{F,S},
                                         index::VI,
                                         value) where {T, F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}}, S}
    MOI.modify(backend(structure.first_stage), ci,
               DecisionCoefficientChange(index, convert(T, value)))
    return nothing
end
function JuMP.set_normalized_coefficient(structure::VerticalStructure{N},
                                         ci::CI{F,S},
                                         index::VI,
                                         stage::Integer,
                                         scenario_index::Integer,
                                         value) where {N, T, F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}}, S}
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    set_normalized_coefficient(scenarioproblems(structure, stage), ci, index, scenario_index, value)
    return nothing
end

function JuMP.normalized_coefficient(structure::VerticalStructure{N},
                                     ci::CI{F,S},
                                     index::VI,
                                     stage::Integer,
                                     scenario_index::Integer) where {N, T, F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}}, S}
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    return normalized_coefficient(scenarioproblems(structure, stage), ci, index, scenario_index)
end

function JuMP.set_normalized_rhs(structure::VerticalStructure,
                                 ci::CI{F,S},
                                 value) where {T,
                                               F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}},
                                               S <: MOIU.ScalarLinearSet{T}}
    MOI.set(backend(structure.first_stage), MOI.ConstraintSet(), ci,
            S(convert(T, value)))
    return nothing
end
function JuMP.set_normalized_rhs(structure::VerticalStructure{N},
                                 ci::CI{F,S},
                                 stage::Integer,
                                 scenario_index::Integer,
                                 value) where {N,
                                               T,
                                               F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}},
                                               S <: MOIU.ScalarLinearSet{T}}
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    stage > 1 || error("There are no scenarios in the first stage.")
    n = num_scenarios(structure, stage)
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    set_normalized_rhs(scenarioproblems(structure, stage), ci, scenario_index, value)
    return nothing
end

function DecisionRef(proxy::JuMP.Model, structure::VerticalStructure, index::VI)
    return DecisionRef(structure.first_stage, index)
end

function JuMP.jump_function(structure::VerticalStructure{N},
                            proxy::JuMP.Model,
                            stage::Integer,
                            f::MOI.AbstractFunction) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    if stage == 1
        return JuMP.jump_function(structure.first_stage, f)
    else
        return JuMP.jump_function(proxy, f)
    end
end

# Getters #
# ========================== #
function structure_name(structure::VerticalStructure)
    return "Vertical"
end

deferred_first_stage(structure::VerticalStructure, ::Val{1}) = num_variables(first_stage(structure)) == 0

function decision(structure::VerticalStructure{N}, index::MOI.VariableIndex, stage::Integer) where N
    stage == 1 || error("No scenario index specified.")
    return decision(structure.decisions[stage], index)
end
# ========================== #

# Setters
# ========================== #
function update_known_decisions!(structure::VerticalStructure)
    update_known_decisions!(structure.first_stage)
    return nothing
end


text:
# Generative model using the change of variables formula
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using LinearAlgebra, InvertibleNetworks, Flux, Test, Random
import Flux.Optimise.update!
Random.seed!(11)

# Target distribution
function swirl(batchsize; noise=.5f0)
    n = sqrt.(rand(Float32, batchsize)) * 1000f0 * 2f0*pi/360f0
    d1x = -cos.(n).*n + rand(Float32, batchsize) .* noise
    d1y = sin.(n).*n + rand(Float32, batchsize) .* noise
    X = cat(d1x', d1y'; dims=1)
    return reshape(X, 1, 1, 2, batchsize) / norm(X, Inf)
end

# Define network
nx = 1
ny = 1
n_in = 2
n_hidden = 32
batchsize = 10
depth = 2
AN = Array{ActNorm}(undef, depth)
L = Array{CouplingLayerGlow}(undef, depth)
AN0 = Array{ActNorm}(undef, depth)
L0 = Array{CouplingLayerGlow}(undef, depth)

for j=1:depth
    AN[j] = ActNorm(n_in; logdet=true)
    L[j] = CouplingLayerGlow(n_in, n_hidden; k1=1, k2=1, p1=0, p2=0, logdet=true)
    AN0[j] = ActNorm(n_in; logdet=true)
    L0[j] = CouplingLayerGlow(n_in, n_hidden; k1=1, k2=1, p1=0, p2=0, logdet=true)
end

# Forward pass
function forward(AN, L, X)
    logdet = 0f0
    for j=1:depth
        X_, logdet1 = AN[j].forward(X)
        X, logdet2 = L[j].forward(X_)
        logdet += (logdet1 + logdet2)
    end
    return X, logdet
end

# Backward pass
function backward(AN, L, ΔX, X)
    logdet = 0f0
    for j=depth:-1:1
        ΔX_, X_ = L[j].backward(ΔX, X)
        ΔX, X = AN[j].backward(ΔX_, X_)
    end
    return ΔX, X
end

###################################################################################################

# Invertibility
X = swirl(1000)
X0 = swirl(1000)
dX = X - X0
Y = forward(AN, L, X)[1]
X_ = backward(AN, L, Y, Y)[2]
@test isapprox(norm(X_ - X)/norm(X), 0f0; atol=1f-3)

###################################################################################################

# Loss
function loss(AN, L, X)
    Y_, logdet = forward(AN, L, X)
    f = .5f0/batchsize*norm(Y_)^2 - logdet
    ΔX, X = backward(AN, L, 1f0/batchsize*Y_, Y_)
    return f, ΔX, L[1].RB.W1.grad, AN[1].s.grad
end

# Gradient test for input
f0, g = loss(AN, L, X0)[1:2]
h = 0.1f0
maxiter = 4
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)

print("\nGradient test loop unrolling\n")
for j=1:maxiter
    f = loss(AN, L, X0 + h*dX)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(dX, g))
    print(err1[j], "; ", err2[j], "\n")
    global h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test for weights
forward(AN0, L0, X) # initialize parameters
ANini = deepcopy(AN0)
Lini  = deepcopy(L0)
dW = L[1].RB.W1.data - L0[1].RB.W1.data
ds = AN[1].s.data - AN0[1].s.data

f0, gX, gW, gs = loss(AN0, L0, X)
h = 0.1f0
maxiter = 4
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)

print("\nGradient test loop unrolling\n")
for j=1:maxiter
    L0[1].RB.W1.data = Lini[1].RB.W1.data + h*dW
    AN0[1].s.data = ANini[1].s.data + h*ds
    f = loss(AN0, L0, X)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(ds, gs) - h*dot(dW, gW))
    print(err1[j], "; ", err2[j], "\n")
    global h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)


text:
#=
Parallel WE estimation of the probability for the Schlögl chemical reaction network
with X(0)= x₀ to reach X(T) ∈ A, the target set.

This reaction network corresponds to:
A + 2S → 3S
3S → A + 2S
B → S
S → B

Parameters were chosen so that the system is bistable, as in:
https://doi.org/10.1063/1.5017955
=#

using Distributed
using Statistics
using HypothesisTests
using Printf

# addprocs(4);

@everywhere using NearestNeighbors

@everywhere include("schlogl_setup.jl");
@everywhere using WeightedEnsemble

# number of coarse steps in WE
@everywhere n_we_steps = 20;
# number of time steps during mutation step
@everywhere T_coarse = Tmax/n_we_steps;
# number of samples in coarse matrix
n_samples_per_bin = 10^2;
# ensemble size
n_particles = 10^2;

# define bin structure
@everywhere voronoi_pts = [[x] for x in 5:5:100];
B₀ = WeightedEnsemble.Voronoi_to_Bins(voronoi_pts);
@everywhere tree = KDTree(hcat(float.(voronoi_pts)...));

# define bin id mapping
@everywhere bin_id = x-> WeightedEnsemble.Voronoi_bin_id(x,tree);
# define the rebinning function
function rebin!(E, B, t)
    @. E.b = bin_id(E.ξ);
    WeightedEnsemble.update_bin_weights!(B, E);
    E, B
end

@everywhere mutation = x -> integrate_schlogl(x, T_coarse);
@everywhere function mutation!(x)
    y = integrate_schlogl(x, T_coarse)
    @. x = y;
    x
end

Random.seed!(100);
x0_vals = copy(voronoi_pts);
bin0_vals = bin_id.(voronoi_pts);
n_bins = length(B₀);
K̃ = WeightedEnsemble.pbuild_coarse_transition_matrix(mutation!, bin_id, x0_vals,bin0_vals, n_bins, n_samples_per_bin);

# define coarse observable as a bin function
f̃ = f.(voronoi_pts);
_,v²_vectors = WeightedEnsemble.build_coarse_vectors(n_we_steps,K̃,float.(f̃));
v² = (x,t)-> v²_vectors[t+1][bin_id(x)]
# define selection function
selection! = (E, B, t)-> WeightedEnsemble.optimal_allocation!(E, B, v², t)

# set up ensemble
E₀ = WeightedEnsemble.Dirac_to_EnsembleWithBins(x₀, n_particles);
rebin!(E₀, B₀, 0);

# run
E = deepcopy(E₀);
B = deepcopy(B₀);
Random.seed!(200)
WeightedEnsemble.prun_we!(E, B, mutation,selection!, rebin!, n_we_steps);
p_est = f.(E.ξ) ⋅ E.ω
@printf("WE Estimate = %g\n", p_est)


text:
##############################################################################
# predict

# if x happens to be a scalar, do we want this to work?
# predict(M::Model, x::Number,            theta=thetaopt(M))   = [x*theta]
# "predict" could be called "predict_y_from_x"
# for each record, x,y,u,v are always vectors
predict(M::Model, x::Array{Float64, 1}, theta=thetaopt(M))   = [dot(x, theta)]
# following could be defined using rowwise, but efficiency might dictate otherwise
predict(M::Model, X::Array{Float64, 2}, theta=thetaopt(M))   = X*theta

"""`predict_y_from_test(M [, theta])`

Allows you compute embedded predictions (i.e., y values) based on a trained ERM model M on test data. Option to
specify a choice of theta. It defaults to `theta=thetaopt(M)`
"""
predict_y_from_test(M::Model,                   theta=thetaopt(M))   = predict(M, Xtest(M), theta)

"""`predict_y_from_train(M [, theta])`

Allows you compute embedded predictions (i.e., y values)  based on a trained ERM model M on train data. Option to
specify a choice of theta. It defaults to `theta=thetaopt(M)`
"""
predict_y_from_train(M::Model,                  theta=thetaopt(M))   = predict(M, Xtrain(M), theta)

"""`predict_v_from_test(M [, theta])`

Allows you compute unembedded predictions (i.e., in V space)
based on a trained ERM model M on test data. Option to
specify a choice of theta. It defaults to `theta=thetaopt(M)`
"""
predict_v_from_test(M::Model,                   theta=thetaopt(M))   = unembedY(M.S, predict(M, Xtest(M), theta))

"""`predict_v_from_train(M [, theta])`

Allows you compute unembedded predictions (i.e., in V space)
based on a trained ERM model M on train data. Option to
specify a choice of theta. It defaults to `theta=thetaopt(M)`
"""
predict_v_from_train(M::Model,                  theta=thetaopt(M))   = unembedY(M.S, predict(M, Xtrain(M), theta))

"""`predict_y_from_u(M, U [, theta])`

Allows you compute embedded predictions (i.e., y values)
based on a trained ERM model M on one or many raw inputs, `U`. Option to
specify a choice of theta. It defaults to `theta=thetaopt(M)`
"""
predict_y_from_u(M::Model, u::Array{T,1}, theta=thetaopt(M)) where {T<:Any} =  predict(M::Model, embedU(M.S, u), theta)
predict_y_from_u(M::Model, U::Array{T,2}, theta=thetaopt(M)) where {T<:Any} =  rowwise(u -> predict_y_from_u(M, u, theta), U)

"""`predict_y_from_u(M, U [, theta])`

Allows you compute unembedded predictions (i.e., in V space)
based on a trained ERM model M on one or many raw inputs, `U`. Option to
specify a choice of theta. It defaults to `theta=thetaopt(M)`
"""
predict_v_from_u(M::Model, u::Array{T,1}, theta=thetaopt(M)) where {T<:Any} =  unembedY(M.S, predict_y_from_u(M, u, theta))
predict_v_from_u(M::Model, U::Array{T,2}, theta=thetaopt(M)) where {T<:Any} =  rowwise(u -> predict_v_from_u(M, u, theta), U)

function uniquevalues(x::Array{T,1}) where {T<:Any}
    n = length(x)
    valtonum = Dict()
    vals = Any[]
    num = 1
    possiblevalues = sort(collect(Set(x)))
    for val in possiblevalues
        valtonum[val] = num
        push!(vals, val)
        num += 1
    end
    return valtonum, vals
end

uniquevalues(x::Array{T,2}) where {T<:Any} = uniquevalues(x[:])

function confusionx(vhat, v)
    K = length(Set(v))
    C = zeros(K,K)
    n = size(v,1)
    valtonum, vals = uniquevalues(v)
    for k=1:n
        i = convert(Int64, vhat[k])
        j = convert(Int64, v[k])
        C[valtonum[i], valtonum[j]] += 1
    end
    return C
end

"""`confusion_train(M)`

Returns the confusion matrix based on training data
for a trained ERM model"""
function confusion_train(M::Model)
    vhat = predict_v_from_train(M)
    return confusionx(vhat, Vtrain(M))
end

"""`confusion_test(M)`

Returns the confusion matrix based on test data
for a trained ERM model"""
function confusion_test(M::Model)
    vhat = predict_v_from_test(M)
    return confusionx(vhat, Vtest(M))
end


text:
# Single decision #
# ========================== #
const FixingConstraint{T} = CI{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}
mutable struct DecisionBridge{T} <: MOIB.Variable.AbstractBridge
    decision::Decision{T}
    variable::MOI.VariableIndex
    fixing_constraint::FixingConstraint{T}

    function DecisionBridge(decision::Decision{T},
                            variable::MOI.VariableIndex,
                            fixing_constraint::FixingConstraint{T}) where T
        return new{T}(decision, variable, fixing_constraint)
    end
end

function MOIB.Variable.bridge_constrained_variable(::Type{DecisionBridge{T}},
                                                   model::MOI.ModelLike,
                                                   set::SingleDecisionSet{T}) where T
    if set.constraint isa NoSpecifiedConstraint
        variable = MOI.add_variable(model)
    else
        variable, constraint = MOI.add_constrained_variable(model, set.constraint)
    end
    # Check state of decision
    if state(set.decision) == NotTaken
        fixing_constraint = FixingConstraint{T}(0)
    else
        # Decision initially fixed
        fixing_constraint =
            MOI.add_constraint(model,
                               MOI.ScalarAffineFunction{T}([MOI.ScalarAffineTerm(one(T), variable)], zero(T)),
                               MOI.EqualTo(set.decision.value))
    end
    return DecisionBridge(set.decision, variable, fixing_constraint)
end

function MOIB.Variable.supports_constrained_variable(
    ::Type{DecisionBridge{T}}, ::Type{SingleDecisionSet{T}}) where T
    return true
end

function MOIB.added_constrained_variable_types(::Type{<:DecisionBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{<:DecisionBridge})
    return Tuple{DataType, DataType}[]
end

# Attributes, Bridge acting as a model
MOI.get(bridge::DecisionBridge, ::MOI.NumberOfVariables) = 1

function MOI.get(bridge::DecisionBridge, ::MOI.ListOfVariableIndices)
    return [bridge.variable]
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::DecisionBridge{T}) where T
    if bridge.fixing_constraint.value != 0
        # Remove the fixing constraint
        MOI.delete(model, bridge.fixing_constraint)
        bridge.fixing_constraint = FixingConstraint{T}(0)
    end
    MOI.delete(model, bridge.variable)
end

# Attributes, Bridge acting as a constraint
function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet,
                 bridge::DecisionBridge)
    return SingleDecisionSet(1, bridge.decision)
end

function MOI.get(model::MOI.ModelLike, ::MOI.ConstraintPrimal,
                 bridge::DecisionBridge{T}) where T
    if bridge.decision.state == Taken
        return bridge.decision.value
    end
    return MOI.get(model, MOI.VariablePrimal(), bridge.variable)
end

function MOI.get(model::MOI.ModelLike, attr::Union{MOI.VariablePrimal, MOI.VariablePrimalStart},
                 bridge::DecisionBridge{T}) where T
    if bridge.decision.state == Taken
        return bridge.decision.value
    end
    return MOI.get(model, attr, bridge.variable)
end

function MOI.get(model::MOI.ModelLike, ::DecisionIndex,
                 bridge::DecisionBridge{T}) where T
    return bridge.variable
end

function MOI.set(model::MOI.ModelLike, attr::MOI.VariablePrimalStart,
                 bridge::DecisionBridge{T}, value) where T
    MOI.set(model, attr, bridge.variable, value)
end

function MOI.set(model::MOI.ModelLike, attr::MOI.ConstraintSet,
                 bridge::DecisionBridge{T}, new_set::SingleDecisionSet{T}) where T
    bridge.decision = new_set.decision
end

function MOI.modify(model::MOI.ModelLike, bridge::DecisionBridge{T}, change::DecisionStateChange) where T
    # Switch on state transition
    if change.new_state == NotTaken
        if bridge.fixing_constraint.value != 0
            # Remove the fixing constraint
            MOI.delete(model, bridge.fixing_constraint)
            bridge.fixing_constraint = FixingConstraint{T}(0)
        end
    else
        set = MOI.EqualTo(bridge.decision.value)
        if bridge.fixing_constraint.value != 0
            # Update existing
            MOI.set(model.model, MOI.ConstraintSet(), bridge.fixing_constraint, set)
        else
            # Add new fixing constraint
            bridge.fixing_constraint = MOI.add_constraint(model,
                                                          MOI.ScalarAffineFunction{T}([MOI.ScalarAffineTerm(one(T), bridge.variable)], zero(T)),
                                                          set)
        end
    end
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::DecisionBridge{T}, ::KnownValuesChange) where T
    # Update fixing constraint if decision is known
    if state(bridge.decision) == Known
        set = MOI.EqualTo(bridge.decision.value)
        MOI.set(model.model, MOI.ConstraintSet(), bridge.fixing_constraint, set)
    end
    return nothing
end

function MOIB.bridged_function(bridge::DecisionBridge{T}) where T
    # Return mapped variable
    return MOI.ScalarAffineFunction{T}(MOI.SingleVariable(bridge.variable))
end

function MOIB.Variable.unbridged_map(bridge::DecisionBridge, vi::MOI.VariableIndex)
    return (bridge.variable => MOI.SingleVariable(vi),)
end

# Multiple decisions #
# ========================== #
struct DecisionsBridge{T} <: MOIB.Variable.AbstractBridge
    decisions::Vector{Decision{T}}
    variables::Vector{MOI.VariableIndex}
    fixing_constraints::Vector{FixingConstraint{T}}
end

function MOIB.Variable.bridge_constrained_variable(::Type{DecisionsBridge{T}},
                                                   model::MOI.ModelLike,
                                                   set::MultipleDecisionSet{T}) where T
    if set.constraint isa NoSpecifiedConstraint
        variables = MOI.add_variables(model, length(set.decisions))
    else
        variables, constraints = MOI.add_constrained_variables(model, set.constraint)
    end
    fixing_constraints = FixingConstraint{T}[]
    # Check decision states
    for (variable, decision) in zip(variables, set.decisions)
        if state(decision) == NotTaken
            push!(fixing_constraints, FixingConstraint{T}(0))
        else
            # Decision initially fixed
            fixing_constraint =
                MOI.add_constraint(model,
                                   MOI.ScalarAffineFunction{T}([MOI.ScalarAffineTerm(one(T), variable)], zero(T)),
                                   MOI.EqualTo(decision.value))
            push!(fixing_constraints, fixing_constraint)
        end
    end
    return DecisionsBridge(set.decisions, variables, fixing_constraints)
end

function MOIB.Variable.supports_constrained_variable(
    ::Type{DecisionsBridge{T}}, ::Type{MultipleDecisionSet{T}}) where T
    return true
end

function MOIB.added_constrained_variable_types(::Type{<:DecisionsBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{<:DecisionsBridge})
    return Tuple{DataType, DataType}[]
end

# Attributes, Bridge acting as a model
MOI.get(bridge::DecisionsBridge, ::MOI.NumberOfVariables) = length(bridge.variables)

function MOI.get(bridge::DecisionsBridge, ::MOI.ListOfVariableIndices)
    return bridge.variables
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::DecisionsBridge{T}) where T
    for (i,fixing_constraint) in enumerate(bridge.fixing_constraints)
        if fixing_constraint.value != 0
            # Remove the fixing constraint
            MOI.delete(model, fixing_constraint)
            bridge.fixing_constraints[i] = FixingConstraint{T}(0)
        end
    end
    MOI.delete(model, bridge.variables)
end

function MOI.delete(model::MOI.ModelLike, bridge::DecisionsBridge{T}, i::MOIB.Variable.IndexInVector) where T
    if bridge.fixing_constraints[i].value != 0
        # Remove the fixing constraint
        MOI.delete(model, bridge.fixing_constraints[i])
        bridge.fixing_constraints[i] = FixingConstraint{T}(0)
    end
    MOI.delete(model, bridge.variables[i.value])
    deleteat!(bridge.variables, i.value)
end

# Attributes, Bridge acting as a constraint
function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet,
                 bridge::DecisionsBridge)
    return MultipleDecisionSet(bridge.decisions)
end

function MOI.get(model::MOI.ModelLike, ::MOI.ConstraintPrimal,
                 bridge::DecisionsBridge{T}) where T
    return [decision.state == Taken ? decision_value.(bridge.decisions) :
            MOI.get(model,
                    MOI.VariablePrimal(),
                    bridge.variables[i]) for (i,decision) in bridge.decisions]
end

function MOI.get(model::MOI.ModelLike,
                 attr::Union{MOI.VariablePrimal, MOI.VariablePrimalStart},
                 bridge::DecisionsBridge, i::MOIB.Variable.IndexInVector)
    decision = bridge.decisions[i.value]
    if decision.state == Taken
        return decision.value
    else
        return MOI.get(model, attr, bridge.variables[i.value])
    end
end

function MOI.get(model::MOI.ModelLike, ::DecisionIndex,
                 bridge::DecisionBridge{T}, i::MOIB.Variable.IndexInVector) where T
    return bridge.variables[i]
end

function MOI.set(model::MOI.ModelLike, attr::MOI.VariablePrimalStart,
                 bridge::DecisionsBridge{T}, val, i::MOIB.Variable.IndexInVector) where T
    MOI.set(model, attr, bridge.variables[i.value], val)
end

function MOI.set(model::MOI.ModelLike, attr::MOI.ConstraintSet,
                 bridge::DecisionsBridge{T}, new_set::MultipleDecisionSet{T}) where T
    bridge.decisions .= new_set.decisions
end

function MOI.modify(model::MOI.ModelLike, bridge::DecisionsBridge{T}, change::DecisionStateChange) where T
    # Switch on state transition
    if change.new_state == NotTaken
        if bridge.fixing_constraints[change.index].value != 0
            # Remove the fixing constraint
            MOI.delete(model, bridge.fixing_constraints[change.index])
            bridge.fixing_constraints[change.index] = FixingConstraint{T}(0)
        end
    else
        set = MOI.EqualTo(bridge.decisions[change.index].value)
        if bridge.fixing_constraints[change.index].value != 0
            # Update existing
            MOI.set(model.model, MOI.ConstraintSet(), bridge.fixing_constraints[change.index], set)
        else
            # Add new fixing constraint
            bridge.fixing_constraints[change.index] = MOI.add_constraint(model,
                                                                         MOI.ScalarAffineFunction{T}([MOI.ScalarAffineTerm(one(T), bridge.variables[change.index])], zero(T)),
                                                                         set)
        end
    end
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::DecisionsBridge{T}, ::KnownValuesChange) where T
    # Update fixing constraints of known decisions
    for (decision,fixing_constraint) in zip(bridge.decisions, bridge.fixing_constraints)
        if state(decision) == Known
            set = MOI.EqualTo(decision.value)
            MOI.set(model.model, MOI.ConstraintSet(), fixing_constraint, set)
        end
    end
    return nothing
end

function MOIB.bridged_function(bridge::DecisionsBridge{T}, i::MOIB.Variable.IndexInVector) where T
    # Return mapped variable
    return MOI.ScalarAffineFunction{T}(MOI.SingleVariable(bridge.variables[i.value]))
end
function MOIB.Variable.unbridged_map(bridge::DecisionsBridge, vi::MOI.VariableIndex, i::MOIB.Variable.IndexInVector)
    return (bridge.variables[i.value] => MOI.SingleVariable(vi),)
end

# Modifications #
# ========================== #
# To enable certain modifications (decision coefficient changes and decision state updates)
# a few methods must be added to MathOptInterface
function MOIB.is_bridged(b::MOIB.AbstractBridgeOptimizer,
                         change::Union{DecisionCoefficientChange, DecisionMultirowChange})
    return MOIB.is_bridged(b, change.decision)
end

function MOIB.is_bridged(b::MOIB.AbstractBridgeOptimizer,
                         change::Union{DecisionStateChange, KnownValuesChange})
    # These modifications should not be handled by the variable bridges
    return false
end

function MOIB.modify_bridged_change(b::MOIB.AbstractBridgeOptimizer, obj,
                                    change::DecisionCoefficientChange)
    f = MOIB.bridged_variable_function(b, change.decision)
    # Continue modification with mapped variable
    MOI.modify(b, obj, DecisionCoefficientChange(only(f.terms).variable_index, change.new_coefficient))
    return nothing
end

function MOIB.modify_bridged_change(b::MOIB.AbstractBridgeOptimizer, obj,
                                    change::DecisionMultirowChange)
    f = MOIB.bridged_variable_function(b, change.decision)
    # Continue modification with mapped variable
    MOI.modify(b, obj, DecisionMultirowChange(only(f.terms).variable_index, change.new_coefficients))
    return nothing
end


text:
# Continuous quadratic problems

function qp1test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # homogeneous quadratic objective
    # Min x^2 + xy + y^2 + yz + z^2
    # st  x + 2y + 3z >= 4 (c1)
    #     x +  y      >= 1 (c2)
    #     x,y \in R

    @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
    MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())
    MOI.supports_constraint(model, MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64})

    MOI.empty!(model)
    @test MOI.is_empty(model)

    v = MOI.add_variables(model, 3)
    @test MOI.get(model, MOI.NumberOfVariables()) == 3

    cf1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,2.0,3.0], v), 0.0)
    c1 = MOI.add_constraint(model, cf1, MOI.GreaterThan(4.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}()) == 1

    c2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,1.0], [v[1],v[2]]), 0.0), MOI.GreaterThan(1.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}()) == 2

    obj = MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Float64}[], MOI.ScalarQuadraticTerm.([2.0, 1.0, 2.0, 1.0, 2.0], v[[1,1,2,2,3]], v[[1,2,2,3,3]]), 0.0)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), obj)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE

    if config.query
        @test obj ≈ MOI.get(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())

        @test cf1 ≈ MOI.get(model, MOI.ConstraintFunction(), c1)

        @test MOI.GreaterThan(4.0) == MOI.get(model, MOI.ConstraintSet(), c1)
    end

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 13/7 atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [4/7,3/7,6/7] atol=atol rtol=rtol
    end
end

function qp2test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # same as QP0 but with duplicate terms
    # then change the objective and sense
    # simple quadratic objective
    # Min x^2 + xy + y^2 + yz + z^2
    # st  x + 2y + 3z >= 4 (c1)
    #     x +  y      >= 1 (c2)
    #     x,y \in R

    @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())
    @test MOI.supports_constraint(model, MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64})

    MOI.empty!(model)
    @test MOI.is_empty(model)

    v = MOI.add_variables(model, 3)
    @test MOI.get(model, MOI.NumberOfVariables()) == 3

    c1f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,2.0,3.0], v), 0.0)
    c1 = MOI.add_constraint(model, c1f, MOI.GreaterThan(4.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}()) == 1

    c2 = MOI.add_constraint(model, MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,1.0], [v[1],v[2]]), 0.0), MOI.GreaterThan(1.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}()) == 2

    obj = MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm.(0.0, v), MOI.ScalarQuadraticTerm.([2.0, 0.5, 0.5, 2.0, 1.0, 1.0, 1.0], [v[1], v[1], v[1], v[2], v[2], v[3], v[3]], [v[1], v[2], v[2], v[2], v[3], v[3], v[3]]), 0.0)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), obj)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE

    if config.query
        @test obj ≈ MOI.get(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())

        @test c1f ≈ MOI.get(model, MOI.ConstraintFunction(), c1)

        @test MOI.GreaterThan(4.0) == MOI.get(model, MOI.ConstraintSet(), c1)
    end

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 13/7 atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [4/7,3/7,6/7] atol=atol rtol=rtol
    end

    # change objective to Max -2(x^2 + xy + y^2 + yz + z^2)
    obj2 = MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm.(0.0, v), MOI.ScalarQuadraticTerm.([-4.0, -1.0, -1.0, -4.0, -2.0, -2.0, -2.0], [v[1], v[1], v[1], v[2], v[2], v[3], v[3]], [v[1], v[2], v[2], v[2], v[3], v[3], v[3]]), 0.0)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), obj2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE

    if config.query
        @test obj2 ≈ MOI.get(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())
    end

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ -2*13/7 atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), v) ≈ [4/7,3/7,6/7] atol=atol rtol=rtol
    end
end

function qp3test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # non-homogeneous quadratic objective
    #    minimize 2 x^2 + y^2 + xy + x + y + 1
    #       s.t.  x, y >= 0
    #             x + y = 1

    @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())
    MOI.supports_constraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    MOI.supports_constraint(model, MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})

    MOI.empty!(model)
    @test MOI.is_empty(model)

    x = MOI.add_variable(model)
    y = MOI.add_variable(model)

    MOI.add_constraint(model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,1.0], [x,y]), 0.0),
        MOI.EqualTo(1.0)
    )

    vc1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    # We test this after the creation of every `SingleVariable` constraint
    # to ensure a good coverage of corner cases.
    @test vc1.value == x.value
    vc2 = MOI.add_constraint(model, MOI.SingleVariable(y), MOI.GreaterThan(0.0))
    @test vc2.value == y.value

    obj = MOI.ScalarQuadraticFunction(
            MOI.ScalarAffineTerm.([1.0,1.0], [x,y]),
            MOI.ScalarQuadraticTerm.([4.0, 2.0, 1.0], [x,y,x], [x,y,y]),
            1.0
          )
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(), obj)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 2.875 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), [x, y]) ≈ [0.25, 0.75] atol=atol rtol=rtol
        @test MOI.get(model, MOI.ConstraintPrimal(), vc1) ≈ 0.25 atol=atol rtol=rtol
        @test MOI.get(model, MOI.ConstraintPrimal(), vc2) ≈ 0.75 atol=atol rtol=rtol
    end

    # change back to linear
    #        max 2x + y + 1
    #       s.t.  x, y >= 0
    #             x + y = 1
    # (x,y) = (1,0), obj = 3
    objf = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0,1.0], [x,y]), 1.0)
    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objf)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    if config.solve
        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 3.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), [x,y]) ≈ [1.0, 0.0] atol=atol rtol=rtol
        @test MOI.get(model, MOI.ConstraintPrimal(), vc1) ≈ 1.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.ConstraintPrimal(), vc2) ≈ 0.0 atol=atol rtol=rtol
    end
end

const qptests = Dict("qp1" => qp1test,
                     "qp2" => qp2test,
                     "qp3" => qp3test)

@moitestset qp

#=
    Quadratically constrained (convex) programs
=#

function qcp1test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # quadratic constraint
    # Max x  + y
    # st -x  + y >= 0 (c1[1])
    #     x  + y >= 0 (c1[2])
    #     x² + y <= 2 (c2)
    # Optimal solution
    # x = 1/2, y = 7/4

    @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    @test MOI.supports_constraint(model, MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives)
    @test MOI.supports_constraint(model, MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64})

    MOI.empty!(model)
    @test MOI.is_empty(model)

    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2

    c1 = MOI.add_constraint(model, MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1,1,2,2], MOI.ScalarAffineTerm.([-1.0,1.0,1.0,1.0], [x,y,x,y])), [0.0,0.0]), MOI.Nonnegatives(2))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives}()) == 1

    c2f = MOI.ScalarQuadraticFunction([MOI.ScalarAffineTerm(1.0, y)], [MOI.ScalarQuadraticTerm(2.0, x, x)], 0.0)
    c2 = MOI.add_constraint(model, c2f, MOI.LessThan(2.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}}()) == 1

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,1.0], [x,y]), 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE

    if config.query
        @test c2f ≈ MOI.get(model, MOI.ConstraintFunction(), c2)
    end

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 2.25 atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), [x,y]) ≈ [0.5,1.75] atol=atol rtol=rtol

        @test MOI.get(model, MOI.ConstraintPrimal(), c1) ≈ [5/4, 9/4] atol=atol rtol=rtol
        @test MOI.get(model, MOI.ConstraintPrimal(), c2) ≈ 2 atol=atol rtol=rtol
    end

    # try delete quadratic constraint and go back to linear

    # FIXME why is this commented ?
    # MOI.delete(model, c2)
    #
    # MOI.optimize!(model)
    #
    # @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status
    #
    # @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    #
    # @test MOI.get(model, MOI.ObjectiveValue()) ≈ 0.0 atol=atol rtol=rtol
end

function qcp2test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # Max x
    # s.t. x^2 <= 2 (c)

    @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    @test MOI.supports_constraint(model, MOI.ScalarQuadraticFunction{Float64},MOI.LessThan{Float64})

    MOI.empty!(model)
    @test MOI.is_empty(model)

    x = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 1

    cf = MOI.ScalarQuadraticFunction([MOI.ScalarAffineTerm(0.0, x)], [MOI.ScalarQuadraticTerm(2.0, x, x)], 0.0)
    c = MOI.add_constraint(model, cf, MOI.LessThan(2.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}}()) == 1

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE

    if config.query
        @test cf ≈ MOI.get(model, MOI.ConstraintFunction(), c)
    end

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        if config.duals
            @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
        end

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ sqrt(2) atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), x) ≈ sqrt(2) atol=atol rtol=rtol

        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 2 atol=atol rtol=rtol

        # TODO - duals
        # @test MOI.get(model, MOI.ConstraintDual(), c) ≈ 0.5/sqrt(2) atol=atol rtol=rtol
    end
end

function qcp3test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # Min -x
    # s.t. x^2 <= 2

    @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    @test MOI.supports_constraint(model, MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64})

    MOI.empty!(model)
    @test MOI.is_empty(model)

    x = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 1

    cf = MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Float64}[], [MOI.ScalarQuadraticTerm(2.0, x, x)], 0.0)
    c = MOI.add_constraint(model, cf, MOI.LessThan(2.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}}()) == 1

    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)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE

    if config.query
        @test cf ≈ MOI.get(model, MOI.ConstraintFunction(), c)
    end

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        if config.duals
            @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT
        end

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ -sqrt(2) atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), x) ≈ sqrt(2) atol=atol rtol=rtol

        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 2 atol=atol rtol=rtol

        # TODO - duals
        # @test MOI.get(model, MOI.ConstraintDual(), c) ≈ -0.5/sqrt(2) atol=atol rtol=rtol
    end
end

function _qcp4test(model::MOI.ModelLike, config::TestConfig, less_than::Bool)
    atol = config.atol
    rtol = config.rtol
    # Max  x
    # s.t. x^2 + x * y + y^2 <= 3 (c)
    #      y == 1

    @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    quad_set = less_than ? MOI.LessThan(3.0) : MOI.GreaterThan(-3.0)
    @test MOI.supports_constraint(model, MOI.ScalarQuadraticFunction{Float64}, typeof(quad_set))

    MOI.empty!(model)
    @test MOI.is_empty(model)

    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2

    vc = MOI.add_constraint(model, MOI.SingleVariable(y), MOI.EqualTo(1.0))
    @test vc.value == y.value

    cf = MOI.ScalarQuadraticFunction([MOI.ScalarAffineTerm(0.0, x)],
                                     MOI.ScalarQuadraticTerm.([2.0, 1.0, 2.0], [x, x, y], [x, y, y]), 0.0)
    if !less_than
        MOIU.operate!(-, Float64, cf)
    end
    c = MOI.add_constraint(model, cf, quad_set)
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarQuadraticFunction{Float64}, typeof(quad_set)}()) == 1

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE

    if config.query
        @test cf ≈ MOI.get(model, MOI.ConstraintFunction(), c)
        @test quad_set == MOI.get(model, MOI.ConstraintSet(), c)
    end

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 1.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 1.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 1.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ MOI.constant(quad_set) atol=atol rtol=rtol
        @test MOI.get(model, MOI.ConstraintPrimal(), vc) ≈ 1.0 atol=atol rtol=rtol
    end
end

qcp4test(model::MOI.ModelLike, config::TestConfig) = _qcp4test(model, config, true)
qcp5test(model::MOI.ModelLike, config::TestConfig) = _qcp4test(model, config, false)

const qcptests = Dict("qcp1" => qcp1test,
                      "qcp2" => qcp2test,
                      "qcp3" => qcp3test,
                      "qcp4" => qcp4test,
                      "qcp5" => qcp5test)

@moitestset qcp

#=
    Quadratically constrained (non-convex) programs
=#

function ncqcp1test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # Max 2x + y
    # s.t. x * y <= 4 (c)
    #      x, y >= 1

    @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    @test MOI.supports_constraint(model, MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64})

    MOI.empty!(model)
    @test MOI.is_empty(model)

    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2

    vc1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(1.0))
    @test vc1.value == x.value
    vc2 = MOI.add_constraint(model, MOI.SingleVariable(y), MOI.GreaterThan(1.0))
    @test vc2.value == y.value

    cf = MOI.ScalarQuadraticFunction([MOI.ScalarAffineTerm(0.0, x)], [MOI.ScalarQuadraticTerm(1.0, x, y)], 0.0)
    c = MOI.add_constraint(model, cf, MOI.LessThan(4.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}}()) == 1

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], [x, y]), 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE

    if config.query
        @test cf ≈ MOI.get(model, MOI.ConstraintFunction(), c)
    end

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 9.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 4.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 1.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 4.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.ConstraintPrimal(), vc1) ≈ 4.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.ConstraintPrimal(), vc2) ≈ 1.0 atol=atol rtol=rtol
    end
end

function ncqcp2test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # Find x, y
    # s.t. x * y == 4 (c)
    #      x * x == 4 (c2)
    #      x, y >= 0

    @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
    @test MOI.supports_constraint(model, MOI.ScalarQuadraticFunction{Float64}, MOI.EqualTo{Float64})

    MOI.empty!(model)
    @test MOI.is_empty(model)

    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2

    vc1 = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    @test vc1.value == x.value
    vc2 = MOI.add_constraint(model, MOI.SingleVariable(y), MOI.GreaterThan(0.0))
    @test vc2.value == y.value

    cf = MOI.ScalarQuadraticFunction([MOI.ScalarAffineTerm(0.0, x)], [MOI.ScalarQuadraticTerm(1.0, x, y)], 0.0)
    c = MOI.add_constraint(model, cf, MOI.EqualTo(4.0))

    cf2 = MOI.ScalarQuadraticFunction([MOI.ScalarAffineTerm(0.0, x)], [MOI.ScalarQuadraticTerm(2.0, x, x)], 0.0)
    c2 = MOI.add_constraint(model, cf2, MOI.EqualTo(4.0))

    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarQuadraticFunction{Float64}, MOI.EqualTo{Float64}}()) == 2

    if config.query
        @test cf ≈ MOI.get(model, MOI.ConstraintFunction(), c)
        @test cf2 ≈ MOI.get(model, MOI.ConstraintFunction(), c2)
    end

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 2.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 2.0 atol=atol rtol=rtol

        @test MOI.get(model, MOI.ConstraintPrimal(), c) ≈ 4.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.ConstraintPrimal(), c2) ≈ 4.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.ConstraintPrimal(), vc1) ≈ 2.0 atol=atol rtol=rtol
        @test MOI.get(model, MOI.ConstraintPrimal(), vc2) ≈ 2.0 atol=atol rtol=rtol
    end
end

const ncqcptests = Dict("ncqcp1" => ncqcp1test,
                        "ncqcp2" => ncqcp2test)

@moitestset ncqcp

#=
    SOCP
=#

function socp1test(model::MOI.ModelLike, config::TestConfig)
    atol = config.atol
    rtol = config.rtol
    # min t
    # s.t. x + y >= 1 (c1)
    #      x^2 + y^2 <= t^2 (c2)
    #      t >= 0 (bound)

    @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    @test MOI.supports_constraint(model, MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64})
    @test MOI.supports_constraint(model, MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64})
    @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})

    MOI.empty!(model)
    @test MOI.is_empty(model)

    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    t = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 3

    c1f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x,y]), 0.0)
    c1 = MOI.add_constraint(model, c1f, MOI.GreaterThan(1.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}()) == 1

    c2f = MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Float64}[], MOI.ScalarQuadraticTerm.([1.0,1.0,-1.0], [x,y,t], [x,y,t]), 0.0)
    c2 = MOI.add_constraint(model, c2f, MOI.LessThan(0.0))
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}}()) == 1

    bound = MOI.add_constraint(model, MOI.SingleVariable(t), MOI.GreaterThan(0.0))
    @test bound.value == t.value
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable, MOI.GreaterThan{Float64}}()) == 1

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, t)], 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE

    if config.query
        @test c1f ≈ MOI.get(model, MOI.ConstraintFunction(), c1)

        @test c2f ≈ MOI.get(model, MOI.ConstraintFunction(), c2)
    end

    if config.solve
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

        MOI.optimize!(model)

        @test MOI.get(model, MOI.TerminationStatus()) == config.optimal_status

        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ObjectiveValue()) ≈ sqrt(1/2) atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), [x,y,t]) ≈ [0.5,0.5,sqrt(1/2)] atol=atol rtol=rtol

        @test MOI.get(model, MOI.VariablePrimal(), [t,x,y,t]) ≈ [sqrt(1/2),0.5,0.5,sqrt(1/2)] atol=atol rtol=rtol
    end
end

const socptests = Dict("socp1" => socp1test)

@moitestset socp

const contquadratictests = Dict("qp" => qptest,
                                "qcp" => qcptest,
                                "ncqcp" => ncqcptest,
                                "socp" => socptest)

@moitestset contquadratic true


text:
#-----------------------------------------------------------------------
#       BEGIN   atmios.jl
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
#   PSEUDO-EPSTEIN FUNCTION
#-----------------------------------------------------------------------
"""
    pEp( h , hc , B , Nmax )

The function `pEp` implements the pseudo-Epstein line shape function
used to approximate the Epstein function used to construct the
ionospheric electron density profiles.

"""
function pEp( h , hc , B , Nmax )
    return  (Nmax/16.0)*(   1/(1 + ((h - hc)/(2.0*B))^2) +
                            1/(1 + ((h - hc)/(4.0*B))^2)^2 +
                            1/(1 + ((h - hc)/(6.0*B))^2)^3 +
                            1/(1 + ((h - hc)/(7.0*B))^2)^4
                        )^2
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#   LINE SHAPE FUNCTION
#-----------------------------------------------------------------------
"""
    LSF( x , xo , σ )

The function `LSF` implements a line shape function employed in the
perturbation models; it has a faster falloff than the Lorentzian.

"""
function LSF( x , xo , σ )
    return  ( σ^2 / ( σ^2 + (x-xo)^2 ) )*
            ( σ^4 / ( σ^4 + (x-xo)^4 ) )
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#   INDEX OF REFRACTION FOR IONOSPHERE
#-----------------------------------------------------------------------
"""
    Δnios( h::Real , θ::Real , ϕ::Real )

The function `Δnios` provides a crude model for the ionospheric
electron density profile. The variable `h` represents height from the
Earth's surface in units of km.

"""
function Δnios( h::Real , θ::Real , ϕ::Real )
    # This is a crude model which emulates ionospheric effects
    # h is height from surface in km
    tpfl=typeof(h)
    return  (
            pEp(h,300,50,4.02411e-5) +    # F layer
            pEp(h,130,30,1.00603e-5) +    # E layer
            pEp(h,75,5,4.02411e-6)        # D layer
            )
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#   TOTAL INDEX OF REFRACTION
#-----------------------------------------------------------------------
"""
    Δntot( h::Real , θ::Real , ϕ::Real )

The function `Δntot` is the total index of refraction profile, which
is a sum of atmospheric and ionospheric index of refraction profiles.
The variable `h` represents height from the Earth's surface in units of
km.

"""
function Δntot( h::Real , θ::Real , ϕ::Real )
    # h is height from surface in km
    tpfl=typeof(h)
    return  ΔnatmStd( h ) + Δnios( h , θ , ϕ )
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#   PERTURBATION MODELS
#-----------------------------------------------------------------------
"""
    P( h::Real , h0::RealVec , σ::RealVec )

The function `P` constructs the perturbation model, which consists of
a sum of line shape functions `LSF`, with the parameters provided in
the vector arguments `h0` and `σ`.

"""
function P( h::Real , h0::RealVec , σ::RealVec )
    tpfl=typeof(h)
    n = length(h0)
    return sum(LSF(h,h0[i],σ[i]) for i=1:n)
end

#-----------------------------------------------------------------------
#   TOTAL INDEX OF REFRACTION WITH PERTURBATION
#-----------------------------------------------------------------------
"""
    Δntp( h::Real , θ::Real , ϕ::Real , δ1::Real=0.001
               , δ2::Real=0.010 , Patm::Function=h->1.0
               , Pion::Function=h->1.0 )

The function `Δntp` adds the perturbations provided in the perturbation
model to the atmospheric and ionospheric profiles defined earlier.

"""
function Δntp( h::Real , θ::Real , ϕ::Real , δ1::Real=0.001
               , δ2::Real=0.010 , Patm::Function=h->1.0
               , Pion::Function=h->1.0 )
    # h is height from surface in km
    tpfl=typeof(h)

    return  ΔnatmStd( h ) * ( 1 + δ1*Patm(h) ) +
            Δnios(h,θ,ϕ) * ( 1 + δ2*Pion(h) )
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#   COMPUTE INDEX OF REFRACTION (SPHERICAL COORDINATES)
#-----------------------------------------------------------------------
"""
    nIRs( RE::Real , r::Real , θ::Real , ϕ::Real , Δnf::Function )

The function `nIRs` computes the index of refraction in spherical
coordinates from the index of refraction profile function `Δnf` in
spherical coordinates.

"""
function nIRs( RE::Real , r::Real , θ::Real , ϕ::Real , Δnf::Function )
    # RE and r in units of Earth mass
    tpfl=typeof(RE)
    um = tpfl(4.435e-6)                     # Conversion factor to km
    h = (r - RE)*um                         # h in units of km
    return  one(tpfl) + Δnf(h,θ,ϕ)
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#   COMPUTE INDEX OF REFRACTION
#-----------------------------------------------------------------------
"""
    nIR( X::RealVec , Δnf::Function=Δntot )

The function `nIR` computes the index of refraction in Cartesian
coordinates; the profile is anchored to the WGS84 ellipsoid defined by
the `rell` function.

"""
function nIR( X::RealVec , Δnf::Function=Δntot )
    # RE and r in units of Earth mass

    xs = SphericalFromCartesian()(X[2:4])

    return  nIRs( rell(X[2:4]) , norm(X[2:4]) , xs.θ , xs.ϕ
                  , Δnf )
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#       END     atmios.jl
#-----------------------------------------------------------------------


text:
module NamedNamedTuples
export @NamedNamedTuple

import NamedTupleTools
import ConstructionBase

macro NamedNamedTuple(Type)
    quote
      struct $(esc(Type)){T}
          x::T
      end

      Base.getproperty(nnt::$(esc(Type)), p::Symbol) = getproperty(getfield(nnt, :x), p)
      Base.propertynames(nnt::$(esc(Type))) = propertynames(getfield(nnt, :x))
      Base.NamedTuple(nnt::$(esc(Type))) = getfield(nnt, :x)
      Base.convert(::Type{NamedTuple}, nnt::$(esc(Type))) = getfield(nnt, :x)::NamedTuple
      Base.convert(::Type{$(esc(Type))}, nnt::NamedTuple) = $(esc(Type))(nnt)
      Base.merge(nnt::$(esc(Type)), fields::NamedTuple) = convert($(esc(Type)), merge(convert(NamedTuple, nnt), fields))
      Base.merge(fields::NamedTuple, nnt::$(esc(Type))) = convert($(esc(Type)), merge(fields, convert(NamedTuple, nnt)))

      NamedTupleTools.delete(nnt::$(esc(Type)), args...) = convert($(esc(Type)), NamedTupleTools.delete(convert(NamedTuple, nnt), args...))

      function ConstructionBase.setproperties(nnt::$(esc(Type)), patch::NamedTuple)
        $(esc(Type))(ConstructionBase.setproperties(getfield(nnt, :x), patch))
      end
  end
end

end


text:
using Unicode
const TEST_GRAPHEMES = true

myreverse(s) = foldl((a,b) -> Base.string(b, a), graphemes(s); init="")

max_stars_repo_name: koladilip/julia-exercism

max_issues_repo_name: koladilip/julia-exercism

max_forks_repo_name: koladilip/julia-exercism


text:
using Base: depwarn
using Flux: Params

check_decay(opt, decay) = decay == 0 ? opt : Optimiser(opt, InvDecay(decay))

# legacy update rule
updaterule(opt, ps) = () -> update!(opt, ps)

function SGD(params::Union{AbstractArray, Params}, η = 0.1; decay = 0.)
  depwarn("SGD(params) is deprecated; use Descent(η::Float64) instead", :SGD)

  ps = params
  opt = Descent(η)
  opt = check_decay(opt, decay)
  updaterule(opt, ps)
end

function Momentum(params::Union{AbstractArray, Params}, η = 0.01; ρ = 0.9, decay = 0.)
  depwarn("Momentum(params) is deprecated; use Momentum(η::Float64) instead", :Momentum)

  ps = params
  opt = Momentum(η, ρ)
  opt = check_decay(opt, decay)
  updaterule(opt, ps)
end

function Nesterov(params::Union{AbstractArray, Params}, η = 0.001; ρ = 0.9, decay = 0.)
  depwarn("Nesterov(params) is deprecated; use Nesterov(η::Float64) instead", :Nesterov)

  ps = params
  opt = Nesterov(η, ρ)
  opt = check_decay(opt, decay)
  updaterule(opt, ps)
end

function RMSProp(params::Union{AbstractArray, Params}, η = 0.001; ρ = 0.9, decay = 0.)
  depwarn("RMSProp(params) is deprecated; use RMSProp(η::Float64) instead", :RMSProp)

  ps = params
  opt = RMSProp(η, ρ)
  opt = check_decay(opt, decay)
  updaterule(opt, ps)
end

function ADAM(params::Union{AbstractArray, Params}, η = 0.001; β1 = 0.9, β2 = 0.999, decay = 0.)
  depwarn("ADAM(params) is deprecated; use ADAM(η::Float64) instead", :ADAM)

  ps = params
  β = (β1, β2)
  opt = ADAM(η, β)
  opt = check_decay(opt, decay)
  updaterule(opt, ps)
end

function ADAGrad(params::Union{AbstractArray, Params}, η::Float64 = 0.1; decay = 0.)
  depwarn("ADAGrad(params) is deprecated; use ADAGrad(η::Float64) instead", :ADAGrad)

  ps = params
  opt = ADAGrad(η)
  opt = check_decay(opt, decay)
  updaterule(opt, ps)
end

function ADADelta(params::Union{AbstractArray, Params}, ρ::Float64 = 0.9; decay = 0.)
  depwarn("ADADelta(params) is deprecated; use ADADelta(η::Float64) instead", :ADADelta)

  ps = params
  opt = ADADelta(ρ)
  opt = check_decay(opt, decay)
  updaterule(opt, ps)
end

function AdaMax(params::Union{AbstractArray, Params}, η = 0.001; β1 = 0.9, β2 = 0.999, decay = 0.)
  depwarn("AdaMax(params) is deprecated; use AdaMax(η::Float64) instead", :AdaMax)

  ps = params
  β = (β1, β2)
  opt = AdaMax(η, β)
  opt = check_decay(opt, decay)
  updaterule(opt, ps)
end

function AMSGrad(params::Union{AbstractArray, Params}, η = 0.001; β1 = 0.9, β2 = 0.999, decay = 0.)
  depwarn("AMSGrad(params) is deprecated; use AMSGrad(η::Float64) instead", :AMSGrad)

  ps = params
  β = (β1, β2)
  opt = AMSGrad(η, β)
  opt = check_decay(opt, decay)
  updaterule(opt, ps)
end

function NADAM(params::Union{AbstractArray, Params}, η = 0.001; β1 = 0.9, β2 = 0.999, decay = 0.)
  depwarn("NADAM(params) is deprecated; use NADAM(η::Float64) instead", :NADAM)

  ps = params
  β = (β1, β2)
  opt = NADAM(η, β)
  opt = check_decay(opt, decay)
  updaterule(opt, ps)
end

function ADAMW(params::Union{AbstractArray, Params}, η = 0.001; β1 = 0.9, β2 = 0.999, decay = 0.)
  depwarn("ADAMW(params) is deprecated; use ADAMW(η::Float64) instead", :ADAMW)

  ps = params
  β = (β1, β2)
  opt = ADAMW(η, β)
  opt = check_decay(opt, decay)
  decay != 0 && (opt = Optimiser(opt, WeightDecay(decay)))
  updaterule(opt, ps)
end

# Old training loop

struct OldOptimiser
  func
end

update!(opt::OldOptimiser, ps) = opt.func()

# Train function
function train!(loss, data, opt; cb = () -> ())
  depwarn("train!(loss, data, opt) is deprecated; use train!(loss, params, data, opt) instead", :train!)
  train!(loss, (), data, OldOptimiser(opt); cb = cb)
end


text:
using QuantEcon
using PyPlot

# == Plot functions == #

function plot_spectral_density(arma::ARMA; ax=None, show=true)
    (w, spect) = spectral_density(arma, two_pi=false)
    if show
        fig, ax = plt.subplots()
    end
    ax[:set_xlim]([0, pi])
    ax[:set_title]("Spectral density")
    ax[:set_xlabel]("frequency")
    ax[:set_ylabel]("spectrum")
    ax[:semilogy](w, spect, axes=ax, color="blue", lw=2, alpha=0.7)
    if show
        plt.show()
    else
        return ax
    end
end

function plot_autocovariance(arma::ARMA; ax=None, show=true)
    acov = autocovariance(arma)
    n = length(acov)
    if show
        fig, ax = plt.subplots()
    end
    ax[:set_title]("Autocovariance")
    ax[:set_xlim](-0.5, n - 0.5)
    ax[:set_xlabel]("time")
    ax[:set_ylabel]("autocovariance")
    ax[:stem](0:(n-1), acov)
    if show
        plt.show()
    else
        return ax
    end
end

function plot_impulse_response(arma::ARMA; ax=None, show=true)
    psi = impulse_response(arma)
    n = length(psi)
    if show
        fig, ax = plt.subplots()
    end
    ax[:set_title]("Impulse response")
    ax[:set_xlim](-0.5, n - 0.5)
    ax[:set_xlabel]("time")
    ax[:set_ylabel]("response")
    ax[:stem](0:(n-1), psi)
    if show
        plt.show()
    else
        return ax
    end
end

function plot_simulation(arma::ARMA; ax=None, show=true)
    X = simulation(arma)
    n = length(X)
    if show
        fig, ax = plt.subplots()
    end
    ax[:set_title]("Sample path")
    ax[:set_xlim](0.0, n)
    ax[:set_xlabel]("time")
    ax[:set_ylabel]("state space")
    ax[:plot](0:(n-1), X, color="blue", lw=2, alpha=0.7)
    if show
        plt.show()
    else
        return ax
    end
end

function quad_plot(arma::ARMA)
    (num_rows, num_cols) = (2, 2)
    fig, axes = plt.subplots(num_rows, num_cols, figsize=(12, 8))
    plt.subplots_adjust(hspace=0.4)
    plot_functions = [plot_impulse_response,
                      plot_spectral_density,
                      plot_autocovariance,
                      plot_simulation]
    for (plot_func, ax) in zip(plot_functions, reshape(axes, 1, 4))
        plot_func(arma, ax=ax, show=false)
    end
    plt.show()
end


text:
# https://www.juliabox.com/notebook/notebooks/tutorials/intro-to-julia/AutoDiff.ipynb
# Auto Diff

using InteractiveUtils
versioninfo()

# The example is the Babylonian algorithm, known to man for millenia, to compute sqrt(x)
function Babylonian(x; N = 10)
    t = (1+x)/2
    for i = 2:N # For illustration purposes, 10 iterations suffice.
        t=(t+x/t)/2
    end
    t # should be square root of x
end

# check that it works
x=2; @show Babylonian(x), √x
@show Babylonian(π), √π # \pi + <tab> , and \sqrt + <tab>

# and now the derivative, almost by magic
struct D <: Number #D is a function-derivative pair
    f::Tuple{Float64,Float64}# f(x), df/dx(x)
end

import Base: +, -, *, /, convert, promote_rule
+(x::D, y::D) = D(x.f .+ y.f)
-(x::D, y::D) = D(x.f .- y.f)
*(x::D, y::D) = D((x.f[1] .* y.f[1], x.f[2] * y.f[1] + x.f[1] * y.f[2]))
/(x::D, y::D) = D((x.f[1] / y.f[1], (y.f[1] * x.f[2] - x.f[1] * y.f[2]) / y.f[1]^2))

convert(::Type{D}, x::Real) = D((x, zero(x)))
promote_rule(::Type{D}, ::Type{<:Number}) = D

x=2.0; @show Babylonian(D((x,1))),(√x,.5/√x)

square = x-> x * x
x=3.0; @show square(D((x, 1))), x*x, 2x

#= output
Julia Version 0.7.0
Commit a4cb80f3ed (2018-08-08 06:46 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) M-5Y51 CPU @ 1.10GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, broadwell)
(Babylonian(x), √x) = (1.414213562373095, 1.4142135623730951)
(Babylonian(π), √π) = (1.7724538509055159, 1.7724538509055159)
(Babylonian(D((x, 1))), (√x, 0.5 / √x)) = (D((1.414213562373095, 0.35355339059327373)), (1.4142135623730951, 0.35355339059327373))
(square(D((x, 1))), x * x, 2x) = (D((9.0, 6.0)), 9.0, 6.0)
=#

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
"""
    Timer

Holds a count that is incremented in units of 1.
"""
mutable struct EvaluationTimer
    Δt::Any

    function EvaluationTimer()
        new(0.0)
    end
end

Δt(t::EvaluationTimer) = getfield(t, :Δt)

on!(t::EvaluationTimer) = setfield!(t, :Δt, Δt(t) - time())

off!(t::EvaluationTimer) = setfield!(t, :Δt, time() + Δt(t))

Base.getproperty(t::EvaluationTimer, s::Symbol) = @restrict EvaluationTimer

Base.propertynames(t::EvaluationTimer) = ()

"""
    EvaluationCounter

Holds a count that is incremented in units of 1.
"""
mutable struct EvaluationCounter
    evaluations::Int

    function EvaluationCounter()
        new(0)
    end
end

evaluations(c::EvaluationCounter) = getfield(c, :evaluations)

increment!(c::EvaluationCounter, units = 1) = begin
    units > 0 && setfield!(c, :evaluations, evaluations(c) + units)
    return evaluations(c)
end

Base.getproperty(t::EvaluationCounter, s::Symbol) = @restrict EvaluationCounter

Base.propertynames(t::EvaluationCounter) = ()

"""
    BlockOptProfile

Holds a count of the trust region subproblem solves and gHS evaluations.

## Development
```md
- Profile should include allocations, increment! by units of allocations
- Generate fields, i.e. counters, timers from simulation units of algorithm 7.1
```
"""
struct BlockOptProfile
    trs_timer::EvaluationTimer
    trs_counter::EvaluationCounter
    ghs_timer::EvaluationTimer
    ghs_counter::EvaluationCounter

    function BlockOptProfile()
        new(EvaluationTimer(), EvaluationCounter(), EvaluationTimer(), EvaluationCounter())
    end
end

trs_timer(p::BlockOptProfile) = getfield(p, :trs_timer)

trs_counter(p::BlockOptProfile) = getfield(p, :trs_counter)

ghs_timer(p::BlockOptProfile) = getfield(p, :ghs_timer)

ghs_counter(p::BlockOptProfile) = getfield(p, :ghs_counter)

Base.getproperty(p::BlockOptProfile, s::Symbol) = @restrict BlockOptProfile

Base.propertynames(p::BlockOptProfile) = ()

"""
    Weaver

Interfaces with Weave.jl package, to construct `.html`
report containing iteration trace plots.
"""
struct Weaver
    # arguments used in *.jmd files to weave report
    f_vals::Vector
    ∇f_norms::Vector
    Δ_vals::Vector
    p_norms::Vector
    ρ_vals::Vector
    weave_level::WeaveLevel

    # weave stuff
    function Weaver(weave_level::WeaveLevel = ALL)
        return new(
            Vector{Float64}(),
            Vector{Float64}(),
            Vector{Float64}(),
            Vector{Float64}(),
            Vector{Float64}(),
            weave_level,
        )
    end
end

f_vals(w::Weaver) = getfield(w, :f_vals)

∇f_norms(w::Weaver) = getfield(w, :∇f_norms)

Δ_vals(w::Weaver) = getfield(w, :Δ_vals)

p_norms(w::Weaver) = getfield(w, :p_norms)

ρ_vals(w::Weaver) = getfield(w, :ρ_vals)

weave!(w::Weaver, field::Symbol, val::Float64) = append!(getfield(w, field), val)

weave!(w::Weaver, accessor::Function, val::Float64) = append!(accessor(w), val)

weave_level(w::Weaver) = getfield(w, :weave_level)

Base.getproperty(w::Weaver, s::Symbol) = @restrict Weaver

Base.propertynames(w::Weaver) = ()

"""
BlockOptTrace

Composite type responsible for logging, generating Weave
reports, and storing a Profile of the simulation.
"""
struct BlockOptTrace
    model::Model
    driver::Driver
    profile::BlockOptProfile
    weaver::Weaver
    log_level::LogLevel
    io::IO

    function BlockOptTrace(model, driver)
        new(
            model,
            driver,
            BlockOptProfile(),
            Weaver(weave_level(driver)),
            log_level(driver),
            open(joinpath(directory(model), "trace.log"), "a"),
        )
    end

    function BlockOptTrace(args...)
        new(args...)
    end
end

# Forwarded Model Methods
model(t::BlockOptTrace) = getfield(t, :model)

directory(t::BlockOptTrace) = directory(model(t))

# Forwarded Driver Methods
driver(t::BlockOptTrace) = getfield(t, :driver)

# Forwarded Profile Methods
profile(t::BlockOptTrace) = getfield(t, :profile)

trs_timer(t::BlockOptTrace) = trs_timer(profile(t))

trs_counter(t::BlockOptTrace) = trs_counter(profile(t))

ghs_timer(t::BlockOptTrace) = ghs_timer(profile(t))

ghs_counter(t::BlockOptTrace) = ghs_counter(profile(t))

# Forwarded Weaver Methods
weaver(t::BlockOptTrace) = getfield(t, :weaver)

weave!(t::BlockOptTrace, field, val) = weave!(weaver(t), field, val)

f_vals(t::BlockOptTrace) = f_vals(weaver(t))

∇f_norms(t::BlockOptTrace) = ∇f_norms(weaver(t))

Δ_vals(t::BlockOptTrace) = Δ_vals(weaver(t))

p_norms(t::BlockOptTrace) = p_norms(weaver(t))

ρ_vals(t::BlockOptTrace) = ρ_vals(weaver(t))

weave_level(t::BlockOptTrace) = weave_level(weaver(t))

log_level(t::BlockOptTrace) = getfield(t, :log_level)

io(t::BlockOptTrace) = getfield(t, :io)

Base.getproperty(t::BlockOptTrace, s::Symbol) = @restrict BlockOptTrace

Base.propertynames(t::BlockOptTrace) = ()

Base.copy(t::BlockOptTrace) = BlockOptTrace(model(t), driver(t), profile(t), weaver(t), log_level(t), io(t))

# generates logging methods info!, warn!, debug!, error!
for level in (:info, :debug, :warn, :error)
    lower_level_str = String(level)
    upper_level_str = uppercase(lower_level_str)
    upper_level_sym = Symbol(upper_level_str)
    fn = Symbol(lower_level_str * "!")
    label = " [" * upper_level_str * "] "
    @eval function $fn(trace::BlockOptTrace, args...)
        if log_level(trace) <= $upper_level_sym
            let io = io(trace)
                println(io, trunc(now(), Second), $label)
                for (idx, arg) in enumerate(args)
                    idx > 0 && show(io, arg)
                    println(io)
                end
                println(io, "")
                flush(io)
            end
        end
    end
end


text:
# ---
# title: 917. Reverse Only Letters
# id: problem917
# author: zhwang
# date: 2022-02-23
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/reverse-only-letters/description/>
# hidden: true
# ---
#
# Given a string `S`, return the "reversed" string where all characters that are
# not a letter stay in the same place, and all letters reverse their positions.
#
#
#
# **Example 1:**
#
#
#
#     Input: "ab-cd"
#     Output: "dc-ba"
#
#
# **Example 2:**
#
#
#
#     Input: "a-bC-dEf-ghIj"
#     Output: "j-Ih-gfE-dCba"
#
#
# **Example 3:**
#
#
#
#     Input: "Test1ng-Leet=code-Q!"
#     Output: "Qedo1ct-eeLg=ntse-T!"
#
#
#
#
# **Note:**
#
#   1. `S.length <= 100`
#   2. `33 <= S[i].ASCIIcode <= 122`
#   3. `S` doesn't contain `\` or `"`
#
#
## @lc code=start
using LeetCode

## Double Pointers
function reverse_only_letters(s::String)::String
    right, res = length(s) + 1, Char[]
    for i in s
        if isletter(i)
            right = findlast(i->isletter(s[i]), 1:(right - 1))
            push!(res, s[right])
        else
            push!(res, i)
        end
    end
    join(res)
end
## 参考 Indigo
function reverse_only_letters(s::String)
    chars = codeunits(s)[:]
    letters = @view(chars[isletter.(Char.(chars))])
    reverse!(letters)
    return String(chars)
end

## @lc code=end
## @lc test=start
@testset "917.reverse-only-letters.jl" begin
    @test reverse_only_letters("ab-cd") == "dc-ba"
    @test reverse_only_letters("a-bC-dEf-ghIj") == "j-Ih-gfE-dCba"
    @test reverse_only_letters("Test1ng-Leet=code-Q!") == "Qedo1ct-eeLg=ntse-T!"
    @test reverse_only_letters("7_28]") == "7_28]"
end
## @lc test=end

## @lc info=start
# link: [solution to 917](https://leetcode-cn.com/problems/reverse-only-letters/solution/python-shuang-zhi-zhen-by-rex-rs-6wmf/)
## @lc info=end


text:
inpt = """
Alice would gain 2 happiness units by sitting next to Bob.
Alice would gain 26 happiness units by sitting next to Carol.
Alice would lose 82 happiness units by sitting next to David.
Alice would lose 75 happiness units by sitting next to Eric.
Alice would gain 42 happiness units by sitting next to Frank.
Alice would gain 38 happiness units by sitting next to George.
Alice would gain 39 happiness units by sitting next to Mallory.
Bob would gain 40 happiness units by sitting next to Alice.
Bob would lose 61 happiness units by sitting next to Carol.
Bob would lose 15 happiness units by sitting next to David.
Bob would gain 63 happiness units by sitting next to Eric.
Bob would gain 41 happiness units by sitting next to Frank.
Bob would gain 30 happiness units by sitting next to George.
Bob would gain 87 happiness units by sitting next to Mallory.
Carol would lose 35 happiness units by sitting next to Alice.
Carol would lose 99 happiness units by sitting next to Bob.
Carol would lose 51 happiness units by sitting next to David.
Carol would gain 95 happiness units by sitting next to Eric.
Carol would gain 90 happiness units by sitting next to Frank.
Carol would lose 16 happiness units by sitting next to George.
Carol would gain 94 happiness units by sitting next to Mallory.
David would gain 36 happiness units by sitting next to Alice.
David would lose 18 happiness units by sitting next to Bob.
David would lose 65 happiness units by sitting next to Carol.
David would lose 18 happiness units by sitting next to Eric.
David would lose 22 happiness units by sitting next to Frank.
David would gain 2 happiness units by sitting next to George.
David would gain 42 happiness units by sitting next to Mallory.
Eric would lose 65 happiness units by sitting next to Alice.
Eric would gain 24 happiness units by sitting next to Bob.
Eric would gain 100 happiness units by sitting next to Carol.
Eric would gain 51 happiness units by sitting next to David.
Eric would gain 21 happiness units by sitting next to Frank.
Eric would gain 55 happiness units by sitting next to George.
Eric would lose 44 happiness units by sitting next to Mallory.
Frank would lose 48 happiness units by sitting next to Alice.
Frank would gain 91 happiness units by sitting next to Bob.
Frank would gain 8 happiness units by sitting next to Carol.
Frank would lose 66 happiness units by sitting next to David.
Frank would gain 97 happiness units by sitting next to Eric.
Frank would lose 9 happiness units by sitting next to George.
Frank would lose 92 happiness units by sitting next to Mallory.
George would lose 44 happiness units by sitting next to Alice.
George would lose 25 happiness units by sitting next to Bob.
George would gain 17 happiness units by sitting next to Carol.
George would gain 92 happiness units by sitting next to David.
George would lose 92 happiness units by sitting next to Eric.
George would gain 18 happiness units by sitting next to Frank.
George would gain 97 happiness units by sitting next to Mallory.
Mallory would gain 92 happiness units by sitting next to Alice.
Mallory would lose 96 happiness units by sitting next to Bob.
Mallory would lose 51 happiness units by sitting next to Carol.
Mallory would lose 81 happiness units by sitting next to David.
Mallory would gain 31 happiness units by sitting next to Eric.
Mallory would lose 73 happiness units by sitting next to Frank.
Mallory would lose 89 happiness units by sitting next to George."""

names = Vector{ASCIIString}()
for line in split(inpt, "\n")
	items = split(line[1:end-1]) # get rid of the . at the end
	if !(items[1] in names)
		push!(names, items[1])
	end
	if !(items[end] in names)
		push!(names, items[end])
	end
end
push!(names, "Jan")
happiness = zeros(Int64, length(names), length(names))

for line in split(inpt, "\n")
	items = split(line[1:end-1]) # get rid of the . at the end
	happy = parse(Int, items[4])
	if items[3] == "lose"
		happy *= -1
	end
	i = findfirst(names, items[1])
	j = findfirst(names, items[end])
	happiness[i, j] = happy
end

function tableHappiness(arrangement)
	sum = 0
	# first wrap around
	i = findfirst(names, arrangement[1])
	j = findfirst(names, arrangement[end])
	sum += happiness[i,j] + happiness[j,i]
	for a in 1:length(arrangement)-1
		i = findfirst(names, arrangement[a])
		j = findfirst(names, arrangement[a+1])
		sum += happiness[i,j] + happiness[j,i]
	end
	return sum
end

best = (0, names)
for p in permutations(names)
	h = tableHappiness(p)
	if h > best[1]
		best = (h, p)
	end
end
println(best)

max_stars_repo_name: jstrube/adventofcode

max_issues_repo_name: jstrube/adventofcode

max_forks_repo_name: jstrube/adventofcode


text:
struct Ticker
    currency::String
    last::Real
    timestamp::DateTime
    volume::Dict{String,Any}
    bid::Real
    ask::Real
    high::Real
    low::Real
    Ticker(currency::String,
           last::Real,
           timestamp::DateTime,
           volume::Dict{String,Any},
           bid::Real,
           ask::Real,
           high::Real,
           low::Real) = new(currency, last,
                                     timestamp, volume,
                                     bid, ask, high, low)
end


text:
# This file is a part of NeighborhoodApproximationIndex.jl

module NeighborhoodApproximationIndex

using SimilaritySearch
using Dates

#include("knr.jl")
include("delone.jl")

end