text:
export BrickletColorV2Color
struct BrickletColorV2Color
    r::Integer
    g::Integer
    b::Integer
    c::Integer
end

export BrickletColorV2ColorCallbackConfiguration
struct BrickletColorV2ColorCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

export BrickletColorV2IlluminanceCallbackConfiguration
struct BrickletColorV2IlluminanceCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
    option::Char
    min::Integer
    max::Integer
end

export BrickletColorV2ColorTemperatureCallbackConfiguration
struct BrickletColorV2ColorTemperatureCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
    option::Char
    min::Integer
    max::Integer
end

export BrickletColorV2Configuration
struct BrickletColorV2Configuration
    gain::Integer
    integration_time::Integer
end

export BrickletColorV2SPITFPErrorCount
struct BrickletColorV2SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

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

export BrickletColorV2
"""
Measures color (RGB value), illuminance and color temperature
"""
mutable struct BrickletColorV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

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

        return new(ipcon, deviceInternal)
    end
end

export get_color
"""
    $(SIGNATURES)

Returns the measured color of the sensor.

The red (r), green (g), blue (b) and clear (c) colors are measured
with four different photodiodes that are responsive at different
wavelengths:

.. image:: /Images/Bricklets/bricklet_color_wavelength_chart_600.jpg
   :scale: 100 %
   :alt: Chart Responsivity / Wavelength
   :align: center
   :target: ../../_images/Bricklets/bricklet_color_wavelength_chart_600.jpg

If you want to get the color periodically, it is recommended
to use the :cb:`Color` callback and set the period with
:func:`Set Color Callback Configuration`.
"""
function get_color(device::BrickletColorV2)
    return device.deviceInternal.get_color()
end

export set_color_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Color`
callback is triggered periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change within the
period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.
"""
function set_color_callback_configuration(device::BrickletColorV2, period, value_has_to_change)
    device.deviceInternal.set_color_callback_configuration(period, value_has_to_change)
end

export get_color_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by
:func:`Set Color Callback Configuration`.
"""
function get_color_callback_configuration(device::BrickletColorV2)
    return device.deviceInternal.get_color_callback_configuration()
end

export get_illuminance
"""
    $(SIGNATURES)

Returns the illuminance affected by the gain and integration time as
set by :func:`Set Configuration`. To get the illuminance in Lux apply this formula::

 lux = illuminance * 700 / gain / integration_time

To get a correct illuminance measurement make sure that the color
values themselves are not saturated. The color value (R, G or B)
is saturated if it is equal to the maximum value of 65535.
In that case you have to reduce the gain, see :func:`Set Configuration`.

If you want to get the value periodically, it is recommended to use the
:cb:`Illuminance` callback. You can set the callback configuration
with :func:`Set Illuminance Callback Configuration`.
"""
function get_illuminance(device::BrickletColorV2)
    return device.deviceInternal.get_illuminance()
end

export set_illuminance_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Illuminance` callback is triggered
periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change
within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.

It is furthermore possible to constrain the callback with thresholds.

The `option`-parameter together with min/max sets a threshold for the :cb:`Illuminance` callback.

The following options are possible:

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

 "'x'",    "Threshold is turned off"
 "'o'",    "Threshold is triggered when the value is *outside* the min and max values"
 "'i'",    "Threshold is triggered when the value is *inside* or equal to the min and max values"
 "'<'",    "Threshold is triggered when the value is smaller than the min value (max is ignored)"
 "'>'",    "Threshold is triggered when the value is greater than the min value (max is ignored)"

If the option is set to 'x' (threshold turned off) the callback is triggered with the fixed period.
"""
function set_illuminance_callback_configuration(device::BrickletColorV2, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_illuminance_callback_configuration(period, value_has_to_change, option, min, max)
end

export get_illuminance_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Illuminance Callback Configuration`.
"""
function get_illuminance_callback_configuration(device::BrickletColorV2)
    return device.deviceInternal.get_illuminance_callback_configuration()
end

export get_color_temperature
"""
    $(SIGNATURES)

Returns the color temperature.

To get a correct color temperature measurement make sure that the color
values themselves are not saturated. The color value (R, G or B)
is saturated if it is equal to the maximum value of 65535.
In that case you have to reduce the gain, see :func:`Set Configuration`.

If you want to get the value periodically, it is recommended to use the
:cb:`Color Temperature` callback. You can set the callback configuration
with :func:`Set Color Temperature Callback Configuration`.
"""
function get_color_temperature(device::BrickletColorV2)
    return device.deviceInternal.get_color_temperature()
end

export set_color_temperature_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Color Temperature` callback is triggered
periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change
within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.

It is furthermore possible to constrain the callback with thresholds.

The `option`-parameter together with min/max sets a threshold for the :cb:`Color Temperature` callback.

The following options are possible:

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

 "'x'",    "Threshold is turned off"
 "'o'",    "Threshold is triggered when the value is *outside* the min and max values"
 "'i'",    "Threshold is triggered when the value is *inside* or equal to the min and max values"
 "'<'",    "Threshold is triggered when the value is smaller than the min value (max is ignored)"
 "'>'",    "Threshold is triggered when the value is greater than the min value (max is ignored)"

If the option is set to 'x' (threshold turned off) the callback is triggered with the fixed period.
"""
function set_color_temperature_callback_configuration(device::BrickletColorV2, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_color_temperature_callback_configuration(period, value_has_to_change, option, min, max)
end

export get_color_temperature_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Color Temperature Callback Configuration`.
"""
function get_color_temperature_callback_configuration(device::BrickletColorV2)
    return device.deviceInternal.get_color_temperature_callback_configuration()
end

export set_light
"""
    $(SIGNATURES)

Turns the white LED on the Bricklet on/off.
"""
function set_light(device::BrickletColorV2, enable)
    device.deviceInternal.set_light(enable)
end

export get_light
"""
    $(SIGNATURES)

Returns the value as set by :func:`Set Light`.
"""
function get_light(device::BrickletColorV2)
    return device.deviceInternal.get_light()
end

export set_configuration
"""
    $(SIGNATURES)

Sets the configuration of the sensor. Gain and integration time
can be configured this way.

For configuring the gain:

* 0: 1x Gain
* 1: 4x Gain
* 2: 16x Gain
* 3: 60x Gain

For configuring the integration time:

* 0: 2.4ms
* 1: 24ms
* 2: 101ms
* 3: 154ms
* 4: 700ms

Increasing the gain enables the sensor to detect a
color from a higher distance.

The integration time provides a trade-off between conversion time
and accuracy. With a longer integration time the values read will
be more accurate but it will take longer to get the conversion
results.
"""
function set_configuration(device::BrickletColorV2, gain, integration_time)
    device.deviceInternal.set_configuration(gain, integration_time)
end

export get_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Configuration`.
"""
function get_configuration(device::BrickletColorV2)
    return device.deviceInternal.get_configuration()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletColorV2)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletColorV2, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletColorV2)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletColorV2, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletColorV2, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletColorV2, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletColorV2)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletColorV2)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletColorV2)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletColorV2, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletColorV2)
    return device.deviceInternal.read_uid()
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::BrickletColorV2)
    return device.deviceInternal.get_identity()
end

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


text:
abstract type Value end

struct PValue <: Value
    contribs::Int # current numbers of contributors
    expected::Int # total number of expected contributors
    value::Float64
    updated::Int64 # time of evaluation
end

mutable struct SValue <: Value
    value::Float64
    updated::Int64 # time of arrival
    SValue(val::Real) = new(val, nowts())
end

mutable struct LValue <: Value
    contribs::Int
    value::Vector{Float64}
    updated::Int64 # greater updated value of the list of values
    LValue(contribs, values) = new(contribs, values, nowts())
    LValue(values::Vector{<:Real}) = new(1, values, nowts())
end

isready(v::SValue) = !isnan(v.value)

isready(v::LValue) = length(v.value) == v.contribs

mutable struct Formula
    output::String
    inputs::Set{String}
    expr::Union{Expr, Symbol}
    progressive::Bool
    Formula(definition, ast) = new(definition, Set(), ast, false)
end

abstract type AbacoError <: Exception end

"""
Wrong formula definition.

[`formula`](@ref) throws `WrongFormula` when a formula is malformed,
for example:

`formula(abaco, "myformula = x + ")`
"""
struct WrongFormula <: AbacoError
    formula::String
end

"""
Formula evaluation failure.

[`ingest`] throws EvalError when a runtime formula evaluation fails,
for example for a wrong numbers of method args:

    formula(abaco, "div(x,y,z")
    ingest(abaco, ts, ne, Dict("x"=>10, "y"=>1, "z"=1))

"""
struct EvalError <: AbacoError
    formula::String
    cause::Exception
end

"""
Attempt to get a value with an invalid index.
"""
struct ValueNotFound <: AbacoError
    var::String
    index::Int
end

mutable struct FormulaState
    done::Bool
    f::String
end

"""
Maintains the state of the abaco.

Before adding formulas and values an abaco `MonoContext`
must be initialized by [`abaco_init`](@ref).
"""
mutable struct Snap
    ts::Int64
    vals::Dict{String, SValue} # variable name => value
    outputs::Dict{String, FormulaState} # formula name => state
    Snap(formula) = begin
        outputs = Dict()
        for fname in keys(formula)
            outputs[fname] = FormulaState(false, fname)
        end
        new(0, Dict(), outputs)
    end
end

"""
The settings of snapshots.

Before adding formulas and values the `SnapsSetting`
must be initialized by [`abaco_init`](@ref).
"""
mutable struct SnapsSetting
    handle::Any
    emitone::Bool
    formula::Dict{String, Formula} # formula name => formula
    dependents::Dict{String, Set{String}} # variable => array of formula names where variable is used
    oncomplete::Union{Function, Nothing}
    SnapsSetting(handle,
            emitone,
            oncomplete) = new(handle,
                              emitone,
                              Dict(),
                              Dict(),
                              oncomplete)
end

mutable struct Node
    ne::String
    tag::String
    snap::Dict{Int, Snap}
    currsnap::Int
    Node(ne, tag) = new(ne, tag, Dict(1 => Snap(Dict())), 1)
    Node(ne, tag, ages) = new(ne, tag, Dict(i => Snap(Dict()) for i in 1:ages), 1)
end

"""
The abaco registry.
"""
mutable struct Context
    interval::Int64
    ages::Int64
    cfg::Dict{String, SnapsSetting}  # tag => snaps setting
    node::Dict{String, Node} # ne => node
    origins::Dict{String, Set{Node}} # sys_2 => [ne_21, ne_22]
    target::Dict{String, Tuple{String, String}} # ne_21 => (ne_21 tag, sys_2)
    oncompletedefault::Union{Function, Nothing}
    Context(interval, ages, oncomplete=nothing) = begin
        new(interval, ages, Dict(), Dict(), Dict(), Dict(), oncomplete)
    end
    Context(interval, ages, settings, oncomplete) = begin
        new(interval, ages, settings, Dict(), Dict(), Dict(), oncomplete)
    end
end

ext: jl


text:
import Cairo,Fontconfig
using Graphs
using GraphPlot, Compose
g = SimpleGraph(3);
add_edge!(g, 1, 2);
add_edge!(g, 2, 3);
draw(PNG("./simple.png", 16cm, 16cm), gplot(g))


text:
# automatically generated by the FlatBuffers compiler, do not modify

UnionVector.eval(quote

FlatBuffers.@STRUCT struct Rapunzel
    hair_length::Int32
end
FlatBuffers.@ALIGN(Rapunzel, 4)

Rapunzel(buf::AbstractVector{UInt8}) = FlatBuffers.read(Rapunzel, buf)
Rapunzel(io::IO) = FlatBuffers.deserialize(io, Rapunzel)

end)

max_stars_repo_name: rjkat/flatbuffers-julia

max_issues_repo_name: rjkat/flatbuffers-julia

max_forks_repo_name: rjkat/flatbuffers-julia


text:
const pkglist = [
    "BioSequences",
    "GenomicFeatures",
    "BioAlignments",
    "BioStructures",
    "GeneticVariation",
    "BioServices"
    "BioTools"
]

import Pkg
for dep in pkglist
    Pkg.test(dep, coverage = false)
end


text:
#=
    Used from https://github.com/lanl-ansi/Juniper.jl
=#

"""
moi_wrapper for variables start
"""

"""
MOI variables
"""

"""
Getters for variables
"""
MOI.get(model::Optimizer, ::MOI.NumberOfVariables) = length(model.variable_info)

function MOI.get(model::Optimizer, ::MOI.ListOfVariableIndices)
    return [MOI.VariableIndex(i) for i in 1:length(model.variable_info)]
end

"""
``MOI.add_variable()`` overloads and safety functions
"""
function MOI.add_variable(model::Optimizer)
    push!(model.variable_info, VariableInfo())
    return MOI.VariableIndex(length(model.variable_info))
end

function MOI.add_variables(model::Optimizer, n::Int)
    return [MOI.add_variable(model) for i in 1:n]
end

function check_inbounds(model::Optimizer, vi::VI)
	num_variables = length(model.variable_info)
	if !(1 <= vi.value <= num_variables)
	    @error "Invalid variable index $vi. ($num_variables variables in the model.)"
	end
	return
end

check_inbounds(model::Optimizer, var::SVF) = check_inbounds(model, var.variable)

function check_inbounds(model::Optimizer, aff::SAF)
	for term in aff.terms
	    check_inbounds(model, term.variable_index)
	end
	return
end

function check_inbounds(model::Optimizer, quad::SQF)
	for term in quad.affine_terms
	    check_inbounds(model, term.variable_index)
	end
	for term in quad.quadratic_terms
	    check_inbounds(model, term.variable_index_1)
	    check_inbounds(model, term.variable_index_2)
	end
	return
end

has_upper_bound(model::Optimizer, vi::VI) =
    model.variable_info[vi.value].has_upper_bound

has_lower_bound(model::Optimizer, vi::VI) =
    model.variable_info[vi.value].has_lower_bound

is_fixed(model::Optimizer, vi::VI) =
    model.variable_info[vi.value].is_fixed

"""
Primal-start support for HybridMINLPSolver
"""
MOI.supports(::Optimizer, ::MOI.VariablePrimalStart, ::Type{VI}) = true

function MOI.set(model::Optimizer, ::MOI.VariablePrimalStart, vi::MOI.VariableIndex, value::Union{Real, Nothing})
    check_inbounds(model, vi)
    if value === nothing
        value = 0.0
    end
    model.variable_info[vi.value].start = value
    return
end

"""
MOI.VariableName attribute support
"""
MOI.supports(::Optimizer, ::MOI.VariableName, ::Type{VI}) = true

function MOI.set(model::Optimizer, ::MOI.VariableName, vi::VI, name::String)
	model.variable_info[vi.value].name = name
	return
end

function MOI.get(model::Optimizer, ::MOI.VariableName, vi::VI)
	return model.variable_info[vi.value].name
end

"""
Populating Variable bounds
"""
function MOI.add_constraint(model::Optimizer, v::SVF, lt::MOI.LessThan{Float64})
    vi = v.variable
    check_inbounds(model, vi)
    if isnan(lt.upper)
        @error "Invalid upper bound value $(lt.upper)."
    end
    if has_upper_bound(model, vi)
        @error "Upper bound on variable $vi already exists."
    end
    if is_fixed(model, vi)
        @error "Variable $vi is fixed. Cannot also set upper bound."
    end
    model.variable_info[vi.value].upper_bound = lt.upper
    model.variable_info[vi.value].has_upper_bound = true
    return MOI.ConstraintIndex{SVF, MOI.LessThan{Float64}}(vi.value)
end

function MOI.add_constraint(model::Optimizer, v::SVF, gt::MOI.GreaterThan{Float64})
    vi = v.variable
    check_inbounds(model, vi)
    if isnan(gt.lower)
        @error "Invalid lower bound value $(gt.lower)."
    end
    if has_lower_bound(model, vi)
        @error "Lower bound on variable $vi already exists."
    end
    if is_fixed(model, vi)
        @error "Variable $vi is fixed. Cannot also set lower bound."
    end
    model.variable_info[vi.value].lower_bound = gt.lower
    model.variable_info[vi.value].has_lower_bound = true
    return MOI.ConstraintIndex{SVF, MOI.GreaterThan{Float64}}(vi.value)
end

function MOI.add_constraint(model::Optimizer, v::SVF, eq::MOI.EqualTo{Float64})
    vi = v.variable
    check_inbounds(model, vi)
    if isnan(eq.value)
        @error "Invalid fixed value $(eq.value)."
    end
    if has_lower_bound(model, vi)
        @error "Variable $vi has a lower bound. Cannot be fixed."
    end
    if has_upper_bound(model, vi)
        @error "Variable $vi has an upper bound. Cannot be fixed."
    end
    if is_fixed(model, vi)
        @error "Variable $vi is already fixed."
    end
    model.variable_info[vi.value].lower_bound = eq.value
    model.variable_info[vi.value].upper_bound = eq.value
    model.variable_info[vi.value].is_fixed = true
    return MOI.ConstraintIndex{SVF, MOI.EqualTo{Float64}}(vi.value)
end

"""
moi_wrapper for variables end
"""

"""
moi_wrapper for constraints start
"""

"""
MOI constraints
"""

"""
Single variable bound constraints
"""
MOI.supports_constraint(::Optimizer, ::Type{SVF}, ::Type{MOI.LessThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{SVF}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{SVF}, ::Type{MOI.EqualTo{Float64}}) = true

"""
Binary/Integer variable support
"""
MOI.supports_constraint(::Optimizer, ::Type{SVF}, ::Type{<:VAR_TYPES}) = true

"""
Linear constraints
"""
MOI.supports_constraint(::Optimizer, ::Type{SAF}, ::Type{<:BOUNDS}) = true

"""
Quadratic constraints (scalar i.e., vectorized constraints are not supported)
"""
MOI.supports_constraint(::Optimizer, ::Type{SQF}, ::Type{MOI.LessThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{SQF}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{SQF}, ::Type{MOI.EqualTo{Float64}}) = true

"""
``define_add_constraint`` macro - from Ipopt.jl
"""
macro define_add_constraint(function_type, set_type, array_name)
    quote
        function MOI.add_constraint(model::Optimizer, func::$function_type, set::$set_type)
            check_inbounds(model, func)
            push!(model.$(array_name), (func, set))
            return MOI.ConstraintIndex{$function_type, $set_type}(length(model.$(array_name)))
        end
    end
end

"""
``MOI.add_constraint()`` overloads for all the supported constraint types
"""
@define_add_constraint(SAF, MOI.LessThan{Float64}, linear_le_constraints)
@define_add_constraint(SAF, MOI.GreaterThan{Float64}, linear_ge_constraints)
@define_add_constraint(SAF, MOI.EqualTo{Float64}, linear_eq_constraints)
@define_add_constraint(SQF, MOI.LessThan{Float64}, quadratic_le_constraints)
@define_add_constraint(SQF, MOI.GreaterThan{Float64}, quadratic_ge_constraints)
@define_add_constraint(SQF, MOI.EqualTo{Float64}, quadratic_eq_constraints)

"""
Binary variable support
"""
function MOI.add_constraint(model::Optimizer, v::SVF, ::MOI.ZeroOne)
    vi = v.variable
    check_inbounds(model, vi)
    # the bounds are set using info_array_of_variables
    # according to mlubin the bounds should not be set here as the bounds should stay when
    # the binary constraint is deleted
    model.variable_info[vi.value].is_binary = true

    return MOI.ConstraintIndex{SVF, MOI.ZeroOne}(vi.value)
end

"""
Integer variable support
"""
function MOI.add_constraint(model::Optimizer, v::SVF, ::MOI.Integer)
    vi = v.variable
	check_inbounds(model, vi)
    model.variable_info[vi.value].is_integer = true

    return MOI.ConstraintIndex{SVF, MOI.Integer}(vi.value)
end

"""
ConstraintIndex support
"""
MOI.supports(::Optimizer, ::MOI.ConstraintName, ::Type{CI}) = true

"""
moi_wrapper for constraints end
"""

"""
moi_wrapper for objective start
"""
"""
MOI objective
"""

"""
Supported objective types and sense
"""
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{SVF}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{SAF}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{SQF}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true

"""
set and get function overloads
"""
MOI.get(model::Optimizer, ::MOI.ObjectiveSense) = model.sense

function MOI.set(model::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    model.sense = sense
    return
end

function MOI.set(model::Optimizer, ::MOI.ObjectiveFunction, func::Union{SVF, SAF, SQF})
    check_inbounds(model, func)
    model.objective = func
    return
end

"""
moi_wrapper for objective end
"""

"""
moi_wrapper for nlp start
"""
"""
MOI NLPBlock
"""

MOI.supports(::Optimizer, ::MOI.NLPBlock) = true

function MOI.set(model::Optimizer, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData)
    model.nlp_data = nlp_data
    return
end

"""
moi_wrapper for nlp end
"""

"""
moi_wrapper for results start
"""
# MathOptInterface results

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
	return model.inner.status
end

function MOI.get(model::Optimizer, ::MOI.RawStatusString)
    return string(model.inner.status)
end

function MOI.get(model::Optimizer, ::MOI.PrimalStatus)
    if state_is_optimal(model.inner.status; allow_almost=model.inner.options.allow_almost_solved)
        return MOI.FEASIBLE_POINT
    else
        return MOI.INFEASIBLE_POINT
    end
end

function MOI.get(model::Optimizer, ::MOI.DualStatus)
    if state_is_optimal(model.inner.status; allow_almost=model.inner.options.allow_almost_solved)
        return MOI.FEASIBLE_POINT
    else
        return MOI.INFEASIBLE_POINT
    end
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveValue)
    if model.inner.status == MOI.OPTIMIZE_NOT_CALLED
        @error "optimize! not called"
    end
    return model.inner.objval
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveBound)
    if model.inner.status == MOI.OPTIMIZE_NOT_CALLED
        @error "optimize! not called"
    end
    return model.inner.best_bound
end

function MOI.get(model::Optimizer, ::MOI.RelativeGap)
    if model.inner.status == MOI.OPTIMIZE_NOT_CALLED
        @error "optimize! not called"
    end
    if isnan(model.inner.objval) || isnan(model.inner.best_bound)
        return NaN
    end
    return abs(model.inner.best_bound-model.inner.objval)/abs(model.inner.objval)
end

function MOI.get(model::Optimizer, ::MOI.SolveTime)
    if model.inner.status == MOI.OPTIMIZE_NOT_CALLED
        @error "optimize! not called"
    end
    return model.inner.soltime
end

function MOI.get(model::Optimizer, ::MOI.VariablePrimal, vi::MOI.VariableIndex)
    if model.inner.status == MOI.OPTIMIZE_NOT_CALLED
        @error "optimize! not called"
    end
    check_inbounds(model, vi)
    return model.inner.solution[vi.value]
end

"""
moi_wrapper for results end
"""

max_stars_repo_name: MdMahbubarRahman/Duran

max_issues_repo_name: MdMahbubarRahman/Duran

max_forks_repo_name: MdMahbubarRahman/Duran


text:
using SparseArrays, LinearAlgebra, BenchmarkTools, MKL

nobs = 200
nlat = 30
A = rand(nobs+nlat, nobs+nlat)
F = zeros(nobs, nobs+nlat)
for i = 1:nobs F[i, i] = 1.0 end

@benchmark F*A

pre = similar(F)

@benchmark copyto!(pre, CartesianIndices(pre), A, CartesianIndices(pre))


text:
const ReadjustT = NamedTuple{(:step, :scale, :min, :max, :trgt, :offset),
                             Tuple{Int64, Float64, Float64, Float64, Float64, Int64}}

function named_readjust(p)
    @assert length(p) == 6
    (step=p[1], scale=p[2], min=p[3], max=p[4], trgt=p[5], offset=p[6])
end

compute_δ(p, mcmc_iter) = p.scale/sqrt(max(1.0, mcmc_iter/p.step-p.offset))

function compute_ϵ(ϵ_old, p, a_r, δ, flip=1.0, f=identity, finv=identity)
    ϵ = finv(f(ϵ_old) + flip*(2*(a_r > p.trgt)-1)*δ)
    ϵ = max(min(ϵ,  p.max), p.min)    # trim excessive updates
end

sigmoid(x, a=1.0) = 1.0 / (1.0 + exp(-a*x))
logit(x, a=1.0) = (log(x) - log(1-x))/a


text:
export Obstacle, Disk, Antidot, RandomDisk, Wall, Circular,
InfiniteWall, PeriodicWall, RandomWall, SplitterWall, FiniteWall,
Semicircle, Ellipse
export translate

using InteractiveUtils
#######################################################################################
## Circles
#######################################################################################
"""
    Obstacle{<:AbstractFloat}
Obstacle supertype.
"""
abstract type Obstacle{T<:AbstractFloat} end
eltype(o::Obstacle{T}) where {T} = T

"""
    Circular{T<:AbstractFloat} <: Obstacle{T}
Circular obstacle supertype.
"""
abstract type Circular{T<:AbstractFloat} <: Obstacle{T} end

"""
    Disk{T<:AbstractFloat}  <: Circular{T}
Disk-like obstacle with propagation allowed outside of the circle (immutable type).
### Fields:
* `c::SVector{2,T}` : Center.
* `r::T` : Radius.
* `name::String` : Some name given for user convenience. Defaults to "Disk".
"""
struct Disk{T<:AbstractFloat} <: Circular{T}
    c::SVector{2,T}
    r::T
    name::String
end

function Disk(c::AbstractVector{T}, r::Real, name::String = "Disk") where {T<:Real}
    S = T <: Integer ? Float64 : T
    return Disk{S}(SVector{2,S}(c), convert(S, abs(r)), name)
end
Disk{T}(args...) where {T} = Disk(args...)

"""
    RandomDisk{T<:AbstractFloat} <: Circular{T}
Disk-like obstacle that randomly (and uniformly) reflects colliding particles.
The propagation is allowed outside of the circle.
### Fields:
* `c::SVector{2,T}` : Center.
* `r::T` : Radius.
* `name::String` : Some name given for user convenience. Defaults to "Random disk".
"""
struct RandomDisk{T<:AbstractFloat} <: Circular{T}
    c::SVector{2,T}
    r::T
    name::String
end
function RandomDisk(
    c::AbstractVector{T}, r::Real, name::String = "RandomDisk") where {T<:Real}
    S = T <: Integer ? Float64 : T
    return RandomDisk{S}(SVector{2,S}(c), convert(S, abs(r)), name)
end
RandomDisk{T}(args...) where {T} = RandomDisk(args...)

"""
    Antidot{T<:AbstractFloat} <: Circular{T}
Disk-like obstacle that allows propagation both inside and outside of the disk
(mutable type). Used in ray-splitting billiards.
### Fields:
* `c::SVector{2,T}` : Center.
* `r::T` : Radius.
* `pflag::Bool` : Flag that keeps track of where the particle is currently
  propagating (`pflag` = propagation-flag).
  `true` stands for *outside* the disk, `false` for *inside* the disk.
  Defaults to `true`.
* `name::String` : Name of the obstacle given for user convenience.
  Defaults to "Antidot".
"""
mutable struct Antidot{T<:AbstractFloat} <: Circular{T}
    c::SVector{2,T}
    r::T
    pflag::Bool
    name::String
end
function Antidot(c::AbstractVector{T}, r::Real,
    pflag::Bool = true, name::String = "Antidot") where {T<:Real}
    S = T <: Integer ? Float64 : T
    return Antidot{S}(SVector{2,S}(c), convert(S, abs(r)), pflag, name)
end
Antidot(c, r, name::String = "Antidot") = Antidot(c, r, true, name)
Antidot{T}(args...) where {T} = Antidot(args...)

"""
    Semicircle{T<:AbstractFloat} <: Circular{T}
Obstacle that represents half a circle. Propagation is allowed only inside
the semicircle.
### Fields:
* `c::SVector{2,T}` : Center.
* `r::T` : Radius.
* `facedir::SVector{2,T}` : Direction where the open face of the Semicircle is facing.
* `name::String` : Name of the obstacle given for user convenience.
  Defaults to "Semicircle".
"""
struct Semicircle{T<:AbstractFloat} <: Circular{T}
    c::SVector{2,T}
    r::T
    facedir::SVector{2,T}
    name::String
end
function Semicircle(
    c::AbstractVector{T}, r::Real, facedir, name = "Semicircle") where {T<:Real}
    S = T <: Integer ? Float64 : T
    return Semicircle{S}(
    SV(c), convert(S, abs(r)), SV(normalize(facedir)), name)
end

show(io::IO, w::Circular{T}) where {T} =
print(io, "$(w.name) {$T}\n", "center: $(w.c)\nradius: $(w.r)")

show(io::IO, w::Semicircle{T}) where {T} =
print(io, "$(w.name) {$T}\n", "center: $(w.c)\nradius: $(w.r)\nfacedir: $(w.facedir)")

#######################################################################################
## Walls
#######################################################################################
"""
    Wall{T<:AbstractFloat} <: Obstacle{T}
Wall obstacle supertype.
"""
abstract type Wall{T<:AbstractFloat} <: Obstacle{T} end

"""
    InfiniteWall{T<:AbstractFloat} <: Wall{T}
Wall obstacle imposing specular reflection during collision (immutable type).
Faster than [`FiniteWall`](@ref), meant to be used for convex billiards.
### Fields:
* `sp::SVector{2,T}` : Starting point of the Wall.
* `ep::SVector{2,T}` : Ending point of the Wall.
* `normal::SVector{2,T}` : Normal vector to the wall, pointing to where the
  particle *will come from before a collision* (pointing towards the inside of the
  billiard). The size of the vector is irrelevant
  since it is internally normalized.
* `name::String` : Name of the obstacle, given for user convenience.
  Defaults to "Wall".
"""
struct InfiniteWall{T<:AbstractFloat} <: Wall{T}
    sp::SVector{2,T}
    ep::SVector{2,T}
    normal::SVector{2,T}
    name::String
end
function InfiniteWall(sp::AbstractVector, ep::AbstractVector,
    n::AbstractVector, name::String = "Wall")
    T = eltype(sp)
    n = normalize(n)
    d = dot(n, ep-sp)
    if abs(d) > 10eps(T)
        error("Normal vector is not actually normal to the wall")
    end
    T = eltype(sp) <: Integer ? Float64 : eltype(sp)
    return InfiniteWall{T}(SVector{2,T}(sp), SVector{2,T}(ep), SVector{2,T}(n), name)
end

"""
    FiniteWall{T<:AbstractFloat} <: Wall{T}
Wall obstacle imposing specular reflection during collision (immutable type).
Slower than [`InfiniteWall`](@ref), meant to be used for non-convex billiards.

Giving a `true` value to the field `isdoor` designates this obstacle to be a `Door`.
This is used in [`escapetime`](@ref) function. A `Door` is a obstacle of the
billiard that the particle can escape from, thus enabling calculations
of escape times.

### Fields:
* `sp::SVector{2,T}` : Starting point of the Wall.
* `ep::SVector{2,T}` : Ending point of the Wall.
* `normal::SVector{2,T}` : Normal vector to the wall, pointing to where the
  particle *will come from before a collision* (pointing towards the inside of the
  billiard). The size of the vector is irrelevant
  since it is internally normalized.
* `isdoor::Bool` : Flag of whether this `FiniteWall` instance is a "Door".
  Defaults to `false`.
* `name::String` : Name of the obstacle, given for user convenience.
  Defaults to "Finite Wall".
"""
struct FiniteWall{T<:AbstractFloat} <: Wall{T}
    sp::SVector{2,T}
    ep::SVector{2,T}
    normal::SVector{2,T}
    width::T
    center::SVector{2,T}
    isdoor::Bool
    name::String
end
function FiniteWall(sp::AbstractVector, ep::AbstractVector,
    n::AbstractVector, isdoor::Bool = false, name::String = "Finite Wall")
    T = eltype(sp)
    n = normalize(n)
    d = dot(n, ep-sp)
    if abs(d) > 10eps(T)
        error("Normal vector is not actually normal to the wall: dot = $d")
    end
    T = eltype(sp) <: Integer ? Float64 : eltype(sp)
    w = norm(ep - sp)
    center = @. (ep+sp)/2
    return FiniteWall{T}(SVector{2,T}(sp), SVector{2,T}(ep), SVector{2,T}(n),
    w, SVector{2,T}(center), isdoor, name)
end
FiniteWall(a, b, c, n::String) = FiniteWall(a, b, c, false, n)

isdoor(w) = w.isdoor

"""
    RandomWall{T<:AbstractFloat} <: Wall{T}
Wall obstacle imposing (uniformly) random reflection during collision (immutable type).
### Fields:
* `sp::SVector{2,T}` : Starting point of the Wall.
* `ep::SVector{2,T}` : Ending point of the Wall.
* `normal::SVector{2,T}` : Normal vector to the wall, pointing to where the
  particle *is expected to come from* (pointing towards the inside of the
  billiard).
* `name::String` : Name of the obstacle, given for user convenience.
  Defaults to "Random wall".
"""
struct RandomWall{T<:AbstractFloat} <: Wall{T}
    sp::SVector{2,T}
    ep::SVector{2,T}
    normal::SVector{2,T}
    name::String
end
function RandomWall(sp::AbstractVector, ep::AbstractVector,
    n::AbstractVector, name::String = "Random wall")
    T = eltype(sp) <: Integer ? Float64 : eltype(sp)
    n = normalize(n)
    d = dot(n, ep-sp)
    if abs(d) > 10eps(T)
        error("Normal vector is not actually normal to the wall")
    end
    return RandomWall{T}(SVector{2,T}(sp), SVector{2,T}(ep), SVector{2,T}(n), name)
end

"""
    PeriodicWall{T<:AbstractFloat} <: Wall{T}
Wall obstacle that imposes periodic boundary conditions upon collision (immutable type).
### Fields:
* `sp::SVector{2,T}` : Starting point of the Wall.
* `ep::SVector{2,T}` : Ending point of the Wall.
* `normal::SVector{2,T}` : Normal vector to the wall, pointing to where the
  particle *will come from* (to the inside the billiard).
  The size of the vector is **important**!
  This vector is added to a particle's `pos` during collision. Therefore the
  size of the normal vector must be correctly associated with the size of the
  periodic cell.
* `name::String` : Name of the obstacle, given for user convenience.
  Defaults to "Periodic wall".
"""
struct PeriodicWall{T<:AbstractFloat} <: Wall{T}
    sp::SVector{2,T}
    ep::SVector{2,T}
    normal::SVector{2,T}
    name::String
end
function PeriodicWall(sp::AbstractVector, ep::AbstractVector,
    n::AbstractVector, name::String = "Periodic wall")
    T = eltype(sp)
    d = dot(n, ep-sp)
    if abs(d) > 10eps(T)
        error("Normal vector is not actually normal to the wall")
    end
    T = eltype(sp) <: Integer ? Float64 : eltype(sp)
    return PeriodicWall{T}(SVector{2,T}(sp), SVector{2,T}(ep), SVector{2,T}(n), name)
end

"""
    SplitterWall{T<:AbstractFloat} <: Wall{T}
Wall obstacle imposing allowing for ray-splitting (mutable type).
### Fields:
* `sp::SVector{2,T}` : Starting point of the Wall.
* `ep::SVector{2,T}` : Ending point of the Wall.
* `normal::SVector{2,T}` : Normal vector to the wall, pointing to where the
  particle *will come from before a collision*.
  The size of the vector is irrelevant.
* `pflag::Bool` : Flag that keeps track of where the particle is currently
  propagating (`pflag` = propagation flag).
  `true` is associated with the `normal` vector the wall is
  instantiated with. Defaults to `true`.
* `name::String` : Name of the obstacle, given for user convenience.
  Defaults to "Splitter wall".
"""
mutable struct SplitterWall{T<:AbstractFloat} <: Wall{T}
    sp::SVector{2,T}
    ep::SVector{2,T}
    normal::SVector{2,T}
    pflag::Bool
    name::String
end
function SplitterWall(sp::AbstractVector, ep::AbstractVector,
    normal::AbstractVector, pflag::Bool = true, name::String = "Splitter wall")
    T = eltype(sp)
    n = normalize(normal)
    d = dot(n, ep-sp)
    if abs(d) > 10eps(T)
        error("Normal vector is not actually normal to the wall")
    end
    T = eltype(sp) <: Integer ? Float64 : eltype(sp)
    return SplitterWall{T}(
    SVector{2,T}(sp), SVector{2,T}(ep), SVector{2,T}(n), pflag, name)
end
SplitterWall(sp, ep, n, name::String = "Splitter wall") =
SplitterWall(sp, ep, n, true, name)
#pretty print:
show(io::IO, w::Wall{T}) where {T} = print(io, "$(w.name) {$T}\n",
"start point: $(w.sp)\nend point: $(w.ep)\nnormal vector: $(w.normal)")

"""
    Ellipse{T<:AbstractFloat}  <: Obstacle{T}
Ellipse obstacle that also allows ray-splitting. The ellipse is always oriented
on the x and y axis (although you can make whichever you want the major one).
### Fields:
* `c::SVector{2,T}` : Center.
* `a::T` : x semi-axis.
* `b::T` : y semi-axis.
* `pflag::Bool` : Flag that keeps track of where the particle is currently
  propagating. `true` (default) is associated with being outside the ellipse.
* `name::String` : Some name given for user convenience. Defaults to `"Ellipse"`.

The ellipse equation is given by
```math
\\left(\\frac{x - c[1]}{a} \\right)^2+ \\left(\\frac{y - c[2]}{b}\\right)^2 = 1
```
"""
mutable struct Ellipse{T<:AbstractFloat} <: Obstacle{T}
    c::SVector{2,T}
    a::T
    b::T
    pflag::Bool
    name::String
    arc::T # arclength of one quadrant of the ellipse
end

"""
    ellipse_arclength(θ, e::Ellipse)
Return the arclength of the ellipse that
spans angle `θ` (in normal coordinates, not in the ellipse parameterization).
Expects `θ` to be in `[0, 2π]`.

After properly calculating the
```math
d=b\\,E\\bigl(\\tan^{-1}(a/b\\,\\tan(\\theta))\\,\\big|\\,1-(a/b)^2\\bigr)
```
"""
function ellipse_arclength(θ, e::Ellipse)
    n, θ = divrem(θ, π/2)
    a = e.a; b = e.b
    if a/b > 1.0
        θ = π/2 - θ
        return (n + 1.0)*e.arc - proper_ellipse_arclength(θ, b, a)
    else
        return n*e.arc + proper_ellipse_arclength(θ, a, b)
    end
end

proper_ellipse_arclength(θ, a, b)  = b*Elliptic.E(atan(a*tan(θ)/b), 1.0 - (a/b)^2)

export ellipse_arclength

function Ellipse(c::AbstractVector{T}, a, b, pflag = true,
                 name::String = "Ellipse") where {T<:Real}
    S = T <: Integer ? Float64 : T
    return Ellipse{S}(SVector{2,S}(c),
            convert(S, abs(a)), convert(S, abs(b)), pflag, name,
            proper_ellipse_arclength(π/2, min(a, b), max(a, b))
        )
end
Ellipse{T}(args...) where {T} = Ellipse(args...)

#######################################################################################
## Normal vectors
#######################################################################################
"""
    normalvec(obst::Obstacle, position)
Return the vector normal to the obstacle's boundary at the given position (which is
assumed to be very close to the obstacle's boundary).
"""
@inline normalvec(wall::Wall, pos) = wall.normal
@inline normalvec(w::PeriodicWall, pos) = normalize(w.normal)
@inline normalvec(w::SplitterWall, pos) = w.pflag ? w.normal : -w.normal
@inline normalvec(disk::Circular, pos) = normalize(pos - disk.c)
@inline normalvec(a::Antidot, pos) =
    a.pflag ? normalize(pos - a.c) : -normalize(pos - a.c)
@inline normalvec(d::Semicircle, pos) = normalize(d.c - pos)

@inline function normalvec(e::Ellipse{T}, pos) where {T}
    # from https://www.algebra.com/algebra/homework/
    # Quadratic-relations-and-conic-sections/Tangent-lines-to-an-ellipse.lesson
    x₀, y₀ = pos
    h, k = e.c
    s = e.pflag ? one(T) : -one(T)
    return s*normalize(SV((x₀-h)/(e.a*e.a), (y₀-k)/(e.b*e.b)))
end

#######################################################################################
## Distances
#######################################################################################
"""
    project_to_line(point, c, n)
Project given `point` to line that contains point `c` and has **normal vector** `n`.
"""
@inline function project_to_line(point, c, n)
    posdot = dot(c - point, n)
    intersection = point + posdot*n
end

"""
    distance(p::AbstractParticle, o::Obstacle)
Return the **signed** distance between particle `p` and obstacle `o`, based on
`p.pos`. Positive distance corresponds to the particle being on the *allowed* region
of the `Obstacle`. E.g. for a `Disk`, the distance is positive when the particle is
outside of the disk, negative otherwise.

    distance(p::AbstractParticle, bd::Billiard)
Return minimum `distance(p, obst)` for all `obst` in `bd`.
If the `distance(p, bd)` is negative and `bd` is convex,
this means that the particle is outside the billiard.

**WARNING** : `distance(p, bd)` may give negative values for non-convex
billiards, or billiards that are composed of several connected sub-billiards.

All `distance` functions can also be given a position (vector) instead of a particle.
"""
(distance(p::AbstractParticle{T}, obst::Obstacle{T})::T) where {T} =
distance(p.pos, obst)

@inline function distance(pos::AbstractVector{T}, w::Wall{T})::T where {T}
    v1 = pos - w.sp
    dot(v1, normalvec(w, pos))
end

# no new distance needed for SplitterWall because the `pflag` field
# has the necessary information to give the correct dinstance,
# since the distance is calculated through the normalvec.

@inline distance(pos::AbstractVector{T}, d::Circular{T}) where {T} =
    norm(pos - d.c) - d.r

@inline function distance(
    pos::AbstractVector{T}, a::Antidot{T})::T where {T}
    d = norm(pos - a.c) - a.r
    a.pflag ? d : -d
end

function distance(pos::AbstractVector{T}, s::Semicircle{T}) where {T}
    # Check on which half of circle is the particle
    v1 = pos - s.c
    nn = dot(v1, s.facedir)
    if nn ≤ 0 # I am "inside" semicircle
        return s.r - norm(pos - s.c)
    else # I am on the "other side"
        # end1 = SV(s.c[1] + s.r*s.facedir[2], s.c[2] - s.r*s.facedir[1])
        # end2 = SV(s.c[1] - s.r*s.facedir[2], s.c[2] + s.r*s.facedir[1])
        # return min(norm(pos - end1), norm(pos - end2))
        return one(T)
    end
end

function distance(pos::SV, e::Ellipse{T})::T where {T}
    d = ((pos[1] - e.c[1])/e.a)^2 + ((pos[2]-e.c[2])/e.b)^2 - 1.0
    e.pflag ? d : -d
end

# The entire functionality of `distance_init` is necessary only for
# FiniteWall !!!
distance_init(p::AbstractParticle, a::Obstacle) = distance_init(p.pos, a)
distance_init(pos::SVector, a::Obstacle) = distance(pos, a)

function distance_init(pos::SVector{2,T}, w::FiniteWall{T})::T where {T}

    n = normalvec(w, pos)
    posdot = dot(w.sp - pos, n)
    if posdot ≥ 0 # I am behind wall
        intersection = project_to_line(pos, w.center, n)
        dfc = norm(intersection - w.center)
        if dfc > w.width/2
            return +1.0 # but not directly behind
        else
            return -1.0
        end
    end
    v1 = pos - w.sp
    dot(v1, n)
end

####################################################
## Initial Conditions
####################################################
"""
    cellsize(bd)
Return the delimiters `xmin, ymin, xmax, ymax` of the given obstacle/billiard.

Used in `randominside()`, error checking and plotting.
"""
function cellsize(d::Circular{T}) where {T}
    xmin = ymin = T(Inf)
    xmax = ymax = T(-Inf)
    return xmin, ymin, xmax, ymax
end

function cellsize(w::Wall)
    xmin = min(w.sp[1], w.ep[1])
    xmax = max(w.sp[1], w.ep[1])
    ymin = min(w.sp[2], w.ep[2])
    ymax = max(w.sp[2], w.ep[2])
    return xmin, ymin, xmax, ymax
end

function cellsize(a::Antidot{T}) where {T}
    if a.pflag
        xmin = ymin = T(Inf)
        xmax = ymax = T(-Inf)
    else
        xmin, ymin = a.c - a.r
        xmax, ymax = a.c + a.r
    end
    return xmin, ymin, xmax, ymax
end

function cellsize(e::Ellipse{T}) where {T}
    if e.pflag
        xmin = ymin = T(Inf)
        xmax = ymax = T(-Inf)
    else
        xmin = e.c[1] - e.a; ymin = e.c[2] - e.b
        xmax = e.c[1] + e.a; ymax = e.c[2] + e.b
    end
    return xmin, ymin, xmax, ymax
end

function cellsize(a::Semicircle{T}) where {T}
    xmin, ymin = a.c - a.r
    xmax, ymax = a.c + a.r
    return xmin, ymin, xmax, ymax
end

####################################################
## Translations
####################################################
"""
    translate(obst::Obstacle, vector)
Create a copy of the given obstacle with its position
translated by `vector`.
"""
function translate end

for T in subtypes(Circular)
  @eval translate(d::$T, vec) = ($T)(d.c .+ vec, d.r)
end

for T in subtypes(Wall)
  @eval translate(w::$T, vec) = ($T)(w.sp + vec, w.ep + vec, w.normal)
end

translate(e::Ellipse, vec) = Ellipse(e.c + vec, e.a, e.b)


text:
function part1(input)
    deck1, deck2 = parse_input(input)
    play_combat!(deck1, deck2)
    return max(score_deck(deck1), score_deck(deck2))
end

function part2(input)
    deck1, deck2 = parse_input(input)
    winner = play_recursive_combat!(deck1, deck2)
    return score_deck((deck1, deck2)[winner])
end

function parse_input(input)
    return parse_deck.(split(readchomp(input), "\n\n"))
end

function parse_deck(deck)
    return parse.(Int, split(deck, "\n")[2:end])
end

function play_combat!(deck1, deck2)
    while !isempty(deck1) && !isempty(deck2)
        card1 = popfirst!(deck1)
        card2 = popfirst!(deck2)
        if card1 > card2
            push!(deck1, card1, card2)
        else
            push!(deck2, card2, card1)
        end
    end
end

function play_recursive_combat!(deck1, deck2)
    played_states = Set{Tuple{Vector{Int}, Vector{Int}}}()
    while true
        isempty(deck1) && return 2
        isempty(deck2) && return 1
        (deck1, deck2) in played_states && return 1
        push!(played_states, (copy(deck1), copy(deck2)))
        card1 = popfirst!(deck1)
        card2 = popfirst!(deck2)
        if card1 <= length(deck1) && card2 <= length(deck2)
            winner = play_recursive_combat!(deck1[1:card1], deck2[1:card2])
        else
            winner = card1 > card2 ? 1 : 2
        end
        if winner == 1
            push!(deck1, card1, card2)
        else
            push!(deck2, card2, card1)
        end
    end
end

function score_deck(deck)
    return sum(reverse(deck) .* (1:length(deck)))
end


text:
# Base.eltype
eltype(::Type{BondGraph}) = AbstractNode
eltype(bg::BondGraph) = AbstractNode

# edgetype
lg.edgetype(::Type{BondGraph}) = lg.AbstractSimpleEdge{Integer}
lg.edgetype(bg::BondGraph) = lg.AbstractSimpleEdge{Integer}

# edges
lg.edges(bg::BondGraph) = bg.bonds
lg.ne(bg::BondGraph) = length(bg.bonds)
lg.has_edge(bg::BondGraph, bond::Bond) =  any(b -> b === bond, bg.bonds) # strong equality
lg.has_edge(bg::BondGraph, n1::AbstractNode, n2::AbstractNode) = lg.has_edge(bg, vertex(n1), vertex(n2))
lg.has_edge(bg::BondGraph, s::Int, d::Int) =
    any(b -> lg.src(b) === s && lg.dst(b) === d, bg.bonds)

# vertices
lg.vertices(bg::BondGraph) = vertex.(bg.nodes)
lg.nv(bg::BondGraph) = length(bg.nodes)
lg.has_vertex(bg::BondGraph, node::AbstractNode) = any(n -> n === node, bg.nodes) # strong equality
lg.has_vertex(bg::BondGraph, v::Int) = v <= length(bg.nodes)

# inneighbors, outneighbors
lg.inneighbors(bg::BondGraph, n::AbstractNode) = bg.nodes[lg.inneighbors(bg, vertex(n))]
lg.inneighbors(bg::BondGraph, v::Int) = [lg.src(b) for b in bg.bonds if lg.dst(b) == v]
lg.outneighbors(bg::BondGraph, n::AbstractNode) = bg.nodes[lg.outneighbors(bg, vertex(n))]
lg.outneighbors(bg::BondGraph, v::Int) = [lg.dst(b) for b in bg.bonds if lg.src(b) == v]
lg.all_neighbors(bg::BondGraph, n::AbstractNode) = bg.nodes[lg.all_neighbors(bg, vertex(n))]

# is_directed
lg.is_directed(::Type{BondGraph}) = true
lg.is_directed(bg::BondGraph) = true

# zero
lg.zero(::Type{BondGraph}) = BondGraph()
lg.zero(bg::BondGraph) = BondGraph()

# src, dst
lg.src(b::Bond) = vertex(srcnode(b))
lg.dst(b::Bond) = vertex(dstnode(b))

# weights
# TODO

# Mutations
function lg.add_vertex!(bg::BondGraph, node::AbstractNode)
    lg.has_vertex(bg, node) && return false
    push!(bg.nodes, node)
    set_vertex!(node, lg.nv(bg))
    return true
end

function lg.rem_vertex!(bg::BondGraph, node::AbstractNode)
    lg.has_vertex(bg, node) || return false
    index = vertex(node)
    deleteat!(bg.nodes, index)
    for n in bg.nodes[index:end]
        n.vertex[] -= 1
    end
    return true
end

function lg.add_edge!(bg::BondGraph, srcport::Port, dstport::Port)
    srcnode = srcport.node
    dstnode = dstport.node
    lg.has_edge(bg, srcnode, dstnode) && error("Bond already exists between $srcnode and $dstnode")
    lg.has_edge(bg, dstnode, srcnode) && error("Bond already exists between $dstnode and $srcnode")
    new_bond = Bond(srcport, dstport)
    push!(bg.bonds, new_bond)
    updateport!(srcnode, srcport.index)
    updateport!(dstnode, dstport.index)
    return new_bond
end
function lg.add_edge!(bg::BondGraph, srcnode::AbstractNode, dstnode::AbstractNode)
    lg.add_edge!(bg, Port(srcnode), Port(dstnode))
end

function lg.rem_edge!(bg::BondGraph, srcnode::AbstractNode, dstnode::AbstractNode)
    lg.has_edge(bg, srcnode, dstnode) || lg.has_edge(bg, dstnode, srcnode) || return
    index = findfirst(b -> srcnode in b && dstnode in b, bg.bonds)
    deleted_bond = bg.bonds[index]
    deleteat!(bg.bonds, index)
    updateport!(srcnode, deleted_bond.src.index)
    updateport!(dstnode, deleted_bond.dst.index)
    return deleted_bond
end

max_stars_repo_name: jedforrest/BondGraphs

max_issues_repo_name: jedforrest/BondGraphs

max_forks_repo_name: jedforrest/BondGraphs


text:
module Elliptic
# elliptic integrals of 1st/2nd/3rd kind
export E, F, K, Pi

# jacobi elliptic functions
export Jacobi

# matlab compatible
export ellipj, ellipke

include("jacobi.jl")
include("slatec.jl")

function E(phi::Float64, m::Float64)
    if isnan(phi) || isnan(m) return NaN end
    if m < 0. || m > 1. throw(DomainError(m, "argument m not in [0,1]")) end
    if abs(phi) > pi/2
        phi2 = phi + pi/2
        return 2*fld(phi2,pi)*E(m) - _E(cos(mod(phi2,pi)), m)
    end
    _E(sin(phi), m)
end
function _E(sinphi::Float64, m::Float64)
    sinphi2 = sinphi^2
    cosphi2 = 1. - sinphi2
    y = 1. - m*sinphi2
    drf,ierr1 = SLATEC.DRF(cosphi2, y, 1.)
    drd,ierr2 = SLATEC.DRD(cosphi2, y, 1.)
    if ierr1 == ierr2 == 0
        return sinphi*(drf - m*sinphi2*drd/3)
    elseif ierr1 == ierr2 == 2
        # 2 - (1+m)*sinphi2 < tol
        return sinphi
    end
    NaN
end
E(phi::Real, m::Real) = E(Float64(phi), Float64(m))

"""
`ellipke(m::Real)`
returns `(K(m), E(m))` for scalar `0 ≤ m ≤ 1`
"""
function ellipke(m::Float64)
    if m < 1.
        y = 1. - m
        drf,ierr1 = SLATEC.DRF(0., y, 1.)
        drd,ierr2 = SLATEC.DRD(0., y, 1.)
        @assert ierr1 == 0 && ierr2 == 0
        return drf, drf - m*drd/3
    elseif m == 1.
        return Inf, 1.
    elseif isnan(m)
        return NaN, NaN
    else
        throw(DomainError(m, "argument m not <= 1"))
    end
end
ellipke(x::Real) = ellipke(Float64(x))

E(m::Float64) = ellipke(m)[2]
E(x::Float32) = Float32(E(Float64(x)))
E(x::Real) = E(Float64(x))

# assumes 0 ≤ m ≤ 1
function rawF(sinphi::Float64, m::Float64)
    if abs(sinphi) == 1. && m == 1. return sign(sinphi)*Inf end
    sinphi2 = sinphi^2
    drf,ierr = SLATEC.DRF(1. - sinphi2, 1. - m*sinphi2, 1.)
    @assert ierr == 0
    sinphi*drf
end

function F(phi::Float64, m::Float64)
    if isnan(phi) || isnan(m) return NaN end
    if m < 0. || m > 1. throw(DomainError(m, "argument m not in [0,1]")) end
    if abs(phi) > pi/2
        # Abramowitz & Stegun (17.4.3)
        phi2 = phi + pi/2
        return 2*fld(phi2,pi)*K(m) - rawF(cos(mod(phi2,pi)), m)
    end
    rawF(sin(phi), m)
end
F(phi::Real, m::Real) = F(Float64(phi), Float64(m))

function K(m::Float64)
    if m < 1.
        drf,ierr = SLATEC.DRF(0., 1. - m, 1.)
        @assert ierr == 0
        return drf
    elseif m == 1.
        return Inf
    elseif isnan(m)
        return NaN
    else
        throw(DomainError(m, "argument m not <= 1"))
    end
end
K(x::Float32) = Float32(K(Float64(x)))
K(x::Real) = K(Float64(x))

function Pi(n::Float64, phi::Float64, m::Float64)
    if isnan(n) || isnan(phi) || isnan(m) return NaN end
    if m < 0. || m > 1. throw(DomainError(m, "argument m not in [0,1]")) end
    sinphi = sin(phi)
    sinphi2 = sinphi^2
    cosphi2 = 1. - sinphi2
    y = 1. - m*sinphi2
    drf,ierr1 = SLATEC.DRF(cosphi2, y, 1.)
    drj,ierr2 = SLATEC.DRJ(cosphi2, y, 1., 1. - n*sinphi2)
    if ierr1 == 0 && ierr2 == 0
        return sinphi*(drf + n*sinphi2*drj/3)
    elseif ierr1 == 2 && ierr2 == 2
        # 2 - (1+m)*sinphi2 < tol
        return Inf
    elseif ierr1 == 0 && ierr2 == 2
        # 1 - n*sinphi2 < tol
        return Inf
    end
    NaN
end
Pi(n::Real, phi::Real, m::Real) = Pi(Float64(n), Float64(phi), Float64(m))
Π = Pi

function ellipj(u::Float64, m::Float64, tol::Float64)
    phi = Jacobi.am(u, m, tol)
    s = sin(phi)
    c = cos(phi)
    d = sqrt(1. - m*s^2)
    s, c, d
end
ellipj(u::Float64, m::Float64) = ellipj(u, m, eps(Float64))
ellipj(u::Real, m::Real) = ellipj(Float64(u), Float64(m))

end # module


text:
# p3.jl - band-limited interpolation

h = 1; xmax = 10; clf();
x = -xmax:h:xmax;                     # computational grid
xx = -xmax-h/20:h/10:xmax+h/20;       # plotting grid
v = zeros(length(x),3);
v[:,1] = @. float(x==0);
v[:,2] = @. float(abs(x) ≤ 3);
v[:,3] = @. max(0,1-abs(x)/3);
for plt = 1:3
    subplot(4,1,plt);
    plot(x,v[:,plt],".",markersize=6);
    p = 0;
    for i = 1:length(x)
        p = @. p + v[i,plt]*sin(pi*(xx-x[i])/h)/(pi*(xx-x[i])/h);
    end
    plot(xx,p,"-")
    axis([-xmax,xmax,-0.5,1.5]);
    xticks(1:0); yticks(0:1);
end


text:
# helps with static types that can't be easily inferred as same parametrically
same_type(::X, ::Y) where {X,Y} =  same_type(X, Y)
same_type(::Type{X}, ::Type{Y}) where {X<:OneToSRange,Y<:OneToSRange} = true
same_type(::Type{X}, ::Type{Y}) where {X<:UnitSRange,Y<:UnitSRange} = true
same_type(::Type{X}, ::Type{Y}) where {X<:StepSRange,Y<:StepSRange} = true
same_type(::Type{X}, ::Type{Y}) where {X<:LinSRange,Y<:LinSRange} = true
same_type(::Type{X}, ::Type{Y}) where {X<:StepSRangeLen,Y<:StepSRangeLen} = true
same_type(::Type{X}, ::Type{X}) where {X} = true
same_type(::Type{X}, ::Type{Y}) where {X,Y} = false

promote_values_rule(::X, ::Y) where {X,Y} = promote_values_rule(X, Y)
function promote_values_rule(::Type{X}, ::Type{Y}) where {X<:AbstractIndex,Y<:AbstractIndex}
    xv = values_type(X)
    yv = values_type(Y)
    return same_type(xv, yv) ? xv : _promote_rule(xv, yv)
end

promote_keys_rule(::X, ::Y) where {X,Y} = promote_keys_rule(X, Y)
function promote_keys_rule(::Type{X}, ::Type{Y}) where {X<:AbstractIndex,Y<:AbstractIndex}
    xv = keys_type(X)
    yv = keys_type(Y)
    return same_type(xv, yv) ? xv : _promote_rule(xv, yv)
end

function _promote_rule(::Type{X}, ::Type{Y}) where {X,Y}
    out = promote_rule(X, Y)
    return out <: Union{} ? promote_rule(Y, X) : out
end

function Base.promote_rule(::Type{X}, ::Type{Y}) where {X<:AbstractIndex,Y<:AbstractIndex}
    return similar_type(X, promote_keys_rule(X, Y), promote_values_rule(X, Y))
end

function Base.promote_rule(::Type{X}, ::Type{Y}) where {X<:AbstractIndex,Y<:AbstractUnitRange}
    return promote_rule(X, similar_type(X, index_keys_type(Y),Y))
end

function Base.promote_rule(::Type{X}, ::Type{Y}) where {X<:UnitRange,Y<:AbstractIndex}
    return promote_rule(Y, X)
end

function Base.promote_rule(::Type{X}, ::Type{Y}) where {X<:AbstractIndex,Y<:AbstractVector}
    return promote_rule(values_type(X), Y)
end
function Base.promote_rule(::Type{X}, ::Type{Y}) where {X<:AbstractVector,Y<:AbstractIndex}
    return promote_rule(Y, X)
end

#function promote_axes(::Type{X}, ::Type{Y}) where {X<:AbstractIndex,Y<:AbstractIndex}
#end

# TODO
#function Base.promote_rule(::Type{}, ::Type{}) where {X<:IndicesArray,Y<:AbstractArray}
#end
#=
function Base.promote_rule(::Type{X}, ::Type{Y}) where {X<:IndicesArray,Y<:IndicesArray}
    return similar_type(X,
                 promote_rule(eltype(X), eltype(Y)),
                 promote_rule(parent_type(X), parent_type(X)),
                 promote_axes(axes_type(X), axes_type(Y))
                )
end

function Base.promote_rule(::Type{X}, ::Type{Y}) where {X<:IndicesArray,Y<:AbstractArray}
end
=#


text:
function mstprimjarnik(G::DictGraph; source = 0)::Tuple{SortedDict{Int, Array{Int,1}},PriorityQueue{Int, Int}}
    𝐷 = SortedDict{Int, Array{Int,1}}(i => [typemax(Int), -1] for i in keys(G.relation))
    𝒬 = PriorityQueue{Int, Int}(i => typemax(Int) for i in vertices(G))

    𝐷[source][1] = 0
    𝒬[source] = 0

    while !isempty(𝒬)
        𝑢 = dequeue_pair!(𝒬)
        for ℒ in G[𝑢.first]                    #for each source adj to 𝑢
            if ℒ[1] ∈ keys(𝒬) && ℒ[2] < 𝒬[ℒ[1]]
                𝐷[ℒ[1]][2] = 𝑢.first;
                𝐷[ℒ[1]][1] = ℒ[2]
                𝒬[ℒ[1]] = ℒ[2]
            end
        end
    end

    #create a new digraph T where T is a spanning tree.
    𝑇 = DictGraph(directed = true)

    for p in 𝐷
        p.first != source && addedge!(𝑇, p.second[2], p.first; weight = p.second[1])
    end

    (𝐷,𝑇)
end


text:
import DataFrames

"""
"""
function DataFrameFeatureContrasts(
        df::DataFrames.AbstractDataFrame,
        columns::AbstractVector{Symbol},
        )
    columns_deepcopy = deepcopy(columns)
    if length(columns_deepcopy) != length(unique(columns_deepcopy))
        error("length(columns_deepcopy) != length(unique(columns_deepcopy))")
    end
    num_df_columns_deepcopy = length(columns_deepcopy)
    formula_wo_intercept = generate_formula(
        columns_deepcopy[1],
        columns_deepcopy;
        intercept = false,
        )
    formula_with_intercept = generate_formula(
        columns_deepcopy[1],
        columns_deepcopy;
        intercept = true,
        )
    schema_wo_intercept = StatsModels.schema(formula_wo_intercept, df)
    schema_with_intercept = StatsModels.schema(formula_with_intercept, df)
    formula_wo_intercept = StatsModels.apply_schema(formula_wo_intercept,
                                                    schema_wo_intercept)
    formula_with_intercept = StatsModels.apply_schema(formula_with_intercept,
                                                      schema_with_intercept)
    response_wo_intercept, predictors_wo_intercept =
        StatsModels.modelcols(formula_wo_intercept, df)
    response_with_intercept, predictors_with_intercept =
        StatsModels.modelcols(formula_with_intercept, df)
    num_array_columns_wo_intercept = size(predictors_wo_intercept, 2)
    num_array_columns_with_intercept = size(predictors_with_intercept, 2)
    result = DataFrameFeatureContrasts(
        columns_deepcopy,
        num_df_columns_deepcopy,
        schema_wo_intercept,
        formula_wo_intercept,
        num_array_columns_wo_intercept,
        schema_with_intercept,
        formula_with_intercept,
        num_array_columns_with_intercept,
        )
    return result
end

"""
"""
function generate_feature_contrasts(
        df::DataFrames.AbstractDataFrame,
        columns::AbstractVector{Symbol},
        )
    result = DataFrameFeatureContrasts(df, columns)
    return result
end


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

include("testkvp.jl")
include("testknr.jl")


text:
#!/usr/bin/env julia6
using ReadGlobal

function model2_part!(r1,r2,r3,t11,t12,t13,t22,t23,t33,out1,out2,out3)
    Threads.@threads for i in linearindices(r1)
        @inbounds out1[i] = t11[i]*r1[i] + t12[i]*r2[i] + t13[i]*r3[i]
        @inbounds out2[i] = t12[i]*r1[i] + t22[i]*r2[i] + t23[i]*r3[i]
        @inbounds out3[i] = t13[i]*r1[i] + t23[i]*r2[i] + t33[i]*r3[i]
    end
end

function main()
    N,Fil = getnfilter()
    r1 = "$(Fil)drhodx_N$N"
    r2 = "$(Fil)drhody_N$N"
    r3 = "$(Fil)drhodz_N$N"
    run(`grad.jl $(Fil)rho_N$N $r1 $r2 $r3`)

    t11 = "$(Fil)T11_N$N"
    t12 = "$(Fil)T12_N$N"
    t13 = "$(Fil)T13_N$N"
    t22 = "$(Fil)T22_N$N"
    t23 = "$(Fil)T23_N$N"
    t33 = "$(Fil)T33_N$N"
    doinchunks(model2_part!, input=(r1,r2,r3,t11,t12,t13,t22,t23,t33), output=("2modelflux1", "2modelflux2", "2modelflux3"))
end

main()


text:
function JuMP.optimize!(m::SJ.StructuredModel; options...)

    # free any existing model pointer
    freeModel(dspenv)

    # set options
    setoptions!(options)

    # classify the second stage constraints into linear ones and quadratic ones
    classifyConstrs(m)

    # load problem
    load_problem!(m)

    # solve DSP problem
    solve!()
    post_solve!()

    return termination_status(m)
end

function JuMP.termination_status(m::StructuredModel)
    status = MOI.OPTIMIZE_NOT_CALLED
    if dspenv.status == 3000
        status = MOI.OPTIMAL
    elseif dspenv.status == 3005
        status = MOI.ALMOST_OPTIMAL
    elseif dspenv.status == 3001
        status = MOI.INFEASIBLE
    elseif dspenv.status == 3002
        status = MOI.DUAL_INFEASIBLE
    elseif dspenv.status in [3004,3007]
        status = MOI.TIME_LIMIT
    elseif dspenv.status == 3010
        status = MOI.ITERATION_LIMIT
    elseif dspenv.status in [3014,3015]
        status = MOI.OBJECTIVE_LIMIT
    elseif dspenv.status == 3006
        status = MOI.NODE_LIMIT
    elseif dspenv.status in [3008,3009,3016]
        status = MOI.OTHER_LIMIT
    else
        status = MOI.OTHER_ERROR
    end
    return status
end

JuMP.objective_value(m::StructuredModel) = dspenv.primVal
JuMP.dual_objective_value(m::StructuredModel) = dspenv.dualVal

function JuMP.value(v::SJ.StructuredVariableRef; result = 1)::Float64
    if !isnothing(v.model.parent)
        @warn("Solution values are available for the parent model only. Please use value() to access all children solutions in vector form.")
        return NaN
    else
        return dspenv.colVal[0][v.idx]
    end
end
JuMP.value() = dspenv.colVal
JuMP.dual() = dspenv.rowVal

JuMP.solve_time(m::StructuredModel) = dspenv.solve_time


text:
using SimpleTraits

# define asset types
include("../src/asset_types.jl")

@traitdef IsLiquid{T}
@traitimpl IsLiquid{Cash}
@traitimpl IsLiquid{Investment}

@traitfn marketprice(x::::IsLiquid) =
    error("Please implement pricing function for ", typeof(x))

@traitfn marketprice(x::::(!IsLiquid)) =
    error("Price for illiquid asset $x is not available.")

marketprice(Stock("AAPL", "Apple"))
marketprice(Residence("Los Angeles"))

marketprice(x::Stock) = 123
marketprice(Stock("AAPL", "Apple"))


text:
"""
    R2BPdynamics(rv, μ, t)

Compute time derivative of state vector in the restricted two-body system. `rv` is the state
vector `[r; v]` {km; km/s}, `μ` is the gravitational parameter {km³/s²}, and `t` is time
{s}.
"""
function R2BPdynamics(rv, μ, t)  #make sure rv and μ are in km and km³/s²
    r,v = rv[1:3], rv[4:6]
    rvdot = zeros(6)
    rvdot[1:3] = v
    rvdot[4:6] = -μ/norm(r)^3 * r   #acceleration [km/s]
    return rvdot
end

"""
    R2BPdynamics(rv, prim::Body, t)

Compute time derivative of state vector in the restricted two-body system. `rv` is the state
vector `[r; v]` {km; km/s}, `prim` is the central body, and `t` is time {s}.
"""
function R2BPdynamics(rv, prim::Body, t)  #make sure rv and μ are in km and km³/s²
    return R2BPdynamics(rv, G*prim.m, t)
end

"""
    R2BPdynamics!(rvdot, rv, μ, t)

In-place version of `R2BPdynamics(rvdot, rv, μ, t)`.
"""
function R2BPdynamics!(rvdot, rv, μ, t)  #make sure rv and μ are in km and km³/s²
    rvdot[:] = R2BPdynamics(rv,μ,t)
    return nothing
end

"""
    R2BPdynamics!(rvdot, rv, prim::Body, t)

In-place version of `R2BPdynamics!(rvdot, rv, prim::Body, t)`.
"""
function R2BPdynamics!(rvdot, rv, prim::Body, t)  #make sure rv and μ are in km and km³/s²
    rvdot[:] = R2BPdynamics(rv, prim, t)
    return nothing
end

"""
    CR3BPdynamics(rv, μ, t)

Compute time derivative of state vector `rv = [r; v]` {NON, NON} in the rotating frame of
the normalized CR3BP where `μ` is the CR3BP mass parameter μ₂/(μ₁+μ₂) {NON} and `t` is time
{NON}.
"""
function CR3BPdynamics(rv,μ,t) #Three body dynamics in Earth/Moon System
    x,y,z,vx,vy,vz = rv
    r₁³= ((x + μ)^2     + y^2 + z^2)^1.5; # distance to m1, LARGER MASS
    r₂³= ((x - 1 + μ)^2 + y^2 + z^2)^1.5; # distance to m2, smaller mass

    rvdot = zeros(6)
    rvdot[1:3] = [vx;vy;vz]
    rvdot[4] = -((1 - μ)*(x + μ)/r₁³) - (μ*(x - 1 + μ)/r₂³) + 2*vy + x;
    rvdot[5] = -((1 - μ)*y      /r₁³) - (μ*y          /r₂³) - 2*vx + y;
    rvdot[6] = -((1 - μ)*z      /r₁³) - (μ*z          /r₂³);
    return rvdot
end

"""
    CR3BPdynamics(rv, sys::System, t)

Compute time derivative of state vector `rv = [r; v]` {NON, NON} in the rotating frame of
the normalized CR3BP where `sys` is the CR3BP system and `t` is time {NON}.
"""
function CR3BPdynamics(rv,sys::System,t) #Three body dynamics in Earth/Moon System
    return CR3BPdynamics(rv,sys.μ,t)
end

"""
    CR3BPdynamics(rv, p::Array, t)

Compute time derivative of state vector `rv = [r; v]` {km, km/s} in the rotating frame of
the non-normalized CR3BP where `p = [μ₁;μ₂;d]` {km³/s²; km³/s²; km} contains the
gravitational parameters of the first and second primary bodies as well as the distance
between them. `t` is time {s}.
"""
function CR3BPdynamics(rv,p::Array,t) #Three body dynamics in Earth/Moon System
    x,y,z,vx,vy,vz = rv
    μ₁,μ₂,d = p # parameters
    R₁ = d*μ₂/(μ₁+μ₂)
    R₂ = d*μ₁/(μ₁+μ₂)
    ωₛ = sqrt((μ₁ + μ₂)/d^3) #rotation rate of system
    r₁³= ((x+R₁)^2 + y^2 + z^2)^1.5; # distance to m1, LARGER MASS
    r₂³= ((x-R₂)^2 + y^2 + z^2)^1.5; # distance to m2, smaller mass

    rvdot = zeros(6)
    rvdot[1:3] = [vx;vy;vz]
    rvdot[4]   = -(μ₁*(x+R₁)/r₁³) - (μ₂*(x-R₂)/r₂³) + 2*ωₛ*vy + ωₛ^2*x;
    rvdot[5]   = -(μ₁*y     /r₁³) - (μ₂*y     /r₂³) - 2*ωₛ*vx + ωₛ^2*y;
    rvdot[6]   = -(μ₁*z     /r₁³) - (μ₂*z     /r₂³);
    return rvdot
end

"""
    CR3BPdynamics!(rvdot, rv, μ, t)

In-place version of `CR3BPdynamics(rv, μ, t)`.
"""
function CR3BPdynamics!(rvdot,rv,μ,t) #Three body dynamics in Earth/Moon System
    rvdot[:] = CR3BPdynamics(rv,μ,t)
    return nothing
end

"""
    CR3BPdynamics!(rvdot, rv, sys::System, t)

In-place version of `CR3BPdynamics(rv, sys::System, t)`.
"""
function CR3BPdynamics!(rvdot,rv,sys::System,t) #Three body dynamics in Earth/Moon System
    rvdot[:] = CR3BPdynamics(rv,sys,t)
    return nothing
end

"""
    CR3BPdynamics!(rvdot, rv, p::Array, t)

In-place version of `CR3BPdynamics(rv, p::Array, t)`.
"""
function CR3BPdynamics!(rvdot,rv,p::Array,t) #Three body dynamics in Earth/Moon System
    rvdot[:] = CR3BPdynamics(rv,p,t)
    return nothing
end

"""
    CR3BPstm(w, μ, t)

Compute time derivative of state vector `w = [r; v; vec(Φ)]` {NON; NON; NON} in the rotating
frame of the normalized CR3BP. `vec(Φ)` is the vectorized state transition matrix while `μ`
is the CR3BP mass parameter μ₂/(μ₁+μ₂) {NON} and `t` is time {NON}.
"""
function CR3BPstm(w,μ,t) #Three body dynamics in Earth/Moon System
    rv = w[1:6]
    Φ = reshape(w[7:42],6,6)
    x,y,z,vx,vy,vz = rv

    r₁³= ((x + μ)^2     + y^2 + z^2)^1.5; # distance to m1, LARGER MASS
    r₂³= ((x - 1 + μ)^2 + y^2 + z^2)^1.5; # distance to m2, smaller mass
    r₁⁵= ((x + μ)^2     + y^2 + z^2)^2.5; # distance to m1, LARGER MASS
    r₂⁵= ((x - 1 + μ)^2 + y^2 + z^2)^2.5; # distance to m2, smaller mass

    omgxx = 1 - (1-μ)*(1/r₁³ - 3*(x + μ)^2/r₁⁵) - μ*(1/r₂³ - 3*(x - 1 + μ)^2/r₂⁵);
    omgxy = 3*(1-μ)*(x + μ)*y/r₁⁵ + 3*μ*(x - 1 + μ)*y/r₂⁵;
    omgxz = 3*(1-μ)*(x + μ)*z/r₁⁵ + 3*μ*(x - 1 + μ)*z/r₂⁵;
    omgyy = 1 - (1-μ)*(1/r₁³ - 3*y^2/r₁⁵) - μ*(1/r₂³ - 3*y^2/r₂⁵);
    omgyz = 3*(1-μ)*y*z/r₁⁵ + 3*μ*y*z/r₂⁵;
    omgzz = - (1-μ)*(1/r₁³ - 3*z^2/r₁⁵) - μ*(1/r₂³ - 3*z^2/r₂⁵);

    F = [   0     0     0     1     0	 0 ;
            0     0     0     0 	1 	 0 ;
            0	  0     0     0     0    1 ;
        omgxx omgxy omgxz     0     2 	 0 ;
        omgxy omgyy omgyz    -2     0 	 0 ;
        omgxz omgyz omgzz     0	    0	 0 ];

    Φdot = F*Φ;
    wdot = zeros(42)
    wdot[1:6] = CR3BPdynamics(rv,μ,t)
    wdot[7:42] = reshape(Φdot, 36, 1);
    return wdot
end

"""
    CR3BPstm(w, sys, t)

Compute time derivative of state vector `w = [r; v; vec(Φ)]` {NON; NON; NON} in the rotating
frame of the normalized CR3BP. `vec(Φ)` is the vectorized state transition matrix while
`sys` is the CR3BP system and `t` is time {NON}.
"""
function CR3BPstm(w,sys::System,t)
    return CR3BPstm(w,sys.μ,t)
end

"""
    CR3BPstm!(wdot, w, μ, t)

In-place version of `CR3BPstm(w, μ, t)`.
"""
function CR3BPstm!(wdot,w,μ,t)
    wdot[:] = CR3BPstm(w,μ,t)
    return nothing
end

"""
    CR3BPstm!(wdot, w, sys::System, t)

In-place version of `CR3BPstm(w, sys::System, t)`.
"""
function CR3BPstm!(wdot,w,sys::System,t) #Three body dynamics in Earth/Moon System
    wdot[:] = CR3BPstm(w,sys,t)
    return nothing
end

"""
    CR3BPinert(rv,μ,t)

Compute time derivative of state vector `rv = [r; v]` {NON, NON} in the inertial frame of
the normalized CR3BP where `μ` is the CR3BP mass parameter μ₂/(μ₁+μ₂) {NON} and `t` is time
{NON}.
"""
function CR3BPinert(rv,μ,t)
    x,y,z,vx,vy,vz = rv
    r₁ = [x +       μ*cos(t); y +       μ*sin(t); 0]
    r₂ = [x - (1 - μ)*cos(t); y - (1 - μ)*cos(t); 0]
    r₁³ = norm(r₁)^3
    r₂³ = norm(r₂)^3
    rvdot = zeros(6)
    rvdot[1:3] = [vx;vy;vz]
    rvdot[4:6] = -(1 - μ)*r₁/r₁³ - μ*r₂/r₂³
    return rvdot
end

"""
    CR3BPinert(rv, sys::System, t)

Compute time derivative of state vector `rv = [r; v]` {NON, NON} in the inertial frame of
the normalized CR3BP where `sys` is the CR3BP system and `t` is time {NON}.
"""
function CR3BPinert(rv,sys::System,t)
    return CR3BPinert(rv,sys.μ,t)
end

"""
    CR3BPinert(rvdot, rv, p::Array, t)

Compute time derivative of state vector `rv = [r; v]` {NON, NON} in the inertial frame of
the non-normalized CR3BP where `p = [μ₁;μ₂;d]` {km³/s²; km³/s²; km} contains the
gravitational parameters of the first and second primary bodies as well as the distance
between them. `t` is time {s}.
"""
function CR3BPinert(rv,p::Array,t)
    x,y,z,vx,vy,vz = rv
    μ₁,μ₂,d = p # parameters
    R₁,R₂ = computeR1R2(p)
    ωₛ = sqrt((μ₁ + μ₂)/d^3);
    r₁ = [x + R₁*cos(ωₛ*t); y + R₁*sin(ωₛ*t); 0] # distance to m1, LARGER MASS
    r₂ = [x - R₂*cos(ωₛ*t); y - R₂*sin(ωₛ*t); 0] # distance to m2, smaller mass
    r₁³= norm(r₁)^3;
    r₂³= norm(r₂)^3;
    rvdot = zeros(6)
    rvdot[1:3] = rv[4:6];
    rvdot[4:6] = -μ₁*r₁/r₁³ - μ₂*r₂/r₂³;
    return rvdot
end

"""
    CR3BPinert!(rvdot, rv, μ, t)

In-place version of `CR3BPinert(rv, μ, t)`.
"""
function CR3BPinert!(rvdot,rv,μ,t)
    rvdot[:] = CR3BPinert(rv,μ,t)
    return nothing
end

"""
    CR3BPinert!(rvdot, rv, sys::System, t)

In-place version of `CR3BPinert(rv, sys::System, t)`.
"""
function CR3BPinert!(rvdot,rv,sys::System,t)
    rvdot[:] = CR3BPinert(rv,sys,t)
    return nothing
end

"""
    CR3BPinert!(rvdot, rv, p::Array, t)

In-place version of `CR3BPinert(rv, p::Array, t)`.
"""
function CR3BPinert!(rvdot,rv,p::Array,t)
    rvdot[:] = CR3BPinert(rv,p,t)
    return nothing
end

"""
    CWdynamics(rv, n, t)

Clohessy-Wiltshire equations. Compute time derivative of state vector `rv = [δr; δv]`
{km; km/s} where `n` {rad/s} is the mean motion of the chief and `t` is time {s}.
"""
function CWdynamics(rv,n,t)
    x,y,z,vx,vy,vz = rv
    rvdot = zeros(6)
    rvdot[1:3] = rv[4:6]
    rvdot[4] = 2*n*vy + 3*n^2*x
    rvdot[5] = -2*n*vx
    rvdot[6] = -n^2*z
    return rvdot
end

"""
    CWdynamics!(rvdot,rv,n,t)

Clohessy-Wiltshire equations

Inputs: n (scalar) mean motion
"""
function CWdynamics!(rvdot,rv,n,t)
    rvdot[:] = CWdynamics(rv,n,t)
    return nothing
end

"""
See G. Gómez, C. Simó, J. Llibre, and R. Martínez, Dynamics and mission design near
libration points. Vol. II, vol. 3. 2001.

    BCPdynamics(rv, μ, m₃, n₃, t)

Compute time derivative of state vector `rv = [r; v]` {km; km/s} in the normalized
Bicircular Four-Body Problem (BCP). `μ` {NON} is the BCP mass parameter and `m₃` {NON} and
`n₃` {NON} are the normalized mass and mean motion of the tertiary body. `t` is time {NON}.
"""
function BCPdynamics(rv, μ, m₃, n₃, t)
    x,y,z,vx,vy,vz = rv

    # Xₑ = μ*cos(t);      Yₑ = μ*sin(t);      Zₑ = 0
    # Xₘ = (μ-1)*cos(t);  Yₘ = (μ-1)*sin(t);  Zₘ = 0
    # Xₛ = a₃*cos(n₃*t);  Yₛ = a₃*sin(n₃*t);  Zₛ = 0

    a₃ = (1+m₃)^(1/3)/n₃^(2/3)
    θ = (1-n₃)*t
    x₃ =  a₃*cos(θ)
    y₃ = -a₃*sin(θ)

    r₁³ = (  (x+μ)^2 +      y^2 + z^2)^1.5 # distance to m1, LARGER MASS
    r₂³ = ((x-1+μ)^2 +      y^2 + z^2)^1.5 # distance to m2, smaller mass
    r₃³ = ( (x-x₃)^2 + (y-y₃)^2 + z^2)^1.5

    rvdot = zeros(6)
    rvdot[1:3] = [vx;vy;vz]
    rvdot[4] = -(1-μ)*(x+μ)/r₁³ - μ*(x-1+μ)/r₂³ - m₃*(x-x₃)/r₃³ - m₃*cos(θ)/a₃^2 + 2*vy + x
    rvdot[5] = -(1-μ)  *  y/r₁³ - μ   *   y/r₂³ - m₃*(y-y₃)/r₃³ + m₃*sin(θ)/a₃^2 - 2*vx + y
    rvdot[6] = -(1-μ)  *  z/r₁³ - μ   *   z/r₂³ - m₃  *   z/r₃³
    return rvdot
end

"""
    BCPdynamics(rv, sys::BicircularSystem, t)

Compute time derivative of state vector `rv = [r; v]` {km; km/s} in the normalized
Bicircular Four-Body Problem (BCP). `sys` is the BCP system and `t` is time {NON}.
"""
function BCPdynamics(rv, sys::BicircularSystem, t)
    return BCPdynamics(rv, sys.μ, sys.m₃, sys.n₃, t)
end

"""
    BCPdynamics!(rvdot, rv, μ, m₃, n₃, t)

In-place version of `BCPdynamics(rv, μ, m₃, n₃, t)`.
"""
function BCPdynamics!(rvdot, rv, μ, m₃, n₃, t)
    rvdot[:] = BCPdynamics(rv, μ, m₃, n₃, t)
    return nothing
end

"""
    BCPdynamics!(rvdot, rv, sys::BicircularSystem, t)

In-place version of `BCPdynamics(rv, sys::BicircularSystem, t)`.
"""
function BCPdynamics!(rvdot, rv, sys::BicircularSystem, t)
    rvdot[:] = BCPdynamics(rv, sys, t)
    return nothing
end

"""
    BCPstm(wdot, w, μ, m₃, n₃, t)

Compute time derivative of state vector `w = [r; v; vec(Φ)]` {NON; NON; NON} in the
normalized Bicircular Four-Body Problem (BCP). `vec(Φ)` is the vectorized state transition
matrix. `μ` {NON} is the BCP mass parameter and `m₃` {NON} and `n₃` {NON} are the normalized
mass and mean motion of the tertiary body. `t` is time {NON}.
"""
function BCPstm(w, μ, m₃, n₃, t) #Three body dynamics in Earth/Moon System
    rv = w[1:6]
    Φ = reshape(w[7:42],6,6)
    x,y,z,vx,vy,vz = rv

    a₃ = (1+m₃)^(1/3)/n₃^(2/3)
    θ = (1-n₃)*t
    x₃ =  a₃*cos(θ)
    y₃ = -a₃*sin(θ)

    r₁ = sqrt(  (x+μ)^2 +      y^2 + z^2) # distance to m1, Larger Mass
    r₂ = sqrt((x-1+μ)^2 +      y^2 + z^2) # distance to m2, smaller mass
    r₃ = sqrt( (x-x₃)^2 + (y-y₃)^2 + z^2) # distance to m3, LARGEST MASS
    r₁³ = r₁^3; r₂³ = r₂^3; r₃³ = r₃^3
    r₁⁵ = r₁^5; r₂⁵ = r₂^5; r₃⁵ = r₃^5

    g11 = 1 - (1-μ)*(1/r₁³ - 3*(x + μ)^2/r₁⁵) - μ*(1/r₂³ - 3*(x - 1 + μ)^2/r₂⁵) - m₃*(1/r₃³ - 3*(x-x₃)^2/r₃⁵)
    g12 = 3*(1-μ)*(x + μ)*y/r₁⁵ + 3*μ*(x - 1 + μ)*y/r₂⁵ + 3*m₃*(x-x₃)*(y-y₃)/r₃⁵
    # g13 = 0 # planar assumption
    g13 = 3*(1-μ)*(x + μ)*z/r₁⁵ + 3*μ*(x - 1 + μ)*z/r₂⁵ + 3*m₃*(x-x₃)*z/r₃⁵
    g22 = 1 - (1-μ)*(1/r₁³ - 3*y^2/r₁⁵) - μ*(1/r₂³ - 3*y^2/r₂⁵) - m₃*(1/r₃³ - 3*(y-y₃)^2/r₃⁵)
    # g23 = 0 # planar assumption
    g23 = 3*(1-μ)*y*z/r₁⁵ + 3*μ*y*z/r₂⁵ + 3*m₃*(y-y₃)*z/r₃⁵
    # g33 = -(1-μ)/r₁³ - μ/r₂³ - m₃/r₃³ # planar assumption
    g33 = -(1-μ)*(1/r₁³ - 3*z^2/r₁⁵) - μ*(1/r₂³ - 3*z^2/r₂⁵) - m₃*(1/r₃³ - 3*z^2/r₃⁵)

    F = [  0   0   0  1   0   0 ;
           0   0   0  0   1   0 ;
           0   0   0  0   0   1 ;
         g11 g12 g13  0   2   0 ;
         g12 g22 g23 -2   0   0 ;
         g13 g23 g33  0	  0   0 ]

    Φdot = F*Φ
    wdot[1:6] = BCPdynamics(rv,μ,m₃,n₃,t)
    wdot[7:42] = reshape(Φdot, 36, 1);
    return wdot
end

"""
    BCPstm(wdot, w, μ, m₃, n₃, t)

Compute time derivative of state vector `w = [r; v; vec(Φ)]` {NON; NON; NON} in the
normalized Bicircular Four-Body Problem (BCP). `vec(Φ)` is the vectorized state transition
matrix, `sys` is the BCP system and `t` is time {NON}.
"""
function BCPstm(w,sys::BicircularSystem,t) #Three body dynamics in Earth/Moon System
    return BCPstm(w, sys.mu, sys.m3, sys.n3 ,t)
end

"""
    BCPstm!(wdot, w, μ, m₃, n₃, t)

In-place version of `BCPstm(w, μ, m₃, n₃, t)`.
"""
function BCPstm!(wdot,w,μ,m₃,n₃,t) #Three body dynamics in Earth/Moon System
    wdot[:] = BCPstm(w,μ,m₃,n₃,t)
    return nothing
end

"""
    BCPstm!(wdot, w, sys::BicircularSystem, t)

In-place version of `BCPstm(w, sys::BicircularSystem, t)`.
"""
function BCPstm!(wdot,w,sys::BicircularSystem,t) #Three body dynamics in Earth/Moon System
    wdot[:] = BCPstm(w,sys,t)
    return nothing
end


text:
uri_pages = ["index", "home", "search", "main", "post", "homepage", "category", "register", "login", "faq", "about", "terms","privacy", "author"]
uri_paths = ["app", "main", "wp-content", "search", "category", "tag", "categories", "tags", "blog", "posts", "list", "explore"]
uri_extensions = [".html", ".html", ".html", ".htm", ".htm", ".php", ".php", ".jsp", ".asp"]

pixels()= rand(1:1024)

user_name()= bothify(executor(data["en"]["faker"]["internet"]["user_name_formats"]))

domain_suffix()= executor(data["en"]["faker"]["internet"]["domain_suffix"])
safe_email()= "$(user_name())@example.$(domain_suffix())"

free_email_domain() = executor(data["en"]["faker"]["internet"]["free_email"])
free_email()= "$(user_name())@$(free_email_domain())"

domain_word()=
   ( company_elements = split(company(), [',', ' ','-']);
    sal = "";
    for i in company_elements
               if length(i)>0
                   sal*="$(i[1])"
               end
          end;

    lowercase(sal) )

domain_name()="$(domain_word()).$(domain_suffix())"

company_email()= "$(user_name())@$(domain_name())"

email()= executor(data["en"]["faker"]["internet"]["email_formats"])

url()= executor(data["en"]["faker"]["internet"]["url_formats"])

ipv4()= "$(rand( 172:255)).$(rand( 0:255)).$(rand(0:255)).$(rand(0:255))"

ipv6()= (res = "";
    for i = 1:8
        res*= "$(string(rand( 4100:65535),base=16))"
        if i < 8
            res*= ":"
        end
    end;res)

mac_address()=
    (res = "";
    for i = 1:6
        res*= "$(Base.hex(unsigned(rand( 0x00:0xff)),2,false))"
        i < 6 && (res*= ":")
    end;res)

uri_page()= executor(uri_pages)

uri_path(deep="None")=
    (if deep =="None"
        deep = rand(1:3)
    end; sa="";
        for i =0:deep
           sa*= uri_paths[rand( 1:length(uri_paths))]
            i < deep && (sa*= "/")
        end; sa )

uri_extension()= executor(uri_extensions)
uri()= executor(data["en"]["faker"]["internet"]["uri_formats"])
image_url()= executor(data["en"]["faker"]["internet"]["image_placeholder_services"])


text:
using NaturalES

sols=[]
cnt=0
function rosenbrock(x::AbstractVector{T}) where T
    global cnt=(cnt+1)%20
    cnt==10 && push!(sols,copy(x))
    s=(1.0 - x[1])^2
    for i in 1:(length(x)-1)
        s+=100.0 * (x[i+1] - x[i]^2)^2
    end
    return s
end

exponential_nes(rosenbrock,[-0.41,0.17],0.01)
sols
A=hcat(sols...)'|>collect

using PyPlot

pygui(true)
function rosenbrock2d(x,y)
  return (1.0 - x)^2 + 100.0 * (y - x^2)^2
end

x=-0.5:0.001:1.5
y=-0.5:0.001:1.5
M=rosenbrock2d.(x,y')
plt.style.use("dark_background")
plot(A[:,2],A[:,1],color="red")
contour(x,y,M,log.(1.01:0.5:10))
colorbar()
scatter(1,1,color="blue")

using NaturalES
function rosenbrock2d(x::AbstractVector{T}) where T
    s=(1.0 - x[1])^2
    for i in 1:(length(x)-1)
        s+=100.0 * (x[i+1] - x[i]^2)^2
    end
    return s
end

using BenchmarkTools
@btime exponential_nes($rosenbrock2d,$[0.0,0.0],$1.0)


text:
####################################################
# Script to train a spiking (LIF, STDP), patchy RP/RG network with 1 hidden layer on MNIST.
# Event-based and euler-forward integration are available, the latter being faster
# depending on the network size, firing rates and of course accepted error made by
# the time discretization (Δt keyword in EulerNet function).
# Simulations might take several days (!) for achieving good accuracy (test acc. > 98%).

###################################################
# Include libraries and core

using Pkg; Pkg.activate("./../core/"); Pkg.instantiate()
push!(LOAD_PATH, string(pwd(),"./../core/src/lifintegrator/core/"))
using EventBasedIntegrator, LinearAlgebra, ProgressMeter

###################################################
# Define task and network

# This keyword decides whether Random Projections ("RP")
# or Random Gabor filters ("RG") are used for feature extraction:
mode = "RG"
# OR "RP"

n_hidden = 100 # number of hidden neurons
patchsize = 10 # (linear) patch size of receptive fields in first layer p ∈ [1,28]
n_iterations = 10^6 # number of iterations (images) for training
learningrate = 2e-4 # learning rate for STDP in second layer
n_eval_trn = n_eval_tst = 10^4 # number of patterns used for evaluation (train & test).
# This is done because evaluating the whole train and test sets might take a long time.
# For the paper the whole data sets were evaluated: (n_eval_trn = 6 * 10^4, n_eval_tst = 10^4)
doEuler = true # decides whether Euler forward integration (true) or event based integration is used (false)

###################################################
# Load and preprocess data

using Knet; include(Knet.dir("data", "mnist.jl"))
function import_mnist()
    smallimgs, labels, smallimgstest, labelstest = mnist()
    smallimgs = reshape(smallimgs, 28^2, 60000)
    smallimgstest = reshape(smallimgstest, 28^2, 10000)
	return smallimgs, labels .- 1, smallimgstest, labelstest .- 1
end

smallimgs, labels, smallimgstest, labelstest = import_mnist()
subtractmean!(smallimgs) # substract pixel-wise mean
subtractmean!(smallimgstest)

###################################################
# Setup network and plasticity

# set additional parameters
n_in = size(smallimgs)[1] # number of input neurons
n_out = 10 # number of output neurons
n_total = n_in + n_hidden + n_out # total number of neurons
threshold = 20. # mV
in_amplitude = 500 # mA
teacher_amplitude = 0.5 # kHz
pat_duration = 5e1 # learning rate should be scaled inversely to this!
transient_duration = 1e2 # ms
pat_duration_test = 2e2 # ms

## Set plasticity rule
plasticityrule = TeacherSignalPlasticity(falses(n_hidden * (n_in + n_out)), zeros(n_total), learningrate)
plasticityrule.isweightplastic[end - n_hidden * n_out + 1:end] .= true

## Setup network
layer_sizes = [n_in,n_hidden,n_out]
network = FeedForwardNet(layer_sizes, LIFParams(threshold = threshold .+ 1. .* rand(n_total));
                        recorder = NoRecorder, plasticityrule = plasticityrule)
if mode == "RP"
    set_connectivity!(network, n_in, n_hidden, patchsize)
elseif mode == "RG"
    set_connectivity_gabor!(network, n_in, n_hidden, patchsize)
else
    error("mode must be either 'RP' or 'RG'")
end
## Set forward euler integration if applicable
doEuler && (network = EulerNet(network; Δt = 5e-2))

###################################################
# Train network

print("Train network... \n")
let endtime = 0
elapsed = @elapsed @showprogress for i in 1:n_iterations
  input, label = getMNISTinput(smallimgs, labels, threshold, n_hidden, n_out;
    in_amplitude = in_amplitude)
  plasticityrule.islearning = false # learning off during transient
  endtime += transient_duration
  integratenet!(network, input, endtime) # integrate transient phase

  plasticityrule.islearning = true # enable learning and define target
  plasticityrule.targets[end-n_out+1:end] = [Int(i == label+1) for i in 1:n_out]*teacher_amplitude
  endtime += pat_duration
  integratenet!(network, input, endtime) # integrate learning period
end
end

###################################################
# Get errors

print("Test phase... \n")
error_train, endtime = geterror!(network, smallimgs, labels, threshold,
                    n_hidden, n_out, in_amplitude, n_eval_trn, pat_duration_test)
error_test, endtime = geterror!(network, smallimgstest, labelstest, threshold,
                    n_hidden,n_out, in_amplitude, n_eval_tst, pat_duration_test)
print("\n Train error:")
print(error_train)
print("\n Test error:")
print(error_test)
print("\n")


text:
import OrdinaryDiffEq, DiffEqCallbacks, DelayDiffEq

function recast! end

"""
    integrate(tspan, df::Function, x0::Vector{ComplexF64},
            state::T, dstate::T, fout::Function; kwargs...)

Integrate using OrdinaryDiffEq
"""
function integrate(
    tspan,
    df::Function,
    x0::X,
    state::T,
    dstate::T,
    fout::Function;
    reltol::Float64 = 1.0e-12,
    abstol::Float64 = 1.0e-12,
    alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = OrdinaryDiffEq.DP5(),
    steady_state = false,
    tol = 1e-3,
    save_everystep = false,
    saveat = tspan,
    callback = nothing,
    kwargs...,
) where {T,X}

    function df_(dx::T, x::T, p, t) where {T}
        recast!(x, state)
        recast!(dx, dstate)
        df(t, state, dstate)
        recast!(dstate, dx)
    end
    function fout_(x, t, integrator)
        recast!(x, state)
        fout(t, state)
    end

    out_type = pure_inference(fout, Tuple{eltype(tspan),typeof(state)})

    out = DiffEqCallbacks.SavedValues(eltype(tspan), out_type)

    scb = DiffEqCallbacks.SavingCallback(
        fout_,
        out,
        saveat = saveat,
        save_everystep = save_everystep,
        save_start = false,
    )

    prob = OrdinaryDiffEq.ODEProblem{true}(df_, x0, (tspan[1], tspan[end]))

    if steady_state
        affect! = function (integrator)
            !save_everystep && scb.affect!(integrator, true)
            OrdinaryDiffEq.terminate!(integrator)
        end
        _cb = OrdinaryDiffEq.DiscreteCallback(
            SteadyStateCondtion(copy(state), tol, state),
            affect!;
            save_positions = (false, false),
        )
        cb = OrdinaryDiffEq.CallbackSet(_cb, scb)
    else
        cb = scb
    end

    full_cb = OrdinaryDiffEq.CallbackSet(callback, cb)

    sol = OrdinaryDiffEq.solve(
        prob,
        alg;
        reltol = reltol,
        abstol = abstol,
        save_everystep = false,
        save_start = false,
        save_end = false,
        callback = full_cb,
        kwargs...,
    )
    out.t, out.saveval
end

function integrate(
    tspan,
    df::Function,
    x0::X,
    state::T,
    dstate::T,
    ::Nothing;
    kwargs...,
) where {T,X}
    function fout(t, state::T)
        copy(state)
    end
    integrate(tspan, df, x0, state, dstate, fout; kwargs...)
end

function integrate_delayed(
    tspan,
    df::Function,
    h::Function,
    x0::X,
    state::T,
    dstate::T,
    fout::Function;
    p = nothing,
    reltol::Float64 = 1.0e-6,
    abstol::Float64 = 1.0e-6,
    alg::Any = DelayDiffEq.MethodOfSteps(DelayDiffEq.Vern6()),
    steady_state = false,
    tol = 1e-3,
    save_everystep = false,
    saveat = tspan,
    callback = nothing,
    kwargs...,
) where {T,X}

    function df_(dx::T, x::T, h, p, t) where {T}
        recast!(x, state)
        recast!(dx, dstate)
        df(t, state, dstate, h, p)
        recast!(dstate, dx)
    end
    function fout_(x, t, integrator)
        recast!(x, state)
        fout(t, state)
    end

    out_type = pure_inference(fout, Tuple{eltype(tspan),typeof(state)})

    out = DiffEqCallbacks.SavedValues(eltype(tspan), out_type)

    scb = DiffEqCallbacks.SavingCallback(
        fout_,
        out,
        saveat = saveat,
        save_everystep = save_everystep,
        save_start = false,
    )

    prob = DelayDiffEq.DDEProblem{true}(df_, x0, h, (tspan[1], tspan[end]), p)

    if steady_state
        affect! = function (integrator)
            !save_everystep && scb.affect!(integrator, true)
            DelayDiffEq.terminate!(integrator)
        end
        _cb = DelayDiffEq.DiscreteCallback(
            SteadyStateCondtion(copy(state), tol, state),
            affect!;
            save_positions = (false, false),
        )
        cb = DelayDiffEq.CallbackSet(_cb, scb)
    else
        cb = scb
    end

    full_cb = DelayDiffEq.CallbackSet(callback, cb)

    sol = DelayDiffEq.solve(
        prob,
        alg;
        reltol = reltol,
        abstol = abstol,
        save_everystep = false,
        save_start = false,
        save_end = false,
        callback = full_cb,
        kwargs...,
    )
    out.t, out.saveval
end

function integrate_delayed(
    tspan,
    df::Function,
    h::Function,
    x0::X,
    state::T,
    dstate::T,
    ::Nothing;
    kwargs...,
) where {T,X}
    function fout(t, state::T)
        copy(state)
    end
    integrate_delayed(tspan, df, h::Function, x0, state, dstate, fout; kwargs...)
end

struct SteadyStateCondtion{T,T2,T3}
    rho0::T
    tol::T2
    state::T3
end
function (c::SteadyStateCondtion)(rho, t, integrator)
    recast!(rho, c.state)
    dt = integrator.dt
    drho = tracedistance(c.rho0, c.state)
    c.rho0.data[:] = c.state.data
    drho / dt < c.tol
end

const QO_CHECKS = Ref(true)
#=
"""
    @skiptimechecks

Macro to skip checks during time-dependent problems.
Useful for [`master_dynamic`](@ref) and similar functions.
"""
=#
macro skiptimechecks(ex)
    return quote
        QO_CHECKS.x = false
        local val = $(esc(ex))
        QO_CHECKS.x = true
        val
    end
end

Base.@pure pure_inference(fout, T) = Core.Compiler.return_type(fout, T)

function recast!(x::Union{Vector,SubArray}, rho::Operator{B,B,T}) where {B<:Basis,T}
    rho.data = reshape(x, size(rho.data))
end
recast!(state::Operator{B,B}, x::SubArray) where {B<:Basis} = (x[:] = state.data)
recast!(state::Operator{B,B}, x::Vector) where {B<:Basis} = nothing

# Recasting needed for the ODE solver is just providing the underlying data
function recast!(x::T, rho::Operator{B,B,T}) where {B<:Basis,T}
    rho.data = x
end
recast!(rho::Operator{B,B,T}, x::T) where {B<:Basis,T} = nothing


text:
"""
$(SIGNATURES)
returns the genotype data. Will skip every other column because genotype probability is duplicated.

"""
function get_geno_data(file, datatype)

    geno_prob = convert(Array{datatype,2},readdlm(file, ','; skipstart=1))
    return geno_prob[:,1:2:end]
end

# function try_string2num(num)
#     return tryparse(Float64,num) != nothing
# end

"""
$(SIGNATURES)
returns the phenotype data. If transposed=true, then the data will be transposed.

"""
function get_pheno_data(file, datatype; transposed=true)

    #first column is individual ID such as : BXD1 , need to be removed.
    pheno = readdlm(file, ','; skipstart=1)[:, 2:end]

    if pheno[1,end] == "f" || pheno[1,end] == "m"
        @info "Removing sex column of phenotype. "
        pheno = pheno[:, 1:end-1]

    end

    pheno = convert(Array{datatype,2}, pheno)

    if transposed
        return transpose(pheno) |> collect
    else
        return pheno
    end

end

"""
$(SIGNATURES)
Filter genotype data with a minor allele frequency threshold.

# Arguments
- `genotype` : genotype matrix
- `maf_threshold` : default value is 0.05.

# Output:
returns filtered genotype matrix.

"""
function filter_maf(genotype::Array{<:Real, 2}; maf_threshold=0.05)
    alleles = 2

    af = sum(genotype, dims=1) ./ (alleles * size(genotype, 1))
    maf = replace(x -> x > 0.5 ? 1-x : x, af)

    if maf_threshold > 0
        mask = vec(maf .>= maf_threshold)
        if size(mask,1) != size(genotype, 2)
            error("Mask dimention does not match original matrix. Mask size: $(size(mask)), Matrix size: $(size(genotype))")
        end
        newgenotype = genotype[:, mask]
    end
    return newgenotype
end


text:
isdefined(Base, :__precompile__) && __precompile__()

module PowerModels

using JSON
using MathProgBase
using JuMP
using Compat
using Memento

import Compat: @__MODULE__

# Create our module level logger (this will get precompiled)
const LOGGER = getlogger(@__MODULE__)

# Register the module level logger at runtime so that folks can access the logger via `getlogger(PowerModels)`
# NOTE: If this line is not included then the precompiled `PowerModels.LOGGER` won't be registered at runtime.
__init__() = Memento.register(LOGGER)

include("io/matlab.jl")
include("io/matpower.jl")
include("io/json.jl")
include("io/common.jl")
include("io/pti.jl")

include("core/base.jl")
include("core/data.jl")
include("core/variable.jl")
include("core/constraint_template.jl")
include("core/constraint.jl")
include("core/relaxation_scheme.jl")
include("core/objective.jl")
include("core/solution.jl")

include("form/acp.jl")
include("form/acr.jl")
include("form/act.jl")
include("form/dcp.jl")
include("form/wr.jl")
include("form/wrm.jl")
include("form/shared.jl")

include("prob/pf.jl")
include("prob/opf.jl")
include("prob/ots.jl")
include("prob/tnep.jl")

end


text:
using Test
using Pitchjx
using DataFrames

# test for no game day.
empty = DataFrame(
  date = String[],
  pitcherid = Int64[],
  pitcher_teamid = Int64[],
  pitcher_firstname = String[],
  pitcher_lastname = String[],
  pitcher_teamname = String[],
  pitcherthrow = String[],
  batterid = Int64[],
  batter_teamid = Int64[],
  batter_firstname = String[],
  batter_lastname = String[],
  batter_teamname = String[],
  batterstand = String[],
  eventdesc = String[],
  pitchresult = String[],
  x = Float64[],
  y = Float64[],
  px = Float64[],
  pz = Float64[],
  pfxx = Float64[],
  pfxz = Float64[],
  zone = Int64[],
  sztop = Float64[],
  szbottom = Float64[],
  pitchtype = String[],
  startspeed = Float64[],
  endspeed = Float64[],
  spindir = Float64[],
  spinrate = Float64[],
  nasty = Float64[]
)

@test pitchjx("2018-01-01") == empty