text:
export BrickletVoltageCurrentConfiguration
struct BrickletVoltageCurrentConfiguration
    averaging::Integer
    voltage_conversion_time::Integer
    current_conversion_time::Integer
end

export BrickletVoltageCurrentCalibration
struct BrickletVoltageCurrentCalibration
    gain_multiplier::Integer
    gain_divisor::Integer
end

export BrickletVoltageCurrentCurrentCallbackThreshold
struct BrickletVoltageCurrentCurrentCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletVoltageCurrentVoltageCallbackThreshold
struct BrickletVoltageCurrentVoltageCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletVoltageCurrentPowerCallbackThreshold
struct BrickletVoltageCurrentPowerCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

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

export BrickletVoltageCurrent
"""
Measures power, DC voltage and DC current up to 720W/36V/20A
"""
mutable struct BrickletVoltageCurrent <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

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

        return new(ipcon, deviceInternal)
    end
end

export get_current
"""
    $(SIGNATURES)

Returns the current.

If you want to get the current periodically, it is recommended to use the
:cb:`Current` callback and set the period with
:func:`Set Current Callback Period`.
"""
function get_current(device::BrickletVoltageCurrent)
    return device.deviceInternal.get_current()
end

export get_voltage
"""
    $(SIGNATURES)

Returns the voltage.

If you want to get the voltage periodically, it is recommended to use the
:cb:`Voltage` callback and set the period with
:func:`Set Voltage Callback Period`.
"""
function get_voltage(device::BrickletVoltageCurrent)
    return device.deviceInternal.get_voltage()
end

export get_power
"""
    $(SIGNATURES)

Returns the power.

If you want to get the power periodically, it is recommended to use the
:cb:`Power` callback and set the period with
:func:`Set Power Callback Period`.
"""
function get_power(device::BrickletVoltageCurrent)
    return device.deviceInternal.get_power()
end

export set_configuration
"""
    $(SIGNATURES)

Sets the configuration of the Voltage/Current Bricklet. It is
possible to configure number of averages as well as
voltage and current conversion time.
"""
function set_configuration(device::BrickletVoltageCurrent, averaging, voltage_conversion_time, current_conversion_time)
    device.deviceInternal.set_configuration(averaging, voltage_conversion_time, current_conversion_time)
end

export get_configuration
"""
    $(SIGNATURES)

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

export set_calibration
"""
    $(SIGNATURES)

Since the shunt resistor that is used to measure the current is not
perfectly precise, it needs to be calibrated by a multiplier and
divisor if a very precise reading is needed.

For example, if you are expecting a measurement of 1000mA and you
are measuring 1023mA, you can calibrate the Voltage/Current Bricklet
by setting the multiplier to 1000 and the divisor to 1023.
"""
function set_calibration(device::BrickletVoltageCurrent, gain_multiplier, gain_divisor)
    device.deviceInternal.set_calibration(gain_multiplier, gain_divisor)
end

export get_calibration
"""
    $(SIGNATURES)

Returns the calibration as set by :func:`Set Calibration`.
"""
function get_calibration(device::BrickletVoltageCurrent)
    return device.deviceInternal.get_calibration()
end

export set_current_callback_period
"""
    $(SIGNATURES)

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

The :cb:`Current` callback is only triggered if the current has changed since
the last triggering.
"""
function set_current_callback_period(device::BrickletVoltageCurrent, period)
    device.deviceInternal.set_current_callback_period(period)
end

export get_current_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Current Callback Period`.
"""
function get_current_callback_period(device::BrickletVoltageCurrent)
    return device.deviceInternal.get_current_callback_period()
end

export set_voltage_callback_period
"""
    $(SIGNATURES)

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

The :cb:`Voltage` callback is only triggered if the voltage has changed since
the last triggering.
"""
function set_voltage_callback_period(device::BrickletVoltageCurrent, period)
    device.deviceInternal.set_voltage_callback_period(period)
end

export get_voltage_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Voltage Callback Period`.
"""
function get_voltage_callback_period(device::BrickletVoltageCurrent)
    return device.deviceInternal.get_voltage_callback_period()
end

export set_power_callback_period
"""
    $(SIGNATURES)

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

The :cb:`Power` callback is only triggered if the power has changed since the
last triggering.
"""
function set_power_callback_period(device::BrickletVoltageCurrent, period)
    device.deviceInternal.set_power_callback_period(period)
end

export get_power_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Get Power Callback Period`.
"""
function get_power_callback_period(device::BrickletVoltageCurrent)
    return device.deviceInternal.get_power_callback_period()
end

export set_current_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Current Reached` callback.

The following options are possible:

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

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

export get_current_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Current Callback Threshold`.
"""
function get_current_callback_threshold(device::BrickletVoltageCurrent)
    return device.deviceInternal.get_current_callback_threshold()
end

export set_voltage_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Voltage Reached` callback.

The following options are possible:

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

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

export get_voltage_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Voltage Callback Threshold`.
"""
function get_voltage_callback_threshold(device::BrickletVoltageCurrent)
    return device.deviceInternal.get_voltage_callback_threshold()
end

export set_power_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Power Reached` callback.

The following options are possible:

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

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

export get_power_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Power Callback Threshold`.
"""
function get_power_callback_threshold(device::BrickletVoltageCurrent)
    return device.deviceInternal.get_power_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callbacks

* :cb:`Current Reached`,
* :cb:`Voltage Reached`,
* :cb:`Power Reached`

are triggered, if the thresholds

* :func:`Set Current Callback Threshold`,
* :func:`Set Voltage Callback Threshold`,
* :func:`Set Power Callback Threshold`

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

export get_debounce_period
"""
    $(SIGNATURES)

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

export get_identity
"""
    $(SIGNATURES)

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

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

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

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


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

function _add_ramp_eqs!(
    model::JuMP.Model,
    g::Unit,
    formulation_prod_vars::Gar1962.ProdVars,
    formulation_ramping::MorLatRam2013.Ramping,
    formulation_status_vars::Gar1962.StatusVars,
)::Nothing
    # TODO: Move upper case constants to model[:instance]
    RESERVES_WHEN_START_UP = true
    RESERVES_WHEN_RAMP_UP = true
    RESERVES_WHEN_RAMP_DOWN = true
    RESERVES_WHEN_SHUT_DOWN = true
    is_initially_on = (g.initial_status > 0)
    SU = g.startup_limit
    SD = g.shutdown_limit
    RU = g.ramp_up_limit
    RD = g.ramp_down_limit
    gn = g.name
    eq_ramp_down = _init(model, :eq_ramp_down)
    eq_ramp_up = _init(model, :eq_str_ramp_up)
    reserve = model[:reserve]

    # Gar1962.ProdVars
    prod_above = model[:prod_above]

    # Gar1962.StatusVars
    is_on = model[:is_on]
    switch_off = model[:switch_off]
    switch_on = model[:switch_on]

    for t in 1:model[:instance].time
        time_invariant =
            (t > 1) ? (abs(g.min_power[t] - g.min_power[t-1]) < 1e-7) : true

        # Ramp up limit
        if t == 1
            if is_initially_on
                eq_ramp_up[gn, t] = @constraint(
                    model,
                    g.min_power[t] +
                    prod_above[gn, t] +
                    (RESERVES_WHEN_RAMP_UP ? reserve[gn, t] : 0.0) <=
                    g.initial_power + RU
                )
            end
        else
            # amk: without accounting for time-varying min power terms,
            #      we might get an infeasible schedule, e.g. if min_power[t-1] = 0, min_power[t] = 10
            #      and ramp_up_limit = 5, the constraint (p'(t) + r(t) <= p'(t-1) + RU)
            #      would be satisfied with p'(t) = r(t) = p'(t-1) = 0
            #      Note that if switch_on[t] = 1, then eqns (20) or (21) go into effect
            if !time_invariant
                # Use equation (24) instead
                SU = g.startup_limit
                max_prod_this_period =
                    g.min_power[t] * is_on[gn, t] +
                    prod_above[gn, t] +
                    (
                        RESERVES_WHEN_START_UP || RESERVES_WHEN_RAMP_UP ?
                        reserve[gn, t] : 0.0
                    )
                min_prod_last_period =
                    g.min_power[t-1] * is_on[gn, t-1] + prod_above[gn, t-1]
                eq_ramp_up[gn, t] = @constraint(
                    model,
                    max_prod_this_period - min_prod_last_period <=
                    RU * is_on[gn, t-1] + SU * switch_on[gn, t]
                )
            else
                # Equation (26) in Kneuven et al. (2020)
                # TODO: what if RU < SU? places too stringent upper bound
                # prod_above[gn, t] when starting up, and creates diff with (24).
                eq_ramp_up[gn, t] = @constraint(
                    model,
                    prod_above[gn, t] +
                    (RESERVES_WHEN_RAMP_UP ? reserve[gn, t] : 0.0) -
                    prod_above[gn, t-1] <= RU
                )
            end
        end

        # Ramp down limit
        if t == 1
            if is_initially_on
                # TODO If RD < SD, or more specifically if
                #        min_power + RD < initial_power < SD
                #      then the generator should be able to shut down at time t = 1,
                #      but the constraint below will force the unit to produce power
                eq_ramp_down[gn, t] = @constraint(
                    model,
                    g.initial_power - (g.min_power[t] + prod_above[gn, t]) <= RD
                )
            end
        else
            # amk: similar to ramp_up, need to account for time-dependent min_power
            if !time_invariant
                # Revert to (25)
                SD = g.shutdown_limit
                max_prod_last_period =
                    g.min_power[t-1] * is_on[gn, t-1] +
                    prod_above[gn, t-1] +
                    (
                        RESERVES_WHEN_SHUT_DOWN || RESERVES_WHEN_RAMP_DOWN ?
                        reserve[gn, t-1] : 0.0
                    )
                min_prod_this_period =
                    g.min_power[t] * is_on[gn, t] + prod_above[gn, t]
                eq_ramp_down[gn, t] = @constraint(
                    model,
                    max_prod_last_period - min_prod_this_period <=
                    RD * is_on[gn, t] + SD * switch_off[gn, t]
                )
            else
                # Equation (27) in Kneuven et al. (2020)
                # TODO: Similar to above, what to do if shutting down in time t
                # and RD < SD? There is a difference with (25).
                eq_ramp_down[gn, t] = @constraint(
                    model,
                    prod_above[gn, t-1] +
                    (RESERVES_WHEN_RAMP_DOWN ? reserve[gn, t-1] : 0.0) -
                    prod_above[gn, t] <= RD
                )
            end
        end
    end
end

ext: jl


text:
# ---
# title: 1033. Moving Stones Until Consecutive
# id: problem1033
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Brainteaser
# link: <https://leetcode.com/problems/moving-stones-until-consecutive/description/>
# hidden: true
# ---
#
# Three stones are on a number line at positions `a`, `b`, and `c`.
#
# Each turn, you pick up a stone at an endpoint (ie., either the lowest or
# highest position stone), and move it to an unoccupied position between those
# endpoints.  Formally, let's say the stones are currently at positions `x, y,
# z` with `x < y < z`.  You pick up the stone at either position `x` or position
# `z`, and move that stone to an integer position `k`, with `x < k < z` and `k
# != y`.
#
# The game ends when you cannot make any more moves, ie. the stones are in
# consecutive positions.
#
# When the game ends, what is the minimum and maximum number of moves that you
# could have made?  Return the answer as an length 2 array: `answer =
# [minimum_moves, maximum_moves]`
#
#
#
# **Example 1:**
#
#
#
#     Input: a = 1, b = 2, c = 5
#     Output: [1,2]
#     Explanation: Move the stone from 5 to 3, or move the stone from 5 to 4 to 3.
#
#
# **Example 2:**
#
#
#
#     Input: a = 4, b = 3, c = 2
#     Output: [0,0]
#     Explanation: We cannot make any moves.
#
#
# **Example 3:**
#
#
#
#     Input: a = 3, b = 5, c = 1
#     Output: [1,2]
#     Explanation: Move the stone from 1 to 4; or move the stone from 1 to 2 to 4.
#
#
#
#
# **Note:**
#
#   1. `1 <= a <= 100`
#   2. `1 <= b <= 100`
#   3. `1 <= c <= 100`
#   4. `a != b, b != c, c != a`
#
#
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# reach, varying block size, fixed delta, all variables
using Reachability

reach_kD_all(model::String, create_plots::Bool=false) =
    reach_kD_all([model], create_plots)

function reach_kD_all(models::Vector{String}, create_plots::Bool=false)
    # load models
    models_loaded = false
    for model in models
        if !isdefined(Main, Symbol(model))
            if !models_loaded
                model_library = include(@relpath "model_library.jl")
            end
            path = (@relpath "../") * model_library[model]
            include(path)
            models_loaded = true
        end
    end
    plots_loaded = false
    if create_plots && !isdefined(Main, :plot_reach)
        include(@relpath "plot_reach.jl")
        plots_loaded = true
    end
    if models_loaded
        println("models have been loaded; please run again for analysis")
        return
    elseif plots_loaded
        println("plot backend has been loaded; please run again for analysis")
        return
    end

    println("-- benchmark suite 'reach_kD_all' --")

    for model in models
        println("- analyzing model '$model' -")

        # raw model (check mode to obtain the property)
        func = getfield(Main, Symbol(model))
        S, options_raw = func(:mode => "check")
        n = MathematicalSystems.statedim(S)

        # extract output function from property
        property = options_raw[:property]
        if !(property isa LinearConstraintProperty)
            println("cannot construct output function from properties of " *
                "type $(typeof(property))... skipping model")
            break
        end
        projection_matrix = reshape(property.clauses[1].atoms[1].a, (1, n))

        # options
        dict_raw = options_raw.dict
        dict_raw[:mode] = "reach"
        dict_raw[:verbosity] = "info"
        dict_raw[:δ] = 1e-3
        dict_raw[:vars] = 1:n
        dict_raw[:plot_vars] = [0, n]
        dict_raw[:lazy_X0] = true
        dict_raw[:lazy_inputs_interval] = -1
        dict_raw[:project_reachset] = false
        dict_raw[:projection_matrix] = projection_matrix

        # create uniform partitions
        partitions = Vector{Vector{AbstractVector{Int}}}(ceil(Int, log(2, n)) + 1)
        k = 0
        while true
            m = 2^k
            k += 1
            if m >= n
                m = n
            end
            uniform = [m*i-(m-1):m*i for i in 1:div(n, m)]
            partitions[k] = n%m == 0 ? uniform : vcat(uniform, [div(n, m) *m+1:n])
            if m == n
                break
            end
        end

        for i in 1:(length(partitions) + 1)
            dict = copy(dict_raw)
            if i == 1
                # warm-up run
                dict[:N] = 3
                dict_raw[:partition] = partitions[1]
            else
                # benchmark settings
                k = i == length(partitions) + 1 ? n : 2^(i-2)
                dict[:N] = 20
                dict[:logfile] = "$model-reach-$(k)D-varying-fixedstep-allvars.txt"
                dict_raw[:partition] = partitions[i-1]
            end
            result = solve(S, Options(dict))
            if create_plots && i > 1
                plot_reach(result, "$model-$(k)D-varying-all")
            end
        end
        println()
    end
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
# Types and functions related to (sets of) points

"Set of points with indices."
struct PointSet
    points::Array{Point}          # coordinates, insertion order
    indices::Dict{Point, Int64}   # find index from point
end
PointSet() = PointSet([], Dict())

length(set::PointSet) = length(set.points)

"Add single point to set."
function push!(set::PointSet, point::Point)
    if !haskey(set.indices, point)
        push!(set.points, point)
        set.indices[point] = length(set.points)
    end
    return set
end

"Add single point (from position) to set."
push!(set::PointSet, pos::Position) = push!(set, Point(pos))

"Add all points of a geometry to set recursively."
function push!(set::PointSet, coords::Vector{T}) where T
    for c in coords
        push!(set, c)
    end
    return set
end

"Add all points of a geometry to set recursively."
push!(set::PointSet, geo::AbstractGeometry) = push!(set, coordinates(geo))

"Get index of given point."
getindex(set::PointSet, point::Point) = set.indices[point]

"Get index for a position."
getindex(set::PointSet, pos::Position) = set.indices[Point(pos)]

"Get point at index."
getindex(set::PointSet, idx::Int64) = set.points[idx]


text:
module MutableStringTestAscii
using MutableStrings

const nstr = 10^5
const lstr = 6

function makestr()
    iob = IOBuffer()
    for idx in 1:nstr
        write(iob, randstring(lstr))
        write(iob, ' ')
    end
    takebuf_string(iob)
end

#mk_blank(x::MutableASCIIString, y::Range1{Int64}) = (x.data[y] = uint8(' '))
mk_blank(x) = ""
mk_space(x) = " "

function warmup()
    _mstr = MutableASCIIString(copy("asdsdsf".data))
    l = length(_mstr)
    lowercase!(_mstr);
    uppercase!(_mstr);
    reverse!(_mstr);
    ucfirst!(_mstr);
    lcfirst!(_mstr);
    replace!(_mstr, "abc", mk_blank);
    replace!(_mstr, r"bcd", mk_blank);
    rx = r"c[a-z]+\s"
    search(_mstr, rx)
    search(_mstr, "abc")
    matchall(rx, _mstr)
end

function test_set_get()
    ms = MutableASCIIString("Lorem ipsum dolor sit amet, consectetur adipisicing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat. Duis aute irure dolor in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur. Excepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mollit anim id est laborum.")
    @assert ms[1] == 'L'
    @assert ms[1:10] == "Lorem ipsu"
    @assert ms[2:10] == "orem ipsu"
    @assert typeof(ms[1:10]) == ASCIIString
    ms[1] = 'X'
    @assert ms[1:10] == "Xorem ipsu"
    ms[1:10] = 'X'
    @assert ms[1:10] == "X"^10
    ms[1:10] = "Lorem ipsu"
    @assert ms[1:10] == "Lorem ipsu"

    replace!(ms, "Lorem", "O RUM")
    @assert ms[1:5] == "O RUM"
    replace!(ms, r"O RUM", "Lorem")
    @assert ms[1:10] == "Lorem ipsu"

    s = replace(ms, "Lorem", "O RUM")
    @assert ms[1:5] != "O RUM"
    @assert beginswith(s, "O RUM")
end

function test_time()
    master_str = makestr();

    str = master_str;
    mstr = MutableASCIIString(str);
    @printf("%25s%30s%40s\n", "", "string", "mutable string")
    @printf("%25s%30s%40s\n", "length", "$(length(str))", "$(length(mstr))")
    @printf("%25s%20s%20s%20s%20s\n", "", "time", "bytes", "time", "bytes")

    gc(); _ret, t, b = @timed str = lowercase(str);
    gc(); _mret, mt, mb = @timed lowercase!(mstr);
    @printf("%25s%20s%20s%20s%20s\n", "lowercase", "$t", "$b", "$mt", "$mb")

    gc(); _ret, t, b = @timed str = reverse(str);
    gc(); _mret, mt, mb = @timed reverse!(mstr);
    @printf("%25s%20s%20s%20s%20s\n", "reverse", "$t", "$b", "$mt", "$mb")

    gc(); _ret, t, b = @timed str = replace(str, "ba", "  ");
    gc(); _mret, mt, mb = @timed replace!(mstr, "ba", "  ");
    @printf("%25s%20s%20s%20s%20s\n", "repl spaces", "$t", "$b", "$mt", "$mb")

    gc(); _ret, t, b = @timed str = replace(str, "cd", "c");
    gc(); _mret, mt, mb = @timed replace!(mstr, "cd", "c");
    @printf("%25s%20s%20s%20s%20s\n", "repl str", "$t", "$b", "$mt", "$mb")

    gc(); _ret, t, b = @timed str = replace(str, "cd", 'c');
    gc(); _mret, mt, mb = @timed replace!(mstr, "cd", 'c');
    @printf("%25s%20s%20s%20s%20s\n", "repl char", "$t", "$b", "$mt", "$mb")

    gc(); _ret, t, b = @timed str = replace(str, "ab", "");
    gc(); _mret, mt, mb = @timed replace!(mstr, "ab", "");
    @printf("%25s%20s%20s%20s%20s\n", "repl nothing", "$t", "$b", "$mt", "$mb")

    rx = r"c[a-z]+\s"
    gc(); _ret, t, b = @timed matchall(rx, str);
    gc(); _mret, mt, mb = @timed matchall(rx, mstr);
    @printf("%25s%20s%20s%20s%20s\n", "matchall", "$t", "$b", "$mt", "$mb")

    str = lowercase(master_str);
    mstr = MutableASCIIString(str);
    _cret = _cmret = 0
    gc(); _ret, t, b = @timed begin
        pos = 1
        while pos <= length(str)
            res = search(str, rx, pos)
            (res.start == 0) && break
            pos = res.start + length(res)
            _cret += 1
        end
    end
    gc(); _mret, mt, mb = @timed begin
        pos = 1
        while pos <= length(mstr)
            res = search(mstr, rx, pos)
            (res.start == 0) && break
            pos = res.start + length(res)
            _cmret += 1
        end
    end
    @assert _cret == _cmret
    @printf("%25s%20s%20s%20s%20s\n", "searched $(_cret) times", "$t", "$b", "$mt", "$mb")

    gc(); _ret, t, b = @timed replace(str, rx, "")
    gc(); _mret, mt, mb = @timed replace!(mstr, rx, mk_blank)
    @printf("%25s%20s%20s%20s%20s\n", "repl regex nothing fn", "$t", "$b", "$mt", "$mb")

    rx = Regex(str[1:1] * "[a-z]+\\s")
    gc(); _ret, t, b = @timed replace(str, rx, " ")
    gc(); _mret, mt, mb = @timed replace!(mstr, rx, mk_space)
    @printf("%25s%20s%20s%20s%20s\n", "repl regex space fn", "$t", "$b", "$mt", "$mb")

    str = lowercase(master_str);
    mstr = MutableASCIIString(str);
    l1 = length(mstr)
    rx = Regex(str[1:1] * "[a-z]+\\s")
    gc(); _ret, t, b = @timed replace(str, rx, "")
    gc(); _mret, mt, mb = @timed replace!(mstr, rx, "")
    @assert length(mstr) < l1
    @printf("%25s%20s%20s%20s%20s\n", "repl regex nothing", "$t", "$b", "$mt", "$mb")

    str = lowercase(master_str);
    mstr = MutableASCIIString(str);
    l1 = length(mstr)
    rx = Regex(str[1:1] * "[a-z]+\\s")
    gc(); _ret, t, b = @timed replace(str, rx, " ")
    gc(); _mret, mt, mb = @timed replace!(mstr, rx, " ")
    @assert length(mstr) < l1
    @printf("%25s%20s%20s%20s%20s\n", "repl regex space", "$t", "$b", "$mt", "$mb")

    str = lowercase(master_str);
    mstr = MutableASCIIString(str);
    l1 = length(mstr)
    rx = Regex(str[1:1] * "[a-z]+\\s")
    gc(); _ret, t, b = @timed replace(str, rx, ' ')
    gc(); _mret, mt, mb = @timed replace!(mstr, rx, ' ')
    @assert length(mstr) == l1
    @printf("%25s%20s%20s%20s%20s\n", "repl regex char", "$t", "$b", "$mt", "$mb")
end
end # module MutableStringTestAscii

MutableStringTestAscii.warmup()
MutableStringTestAscii.test_set_get()
MutableStringTestAscii.test_time()


text:
################################################################################
#
#   FUNCTIONS FOR CONSTRUCTION OF K SPACE POINTS
#
#   FILE CONTAINS
#   - general building of points
#   - high symmetry mesh
#
#   ALL FUNCTIONS ARE NOT EXPORTED
#
################################################################################

################################################################################
#
#   FUNCTIONS FOR GENERAL CONSTRUCTION OF K SPACE POINTS
#
################################################################################

# contruct line between vertices with a certain resolution
function getPointsOnLine(
            point1 :: Vector{<:Real},
            point2 :: Vector{<:Real},
            resolution :: Integer
        ) :: Vector{Vector{Float64}}

    # the array of points which is later returned
    line = Vector{Vector{Float64}}(undef, resolution)
    # the multipliers between the two points
    alphas = range(0.0, stop=1.0, length=resolution)

    # set all line points
    for i in 1:resolution
        line[i] = (point1.*(1-alphas[i])) .+ (point2.*alphas[i])
    end

    # return the array of points
    return line
end
function getPointsOnLine(
            point1 :: P1,
            point2 :: P2,
            resolution :: Integer
        ) :: Vector{Vector{Float64}} where {D,P1<:AbstractReciprocalPoint{D}, P2<:AbstractReciprocalPoint{D}}

    # parse to k space points
    return getPointsOnLine(point(point1), point(point2), resolution)
end
export getPointsOnLine

# get random points inside a box
function getPointsRandomInBox(
            N           :: Integer,
            center      :: Vector{<:Real},
            dimensions  :: Vector{<:Real}
        )

    # the list of points to return later
    points = Vector{Vector{Float64}}(undef, N)
    # the dimension of points
    d = length(dimensions)

    # set all points
    for i in 1:N
        # set the point to a random point
        points[i] = rand(Float64,d)
        # shift the point into the box
        points[i] .*= 2
        points[i] .-= 1
        points[i] .*= dimensions
        points[i] .+= center
    end

    # return the array of points
    return points
end
export getPointsRandomInBox

# get random points inside a convex hull of some points
function getPointsRandomInConvexHull(
            N             :: Integer,
            anchor_points :: Vector{<:Real}...
        )

    # the list of points to return later
    points = Vector{Vector{Float64}}(undef, N)

    # set all points
    for i in 1:N
        # set the point to a random point
        alphas = rand(Float64,length(anchor_points))
        alphas ./= sum(alphas)
        # shift the point into the box
        points[i] = sum(anchor_points .* alphas)
    end

    # return the array of points
    return points
end
export getPointsRandomInConvexHull

# compute the center as half the diagonal from a set of vertices forming a BZ surface
function getPointCenter(
            points :: Vector{<:Vector{<:Real}}
        ) :: Vector{Float64}

    # return sum of all points and divide by number of points
    return sum(points) ./ length(points)
end
function getPointCenter(
            points :: Vector{<:Real} ...
        ) :: Vector{Float64}

    # return sum of all points and divide by number of points
    return sum(points) ./ length(points)
end
export getPointCenter


text:
# Power grid/braess paradox
# From https://link.springer.com/article/10.1140/epjb/e2013-40469-4
# Nonlocal failures in complex supply networks by single link additions
# Dirk Witthaut & Marc Timme
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random, OrdinaryDiffEq

using LinearAlgebra
P = 1.0
K = 1.03*P
α = 1.0*P
# Coefficient of power of the nodes, according to figure
power = Float64[-1, +1, +1, +1, -1, -1, +1, -1]
N = length(power) # number of nodes
# Weighted adjacency matrix, in units of K:
links = zeros(N, N)
links[1,2] = links[1,6] = links[1,7] = 1
links[2,1] = links[2,3] = 1
links[3,2] = links[3,4] = links[3,7] = 1
links[4,3] = links[4,5] = links[4,8] = 1
links[5,4] = links[5,6] = 1
links[6,5] = links[6,1] = links[6,8] = 1
links[7,1] = links[7,3] = 1
links[8,4] = links[8,6] = 1
@assert issymmetric(links)
# Colors of consumer/producer
PLUS = "C1"
MINUS = "C3"

function power_grid_f!(du, u, p, t)
    P, K, α, power, links = p
    N = size(u, 1)
    # φ is first column, ω is second column of `u`
    @inbounds for i in 1:N
        du[i,1] = u[i,2] # φdot = ω
        du[i,2] = power[i]*P - α*u[i,2] +
            sum(K*links[i,j]*sin(u[j,1] - u[i,1]) for j in 1:N)
    end
end

u0 = zeros(N,2)
T = 20.0
saveat = 0:0.1:T

function rungrid!(ax, links)
    p = (P, K, α, power, links)
    pg = ODEProblem(power_grid_f!, u0, (0.0, T), p)
    sol = solve(pg, Tsit5(); saveat)
    # Phases timeseries
    φs = []
    for i in 1:N
        push!(φs, [u[i,1] for u in sol.u])
    end
    for i in 1:N
        c = power[i] > 0 ? PLUS : MINUS
        # ax.plot(sol.t, (φs[i] .- φs[1]), alpha = 0.9, color = c, lw = 2)
        ax.plot(sol.t, φs[i], alpha = 0.9, color = c, lw = 2)
    end
end

fig = figure()
rungrid!(gca(), links)

# %% plotting
# node coordinates
nx = [2.5, 2, 2.5, 0.5, 0, 0.5, 3, 1]
ny = [-1, 0, 1, 1, 0, -1, 0, 0]

fig, axs = subplots(2,3)

function plot_graph!(ax, nx, ny, power, links)
    colors = [p>0 ? "C1" : "C3" for p in power]
    colors = darken_color.(colors, 1.2)
    d = darken_color.(colors)
    ax.scatter(nx, ny; c = colors, s = 600, zorder = 98, edgecolors = d, linewidths = 2)
    for i in 1:N # loop over nodes
        t = power[i] > 0 ? "+" : "-"

        ax.annotate(t, (nx[i], ny[i]), fontsize = 24, ha = :center, va = :center, zorder = 99)
        # ax.annotate(i, (nx[i], ny[i]), fontsize = 24)
        # links
        for j in 1:N
            k = links[i,j]
            k == 0 && continue
            if abs(k) == 1
                color = BLACK
                lw = 3
            else
                color = "C0"
                lw = 9
            end
            ax.plot([nx[i],nx[j]], [ny[i],ny[j]]; color, lw, zorder = 1)
        end
    end
end

links2 = copy(links)
links2[3,4] = links2[4,3] = 2
links3 = copy(links)
links3[4,2] = links3[2,4] = 1
alllinks = (links, links2, links3)

for j in 1:3
    rungrid!(axs[2,j], alllinks[j])
    plot_graph!(axs[1,j], nx, ny, power, alllinks[j])
    axs[1,j].axis("off")
    axs[1,j].set_ylim(-1.4,1.4)
    axs[1,j].set_xlim(-0.4,3.4)
    y = 6
    axs[2,j].set_ylim(-y,y)
    axs[2,j].set_yticks([-y,y])
    if j > 1
        plt.setp(axs[2,j].get_yticklabels(), visible=false)
    end
    axs[2,j].set_xticks(0:6:18)
    axs[2,j].set_xlim(0,18)
    axs[2,j].set_xlabel("\$t\$",labelpad = -20)
    t = j == 1 ? "sync." : "no sync."
    axs[2,j].text(0.05, 0.8, t, transform = axs[2,j].transAxes)
end
axs[1,1].set_title("original")
axs[1,2].set_title("added capacity")
axs[1,3].set_title("added link")
axs[2,1].set_ylabel("\$\\phi_i\$", labelpad = -30)
add_identifiers!(fig, axs[2,:])
fig.tight_layout(pad=0.3)
wsave(plotsdir("10", "powergrid"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using Compat
using ProtoBuf
import ProtoBuf.meta
import Base: hash, isequal, ==

type __enum_XAttrSetFlagProto <: ProtoEnum
    XATTR_CREATE::Int32
    XATTR_REPLACE::Int32
    __enum_XAttrSetFlagProto() = new(1,2)
end #type __enum_XAttrSetFlagProto
const XAttrSetFlagProto = __enum_XAttrSetFlagProto()

type __enum_XAttrProto_XAttrNamespaceProto <: ProtoEnum
    USER::Int32
    TRUSTED::Int32
    SECURITY::Int32
    SYSTEM::Int32
    RAW::Int32
    __enum_XAttrProto_XAttrNamespaceProto() = new(0,1,2,3,4)
end #type __enum_XAttrProto_XAttrNamespaceProto
const XAttrProto_XAttrNamespaceProto = __enum_XAttrProto_XAttrNamespaceProto()

type XAttrProto
    namespace::Int32
    name::AbstractString
    value::Array{UInt8,1}
    XAttrProto() = (o=new(); fillunset(o); o)
end #type XAttrProto
const __req_XAttrProto = Symbol[:namespace,:name]
meta(t::Type{XAttrProto}) = meta(t, __req_XAttrProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::XAttrProto) = ProtoBuf.protohash(v)
isequal(v1::XAttrProto, v2::XAttrProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::XAttrProto, v2::XAttrProto) = ProtoBuf.protoeq(v1, v2)

type XAttrEditLogProto
    src::AbstractString
    xAttrs::Array{XAttrProto,1}
    XAttrEditLogProto() = (o=new(); fillunset(o); o)
end #type XAttrEditLogProto
hash(v::XAttrEditLogProto) = ProtoBuf.protohash(v)
isequal(v1::XAttrEditLogProto, v2::XAttrEditLogProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::XAttrEditLogProto, v2::XAttrEditLogProto) = ProtoBuf.protoeq(v1, v2)

type SetXAttrRequestProto
    src::AbstractString
    xAttr::XAttrProto
    flag::UInt32
    SetXAttrRequestProto() = (o=new(); fillunset(o); o)
end #type SetXAttrRequestProto
const __req_SetXAttrRequestProto = Symbol[:src]
meta(t::Type{SetXAttrRequestProto}) = meta(t, __req_SetXAttrRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::SetXAttrRequestProto) = ProtoBuf.protohash(v)
isequal(v1::SetXAttrRequestProto, v2::SetXAttrRequestProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::SetXAttrRequestProto, v2::SetXAttrRequestProto) = ProtoBuf.protoeq(v1, v2)

type SetXAttrResponseProto
    SetXAttrResponseProto() = (o=new(); fillunset(o); o)
end #type SetXAttrResponseProto
hash(v::SetXAttrResponseProto) = ProtoBuf.protohash(v)
isequal(v1::SetXAttrResponseProto, v2::SetXAttrResponseProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::SetXAttrResponseProto, v2::SetXAttrResponseProto) = ProtoBuf.protoeq(v1, v2)

type GetXAttrsRequestProto
    src::AbstractString
    xAttrs::Array{XAttrProto,1}
    GetXAttrsRequestProto() = (o=new(); fillunset(o); o)
end #type GetXAttrsRequestProto
const __req_GetXAttrsRequestProto = Symbol[:src]
meta(t::Type{GetXAttrsRequestProto}) = meta(t, __req_GetXAttrsRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::GetXAttrsRequestProto) = ProtoBuf.protohash(v)
isequal(v1::GetXAttrsRequestProto, v2::GetXAttrsRequestProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::GetXAttrsRequestProto, v2::GetXAttrsRequestProto) = ProtoBuf.protoeq(v1, v2)

type GetXAttrsResponseProto
    xAttrs::Array{XAttrProto,1}
    GetXAttrsResponseProto() = (o=new(); fillunset(o); o)
end #type GetXAttrsResponseProto
hash(v::GetXAttrsResponseProto) = ProtoBuf.protohash(v)
isequal(v1::GetXAttrsResponseProto, v2::GetXAttrsResponseProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::GetXAttrsResponseProto, v2::GetXAttrsResponseProto) = ProtoBuf.protoeq(v1, v2)

type ListXAttrsRequestProto
    src::AbstractString
    ListXAttrsRequestProto() = (o=new(); fillunset(o); o)
end #type ListXAttrsRequestProto
const __req_ListXAttrsRequestProto = Symbol[:src]
meta(t::Type{ListXAttrsRequestProto}) = meta(t, __req_ListXAttrsRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::ListXAttrsRequestProto) = ProtoBuf.protohash(v)
isequal(v1::ListXAttrsRequestProto, v2::ListXAttrsRequestProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::ListXAttrsRequestProto, v2::ListXAttrsRequestProto) = ProtoBuf.protoeq(v1, v2)

type ListXAttrsResponseProto
    xAttrs::Array{XAttrProto,1}
    ListXAttrsResponseProto() = (o=new(); fillunset(o); o)
end #type ListXAttrsResponseProto
hash(v::ListXAttrsResponseProto) = ProtoBuf.protohash(v)
isequal(v1::ListXAttrsResponseProto, v2::ListXAttrsResponseProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::ListXAttrsResponseProto, v2::ListXAttrsResponseProto) = ProtoBuf.protoeq(v1, v2)

type RemoveXAttrRequestProto
    src::AbstractString
    xAttr::XAttrProto
    RemoveXAttrRequestProto() = (o=new(); fillunset(o); o)
end #type RemoveXAttrRequestProto
const __req_RemoveXAttrRequestProto = Symbol[:src]
meta(t::Type{RemoveXAttrRequestProto}) = meta(t, __req_RemoveXAttrRequestProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::RemoveXAttrRequestProto) = ProtoBuf.protohash(v)
isequal(v1::RemoveXAttrRequestProto, v2::RemoveXAttrRequestProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::RemoveXAttrRequestProto, v2::RemoveXAttrRequestProto) = ProtoBuf.protoeq(v1, v2)

type RemoveXAttrResponseProto
    RemoveXAttrResponseProto() = (o=new(); fillunset(o); o)
end #type RemoveXAttrResponseProto
hash(v::RemoveXAttrResponseProto) = ProtoBuf.protohash(v)
isequal(v1::RemoveXAttrResponseProto, v2::RemoveXAttrResponseProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::RemoveXAttrResponseProto, v2::RemoveXAttrResponseProto) = ProtoBuf.protoeq(v1, v2)

export XAttrSetFlagProto, XAttrProto_XAttrNamespaceProto, XAttrProto, XAttrEditLogProto, SetXAttrRequestProto, SetXAttrResponseProto, GetXAttrsRequestProto, GetXAttrsResponseProto, ListXAttrsRequestProto, ListXAttrsResponseProto, RemoveXAttrRequestProto, RemoveXAttrResponseProto


text:
import ForwardDiff
import Optim

@inline function geometric_median(a::AbstractVector{V};
                                  x0 = zeros(T, length(first(a)))) where V <: AbstractVector{T} where T
    function _sum_of_distances(x)
        return sum(l2_distance.(a, Ref(x)))
    end
    result = Optim.optimize(_sum_of_distances,
                            x0,
                            Optim.LBFGS(); autodiff = :forward)
    minimizer = Optim.minimizer(result)
    return minimizer
end


text:
include("data_types.jl")
include("data_tree.jl")

include("util.jl")

include("point_cloud.jl")
include("integration_volume.jl")

include("hyper_rectangle.jl")
include("whitening_transformation.jl")
include("harmonic_mean_integration.jl")
include("hm_integration_rectangle.jl")

include("uncertainty.jl")

# ToDo: rename hm_... functions to ahmi_...


text:
using DifferentialEquations
using PyCall
pygui(:qt)
using PyPlot
pygui(true)

function condition(u,t,integrator) # Event when event_f(u,t) == 0
 u[3]
end

function affect!(integrator)
end
#cb1 = SavingCallback((u,t,integrator)->(u[4],u[2]), saved_values)

cb2= ContinuousCallback(condition, affect!)
Energy=0.112012;
H=(2*Energy)^2
h=sqrt(H)
n_iter=10
Y_start=8e-5
Y_end=1e-2
t_end=10000
Y=range(Y_start,stop=Y_end,length=n_iter)
P=zeros(n_iter)
PofY_find(YY)=sqrt((h+YY.^2-h*YY.^2 )./(1+h-h*YY.^2))
H1(p,q,t)= ( (q[2]-q[1])^2+(p[2]-p[1])^2 )*( (q[2]+q[1])^2+(p[2]+p[1])^2 )
H2(p,q,t)= (p[2]^4+2*p[2]^2*(q[1]^2-1)+(1+q[1]^2)^2 )*(q[2]^4+2*q[2]^2*(p[1]^2+1)+(p[1]^2-1)^2 )
HAM(p,q,t)=H1(p,q,t)/H2(p,q,t)
PofY_find(YY)=sqrt((h+YY.^2-h*YY.^2 )./(1+h-h*YY.^2))

for k=1:n_iter

    q0,p0=[zeros(2) for i in 1:2]
    P[k]=PofY_find(Y[k])
    q0[1]=0
    p0[1]=P[k]
    q0[2]=0
    p0[2]=Y[k]
    prob= HamiltonianProblem{true}(HAM, p0, q0, (0., t_end))
    t,A=solve(prob, Vern9(),reltol=1e-13,abstol=1e-15,maxiters=1e15, callback=cb2,save_end=false,save_everystep=false)
    plot(A[2,:],A[4,:],",")
#     plot!(A[2,:],-A[4,:],seriestype=:scatter,marker=:circle,msw = 0,ms=.1,seriescolor = :black,legend=false)
#     plot!(-A[2,:],A[4,:],seriestype=:scatter,marker=:circle,msw = 0,ms=.1,seriescolor = :black,legend=false)
#     plot(-A[2,:],-A[4,:],seriestype=:scatter,marker=:circle,msw = 0,ms=.1,seriescolor = :black,legend=false)
end
 savefig("fn.eps")
# N=length(saved_values.saveval)
# QList=zeros(N)
# PList=zeros(N)
#
#
# for k in 1:N
#     QList[k],PList[k]=saved_values.saveval[k]
# end
#     plot(QList,PList,seriestype=:scatter,marker=:circle,msw = 0,ms=.1)


text:
# Here we calculate the effective wavenumber and effective wave amplitudes, without any restriction to the volume fraction of particles and incident wave frequency.

# include depricated function to find a single effective wavenumber, when in fact there are many. The code is still used in tests and gives many correct results
# include("wavenumber_single.jl")

" Returns all the transmitted effective wavenumbers"
wavenumbers(ω::T, medium::Medium{T}, specie::Specie{T}; kws...) where T<:Number = wavenumbers(ω, medium, [specie]; kws...)

function wavenumbers(ω::T, medium::Medium{T}, species::Vector{Specie{T}};
        num_wavenumbers = 8,
        apply_meshing::Bool = true,
        kws...) where T<:Number

    k_effs = wavenumbers_path(ω, medium, species; num_wavenumbers=num_wavenumbers, kws...)
    num_wavenumbers = min(length(k_effs),num_wavenumbers)

    if num_wavenumbers > 1 && apply_meshing
        num_wavenumbers = min(length(k_effs),num_wavenumbers)
        k_effs = wavenumbers_mesh(ω, k_effs[1:num_wavenumbers], medium, species; kws...)
    end

    return k_effs
end


text:
using Dates
println(today() - Date(ARGS[1]))


text:
using ParamLevelSet
using jInv.Mesh

println("test for SPD regularizer");

plotting = false;
if plotting
	using jInvVis
	using PyPlot
end

n = [64,64,64];
Mesh = getRegularMesh([0.0;5.0;0.0;5.0;0.0;5.0],n);
midmesh = [(Mesh.domain[2]+Mesh.domain[1])/2 ; (Mesh.domain[4]+Mesh.domain[3])/2; (Mesh.domain[5]+Mesh.domain[6])/2];
lenmesh = [(Mesh.domain[2]-Mesh.domain[1]) ; (Mesh.domain[4]-Mesh.domain[3]); (Mesh.domain[6]-Mesh.domain[5])];
numParamOfRBF = 10;
nRBF = 5
m = zeros(nRBF*numParamOfRBF);
dm = zeros(length(m));

for k=1:nRBF
	offset = (k-1)*numParamOfRBF + 1;
	m[offset] = randn();
	A = randn(3,3); A = 0.3*A'*A + 0.5*Matrix(I,3,3);
	m[(offset+1):(offset+6)] = A[idx_trilA()];
	m[(offset+7):(offset+9)] = midmesh + 0.05*lenmesh.*randn(size(lenmesh));
	dm[offset] = 0.01*randn();
	dm[(offset+1):(offset+6)] = 0.01*randn(6);# 0.01*[randn(); randn(); 0.0; randn(); 0.0; randn()];#
	dm[(offset+7):(offset+9)] = 0.01*randn(3)
end

regfun = (m)-> RBF_SPD_regularization(m,0,nRBF);
hh = 1.0;
R0,dR,d2R = regfun(m);
for k=0:8
	hhh = (0.5^k)*hh;
	Rt,dRt = regfun(m+ hhh*dm);
	println("norm(Rt-R0): ",norm(Rt-R0),", norm(Rt - R0 - dR*dm): ",norm(Rt - R0 - dot(dR,hhh*dm)));
end

println("Testing reg Hessian")
for k=0:8
	hhh = (0.5^k)*hh;
	Rt,dRt = regfun(m+ hhh*dm);
	println("norm(dRt-dR0): ",norm(dRt-dR),", norm(dRt - dR - d2R*dm): ",norm(dRt - dR - d2R*(hhh*dm)));
end

## Below are some checks that we had to do along the way. May be ignored.

# regfun = (A)->(-log(det(A)),-inv(A))
# m = randn(3,3); m = 0.3*m'*m + 0.5*eye(3);
# dm = 0.01*randn(3,3); dm = dm + dm';
# R0,dR = regfun(m);
# for k=0:8
	# hhh = (0.5^k);
	# Rt, = regfun(m+ hhh*dm);
	# println("norm(Rt-R0): ",norm(Rt-R0),", norm(Rt - R0 - dR*dm): ",norm(Rt - R0 - dot(dR[:],hhh*dm[:])));
# end

# regfun = (A)->(-inv(A),kron(inv(A),inv(A)))
# m = randn(3,3); m = 0.3*m'*m + 0.5*eye(3);
# dm = 0.01*randn(3,3); dm = dm + dm';
# R0,dR = regfun(m);
# for k=0:8
	# hhh = (0.5^k);
	# Rt, = regfun(m+ hhh*dm);
	# println("norm(Rt-R0): ",vecnorm(Rt-R0),", vecnorm(Rt - R0 - dR*dm): ",vecnorm(Rt[:] - R0[:] - dR*(hhh*dm[:])));
# end

# itrilA = idx_trilA();
# P = zeros(9,6);
# P[itrilA,1:6] = eye(6);
# P[4,2] = 1.0; P[7,3] = 1.0; P[8,5] = 1.0;
# regfun = (a)->(A = [a[1] a[2] a[3] ; a[2] a[4] a[5] ; a[3] a[5] a[6]]; invA = inv(A); return -invA[itrilA],(kron(invA,invA)*P)[itrilA,:]);
# m = randn(3,3); m = 0.3*m'*m + 0.5*eye(3);
# m = m[itrilA];
# dm = 0.01*randn(6);
# R0,dR = regfun(m);
# for k=0:8
	# hhh = (0.5^k);
	# Rt, = regfun(m+ hhh*dm);
	# println("norm(Rt-R0): ",vecnorm(Rt-R0),", vecnorm(Rt - R0 - dR*dm): ",vecnorm(Rt[:] - R0[:] - dR*(hhh*dm[:])));
# end


text:
export ProductMeasure
using Base: eltype

struct ProductMeasure{T} <: AbstractMeasure{T}
    # TODO: Type annotation
    components

    ProductMeasure(μs...) = new{Tuple{eltype.(μs)...}}(μs)
end

# ProductMeasure(m::NTuple{N, Measure{X}}) where {N,X} = ProductMeasure(m...)

Base.length(m::ProductMeasure{T}) where {T} = length(m.components)

function Base.:*(μ::ProductMeasure{X}, ν::ProductMeasure{Y}) where {X,Y}
    components = (μ.components..., ν.components...)
    ProductMeasure(components...)
end

function Base.:*(μ::AbstractMeasure{X}, ν::ProductMeasure{Y}) where {X,Y}
    components = (μ, ν.components...)
    ProductMeasure(components...)
end

function Base.:*(μ::ProductMeasure{X}, ν::AbstractMeasure{Y}) where {X,Y}
    components = (μ.components..., ν)
    ProductMeasure(components...)
end

function Base.:*(μ::AbstractMeasure{X}, ν::AbstractMeasure{Y}) where {X,Y}
    components = (μ, ν)
    ProductMeasure(components...)
end

export ×
function ×(μ::AbstractMeasure{X}, ν::AbstractMeasure{Y}) where {X,Y}
    return μ*ν
end

function Base.rand(μ::ProductMeasure{T}) where T
    return rand.(μ.components)
end


text:
#!/usr/bin/env julia

using EAGODomainReduction

println("Testing Duality-Based Bound Tightening...")
t = @elapsed include("TestDual.jl")
println("done (took $t seconds).")

println("Testing DAG Contractor ...")
t = @elapsed include("TestDAGcntr.jl")
println("done (took $t seconds).")

println("Testing Range Reduction...")
t = @elapsed include("TestRR.jl")
println("done (took $t seconds).")


text:
# Geometric Utilities
# Areas
export triangle_area, circle_area, annulus_area

triangle_area(a, b, c) =
  let s = (a + b + c)/2.0
    sqrt(max(0.0, s*(s - a)*(s - b)*(s - c)))
  end
circle_area(r) = π*r^2
annulus_area(rₒ, rᵢ) = π*(rₒ^2 - rᵢ^2)

project_to_world(surf) =
    transform(surf, inverse_transformation(frame_at(surf, 0, 0)))

#=

Given a poligonal line described by its vertices, we need to compute another
polygonal line that is parallel to the first one.

=#

v_in_v(v0, v1) =
  let v = v0 + v1
    v*dot(v0, v0)/dot(v, v0)
  end

rotated_v(v, alpha) =
  vpol(pol_rho(v), pol_phi(v) + alpha)

centered_rectangle(p0, w, p1) =
  let v0 = p1 - p0,
      v1 = rotated_v(v0, pi/2),
      c = loc_from_o_vx_vy(p0, v0, v1)
    rectangle(c-vy(w/2, c.cs), distance(p0, p1), w)
  end

offset_vertices(ps::Locs, d::Real, closed) =
  let qs = closed ? [ps[end], ps..., ps[1]] : ps,
      vs = map((p0, p1) -> rotated_v(unitized(p1 - p0)*d, pi/2), qs[1:end-1], qs[2:end]),
      ws = map(v_in_v, vs[1:end-1], vs[2:end])
    map(+, ps, closed ? ws : [vs[1], ws..., vs[end]])
  end

offset(path, d::Real) = d == 0 ? path : nonzero_offset(path, d)
nonzero_offset(path::RectangularPath, d::Real) =
  rectangular_path(add_xy(path.corner, d, d), path.dx - 2d, path.dy - 2d)
nonzero_offset(path::OpenPolygonalPath, d::Real) =
  open_polygonal_path(offset_vertices(path.vertices, d, false))
nonzero_offset(path::ClosedPolygonalPath, d::Real) =
  closed_polygonal_path(offset_vertices(path.vertices, d, true))
nonzero_offset(path::CircularPath, d::Real) =
  circular_path(path.center, path.radius - d)
nonzero_offset(path::ArcPath, d::Real) =
  arc_path(path.center, path.radius - d, path.start_angle, path.amplitude)

export offset

# Polygon combination

closest_vertices_indexes(pts1, pts2) =
  # This is a brute force method. There are better algorithms to do this.
  let min_dist = Inf,
      min_i = nothing,
      min_j = nothing
    for (i, pt1) in enumerate(pts1), (j, pt2) in enumerate(pts2)
      let dist = distance(pt1, pt2)
        if dist < min_dist
          min_dist = dist
          min_i = i
          min_j = j
        end
      end
    end
    min_i, min_j
  end

point_in_segment(r, p, q) =
  let pr = r-p,
      pq = q-p,
      rx = cx(pr)/cx(pq),
      ry = cy(pr)/cy(pq)
    isapprox(rx, ry) && isapprox(ry, cz(pr)/cz(pq))
  end

collinear_segments(p1, p2, q1, q2) =
  point_in_segment(q1, p1, p2) &&
  point_in_segment(q2, p1, p2)

collinear_vertices_indexes(pts1, pts2) =
  for (i1, p1) in enumerate(pts1)
    let i2 = i1%length(pts1)+1,
        p2 = pts1[i2]
      for (j1, q1) in enumerate(pts2)
        let j2 = j1%length(pts2)+1,
            q2 = pts2[j2]
          if collinear_segments(p1, p2, q1, q2)
            return (i1, j1)
          end
        end
      end
    end
  end

subtract_polygon_vertices(pts1, pts2) =
  let ij = collinear_vertices_indexes(pts1, pts2)
    isnothing(ij) ?
      inject_polygon_vertices_at_indexes(pts1, pts2, closest_vertices_indexes(pts1, pts2)) :
      splice_polygon_vertices_at_indexes(pts1, pts2, ij)
  end

inject_polygon_vertices_at_indexes(pts1, pts2, (i, j)) =
  [pts1[1:i]..., reverse([pts2[j:end]..., pts2[1:j]...])..., pts1[i:end]...]

export closest_vertices_indexes, inject_polygon_vertices_at_indexes, subtract_polygon_vertices

# Intersection

segments_intersection(p0, p1, p2, p3) =
  let denom = (p3.y - p2.y)*(p1.x - p0.x) - (p3.x - p2.x)*(p1.y - p0.y)
    if denom == 0
      nothing
    else
      let u = ((p3.x - p2.x)*(p0.y - p2.y) - (p3.y - p2.y)*(p0.x - p2.x))/denom,
          v = ((p1.x - p0.x)*(p0.y - p2.y) - (p1.y - p0.y)*(p0.x - p2.x))/denom
        if 0 <= u <= 1 && 0 <= v <= 1
          xy(p0.x + u*(p1.x - p0.x), p0.y + u*(p1.y - p0.y))
        else
          nothing
        end
      end
    end
  end

lines_intersection(p0, p1, p2, p3) =
  let denom = (p3.y - p2.y)*(p1.x - p0.x) - (p3.x - p2.x)*(p1.y - p0.y)
    if denom == 0
      nothing
    else
      let u = ((p3.x - p2.x)*(p0.y - p2.y) - (p3.y - p2.y)*(p0.x - p2.x))/denom,
          v = ((p1.x - p0.x)*(p0.y - p2.y) - (p1.y - p0.y)*(p0.x - p2.x))/denom
        xy(p0.x + u*(p1.x - p0.x), p0.y + u*(p1.y - p0.y))
      end
    end
  end

export epsilon, nearest_point_from_lines, circle_from_three_points
const epsilon = Parameter(1e-8)

nearest_point_from_lines(l0p0::Loc, l0p1::Loc, l1p0::Loc, l1p1::Loc) =
  let u = l0p1-l0p0,
      v = l1p1-l1p0,
      w = l0p0-l1p0,
      a = dot(u, u),
      b = dot(u, v),
      c = dot(v, v),
      d = dot(u, w),
      e = dot(v, w),
      D = a*c-b*b,
      (sc, tc) = D < epsilon() ?
                #almost parallel
                (0.0, b > c ? d/b : e/c) :
                ((b*e-c*d)/D, (a*e-b*d)/D),
      (p0, p1) = (l0p0+u*sc, l1p0+v*tc)
    intermediate_loc(p0, p1)
  end

circle_from_three_points_2d(v0::Loc, v1::Loc, v2::Loc) =
  let v1sv0 = v1-v0,
      v2sv0 = v2-v0,
      v2sv1 = v2-v1,
      v1pv0 = v1+(v0-u0()),
      v2pv0 = v2+(v0-u0()),
      a = v1sv0.x,
      b = v1sv0.y,
      c = v2sv0.x,
      d = v2sv0.y,
      e = a*v1pv0.x+b*v1pv0.y,
      f = c*v2pv0.x+d*v2pv0.y,
      g = 2.0*(a*v2sv1.y-b*v2sv1.x),
      iscolinear = abs(g) < 1e-8,
      (cx, cy, dx, dy) = iscolinear ?
                         let minx = min(v0.x, v1.x, v2.x),
                             miny = min(v0.y, v1.y, v2.y),
                             maxx = max(v0.x, v1.x, v2.x),
                             maxy = max(v0.y, v1.y, v2.y),
                             x = (minx+maxx)/2,
                             y = (miny+maxy)/2
                           (x, y, x-minx, y-miny)
                         end :
                         let x = (d*e-b*f)/g,
                             y = (a*f-c*e)/g
                           (x, y, x-v0.x, y-v0.y)
                         end,
      radius_squared = dx*dx+dy*dy,
      radius = sqrt(radius_squared)
    (xy(cx, cy), radius)
  end

circle_from_three_points(p0::Loc, p1::Loc, p2::Loc) =
  let cs = cs_from_o_vx_vy(p0, p1-p0, p2-p0)
    with(current_cs, cs) do
      c, r = circle_from_three_points_2d(in_cs(p0, cs),
                                         in_cs(p1, cs),
                                         in_cs(p2, cs))
      (c, r)
    end
  end

const collinearity_tolerance = Parameter(1e-2)
# are the three points sufficiently collinear?
collinear_points(p0, pm, p1, epsilon=collinearity_tolerance()) =
  let a = distance(p0, pm),
      b = distance(pm, p1),
      c = distance(p1, p0)
    triangle_area(a, b, c) < epsilon
  end

#=
export sweep_path_with_path
sweep_path_with_path(path, profile) =
  let vertices = in_world.(path_frames(profile)),
      frames = path_frames(path)
      #show_cs.(frames, 0.1)
    surface_grid([[xyz(cx(p), cy(p), cz(p), frame.cs) for p in vertices]
                  for frame in frames],
                 is_closed_path(profile),
                 is_closed_path(path),
                 is_smooth_path(profile),
                 is_smooth_path(path))
  end
=#
export quad_grid, quad_grid_indexes
quad_grid(quad, points, closed_u, closed_v) =
  let pts = in_world.(points),
      si = size(pts, 1),
      sj = size(pts, 2)
    for i in 1:si-1
      for j in 1:sj-1
        quad(pts[i,j], pts[i+1,j], pts[i+1,j+1], pts[i,j+1])
      end
      if closed_v
        quad(pts[i,sj], pts[i+1,sj], pts[i+1,1], pts[i,1])
      end
    end
    if closed_u
      for j in 1:sj-1
        quad(pts[si,j], pts[1,j], pts[si,j+1], pts[si,j+1])
      end
      if closed_v
        quad(pts[si,sj], pts[1,sj], pts[1,1], pts[si,1])
      end
    end
  end

quad_grid_indexes(si, sj, closed_u, closed_v) =
  let idx(i,j) = (i-1)*sj+(j-1),
      idxs = Vector{Int}[],
      quad(a,b,c,d) = (push!(idxs, [a, b, d]); push!(idxs, [d, b, c]))
    for i in 1:si-1
      for j in 1:sj-1
        quad(idx(i,j), idx(i+1,j), idx(i+1,j+1), idx(i,j+1))
      end
      if closed_v
        quad(idx(i,sj), idx(i+1,sj), idx(i+1,1), idx(i,1))
      end
    end
    if closed_u
      for j in 1:sj-1
        quad(idx(si,j), idx(1,j), idx(1,j+1), idx(si,j+1))
      end
      if closed_v
        quad(idx(si,sj), idx(1,sj), idx(1,1), idx(si,1))
      end
    end
    idxs
  end

export illustrate_path
illustrate_path(path) =
  begin
    for (i, v) in enumerate(path_vertices(path))
      sphere(v, 0.01)
      text(string(i), v+vxy(0.05, 0.05), 0.1)
    end
    stroke(path)
  end

illustrate_expr(expr) =
  :(let p = $(esc(expr)); text($(string(expr)), p); p end)

export illustrate
macro illustrate(exprs...)
  :(tuple($([illustrate_expr(expr) for expr in exprs]...)))
end


text:
include("euler/euler.jl")
using .Numbers: is_s_number

function compute(n::Int)::Int
    result = 0
    for i ∈ 2:n if is_s_number(i, i * i) result += i * i end end
    return result
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using SimpleDirectMediaLayer
const SDL = SimpleDirectMediaLayer

chip82key = Dict(0x01 => SDL.SCANCODE_1, 0x02 => SDL.SCANCODE_2, 0x03 => SDL.SCANCODE_3,
                 0x0C => SDL.SCANCODE_4, 0x04 => SDL.SCANCODE_Q, 0x05 => SDL.SCANCODE_W,
                 0x06 => SDL.SCANCODE_E, 0x0D => SDL.SCANCODE_R, 0x07 => SDL.SCANCODE_A,
                 0x08 => SDL.SCANCODE_S, 0x09 => SDL.SCANCODE_D, 0x0E => SDL.SCANCODE_F,
                 0x0A => SDL.SCANCODE_Z, 0x00 => SDL.SCANCODE_X, 0x0B => SDL.SCANCODE_C,
                 0x0F => SDL.SCANCODE_V)

key2chip8 = Dict(b => a for (a,b) in chip82key)


text:
flip(X::Matrix{T}) where {T} = flipdim(flipdim(X, 1), 2)

sum_1(X) = sum(X, dims=1)
sum_2(X) = sum(X, dims=2)

squeeze_sum(A::Array{T,N}, dim::Integer) where {T,N} = squeeze(sum(A, dims=dim), dims=dim)
squeeze_sum_1(A) = squeeze_sum(A, 1)
squeeze_sum_2(A) = squeeze_sum(A, 2)

## constructor

function struct_like(m)
    try
        return typeof(m)()
    catch e
        if e isa MethodError
            println("ERROR: Trying to instantiate mutable type $(typeof(m)), " *
                    "but it doesn't have a default constructor")
        end
        throw(e)
    end
end

function __construct(mem::Dict, dzdm_v::Symbol, m::T; fields...) where T
    if isimmutable(m)
        return __construct_immutable(T; fields...)
    else
        dz!dm = @get_or_create(mem, dzdm_v, struct_like(m))
        for (f, val) in fields
            setfield!(dz!dm, f, val)
        end
        return dz!dm
    end
end

function __construct(m::T; fields...) where T
    if isimmutable(m)
        return __construct_immutable(T; fields...)
    else
        dz!dm = struct_like(m)
        for (f, val) in fields
            setfield!(dz!dm, f, val)
        end
        return dz!dm
    end
end

function __construct_immutable(::Type{T}; fields...) where T
    f2i = Dict((f, i) for (i, f) in enumerate(fieldnames(T)))
    vals = Array{Any}(length(f2i))
    for (f, val) in fields
        vals[f2i[f]] = val
    end
    return T(vals...)
end


text:
#
# Parrot.jl/src/Parrot.jl
#
# By Elijah Mathews (me@elijahmathews.com)
# openpgp4fpr:ac1d3fb1e8a5eb7d14bd587b2932c725055a90d8
#

module Parrot

    #
    # State packages we're using.
    #
    using Flux, MultivariateStats, OnlineStats, Random, Statistics

    #
    # Export really useful things from Parrot.
    #
    export Alsing, TransformPCA, ReconstructPCA, Normalize, Denormalize, convert

    #
    # Include source code files.
    #
    include("layers.jl")
    # include("prospect.jl")
    include("utils.jl")

end


text:
include("run-gf-solver.jl")

using Logging

logger = ConsoleLogger(stdout, Logging.Error)

function run_dim_vs_nondim_study()
    eos_list = [:full_cnga]
    cases = ["GasLib-11", "GasLib-24", "GasLib-40", "GasLib-134"]

    output_folder = "./output/"
    data_folder = "./data/"

    num_runs = 500

    for case in cases
        for eos in eos_list
            for dimensional in [false, true]
                folder = data_folder * case * "/"
                eos_string = (eos == :ideal) ? "ideal" : "cnga"

                output_file = if dimensional
                    case * "-" * eos_string * "-d.csv"
                else
                    case * "-" * eos_string * ".csv"
                end
                paper_runs(folder::AbstractString, eos::Symbol; num_runs = num_runs,
                    output_folder = output_folder,
                    output_file = output_file,
                    case_name = case, dimensional = dimensional)
            end
        end
    end
end

function run_kekatos_comparison_study()
    eos_list = [:ideal, :full_cnga]
    cases = ["GasLib-11", "GasLib-24", "GasLib-40", "GasLib-134"]

    output_folder = "./output/comparison-study/"
    data_folder = "./data/"

    num_runs = 500

    for case in cases
        for eos in eos_list
            for dimensional in [false]
                folder = data_folder * case * "/"
                eos_string = (eos == :ideal) ? "ideal" : "cnga"

                output_file = if dimensional
                    case * "-" * eos_string * "-d.csv"
                else
                    case * "-" * eos_string * ".csv"
                end
                paper_runs(folder::AbstractString, eos::Symbol; num_runs = num_runs,
                    output_folder = output_folder,
                    output_file = output_file,
                    case_name = case, dimensional = dimensional)
            end
        end
    end
end

run_deviation_study() = ideal_vs_non_ideal_runs()

max_stars_repo_name: kaarthiksundar/GasFlowRuns

max_issues_repo_name: kaarthiksundar/GasFlowRuns

max_forks_repo_name: kaarthiksundar/GasFlowRuns


text:
module PDMesh

# Write your package code here.

include("./shapes/shape.jl")
include("./operations/ops.jl")
include("./io.jl")

end

max_stars_repo_name: sajidmannan/PDMesh

max_issues_repo_name: sajidmannan/PDMesh

max_forks_repo_name: sajidmannan/PDMesh


text:
## Exercise 9-2
# Write a function called hasno_e that returns true if the given word doesn’t have the letter e in it.
using Printf

println("Ans: [Part 1]")

function hasno_e(string)
    return !('e' in string)
end

println(hasno_e("Hello"))
println(hasno_e("print"))

# Modify your program from the previous section to print only the words that have no e and compute the percentage of the words in the list that have no e.
println("[Part 2]")

total = 0
have_no_e = 0

for line in eachline("words.txt")
    global total += 1
    if hasno_e(line)
        global have_no_e += 1
        println(line)
    end
end

@printf("%0.2f%% words have no e.\n", ((have_no_e)/total)*100)

println("End.")


text:
ALL_TASK_OFFSETS[("windows-x86", v"1.6.1")] = Dict(
    :END => 132,
    :_isexception => 30,
    :_state => 28,
    :bufsz => 120,
    :copy_stack => 124,
    :ctx => 44,
    :donenotify => 12,
    :eh => 40,
    :excstack => 36,
    :gcstack => 128,
    :logstate => 20,
    :next => 0,
    :prio => 34,
    :queue => 4,
    :result => 16,
    :sizeof_ctx => 72,
    :start => 24,
    :sticky => 29,
    :stkbuf => 116,
    :tid => 32,
    :tls => 8,
    :tls_base_context => 3360,
)


text:
""" 2-category of finitely presented categories.

This module is for the 2-category **Cat** what the module [`FinSets](@ref) is
for the category **Set**: a finitary, combinatorial setting where explicit
calculations can be carried out. We emphasize that the prefix `Fin` means
"finitely presented," not "finite," as finite categories are too restrictive a
notion for many purposes. For example, the free category on a graph is finite if
and only if the graph is DAG, which is a fairly special condition. This usage of
`Fin` is also consistent with `FinSet` because for sets, being finite and being
finitely presented are equivalent.
"""
module FinCats
export FinCat, Path, ob_generators, hom_generators, equations, is_free,
  graph, edges, src, tgt, presentation,
  FinFunctor, FinDomFunctor, is_functorial, collect_ob, collect_hom,
  FinTransformation, components, is_natural

using AutoHashEquals
using Reexport
using StaticArrays: SVector

@reexport using ..Categories
using ...GAT, ...Present, ...Syntax
import ...Present: equations
using ...Theories: Category, ObExpr, HomExpr
import ...Theories: dom, codom, id, compose, ⋅, ∘
using ...CSetDataStructures, ...Graphs, ..FreeDiagrams, ..FinSets
import ...Graphs: edges, src, tgt
import ..FreeDiagrams: FreeDiagram, diagram_type, cone_objects, cocone_objects
import ..Limits: limit, colimit
import ..Categories: Ob, ob, hom, ob_map, hom_map, component

# Categories
############

""" Abstract type for finitely presented category.
"""
abstract type FinCat{Ob,Hom} <: Cat{Ob,Hom} end

FinCat(g::HasGraph, args...; kw...) = FinCatGraph(g, args...; kw...)
FinCat(pres::Presentation, args...; kw...) =
  FinCatPresentation(pres, args...; kw...)

""" Object generators of finitely presented category.

The object generators are almost always the same as the objects. In principle,
however, it is possible to have equations between objects, so that there are
fewer objects than object generators.
"""
function ob_generators end

""" Morphism generators of finitely presented category.
"""
function hom_generators end

""" Is the category freely generated?
"""
is_free(C::FinCat) = isempty(equations(C))

Ob(C::FinCat{Int}) = FinSet(length(ob_generators(C)))

# Discrete categories
#####################

""" Discrete category on a finite set.

The only morphisms in a discrete category are the identities, which in this
implementation are identified with the objects.
"""
@auto_hash_equals struct DiscreteCat{Ob,S<:FinSet{<:Any,Ob}} <: FinCat{Ob,Ob}
  set::S
end
DiscreteCat(n::Integer) = DiscreteCat(FinSet(n))

FinCat(s::Union{FinSet,Integer}) = DiscreteCat(s)

ob_generators(C::DiscreteCat) = C.set
hom_generators(C::DiscreteCat) = ()

dom(C::DiscreteCat{T}, f) where T = f::T
codom(C::DiscreteCat{T}, f) where T = f::T
id(C::DiscreteCat{T}, x) where T = x::T
compose(C::DiscreteCat{T}, f, g) where T = (f::T == g::T) ? f :
  error("Nontrivial composite in discrete category: $f != $g")

# Categories on graphs
######################

""" Abstract type for category backed by finite generating graph.
"""
abstract type FinCatGraph{Ob,Hom} <: FinCat{Ob,Hom} end

""" Generating graph for a finitely presented category.
"""
graph(C::FinCatGraph) = C.graph

ob_generators(C::FinCatGraph) = vertices(graph(C))
hom_generators(C::FinCatGraph) = edges(graph(C))

# Paths in graphs
#----------------

""" Path in a graph.

The path is allowed to be empty but always has definite start and end points
(source and target vertices).
"""
@auto_hash_equals struct Path{V,E,Edges<:AbstractVector{E}}
  edges::Edges
  src::V
  tgt::V
end
edges(path::Path) = path.edges
src(path::Path) = path.src
tgt(path::Path) = path.tgt

function Path(g::HasGraph, es::AbstractVector)
  !isempty(es) || error("Nonempty edge list needed for nontrivial path")
  all(e -> has_edge(g, e), es) || error("Path has edges not contained in graph")
  Path(es, src(g, first(es)), tgt(g, last(es)))
end

function Path(g::HasGraph, e)
  has_edge(g, e) || error("Edge $e not contained in graph $g")
  Path(SVector(e), src(g,e), tgt(g,e))
end

function Base.empty(::Type{Path}, g::HasGraph, v::T) where T
  has_vertex(g, v) || error("Vertex $v not contained in graph $g")
  Path(SVector{0,T}(), v, v)
end

function Base.vcat(p1::Path, p2::Path)
  tgt(p1) == src(p2) ||
    error("Path start/end points do not match: $(tgt(p1)) != $(src(p2))")
  Path(vcat(edges(p1), edges(p2)), src(p1), tgt(p2))
end

""" Abstract type for category whose morphisms are paths in a graph.

(Or equivalence classes of paths in a graph, but we compute with
"""
const FinCatPathGraph{V,E} = FinCatGraph{V,Path{V,E}}

dom(C::FinCatPathGraph, e) = src(graph(C), e)
dom(C::FinCatPathGraph, path::Path) = src(path)
codom(C::FinCatPathGraph, e) = tgt(graph(C), e)
codom(C::FinCatPathGraph, path::Path) = tgt(path)

id(C::FinCatPathGraph, x) = empty(Path, graph(C), x)
compose(C::FinCatPathGraph, fs...) =
  reduce(vcat, coerce_path(graph(C), f) for f in fs)

ob(C::FinCatPathGraph, x) = has_vertex(graph(C), x) ? x :
  error("Vertex $x not contained in graph $(graph(C))")
hom(C::FinCatPathGraph, f) = coerce_path(graph(C), f)

coerce_path(g::HasGraph, path::Path) = path
coerce_path(g::HasGraph, x) = Path(g, x)

# Free category on graph
#-----------------------

""" Free category generated by a finite graph.

The objects of the free category are vertices in the graph and the morphisms are
(possibly empty) paths.
"""
@auto_hash_equals struct FreeCatGraph{G<:HasGraph} <: FinCatPathGraph{Int,Int}
  graph::G
end

FinCatGraph(g::HasGraph) = FreeCatGraph(g)

is_free(::FreeCatGraph) = true

# Category on graph with equations
#---------------------------------

""" Category presented by a finite graph together with path equations.

The objects of the category are vertices in the graph and the morphisms are
paths, quotiented by the congruence relation generated by the path equations.
See (Spivak, 2014, *Category theory for the sciences*, §4.5).
"""
@auto_hash_equals struct FinCatGraphEq{G<:HasGraph,Eqs<:AbstractVector{<:Pair}} <:
    FinCatPathGraph{Int,Int}
  graph::G
  equations::Eqs
end

equations(C::FinCatGraphEq) = C.equations

function FinCatGraph(g::HasGraph, eqs::AbstractVector)
  eqs = map(eqs) do eq
    lhs, rhs = coerce_path(g, first(eq)), coerce_path(g, last(eq))
    (src(lhs) == src(rhs) && tgt(lhs) == tgt(rhs)) ||
      error("Paths $lhs and $rhs in equation do not have equal (co)domains")
    lhs => rhs
  end
  FinCatGraphEq(g, eqs)
end

# Symbolic categories
#####################

""" Category defined by a `Presentation` object.
"""
struct FinCatPresentation{Ob,Hom} <: FinCat{Ob,Hom}
  presentation::Presentation

  function FinCatPresentation(pres::Presentation)
    new{pres.syntax.Ob,pres.syntax.Hom}(pres)
  end
end

presentation(C::FinCatPresentation) = C.presentation
ob_generators(C::FinCatPresentation) = generators(presentation(C), :Ob)
hom_generators(C::FinCatPresentation) = generators(presentation(C), :Hom)
equations(C::FinCatPresentation) = equations(presentation(C))

ob(C::FinCatPresentation, x) = ob(C, presentation(C)[x])
ob(C::FinCatPresentation, x::GATExpr) = x
  # FIXME: Commented for now because `DataMigration` uses for attr types.
  # gat_typeof(x) == :Ob ? x : error("Expression $x is not an object")

hom(C::FinCatPresentation, f) = hom(C, presentation(C)[f])
hom(C::FinCatPresentation, f::GATExpr) = f
  # FIXME: Commented for now because `DataMigration` uses for attributes.
  # gat_typeof(f) == :Hom ? f : error("Expression $f is not a morphism")
hom(C::FinCatPresentation, fs::AbstractVector) =
  mapreduce(f -> hom(C, f), compose, fs)

# Functors
##########

""" A functor out of a finitely presented category.
"""
const FinDomFunctor{Dom<:FinCat,Codom<:Cat} = Functor{Dom,Codom}

FinDomFunctor(maps::NamedTuple{(:V,:E)}, dom::FinCatGraph, codom::Cat) =
  FinDomFunctor(maps.V, maps.E, dom, codom)

# Diagram interface. See `FreeDiagrams` module.
diagram_type(F::FinDomFunctor{Dom,Codom}) where {Ob,Hom,Dom,Codom<:Cat{Ob,Hom}} =
  Tuple{Ob,Hom}
cone_objects(F::FinDomFunctor) = collect_ob(F)
cocone_objects(F::FinDomFunctor) = collect_ob(F)

function hom_map(F::FinDomFunctor{<:FinCatPathGraph}, path::Path)
  D = codom(F)
  mapreduce(e -> hom_map(F, e), (gs...) -> compose(D, gs...),
            edges(path), init=id(D, ob_map(F, src(path))))
end

ob_map(F::FinDomFunctor, x::GATExpr{:generator}) = ob_map(F, first(x))
hom_map(F::FinDomFunctor, f::GATExpr{:generator}) = hom_map(F, first(f))
hom_map(F::FinDomFunctor, f::GATExpr{:id}) = id(codom(F), ob_map(F, dom(f)))

function hom_map(F::FinDomFunctor, f::GATExpr{:compose})
  D = codom(F)
  mapreduce(f -> hom_map(F, f), (gs...) -> compose(D, gs...), args(f))
end

(F::FinDomFunctor)(expr::ObExpr) = ob_map(F, expr)
(F::FinDomFunctor)(expr::HomExpr) = hom_map(F, expr)

""" Is the purported functor on a presented category functorial?

This function checks that functor preserves domains and codomains. When
`check_equations` is `true` (the default), it also checks that the functor
preserves all equations in the domain category. Note that in some cases this may
not be possible.

See also: [`is_natural`](@ref).
"""
function is_functorial(F::FinDomFunctor; check_equations::Bool=false)
  C, D = dom(F), codom(F)
  all(hom_generators(C)) do f
    g = hom_map(F, f)
    dom(D,g) == ob_map(F, dom(C,f)) && codom(D,g) == ob_map(F, codom(C,f))
  end || return false

  if check_equations
    all(equations(C)) do (lhs, rhs)
      is_hom_equal(D, hom_map(F, lhs), hom_map(F, rhs))
    end || return false
  end

  true
end

""" A functor between finitely presented categories.
"""
const FinFunctor{Dom<:FinCat,Codom<:FinCat} = FinDomFunctor{Dom,Codom}

FinFunctor(maps, dom::FinCat, codom::FinCat) = FinDomFunctor(maps, dom, codom)
FinFunctor(ob_map, hom_map, dom::FinCat, codom::FinCat) =
  FinDomFunctor(ob_map, hom_map, dom, codom)
FinFunctor(ob_map, hom_map, dom::Presentation, codom::Presentation) =
  FinDomFunctor(ob_map, hom_map, FinCat(dom), FinCat(codom))

# Mapping-based functors
#-----------------------

""" Functor out of a finitely presented category given by explicit mappings.
"""
@auto_hash_equals struct FinDomFunctorMap{Dom<:FinCat, Codom<:Cat,
    ObMap, HomMap} <: FinDomFunctor{Dom,Codom}
  ob_map::ObMap
  hom_map::HomMap
  dom::Dom
  codom::Codom
end

function FinDomFunctor(ob_map::Union{AbstractVector,AbstractDict},
                       hom_map::Union{AbstractVector,AbstractDict},
                       dom::FinCat, codom::Cat)
  # FIXME: Disable checks for dicts due to `DataMigration` issue above.
  ob_map isa AbstractDict || length(ob_map) == length(ob_generators(dom)) ||
    error("Length of object map $ob_map does not match domain $dom")
  hom_map isa AbstractDict || length(hom_map) == length(hom_generators(dom)) ||
    error("Length of morphism map $hom_map does not match domain $dom")

  ob_map = mappairs(x -> functor_key(dom, x), y -> ob(codom, y), ob_map)
  hom_map = mappairs(f -> functor_key(dom, f), g -> hom(codom, g), hom_map)
  FinDomFunctorMap(ob_map, hom_map, dom, codom)
end
FinDomFunctor(ob_map::Union{AbstractVector{Ob},AbstractDict{<:Any,Ob}},
              hom_map::Union{AbstractVector{Hom},AbstractDict{<:Any,Hom}},
              dom::FinCat) where {Ob,Hom} =
  FinDomFunctor(ob_map, hom_map, dom, TypeCat(Ob, Hom))

functor_key(C::FinCat, x) = x
functor_key(C::FinCat, expr::GATExpr) = head(expr) == :generator ?
  first(expr) : error("Functor must be defined on generators")

function Categories.do_compose(F::FinDomFunctorMap, G::FinDomFunctorMap)
  FinDomFunctorMap(mapvals(x -> ob_map(G, x), F.ob_map),
                   mapvals(f -> hom_map(G, f), F.hom_map), dom(F), codom(G))
end

""" Functor object and morphism maps given as vectors.
"""
const FinDomFunctorVector{Dom<:FinCat,Codom<:Cat,
                          ObMap<:AbstractVector,HomMap<:AbstractVector} =
  FinDomFunctorMap{Dom,Codom,ObMap,HomMap}

ob_map(F::FinDomFunctorVector, x::Integer) = F.ob_map[x]
hom_map(F::FinDomFunctorVector, f::Integer) = F.hom_map[f]

collect_ob(F::FinDomFunctorVector) = F.ob_map
collect_hom(F::FinDomFunctorVector) = F.hom_map

Ob(F::FinDomFunctorVector) = FinDomFunction(F.ob_map, Ob(codom(F)))

""" Functor with object and morphism maps given as dictionaries.
"""
const FinDomFunctorDict{Dom<:FinCat,Codom<:Cat,ObKey,HomKey,
                        ObMap<:AbstractDict{ObKey},HomMap<:AbstractDict{HomKey}} =
  FinDomFunctorMap{Dom,Codom,ObMap,HomMap}

ob_map(F::FinDomFunctorDict{Dom,Codom,ObKey}, x::ObKey) where
  {Dom,Codom,ObKey} = F.ob_map[x]
hom_map(F::FinDomFunctorDict{Dom,Codom,ObKey,HomKey}, f::HomKey) where
  {Dom,Codom,ObKey,HomKey} = F.hom_map[f]

# C-set interop
#--------------

function FinDomFunctor(pres::Presentation, X::ACSet)
  ob_map = Dict(c => FinSet(X, nameof(c)) for c in generators(pres, :Ob))
  hom_map = Dict(f => FinFunction(X, nameof(f)) for f in generators(pres, :Hom))
  FinDomFunctor(ob_map, hom_map,
                FinCat(pres), TypeCat(FinSet{Int}, FinFunction{Int}))
end

# Free diagram interop
#---------------------

function FreeDiagram(F::FinDomFunctor{<:FreeCatGraph,<:TypeCat{Ob,Hom}}) where {Ob,Hom}
  diagram = FreeDiagram{Ob,Hom}()
  copy_parts!(diagram, graph(dom(F)))
  diagram[:ob] = collect_ob(F)
  diagram[:hom] = collect_hom(F)
  diagram
end

function FinDomFunctor(diagram::FreeDiagram)
  g = Graph()
  copy_parts!(g, diagram)
  FinDomFunctor(ob(diagram), hom(diagram), FinCat(g))
end

limit(F::FinDomFunctor) = limit(FreeDiagram(F))
colimit(F::FinDomFunctor) = colimit(FreeDiagram(F))

# Natural transformations
#########################

""" A natural transformation whose domain category is finitely generated.

This type is for natural transformations ``α: F ⇒ G: C → D`` such that the
domain category ``C`` is finitely generated. Such a natural transformation is
given by a finite amount of data (one morphism in ``D`` for each generating
object of ``C``) and its naturality is verified by finitely many equations (one
equation for each generating morphism of ``C``).
"""
const FinTransformation{C<:FinCat,D<:Cat,Dom<:FinDomFunctor{C,D},Codom<:FinDomFunctor{C,D}} =
  Transformation{C,D,Dom,Codom}

FinTransformation(F, G; components...) = FinTransformation(components, F, G)

""" Components of a natural transformation.
"""
components(α::FinTransformation) = α.components

""" Is the transformation between `FinDomFunctors` a natural transformation?

This function uses the fact that to check whether a transformation is natural,
it suffices to check the naturality equations on a generating set of morphisms
of the domain category. In some cases, checking the equations may be expensive
or impossible. When the keyword argument `check_equations` is `false`, only the
domains and codomains of the components are checked.

See also: [`is_functorial`](@ref).
"""
function is_natural(α::FinTransformation; check_equations::Bool=true)
  F, G = dom(α), codom(α)
  C, D = dom(F), codom(F) # == dom(G), codom(G)
  all(ob_generators(C)) do c
    α_c = α[c]
    dom(D, α_c) == ob_map(F,c) && codom(D, α_c) == ob_map(G,c)
  end || return false

  if check_equations
    all(hom_generators(C)) do f
      Ff, Gf = hom_map(F,f), hom_map(G,f)
      α_c, α_d = α[dom(C,f)], α[codom(C,f)]
      is_hom_equal(D, compose(D, α_c, Gf), compose(D, Ff, α_d))
    end || return false
  end

  true
end

function check_transformation_domains(F::Functor, G::Functor)
  (C = dom(F)) == dom(G) ||
    error("Mismatched domains in functors $F and $G")
  (D = codom(F)) == codom(G) ||
    error("Mismatched codomains in functors $F and $G")
  (C, D)
end

# Mapping-based transformations
#------------------------------

""" Natural transformation with components given by explicit mapping.
"""
@auto_hash_equals struct FinTransformationMap{C<:FinCat,D<:Cat,
    Dom<:FinDomFunctor{C,D},Codom<:FinDomFunctor{C,D},Comp} <:
    FinTransformation{C,D,Dom,Codom}
  components::Comp
  dom::Dom
  codom::Codom
end

function FinTransformation(components::Union{AbstractVector,AbstractDict},
                           F::FinDomFunctor, G::FinDomFunctor)
  C, D = check_transformation_domains(F, G)
  length(components) == length(ob_generators(C)) ||
    error("Incorrect number of components in $components for domain category $C")
  components = mappairs(x -> transformation_key(C,x), f -> hom(D,f), components)
  FinTransformationMap(components, F, G)
end

transformation_key(C::FinCat, x) = x
transformation_key(C::FinCat, expr::GATExpr) = head(expr) == :generator ?
  first(expr) : error("Natural transformation must be defined on generators")

component(α::FinTransformationMap{C,D,F,G,Comp}, c::Integer) where
  {C,D,F,G,Comp<:AbstractVector} = α.components[c]
component(α::FinTransformationMap{C,D,F,G,Comp}, c::Key) where
  {Key,C,D,F,G,Comp<:AbstractDict{Key}} = α.components[c]
component(α::FinTransformationMap, expr::GATExpr) =
  component(α, first(expr))

function Categories.do_compose(α::FinTransformationMap, β::FinTransformation)
  F = dom(α)
  D = codom(F)
  FinTransformationMap(mapvals(α.components, keys=true) do c, f
                         compose(D, f, component(β, c))
                       end, F, codom(β))
end

function Categories.do_composeH(F::FinDomFunctorMap, β::Transformation)
  G, H = dom(β), codom(β)
  FinTransformationMap(mapvals(c -> component(β, c), F.ob_map),
                       compose(F, G), compose(F, H))
end

function Categories.do_composeH(α::FinTransformationMap, H::Functor)
  F, G = dom(α), codom(α)
  FinTransformationMap(mapvals(f -> hom_map(H, f), α.components),
                       compose(F, H), compose(G, H))
end

# Dict utilities
################

# Something like this should be built into Julia, but unfortunately isn't.

function mappairs(kmap, vmap, pairs::T) where T<:AbstractDict
  (dicttype(T))(kmap(k) => vmap(v) for (k,v) in pairs)
end
mappairs(kmap, vmap, vec::AbstractVector) = map(vmap, vec)

function mapvals(f, pairs::T; keys::Bool=false) where T<:AbstractDict
  (dicttype(T))(if keys
    (k => f(k,v) for (k,v) in pairs)
  else
    (k => f(v) for (k,v) in pairs)
  end)
end
mapvals(f, vec::AbstractVector; keys::Bool=false) =
  keys ? map(f, eachindex(vec), vec) : map(f, vec)

dicttype(::Type{T}) where T <: AbstractDict = T.name.wrapper
dicttype(::Type{<:Iterators.Pairs}) = Dict

end


text:
include("../src/lstm_g/network.jl")
include("../src/lstm_g/viz.jl")

nin = 3
nhidden = 2
nout = 1

# basic ANN
inputlayer  = gatedlayer(3, tag=:input)
hiddenlayer = gatedlayer(2, tag=:hidden)
outputlayer = gatedlayer(1, tag=:output)
@show inputlayer hiddenlayer outputlayer

ci = connect(inputlayer, hiddenlayer)
co = connect(hiddenlayer, outputlayer)
@show ci co

# now add a gating gatedlayer
gatinglayer = gatedlayer(2, tag=:gater)
cg = connect(inputlayer, gatinglayer)

# gate ci
gate!(ci, gatinglayer)

# peephole uses SAME type... node i maps to node i
peep = connect(hiddenlayer, gatinglayer, SAME)

@show gatinglayer

net = GatedNetwork(inputlayer, hiddenlayer, GatedLayer{Float64}[inputlayer, gatinglayer, hiddenlayer, outputlayer])


text:
module TestExamples

import Modia3D

# Test Modia3D - Examples

collisionPath = joinpath(Modia3D.path, "examples", "collisions")
include(joinpath(collisionPath, "Simulate_Billiards_OneBall.jl"))
include(joinpath(collisionPath, "Simulate_BouncingBall.jl"))
include(joinpath(collisionPath, "Simulate_NewtonsCradle.jl"))
include(joinpath(collisionPath, "Simulate_SlidingAndRollingBall.jl"))
include(joinpath(collisionPath, "Simulate_TwoCollidingBalls.jl"))
include(joinpath(collisionPath, "Simulate_YouBot.jl"))

# test/collision
collisionPath = joinpath(Modia3D.path, "test", "collision")
include(joinpath(collisionPath, "Collision_3Elements.jl"))
# include(joinpath(collisionPath, "Collision_Bars.jl")) # is not working needs to be fixed
# include(joinpath(collisionPath, "Plot_cor.jl"))             # not included since PyPlot required
#include(joinpath(collisionPath, "Plot_SlidingFriction.jl"))  # not included since PyPlot required
include(joinpath(collisionPath, "Test_Collision.jl"))
include(joinpath(collisionPath, "Test_Collision_moreRevolutes.jl"))
include(joinpath(collisionPath, "Test_Collision_StarSetting.jl"))
include(joinpath(collisionPath, "Test_MiniBsp.jl"))
include(joinpath(collisionPath, "Test_Solids.jl"))

# test/collision/boxOnTable
collisionPath = joinpath(Modia3D.path, "test", "collision","boxOnTable")
include(joinpath(collisionPath, "Simulate_ContactBoxOnTable.jl"))
include(joinpath(collisionPath, "Simulate_YouBotBoxOnTable.jl"))

# test/collision/dynamic_collision
collisionPath = joinpath(Modia3D.path, "test", "collision","dynamic_collision")
include(joinpath(collisionPath, "collision_2_boxes.jl"))
include(joinpath(collisionPath, "collision_ballWithBall.jl"))
include(joinpath(collisionPath, "collision_ballWithBox.jl"))
include(joinpath(collisionPath, "collision_ballWithBox_45Deg.jl"))
include(joinpath(collisionPath, "collision_BallWithBox_Prismatic.jl"))
include(joinpath(collisionPath, "collision_ballWithBox_45Deg.jl"))
include(joinpath(collisionPath, "collision_newtons_cradle.jl"))

# test/collision/contactForceLaw
collisionPath = joinpath(Modia3D.path, "test", "collision","contactForceLaw")
# include(joinpath(collisionPath, "bouncingBallComparedWithImpulses.jl"))  # not included since PyPlot required
include(joinpath(collisionPath, "contactForceLaw_2_boxes.jl"))
include(joinpath(collisionPath, "contactForceLaw_2_boxes_Prismatic.jl"))
include(joinpath(collisionPath, "contactForceLaw_2_boxes2.jl"))
include(joinpath(collisionPath, "contactForceLaw_Ball.jl"))
include(joinpath(collisionPath, "contactForceLaw_ballWithBall.jl"))
include(joinpath(collisionPath, "contactForceLaw_ballWithBox_45Deg.jl"))
include(joinpath(collisionPath, "contactForceLaw_ballWithBox_Prismatic.jl"))
include(joinpath(collisionPath, "contactForceLaw_newtons_cradle.jl"))

# test/collision/Billard_Simulation
collisionPath = joinpath(Modia3D.path, "test", "collision","Billard_Simulations", "BallCushion")
include(joinpath(collisionPath, "BillardBall1_Cushion1_directHit.jl"))
include(joinpath(collisionPath, "BillardBall1_Cushion4_arbitraryHit.jl"))
# include(joinpath(collisionPath, "BillardBall1_Cushion4_directHit.jl"))  # Successful on Windows and Linux, buf fails on Mac with
                                                                          # mxstep steps taken before reaching tout.
                                                                          # Try to set maxNumberOfSteps to a value > 500

# collisionPath = joinpath(Modia3D.path, "test", "collision","Billard_Simulations", "rollingBall")
# include(joinpath(collisionPath, "colliding2rollingBalls_plot1.jl"))   # not included since PyPlot required
# include(joinpath(collisionPath, "colliding2rollingBalls_plot2.jl"))   # not included since PyPlot required
# include(joinpath(collisionPath, "rollingBall1.jl"))                   # not included since PyPlot required

collisionPath = joinpath(Modia3D.path, "test", "collision","Billard_Simulations", "billard")
# include(joinpath(collisionPath, "Billard_Balls4.jl")) # takes not too long, but still to long for an example test
# include(joinpath(collisionPath, "Billard_Balls4_Prismatic.jl")) # takes not too long, but still to long for an example test
# include(joinpath(collisionPath, "Billard_Balls16.jl")) # takes too long

println("\n ...test_Examples_Collision finished!")
end


text:
using Krylov, Stopping

include("random_coordinate_descent_method.jl")
include("stop_random_coordinate_descent_method.jl")
include("stop_random_coordinate_descent_method_LS.jl")

n=5000;
A=rand(n,n)
sol=rand(n)
b=A*sol

#@time stp = StopRandomizedCD(A,b, max_iter = 1000)

x0   = zeros(size(A,2))
pb   = LinearSystem(A,b)
s0   = GenericState(x0, similar(b))
mcnt = Stopping._init_max_counters_linear_operators()
#
@time meta = StoppingMeta(max_cntrs = mcnt,
                        atol = 1e-7, rtol = 1e-15, max_iter = 99,
                        retol = false,
                        tol_check = (atol, rtol, opt0)->(atol + rtol * opt0),
                        optimality_check = (pb, state) -> state.res)
s0   = GenericState(x0, similar(b))
@time stp = LAStopping(pb, meta, cheap_stop_remote_control(), s0)
@time StopRandomizedCD(stp)

@time x, OK = RandomizedCD(A,b, max_iter = 100)

###############################################################################
#
# We compare the stopping_randomized_CD with the same version with an
# artificial substopping created and called at each iteration.
# It appears that the expansive part is to create the SubStopping, and in
# particular create the StoppingMeta.
#
@time meta = StoppingMeta(max_cntrs = mcnt,
                        atol = 1e-7, rtol = 1e-15, max_iter = 99,
                        retol = false,
                        tol_check = (atol, rtol, opt0)->(atol + rtol * opt0),
                        optimality_check = (pb, state) -> state.res)
s0   = GenericState(x0, similar(b))
@time stp2 = LAStopping(pb, meta, cheap_stop_remote_control(), s0)
@time StopRandomizedCD_LS(stp2)

;


text:
using SynchrotronKernel, Test

@testset "Synchrotron Kernel" begin

    @testset "F" begin
        @test SynchrotronKernel.F(1.e-8) ≈ 0.0046310000727714734
        @test SynchrotronKernel.F(  1.0) ≈ 0.6514228153553652
        @test SynchrotronKernel.F(  3.5) ≈ 0.08268719536694746
        @test SynchrotronKernel.F(  4.5) ≈ 0.03357441971502366
        @test SynchrotronKernel.F( 10.0) ≈ 0.00019223826430086885
        @test SynchrotronKernel.F( 50.0) ≈ 1.734785203976593e-21
        @test SynchrotronKernel.F(100.0) ≈ 4.697593665922202e-43
        @test SynchrotronKernel.F( 1.e6) == 0.0

        @test ℱ(1.0) == SynchrotronKernel.F(1.0)
    end

    @testset "G" begin
        @test SynchrotronKernel.G(1.0) ≈ 0.4944750621042115

        @test 𝒢(1.0) == SynchrotronKernel.G(1.0)
    end

    @testset "Synchrotron Kernel" begin
        K_tot, K_ort, K_par = synchrotron_kernel(1.0)

        @test K_tot == SynchrotronKernel.F(1.0)
        @test K_ort ≈ 0.5729489387297884
        @test K_par ≈ 0.07847387662557684
        @test K_tot == K_ort + K_par
    end

    @testset "Intensity" begin
        @test synchrotron_intensity(1.0) == SynchrotronKernel.F(1.0)
    end

    @testset "Polarisation" begin
        K = SynchrotronKernel.F(1.0)
        sk_ort, sk_par = synchrotron_polarisation(1.0)

        @test sk_ort ≈ 0.5729489387297884
        @test sk_par ≈ 0.07847387662557684
        @test K == sk_ort + sk_par
    end

end


text:
function extract_emails(str::String)
    # this might be better
    # eachmatch_string(regex, str) = (String(i.match) for i in eachmatch(regex, str))
    return unique(
        String.(
            SubString.(
                str, findall(r"([a-zA-Z0-9_.+-]+@[a-zA-Z0-9-]+\.[a-zA-Z0-9-.]+)", str)
            ),
        ),
    )
end

function valid_email(str::String)
    ending = String(last(split(str, '.')))
    if str == ""
        return false
    elseif length(split(str, '.')) < 2
        return false
    elseif !(ending in TLDs)
        return false
    else
        return true
    end
end

function find_emails(body::HTMLElement)::Vector{String}
    emails = Vector{Union{String,Missing}}()
    elements_p = eachmatch(
        Selector("p:matches([a-zA-Z0-9._%+-]+@[a-zA-Z0-9.-]+\\.[a-zA-Z]{2,6})"), body
    )
    elements_a = eachmatch(Selector("a[href^=\"mailto:\"]"), body)
    if length(elements_p) > 0
        emails_p = mapreduce(elem -> extract_emails(nodeText(elem)), vcat, elements_p)
        emails = vcat(emails, emails_p)
        filter!(x -> x isa String, emails)
    end
    if length(elements_a) > 0
        emails_a = mapreduce(
            elem -> extract_emails(getattr(elem, "href")), vcat, elements_a
        )
        emails = vcat(emails, emails_a)
        filter!(x -> x isa String, emails)
    end
    return unique(skipmissing(emails))
end


text:
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/build_literate/Greens_functions.ipynb)

# # Green's Functions
# The Green's functions are given by the Fourier transform of the corresponding lineshapes. For a Lorentzian lineshape, the Green's function is
# ```math
# 	G(t,\tau) = \exp (-R_2^s (t-\tau)) \;\; \forall \;\; t \geq \tau,
# ```
# for a Gaussian lineshape it is
# ```math
# 	G(t,\tau) = \exp(- {R_2^s}^2 (t-\tau)^2 / 2)),
# ```
# and for [super-Lorentzian lineshape](http://dx.doi.org/10.1002/mrm.1910330404) it is
# ```math
# 	G(t,\tau) = \int_0^{1} \exp \left(- {R_2^s}^2 (t - \tau)^2 \cdot  \frac{(3 \zeta^2 - 1)^2}{8} \right) d\zeta.
# ```
# As evident from these equations, the Green's functions are merely a function of ``\kappa = R_2^s \cdot (t - \tau) = (t - \tau) / T_2^s``, and in this package we implemented the functions as such: [`greens_lorentzian(κ)`](@ref), [`greens_gaussian(κ)`](@ref), and [`greens_superlorentzian(κ)`](@ref). These functions can be used to reproduce Fig. 1 in the generalized Bloch paper:

using MRIgeneralizedBloch
using Plots
plotlyjs(bg = RGBA(31/255,36/255,36/255,1.0), ticks=:native); #hide #!nb
 #-
T₂ˢ = 10e-6 # s
t = 0 : 1e-6 : 1e-3
p = plot(yaxis=:log, ylim=(1e-6,1), xlabel="(t-τ) [ms]", ylabel="G((t-τ)/T₂ˢ)")
plot!(p, 1e3t, greens_lorentzian.(t ./ T₂ˢ), label="Lorentzian lineshape")
plot!(p, 1e3t, greens_gaussian.(t ./ T₂ˢ), label="Gaussian lineshape")
plot!(p, 1e3t, greens_superlorentzian.(t ./ T₂ˢ), label="super-Lorentzian l.")
#md Main.HTMLPlot(p) #hide

#src ###############################################################
#src export data for plotting
#src ###############################################################
using Printf #src
io = open(expanduser(string("~/Documents/Paper/2021_MT_IDE/Figures/Greens_Functions.txt")), "w") #src
write(io, "t_ms tau G_Lorentzian G_Gaussian G_superLorentzian \n") #src
for i = 1:length(t) #src
    write(io, @sprintf("%1.3e %1.3e %1.3e %1.3e %1.3e \n",  #src
    t[i] * 1e3, #src
    t[i] / T₂ˢ, #src
    greens_lorentzian(t[i] / T₂ˢ), #src
    greens_gaussian(t[i] / T₂ˢ), #src
    greens_superlorentzian(t[i] / T₂ˢ) #src
    )) #src
end #src
close(io) #src


text:
export addprocs_scyld, ScyldManager

immutable ScyldManager <: ClusterManager
end

function launch(manager::ScyldManager, np::Integer, config::Dict, instances_arr::Array, c::Condition)
    try
        home = config[:dir]
        exename = config[:exename]
        exeflags = config[:exeflags]

        beomap_cmd = `bpsh -1 beomap --no-local --np $np`
        out,beomap_proc = open(beomap_cmd)
        wait(beomap_proc)
        if !success(beomap_proc)
            error("node availability inaccessible (could not run beomap)")
        end
        nodes = split(chomp(readline(out)),':')
        io_objs = cell(np)
        configs = cell(np)
        for (i,node) in enumerate(nodes)
            cmd = `bpsh $node sh -l -c "cd $home && $(exename) $(exeflags)"`
            cmd.detach = true
            configs[i] = merge(config, {:node => node})
            io_objs[i],_ = open(cmd)
            io_objs[i].line_buffered = true
        end

        push!(instances_arr, collect(zip(io_objs, configs)))
        notify(c)
   catch e
        println("Error launching beomap")
        println(e)
   end
end

function manage(manager::ScyldManager, id::Integer, config::Dict, op::Symbol)
    if op == :interrupt
        if haskey(config, :ospid)
            node = config[:node]
            if !success(`bpsh $node kill -2 $(config[:ospid])`)
                println("Error sending Ctrl-C to julia worker $id on node $node")
            end
        else
            # This state can happen immediately after an addprocs
            println("Worker $id cannot be presently interrupted.")
        end
    elseif op == :register
        config[:ospid] = remotecall_fetch(id, getpid)
    end
end

addprocs_scyld(np::Integer) = addprocs(np, manager=ScyldManager())


text:
#=
Plots of bivariate Gaussians to illustrate the Kalman filter.

@author : Spencer Lyon <spencer.lyon@nyu.edu>
=#
using PyPlot

# Quick meshgrid function
meshgrid(x::Vector, y::Vector) = (repmat(x, 1, length(y))',
                                  repmat(y, 1, length(x)))

# bivariate normal function. I could call plt.mlab, but this is more fun!
# See http://mathworld.wolfram.com/BivariateNormalDistribution.html
function bivariate_normal(X::Matrix, Y::Matrix, σ_x::Real=1.0, σ_y::Real=1.0,
                          μ_x::Real=0.0, μ_y::Real=0.0, σ_xy::Real=0.0)
    Xμ = X .- μ_x
    Yμ = Y .- μ_y

    ρ = σ_xy/(σ_x*σ_y)
    z = Xμ.^2/σ_x^2 + Yμ.^2/σ_y^2 - 2*ρ.*Xμ.*Yμ/(σ_x*σ_y)
    denom = 2π*σ_x*σ_y*sqrt(1-ρ^2)
    return exp( -z/(2*(1-ρ^2))) ./ denom
end

# == Set up the Gaussian prior density p == #
Σ = [0.4 0.3
     0.3 0.45]
x_hat = [0.2
         -0.2]''

# == Define the matrices G and R from the equation y = G x + N(0, R) == #
G = eye(2)
R = 0.5 .* Σ

# == The matrices A and Q == #
A = [1.2 0
     0   -0.2]
Q = 0.3 .* Σ

# == The observed value of y == #
y = [2.3, -1.9]''

# == Set up grid for plotting == #
x_grid = linspace(-1.5, 2.9, 100)
y_grid = linspace(-3.1, 1.7, 100)
X, Y = meshgrid(x_grid, y_grid)

function gen_gaussian_plot_vals(mu, C)
    "Z values for plotting the bivariate Gaussian N(mu, C)"
    m_x, m_y = mu[1], mu[2]
    s_x, s_y = sqrt(C[1, 1]), sqrt(C[2, 2])
    s_xy = C[1, 2]
    return bivariate_normal(X, Y, s_x, s_y, m_x, m_y, s_xy)
end

# helper function to prepare axis
function prep_ax()
    fig, ax = subplots()
    ax[:xaxis][:grid](true, zorder=0)
    ax[:yaxis][:grid](true, zorder=0)
    return ax
end

function plot1()
    ax = prep_ax()
    Z = gen_gaussian_plot_vals(x_hat, Σ)
    ax[:contourf](X, Y, Z, 6, alpha=0.6, cmap=ColorMap("jet"))
    cs = ax[:contour](X, Y, Z, 6, colors="black")
    ax[:clabel](cs, inline=1, fontsize=10)
end

function plot2()
    ax = prep_ax()
    Z = gen_gaussian_plot_vals(x_hat, Σ)
    ax[:contourf](X, Y, Z, 6, alpha=0.6, cmap=ColorMap("jet"))
    cs = ax[:contour](X, Y, Z, 6, colors="black")
    ax[:clabel](cs, inline=1, fontsize=10)
    ax[:text](y[1], y[2], L"$y$", fontsize=20, color="black")
end

function plot3()
    ax = prep_ax()
    Z = gen_gaussian_plot_vals(x_hat, Σ)
    cs1 = ax[:contour](X, Y, Z, 6, colors="black")
    ax[:clabel](cs1, inline=1, fontsize=10)
    M = Σ * G' * inv(G * Σ * G' + R)
    x_hat_F = x_hat + M * (y - G * x_hat)
    Sigma_F = Σ - M * G * Σ
    new_Z = gen_gaussian_plot_vals(x_hat_F, Sigma_F)
    cs2 = ax[:contour](X, Y, new_Z, 6, colors="black")
    ax[:clabel](cs2, inline=1, fontsize=10)
    ax[:contourf](X, Y, new_Z, 6, alpha=0.6, cmap=ColorMap("jet"))
    ax[:text](y[1], y[2], L"$y$", fontsize=20, color="black")
end

function plot4()
    ax = prep_ax()
    # Density 1
    Z = gen_gaussian_plot_vals(x_hat, Σ)
    cs1 = ax[:contour](X, Y, Z, 6, colors="black")
    ax[:clabel](cs1, inline=1, fontsize=10)
    # Density 2
    M = Σ * G' * inv(G * Σ * G' + R)
    x_hat_F = x_hat + M * (y - G * x_hat)
    Sigma_F = Σ - M * G * Σ
    Z_F = gen_gaussian_plot_vals(x_hat_F, Sigma_F)
    cs2 = ax[:contour](X, Y, Z_F, 6, colors="black")
    ax[:clabel](cs2, inline=1, fontsize=10)
    # Density 3
    new_x_hat = A * x_hat_F
    new_Sigma = A * Sigma_F * A' + Q
    new_Z = gen_gaussian_plot_vals(new_x_hat, new_Sigma)
    cs3 = ax[:contour](X, Y, new_Z, 6, colors="black")
    ax[:clabel](cs3, inline=1, fontsize=10)
    ax[:contourf](X, Y, new_Z, 6, alpha=0.6, cmap=ColorMap("jet"))
    ax[:text](y[1], y[2], L"$y$", fontsize=20, color="black")
end