text:
# ____________________________________________________________________
# SpecLineAsymmGauss
#
mutable struct SpecLineAsymmGauss <: AbstractComponent
    norm::Parameter
    center::Parameter
    fwhm::Parameter
    voff::Parameter
    asymm::Parameter
    index::Vector{Int}  # optimization

    function SpecLineAsymmGauss(center::Number)
        out = new(Parameter(1),
                  Parameter(center),
                  Parameter(3000),
                  Parameter(0),
                  Parameter(0),
                  Vector{Float64}(), true)

        @assert center > 0
        out.norm.low = 0
        out.center.low = 0
        out.fwhm.low = 0
        out.voff.low = 0
        out.center.fixed = true
        return out
    end
end

function prepare!(comp::SpecLineAsymmGauss, domain::Domain{1})
    comp.index = collect(1:length(domain))
    return fill(NaN, length(domain))
end

function evaluate!(buffer, comp::SpecLineAsymmGauss, x::Domain{1},
                   norm, center, fwhm, voff, asymm)
    buffer[comp.index] .= 0.
    empty!(comp.index)

    x0 = center - (voff / 3.e5) * center
    hwhm = fwhm / 3.e5 * center / 2  # Note: this is in `center` units

    sigma0 = hwhm / (2.355 / 2)
    sigma = 2. * sigma0 ./ (1 .+ exp.(asymm .* (x .- x0) ./ 2 ./ sigma0))
    X = (x .- x0) ./ sigma
    i = findall(abs.(X) .< 4) # optimization
    append!(comp.index, i)
    buffer[i] .= norm * exp.(-X[i].^2 ./ 2) ./ sqrt(2pi) * sigma0
end


text:
include("metrics/seed.jl")
include("metrics/pricer.jl")
include("metrics/pricer_cv.jl")
include("metrics/distribution.jl")
include("metrics/delta.jl")
include("metrics/rho.jl")
include("metrics/confinterval.jl")
include("metrics/variance.jl")

ext: jl


text:
module Aqua

@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), "```julia" => "```jldoctest")
end Aqua

using Base: PkgId
using Test

include("ambiguities.jl")
include("unbound_args.jl")
include("exports.jl")

"""
    test_all(testtarget::Module)

Run following tests in isolated testset:

* [`test_ambiguities([testtarget, Base])`](@ref test_ambiguities)
* [`test_unbound_args(testtarget)`](@ref test_unbound_args)
* [`test_undefined_exports(testtarget)`](@ref test_undefined_exports)

# Keyword Arguments
- `ambiguities`: Keyword arguments passed to [`test_ambiguities`](@ref).
"""
function test_all(
    testtarget::Module;
    ambiguities = (),
    # unbound_args = (),
    # undefined_exports = (),
)
    @testset "Method ambiguity" begin
        test_ambiguities([testtarget, Base]; ambiguities...)
    end
    @testset "Unbound type parameters" begin
        test_unbound_args(testtarget)
    end
    @testset "Undefined exports" begin
        test_undefined_exports(testtarget)
    end
end

end # module


text:
@static if Compat.Sys.isapple()
    const RingBufferSize = Int32
else
    const RingBufferSize = Clong
end

"""
    PaUtilRingBuffer(elementSizeBytes, elementCount)

Mirrors the C-side PaUtilRingBuffer struct that describes the ringbuffer state.
"""
mutable struct PaUtilRingBuffer
    bufferSize::RingBufferSize # Number of elements in FIFO. Power of 2. Set by PaUtil_InitializeRingBuffer.
    writeIndex::RingBufferSize # Index of next writable element. Set by PaUtil_AdvanceRingBufferWriteIndex.
    readIndex::RingBufferSize # Index of next readable element. Set by PaUtil_AdvanceRingBufferReadIndex.
    bigMask::RingBufferSize # Used for wrapping indices with extra bit to distinguish full/empty.
    smallMask::RingBufferSize # Used for fitting indices to buffer.
    elementSizeBytes::RingBufferSize # Number of bytes per element.
    buffer::Ptr{Cchar} # Pointer to the buffer containing the actual data.

    # allow it to be created uninitialized because we need to pass it to PaUtil_InitializeRingBuffer
    function PaUtilRingBuffer(elementSizeBytes, elementCount)
        data = Base.Libc.malloc(elementSizeBytes * elementCount)
        rbuf = new()
        PaUtil_InitializeRingBuffer(rbuf, elementSizeBytes, elementCount, data)

        @compat finalizer(close, rbuf)
        rbuf
    end
end

function close(rbuf::PaUtilRingBuffer)
    if rbuf.buffer != C_NULL
        Base.Libc.free(rbuf.buffer)
        rbuf.buffer = C_NULL
    end
end

isopen(rbuf::PaUtilRingBuffer) = rbuf.buffer != C_NULL

"""
    PaUtil_InitializeRingBuffer(rbuf, elementSizeBytes, elementCount, dataPtr)

Initialize Ring Buffer to empty state ready to have elements written to it.

# Arguments
* `rbuf::PaUtilRingBuffer`: The ring buffer.
* `elementSizeBytes::RingBufferSize`: The size of a single data element in bytes.
* `elementCount::RingBufferSize`: The number of elements in the buffer (must be a power of 2).
* `dataPtr::Ptr{Cvoid}`: A pointer to a previously allocated area where the data
  will be maintained.  It must be elementCount*elementSizeBytes long.
"""
function PaUtil_InitializeRingBuffer(rbuf, elementSizeBytes, elementCount, dataPtr)
    if !ispow2(elementCount)
        throw(ErrorException("elementCount($elementCount) must be a power of 2"))
    end
    status = ccall((:PaUtil_InitializeRingBuffer, libpa_ringbuffer),
                   RingBufferSize,
                   (Ref{PaUtilRingBuffer}, RingBufferSize, RingBufferSize, Ptr{Cvoid}),
                   rbuf, elementSizeBytes, elementCount, dataPtr)
    if status != 0
        throw(ErrorException("PaUtil_InitializeRingBuffer returned status $status"))
    end

    nothing
end

"""
    PaUtil_FlushRingBuffer(rbuf::PaUtilRingBuffer)

Reset buffer to empty. Should only be called when buffer is NOT being read or written.
"""
function PaUtil_FlushRingBuffer(rbuf)
    ccall((:PaUtil_FlushRingBuffer, libpa_ringbuffer),
          Cvoid,
          (Ref{PaUtilRingBuffer}, ),
          rbuf)
end

"""
    PaUtil_GetRingBufferWriteAvailable(rbuf::PaUtilRingBuffer)

Retrieve the number of elements available in the ring buffer for writing.
"""
function PaUtil_GetRingBufferWriteAvailable(rbuf)
    ccall((:PaUtil_GetRingBufferWriteAvailable, libpa_ringbuffer),
          RingBufferSize,
          (Ref{PaUtilRingBuffer}, ),
          rbuf)
end

"""
    PaUtil_GetRingBufferReadAvailable(rbuf::PaUtilRingBuffer)

Retrieve the number of elements available in the ring buffer for reading.
"""
function PaUtil_GetRingBufferReadAvailable(rbuf)
    ccall((:PaUtil_GetRingBufferReadAvailable, libpa_ringbuffer),
          RingBufferSize,
          (Ref{PaUtilRingBuffer}, ),
          rbuf)
end

"""
    PaUtil_WriteRingBuffer(rbuf::PaUtilRingBuffer,
                           data::Ptr{Cvoid},
                           elementCount::RingBufferSize)

Write data to the ring buffer and return the number of elements written.
"""
function PaUtil_WriteRingBuffer(rbuf, data, elementCount)
    ccall((:PaUtil_WriteRingBuffer, libpa_ringbuffer),
          RingBufferSize,
          (Ref{PaUtilRingBuffer}, Ptr{Cvoid}, RingBufferSize),
          rbuf, data, elementCount)
end

"""
    PaUtil_ReadRingBuffer(rbuf::PaUtilRingBuffer,
                          data::Ptr{Cvoid},
                          elementCount::RingBufferSize)

Read data from the ring buffer and return the number of elements read.
"""
function PaUtil_ReadRingBuffer(rbuf, data, elementCount)
    ccall((:PaUtil_ReadRingBuffer, libpa_ringbuffer),
          RingBufferSize,
          (Ref{PaUtilRingBuffer}, Ptr{Cvoid}, RingBufferSize),
          rbuf, data, elementCount)
end

max_stars_repo_name: wookay/TestRingBuffers

max_issues_repo_name: wookay/TestRingBuffers

max_forks_repo_name: wookay/TestRingBuffers


text:
import Random: randstring
import Dates

export Slider, NumberField, Button, CheckBox, TextField, PasswordField, Select, MultiSelect, Radio, FilePicker, DateField, TimeField, ColorStringPicker

struct Slider
    range::AbstractRange
    default::Number
    show_value::Bool
end

"""A Slider on the given `range`.

## Examples
`@bind x Slider(1:10)`

`@bind x Slider(0.00 : 0.01 : 0.30)`

`@bind x Slider(1:10; default=8, show_value=true)`

"""
Slider(range::AbstractRange; default=missing, show_value=false) = Slider(range, (default === missing) ? first(range) : default, show_value)

function show(io::IO, ::MIME"text/html", slider::Slider)
    print(io, """<input
        type="range"
        min="$(first(slider.range))"
        step="$(step(slider.range))"
        max="$(last(slider.range))"
        value="$(slider.default)"
        $(slider.show_value ? "oninput=\"this.nextElementSibling.value=this.value\"" : "")
        >""")

    if slider.show_value
        print(io, """<output>$(slider.default)</output>""")
    end
end

get(slider::Slider) = slider.default

"""A box where you can type in a number, within a specific range.

## Examples
`@bind x NumberField(1:10)`

`@bind x NumberField(0.00 : 0.01 : 0.30)`

`@bind x NumberField(1:10; default=8)`

"""
struct NumberField
    range::AbstractRange
    default::Number
end

NumberField(range::AbstractRange; default=missing) = NumberField(range, (default === missing) ? first(range) : default)

function show(io::IO, ::MIME"text/html", numberfield::NumberField)
    print(io, """<input type="number" min="$(first(numberfield.range))" step="$(step(numberfield.range))" max="$(last(numberfield.range))" value="$(numberfield.default)">""")
end

get(numberfield::NumberField) = numberfield.default

"""A button that sends back the same value every time that it is clicked.

You can use it to _trigger reactive cells_.

## Examples

In one cell:

```julia
@bind go Button("Go!")
```

and in a second cell:

```julia
begin
    # reference the bound variable - clicking the button will run this cell
    go

    md"My favorite number is $(rand())!"
end
```
"""
struct Button
    label::AbstractString
end
Button() = Button("Click")

function show(io::IO, ::MIME"text/html", button::Button)
    print(io, """<input type="button" value="$(htmlesc(button.label))">""")
end

get(button::Button) = button.label

"""A checkbox to choose a Boolean value `true`/`false`.

## Examples

`@bind programming_is_fun CheckBox()`

`@bind julia_is_fun CheckBox(default=true)`

`md"Would you like the thing? \$(@bind enable_thing CheckBox())"`
"""
struct CheckBox
    default::Bool
end

CheckBox(;default::Bool=false) = CheckBox(default)

function show(io::IO, ::MIME"text/html", button::CheckBox)
    print(io, """<input type="checkbox"$(button.default ? " checked" : "")>""")
end

get(checkbox::CheckBox) = checkbox.default

"""A text input (`<input type="text">`) - the user can type text, the text is returned as `String` via `@bind`.

If `dims` is a tuple `(cols::Integer, row::Integer)`, a `<textarea>` will be shown, with the given dimensions

Use `default` to set the initial value.

See the [Mozilla docs about `<input type="text">`](https://developer.mozilla.org/en-US/docs/Web/HTML/Element/input/text) and [`<textarea>`](https://developer.mozilla.org/en-US/docs/Web/HTML/Element/textarea)

# Examples
`@bind poem TextField()`

`@bind poem TextField((30,5); default="Hello\nJuliaCon!")`"""
struct TextField
    dims::Union{Tuple{Integer,Integer},Nothing}
    default::AbstractString
end
TextField(dims::Union{Tuple{Integer,Integer},Nothing}=nothing; default::AbstractString="") = TextField(dims, default)

function show(io::IO, ::MIME"text/html", textfield::TextField)
    if textfield.dims === nothing
        print(io, """<input type="text" value="$(htmlesc(textfield.default))">""")
    else
        print(io, """<textarea cols="$(textfield.dims[1])" rows="$(textfield.dims[2])">$(htmlesc(textfield.default))</textarea>""")
    end
end

get(textfield::TextField) = textfield.default

"""A password input (`<input type="password">`) - the user can type text, the text is returned as `String` via `@bind`.

This does not provide any special security measures, it just renders black dots (•••) instead of the typed characters.

Use `default` to set the initial value.

See the [Mozilla docs about `<input type="password">`](https://developer.mozilla.org/en-US/docs/Web/HTML/Element/input/password)

# Examples
`@bind secret_poem PasswordField()`

`@bind secret_poem PasswordField(default="Te dansen omdat men leeft")`"""
Base.@kwdef struct PasswordField
    default::AbstractString=""
end

function show(io::IO, ::MIME"text/html", passwordfield::PasswordField)
    print(io, """<input type="password" value="$(htmlesc(passwordfield.default))">""")
end

get(passwordfield::PasswordField) = passwordfield.default

"""A dropdown menu (`<select>`) - the user can choose one of the `options`, an array of `String`s.

See [`MultiSelect`](@ref) for a version that allows multiple selected items.

`options` can also be an array of pairs `key::String => value::Any`. The `key` is returned via `@bind`; the `value` is shown.

See the [Mozilla docs about `select`](https://developer.mozilla.org/en-US/docs/Web/HTML/Element/select)

# Examples
`@bind veg Select(["potato", "carrot"])`

`@bind veg Select(["potato" => "🥔", "carrot" => "🥕"])`

`@bind veg Select(["potato" => "🥔", "carrot" => "🥕"], default="carrot")`"""
struct Select
    options::Vector{Pair{<:AbstractString,<:Any}}
    default::Union{Missing, AbstractString}
end
Select(options::AbstractVector{<:AbstractString}; default=missing) = Select([o => o for o in options], default)
Select(options::AbstractVector{<:Pair{<:AbstractString,<:Any}}; default=missing) = Select(options, default)

function show(io::IO, ::MIME"text/html", select::Select)
    withtag(io, :select) do
        for o in select.options
            print(io, """<option value="$(htmlesc(o.first))"$(select.default === o.first ? " selected" : "")>""")
            if showable(MIME"text/html"(), o.second)
                show(io, MIME"text/html"(), o.second)
            else
                print(io, o.second)
            end
            print(io, "</option>")
        end
    end
end

get(select::Select) = ismissing(select.default) ? first(select.options).first : select.default

"""A multi-selector (`<select multi>`) - the user can choose one or more of the `options`, an array of `Strings.

See [`Select`](@ref) for a version that allows only one selected item.

`options` can also be an array of pairs `key::String => value::Any`. The `key` is returned via `@bind`; the `value` is shown.

See the [Mozilla docs about `select`](https://developer.mozilla.org/en-US/docs/Web/HTML/Element/select)

# Examples
`@bind veg MultiSelect(["potato", "carrot"])`

`@bind veg MultiSelect(["potato" => "🥔", "carrot" => "🥕"])`

`@bind veg MultiSelect(["potato" => "🥔", "carrot" => "🥕"], default=["carrot"])`"""
struct MultiSelect
    options::Vector{Pair{<:AbstractString,<:Any}}
    default::Union{Missing, AbstractVector{AbstractString}}
end
MultiSelect(options::AbstractVector{<:AbstractString}; default=missing) = MultiSelect([o => o for o in options], default)
MultiSelect(options::AbstractVector{<:Pair{<:AbstractString,<:Any}}; default=missing) = MultiSelect(options, default)

function show(io::IO, ::MIME"text/html", select::MultiSelect)
    withtag(io, Symbol("select multiple")) do
        for o in select.options
            print(io, """<option value="$(htmlesc(o.first))"$(!ismissing(select.default) && o.first ∈ select.default ? " selected" : "")>""")
            if showable(MIME"text/html"(), o.second)
                show(io, MIME"text/html"(), o.second)
            else
                print(io, o.second)
            end
            print(io, "</option>")
        end
    end
end

get(select::MultiSelect) = ismissing(select.default) ? Any[] : select.default

"""A file upload box. The chosen file will be read by the browser, and the bytes are sent back to Julia.

The optional `accept` argument can be an array of `MIME`s. The user can only select files with these MIME. If only `image/*` MIMEs are allowed, then smartphone browsers will open the camera instead of a file browser.

## Examples

`@bind file_data FilePicker()`

`file_data["data"]`

You can limit the allowed MIME types:

```julia
@bind image_data FilePicker([MIME("image/jpg"), MIME("image/png")])
# and use MIME groups:
@bind image_data FilePicker([MIME("image/*")])
```
"""
struct FilePicker
    accept::Array{MIME,1}
end
FilePicker() = FilePicker(MIME[])

function show(io::IO, ::MIME"text/html", filepicker::FilePicker)
    print(io, """<input type='file' accept='""")
    join(io, string.(filepicker.accept), ",")
    print(io, "'>")
end

get(select::FilePicker) = Dict("name" => "", "data" => UInt8[], "type" => "")

"""A group of radio buttons - the user can choose one of the `options`, an array of `String`s.

`options` can also be an array of pairs `key::String => value::Any`. The `key` is returned via `@bind`; the `value` is shown.

# Examples
`@bind veg Radio(["potato", "carrot"])`

`@bind veg Radio(["potato" => "🥔", "carrot" => "🥕"])`

`@bind veg Radio(["potato" => "🥔", "carrot" => "🥕"], default="carrot")`

"""
struct Radio
    options::Vector{Pair{<:AbstractString,<:Any}}
    default::Union{Missing, AbstractString}
end
Radio(options::AbstractVector{<:AbstractString}; default=missing) = Radio([o => o for o in options], default)
Radio(options::AbstractVector{<:Pair{<:AbstractString,<:Any}}; default=missing) = Radio(options, default)

function show(io::IO, ::MIME"text/html", radio::Radio)
    groupname = randstring('a':'z')
    withtag(io, :form, :id=>groupname) do
        for o in radio.options
            withtag(io, :div) do
                print(io, """<input type="radio" id="$(htmlesc(groupname * o.first))" name="$(groupname)" value="$(htmlesc(o.first))"$(radio.default === o.first ? " checked" : "")>""")

                withtag(io, :label, :for=>(groupname * o.first)) do
                    if showable(MIME"text/html"(), o.second)
                        show(io, MIME"text/html"(), o.second)
                    else
                        print(io, o.second)
                    end
                end
            end
        end
    end
    withtag(io, :script) do
        print(io, """
        const form = document.querySelector('#$(groupname)')

        form.oninput = (e) => {
            form.value = e.target.value
            // and bubble upwards
        }

        // set initial value:
        const selected_radio = form.querySelector('input[checked]')
        if(selected_radio != null){
            form.value = selected_radio.value
        }
        """)
    end
end

get(radio::Radio) = radio.default

"""A date input (`<input type="date">`) - the user can pick a date, the date is returned as `Dates.DateTime` via `@bind`.

Use `default` to set the initial value.

See the [Mozilla docs about `<input type="date">`](https://developer.mozilla.org/en-US/docs/Web/HTML/Element/input/date)

# Examples
`@bind best_day_of_my_live DateField()`

`@bind best_day_of_my_live DateField(default=today())`"""
Base.@kwdef struct DateField
    default::Union{Dates.TimeType,Missing}=missing
end

function show(io::IO, ::MIME"text/html", datefield::DateField)
    withtag(() -> (), io, :input, :type=>"date", :value=>datefield.default === missing ? "" : Dates.format(datefield.default, "Y-mm-dd"))
end
get(datefield::DateField) = datefield.default

"""A time input (`<input type="time">`) - the user can pick a time, the time is returned as `Dates.DateTime` via `@bind`.

Use `default` to set the initial value.

See the [Mozilla docs about `<input type="time">`](https://developer.mozilla.org/en-US/docs/Web/HTML/Element/input/time)

# Examples
`@bind lunch_time TimeField()`

`@bind lunch_time TimeField(default=now())`"""
Base.@kwdef struct TimeField
    default::Union{Dates.TimeType,Missing}=missing
end

function show(io::IO, ::MIME"text/html", timefield::TimeField)
    withtag(() -> (), io, :input, :type=>"time", :value=>timefield.default === missing ? "" : Dates.format(timefield.default, "HH:MM:SS"))
end
get(timefield::TimeField) = timefield.default

"""A color input (`<input type="color">`) - the user can pick an RGB color, the color is returned as color hex `String` via `@bind`. The value is lowercase and starts with `#`.

Use `default` to set the initial value.

See the [Mozilla docs about `<input type="color">`](https://developer.mozilla.org/en-US/docs/Web/HTML/Element/input/color)

# Examples
`@bind color ColorStringPicker()`

`@bind color ColorStringPicker(default="#aabbcc")`
"""
Base.@kwdef struct ColorStringPicker
    default::String="#000000"
end

function show(io::IO, ::MIME"text/html", colorStringPicker::ColorStringPicker)
    withtag(() -> (), io, :input, :type=>"color", :value=>colorStringPicker.default)
end
get(colorStringPicker::ColorStringPicker) = colorStringPicker.default


text:
using Compat.Test
using CodecLz4

@testset "CodecLz4.jl" begin
    include("lz4frame.jl")
    include("stream_compression.jl")
end


text:
if VERSION >= v"0.7.0"
	import Random: randperm, shuffle!, shuffle, seed!
	import Printf.@printf
	import LinearAlgebra
	import LinearAlgebra: norm, Diagonal, dot
	import Statistics: var, mean, std
	using UnicodePlots
	import Base.minimum
end


text:
# make a taxonomy tree of all the birds on NCBI
using Pkg; Pkg.activate(".")
using Taxonomy
using NCBITaxonomy

procs = ncbi"Procellariiformes"
birds = ncbi"Aves"

birdtaxonomy = taxonomy(NCBITaxonWrapper(birds), FixedRankTaxonomy)
proctaxonomy = taxonomy(NCBITaxonWrapper(procs), FixedRankTaxonomy)

# things that should be true
birdtaxonomy ∩ proctaxonomy == proctaxnomy
birdtaxonomy ∪ proctaxonomy == birdtaxonomy
birdtaxonomy ∈ proctaxonomy == false
proctaxonomy ∈ birdtaxonomy == true


text:
module HyperspectraWithNeXL

using DrWatson
using Reexport
@reexport using DataDeps

# Register hyperspectral data and standards
function __init__()
    (!isdir(datadir(),"exp_raw")) && mkdirs(datadir(),"exp_raw")
    ENV["DATADEPS_LOAD_PATH"] = mkpath(joinpath(datadir(),"exp_raw"))
    register(DataDep("MnNodule",
        """
        Dataset: Deep sea manganese nodule electron excited X-ray microanalysis hyperspectral data set
         Author: Nicholas W. M. Ritchie (NIST)
        License: Public Domain
        Website: https://data.nist.gov/od/id/mds2-2467
         Notice: This file is over 600 Mb
        """,
        "https://data.nist.gov/od/ds/mds2-2467/MnNodule.tar.gz",
        "5b5b6623b8f4daca3ff3073708442ac5702ff690aa12668659875ec5642b458d",
        post_fetch_method = unpack
    ))
    register(DataDep("MnNodule_Standards",
        """
        Dataset: Standard spectra for the deep sea manganese nodule electron excited X-ray microanalysis hyperspectral data set
         Author: Nicholas W. M. Ritchie (NIST)
        License: Public Domain
        Website: https://data.nist.gov/od/id/mds2-2467
        """,
        "https://data.nist.gov/od/ds/mds2-2467/MnNodule_Standards.tar.gz",
        "69283ba72146932ba451e679cf02fbd6b350f96f6d012d50f589ed9dd2e35f1a",
        post_fetch_method = unpack
    ))
end

# Nothing to export
include("tohdf5.jl")

end


text:
using Markdown

function lx_escape(com, _)
    # keep this first line
    content = Franklin.content(com.braces[1]) # input string
    lang, code = split(content, "::")
    scode = strip(code)
    esc_code = Markdown.htmlesc(scode)
    io = IOBuffer()
    println(io, "```$lang")
    println(io, esc_code)
    println(io, "```")
    return String(take!(io))
end

max_stars_repo_name: mcognetta/franklindocs

max_issues_repo_name: mcognetta/franklindocs

max_forks_repo_name: mcognetta/franklindocs


text:
"""
Nelder-mead enhanced extreme learning machine.

Philip Reiner and Bogdan M Wilamowski.

2013 IEEE 17th International Conference on Intelligent Engineering Systems
(INES), 2013 pp. 225-230.

http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=6632816

"""
mutable struct NMEELM{TV<:AbstractArray{Float64}} <: AbstractSLFN
    p::Int  # Number of training points
    q::Int  # Dimensionality of function domain
    s::Int  # number of neurons
    ϵ::Float64  # expected learning accuracy

    # Nelder Mead parameters
    ρ::Float64
    χ::Float64
    γ::Float64
    α::Float64

    # Number of simplex operations
    k::Int

    # internal parameters
    Wt::Matrix{Float64}  # transpose of W matrix
    d::Vector{Float64}
    activation::Identity
    neuron_type::RBF{Gaussian}
    β::TV

    function NMEELM(p::Int, q::Int, s::Int, ϵ::Float64,
                    ρ::Float64, χ::Float64, γ::Float64, α::Float64)
        Wt = 2*rand(q, s) - 1
        d = rand(s)
        β = TV(s)
        new(p, q, s, ϵ, ρ, χ, γ, α, Wt, d, Identity(), RBF(Gaussian), β)
    end
end

function NMEELM{TV<:AbstractArray}(x::AbstractArray, y::TV; s::Int=size(y, 1),
                                   ϵ::Float64=1e-6, ρ::Float64=1.0, χ::Float64=2.0,
                                   γ::Float64=1/2, α::Float64=1/2, k::Int=7)
    q = size(x, 2)  # dimensionality of function domain
    p = size(y, 1)  # number of training points
    s = min(p, s)   # Can't have more neurons than training points
    out = NMEELM{TV}(p, q, s, ϵ, ρ, χ, γ, α, k)
    fit!(out, x, y)
end

## API methods
isexact(elm::NMEELM) = false
input_to_node(elm::NMEELM, x::AbstractArray) = input_to_node(elm.neuron_type, x, elm.Wt, elm.d)
hidden_out(elm::NMEELM, x::AbstractArray) = elm.activation(input_to_node(elm, x))

function fit!(elm::NMEELM, x::AbstractArray, y::AbstractArray)
    # ------- Step 1: initialize
    n = 0
    E = copy(y)
    err = Inf
    xt = x'

    # ------- Step 2: Learning
    while n < elm.s && err > elm.ϵ
        # --- Step 2.a: Increase n
        n += 1

        # --- Step 2.b: find index of max error
        j = indmax(E)

        # --- Step 2.c: Assign new center to be input x
        cn = xt[:, j]
        elm.Wt[:, n] = cn

        # --- Step 2.d: Initialze βₙ = Eⱼ
        βn = E[j]

        # --- Step 2.e: Do Nelder-Mead simplex algorithm
        # define objective function
        function obj(σ)
            sse = 0.0
            for i in 1:elm.p
                gi = input_to_node(elm.neuron_type, view(xt, :, i), cn, σ)
                sse += (E[i] - βn * gi)^2
            end
            sse
        end

        # call to elm.k nm simplex iterations (or call Optim??)
        # TODO: pick up here

    end

end

@compat function (elm::NMEELM)(x′::AbstractArray)
    @assert size(x′, 2) == elm.q "wrong input dimension"
    return hidden_out(elm, x′) * elm.v
end

function Base.show{TA}(io::IO, elm::NMEELM{TA})
    s =
    """
    NMEELM with
      - Identity Activation function
      - $(elm.q) input dimension(s)
      - RBF{Gaussian} neuron(s)
      - $(elm.p) training point(s)
    """
    print(io, s)
end


text:
"""
    solve(
        problem::GeneralBTPDE,
        gradient,
        odesolver = QNDF();
        abstol = 1e-6,
        reltol = 1e-4,
        callbacks = [].
    )

Solve the Bloch-Torrey partial differential equation using P1 finite elements in space and
`odesolver` in time.
"""
function solve(
    problem::BTPDE{T},
    gradient,
    odesolver::OrdinaryDiffEqAlgorithm = QNDF(autodiff = false);
    abstol = 1e-6,
    reltol = 1e-4,
    callbacks = AbstractCallback[],
) where {T}
    (; model, matrices) = problem
    (; γ) = model
    (; M, S, R, Mx, Q) = matrices

    ρ = initial_conditions(model)

    # Time dependent ODE function
    function Mdξ!(dξ, ξ, p, t)
        Jac!(p.J, ξ, p, t)
        mul!(dξ, p.J, ξ)
    end

    # Time dependent Jacobian of ODE function with respect to the state `ξ`
    function Jac!(J, ξ, p, t)
        (; S, Q, R, gradient) = p
        g⃗ = gradient(t)
        @. J = -(S + Q + R + im * γ * (g⃗[1] * Mx[1] + g⃗[2] * Mx[2] + g⃗[3] * Mx[3]))
    end

    function func(u, t, integrator)
        for cb ∈ callbacks
            update!(cb, problem, gradient, u, t)
        end
    end

    # Jacobian sparsity pattern
    jac_prototype = -(S + Q + im * sum(Mx))

    # ODE problem
    J = copy(jac_prototype)
    p = (; J, S, Q, R, gradient)
    TE = echotime(gradient)

    for cb ∈ callbacks
        initialize!(cb, problem, gradient, ρ, 0)
    end

    odefunction = ODEFunction(Mdξ!; mass_matrix = M, jac = Jac!, jac_prototype)
    odeproblem = ODEProblem(odefunction, ρ, (0, TE), p, progress = false)

    callback = FunctionCallingCallback(func; func_start = false)

    # Solve ODE problem
    sol = OrdinaryDiffEq.solve(
        odeproblem,
        odesolver;
        reltol,
        abstol,
        callback,
        save_everystep = false,
    )

    finalize!.(callbacks)

    sol.u[end]
end


text:
using TruncatedNormal, Test, SpecialFunctions
using TruncatedNormal: xerfcx_asym, xerfcx_asym_pi

@test isnan(erf(NaN, NaN))
@test isnan(erf(0, NaN))
@test isnan(erf(NaN, 0))

@test erf(0, 0) == 0
@test erf(+1, +1) == 0
@test erf(-1, -1) == 0

@test erf(0, +Inf) == +1
@test erf(-Inf, 0) == +1
@test erf(+Inf, 0) == -1
@test erf(0, -Inf) == -1

@test erf(-1, 1) ≈ 1.68540158589942973868244127017
@test erf(-2, 2) ≈ 1.99064453003790546832413851273
@test erf(-5, 5) ≈ 1.99999999999692508041114393030

@test xerfcx_asym(500, 3) ≈ 500erfcx(500)
@test xerfcx_asym(1e6, 3) ≈ 1e6erfcx(1e6)
@test xerfcx_asym_pi(500, 3) ≈ 1.12837239688821180302500663804e-6
@test xerfcx_asym_pi(9e5, 3) ≈ 3.48265175028834253964811258438e-13
@test xerfcx_asym_pi(1e6, 3) ≈ 2.82094791773455001286379966421e-13


text:
##########################################################################
#
# Benchmark optimization algorithms by tracking:
# * Run time over 100 runs -- after 1 initial run that forces JIT.
# * Number of iterations
# * Euclidean error of solution
# * Memory requirements (TODO)
#
##########################################################################

##########################################################################
#
# Load the software and test functions.
#
##########################################################################

using Optim
include("test_functions.jl")

##########################################################################
#
# Iterate over 5 optimizable functions:
# * parabola
# * powell
# * rosenbrock
# * polynomial
# * exponential
#
##########################################################################

problems = Array(Any, 6)

parabola_problem = Dict()
parabola_problem[:name] = "Parabola"
parabola_problem[:f] = parabola
parabola_problem[:g] = parabola_gradient!
parabola_problem[:h] = parabola_hessian!
parabola_problem[:initial_x] = [0.0, 0.0, 0.0, 0.0, 0.0]
parabola_problem[:solution] = [1.0, 2.0, 3.0, 5.0, 8.0]
problems[1] = parabola_problem

powell_problem = Dict()
powell_problem[:name] = "Powell"
powell_problem[:f] = powell
powell_problem[:g] = powell_gradient!
powell_problem[:h] = powell_hessian!
powell_problem[:initial_x] = [3.0, -1.0, 0.0, 1.0]
powell_problem[:solution] = [0.0, 0.0, 0.0, 0.0]
problems[2] = powell_problem

rosenbrock_problem = Dict()
rosenbrock_problem[:name] = "Rosenbrock"
rosenbrock_problem[:f] = rosenbrock
rosenbrock_problem[:g] = rosenbrock_gradient!
rosenbrock_problem[:h] = rosenbrock_hessian!
rosenbrock_problem[:initial_x] = [0.0, 0.0]
rosenbrock_problem[:solution] = [1.0, 1.0]
problems[3] = rosenbrock_problem

polynomial_problem = Dict()
polynomial_problem[:name] = "Polynomial"
polynomial_problem[:f] = polynomial
polynomial_problem[:g] = polynomial_gradient!
polynomial_problem[:h] = polynomial_hessian!
polynomial_problem[:initial_x] = [0.0, 0.0, 0.0]
polynomial_problem[:solution] = [10.0, 7.0, 108.0]
problems[4] = polynomial_problem

exponential_problem = Dict()
exponential_problem[:name] = "Exponential"
exponential_problem[:f] = exponential
exponential_problem[:g] = exponential_gradient!
exponential_problem[:h] = exponential_hessian!
exponential_problem[:initial_x] = [0.0, 0.0]
exponential_problem[:solution] = [2.0, 3.0]
problems[5] = exponential_problem

large_polynomial_problem = Dict()
large_polynomial_problem[:name] = "Large Polynomial"
large_polynomial_problem[:f] = large_polynomial
large_polynomial_problem[:g] = large_polynomial_gradient!
large_polynomial_problem[:h] = large_polynomial_hessian!
large_polynomial_problem[:initial_x] = zeros(250)
large_polynomial_problem[:solution] = [1.0:1.0:250.0]
problems[6] = large_polynomial_problem

##########################################################################
#
# Iterate over 6 optimization functions:
# * gradient_descent: Gradient Descent
# * newton: Newton's Method
# * bfgs: BFGS
# * l_bfgs: L-BFGS
# * nelder_mead: Nelder-Mead
# * sa: Simulated Annealing
#
##########################################################################

algorithms = [:gradient_descent,
              :newton,
              :bfgs,
              :l_bfgs,
              :nelder_mead,
              :simulated_annealing]

# Print out a header line for the TSV-formatted report.
println(join({"Problem",
              "Algorithm",
              "AverageRunTimeInMilliseconds",
              "Iterations",
              "Error"},
             "\t"))

for problem = problems
  for algorithm = algorithms
    # Force compilation
    results = optimize(problem[:f],
                       problem[:g],
                       problem[:h],
                       problem[:initial_x],
                       method = algorithm,
                       tolerance = 1e-16)

    # Run each algorithm 100 times
    n = 100

    # Estimate run time in milliseconds
    run_time = @elapsed for i = 1:n
      results = optimize(problem[:f],
                         problem[:g],
                         problem[:h],
                         problem[:initial_x],
                         method = algorithm,
                         tolerance = 1e-16)
    end
    run_time = run_time * n

    # Estimate error in discovered solution
    results = optimize(problem[:f],
                       problem[:g],
                       problem[:h],
                       problem[:initial_x],
                       method = algorithm,
                       tolerance = 1e-16)
    errors = norm(results.minimum - problem[:solution])

    # Count iterations
    iterations = results.iterations

    # Print out results.
    println(join({problem[:name],
                  results.method,
                  run_time / n,
                  iterations,
                  errors},
                  "\t"))
  end
end


text:
using Test
using ForwardDiff2: Dual, partials, dualrun, find_dual, Tag

const Tag1 = Tag{Nothing}
const Tag2 = Tag{Tag{Nothing}}
const tag1 = Tag1()
const tag2 = Tag2()

@testset "find_dual" begin
    @test find_dual(tag1, 0, 0) == 0
    @test find_dual(tag1, Dual{Tag1}(1,1), 1) == 1
    @test find_dual(tag1, 1, Dual{Tag1}(1,1)) == 2
    @test find_dual(tag2, 1, Dual{Tag1}(1,1)) == 0
    @test find_dual(tag2, Dual{Tag1}(1,1), 1) == 0
    @test find_dual(tag1, Dual{Tag2}(1,1), 1) == 0
end

using Cassette
Cassette.@context TestCtx

const TaggedTestCtx{T} = Cassette.Context{Cassette.nametype(TestCtx), T}

@inline function find_dual_ctx(::TaggedTestCtx{T}, args...) where T
    find_dual(Tag{T}(), args...)
end

@testset "find_dual_ctx" begin
    ctx = TaggedTestCtx(metadata=nothing)
    @test find_dual_ctx(ctx, 1, Dual{Tag1}(1,1)) == 2
    @test find_dual_ctx(ctx, Dual{Tag1}(1,1), 1) == 1
    @test find_dual_ctx(ctx, Dual{Tag1}(1,1), Dual{Tag2}(1,1)) == 1
    @test find_dual_ctx(ctx, Dual{Tag2}(1,1), Dual{Tag1}(1,1)) == 2

    ctx = TaggedTestCtx(metadata=Tag{Nothing}())
    @test find_dual_ctx(ctx, Dual{Tag1}(1,1), Dual{Tag2}(1,1)) == 2
    @test find_dual_ctx(ctx, Dual{Tag2}(1,1), Dual{Tag1}(1,1)) == 1
end


text:
module ShaleDrillingLikelihood_NEWFlow_Test

using ShaleDrillingLikelihood
using Test
using Calculus
using Random

using InteractiveUtils

using Base.Iterators: product
using ShaleDrillingLikelihood: showtypetree,
    ObservationDrill,
    Observation,
    SimulationDraw,
    check_flow_grad,
    TestDrillModel,
    _nparm,
    flow!,
    theta_ρ,
    vw_revenue,
    vw_cost,
    vw_extend,
    revenue

using Calculus: finite_difference!

const DOPRINT = false

@testset "Drilling Reward" begin

    @testset "Flow gradients" begin
        problem = DrillReward(
            DrillingRevenue(Unconstrained(),NoTrend(),NoTaxes()),
            DrillingCost_constant(),
            ExtensionCost_Constant()
        )

        Random.seed!(1234)

        DOPRINT && showtypetree(AbstractPayoffFunction)
        @test true == true

        println("")

        payoff_functions_to_test = (
            DrillingCost_constant(),
            DrillingCost_dgt1(),
            DrillingCost_TimeFE(2008,2012),
            DrillingCost_TimeFE(2009,2011),
            DrillingCost_TimeFE_rigrate(2008,2012),
            DrillingCost_DoubleTimeFE(2008,2012),

            ScrapValue_Constant(),
            ScrapValue_Price(),

            DrillingRevenue(Constrained(),NoTrend(),NoTaxes(), Learn(), WithRoyalty()),
            DrillingRevenue(Constrained(),NoTrend(),WithTaxes(), Learn(), WithRoyalty()),
            DrillingRevenue(Constrained(),TimeTrend(),NoTaxes(), Learn(), WithRoyalty()),
            DrillingRevenue(Constrained(),TimeTrend(),WithTaxes(), Learn(), WithRoyalty()),
            DrillingRevenue(Unconstrained(),NoTrend(),NoTaxes(), Learn(), WithRoyalty()),
            DrillingRevenue(Unconstrained(),NoTrend(),WithTaxes(), Learn(), WithRoyalty()),
            DrillingRevenue(Unconstrained(),TimeTrend(),NoTaxes(), Learn(), WithRoyalty()),
            DrillingRevenue(Unconstrained(),TimeTrend(),WithTaxes(), Learn(), WithRoyalty()),

            DrillingRevenue(Unconstrained(),TimeTrendDgt0(),NoTaxes(), Learn(), WithRoyalty()),
            DrillingRevenue(Unconstrained(),TimeFE(),WithTaxes(), Learn(), WithRoyalty()),

            DrillingRevenue(Constrained(),NoTrend(),NoTaxes(), Learn(), NoRoyalty()),
            DrillingRevenue(Constrained(),NoTrend(),WithTaxes(), Learn(), NoRoyalty()),
            DrillingRevenue(Constrained(),TimeTrend(),NoTaxes(), Learn(), NoRoyalty()),
            DrillingRevenue(Constrained(),TimeTrend(),WithTaxes(), Learn(), NoRoyalty()),
            DrillingRevenue(Unconstrained(),NoTrend(),NoTaxes(), Learn(), NoRoyalty()),
            DrillingRevenue(Unconstrained(),NoTrend(),WithTaxes(), Learn(), NoRoyalty()),
            DrillingRevenue(Unconstrained(),TimeTrend(),NoTaxes(), Learn(), NoRoyalty()),
            DrillingRevenue(Unconstrained(),TimeTrend(),WithTaxes(), Learn(), NoRoyalty()),

            ExtensionCost_Zero(),
            ExtensionCost_Constant(),
            DrillReward(DrillingRevenue(Unconstrained(),NoTrend(),NoTaxes()), DrillingCost_TimeFE(2009,2011), ExtensionCost_Constant(), ScrapValue_Zero()), # UnconstrainedProblem( StaticDrillingPayoff(DrillingRevenue(Unconstrained(),NoTrend(),NoTaxes()), DrillingCost_TimeFE(2009,2011), ExtensionCost_Constant()), ),
            DrillReward(DrillingRevenue(Unconstrained(),NoTrend(),NoTaxes()), DrillingCost_TimeFE(2009,2011), ExtensionCost_Constant(), ScrapValue_Constant()), # UnconstrainedProblem( StaticDrillingPayoff(DrillingRevenue(Unconstrained(),NoTrend(),NoTaxes()), DrillingCost_TimeFE(2009,2011), ExtensionCost_Constant()), ),
            DrillReward(DrillingRevenue(Unconstrained(),NoTrend(),NoTaxes()), DrillingCost_TimeFE(2009,2011), ExtensionCost_Constant(), ScrapValue_Price()), # UnconstrainedProblem( StaticDrillingPayoff(DrillingRevenue(Unconstrained(),NoTrend(),NoTaxes()), DrillingCost_TimeFE(2009,2011), ExtensionCost_Constant()), ),
            ConstrainedProblem(   DrillReward(DrillingRevenue(Unconstrained(),NoTrend(),NoTaxes()), DrillingCost_TimeFE(2009,2011), ExtensionCost_Constant()), ),
            DrillReward(DrillingRevenue(Constrained(),NoTrend(),NoTaxes()), DrillingCost_constant(), ExtensionCost_Constant()),
        )

        for f in payoff_functions_to_test
            strout = "Testing $f"
            strout = replace(strout, "ShaleDrillingLikelihood." => "")
            DOPRINT && println(strout)
            let z = (2.5, 2.0, 2010), ichars = (4.5, 0.25)

                n = _nparm(f)
                cfs = coefnames(f)
                @test cfs isa Vector{<:Union{String,Symbol}}
                @test n == length(cfs)
                θ0 = rand(n)
                u, v = randn(2)

                if n > 0
                    sim = SimulationDraw(u, v, last(θ0))
                else
                    sim = SimulationDraw(u, v, 0.0)
                end

                # m = TestDrillModel()
                # (d,i) = (1,1)
                # theta = θ0
                # s = sim
                # obs = ObservationDrill(m, ichars, z, d, i)
                # g = zero(theta)
                # @code_warntype flow!(g, f, d, obs, theta, s, true)

                for (d,i) in product(0:2, 1:3)

                    m = TestDrillModel()
                    obs = ObservationDrill(m, ichars, z, d, i)

                    check_flow_grad(f, d, obs, θ0, sim)
                end

            end
        end
    end

    @testset "parameter access" begin

        discount = ((0x1.006b55c832502p+0)^12 / 1.125) ^ (1/4)  # real discount rate

        # set up coefs
        θρ =  0.0
        αψ = 0.33
        αg = 0.56

        # set up coefs
        #            θρ  ψ        g    x_r     κ₁   κ₂
        θ_royalty = [θρ, 1.0,    0.1, -0.3,  -0.6, 0.6]
        #                αψ, αg, γx   σ2η, σ2u   # η is iwt, u is iw
        θ_produce = vcat(αψ, αg, 0.2, 0.3, 0.4)
        #            drill  ext    α0  αg  αψ  θρ
        θ_drill_u = [-5.5, -2.0, -2.8, αg, αψ, θρ]
        # θ_drill_u = [-5.8,       -2.8, αg, αψ, θρ]
        θ_drill_c = vcat(θ_drill_u[1:3], θρ)

        # model
        rwrd_c = DrillReward(DrillingRevenue(Constrained(;log_ogip=αg, α_ψ=αψ),NoTrend(),NoTaxes()), DrillingCost_constant(), ExtensionCost_Constant())
        rwrd_u = UnconstrainedProblem(rwrd_c)
        rwrd = rwrd_u
        @test _nparm(rwrd) == length(θ_drill_u)

        # test drilling coef
        @test theta_ρ(rwrd, θ_drill_u) == θρ
        @test theta_ρ(revenue(rwrd), vw_revenue(rwrd, θ_drill_u)) == θρ
        @test θ_drill_u[1:1] == vw_cost(rwrd, θ_drill_u)
        @test θ_drill_u[2:2] == vw_extend(rwrd, θ_drill_u)
        @test θ_drill_u[end-3:end] == vw_revenue(rwrd, θ_drill_u)

    end

end

end # Drilling Model Flow Payoffs


text:
import NetCDF.nccreate, NetCDF.ncwrite, NetCDF.ncclose
# writes the first 1000 extremes of an extremelist to a nectdf file
function writeExtremes(el::ExtremeList,folder::String;nmax=1000,prefix="e")
    run(`mkdir -p $(folder)`)
    for iex=1:min(nmax,length(el.extremes))
        si=@sprintf("%04d",iex)
        filename="$(folder)/$(prefix)$(si).nc"
        isfile(filename) && rm(filename)
        nccreate(filename,"locs",{"units"=>"Lon,Lat,Time"},"pixel",[1:size(el.extremes[iex].locs,1)],"dimension",[1,2,3,4])
        ncwrite(el.extremes[iex].locs,filename,"locs",start=[1,1],count=[-1,3])
        ncwrite(el.extremes[iex].zvalues,filename,"locs",start=[1,4],count=[-1,1])
        ncclose()
    end
end

function unroll_tuples(f::Vector)
    nwrites=Int[length(fi[1]) for fi in f]
    maximum(nwrites)>5 && error("Do not use the writeFeatures function to write time series, use it only for single indices")
    nEx=length(f[1])
    fnew=Array(Any,sum(nwrites))
    inew=1
    for i=1:length(f)
        l=length(f[i][1])
        if l==1
            #Everything is ok
            fnew[inew]=f[i]
            inew=inew+1
        else
            for j=1:l
                fnew[inew]=[f[i][k][j] for k=1:nEx]
                inew=inew+1
            end
        end
    end
    return fnew
end

function writeFeatures(f::Vector,featnames::Vector,file;overwrite=true)
    fnew=unroll_tuples(f)
    length(fnew)==length(featnames) || error("Number of feature names must match number of written features. Make sure that you named each quantile separately!")
    featDict=Dict{Any,Any}()
    [featDict["$i"]="$(featnames[i])" for i=1:length(featnames)]
    #Write features
    isfile(file) && overwrite && rm(file)
    nccreate(file,"features",{"unit"=>"diverse"},"N_Extreme",[1:length(fnew[1])],"Feat",[1:length(featnames)],featDict)
    for i=1:length(fnew)
            ncwrite(fnew[i],file,"features",start=[1,i],count=[-1,1])
    end
    ncclose()
end

writeTimeSeries(f::Vector,filename,varNames::String)=writeTimeSeries(f,filename,[varNames])
function writeTimeSeries(f::Vector,filename,varNames::Vector)
    isfile(filename) && rm(filename)
    ntstep=length(f[1][1])
    length(varNames)==length(f) || error("number of variable names must equal")
    for i=1:length(f)
        nccreate(filename,varNames[i],{"unit"=>"none"},"timestep",[1:ntstep])
    end
    for i=1:length(f)
        ncwrite(f[i][1],filename,varNames[i])
    end
    ncclose()
end

max_stars_repo_name: CAB-LAB/ExtremeStats

max_issues_repo_name: CAB-LAB/ExtremeStats

max_forks_repo_name: CAB-LAB/ExtremeStats


text:
## Sequential Search / Linear Search

function search(list, num)
for i in list
        if num == i
            println("Found it")
            break
        end
    end
end

search([1,2,3,1,2,2,3,1,56,43,43,4,3,4,3,6,36,43,], 50)


text:
# Uses obs_eq, state_eq

export simulate

function simulate(m, x0, w, v; squeeze=true)
    @assert length(w) == length(v) + 1
    state = similar(w, length(x0), length(w))
    state[:, 1] .= x0
    first_obs = obs_eq(m, view(state, :, 1), w[1])
    obs = similar(w, length(first_obs), length(w))
    obs[:, 1] .= first_obs
    for i in eachindex(v)
        obs[:, i] .= obs_eq(m, view(state, :, i), w[i])
        state[:, i+1] .= state_eq(m, view(state, :, i), v[i])
    end
    obs[:, end] .= obs_eq(m, view(state, :, size(state)[2]), w[end])
    if squeeze && size(obs)[1] == 1
        obs = reshape(obs, :)
    end
    if squeeze && size(state)[1] == 1
        state = reshape(state, :)
    end
    obs, state
end

max_stars_repo_name: venuur/TSExplore

max_issues_repo_name: venuur/TSExplore

max_forks_repo_name: venuur/TSExplore


text:
# This file was generated, do not modify it. # hide
dfs = [wind, temp, raddir, raddif]
col_mean = [:windspeed_mean, :temp_mean, :raddir_mean, :raddif_mean];


text:
using Flux, DiffEqFlux, OrdinaryDiffEq, DiffEqParamEstim, Plots, Optim, Dates

using DiffEqBase, Flux, DiffResults, DiffEqSensitivity, ForwardDiff
using Flux.Tracker: @grad
using DiffEqSensitivity: adjoint_sensitivities_u0
######################################
mutable struct saver
    losses::Array{Float64,1}
    times::Array{Dates.Time,1}
    count_epochs::Int128
end
function saver(n_epochs)
    losses = zeros(n_epochs)
    times = fill(Dates.Time(Dates.now()),n_epochs)
    count_epochs = 0
    return saver(losses,times,count_epochs)
end
function update_saver(saver, loss_i, time_i)
    epoch_i = saver.count_epochs
    saver.losses[epoch_i] = loss_i
    saver.times[epoch_i] = time_i
end
####################################################### Observation ###############################################################
# Start conditions for the two species in the system
u0 = Float32[2.; 0.]
# Number of evaluations of the neural ODE. It relates to the numbers of layers of the neural net (depth of network).
datasize = 30
# Time span in which of evaluation will be and actual timepoints of evaluations
tspan = (0.0f0, 1.5f0)
t = range(tspan[1], tspan[2], length = datasize)
# The true ODE (with the true parameters) which the neural net should learn
function trueODEfunc(du, u, p, t)
  true_A = [-0.1 2.0; -2.0 -0.1]
  du .= ((u.^3)'true_A)'
end
# Construction of the ODEProblem and solving the ODEProblem with Tsit5 solver
prob = ODEProblem(trueODEfunc, u0, tspan)
ode_data = Array(solve(prob,Tsit5(),saveat=t))
################################################### Building a neural ODE ##########################################################
# Derivative is modeled by a neural net. Chain concatinates the functions ode function and two dense layers.
dudt = Chain(x -> x.^3,
       Dense(2,50,tanh),
       Dense(50,2))
# Parameters of the model which are to be learnt. They are: W1 (2x50), b1 (50), W2 (50x2), b2 (2)
ps = Flux.params(dudt)
# Build node
n_ode = x->neural_ode(dudt, x, tspan, Tsit5(), saveat=t, reltol=1e-7, abstol=1e-9)
# To set by user
init_train, number_sets, number_reps, number_epochs_loss1, number_epochs_loss2 = 5000, 1, 1, 1, 1000
# For saving
sa = saver(init_train+number_sets*(number_reps*(number_epochs_loss1+number_epochs_loss2)+init_train))
# Define collocation loss function with callback
function node_two_stage_function(model, x, tspan, saveat, ode_data,
            args...; kwargs...)
  dudt_(du,u,p,t) = du .= model(u)
  prob_fly = ODEProblem(dudt_,x,tspan)
  two_stage_method(prob_fly, saveat, ode_data)
end
loss_n_ode = node_two_stage_function(dudt, u0, tspan, t, ode_data, Tsit5(), reltol=1e-7, abstol=1e-9)
collocation_loss_fct()=loss_n_ode.cost_function(ps)
cb_collocation = function ()
    sa.count_epochs = sa.count_epochs +  1
    update_saver(sa, Tracker.data(collocation_loss_fct()),Dates.Time(Dates.now()))
    #println("\"",Tracker.data(collocation_loss_fct()),"\" \"",Dates.Time(Dates.now()),"\";")
end
# Define L2 loss function with callback
l2_loss_fct() = sum(abs2,ode_data .- n_ode(u0))
cb_l2 = function ()
    sa.count_epochs = sa.count_epochs +  1
    update_saver(sa, Tracker.data(l2_loss_fct()),Dates.Time(Dates.now()))
    #println("\"",Tracker.data(l2_loss_fct()),"\" \"",Dates.Time(Dates.now()),"\";")
end
# Training
data_iter_init_collocation = Iterators.repeated((), init_train)
data_iter_rep_collocation = Iterators.repeated((), number_epochs_loss1)
data_iter_rep_l2 = Iterators.repeated((), number_epochs_loss2)
opt = ADAM(0.1)
opt_l2 = ADAM(0.1)

# init train
Flux.train!(collocation_loss_fct, ps, data_iter_init_collocation, opt, cb = cb_collocation)
for j in 1:number_sets
    # Push in new area
    #for i in 1:number_reps
        Flux.train!(l2_loss_fct, ps, data_iter_rep_l2, opt_l2, cb = cb_l2)
        print("sa.count_epochs ", sa.count_epochs)
        print("\n")
        #Flux.train!(collocation_loss_fct, ps, data_iter_rep_collocation, opt, cb = cb_collocation)
    #end
    # Find local minimum
    #Flux.train!(collocation_loss_fct, ps, data_iter_init_collocation, opt, cb = cb_collocation)
end

# Call n_ode to get first prediction and to show startpoint for training.
pred = n_ode(u0)
scatter(t, ode_data[1,:], label="data")
    scatter!(t, Flux.data(pred[1,:]), label="prediction")
    scatter!(t, ode_data[2,:], label="data")
    scatter!(t, Flux.data(pred[2,:]), label="prediction")

#savefig("sogood.png")
using BSON: @save
name = string(init_train)*"_"*string(number_sets)*"_"*string(number_reps)*"_"*string(number_epochs_loss1)*"_"*string(number_epochs_loss2)
#@save "DiffEqFlux.jl/test/results/mix/model_mix_"*name*"_epochs.bson" dudt


text:
export isiso31661alpha3

validISO31661Alpha3CountriesCodes = Set([
    "AFG", "ALA", "ALB", "DZA", "ASM", "AND", "AGO", "AIA", "ATA", "ATG", "ARG", "ARM", "ABW", "AUS", "AUT", "AZE",
    "BHS", "BHR", "BGD", "BRB", "BLR", "BEL", "BLZ", "BEN", "BMU", "BTN", "BOL", "BES", "BIH", "BWA", "BVT", "BRA",
    "IOT", "BRN", "BGR", "BFA", "BDI", "KHM", "CMR", "CAN", "CPV", "CYM", "CAF", "TCD", "CHL", "CHN", "CXR", "CCK",
    "COL", "COM", "COG", "COD", "COK", "CRI", "CIV", "HRV", "CUB", "CUW", "CYP", "CZE", "DNK", "DJI", "DMA", "DOM",
    "ECU", "EGY", "SLV", "GNQ", "ERI", "EST", "ETH", "FLK", "FRO", "FJI", "FIN", "FRA", "GUF", "PYF", "ATF", "GAB",
    "GMB", "GEO", "DEU", "GHA", "GIB", "GRC", "GRL", "GRD", "GLP", "GUM", "GTM", "GGY", "GIN", "GNB", "GUY", "HTI",
    "HMD", "VAT", "HND", "HKG", "HUN", "ISL", "IND", "IDN", "IRN", "IRQ", "IRL", "IMN", "ISR", "ITA", "JAM", "JPN",
    "JEY", "JOR", "KAZ", "KEN", "KIR", "PRK", "KOR", "KWT", "KGZ", "LAO", "LVA", "LBN", "LSO", "LBR", "LBY", "LIE",
    "LTU", "LUX", "MAC", "MKD", "MDG", "MWI", "MYS", "MDV", "MLI", "MLT", "MHL", "MTQ", "MRT", "MUS", "MYT", "MEX",
    "FSM", "MDA", "MCO", "MNG", "MNE", "MSR", "MAR", "MOZ", "MMR", "NAM", "NRU", "NPL", "NLD", "NCL", "NZL", "NIC",
    "NER", "NGA", "NIU", "NFK", "MNP", "NOR", "OMN", "PAK", "PLW", "PSE", "PAN", "PNG", "PRY", "PER", "PHL", "PCN",
    "POL", "PRT", "PRI", "QAT", "REU", "ROU", "RUS", "RWA", "BLM", "SHN", "KNA", "LCA", "MAF", "SPM", "VCT", "WSM",
    "SMR", "STP", "SAU", "SEN", "SRB", "SYC", "SLE", "SGP", "SXM", "SVK", "SVN", "SLB", "SOM", "ZAF", "SGS", "SSD",
    "ESP", "LKA", "SDN", "SUR", "SJM", "SWZ", "SWE", "CHE", "SYR", "TWN", "TJK", "TZA", "THA", "TLS", "TGO", "TKL",
    "TON", "TTO", "TUN", "TUR", "TKM", "TCA", "TUV", "UGA", "UKR", "ARE", "GBR", "USA", "UMI", "URY", "UZB", "VUT",
    "VEN", "VNM", "VGB", "VIR", "WLF", "ESH", "YEM", "ZMB", "ZWE",
])

function isiso31661alpha3(str::AbstractString)::Bool
    return uppercase(str) in validISO31661Alpha3CountriesCodes
end


text:
@testset "vector1" begin
    env = SymbolicEnv()
    cache = SymbolicCache()
    x = [@expr(x), 1, 1.0]
    println(x)
end

@testset "vector2" begin
    env = SymbolicEnv()
    cache = SymbolicCache()
    x = [@expr(x), 1, 1.0]
    println(convert(AbstractSymbolic{Float64}, x))
end

@testset "SymbolicVector1" begin
    @bind x = 10
    v = [x + i for i = 1:10]
    @test seval(v) == [seval(x)+i for i = 1:10]
end

@testset "SymbolicVector1" begin
    @bind x = 10
    v = [x^i + i for i = 1:10]
    @test seval(v, :x) == [i * seval(x)^(i-1) for i = 1:10]
end

@testset "SymbolicVec3" begin
    m = @expr [1, x, 3]
    @test m[1].val == 1
    @test m[2].var == :x
    @test m[3].val == 3
end

@testset "SymbolicVec4" begin
    test = SymbolicEnv()
    m = @expr [1, x, y]
    @bind test x = 10
    @bind test y = 20
    result = seval(m[1:2], test)
    @test result == [1, 10]
end

@testset "SymbolicVec5" begin
    test = SymbolicEnv()
    m = @expr [1, x, y]
    expr = sum(m)
    println(expr)
    @bind test x = 10
    @bind test y = 20
    result = seval(expr, test)
    @test result == 31.0
    result = seval(expr, :x, test)
    @test result == 1.0
end

@testset "SymbolicDot1" begin
    test = SymbolicEnv()
    e1 = @expr [x^2, y, 10]
    e2 = @expr [x^2, y, 10]
    expr = dot(e1, e2)
    println(expr)
    x = 0.5
    y = 0.8
    @bind test x = x
    @bind test y = y
    @test seval(expr, test, SymbolicCache()) == dot([seval(x,test)^2, seval(y,test), 10], [seval(x,test)^2, seval(y,test), 10])
end

@testset "SymbolicDot2" begin
    test = SymbolicEnv()
    e1 = @expr [x^2, y, 10]
    e2 = @expr [x^2, y, 10]
    expr = dot(e1, e2)
    x0 = 0.5
    y0 = 0.8
    @bind test x = x0
    @bind test y = y0
    h = 0.0001
    test1 = SymbolicEnv()
    @bind test1 begin
        x = x0 + h
        y = y0
    end
    test2 = SymbolicEnv()
    @bind test2 begin
        x = x0 - h
        y = y0
    end
    ex = (seval(expr, test1, SymbolicCache()) - seval(expr, test2, SymbolicCache())) / (2*h)
    @test isapprox(seval(expr, :x, test, SymbolicCache()), ex, atol=1.0e-5)
end

@testset "SymbolicDot3" begin
    test = SymbolicEnv()
    e1 = @expr [x^2, y, 10]
    e2 = @expr [x^2, y, 10]
    expr = dot(e1, e2)
    x0 = 10.0
    y0 = 0.8
    @bind test begin
        x = x0
        y = y0
    end
    # h = 0.00001
    # @env test1 begin
    #     x = x+h
    #     y = y
    # end
    # @env test2 begin
    #     x = x-h
    #     y = y
    # end
    # ex = (seval(expr, test1, SymbolicCache()) - 2*seval(expr, test, SymbolicCache()) + seval(expr, test2, SymbolicCache())) / (h^2)
    # @test isapprox(ex, 4*3*x^2)
    @test isapprox(seval(expr, (:x,:x), test), 4*3*x0^2)
end

@testset "SymbolicPlus1" begin
    test = SymbolicEnv()
    e1 = convert(AbstractVectorSymbolic{Float64}, @expr [x^2, y, 10])
    e2 = convert(AbstractVectorSymbolic{Float64}, @expr [x, y^2, 10])
    expr = e1 + e2
    println(expr)
    x0 = 10.0
    y0 = 0.8
    @bind test begin
        x = x0
        y = y0
    end
    println(seval(expr, test))
    println(seval(expr[1], test))
    println(seval(expr[1:2], test))
    println(seval(expr, :x, test))
    println(seval(expr[1], :x, test))
    println(seval(expr[1:2], :x, test))
    println(seval(expr, (:x,:y), test))
    println(seval(expr[1], (:x,:y), test))
    println(seval(expr[1:2], (:x,:y), test))
    # @test isapprox(seval(expr, (:x,:x), test), 4*3*x0^2)
end

@testset "SymbolicDot5" begin
    @bind begin
        x = 4
        y = 3
    end
    v1 = [x, y]
    v2 = [x, 8]
    println(seval(dot(v1, v2)))
    println(seval(dot(v1, v2), :x))
    println(seval(dot(v1, v2), :y))
    println(seval(dot(v1, v2), (:x,:y)))
end

@testset "SymbolicPlus2" begin
    test = SymbolicEnv()
    e1 = convert(AbstractVectorSymbolic{Float64}, @expr [x^2, y, 10])
    e2 = convert(AbstractVectorSymbolic{Float64}, @expr [x, y^2, 10])
    expr = e1 + e2
    println(expr)
    x0 = 10.0
    y0 = 0.8
    @bind test begin
        x = x0
        y = y0
    end
    println(seval(expr, test))
    println(seval(expr[1], test))
    println(seval(expr[[true, false, false]], test))
end


text:
@testset "Comparisons" begin
    @testset "Equality" begin
        function check_seq_kmer_equality(len)
            a = DNAKmer(random_dna_kmer(len))
            b = LongDNA{4}(a)
            c = LongDNA{2}(a)
            return a == b == c && c == b == a
        end

        for len in [1, 10, 32, 64, 128]
            @test all(Bool[check_seq_kmer_equality(len) for _ in 1:reps])
        end

        # True negatives
        @test DNAKmer("ACG") != RNAKmer("ACG")
        @test DNAKmer("T")   != RNAKmer("U")
        @test DNAKmer("AC")  != DNAKmer("AG")
        @test RNAKmer("AC")  != RNAKmer("AG")
        @test AAKmer("MV")   != AAKmer("NM")

        @test DNAKmer("ACG") != rna"ACG"
        @test DNAKmer("T")   != rna"U"
        @test DNAKmer("AC")  != dna"AG"
        @test RNAKmer("AC")  != rna"AG"
        @test AAKmer("MV")   != aa"NM"

        @test rna"ACG" != DNAKmer("ACG")
        @test rna"U"   != DNAKmer("T")
        @test dna"AG"  != DNAKmer("AC")
        @test rna"AG"  != RNAKmer("AC")
        @test aa"MV"   != AAKmer("NM")
    end

    @testset "Inequality" begin
        for len in [1, 10, 32, 64]
            if len <= 32
                @test  isless(DNAKmer{1}((UInt64(0),)), DNAKmer{1}((UInt64(1),)))
                @test !isless(DNAKmer{1}((UInt64(0),)), DNAKmer{1}((UInt64(0),)))
                @test !isless(DNAKmer{1}((UInt64(1),)), DNAKmer{1}((UInt64(0),)))

                @test  isless(RNAKmer{1}((UInt64(0),)), RNAKmer{1}((UInt64(1),)))
                @test !isless(RNAKmer{1}((UInt64(0),)), RNAKmer{1}((UInt64(0),)))
                @test !isless(RNAKmer{1}((UInt64(1),)), RNAKmer{1}((UInt64(0),)))
            end
        end
    end

    @testset "Hash" begin
        kmers = map(DNAKmer, ["AAAA", "AACT", "ACGT", "TGCA"])
        for x in kmers, y in kmers
            @test (x == y) == (hash(x) == hash(y))
        end

        kmers = map(RNAKmer, ["AAAA", "AACU", "ACGU", "UGCA"])
        for x in kmers, y in kmers
            @test (x == y) == (hash(x) == hash(y))
        end

        kmers = map(AAKmer, ["AMVK", "FPST", "QEGH", "ARND"])
        for x in kmers, y in kmers
            @test (x == y) == (hash(x) == hash(y))
        end
    end
end


text:
#
# Do work here ...
#

v = 1 + 1
println("Julia side: I will return an Int64: ", v)

return convert(Int64, v)

max_stars_repo_name: yundantianchang/dynload-julia

max_issues_repo_name: yundantianchang/dynload-julia

max_forks_repo_name: yundantianchang/dynload-julia


text:
using MPI

function do_hello()
    comm = MPI.COMM_WORLD
    rk = MPI.Comm_rank(comm)
    sz = (MPI.Comm_size(comm))
    println("Hello world, I'm $(rk) of $(sz)")
    MPI.Barrier(comm)
end

function main()
    MPI.Init()
    do_hello()
    MPI.Finalize()
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
module Space

using Omega.Space
using Random
using Test
export test_resolve_repeatable, test_resolve_uniform

const TEST_SEED = 12345

function randvec(ω::Ω{I}, startIdx, nextIdx, n, T::Type) where {I}
  id = startIdx
  vec = []
  for i in 1:n
    push!(vec, memrand(ω, id, T))
    id = nextIdx(id)
  end
  vec
end

function test_resolve_repeatable(ω::Ω{I}, startIdx, nextIdx) where {I}
  Random.seed!(TEST_SEED)
  N = 10000
  vec1 = randvec(ω, startIdx, nextIdx, N, Int)
  vec2 = randvec(ω, startIdx, nextIdx, N, Int)
  vec1 == vec2
end

function close_to(num, targ, sigma)
  abs(num-targ) < 6 * sigma
end

# Extra precondition: ω must be able to take Bool values
function test_resolve_uniform(ω::Ω{I}, startIdx, nextIdx) where {I}
  Random.seed!(TEST_SEED)
  N = 10000
  vec = randvec(ω, startIdx, nextIdx, N, Bool)
  numTrue = length(filter(x -> x, vec))
  bernVariance = sqrt(0.5*(1-0.5))
  close_to(numTrue / N, 0.5, bernVariance / sqrt(N))
end

end


text:
# Generate an array of thetas (useful for custom parallelization)
function generate_thetas(plan::abc_pmc_plan_type, sampler::Distribution, ss_true, epsilon::Float64, num_to_generate::Integer; num_max_attempt = plan.num_max_attempt)
    @assert(epsilon>0.0)
    ex_sample = rand(sampler)
    num_param = length(ex_sample)
    thetas = Array(eltype(ex_sample),num_param,num_to_generate)
    dists = Array(Float64, num_to_generate)  # Warning should make aware of distance return type
    attempts = Array(Int64,  num_to_generate)
    for i in 1:num_to_generate
      thetas[:,i], dists[i], attempts[i] = generate_theta(plan, sampler, ss_true, epsilon, num_max_attempt = num_max_attempt )
    end
    return (thetas, dists, attempts)
end

# Generate initial abc population from prior, aiming for d(ss,ss_true)<epsilon
function init_abc_parallel_custom(plan::abc_pmc_plan_type, ss_true)   # WARNING PARALLEL CODE UNTESTED
  num_param = length(Distributions.rand(plan.prior))
  num_draw = plan.num_part
  nw = nworkers()

  job = Array(RemoteRef, nw)
  a = [1:ifloor(num_draw/nw):num_draw][1:nw]
  b = Array(Int64,nw)
  for i in 1:nw-1
     b[i] = a[i+1]-1
  end
  b[end] = num_draw

  for j in 1:nw
    #println("# Starting jobs ",j," with entries ", a[j], " to ", b[j])
    # Draw initial set of theta's from prior (either with dist<epsilon or best of num_max_attempts)
    job[j] = @spawn generate_thetas(plan, plan.prior, ss_true, plan.epsilon_init, b[j]-a[j]+1)
  end

  # Allocate arrays for results on master
  theta = Array(Float64,(num_param,plan.num_part))
  dist_theta = Array(Float64, plan.num_part)
  attempts = Array(Int64, plan.num_part)
  # Retreive results
  for j in 1:nw
    theta[:,a[j]:b[j]], dist_theta[a[j]:b[j]], attempts[a[j]:b[j]] = fetch(job[j])
  end
  weights = fill(1.0/plan.num_part,plan.num_part)
  return abc_population_type(theta,weights,dist_theta)
end


text:
import ..SamplingConstraints: constrain
import ..UncertainDatasets:
    UncertainIndexValueDataset,
    UncertainIndexDataset,
    ConstrainedUncertainIndexDataset,
    ConstrainedUncertainValueDataset,
    AbstractUncertainValueDataset,
    AbstractUncertainIndexValueDataset

function constrain(x::UncertainIndexDataset,
        resampling::ConstrainedValueResampling{1})

    ConstrainedUncertainIndexDataset(constrain(x.indices, resampling.constraints[1]))
end

function constrain(x::AbstractUncertainValueDataset,
        resampling::ConstrainedValueResampling{1})

    ConstrainedUncertainValueDataset(constrain(x.values, resampling.constraints[1]))
end

function constrain(x::AbstractUncertainIndexValueDataset,
        resampling_inds::ConstrainedValueResampling{N1},
        resampling_vals::ConstrainedValueResampling{N2}) where {N1, N2}

    inds = constrain(x.indices, resampling_inds)
    vals = constrain(x.values, resampling_vals)

    UncertainIndexValueDataset(inds, vals)
end

function constrain(x::AbstractUncertainIndexValueDataset, resampling::ConstrainedIndexValueResampling{2, 1})

    inds = constrain(x.indices, resampling[1, 1])
    vals = constrain(x.values, resampling[1, 2])

    UncertainIndexValueDataset(inds, vals)
end

function constrain(x::AbstractUncertainIndexValueDataset,
        constraints_inds::CONSTRAINT_TYPES,
        constraints_vals::CONSTRAINT_TYPES)

    inds = constrain(x.indices, constraints_inds)
    vals = constrain(x.values, constraints_vals)

    UncertainIndexValueDataset(inds, vals)
end

export constrain


text:
using GPlot; include("_test.jl")

include("utils.jl")

include("figure.jl")
include("ax.jl")
include("ax_elem.jl")
include("legend.jl")

include("drawing.jl")
include("drawing2.jl")

include("style.jl")

include("layout.jl")
include("render.jl")

include("properties.jl")

include("object.jl")

continuous_preview(_bk_cont_preview)
G.GP_ENV["PALETTE"] = _bk_palette
G.GP_ENV["SIZE_PALETTE"] = _bk_palette
true


text:
struct SType

end


text:
function _AudzeEglaisObjective!(dim::Continuous,dist,LHC)
    n,d = size(LHC)
    dist .= 0.0

    # Squared l-2 norm of distances between all (unique) points
    l = 0
    for i = 2:n
        for j = 1:i-1
            l += 1
            for k = 1:d
                dist[l] += (LHC[i,k]-LHC[j,k])^2
            end
            dist[l] = 1/dist[l]
        end
    end
    output = 1/(sum(dist))
end

function _AudzeEglaisObjective!(dim::Categorical,dist,LHC)
    n,d = size(LHC)
    dist .= 0.0

    # Squared l-2 norm of distances between all points of the same category
    l = 0

    for i = 2:n
        for j = 1:i-1
            l += 1
            for k = 1:d
                dist[l] += (LHC[i,k]-LHC[j,k])^2
            end
            dist[l] = 1/dist[l]
        end
    end
    output = 1/(sum(dist))
    if output == Inf
        output = 0
    end
    output
end

"""
    function AudzeEglaisObjective!(dist,LHC::T) where T <: AbstractArray
Return the scalar which should be maximized when using the Audze-Eglais
distance as the objective function. Note this is the inverse of the typical
Audze-Eglais distance which normally is minimized.
"""
function AudzeEglaisObjective!(dist,LHC::T; dims::Array{V,1} =[Continuous() for i in 1:size(LHC,2)],
                                            interSampleWeight::Float64=1.0,
                                            ) where T <: AbstractArray where V <: LHCDimension

    out = 0.0

    #Compute the objective function among all points
    out += _AudzeEglaisObjective!(Continuous(),dist,LHC)*interSampleWeight

    categoricalDimInds = findall(x->typeof(x)==Categorical,dims)
    for i in categoricalDimInds
        for j = 1:dims[i].levels
            subLHC = @view LHC[LHC[:,i] .== j,:]
            out += _AudzeEglaisObjective!(dims[i],dist,subLHC)*dims[i].weight
        end
    end

    return out
end

"""
    function AudzeEglaisObjective(LHC::T) where T <: AbstractArray
Same as AudzeEglaisObjective!(dist,LHC::Array) but creating a new distance array.
"""
function AudzeEglaisObjective(LHC::T;   dims::Array{V,1}=[Continuous() for i in 1:size(LHC,2)],
                                        interSampleWeight::Float64=1.0,
                                        ) where T <: AbstractArray where V <: LHCDimension

    n = size(LHC,1)
    dist = zeros(Float64,Int(n*(n-1)*0.5))

    return AudzeEglaisObjective!(dist,LHC; dims=dims, interSampleWeight=interSampleWeight)

end


text:
#
# Structure to store auxiliary arrays
#
struct VAux{T,F}
    g::T
    xn::T
    gn::T
    fprev::Vector{F}
end
function VAux(x,fx;m=10)
    g = similar(x)
    xn = similar(x)
    gn = similar(x)
    fprev = typeof(fx)[]
    for i in 1:m
        push!(fprev,zero(fx))
    end
    return VAux{typeof(g),typeof(fx)}(g, xn, gn, fprev)
end


text:
#
#      MatrixManifold -- a matrix manifold,
#       i.e. the values/points on the manifold are matrices
#		this is here implemented as traits, such that every manifold can
#		“gain” this feature by setting @traitimlp IsMatrixM{}
#		for the manifold, its point and the tangent vector
#
# Manopt.jl, R. Bergmann, 2018-06-26
import LinearAlgebra: transpose
import Base: +, -, *, /
export transpose, IsMatrixM, IsMatrixP, IsMatrixTV
export +,-,*,/
"""
    IsMatrixM{X}
An abstract Manifold to represent a manifold whose points are matrices.
For these manifolds the usual operators (+,-,*) are overloaded for points.
Furthermore, the `transpose` is also overloaded, though it returns the matrix,
since the dimensions mit be different for rectangular matrices.
"""
@traitdef IsMatrixM{X}
"""
    IsMatrixP{X}
An abstract Manifold Point belonging to a matrix manifold.
"""
@traitdef IsMatrixP{X}
"""
    IsMatrixTV{X}
An abstract Manifold Point belonging to a matrix manifold.
"""
@traitdef IsMatrixTV{X}

# for all that satisfy IsMatrixM -> introduce operators on points and points/TVecs
@traitfn +(x::P, y::P) where {P <: MPoint; IsMatrixP{P}}  = P( getValue(x) + getValue(y) )
@traitfn -(x::P, y::P) where {P <: MPoint; IsMatrixP{P}}  = P( getValue(x) - getValue(y) )
@traitfn *(x::P, y::P) where {P <: MPoint; IsMatrixP{P}}  = P( getValue(x) * getValue(y) )
@traitfn transpose(x::P) where {P <: MPoint; IsMatrixP{P}} = P( Matrix(transpose( getValue(x) ) ) )
@traitfn *(ξ::T, ν::T) where {T <: TVector; IsMatrixTV{T}} = T( getValue(ξ) * getValue(ν) )
@traitfn transpose(ξ::T) where {T <: TVector; IsMatrixTV{T}} = T( Matrix( transpose( getValue(ξ) ) ) )


text:
module FreeCell
using PlayingCards52, ImplicitGraphs

import Base: show, string, length, isempty, hash, (==)

# A FreeCell tableau consists of:
# * four free cells
# * four foundations
# * eight cascades

ace_high(false)

function row_string(list::Vector{Card})::String
    if length(list) == 0
        return "\t--"
    end
    result = "\t"
    for c in list
        result *= string(c)
        if c != list[end]
            result *= " "
        end
    end
    return result
end

export score

include("free_cells.jl")
include("foundation.jl")
include("cascade.jl")
include("tableau.jl")
include("mover.jl")
include("game_graph.jl")
include("game_solver.jl")

export new_game
"""
`new_game()` creates a new, random instance of a FreeCell game.
"""
new_game() = Tableau(deck())

end # module


text:
function check_registry()

    basedir = Registry.environment_path
    mljmodelsdir = joinpath(basedir, "..", "..", ".")
    Pkg.activate(basedir)
    Pkg.develop(PackageSpec(path=mljmodelsdir))
    Pkg.instantiate()
    Pkg.precompile()

    # Read Metadata.toml
    dict = TOML.parsefile(joinpath(basedir, "Metadata.toml"))

    problems = String[]
    for (package, model_dict) in dict
        for (model, meta) in model_dict
            # check if new entry or changed entry, otherwise don't test
            key = "$package.$model"
            program = quote
                @load $model pkg=$package verbosity=-1
            end
            try
                eval(program)
                # add/refresh entry
                print(rpad("Entry for $key was loaded properly ✓", 79)*"\r")
            catch ex
                push!(problems, string(key))
                @error "⚠ there was an issue trying to load $key" exception=ex
            end
        end
    end
    return problems
end


text:
# NOTE: we cannot change subblocks of a primitive block
#      since they are primitive, therefore we return themselves
chsubblocks(x::PrimitiveBlock, it) = x
subblocks(x::PrimitiveBlock) = ()
# NOTE: all primitive block should name with postfix Gate
#       and each primitive block should stay in a single
#       file whose name is in lowercase and underscore.
include("const_gate.jl")
include("identity_gate.jl")
include("phase_gate.jl")
include("shift_gate.jl")
include("rotation_gate.jl")
include("time_evolution.jl")
include("general_matrix_gate.jl")
include("measure.jl")


text:
using Roots, NLopt, PyPlot, LsqFit, JLD, Polynomials, JFVM, Dierckx
include("DME_solver.jl")

function rec_fact(param, G, K, MT, T, BC, IC, k_w, t_exp, R_exp)
  k_o = param[1]
  (t, R) = DME_solver(G, K, MT, T, BC, IC, k_o, k_w, false)
  R_int = Spline1D(t, R, k=1)
  return R_int(t_exp)
end

model = (t_exp, param) -> rec_fact(param, G, K, MT, T, BC, IC, k_w, t_exp, R_exp)

function error_calc(param, G, K, MT, T, BC, IC, k_w, t_exp, R_exp)
  return rec_fact(param, G, K, MT, T, BC, IC, k_w, t_exp, R_exp) - R_exp
end

function error_calculation(param, param_ind, G, K, MT, T, BC, IC, k_w, t_exp, R_exp; grad_calc=true)
  eps1 = 1e-8
  error_vals = error_calc(param, G, K, MT, T, BC, IC, k_w, t_exp, R_exp)
  if grad_calc
    jac  = zeros(length(error_vals), length(param_ind))
  else
    jac = zeros(0)
  end
  if grad_calc
    for j in eachindex(param_ind)
      param2 = copy(param)
      param2[param_ind[j]]+=eps1
      error_val2 = error_calc(param2, G, K, MT, T, BC, IC, k_w, t_exp, R_exp)
      jac[:, j] = (error_val2 - error_vals) / eps1
    end
  end
  return error_vals, jac
end

function objective_function(param::Vector{Float64}, grad::Vector{Float64}, param_all, param_ind, G, K, MT, T, BC, IC, k_w, t_exp, R_exp; w=1.0)
  param2 = copy(param_all)
  param2[param_ind] = param
  grad_calc = false
  if length(grad)>0
    grad_calc = true
  end
  error_val, jac = error_calculation(param2, param_ind, G, K, MT, T, BC, IC, k_w, t_exp, R_exp, grad_calc = grad_calc)
  obj_func_value = sum(abs2, w.*error_val)
  if length(grad)>0
      grad[:] = 2.0*sum(w.*error_val.*jac, 1)
  end
  return obj_func_value
end

obj_fun = (param, grad) -> objective_function(param::Vector{Float64}, grad::Vector{Float64}, param_all, param_ind, G, K, MT, T, BC, IC, k_w, t_sec_cor, R_oil, w=w)

  perm_val  = 0.01e-12 # [m^2] permeability
  poros_val = 0.40     # [-] porosity
  ρ_rock    = 2700     # [kg/m^3]

  u_inj = 1.0/(3600*24) # 1 m/day to m/s injection velocity
  c_inj1 = 0.0 # mol/m^3 DME in injection water at initial water flooding
  c_inj2 = 10.0 # mol/m^3 DME in injection water at DME flooding
  c_tracer1 = 0.0
  c_tracer2 = 1.0

  Kval = 2.0
  a = 2.0*2700.0*1000.0*(1.0-poros_val) # m^2/m^3 = 2 m^2/g * 2700 kg/m^3 * 1000 g/kg * (1-phi)
  D_oil = 0.5*10^(-9.0) # m^2/s
  D_water = 2*10^(-9.0) # m^2/s

  ## Creating DME parameters
  ρ_oil = c_o -> 748.0015057043045-0.0016431963703508173.*c_o+2.9364412314092913.*10.0.^(-7.0).*c_o.^2.0-4.107225978803957.*10.0.^(-11.0).*c_o.^3.0
  dρ_oil = c_o -> -0.0016431963703508173+2.0.*2.9364412314092913.*10.0.^(-7.0).*c_o-3.0.*4.107225978803957.*10.0.^(-11.0).*c_o.^2.0
  ρ_water = 983.6319389796079
  μ_oil = c_o -> 0.0010131130815849523-2.528969514668829.*10.0.^(-5.0).*c_o+2.3873417406201346.*10.0.^(-11.0).*c_o.^2.0
  dμ_oil = c_o -> -2.528969514668829.*10.0.^(-5.0)+2*2.3873417406201346.*10.0.^(-11.0).*c_o
  μ_water = c_w -> 0.00040706344808758737+1.47269029810352587.*10.0.^(-7.0).*c_w
  dμ_water = c_w -> 1.47269029810352587.*10.0.^(-7.0).*ones(size(c_w))

  MT = FluidParametersDME(ρ_oil, dρ_oil, ρ_water, μ_oil, dμ_oil, μ_water, dμ_water, Kval, a, D_oil, D_water)

  # Creating realitive permeability
  kro0 = 0.8
  krw0 = 0.2
  n_o  = 2.0
  n_w  = 2.0
  swc  = 0.08

  sor_min  = 0.05
  sor_max  = 0.3
  c_min = 0.0 # min DME mol/m^3 in oil
  c_max = Kval*c_inj2 # max DME mol/m^3 fraction in oil
  m_sor = -(sor_max-sor_min)/(c_max-c_min)
  sor_c = c -> m_sor*c+sor_max

  K = RelativePermeabilityDME(kro0, krw0, n_o, n_w, sor_c, swc)

  # =============================================================================
  # Preparing geology/mesh
  # =============================================================================
  dim = 1

  srand(123456543)
  λ_dp    = 0.8
  cor_x   = 0.1
  cor_y   = 0.1
  cor_z   = 0.1
  Lx      = 80 # [m]
  Ly      = 20 # [m]
  Lz      = 10 # [m]
  Nx      = 50  # number of grids in the x direction
  Ny      = 20  # number of grids in the y direction
  Nz      = 10 # number of grids in the z direction

  if dim == 1
      m           = createMesh1D(Nx, Lx)
      perm_mat    = permfieldlogrnde(Nx, perm_val, λ_dp, cor_x)
  elseif dim == 2
      m           = createMesh2D(Nx, Ny, Lx, Ly)
      perm_mat    = permfieldlogrnde(Nx, Ny, perm_val, λ_dp, cor_x, cor_y)
  else
      m           = createMesh3D(Nx, Ny, Nz, Lx, Ly, Lz)
      perm_mat    = permfieldlogrnde(Nx, Ny, Nz, perm_val, λ_dp, cor_x, cor_y, cor_z)
  end

  k       = createCellVariable(m, perm_mat)
  k_face  = harmonicMean(k)     # permeability on the cell faces
  ϕ       = createCellVariable(m, poros_val)

  G = Geology(m, perm_val, poros_val, k, ϕ)

  # =============================================================================
  # Preparing time-stepping parameters
  # =============================================================================

  ## Simulation Time
  n_pv_water1       = 0.8 # number of injected pore volumes
  t_final_water1    = n_pv_water1*Lx/(u_inj/poros_val) # [s] final time for water flood
  dt_water1         = t_final_water1/n_pv_water1/Nx/5 # [s] time step

  n_pv_DME          = 1.0
  t_final_DME       = t_final_water1+n_pv_DME*Lx/(u_inj/poros_val) # [s] final time for DME flood
  dt_DME            = t_final_DME/(n_pv_water1 + n_pv_DME)/Nx/5 # [s] time step

  n_pv_water2       = 0.5
  t_final_water2    = t_final_DME+n_pv_water2*Lx/(u_inj/poros_val) # [s] final time for DME flood
  dt_water2         = t_final_water2/(n_pv_water1 + n_pv_DME + n_pv_water2)/Nx/5 # [s] time step

  T_water1 = SimulationTime(t_final_water1, dt_water1)
  T_DME    = SimulationTime(t_final_DME, dt_DME)
  T_water2 = SimulationTime(t_final_water2, dt_water2)

  T = [T_water1, T_DME, T_water2]

  ## Boundary Conditions for water flooding
  #T0 = 70 + 273.15 # [K]

  c_oil0      = 0.0 # initial DME concentration in the oil phase
  c_water0    = 0.0 # initial DME concentration in the water phase
  p0          = 2000/14.7*1e5   # [Pa]
  sw0         = swc+0.05 # [vol frac] initial water saturation
  c_t0        = 0.0 # initial tracer concentration

  BCc_oil_water     = createBC(m) # concentration (DME in oil) boundary condition
  BCc_water_water   = createBC(m) # concentration (DME in water) boundary condition
  BCp_water         = createBC(m) # pressure boundary condition
  BCs_water         = createBC(m) # saturation boundary condition
  BCt_water         = createBC(m) # tracer concentration boundary condition

  BCp_water.right.a[:] = 0.0
  BCp_water.right.b[:] = 1.0 #Dirichlet
  BCp_water.right.c[:] = p0

  BCp_water.left.a[:]  = k_face.xvalue[1,:]./μ_water(c_inj1) #Neumann
  BCp_water.left.b[:]  = 0.0
  BCp_water.left.c[:]  = -u_inj

  BCs_water.left.a[:]  = 0.0
  BCs_water.left.b[:]  = 1.0 #Dirichlet
  BCs_water.left.c[:]  = 1.0 #Value S_w=1

  BCc_water_water.left.a[:]  = 0.0 #
  BCc_water_water.left.b[:]  = 1.0 #Dirichlet
  BCc_water_water.left.c[:]  = c_inj1 #injected concentration of DME in water

  BCc_oil_water.left.a[:]  = 0.0 #
  BCc_oil_water.left.b[:]  = 1.0 #Dirichlet
  BCc_oil_water.left.c[:]  = Kval*c_inj1 #injected concentration of DME in water

  BCt_water.left.a[:]  = 0.0
  BCt_water.left.b[:]  = 1.0 #Dirichlet
  BCt_water.left.c[:]  = c_tracer1

  BC_water = BoundaryConditionDME(BCc_oil_water, BCc_water_water, BCp_water, BCs_water, BCt_water)

  c_oil_init      = createCellVariable(m, c_oil0, BCc_oil_water)
  c_water_init    = createCellVariable(m, c_water0, BCc_water_water)
  p_init          = createCellVariable(m, p0, BCp_water)
  sw_init         = createCellVariable(m, sw0, BCs_water)
  c_t_init        = createCellVariable(m, c_t0, BCt_water)

  IC = InitialConditionDME(c_oil_init, c_water_init, p_init, sw_init, c_t_init)

  # Creating Boundary conditions for DME flooding
  BCc_oil_DME     = createBC(m) # concentration (DME in oil) boundary condition
  BCc_water_DME   = createBC(m) # concentration (DME in water) boundary condition
  BCp_DME         = createBC(m) # pressure boundary condition
  BCs_DME         = createBC(m) # saturation boundary condition
  BCt_DME         = createBC(m) # tracer concentration boundary condition

  BCc_water_DME.left.a[:] = 0.0
  BCc_water_DME.left.b[:] = 1.0 #Dirichlet
  BCc_water_DME.left.c[:] = c_inj2 #Injected concentration of DME

  BCc_oil_DME.left.a[:] = 0.0
  BCc_oil_DME.left.b[:] = 1.0 #Dirichlet
  BCc_oil_DME.left.c[:] = Kval*c_inj2 #Created from reaction at boundary

  BCp_DME.right.a[:] = 0.0
  BCp_DME.right.b[:] = 1.0 #Dirichlet
  BCp_DME.right.c[:] = p0

  BCp_DME.left.a[:]  = k_face.xvalue[1,:]/μ_water(c_inj2) #Neumann
  BCp_DME.left.b[:]  = 0.0
  BCp_DME.left.c[:]  = -u_inj

  BCs_DME.left.a[:]  = 0.0
  BCs_DME.left.b[:]  = 1.0 #Dirichlet
  BCs_DME.left.c[:]  = 1.0 #Value S_w=1

  BCt_DME.left.a[:]       = 0.0
  BCt_DME.left.b[:]       = 1.0 #Dirichlet
  BCt_DME.left.c[:]       = c_tracer2

  BC_DME = BoundaryConditionDME(BCc_oil_DME, BCc_water_DME, BCp_DME, BCs_DME, BCt_DME)

  BC = [BC_water, BC_DME, BC_water]

## Experimental results
t_sec_cor = [0, 16.0/6, 16.0/6*2, 16.0/6*3, 16.0/6*4, 16.0/6*5,  16.0, 82.0, 100.0, 112.0].*24.0.*60.0.*60.0
R_exp = [0.0, 10.0/100, 20.0/100, 30.0/100, 40.0/100, 50.0/100, 58.0/100, 62.0/100, 80.0/100, 88.0/100]
R_oil = R_exp

# =============================================================================
# Calculating heuristic for k_o and k_w
# =============================================================================
d = sqrt(perm_val) #Characteristic length of capillary in porous medium
c_o = Kval*c_inj2
c_w = c_inj2
v_o = u_inj/poros_val
v_w = u_inj/poros_val

k_o_ini = 1.1*10.0^(-0.0)*0.023*D_oil/d*((ρ_oil(c_o)*v_o*d)/μ_oil(c_o))^(0.83)*(μ_oil(c_o)/(ρ_oil(c_o)*D_oil))^(1.0/3.0)
k_w = 10.0^(-0.0)*0.023*D_water/d*((ρ_water*v_w*d)/μ_water(c_w))^(0.83)*(μ_water(c_w)/(ρ_water*D_water))^(1.0/3.0)

param_all = [k_o_ini]
param_ind = [1]
w = ones(length(R_exp))

x_init = copy(param_all)
x_lb = [0.0]
x_ub = [1e-5]

# =============================================================================
# Optimizing
# =============================================================================

opt_type = 2

if opt_type == 1
    fit = curve_fit(model, t_sec_cor, R_exp, w, x_init, lower = x_lb, upper = x_ub)
    print("k_o_ini=", k_o_ini)
    print("k_o_opt=", fit.param[1])

elseif opt_type == 2
    # algorithms
    # :LD_MMA
    # :LN_COBYLA
    # :LD_LBFGS
    # :GN_DIRECT
    # :GN_DIRECT_L
    opt_alg=:LD_MMA

    opt1 = Opt(opt_alg, length(x_init)) # choose the algorithm
    lower_bounds!(opt1, x_lb)
    upper_bounds!(opt1, x_ub)
    ftol_rel!(opt1, 1e-5)
    ftol_abs!(opt1, 1e-5)

    min_objective!(opt1, obj_fun)
    (fObjOpt, paramOpt, flag) = optimize(opt1, x_init)
    print("k_o_ini=", k_o_ini, "\n\n")
    print("k_o_opt=", paramOpt[1])
end

max_stars_repo_name: FrederikLehn/coupled-solvers

max_issues_repo_name: FrederikLehn/coupled-solvers

max_forks_repo_name: FrederikLehn/coupled-solvers


text:
const LinsolveSubarray = SubArray{Float64,1,Vector{Float64},Tuple{UnitRange{Int64}},true}

function admm_step!(x::Vector{Float64},
	s::SplitVector{Float64},
	μ::Vector{Float64},
	ν::LinsolveSubarray,
	x_tl::LinsolveSubarray,
	s_tl::Vector{Float64},
	ls::Vector{Float64},
	sol::Vector{Float64},
	kkt_solver::AbstractKKTSolver,
	q::Vector{Float64},
	b::Vector{Float64},
	ρ::Vector{Float64},
	α::Float64,
	σ::Float64,
	m::Int64,
	n::Int64,
	set::CompositeConvexSet{Float64})
	# linear solve
	# Create right hand side for linear system
	# deconstructed solution vector is ls = [x_tl(n+1); ν(n+1)]
	# x_tl and ν are automatically updated, since they are views on sol
	@. ls[1:n] = σ * x - q
	@. ls[(n + 1):end] = b - s + μ / ρ
	s_time = @elapsed solve!(kkt_solver,sol,ls)

	# Over relaxattion
	@. x = α * x_tl + (1.0 - α) * x
	@. s_tl = s - (ν + μ) / ρ
	@. s_tl = α * s_tl + (1.0 - α) * s
	@. s = s_tl + μ / ρ

	# Project onto cone
	p_time = @elapsed project!(s, set)

	# update dual variable μ
	@. μ = μ + ρ .* (s_tl - s)
	return p_time, s_time
end

# SOLVER ROUTINE
# -------------------------------------

"""
optimize!(model)

Attempts to solve the optimization problem defined in `COSMO.Model` object with the user settings defined in `COSMO.Settings`. Returns a `COSMO.Result` object.
"""
function optimize!(ws::COSMO.Workspace)
	solver_time_start = time()
	settings = ws.settings
  # perform chordal decomposition
  if settings.decompose
  	ws.times.graph_time = @elapsed COSMO.chordal_decomposition!(ws)
  end
	# create scaling variables
	# with scaling    -> uses mutable diagonal scaling matrices
	# without scaling -> uses identity matrices
	ws.sm = (settings.scaling > 0) ? ScaleMatrices(ws.p.model_size[1], ws.p.model_size[2]) : ScaleMatrices()

	# perform preprocessing steps (scaling, initial KKT factorization)
	if isnothing(ws.kkt_solver)
		ws.sm = (settings.scaling > 0) ? ScaleMatrices(ws.p.model_size[1], ws.p.model_size[2]) : ScaleMatrices()
		# perform preprocessing steps (scaling, initial KKT factorization)
		ws.times.factor_time = 0
		ws.times.setup_time = @elapsed setup!(ws);
		ws.times.proj_time  = 0. #reset projection time
		ws.times.sol_time  = 0. #reset linear solve time
	else
		ws.times.setup_time = 0.
		ws.times.factor_time = 0
		ws.times.proj_time  = 0. #reset projection time
		ws.times.sol_time  = 0. #reset linear solve time
		reset_iteration_counters!(ws) # this sets lobpcg's iteration = 1, i.e. it is
		# only suitable for warm starting
	end

	# instantiate variables
	ws.iterations = 0
	status = :Unsolved
	cost = Inf
	r_prim = Inf
	r_dual = Inf

	# print information about settings to the screen
	settings.verbose && print_header(ws)
	time_limit_start = time()

	m, n = ws.p.model_size
	δx = zeros(n)
	δy = SplitVector(zeros(m), ws.p.C)

	s_tl = zeros(m) # i.e. sTilde

	ls = zeros(n + m)
	sol = zeros(n + m)
	x_tl = view(sol, 1:n) # i.e. xTilde
	ν = view(sol, (n + 1):(n + m))

	settings.verbose_timing && (iter_start = time())

	for iter = 1:settings.max_iter

		ws.iterations += 1

		@. δx = ws.vars.x
		@. δy = ws.vars.μ

		proj_time, sol_time = admm_step!(
			ws.vars.x, ws.vars.s, ws.vars.μ, ν,
			x_tl, s_tl, ls, sol,
			ws.kkt_solver, ws.p.q, ws.p.b, ws.ρvec,
			settings.alpha, settings.sigma,
			m, n, ws.p.C);

		ws.times.proj_time += proj_time
		ws.times.sol_time += sol_time
		# compute deltas for infeasibility detection
		@. δx = ws.vars.x - δx

		#@. δy = -ws.vars.μ + δy
		@. δy -= ws.vars.μ

		# compute residuals (based on optimality conditions of the problem) to check for termination condition
		# compute them every {settings.check_termination} step
		mod(iter, settings.check_termination)  == 0 && ((r_prim, r_dual) = calculate_residuals!(ws))

		# check convergence with residuals every {settings.checkIteration} steps
		if mod(iter, settings.check_termination) == 0
			# update cost
			cost = ws.sm.cinv[] * (1/2 * ws.vars.x' * ws.p.P * ws.vars.x + ws.p.q' * ws.vars.x)[1]

			if abs(cost) > 1e20
				status = :Unsolved
				break
			end

			# print iteration steps
			settings.verbose && print_iteration(ws, iter, cost, r_prim, r_dual, time() - solver_time_start)

			if has_converged(ws, r_prim, r_dual)
				status = :Solved
				break
			end
		end

		# check infeasibility conditions every {settings.checkInfeasibility} steps
		if mod(iter, settings.check_infeasibility) == 0
			if is_primal_infeasible!(δy, ws)
				status = :Primal_infeasible
				cost = Inf
				break
			end

			if is_dual_infeasible!(δx, ws)
				status = :Dual_infeasible
				cost = -Inf
				break
			end
		end

		# adapt rhoVec if enabled
		if settings.adaptive_rho && (mod(iter, settings.adaptive_rho_interval) == 0) && (settings.adaptive_rho_interval > 0)
			adapt_rho_vec!(ws)
		end

		if settings.time_limit !=0 &&  (time() - time_limit_start) > settings.time_limit
			status = :Time_limit_reached
			break
		end

	end #END-ADMM-MAIN-LOOP

	settings.verbose_timing && (ws.times.iter_time = (time() - iter_start))
	settings.verbose_timing && (ws.times.post_time = time())

	# calculate primal and dual residuals
	if ws.iterations == settings.max_iter
		r_prim, r_dual = calculate_residuals!(ws)
		status = :Max_iter_reached
	end

	# reverse scaling for scaled feasible cases
	if settings.scaling != 0
		reverse_scaling!(ws)
		# FIXME: Another cost calculation is not necessary since cost value is not affected by scaling
		cost =  (1/2 * ws.vars.x' * ws.p.P * ws.vars.x + ws.p.q' * ws.vars.x)[1] #sm.cinv * not necessary anymore since reverseScaling
	end

	# reverse chordal decomposition
	if settings.decompose
	 reverse_decomposition!(ws, settings)
	end

	ws.times.solver_time = time() - solver_time_start
	settings.verbose_timing && (ws.times.post_time = time() - ws.times.post_time)
	# print solution to screen
	settings.verbose && print_result(status, ws.iterations, cost, ws.times.solver_time)

	# create result object
	res_info = ResultInfo(r_prim, r_dual)

	y = -ws.vars.μ

	return Result{Float64}(ws.vars.x, y, ws.vars.s.data, cost, ws.iterations, status, res_info, ws.times);

end

function free_memory!(ws)
	free_memory!(ws.kkt_solver)
end

function reset_iteration_counters!(ws)
	reset_iteration_counters!(ws.kkt_solver)
	for set in ws.p.C.sets
		reset_iteration_counters!(set)
	end
end


text:
abstract type scatteringTransform{Dimension,Depth} end
struct stFlux{Dimension,Depth,ChainType,D,E,F} <: scatteringTransform{Dimension,Depth}
    mainChain::ChainType
    normalize::Bool
    outputSizes::D
    outputPool::E
    settings::F
end

import Base.ndims
ndims(s::scatteringTransform{D}) where {D} = D
nPathDims(ii) = 1 + max(min(ii - 2, 1), 0) # the number of path dimensions at layer ii (zeroth
# is ii=1)
depth(s::scatteringTransform{Dim,Depth}) where {Dim,Depth} = Depth
function Base.show(io::IO, st::stFlux{Dim,Dep}) where {Dim,Dep}
    layers = st.mainChain.layers
    σ = st.settings[:σ]
    Nd = ndims(st)
    nFilters = [size(layers[i].weight, 3) - 1 for i = 1:3:(3*Dim)]
    batchSize = getBatchSize(layers[1])
    print(io, "stFlux{$(Dep), Nd=$(Nd), filters=$(nFilters), σ = " *
              "$(σ), batchSize = $(batchSize), normalize = $(st.normalize)}")
end

# the type T is a type of frame transform that forms the backbone of the transform
# the type Dimension<:Integer gives the dimension of the transform

@doc """
        st(X::Array{T, N}, layers::stParallel, nonlinear::nl; fullOr::fullType=fullType(),# subsam::Sub = bspline(), thin::Bool=true, outputSubsample::Tuple{Int, Int}=(-1,-1), subsam::Bool=true, totalScales = [-1 for i=1:depth(layers)+1], percentage = .9, fftPlans = -1) where {T <: Real, S <: Union, N, nl <: Function, Sub <: resamplingMethod}
1D scattering transform using the stParallel layers. you can switch out the nonlinearity as well as the method of subsampling. Finally, the stType is a string. If it is "full", it will produce all paths. If it is "decreasing", it will only keep paths of increasing scale. If it is "collating", you must also include a vector of matrices.
# Arguments
- `nonlinear` : a type of nonlinearity. Should be a function that acts on Complex numbers
- `thin` : determines whether to wrap the output into a format that can be indexed using paths. `thin` cannot.
- `totalScales`, if positive, gives the number of non-averaging wavelets.
- `outputSubsample` is a tuple, with the first number indicating a rate, and the second number giving the minimum allowed size. If some of the entries are less than 1, it has different behaviour:
    + `(<1, x)` : subsample to x elements for each path.
    + `(<1, <1)` : no ssubsampling
    + `(x, <1)` : subsample at a rate of x, with at least one element kept in each path
- `fullOr::fullType=fullType()` : the structure of the transform either
       `fullType()`, `collatingType()` or `decreasingType()`. At the moment,
       only `fullType()` is functional.
- `fftPlans = false` if not `false`, it should be a 2D array of `Future`s, where the first index is the layer, and the second index the core. See `createFFTPlans` if you want to do this.
"""
struct stParallel{T,Dimension,Depth,subsampType,outType} <: scatteringTransform{Dimension,Depth}
    n::Tuple{Vararg{Int,Dimension}} # the dimensions of a single entry
    shears::Array{T,1} # the array of the transforms; the final of these is
    # used only for averaging, so it has length m+1
    subsampling::subsampType # for each layer, the rate of
    # subsampling. There is one of these for layer zero as well, since the
    # output is subsampled, so it should have length m+1
    outputSize::outType # a list of the size of a single output example
    # dimensions in each layer. The first index is layer, the second is
    # dimension (e.g. a 3 layer shattering transform is 3×2) TODO: currently
    # unused for the 1D case
end

function Base.show(io::IO, l::stParallel{T,D,Depth}) where {T,D,Depth}
    print(io, "stParallel{$T,$D} depth $(Depth), input size $(l.n), subsampling rates $(l.subsampling), outputsizes = $(l.outputSize)")
end

function eltypes(f::stParallel)
    return (eltype(f.shears), length(f.n))
end

"""
    listVargs = processArgs(m, varargs)
method to go from arguments given to the scattering transform constructor to
those for the frame transform, e.g. shearlets or wavelet. `listVargs` is a list
of length `m` of one argument from each of vargs, with insufficiently long
entries filled in by repeating the last value. An example:
```
julia> varargs
pairs(::NamedTuple) with 3 entries:
  :boundary      => PerBoundary()
  :frameBound    => [1, 1]
  :normalization => (Inf, Inf)

julia> listVargs = processArgs(3,varargs)
(Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol,Symbol},NamedTuple{(:boundary, :frameBound, :normalization),Tuple{PerBoundary,Int64,Float64}}}(:boundary => PerBoundary(),:frameBound => 1,:normalization => Inf), Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol,Symbol},NamedTuple{(:boundary, :frameBound, :normalization),Tuple{PerBoundary,Int64,Float64}}}(:boundary => PerBoundary(),:frameBound => 1,:normalization => Inf), Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol,Symbol},NamedTuple{(:boundary, :frameBound, :normalization),Tuple{PerBoundary,Int64,Float64}}}(:boundary => PerBoundary(),:frameBound => 1,:normalization => Inf))

julia> listVargs[1]
pairs(::NamedTuple) with 3 entries:
  :boundary      => PerBoundary()
  :frameBound    => 1
  :normalization => Inf

julia> listVargs[2]
pairs(::NamedTuple) with 3 entries:
  :boundary      => PerBoundary()
  :frameBound    => 1
  :normalization => Inf

julia> listVargs[3]
pairs(::NamedTuple) with 3 entries:
  :boundary      => PerBoundary()
  :frameBound    => 1
  :normalization => Inf
```
"""
function processArgs(m, varargs)
    keysVarg = keys(varargs)
    valVarg = map(v -> makeTuple(m, v), values(varargs))
    pairedArgs = Iterators.Pairs(valVarg, keysVarg)
    listVargs = ntuple(i -> Iterators.Pairs(map(x -> x[i], valVarg), keysVarg), m)
end

function makeTuple(m, v::Tuple)
    if length(v) >= m
        return v
    end
    return (v..., ntuple(i -> v[end], m - length(v))...)
end
makeTuple(m, v::AbstractArray) = makeTuple(m, (v...,))
makeTuple(m, v) = ntuple(i -> v, m)

# shouldn't be necessary
# function cwtDefaults(varargs)
#     keyLimes = keys(varargs)
#     DefaultList = ((:averagingLength,4),
#                    (:averagingType,ContinuousWavelets.Father()),
#                    (:boundary,DEFAULT_BOUNDARY),
#                    (:frameBound,1),
#                    (:normalization,Inf),
#                    (:decreasing, 4.0))
#     for (key,val) in DefaultList

#     end
# end


text:
#=
PostProcessing:
- Julia version:
- Author: dgleich
- Date: 2020-05-22
=#

#Notes: change b-matching search order to k=1:n, i = 1:n , j = k+1 - i (verify)

"""-----------------------------------------------------------------------------
   Greedily creates a b matching from a rank 1 bipartite matching problem. Only
  adds pairs to the matchings if they have the same sign, and neither edge value
  is zero.

  Output:
  -------
  * Matchings - (2D Int Array):
    The ith row corresponds to the b-matchings of the ith vertex represented by
    the vector u to the vertices in v. unmatched vertices are set to -1.

-----------------------------------------------------------------------------"""
function b_matching(u::Array{Float64,1}, v::Array{Float64,1}, b::Int)

    @assert length(u) <= length(v) #make first vector shorter

    n = length(u)

    #get indices sorted by increasing order
    u_indices = sort(1:length(u),by=i->-u[i])
    v_indices = sort(1:length(v),by=i->-v[i])

    matching_degrees = zeros(length(v))

    Matchings = Array{Tuple{Int,Int,Float64},1}(undef,b*n)  #-1 indicates unset matching
    match_idx = 1
    starting_index = 0b1
    ending_index = b

    for i in 1:n

        if u[i] == 0.0
            continue
        else
            #check matching degree of starting index
            while matching_degrees[v_indices[starting_index]] >= b
                starting_index += 1
            end
            ending_index = minimum((starting_index + b,length(v)))
            #check that sign of ending vertex doesn't change
            while v[v_indices[ending_index]]*u[u_indices[i]] <= 0
                ending_index -= 1
            end

            @assert ending_index != starting_index  #make sure something's added

            for j in 1:(ending_index - starting_index)
                v_j = v_indices[starting_index + (j-1)]
                Matchings[match_idx] = (u_indices[i],v_j,v[v_j]*u[u_indices[i]])
                match_idx += 1
                matching_degrees[v_j] += 1
            end
        end
    end

    return Matchings[1:(match_idx-1)]

end

"""-----------------------------------------------------------------------------
    Computes the cost of swapping the matching of (i,j) & (k,l) to (i,l) & (k,j)
  in terms of the number of triangles. Inputs are the triangles associated with
  each of the vertices

  Inputs:
  -------
  * T_{i,j,k,l} -(Set{Array{Int,1}}):
    All triangles incident with the vertices i,j,k,l .
  * original_matching - (Dict{Int,Int})
    The original matching from the vertices in graph A to the vertices in graph
    B.

-----------------------------------------------------------------------------"""
function compute_triangles_aligned(T_i::Set{Array{Int,1}},T_j::Set{Array{Int,1}},
  #                                  T_k::Set{Array{Int,1}},T_l::Set{Array{Int,1}},
                                    M::Dict{Int,Int})
    #compute the original number of aligned triangles
    matched_tris = 0
    for (v_i,v_j,v_k) in T_i
        if sort([get(M,v_i,-1),get(M,v_j,-1),get(M,v_k,-1)]) in T_j
            matched_tris += 1
        end
    end

    #=
    for (v_i,v_j,v_k) in T_k
        if sort([get(M,v_i,-1),get(M,v_j,-1),get(M,v_k,-1)]) in T_l
            matched_tris += 1
        end
    end

    #compute the new matched_triangles end
    new_matched_tris = 0
    for (v_i,v_j,v_k) in T_i
        if sort([get(M,v_i,-1),get(M,v_j,-1),get(M,v_k,-1)]) in T_l
            new_matched_tris += 1
        end
    end
    for (v_i,v_j,v_k) in T_k
        if sort([get(M,v_i,-1),get(M,v_j,-1),get(M,v_k,-1)]) in T_j
            new_matched_tris += 1
        end
    end
    =#
    return matched_tris#, new_matched_tris
end

"""-----------------------------------------------------------------------------
     Finds all the triangles associated with the vertices and preprocesses them
   into a list of Sets. Used to preprocess the searching for triangles
   associated with a vertex.

   Inputs:
   -------
   * triangles - (Array{Int,2}):
     The triangles contained within the graph. indices in each row are expected
     to be sorted in ascending order.
   * n - (Int):
     The number of vertices in the graph.

   Outputs:
   --------
   * incidence - (Array{Set{Array{Int,1}}},1):
     An array where the ith entry contains a set linking the vertex to the
     triangles it's contained within.
-----------------------------------------------------------------------------"""
function produce_triangle_incidence(triangles::Array{Int,2},n::Int)
    incidence = Array{Set{Array{Int,1}},1}(undef,n)
    for i in 1:n
        incidence[i] = Set{Array{Int,1}}()
    end

    for t in eachrow(triangles)
        for v in t
            push!(incidence[v],t)
        end
    end

    return incidence
end