text:
using PyCall
using PyPlot

"""
safe_axes(axh; further_params...)

If you're going to make axh the current axes, this function
first makes axh's figure the current figure. Some Julias
break without that.

Any optional keyword-value arguments are passed on to axes()
"""

function safe_axes(axh; further_params...)
    figure(axh[:figure][:number])
    axes(axh; Dict(further_params)...)
end

"""
    ax = axisWidthChange(factor; lock="c", ax=nothing)
"""
function axisWidthChange(factor; lock="c", ax=nothing)
    if ax==nothing; ax=gca(); end
    x, y, w, h = ax[:get_position]()[:bounds]

    if lock=="l";
    elseif lock=="c" || lock=="m"; x = x + w*(1-factor)/2;
    elseif lock=="r"; x = x + w*(1-factor);
    else error("I don't know lock type ", lock)
    end

    w = w*factor;
    ax[:set_position]([x, y, w, h])

    return ax
end

"""
ax = axisHeightChange(factor; lock="c", ax=nothing)
"""
function axisHeightChange(factor; lock="c", ax=nothing)
    if ax==nothing; ax=gca(); end
    x, y, w, h = ax[:get_position]()[:bounds]

    if lock=="b";
    elseif lock=="c" || lock=="m"; y = y + h*(1-factor)/2;
    elseif lock=="t"; y = y + h*(1-factor);
    else error("I don't know lock type ", lock)
    end

    h = h*factor;
    ax[:set_position]([x, y, w, h])

    return ax
end

"""
   ax = axisMove(xd, yd; ax=nothing)
"""
function axisMove(xd, yd; ax=nothing)
    if ax==nothing; ax=gca(); end
    x, y, w, h = ax[:get_position]()[:bounds]

    x += xd
    y += yd

    ax[:set_position]([x, y, w, h])
    return ax
end

"""
[] = remove_xtick_labels(ax=NaN)

Given an axis object, or an array of axes objects, replaces each xtick label string with the empty string "".

If no axis is passed, uses gca() to work with the current axis.

"""
function remove_xtick_labels(ax=nothing)

    if ax==nothing
        ax = gca()
    end

    if typeof(ax) <: Array
        for i=1:length(ax)
            remove_xtick_labels(ax[i])
        end
        return
    end

    nlabels = length(ax[:xaxis][:get_ticklabels]())

    newlabels = Array{String,1}(nlabels)
    for i=1:length(newlabels);
        newlabels[i] = ""
    end

    ax[:xaxis][:set_ticklabels](newlabels)
    return
end

"""
[] = remove_ytick_labels(ax=NaN)

Given an axis object, or an array of axes objects, replaces each ytick label string with the empty string "".

If no axis is passed, uses gca() to work with the current axis.

"""
function remove_ytick_labels(ax=nothing)

    if ax==nothing
        ax = gca()
    end

    if typeof(ax) <: Array
        for i=1:length(ax)
            remove_ytick_labels(ax[i])
        end
        return
    end

    nlabels = length(ax[:yaxis][:get_ticklabels]())

    newlabels = Array{String,1}(nlabels)
    for i=1:length(newlabels);
        newlabels[i] = ""
    end

    ax[:yaxis][:set_ticklabels](newlabels)
    return
end

# If the Python path does not already have the local directory in it
if PyVector(pyimport("sys")["path"])[1] != ""
    # Then the following line is PyCall-ese for "add the current directory to the Python path"
    unshift!(PyVector(pyimport("sys")["path"]), "")
end
# We use Python to enable callbacks from the figures:
@pyimport kbMonitorModule

__permanent_BP_store = []   # The user doesn't need to worry about this variable, it is here to ensure that
                            # kbMonitorModule.kb_monitor objects created inside the install_nearest_callback()
                            # function do not get deleted upon exit of that function

"""
BP = install_nearest_point_callback(fighandle, user_callback; user_data=nothing)

This function makes the figure indicated by fighandle interactive: any time the mouse is clicked
while pointing inside any of the axes of the figure, the function user_callback() will be called,
and will be passed parameters that indicate which of the points drawn in the axis is the one
closest to the clicked point (closest in Euclidean data units).

**WARNING** because this runs through PyCall, any errors in your user_callback function will sadly
not show up.  The function will simply fail to work. So be careful and debug with lots of print statements.

# PARAMETERS:

- fighandle       A matplotlib figure handle, e.g. the result of figure(2)

- user_callback   A function, which must take 4 or 5 parameters (see below). These will be passed to it as:

* PARAMETERS OF YOUR FUNCTION USER_CALLBACK:

        - xy          A 2-element tuple, indicating the (x,y) position of the drawn point closest to the clicked point

        - r           A scalar, indicating the Euclidean distance between the clicked location and xy

        - linehandle  A matplotlib handle to a Lines2D object (e.g., the result of plot([1,2], [3, 10])) or a PathCollection object (as returned by scatter()).

        - axhandle    A matplotlib handle to the axis (e.g., the result of gca()) in which the event occurred.

        - user_data   If `install_nearest_point_callback()` was called user_data set to something, then
                      your function will be called with *five* parameters, and the last one will be the contents of
                      the user_data

# OPTIONAL PARAMETERS FOR INSTALL_NEAREST_POINT_CALLBACK():

- user_data       Data to be stored internally and made available to the callback function. Default
                  is nothing, in which case the callback function is called with 4 params

# RETURNS:

- BP     A PyCall.PyObject kbMonitorModule_kb_monitor object. This object contains the underlying engine
linking the figure to the callback function. To disconnect that link, call "remove_BP(BP)". To disconnect
*all* existing BP-function links, call "remove_all_BPs()".

# EXAMPLE:

pygui(true)

```jldoctest
function mycallback(xy, r, h, ax)
    @printf("(%.3f,%.3f), r=%3f ", xy[1], xy[2], r);
    print(h)
    print(ax)
    print("\n")
end

BP = install_nearest_point_callback(figure(2), mycallback)
plot([2,2])
```

"""
function install_nearest_point_callback(fighandle, user_callback; user_data=nothing)

    function point_nearest_to_click(BP)
        bpe = BP[:buttonlist]()
        # Remove any leading clicks that weren't inside an axis:
        while length(bpe)>0 && ((bpe[1][1]==nothing) || (bpe[1][1]==Void))
            bpe = bpe[2:end]
        end
        if length(bpe)>0
            ax = BP[:buttonlist]()[1][1]   # the axis we're working with
            x  = BP[:buttonlist]()[1][2]   # the x of the clicked point
            y  = BP[:buttonlist]()[1][3]   # the y of the clicked point

            ch = ax[:get_children]()       # all children of the axis

            idx    = nothing    # this'll be the index of the data point closest to the clickpoint
            minJ   = nothing    # the smallest squared distance between data point and clickpoint found so far
            handle = nothing    # the matplotlib handle of the line object for which the closes data point is found
            dx     = nothing    # closest data point x position
            dy     = nothing    # closest data point y position

            # Look over all children of the axis:
            for i=1:length(ch)
                # But only consider line objects:
                if contains(pystring(ch[i]), "lines.Line2D")
                    D = ch[i][:get_data]()    # D will be a Tuple with xdata, ydata vectors
                elseif contains(pystring(ch[i]), "PathCollection")
                    D = ch[i][:get_offsets]()    # D will be a matrix with xdata, ydata columns
                    D = (D[:,1], D[:,2])         # Turn it into a Tuple like for Line2D objects
                end
                if contains(pystring(ch[i]), "lines.Line2D") || contains(pystring(ch[i]), "PathCollection")
                    J = (D[1] - x).^2 + (D[2] - y).^2
                    ix = indmin(J)
                    if idx == nothing || J[ix] < minJ   # if we did not yet have a minimum candidate or this one is better
                        idx = ix; minJ = J[ix]; handle = ch[i]   # store our candidate
                        dx = D[1][ix]; dy = D[2][ix]
                    end
                end
            end

            # @printf("install: Am about to call the user callback\n")
            if minJ != nothing
                if BP[:get_userdata]() == nothing
                    user_callback((dx,dy), sqrt(minJ), handle, ax)
                else
                    user_callback((dx,dy), sqrt(minJ), handle, ax, BP[:get_userdata]())
                end
            end
            # @printf("install: Just returned from the user callback\n")

            # After dealing with all the buttonclick callbacks and so on, bring focus back to the figure that was clicked:
            figure(ax[:figure][:number])
        end

        # We've dealt with the buttonclick, clear the buttonlist
        # @printf("Am about to clear the button list on button "); print(BP); print("\n")
        BP[:clear_buttonlist]()
    end

    BP = kbMonitorModule.kb_monitor(fighandle, callback = point_nearest_to_click, userData=user_data)
    global __permanent_BP_store = [__permanent_BP_store ; BP]

    return BP
end

"""
    userdata = install_callback_reporter(xy, r, axhandle, dothandle, userdata)

Useful as a debugging tool for `install_nearest_point_callback()`: can be used
with that function as a callback; when called, simply prints out its parameters.
Returns userdata, does not print it.

"""
function install_callback_reporter(xy, r, linehandle, axhandle, userdata=nothing)
    @printf("xy=(%g,%g), r=%g\n", xy[1], xy[2], r)
    print("Line Handle:\n"); print(linehandle); print("\n")
    print("Axis Handle:\n"); print(axhandle); print("\n")
    # print("User Data:\n"); print(userdata); print("\n")

    return userdata
end

"""
    remove_BP(BP::PyCall.PyObject)

Disconnects a kbMonitorModule.kb_monitor object from its figure
"""
function remove_BP(BP::PyCall.PyObject)
    if contains(pystring(BP), "kbMonitorModule.kb_monitor")
        BP[:__del__]()

        i = find(__permanent_BP_store .== BP)
        if length(i)>0;
            i = i[1];
            global __permanent_BP_store = __permanent_BP_store[[1:(i-1) ; (i+1):end]]
        end
    end
end

"""
    remove_all_BPs()

    Disconnects *all* kbMonitorModule.kb_monitor objects from their figures.
"""
function remove_all_BPs()
    for BP in __permanent_BP_store
        BP[:__del__]()
    end

    global __permanent_BP_store = []
end

"""

    (x, y, w, h) = get_current_fig_position()

Works only when pygui(true) and when the back end is QT. Has been tested only with PyPlot.
"""
function get_current_fig_position()
    # if !contains(pystring(plt[:get_current_fig_manager]()), "FigureManagerQT")
    try
        if contains(pystring(plt[:get_current_fig_manager]()), "Tk")
            g = split(plt[:get_current_fig_manager]()[:window][:geometry](), ['x', '+'])
            w = parse(Int64, g[1])
            h = parse(Int64, g[2])
            x = parse(Int64, g[3])
            y = parse(Int64, g[4])
        elseif contains(pystring(plt[:get_current_fig_manager]()), "QT")
            x = plt[:get_current_fig_manager]()[:window][:pos]()[:x]()
            y = plt[:get_current_fig_manager]()[:window][:pos]()[:y]()
            w = plt[:get_current_fig_manager]()[:window][:width]()
            h = plt[:get_current_fig_manager]()[:window][:height]()
        else
            error("Only know how to work with matplotlib graphics backends that are either Tk or QT")
        end

        return (x, y, w, h)
    catch
        error("Failed to get current figure position. Is pygui(false) or are you using a back end other than QT or Tk?")
    end
end

"""

    set_current_fig_position(x, y, w, h)

Works only when pygui(true) and when the back end is QT. Has been tested only with PyPlot.
"""
function set_current_fig_position(x, y, w, h)
    # if !contains(pystring(plt[:get_current_fig_manager]()), "FigureManagerQT")
    try
        if contains(pystring(plt[:get_current_fig_manager]()), "Tk")
            plt[:get_current_fig_manager]()[:window][:geometry](@sprintf("%dx%d+%d+%d", w, h, x, y))
        elseif contains(pystring(plt[:get_current_fig_manager]()), "QT")
            plt[:get_current_fig_manager]()[:window][:setGeometry](x, y, w, h)
        else
            error("Only know how to work with matplotlib graphics backends that are either Tk or QT")
        end
    catch
        error("Failed to set current figure position. Is pygui(false) or are you using a back end other than QT?")
    end
end

"""
    C = capture_current_figure_configuration()

Collects the positions of all current figures and
prints out to the screen code, that can be copy-pasted,
that would reproduce that positioning configuration.

# PARAMETERS:

None

# RETURNS:

- C    A matrix that is nfigures-by-5 in size. You probably
    don't want this, you probably want the text printed to
    the screen, but here just in case.  Each row will have,
    in order: figure number, x, y, width, height
"""
function capture_current_figure_configuration()
    @printf("The following code will reproduce your current figure placement:\n\n")
    C = []
    for f in sort(plt[:get_fignums]())
        figure(f)
        x, y, w, h = get_current_fig_position()
        @printf("figure(%d); set_current_fig_position(%d, %d, %d, %d)   # x, y, width, height\n",
            f, x, y, w, h)
        C = [C ; [f x y w h]]
    end
    return C
end


text:
module MatrixCoalescent

using LinearAlgebra

export MatCoal, eigenvals!, project!, interval_lengths!, dim

"""
An object of type MatCoal holds the data needed to calculate two
quantities under the matrix coalescent: (1) the probability that there
are k lineages at the ancient end of an epoch, given that there are
nLineages at the recent end, and (2) the expected duration of the
subinterval during which there are k lineages.
"""
mutable struct MatCoal{T<:Real}
    nLineages :: Int # number of lineages at recent end of epoch

    # Vector of dimension nLineages-1. beta[i] is (i+1) choose 2.  For
    # example, if nLineages=3, then beta has two entries: 2 choose 2
    # and 3 choose 2.
    beta :: Vector{T}

    # Vector for calculating E[len] of coalescent intervals
    # Reciprocal of beta.
    rbeta :: Vector{T}

    # Matrix of scaled column eigenvectors (upper triangular)
    gmat :: Matrix{T}

    # Matrix for calculating expected lengths of coalescent intervals.
    # Also upper triangular.
    hmat :: Matrix{T}
end

"""
    MatCoal(el_type::DataType, nLineages)

Construct an object of type `MatCoal`, which describes an epoch of
population history during which the population's size was
constant. The first argument is the type to be used in representing
probabilities and expectations. The constructor will accept any Real
type, but the other functions require a subtype of AbstractFloat. The
second argument is the number of lineages at the recent end of this
epoch of population history.
"""
function MatCoal(el_type::DataType, nLineages)

    nLin = Int(nLineages)
    n = nLin - 1

    # Column eigenvectors and row eigenvectors
    cvec = zeros(Rational{Int}, n, n)
    rvec = zeros(Rational{Int}, n, n)
    for j = 2:nLin
        jj = j-1
        lambda = -j*(j-1) # eigenvalue
        cvec[jj,jj] = rvec[jj,jj] = 1//1
        for i in j-1 : -1 : 2
            ii = i-1
            m = (i*(i+1)) // (i*(i-1) + lambda)
            cvec[ii, jj] = m * cvec[ii+1, jj]
        end
        for i in j+1 : nLin
            ii = i-1
            m = (i*(i-1)) // (i*(i-1) + lambda)
            rvec[jj, ii] = m * rvec[jj, ii-1]
        end
    end

    # Calculate coefficients of exponentials in x(t).  This will be
    # converted to float to produce the the matrix gmat.
    for ii in 1:n
        for jj in ii:n
            cvec[ii,jj] *= rvec[jj, n]
        end
    end

    # rbeta[i] is reciprocal of beta[i] = (i+1) choose 2.
    beta = Vector{Rational{Int}}(undef, n)
    rbeta = Vector{Rational{Int}}(undef, n)
    for i in 1:n
        j = i+1
        beta[i] = (j*(j-1)) // 2
        rbeta[i] = 1/beta[i]
    end

    # Initially hmat == cvec
    hmat = Matrix{Rational{Int}}(cvec)

    # Cumulative sum, so i'th row is sum of i:n initial rows
    for i in n-1 : -1 : 1
        for j in 1:n
            hmat[i,j] += hmat[i+1, j]
        end
    end

    # Weight row i by -1/beta[i]
    for i in 1:n
        for j in 1:n
            hmat[i,j] *= -rbeta[i]
        end
    end

    MatCoal(nLin,
            Vector{el_type}(beta),
            Vector{el_type}(rbeta),
            Matrix{el_type}(cvec),
            Matrix{el_type}(hmat))
end

"""
Dimension of a MatCoal object is 1 less than the number (nLin) of
lineages at the recent end of the epoch.
"""
function dim(mc::MatCoal)
    return length(mc.beta)
end

"""
Write a representation of an MatCoal to io.
"""
function Base.show(io::IO, mc::MatCoal)
    println(io, "nLineages:", mc.nLineages)
    print(io, "beta:")
    for b in mc.beta
        print(io, " ", b)
    end
    println(io)
    print(io, "rbeta:")
    for nrb in mc.rbeta
        print(io, " ", nrb)
    end
    println(io)
    println(io, "gmat:")
    println(io, mc.gmat)
    println(io, "hmat:")
    println(io, mc.hmat)
end

"""
    eigenvals!(eig::Vector{T}, v::T,
                    mc::MatCoal{T}) where {T<:AbstractFloat}

Calculate eigenvalues for time ``v``, where ``v = t/2N`` is time
measured backwards from the recent end of the epoch in units of
``2N`` generations. It represents the length of the epoch and may be
0 or infinite. Time affects the coalescent process only via its effect
on eigenvalues. Consequently, this is the only function that takes
time as an argument.
"""
function eigenvals!(eig::Vector{T}, v::T,
                    mc::MatCoal{T}) where {T<:AbstractFloat}
    n = length(eig)
    @assert n == dim(mc)
    for i in 1:n
        eig[i] = exp(-v*mc.beta[i]);
    end
    return eig
end

"""
    project!(ans::Vector{T}, eig::Vector{T},
                 mc::MatCoal{T}) where {T <: AbstractFloat}

Calculate the probability that there are 2,3,...(dim+1) lines of
descent at the ancient end of the epoch. Call eigenvals! first to
calculate eig. The length of the interval affects the calculation only
through its effect on the eigenvalues in eig, so it does not appear as
an argument to project!.
"""
function project!(ans::Vector{T}, eig::Vector{T},
                 mc::MatCoal{T}) where {T <: AbstractFloat}
    @assert length(ans) == length(eig)
    @assert length(ans) == dim(mc)

    mul!(ans, mc.gmat, eig)
end

"""
    interval_lengths!(ans::Vector{T}, eig::Vector{T},
                 mc::MatCoal{T}) where {T <: AbstractFloat}

Vector of expected lengths of coalescent intervals during which there
were 2,3,...(dim+1) lines of descent. To get the expected length of
the interval with 1 line of descent, subtract the sum of ans from v,
where v is the time argument given to eigenvals!.  Call eigenvals!
first to calculate eig.  The length of the interval affects the
calculation only via the eigenvalues in eig and therefore does not
appear as an argument in interval_lengths;
"""
function interval_lengths!(ans::Vector{T}, eig::Vector{T},
                 mc::MatCoal{T}) where {T <: AbstractFloat}
    @assert length(ans) == length(eig)
    @assert length(ans) == dim(mc)

    # ans = rbeta + hmat*eig
    copy!(ans, mc.rbeta)
    mul!(ans, mc.hmat, eig, 1, 1)
end

end # module

ext: jl


text:
using ITensors

# Define indices
a = Index(2, "a")
b = Index(2, "b")
c = Index(2, "c")

# Note that Index tags don't have
# to match variable names
i = Index(3, "Red,Link")

# Define 3 order 2 tensors (matrices)
Z = ITensor(a, b)
X = ITensor(b, c)
Y = ITensor(b, c)

# Set some elements
Z[a(1), b(1)] = 1.0
Z[a(2), b(2)] = -1.0

X[b(1), c(2)] = 1.0
X[b(2), c(1)] = 1.0

Y[b(1), c(1)] = 1.0
Y[b(2), c(2)] = 1.0

# Operations with tensors
R = Z * X

S = Y + X

T = Y - X

# Print results
println("Z =\n", Z, "\n")
println("X =\n", X, "\n")
println("Y =\n", Y, "\n")
println("R = Z * X =\n", R, "\n")
println("S = Y + X =\n", S, "\n")
println("S = Y - X =\n", T, "\n")

# Check that adding incompatible tensors cause an error
try
  U = Z + X
catch
  println("Cannot add Z and X")
end

# Compare calculations to Julia arrays
jZ = [
  1.0 0.0
  0.0 -1.0
]
jX = [
  0.0 1.0
  1.0 0.0
]
jY = [
  1.0 0.0
  0.0 1.0
]
@assert Array(R, a, c) == jZ * jX
@assert Array(S, b, c) == jY + jX
@assert Array(T, b, c) == jY - jX


text:
module HavenOnDemand

using JSON
using Requests

       # Types
export HODClient,
       HODIndex,
       HODSyncResponse,
       HODAsyncResponse,
       HODException,
       DocumentException,

       # Client
       create_index,
       parse_index,
       has_index,
       get_index,
       delete_index,
       list_indexes,
       post,
       post_async,
       call_API,
       get_async_result,

       # Index
       size,
       push_doc,
       push_docs,
       commit,
       add_doc,
       add_docs,
       delete,

       # Responses
       json,
       get_result

include("types.jl")
include("client.jl")
include("responses.jl")
include("index.jl")
include("exceptions.jl")

end # module


text:
# Function that updates partial force and energy vectors, used for parallel computation

function compute_uf_partials_LJ!(it :: Int64, f :: Forces,
                                 i :: Int64, j :: Int64,
                                 xj :: Float64, yj :: Float64, r :: Float64,
                                 input :: MDInput)
  r6 = r^6
  r12 = r6^2
  r7 = r6*r
  r13 = r12*r

  f.upartial[it] = f.upartial[it] + f.eps4*( f.sig12 / r12 - f.sig6 / r6 )

  drdxi = -xj/r
  drdyi = -yj/r

  fx = -f.eps4*( (f.sig1212/r13) - (f.sig66/r7) )*drdxi
  fy = -f.eps4*( (f.sig1212/r13) - (f.sig66/r7) )*drdyi

  f.fpartial[it][i,1] = f.fpartial[it][i,1] + fx
  f.fpartial[it][i,2] = f.fpartial[it][i,2] + fy
  f.fpartial[it][j,1] = f.fpartial[it][j,1] - fx
  f.fpartial[it][j,2] = f.fpartial[it][j,2] - fy

end


text:
function load_messages(client::OutlookClient, folderId::AbstractString)
    msgs = Array{OutlookMessage, 1}()

    headers = Dict(
        "Authorization"     => build_token(client),
        "Accept"            => "application/json",
        "X-AnchorMailbox"   => emailaddress(client),
    )

    response = get(API * "MailFolders/$folderId/messages", headers=headers)
    json_response = JSON.parse(readall(response))

    for (key, value) in enumerate(json_response["value"])
        msg = OutlookMessage(
            value["Id"],
            value["Subject"],
            value["IsDraft"],
            value["HasAttachments"],
            haskey(value["From"]["EmailAddress"], "Address") ? value["From"]["EmailAddress"]["Address"] : "Email Address was not found.",
            value["From"]["EmailAddress"]["Name"],
            value["Body"]["Content"],
            value["IsRead"],
            value["Importance"],
        )

        push!(msgs, msg)
    end

    msgs
end

messageid(msg::OutlookMessage)      = msg.id
subject(msg::OutlookMessage)        = msg.subject
isdraft(msg::OutlookMessage)        = msg.draft
hasattachments(msg::OutlookMessage) = msg.hasAttachments
fromaddress(msg::OutlookMessage)    = msg.fromAddress
fromname(msg::OutlookMessage)       = msg.fromName
body(msg::OutlookMessage)           = msg.body
isread(msg::OutlookMessage)         = msg.isread
importance(msg::OutlookMessage)     = msg.importance


text:
using ArtifactUtils, Artifacts

url = "https://github.com/polyanskiy/refractiveindex.info-database/archive/v2021-07-18.tar.gz"

add_artifact!(
    joinpath(@__DIR__, "../Artifacts.toml"),
    "refractiveindex.info",
    url,
    force=true,
    clear=false,
)


text:
# ---
# title: ART Options Example
# id: options
# cover: ../assets/art.png
# date: 2021-12-2
# author: "[Sasha Petrenko](https://github.com/AP6YC)"
# julia: 1.6
# description: This demo illustrates how to use options and modify the options for all ART and ARTMAP modules.
# ---

# ## Overview

# The `AdaptiveResonance.jl` package has several ways of handling options for ART modules.
# These methods are meant to give maximum flexibility to the user for sharing and interpreting options, which themselves vary between each module.

# !!! note
#     For more info on options in ART modules, see the guide in the docs on [ART options](@ref art_options).

# ## ART Options

# To get a feel for the ART options system, we will inspect different options and their instantiation methods.

# ### Inspection

# First, we load `AdaptiveResonance`:
using AdaptiveResonance

# Every ART module has a default constructor, which can be instantiated in the usual way:

## Create a FuzzyART module with default options
my_fuzzyart = FuzzyART()
typeof(my_fuzzyart)

# Within every ART module is a [Parameters.jl](https://github.com/mauro3/Parameters.jl) struct named `opts` containing the options for the module

## Check the FuzzyART options
my_fuzzyart.opts

# Note that the options here have the type `opts_FuzzyART`.
# This nomenclature is used throughout the module to indicate an options type associated with an ART module.
# For example, the options for a DDVFA module are `opts_DDVFA`:

## Create a DDVFA module and check the type of the options
my_ddvfa = DDVFA()
typeof(my_ddvfa.opts)

# In fact, we can create an instance of these options with a default constructor:

## Create a separate options struct
my_fuzzyart_opts = opts_FuzzyART()

# In addition to the default constructor, we can construct ART modules by instantiating these options and passing them to the module during construction:

## Instantiate an ART module by passing our options
my_fuzzyart = FuzzyART(my_fuzzyart_opts)
my_other_fuzzyart = FuzzyART(my_fuzzyart_opts)

# ### Specifying Options

# Now to the good stuff: because of the behavior of the `Parameters.jl` type, each option has a default value that we can modify during instantiation with keyword arguments:

## Change some of the default FuzzyART options
my_fuzzyart_opts = opts_FuzzyART(
    rho=0.6,
    gamma_normalization=true
)
my_fuzzyart = FuzzyART(my_fuzzyart_opts)

# As some syntactic sugar, we can pass these keyword arguments directly to the module during instantiation if we have no need to share option structs:

## Pass these keyword arguments to the module directly
my_fuzzyart = FuzzyART(
    rho=0.6,
    gamma_normalization=true
)

# Before training, we can also instantiate the model and alter the options afterward:
my_fuzzyart = FuzzyART()
my_fuzzyart.opts.rho=0.6

# !!! note
#     All ART modules are designed to use this options struct internally when the parameters are needed.
#     It is possible to change these parameters in the middle of training and evaluation, but some algorithmic instability may occur.

# ## Comparison

# To see the effect that changing these parameters has on the modules, we can train and test them side-by-side.

# We begin with importing AdaptiveResonance for the ART modules and MLDatasets for some data utilities.
using MLDatasets        # Iris dataset
using MLDataUtils       # Shuffling and splitting
using Printf            # Formatted number printing
using MultivariateStats # Principal component analysis (PCA)
using Plots             # Plotting

# We will download the Iris dataset for its small size and benchmark use for clustering algorithms.
Iris.download(i_accept_the_terms_of_use=true)
features, labels = Iris.features(), Iris.labels()

# Because the MLDatasets package gives us Iris labels as strings, we will use the `MLDataUtils.convertlabel` method with the `MLLabelUtils.LabelEnc.Indices` type to get a list of integers representing each class:
labels = convertlabel(LabelEnc.Indices{Int}, labels)
unique(labels)

# Next, we will create a train/test split with the `MLDataUtils.stratifiedobs` utility:
(X_train, y_train), (X_test, y_test) = stratifiedobs((features, labels))

# Now we can create several FuzzyART modules with different options.

## Create two FuzzyARTs with different vigilance values and suppressing logging messages
rho_1 = 0.5
rho_2 = 0.7
my_fuzzyart_1 = FuzzyART(rho=rho_1, display=false)
my_fuzzyart_2 = FuzzyART(rho=rho_2, display=false)

# Here, we will train these FuzzyART modules in simple supervised mode by passing the supervised labels as a keyword argument:

## Train in simple supervised mode by passing the labels as a keyword argument.
y_hat_train_1 = train!(my_fuzzyart_1, X_train, y=y_train)
y_hat_train_2 = train!(my_fuzzyart_2, X_train, y=y_train)

# We then classify the test data with both modules:
y_hat_1 = AdaptiveResonance.classify(my_fuzzyart_1, X_test, get_bmu=true)
y_hat_2 = AdaptiveResonance.classify(my_fuzzyart_2, X_test, get_bmu=true)

## Check the shape and type of the output labels
println("FuzzyART 1 labels: ",  size(y_hat_1), " ", typeof(y_hat_1))
println("FuzzyART 2 labels: ",  size(y_hat_2), " ", typeof(y_hat_2))

## Calculate the performance on the test data
perf_test_1 = performance(y_hat_1, y_test)
perf_test_2 = performance(y_hat_2, y_test)

## Format each performance number for comparison
@printf "Testing performance rho=%.1f: %.4f\n" rho_1 perf_test_1
@printf "Testing performance rho=%.1f: %.4f\n" rho_2 perf_test_2

# In addition to having different performances, we can see that there is a subsequent trade-off in the number of categories used:
## Print the number of categories for each vigilance parameter
@printf "Number of categories rho=%.1f: %i\n" rho_1 my_fuzzyart_1.n_categories
@printf "Number of categories rho=%.1f: %i\n" rho_2 my_fuzzyart_2.n_categories

# The variation between vigilance parameter, number of categories created during learning, and testing performance/generalization is a central theme in ART-based algorithms.

# ## Visualization

# Now, to visualize how the two models differ in how they partition the data, we can use principal component analysis (PCA) to compress to two plotting dimensions.
# PCA is a method to represent a dataset in a different number of dimensions while preserving the relative separation between datapoints.
# Though most datasets are not able to be effectively transformed down to two dimensions, this technique is useful to get a general sense of how well separated the classes are and how well your algorithm classifies them.
## Train a PCA model to visually separate the features in two dimensions.
M = fit(PCA, features; maxoutdim=2)

## Apply the PCA model to the testing set
X_test_pca = transform(M, X_test)

# We can now plot the PCA'ed test set and label them according to the two FuzzyART's
p1 = scatter(X_test_pca[1,:], X_test_pca[2,:], group=y_hat_1, title=@sprintf "FuzzyART rho = %.1f" rho_1)
p2 = scatter(X_test_pca[1,:], X_test_pca[2,:], group=y_hat_2, title=@sprintf "FuzzyART rho = %.1f" rho_2)
plot(p1, p2, layout=(1, 2), legend = false, xtickfontsize=6, xguidefontsize=8, titlefont=font(8))


text:
@derived_dimension ElectricField Unitful.𝐈^-1*Unitful.𝐋*Unitful.𝐌*Unitful.𝐓^-3
@derived_dimension VectorPotential Unitful.𝐈^-1*Unitful.𝐋*Unitful.𝐌*Unitful.𝐓^-2
@derived_dimension Intensity Unitful.𝐌*Unitful.𝐓^-3

# Unit arithmetic is performed "outside" of @u_str, to ensure type
# stability and allow precompilation
const base_units = Dict{Symbol,Unitful.FreeUnits}(
    :λ => u"nm",
    :I₀ => u"W"/u"cm"^2,
    :E₀ => u"V"/u"m",
    :A₀ => (u"V"/u"m")/(u"rad"/u"s"),
    :τ => u"fs",
    :σ => u"fs",
    :σ′ => u"fs",
    :σmax => NoUnits,
    :σ′max => NoUnits,
    :toff => u"fs",
    :tmax => u"fs",
    :Tmax => NoUnits,
    :T => u"fs",
    :f => u"THz",
    :ν => u"cm"^-1,
    :ω => u"Trad"/u"s",
    :ħω => u"eV",
    :Uₚ => u"eV"
)

const Iau = 3.5094452e16u"W"/u"cm"^2

# We have to use this hack, since UnitfulAtomic does not consider Iau
# as the atomic unit of intensity (which in itself is a contested name
# for the quantity in question).
oneaunit(x) = aunit(x)
oneaunit(x::Intensity) = Iau

function get_unitful_quantity(field_params::Dict{Symbol,Any}, sym::Symbol)
    v = field_params[sym]
    typeof(v) <: Quantity || sym ∉ keys(base_units) ? v : v*oneaunit(1base_units[sym])
end

function shift_unit(u::U, d, n) where {U<:Unitful.Unit}
    i₀ = floor(Int, d/n)
    d = 3i₀
    iszero(d) && return u,0
    for (i,tens) in enumerate(Unitful.tens(u) .+ (d:(-3*sign(d)):0))
        haskey(Unitful.prefixdict, tens) && return U(tens, u.power),i-1 + i₀
    end
    u,0
end

function shift_unit(u::Unitful.FreeUnits, d)
    tu = typeof(u)
    us,ds,a = tu.parameters
    ui = findfirst(u -> u.power>0, us)

    uu,i = shift_unit(us[ui], d, 3)

    Unitful.FreeUnits{(us[1:ui-1]...,uu,us[ui+1:end]...), ds, a}(),i
end

function si_round(q::Quantity; fspec="{1:.4f} {2:s}")
    v,u = ustrip(q), unit(q)
    if !iszero(v)
        u,i = shift_unit(u, log10(abs(v)))
        q = u(q)
    end
    format(fspec, ustrip(q), unit(q))
end

au2si(v, u) = u(v*aunit(u))
au2si(u::Unitful.Units)= au2si(1, u)

au2si_round(v, u; kwargs...) =
    si_round(au2si(v, u); kwargs...)

I2si_round(I; kwargs...) = si_round(I*Iau)

Iaustrip(I::Intensity) = NoUnits(I/Iau)


text:
function solver_default(X,data::AbstractArray{T,2}) where {T<:Union{Missing, <:Number}}
    # likely much larger matrix, using lsqr
    modelinfo = []
    beta = Array{Float64}(undef,size(data,1),size(X,2))

    @showprogress .1 for ch in 1:size(data,1)
        ix = .!ismissing.(data[ch,:])
        beta[ch,:],h = lsmr(X[ix,:],data[ch,ix],log=true)
        push!(modelinfo,h)
    end

    return beta, modelinfo
end

function solver_default(X,data::AbstractArray{T,3}) where {T<:Union{Missing, <:Number}}
       beta = Array{Union{Missing,Number}}(undef,size(data,1),size(data,2),size(X,2))
       @showprogress .1 for ch in 1:size(data,1)
           for t in 1:size(data,2)
               @debug("$(ndims(data,)),$t,$ch")
               ix = .!ismissing.(data[ch,t,:])
               beta[ch,t,:] = X[ix,:] \ data[ch,t,ix]
           end
       end

       modelinfo = [undef] # no history implemented (yet?)
       return beta, modelinfo
end

function solver_b2b(X,data::AbstractArray{T,3},cross_val_reps = 10) where {T<:Union{Missing, <:Number}}

    X,data = dropMissingEpochs(X,data)

    E = zeros(size(data,2),size(X,2),size(X,2))
    W = Array{Float64}(undef,size(data,2),size(X,2),size(data,1))
    println("n = samples = $(size(X,1)) = $(size(data,3))")
    for t in 1:size(data,2)
        for m in 1:cross_val_reps
            k_ix = collect(Kfold(size(data,3),2))
            Y1 = data[:,t,k_ix[1]]
            Y2 = data[:,t,k_ix[2]]
            X1 = X[k_ix[1],:]
            X2 = X[k_ix[2],:]

            G = (Y1' \ X1)
            H = X2 \ (Y2'*G)

            E[t,:,:] = E[t,:,:]+Diagonal(H[diagind(H)])

        end
        E[t,:,:] = E[t,:,:] ./ cross_val_reps
        W[t,:,:] = (X*E[t,:,:])' / data[:,t,:]

    end

    # extract diagonal
    beta = mapslices(diag,E,dims=[2,3])
    # reshape to conform to ch x time x pred
    beta = permutedims(beta,[3 1 2])
    modelinfo = Dict("W"=>W,"E"=>E,"cross_val_reps"=>cross_val_reps) # no history implemented (yet?)
    return beta, modelinfo
end


text:
function bridgeold!(Y::SamplePath, W::SamplePath, P, scheme!)
    !(W.tt[1] == P.t0 && W.tt[end] == P.t1) && error("Time axis mismatch between bridge P and driving W.") # not strictly an error
    scheme!(Y, P.v0, W, P)
    Y.yy[.., length(W.tt)] = P.v1
    Y
end

function bridgeold!(Y::SamplePath, W::SamplePath, P)
    !(W.tt[1] == P.t0 && W.tt[end] == P.t1) && error("Time axis mismatch between bridge P and driving W.") # not strictly an error
    solve!(Euler(), Y, P.v0, W, P)
    Y.yy[.., length(W.tt)] = P.v1
    Y
end

@deprecate bridge(W::SamplePath, P) solve(Euler(), W, P)
@deprecate bridge!(Y::SamplePath, W::SamplePath, P) solve!(Euler(), Y, W, P)


text:
# Julia wrapper for header: /usr/include/openssl/evp.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@c Int32 EVP_MD_type (Ptr{EVP_MD},) libcrypto
@c Int32 EVP_MD_pkey_type (Ptr{EVP_MD},) libcrypto
@c Int32 EVP_MD_size (Ptr{EVP_MD},) libcrypto
@c Int32 EVP_MD_block_size (Ptr{EVP_MD},) libcrypto
@c Uint32 EVP_MD_flags (Ptr{EVP_MD},) libcrypto
@c Ptr{EVP_MD} EVP_MD_CTX_md (Ptr{EVP_MD_CTX},) libcrypto
@c Int32 EVP_CIPHER_nid (Ptr{EVP_CIPHER},) libcrypto
@c Int32 EVP_CIPHER_block_size (Ptr{EVP_CIPHER},) libcrypto
@c Int32 EVP_CIPHER_key_length (Ptr{EVP_CIPHER},) libcrypto
@c Int32 EVP_CIPHER_iv_length (Ptr{EVP_CIPHER},) libcrypto
@c Uint32 EVP_CIPHER_flags (Ptr{EVP_CIPHER},) libcrypto
@c Ptr{EVP_CIPHER} EVP_CIPHER_CTX_cipher (Ptr{EVP_CIPHER_CTX},) libcrypto
@c Int32 EVP_CIPHER_CTX_nid (Ptr{EVP_CIPHER_CTX},) libcrypto
@c Int32 EVP_CIPHER_CTX_block_size (Ptr{EVP_CIPHER_CTX},) libcrypto
@c Int32 EVP_CIPHER_CTX_key_length (Ptr{EVP_CIPHER_CTX},) libcrypto
@c Int32 EVP_CIPHER_CTX_iv_length (Ptr{EVP_CIPHER_CTX},) libcrypto
@c Int32 EVP_CIPHER_CTX_copy (Ptr{EVP_CIPHER_CTX}, Ptr{EVP_CIPHER_CTX}) libcrypto
@c Ptr{None} EVP_CIPHER_CTX_get_app_data (Ptr{EVP_CIPHER_CTX},) libcrypto
@c None EVP_CIPHER_CTX_set_app_data (Ptr{EVP_CIPHER_CTX}, Ptr{None}) libcrypto
@c Uint32 EVP_CIPHER_CTX_flags (Ptr{EVP_CIPHER_CTX},) libcrypto
@c Int32 EVP_Cipher (Ptr{EVP_CIPHER_CTX}, Ptr{Uint8}, Ptr{Uint8}, Uint32) libcrypto
@c None EVP_MD_CTX_init (Ptr{EVP_MD_CTX},) libcrypto
@c Int32 EVP_MD_CTX_cleanup (Ptr{EVP_MD_CTX},) libcrypto
@c Ptr{EVP_MD_CTX} EVP_MD_CTX_create () libcrypto
@c None EVP_MD_CTX_destroy (Ptr{EVP_MD_CTX},) libcrypto
@c Int32 EVP_MD_CTX_copy_ex (Ptr{EVP_MD_CTX}, Ptr{EVP_MD_CTX}) libcrypto
@c None EVP_MD_CTX_set_flags (Ptr{EVP_MD_CTX}, Int32) libcrypto
@c None EVP_MD_CTX_clear_flags (Ptr{EVP_MD_CTX}, Int32) libcrypto
@c Int32 EVP_MD_CTX_test_flags (Ptr{EVP_MD_CTX}, Int32) libcrypto
@c Int32 EVP_DigestInit_ex (Ptr{EVP_MD_CTX}, Ptr{EVP_MD}, Ptr{ENGINE}) libcrypto
@c Int32 EVP_DigestUpdate (Ptr{EVP_MD_CTX}, Ptr{None}, size_t) libcrypto
@c Int32 EVP_DigestFinal_ex (Ptr{EVP_MD_CTX}, Ptr{Uint8}, Ptr{Uint32}) libcrypto
@c Int32 EVP_Digest (Ptr{None}, size_t, Ptr{Uint8}, Ptr{Uint32}, Ptr{EVP_MD}, Ptr{ENGINE}) libcrypto
@c Int32 EVP_MD_CTX_copy (Ptr{EVP_MD_CTX}, Ptr{EVP_MD_CTX}) libcrypto
@c Int32 EVP_DigestInit (Ptr{EVP_MD_CTX}, Ptr{EVP_MD}) libcrypto
@c Int32 EVP_DigestFinal (Ptr{EVP_MD_CTX}, Ptr{Uint8}, Ptr{Uint32}) libcrypto
@c Int32 EVP_read_pw_string (Ptr{Uint8}, Int32, Ptr{Uint8}, Int32) libcrypto
@c Int32 EVP_read_pw_string_min (Ptr{Uint8}, Int32, Int32, Ptr{Uint8}, Int32) libcrypto
@c None EVP_set_pw_prompt (Ptr{Uint8},) libcrypto
@c Ptr{Uint8} EVP_get_pw_prompt () libcrypto
@c Int32 EVP_BytesToKey (Ptr{EVP_CIPHER}, Ptr{EVP_MD}, Ptr{Uint8}, Ptr{Uint8}, Int32, Int32, Ptr{Uint8}, Ptr{Uint8}) libcrypto
@c None EVP_CIPHER_CTX_set_flags (Ptr{EVP_CIPHER_CTX}, Int32) libcrypto
@c None EVP_CIPHER_CTX_clear_flags (Ptr{EVP_CIPHER_CTX}, Int32) libcrypto
@c Int32 EVP_CIPHER_CTX_test_flags (Ptr{EVP_CIPHER_CTX}, Int32) libcrypto
@c Int32 EVP_EncryptInit (Ptr{EVP_CIPHER_CTX}, Ptr{EVP_CIPHER}, Ptr{Uint8}, Ptr{Uint8}) libcrypto
@c Int32 EVP_EncryptInit_ex (Ptr{EVP_CIPHER_CTX}, Ptr{EVP_CIPHER}, Ptr{ENGINE}, Ptr{Uint8}, Ptr{Uint8}) libcrypto
@c Int32 EVP_EncryptUpdate (Ptr{EVP_CIPHER_CTX}, Ptr{Uint8}, Ptr{Int32}, Ptr{Uint8}, Int32) libcrypto
@c Int32 EVP_EncryptFinal_ex (Ptr{EVP_CIPHER_CTX}, Ptr{Uint8}, Ptr{Int32}) libcrypto
@c Int32 EVP_EncryptFinal (Ptr{EVP_CIPHER_CTX}, Ptr{Uint8}, Ptr{Int32}) libcrypto
@c Int32 EVP_DecryptInit (Ptr{EVP_CIPHER_CTX}, Ptr{EVP_CIPHER}, Ptr{Uint8}, Ptr{Uint8}) libcrypto
@c Int32 EVP_DecryptInit_ex (Ptr{EVP_CIPHER_CTX}, Ptr{EVP_CIPHER}, Ptr{ENGINE}, Ptr{Uint8}, Ptr{Uint8}) libcrypto
@c Int32 EVP_DecryptUpdate (Ptr{EVP_CIPHER_CTX}, Ptr{Uint8}, Ptr{Int32}, Ptr{Uint8}, Int32) libcrypto
@c Int32 EVP_DecryptFinal (Ptr{EVP_CIPHER_CTX}, Ptr{Uint8}, Ptr{Int32}) libcrypto
@c Int32 EVP_DecryptFinal_ex (Ptr{EVP_CIPHER_CTX}, Ptr{Uint8}, Ptr{Int32}) libcrypto
@c Int32 EVP_CipherInit (Ptr{EVP_CIPHER_CTX}, Ptr{EVP_CIPHER}, Ptr{Uint8}, Ptr{Uint8}, Int32) libcrypto
@c Int32 EVP_CipherInit_ex (Ptr{EVP_CIPHER_CTX}, Ptr{EVP_CIPHER}, Ptr{ENGINE}, Ptr{Uint8}, Ptr{Uint8}, Int32) libcrypto
@c Int32 EVP_CipherUpdate (Ptr{EVP_CIPHER_CTX}, Ptr{Uint8}, Ptr{Int32}, Ptr{Uint8}, Int32) libcrypto
@c Int32 EVP_CipherFinal (Ptr{EVP_CIPHER_CTX}, Ptr{Uint8}, Ptr{Int32}) libcrypto
@c Int32 EVP_CipherFinal_ex (Ptr{EVP_CIPHER_CTX}, Ptr{Uint8}, Ptr{Int32}) libcrypto
@c Int32 EVP_SignFinal (Ptr{EVP_MD_CTX}, Ptr{Uint8}, Ptr{Uint32}, Ptr{EVP_PKEY}) libcrypto
@c Int32 EVP_VerifyFinal (Ptr{EVP_MD_CTX}, Ptr{Uint8}, Uint32, Ptr{EVP_PKEY}) libcrypto
@c Int32 EVP_DigestSignInit (Ptr{EVP_MD_CTX}, Ptr{Ptr{EVP_PKEY_CTX}}, Ptr{EVP_MD}, Ptr{ENGINE}, Ptr{EVP_PKEY}) libcrypto
@c Int32 EVP_DigestSignFinal (Ptr{EVP_MD_CTX}, Ptr{Uint8}, Ptr{size_t}) libcrypto
@c Int32 EVP_DigestVerifyInit (Ptr{EVP_MD_CTX}, Ptr{Ptr{EVP_PKEY_CTX}}, Ptr{EVP_MD}, Ptr{ENGINE}, Ptr{EVP_PKEY}) libcrypto
@c Int32 EVP_DigestVerifyFinal (Ptr{EVP_MD_CTX}, Ptr{Uint8}, size_t) libcrypto
@c Int32 EVP_OpenInit (Ptr{EVP_CIPHER_CTX}, Ptr{EVP_CIPHER}, Ptr{Uint8}, Int32, Ptr{Uint8}, Ptr{EVP_PKEY}) libcrypto
@c Int32 EVP_OpenFinal (Ptr{EVP_CIPHER_CTX}, Ptr{Uint8}, Ptr{Int32}) libcrypto
@c Int32 EVP_SealInit (Ptr{EVP_CIPHER_CTX}, Ptr{EVP_CIPHER}, Ptr{Ptr{Uint8}}, Ptr{Int32}, Ptr{Uint8}, Ptr{Ptr{EVP_PKEY}}, Int32) libcrypto
@c Int32 EVP_SealFinal (Ptr{EVP_CIPHER_CTX}, Ptr{Uint8}, Ptr{Int32}) libcrypto
@c None EVP_EncodeInit (Ptr{EVP_ENCODE_CTX},) libcrypto
@c None EVP_EncodeUpdate (Ptr{EVP_ENCODE_CTX}, Ptr{Uint8}, Ptr{Int32}, Ptr{Uint8}, Int32) libcrypto
@c None EVP_EncodeFinal (Ptr{EVP_ENCODE_CTX}, Ptr{Uint8}, Ptr{Int32}) libcrypto
@c Int32 EVP_EncodeBlock (Ptr{Uint8}, Ptr{Uint8}, Int32) libcrypto
@c None EVP_DecodeInit (Ptr{EVP_ENCODE_CTX},) libcrypto
@c Int32 EVP_DecodeUpdate (Ptr{EVP_ENCODE_CTX}, Ptr{Uint8}, Ptr{Int32}, Ptr{Uint8}, Int32) libcrypto
@c Int32 EVP_DecodeFinal (Ptr{EVP_ENCODE_CTX}, Ptr{Uint8}, Ptr{Int32}) libcrypto
@c Int32 EVP_DecodeBlock (Ptr{Uint8}, Ptr{Uint8}, Int32) libcrypto
@c None EVP_CIPHER_CTX_init (Ptr{EVP_CIPHER_CTX},) libcrypto
@c Int32 EVP_CIPHER_CTX_cleanup (Ptr{EVP_CIPHER_CTX},) libcrypto
@c Ptr{EVP_CIPHER_CTX} EVP_CIPHER_CTX_new () libcrypto
@c None EVP_CIPHER_CTX_free (Ptr{EVP_CIPHER_CTX},) libcrypto
@c Int32 EVP_CIPHER_CTX_set_key_length (Ptr{EVP_CIPHER_CTX}, Int32) libcrypto
@c Int32 EVP_CIPHER_CTX_set_padding (Ptr{EVP_CIPHER_CTX}, Int32) libcrypto
@c Int32 EVP_CIPHER_CTX_ctrl (Ptr{EVP_CIPHER_CTX}, Int32, Int32, Ptr{None}) libcrypto
@c Int32 EVP_CIPHER_CTX_rand_key (Ptr{EVP_CIPHER_CTX}, Ptr{Uint8}) libcrypto
@c Ptr{BIO_METHOD} BIO_f_md () libcrypto
@c Ptr{BIO_METHOD} BIO_f_base64 () libcrypto
@c Ptr{BIO_METHOD} BIO_f_cipher () libcrypto
@c Ptr{BIO_METHOD} BIO_f_reliable () libcrypto
@c None BIO_set_cipher (Ptr{BIO}, Ptr{EVP_CIPHER}, Ptr{Uint8}, Ptr{Uint8}, Int32) libcrypto
@c Ptr{EVP_MD} EVP_md_null () libcrypto
@c Ptr{EVP_MD} EVP_md4 () libcrypto
@c Ptr{EVP_MD} EVP_md5 () libcrypto
@c Ptr{EVP_MD} EVP_sha () libcrypto
@c Ptr{EVP_MD} EVP_sha1 () libcrypto
@c Ptr{EVP_MD} EVP_dss () libcrypto
@c Ptr{EVP_MD} EVP_dss1 () libcrypto
@c Ptr{EVP_MD} EVP_ecdsa () libcrypto
@c Ptr{EVP_MD} EVP_sha224 () libcrypto
@c Ptr{EVP_MD} EVP_sha256 () libcrypto
@c Ptr{EVP_MD} EVP_sha384 () libcrypto
@c Ptr{EVP_MD} EVP_sha512 () libcrypto
@c Ptr{EVP_MD} EVP_ripemd160 () libcrypto
@c Ptr{EVP_MD} EVP_whirlpool () libcrypto
@c Ptr{EVP_CIPHER} EVP_enc_null () libcrypto
@c Ptr{EVP_CIPHER} EVP_des_ecb () libcrypto
@c Ptr{EVP_CIPHER} EVP_des_ede () libcrypto
@c Ptr{EVP_CIPHER} EVP_des_ede3 () libcrypto
@c Ptr{EVP_CIPHER} EVP_des_ede_ecb () libcrypto
@c Ptr{EVP_CIPHER} EVP_des_ede3_ecb () libcrypto
@c Ptr{EVP_CIPHER} EVP_des_cfb64 () libcrypto
@c Ptr{EVP_CIPHER} EVP_des_cfb1 () libcrypto
@c Ptr{EVP_CIPHER} EVP_des_cfb8 () libcrypto
@c Ptr{EVP_CIPHER} EVP_des_ede_cfb64 () libcrypto
@c Ptr{EVP_CIPHER} EVP_des_ede3_cfb64 () libcrypto
@c Ptr{EVP_CIPHER} EVP_des_ede3_cfb1 () libcrypto
@c Ptr{EVP_CIPHER} EVP_des_ede3_cfb8 () libcrypto
@c Ptr{EVP_CIPHER} EVP_des_ofb () libcrypto
@c Ptr{EVP_CIPHER} EVP_des_ede_ofb () libcrypto
@c Ptr{EVP_CIPHER} EVP_des_ede3_ofb () libcrypto
@c Ptr{EVP_CIPHER} EVP_des_cbc () libcrypto
@c Ptr{EVP_CIPHER} EVP_des_ede_cbc () libcrypto
@c Ptr{EVP_CIPHER} EVP_des_ede3_cbc () libcrypto
@c Ptr{EVP_CIPHER} EVP_desx_cbc () libcrypto
@c Ptr{EVP_CIPHER} EVP_rc4 () libcrypto
@c Ptr{EVP_CIPHER} EVP_rc4_40 () libcrypto
@c Ptr{EVP_CIPHER} EVP_rc4_hmac_md5 () libcrypto
@c Ptr{EVP_CIPHER} EVP_rc2_ecb () libcrypto
@c Ptr{EVP_CIPHER} EVP_rc2_cbc () libcrypto
@c Ptr{EVP_CIPHER} EVP_rc2_40_cbc () libcrypto
@c Ptr{EVP_CIPHER} EVP_rc2_64_cbc () libcrypto
@c Ptr{EVP_CIPHER} EVP_rc2_cfb64 () libcrypto
@c Ptr{EVP_CIPHER} EVP_rc2_ofb () libcrypto
@c Ptr{EVP_CIPHER} EVP_bf_ecb () libcrypto
@c Ptr{EVP_CIPHER} EVP_bf_cbc () libcrypto
@c Ptr{EVP_CIPHER} EVP_bf_cfb64 () libcrypto
@c Ptr{EVP_CIPHER} EVP_bf_ofb () libcrypto
@c Ptr{EVP_CIPHER} EVP_cast5_ecb () libcrypto
@c Ptr{EVP_CIPHER} EVP_cast5_cbc () libcrypto
@c Ptr{EVP_CIPHER} EVP_cast5_cfb64 () libcrypto
@c Ptr{EVP_CIPHER} EVP_cast5_ofb () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_128_ecb () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_128_cbc () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_128_cfb1 () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_128_cfb8 () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_128_cfb128 () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_128_ofb () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_128_ctr () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_128_gcm () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_128_ccm () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_128_xts () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_192_ecb () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_192_cbc () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_192_cfb1 () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_192_cfb8 () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_192_cfb128 () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_192_ofb () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_192_ctr () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_192_gcm () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_192_ccm () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_256_ecb () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_256_cbc () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_256_cfb1 () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_256_cfb8 () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_256_cfb128 () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_256_ofb () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_256_ctr () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_256_gcm () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_256_ccm () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_256_xts () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_128_cbc_hmac_sha1 () libcrypto
@c Ptr{EVP_CIPHER} EVP_aes_256_cbc_hmac_sha1 () libcrypto
@c Ptr{EVP_CIPHER} EVP_camellia_128_ecb () libcrypto
@c Ptr{EVP_CIPHER} EVP_camellia_128_cbc () libcrypto
@c Ptr{EVP_CIPHER} EVP_camellia_128_cfb1 () libcrypto
@c Ptr{EVP_CIPHER} EVP_camellia_128_cfb8 () libcrypto
@c Ptr{EVP_CIPHER} EVP_camellia_128_cfb128 () libcrypto
@c Ptr{EVP_CIPHER} EVP_camellia_128_ofb () libcrypto
@c Ptr{EVP_CIPHER} EVP_camellia_192_ecb () libcrypto
@c Ptr{EVP_CIPHER} EVP_camellia_192_cbc () libcrypto
@c Ptr{EVP_CIPHER} EVP_camellia_192_cfb1 () libcrypto
@c Ptr{EVP_CIPHER} EVP_camellia_192_cfb8 () libcrypto
@c Ptr{EVP_CIPHER} EVP_camellia_192_cfb128 () libcrypto
@c Ptr{EVP_CIPHER} EVP_camellia_192_ofb () libcrypto
@c Ptr{EVP_CIPHER} EVP_camellia_256_ecb () libcrypto
@c Ptr{EVP_CIPHER} EVP_camellia_256_cbc () libcrypto
@c Ptr{EVP_CIPHER} EVP_camellia_256_cfb1 () libcrypto
@c Ptr{EVP_CIPHER} EVP_camellia_256_cfb8 () libcrypto
@c Ptr{EVP_CIPHER} EVP_camellia_256_cfb128 () libcrypto
@c Ptr{EVP_CIPHER} EVP_camellia_256_ofb () libcrypto
@c Ptr{EVP_CIPHER} EVP_seed_ecb () libcrypto
@c Ptr{EVP_CIPHER} EVP_seed_cbc () libcrypto
@c Ptr{EVP_CIPHER} EVP_seed_cfb128 () libcrypto
@c Ptr{EVP_CIPHER} EVP_seed_ofb () libcrypto
@c None OPENSSL_add_all_algorithms_noconf () libcrypto
@c None OPENSSL_add_all_algorithms_conf () libcrypto
@c None OpenSSL_add_all_ciphers () libcrypto
@c None OpenSSL_add_all_digests () libcrypto
@c Int32 EVP_add_cipher (Ptr{EVP_CIPHER},) libcrypto
@c Int32 EVP_add_digest (Ptr{EVP_MD},) libcrypto
@c Ptr{EVP_CIPHER} EVP_get_cipherbyname (Ptr{Uint8},) libcrypto
@c Ptr{EVP_MD} EVP_get_digestbyname (Ptr{Uint8},) libcrypto
@c None EVP_cleanup () libcrypto
@c None EVP_CIPHER_do_all (Ptr{Void}, Ptr{None}) libcrypto
@c None EVP_CIPHER_do_all_sorted (Ptr{Void}, Ptr{None}) libcrypto
@c None EVP_MD_do_all (Ptr{Void}, Ptr{None}) libcrypto
@c None EVP_MD_do_all_sorted (Ptr{Void}, Ptr{None}) libcrypto
@c Int32 EVP_PKEY_decrypt_old (Ptr{Uint8}, Ptr{Uint8}, Int32, Ptr{EVP_PKEY}) libcrypto
@c Int32 EVP_PKEY_encrypt_old (Ptr{Uint8}, Ptr{Uint8}, Int32, Ptr{EVP_PKEY}) libcrypto
@c Int32 EVP_PKEY_type (Int32,) libcrypto
@c Int32 EVP_PKEY_id (Ptr{EVP_PKEY},) libcrypto
@c Int32 EVP_PKEY_base_id (Ptr{EVP_PKEY},) libcrypto
@c Int32 EVP_PKEY_bits (Ptr{EVP_PKEY},) libcrypto
@c Int32 EVP_PKEY_size (Ptr{EVP_PKEY},) libcrypto
@c Int32 EVP_PKEY_set_type (Ptr{EVP_PKEY}, Int32) libcrypto
@c Int32 EVP_PKEY_set_type_str (Ptr{EVP_PKEY}, Ptr{Uint8}, Int32) libcrypto
@c Int32 EVP_PKEY_assign (Ptr{EVP_PKEY}, Int32, Ptr{None}) libcrypto
@c Ptr{None} EVP_PKEY_get0 (Ptr{EVP_PKEY},) libcrypto
@c Int32 EVP_PKEY_set1_RSA (Ptr{EVP_PKEY}, Ptr{Void}) libcrypto
@c Ptr{Void} EVP_PKEY_get1_RSA (Ptr{EVP_PKEY},) libcrypto
@c Int32 EVP_PKEY_set1_DSA (Ptr{EVP_PKEY}, Ptr{Void}) libcrypto
@c Ptr{Void} EVP_PKEY_get1_DSA (Ptr{EVP_PKEY},) libcrypto
@c Int32 EVP_PKEY_set1_DH (Ptr{EVP_PKEY}, Ptr{Void}) libcrypto
@c Ptr{Void} EVP_PKEY_get1_DH (Ptr{EVP_PKEY},) libcrypto
@c Int32 EVP_PKEY_set1_EC_KEY (Ptr{EVP_PKEY}, Ptr{Void}) libcrypto
@c Ptr{Void} EVP_PKEY_get1_EC_KEY (Ptr{EVP_PKEY},) libcrypto
@c Ptr{EVP_PKEY} EVP_PKEY_new () libcrypto
@c None EVP_PKEY_free (Ptr{EVP_PKEY},) libcrypto
@c Ptr{EVP_PKEY} d2i_PublicKey (Int32, Ptr{Ptr{EVP_PKEY}}, Ptr{Ptr{Uint8}}, Int32) libcrypto
@c Int32 i2d_PublicKey (Ptr{EVP_PKEY}, Ptr{Ptr{Uint8}}) libcrypto
@c Ptr{EVP_PKEY} d2i_PrivateKey (Int32, Ptr{Ptr{EVP_PKEY}}, Ptr{Ptr{Uint8}}, Int32) libcrypto
@c Ptr{EVP_PKEY} d2i_AutoPrivateKey (Ptr{Ptr{EVP_PKEY}}, Ptr{Ptr{Uint8}}, Int32) libcrypto
@c Int32 i2d_PrivateKey (Ptr{EVP_PKEY}, Ptr{Ptr{Uint8}}) libcrypto
@c Int32 EVP_PKEY_copy_parameters (Ptr{EVP_PKEY}, Ptr{EVP_PKEY}) libcrypto
@c Int32 EVP_PKEY_missing_parameters (Ptr{EVP_PKEY},) libcrypto
@c Int32 EVP_PKEY_save_parameters (Ptr{EVP_PKEY}, Int32) libcrypto
@c Int32 EVP_PKEY_cmp_parameters (Ptr{EVP_PKEY}, Ptr{EVP_PKEY}) libcrypto
@c Int32 EVP_PKEY_cmp (Ptr{EVP_PKEY}, Ptr{EVP_PKEY}) libcrypto
@c Int32 EVP_PKEY_print_public (Ptr{BIO}, Ptr{EVP_PKEY}, Int32, Ptr{ASN1_PCTX}) libcrypto
@c Int32 EVP_PKEY_print_private (Ptr{BIO}, Ptr{EVP_PKEY}, Int32, Ptr{ASN1_PCTX}) libcrypto
@c Int32 EVP_PKEY_print_params (Ptr{BIO}, Ptr{EVP_PKEY}, Int32, Ptr{ASN1_PCTX}) libcrypto
@c Int32 EVP_PKEY_get_default_digest_nid (Ptr{EVP_PKEY}, Ptr{Int32}) libcrypto
@c Int32 EVP_CIPHER_type (Ptr{EVP_CIPHER},) libcrypto
@c Int32 EVP_CIPHER_param_to_asn1 (Ptr{EVP_CIPHER_CTX}, Ptr{ASN1_TYPE}) libcrypto
@c Int32 EVP_CIPHER_asn1_to_param (Ptr{EVP_CIPHER_CTX}, Ptr{ASN1_TYPE}) libcrypto
@c Int32 EVP_CIPHER_set_asn1_iv (Ptr{EVP_CIPHER_CTX}, Ptr{ASN1_TYPE}) libcrypto
@c Int32 EVP_CIPHER_get_asn1_iv (Ptr{EVP_CIPHER_CTX}, Ptr{ASN1_TYPE}) libcrypto
@c Int32 PKCS5_PBE_keyivgen (Ptr{EVP_CIPHER_CTX}, Ptr{Uint8}, Int32, Ptr{ASN1_TYPE}, Ptr{EVP_CIPHER}, Ptr{EVP_MD}, Int32) libcrypto
@c Int32 PKCS5_PBKDF2_HMAC_SHA1 (Ptr{Uint8}, Int32, Ptr{Uint8}, Int32, Int32, Int32, Ptr{Uint8}) libcrypto
@c Int32 PKCS5_PBKDF2_HMAC (Ptr{Uint8}, Int32, Ptr{Uint8}, Int32, Int32, Ptr{EVP_MD}, Int32, Ptr{Uint8}) libcrypto
@c Int32 PKCS5_v2_PBE_keyivgen (Ptr{EVP_CIPHER_CTX}, Ptr{Uint8}, Int32, Ptr{ASN1_TYPE}, Ptr{EVP_CIPHER}, Ptr{EVP_MD}, Int32) libcrypto
@c None PKCS5_PBE_add () libcrypto
@c Int32 EVP_PBE_CipherInit (Ptr{ASN1_OBJECT}, Ptr{Uint8}, Int32, Ptr{ASN1_TYPE}, Ptr{EVP_CIPHER_CTX}, Int32) libcrypto
@c Int32 EVP_PBE_alg_add_type (Int32, Int32, Int32, Int32, Ptr{EVP_PBE_KEYGEN}) libcrypto
@c Int32 EVP_PBE_alg_add (Int32, Ptr{EVP_CIPHER}, Ptr{EVP_MD}, Ptr{EVP_PBE_KEYGEN}) libcrypto
@c Int32 EVP_PBE_find (Int32, Int32, Ptr{Int32}, Ptr{Int32}, Ptr{Ptr{EVP_PBE_KEYGEN}}) libcrypto
@c None EVP_PBE_cleanup () libcrypto
@c Int32 EVP_PKEY_asn1_get_count () libcrypto
@c Ptr{EVP_PKEY_ASN1_METHOD} EVP_PKEY_asn1_get0 (Int32,) libcrypto
@c Ptr{EVP_PKEY_ASN1_METHOD} EVP_PKEY_asn1_find (Ptr{Ptr{ENGINE}}, Int32) libcrypto
@c Ptr{EVP_PKEY_ASN1_METHOD} EVP_PKEY_asn1_find_str (Ptr{Ptr{ENGINE}}, Ptr{Uint8}, Int32) libcrypto
@c Int32 EVP_PKEY_asn1_add0 (Ptr{EVP_PKEY_ASN1_METHOD},) libcrypto
@c Int32 EVP_PKEY_asn1_add_alias (Int32, Int32) libcrypto
@c Int32 EVP_PKEY_asn1_get0_info (Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Ptr{Uint8}}, Ptr{Ptr{Uint8}}, Ptr{EVP_PKEY_ASN1_METHOD}) libcrypto
@c Ptr{EVP_PKEY_ASN1_METHOD} EVP_PKEY_get0_asn1 (Ptr{EVP_PKEY},) libcrypto
@c Ptr{EVP_PKEY_ASN1_METHOD} EVP_PKEY_asn1_new (Int32, Int32, Ptr{Uint8}, Ptr{Uint8}) libcrypto
@c None EVP_PKEY_asn1_copy (Ptr{EVP_PKEY_ASN1_METHOD}, Ptr{EVP_PKEY_ASN1_METHOD}) libcrypto
@c None EVP_PKEY_asn1_free (Ptr{EVP_PKEY_ASN1_METHOD},) libcrypto
@c None EVP_PKEY_asn1_set_public (Ptr{EVP_PKEY_ASN1_METHOD}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}) libcrypto
@c None EVP_PKEY_asn1_set_private (Ptr{EVP_PKEY_ASN1_METHOD}, Ptr{Void}, Ptr{Void}, Ptr{Void}) libcrypto
@c None EVP_PKEY_asn1_set_param (Ptr{EVP_PKEY_ASN1_METHOD}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}) libcrypto
@c None EVP_PKEY_asn1_set_free (Ptr{EVP_PKEY_ASN1_METHOD}, Ptr{Void}) libcrypto
@c None EVP_PKEY_asn1_set_ctrl (Ptr{EVP_PKEY_ASN1_METHOD}, Ptr{Void}) libcrypto
@c Ptr{EVP_PKEY_METHOD} EVP_PKEY_meth_find (Int32,) libcrypto
@c Ptr{EVP_PKEY_METHOD} EVP_PKEY_meth_new (Int32, Int32) libcrypto
@c None EVP_PKEY_meth_get0_info (Ptr{Int32}, Ptr{Int32}, Ptr{EVP_PKEY_METHOD}) libcrypto
@c None EVP_PKEY_meth_copy (Ptr{EVP_PKEY_METHOD}, Ptr{EVP_PKEY_METHOD}) libcrypto
@c None EVP_PKEY_meth_free (Ptr{EVP_PKEY_METHOD},) libcrypto
@c Int32 EVP_PKEY_meth_add0 (Ptr{EVP_PKEY_METHOD},) libcrypto
@c Ptr{EVP_PKEY_CTX} EVP_PKEY_CTX_new (Ptr{EVP_PKEY}, Ptr{ENGINE}) libcrypto
@c Ptr{EVP_PKEY_CTX} EVP_PKEY_CTX_new_id (Int32, Ptr{ENGINE}) libcrypto
@c Ptr{EVP_PKEY_CTX} EVP_PKEY_CTX_dup (Ptr{EVP_PKEY_CTX},) libcrypto
@c None EVP_PKEY_CTX_free (Ptr{EVP_PKEY_CTX},) libcrypto
@c Int32 EVP_PKEY_CTX_ctrl (Ptr{EVP_PKEY_CTX}, Int32, Int32, Int32, Int32, Ptr{None}) libcrypto
@c Int32 EVP_PKEY_CTX_ctrl_str (Ptr{EVP_PKEY_CTX}, Ptr{Uint8}, Ptr{Uint8}) libcrypto
@c Int32 EVP_PKEY_CTX_get_operation (Ptr{EVP_PKEY_CTX},) libcrypto
@c None EVP_PKEY_CTX_set0_keygen_info (Ptr{EVP_PKEY_CTX}, Ptr{Int32}, Int32) libcrypto
@c Ptr{EVP_PKEY} EVP_PKEY_new_mac_key (Int32, Ptr{ENGINE}, Ptr{Uint8}, Int32) libcrypto
@c None EVP_PKEY_CTX_set_data (Ptr{EVP_PKEY_CTX}, Ptr{None}) libcrypto
@c Ptr{None} EVP_PKEY_CTX_get_data (Ptr{EVP_PKEY_CTX},) libcrypto
@c Ptr{EVP_PKEY} EVP_PKEY_CTX_get0_pkey (Ptr{EVP_PKEY_CTX},) libcrypto
@c Ptr{EVP_PKEY} EVP_PKEY_CTX_get0_peerkey (Ptr{EVP_PKEY_CTX},) libcrypto
@c None EVP_PKEY_CTX_set_app_data (Ptr{EVP_PKEY_CTX}, Ptr{None}) libcrypto
@c Ptr{None} EVP_PKEY_CTX_get_app_data (Ptr{EVP_PKEY_CTX},) libcrypto
@c Int32 EVP_PKEY_sign_init (Ptr{EVP_PKEY_CTX},) libcrypto
@c Int32 EVP_PKEY_sign (Ptr{EVP_PKEY_CTX}, Ptr{Uint8}, Ptr{size_t}, Ptr{Uint8}, size_t) libcrypto
@c Int32 EVP_PKEY_verify_init (Ptr{EVP_PKEY_CTX},) libcrypto
@c Int32 EVP_PKEY_verify (Ptr{EVP_PKEY_CTX}, Ptr{Uint8}, size_t, Ptr{Uint8}, size_t) libcrypto
@c Int32 EVP_PKEY_verify_recover_init (Ptr{EVP_PKEY_CTX},) libcrypto
@c Int32 EVP_PKEY_verify_recover (Ptr{EVP_PKEY_CTX}, Ptr{Uint8}, Ptr{size_t}, Ptr{Uint8}, size_t) libcrypto
@c Int32 EVP_PKEY_encrypt_init (Ptr{EVP_PKEY_CTX},) libcrypto
@c Int32 EVP_PKEY_encrypt (Ptr{EVP_PKEY_CTX}, Ptr{Uint8}, Ptr{size_t}, Ptr{Uint8}, size_t) libcrypto
@c Int32 EVP_PKEY_decrypt_init (Ptr{EVP_PKEY_CTX},) libcrypto
@c Int32 EVP_PKEY_decrypt (Ptr{EVP_PKEY_CTX}, Ptr{Uint8}, Ptr{size_t}, Ptr{Uint8}, size_t) libcrypto
@c Int32 EVP_PKEY_derive_init (Ptr{EVP_PKEY_CTX},) libcrypto
@c Int32 EVP_PKEY_derive_set_peer (Ptr{EVP_PKEY_CTX}, Ptr{EVP_PKEY}) libcrypto
@c Int32 EVP_PKEY_derive (Ptr{EVP_PKEY_CTX}, Ptr{Uint8}, Ptr{size_t}) libcrypto
@c Int32 EVP_PKEY_paramgen_init (Ptr{EVP_PKEY_CTX},) libcrypto
@c Int32 EVP_PKEY_paramgen (Ptr{EVP_PKEY_CTX}, Ptr{Ptr{EVP_PKEY}}) libcrypto
@c Int32 EVP_PKEY_keygen_init (Ptr{EVP_PKEY_CTX},) libcrypto
@c Int32 EVP_PKEY_keygen (Ptr{EVP_PKEY_CTX}, Ptr{Ptr{EVP_PKEY}}) libcrypto
@c None EVP_PKEY_CTX_set_cb (Ptr{EVP_PKEY_CTX}, Ptr{EVP_PKEY_gen_cb}) libcrypto
@c Ptr{EVP_PKEY_gen_cb} EVP_PKEY_CTX_get_cb (Ptr{EVP_PKEY_CTX},) libcrypto
@c Int32 EVP_PKEY_CTX_get_keygen_info (Ptr{EVP_PKEY_CTX}, Int32) libcrypto
@c None EVP_PKEY_meth_set_init (Ptr{EVP_PKEY_METHOD}, Ptr{Void}) libcrypto
@c None EVP_PKEY_meth_set_copy (Ptr{EVP_PKEY_METHOD}, Ptr{Void}) libcrypto
@c None EVP_PKEY_meth_set_cleanup (Ptr{EVP_PKEY_METHOD}, Ptr{Void}) libcrypto
@c None EVP_PKEY_meth_set_paramgen (Ptr{EVP_PKEY_METHOD}, Ptr{Void}, Ptr{Void}) libcrypto
@c None EVP_PKEY_meth_set_keygen (Ptr{EVP_PKEY_METHOD}, Ptr{Void}, Ptr{Void}) libcrypto
@c None EVP_PKEY_meth_set_sign (Ptr{EVP_PKEY_METHOD}, Ptr{Void}, Ptr{Void}) libcrypto
@c None EVP_PKEY_meth_set_verify (Ptr{EVP_PKEY_METHOD}, Ptr{Void}, Ptr{Void}) libcrypto
@c None EVP_PKEY_meth_set_verify_recover (Ptr{EVP_PKEY_METHOD}, Ptr{Void}, Ptr{Void}) libcrypto
@c None EVP_PKEY_meth_set_signctx (Ptr{EVP_PKEY_METHOD}, Ptr{Void}, Ptr{Void}) libcrypto
@c None EVP_PKEY_meth_set_verifyctx (Ptr{EVP_PKEY_METHOD}, Ptr{Void}, Ptr{Void}) libcrypto
@c None EVP_PKEY_meth_set_encrypt (Ptr{EVP_PKEY_METHOD}, Ptr{Void}, Ptr{Void}) libcrypto
@c None EVP_PKEY_meth_set_decrypt (Ptr{EVP_PKEY_METHOD}, Ptr{Void}, Ptr{Void}) libcrypto
@c None EVP_PKEY_meth_set_derive (Ptr{EVP_PKEY_METHOD}, Ptr{Void}, Ptr{Void}) libcrypto
@c None EVP_PKEY_meth_set_ctrl (Ptr{EVP_PKEY_METHOD}, Ptr{Void}, Ptr{Void}) libcrypto
@c None ERR_load_EVP_strings () libcrypto


text:
import StaticLint: hasfile, canloadfile, loadfile, setfile, getfile, getsymbols, getsymbolextendeds, getenv
import StaticLint: getpath, getroot, setroot, getcst, setcst, semantic_pass, getserver, setserver
hasfile(server::LanguageServerInstance, path::String) = !isempty(path) && hasdocument(server, URI2(filepath2uri(path)))
function canloadfile(server::LanguageServerInstance, path::String)
    try
        return !isempty(path) && !safe_isfile(path)
    catch err
        isa(err, Base.IOError) || isa(err, Base.SystemError) || rethrow()
        return false
    end
end
function loadfile(server::LanguageServerInstance, path::String)
    source = try
        s = read(path, String)
        isvalid(s) || return
        s
    catch err
        isa(err, Base.IOError) || isa(err, Base.SystemError) || rethrow()
        return
    end
    uri = filepath2uri(path)
    doc = Document(uri, source, true, server)
    StaticLint.setfile(server, path, doc)
end
function setfile(server::LanguageServerInstance, path::String, x::Document)
    uri = URI2(filepath2uri(path))
    if hasdocument(server, uri)
        error("StaticLint should not try to set documents that are already tracked.")
    end

    setdocument!(server, uri, x)
end
getfile(server::LanguageServerInstance, path::String) = getdocument(server, URI2(filepath2uri(path)))

function getenv(doc::Document, server::LanguageServerInstance)
    get(server.roots_env_map, doc.root, server.global_env)
end
getenv(doc::Document) = getenv(doc, doc.server)
getenv(server::LanguageServerInstance) = server.global_env

getpath(d::Document) = d._path

getroot(d::Document) = d.root
function setroot(doc::Document, root::Document)
    if isdefined(doc, :root) && doc == doc.root && root !== doc
        # doc is being unset as a root - remove ExternalEnv if there is one
        if doc.server isa LanguageServerInstance && haskey(doc.server.roots_env_map, doc)
            delete!(doc.server.roots_env_map, doc)
        end
    end
    doc.root = root
    if doc == root && doc.server isa LanguageServerInstance
        # doc is being set as it's own root, lets find
        extenv = get_env_for_root(doc, doc.server)
        if extenv !== nothing
            doc.server.roots_env_map[doc] = extenv
        end
    end
    return doc
end

getcst(d::Document) = d.cst
function setcst(d::Document, cst::EXPR)
    d.cst = cst
    return d
end

getserver(file::Document) = file.server
function setserver(file::Document, server::LanguageServerInstance)
    file.server = server
    return file
end

function lint!(doc, server)
    StaticLint.check_all(getcst(doc), server.lint_options, getenv(doc, server))
    empty!(doc.diagnostics)
    mark_errors(doc, doc.diagnostics)
    # TODO Ideally we would not want to acces jr_endpoint here
    publish_diagnostics(doc, server, server.jr_endpoint)
end


text:
@testset "Nextfloat" begin
    @testset for T in (Posit8, Posit16, Posit32,

                    # Doesn't work yet because -PositX1 requires pX1_mul
                    # Posit8_1, Posit16_1, Posit24_1,
                    Posit8_2, Posit16_2, Posit24_2)

        @test -floatmax(T) == nextfloat(nextfloat(floatmax(T)))
        @test -floatmin(T) == prevfloat(prevfloat(floatmin(T)))
    end
end


text:
"""
Abstract supertype for all types describing the MTG coding
for a node.

See [`NodeMTG`](@ref) and [`MutableNodeMTG`](@ref) for examples
of implementation.
"""
abstract type AbstractNodeMTG end

"""
    NodeMTG(link, symbol, index, scale)
    MutableNodeMTG(link, symbol, index, scale)

# NodeMTG structure

Builds an MTG node to hold data about the link to the previous node,
the symbol of the node, and its index.

# Note

- The symbol should match the possible values listed in the `SYMBOL` column of the `CLASSES` section
 in the mtg file if read from a file.

- The index is totaly free, and can be used as a way to *e.g.* keep track of the branching order.

```julia
NodeMTG("<", "Leaf", 2, 0)
```
"""
NodeMTG, MutableNodeMTG

struct NodeMTG <: AbstractNodeMTG
    link::String
    symbol::Union{String,SubString,Char}
    index::Union{Int,Nothing}
    scale::Int
end

mutable struct MutableNodeMTG <: AbstractNodeMTG
    link::String
    symbol::Union{String,SubString,Char}
    index::Union{Int,Nothing}
    scale::Int
end


text:
# Autogenerated wrapper script for OpenSSL_jll for x86_64-apple-darwin
export libcrypto, libssl, openssl

JLLWrappers.@generate_wrapper_header("OpenSSL")
JLLWrappers.@declare_library_product(libcrypto, "@rpath/libcrypto.1.1.dylib")
JLLWrappers.@declare_library_product(libssl, "@rpath/libssl.1.1.dylib")
JLLWrappers.@declare_executable_product(openssl)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcrypto,
        "lib/libcrypto.1.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libssl,
        "lib/libssl.1.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        openssl,
        "bin/openssl",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
# PSFModels

Statistical models for constructing point-spread functions (PSFs). These models act like matrices but without allocating any memory, which makes them efficient to fit and apply.

## Models

The following models are currently implemented
* [`PSFModels.Gaussian`](@ref)/[`PSFModels.Normal`](@ref)
* [`PSFModels.AiryDisk`](@ref)
* [`PSFModels.Moffat`](@ref)

## Parameters

In general, the PSFs have a position, a full-width at half-maximum (FWHM) measure, and an amplitude. The position follows a 1-based pixel coordinate system, where `(1, 1)` represents the *center* of the bottom left pixel. This matches the indexing style of Julia as well as DS9, IRAF, SourceExtractor, and WCS. If a position is not specified, it is set to `(0, 0)`. The FWHM is a consistent scale parameter for the models. All models support a scalar (isotropic) FWHM and a FWHM for each axis (diagonal).

## Usage

Using the models should feel just like an array. In fact, `PSFModels.PSFModel <: AbstractMatrix`. However, no data is stored and no allocations have to be made. In other words, representing the models as matrices is merely a convenience, since typically astronomical data is stored in dense arrays. Another way of thinking of these is a lazy array that applies a function when indexed, rather than returning data stored in memory.

```jldoctest model
julia> m = PSFModels.Gaussian(fwhm=3); # centered at (0, 0)

julia> m[0, 0] # [y, x] for indexing
1.0

julia> m(0, 0) # (x, y) for evaluating
1.0
```

!!! note "axis order"

    It's important to note the difference in the axis ordering between the index-style calls and the function-style calls. The index-style calls are reverse cartesian order (e.g., `(z, y, x)`), while function calls are the typical cartesian order `(x, y, z)`. Regardless, the constructors are always in cartesian order (`(x, y, z)`).

Because the model is a matrix, it needs to have a size. Each model has a bounding box which can be controlled with the `extent` keyword. By default the extent is set by a scalar factor of the FWHM (e.g., `maxsize * FWHM` pixels), centered around the PSF, and rounded up. We can see how this alters the indices from a typical `Matrix`

```jldoctest model
julia> size(m) # default 'stamp' size is fwhm * 3
(11, 11)

julia> axes(m)
(-5:5, -5:5)
```

if we want to collect the model into a dense matrix, regardless of the indexing (e.g. to prepare for cross-correlation), we can simply

```jldoctest model
julia> stamp = collect(m);

```

these axes are merely a convenience for bounding the model, since they accept any real number as input.

```jldoctest model
julia> m[100, 10000] # index-like inputs [y, x]
0.0

julia> m(2.4, 1.7) # valid for any real (x, y)
0.0696156536973086
```

By bounding the model, we get a cutout which can be applied to arrays with much larger dimensions without having to iterate over the whole matrix

```jldoctest
julia> big_mat = ones(1001, 1001);

julia> model = PSFModels.Gaussian(x=51, y=51, fwhm=2);

julia> ax = map(intersect, axes(big_mat), axes(model))
(48:54, 48:54)

julia> cutout = @view big_mat[ax...]
7×7 view(::Matrix{Float64}, 48:54, 48:54) with eltype Float64:
 1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0

julia> stamp = @view model[ax...];

julia> photsum = sum(cutout .* stamp)
4.5322418212890625
```
Nice- we only had to reduce ~50 pixels instead of ~1,000,000 to calculate the aperture sum, all in under a microsecond (on my machine).

Since the models are lazy, that means the type of the output can be specified, as long as it can be converted to from a real number (so no integer types).

```jldoctest
julia> mbig = PSFModels.Gaussian(BigFloat, fwhm=12);

julia> sum(mbig)
163.07467408408593885562554918859656805096847165259532630443572998046875
```

finally, we provide plotting recipes from [RecipesBase.jl](https://github.com/JuliaPlots/RecipesBase.jl), which can be seen in use in the [API/Reference](@ref) section.

```julia
using Plots
model = PSFModels.Gaussian(8)
plot(model)              # default axes
plot(model, :, 1:5)    # custom axes (y, x)
plot(model, axes(other)) # use axes from other array
```
"""
module PSFModels

using CoordinateTransformations
using Distances
using KeywordCalls
using SpecialFunctions
using StaticArrays

include("core.jl")
include("gaussian.jl")
include("moffat.jl")
include("airy.jl")
include("plotting.jl")

end # module PSFModels


text:
module VlasiatorMakie

using Vlasiator, StatsBase, Printf

using Makie

include("typerecipe.jl")
include("fullrecipe.jl")
include("interactive.jl")

export
   vlheatmap, vlslice, vlslices,
   vdfvolume, vdfslice, vdfslices

end


text:
using DrWatson
@quickactivate "PECUZAL_Julia"

# Core functions for the usage in `/scripts`-folder
using DelayEmbeddings
using DynamicalSystemsBase
using Random
using Neighborhood
using StatsBase
using Distances

export pecuzal_embedding_update

"""
    pecuzal_embedding_update(s; kwargs...) → Y, τ_vals, ts_vals, ΔLs ,⟨ε★⟩
A unified approach to properly embed a time series or a set of time series
(`Dataset`) based on the ideas of Pecora et al. [^Pecoral2007] and Uzal et al.
[^Uzal2011].

## Keyword arguments

* `τs = 0:50`: Possible delay values `τs` (in sampling time units). For each of
  the `τs`'s the continuity statistic ⟨ε★⟩ gets computed and further processed
  in order to find optimal delays `τᵢ` for each embedding cycle `i` (read
  algorithm description).
* `w::Int = 1`: Theiler window (neighbors in time with index `w` close to the point,
  that are excluded from being true neighbors). `w=0` means to exclude only the
  point itself, and no temporal neighbors.
* `samplesize::Real = 1`: determine the fraction of all phase space points
  (=`length(s)`) to be considered (fiducial points v) to average ε★, in order to
  produce `⟨ε★⟩`.
* `K::Int = 13`: the amount of nearest neighbors in the δ-ball (read algorithm
  description). Must be at least 8 (in order to gurantee a valid statistic).
  `⟨ε★⟩` is computed taking the minimum result over all `k ∈ K`.
* `KNN::Int = 3`: the amount of nearest neighbors considered, in order to compute
  σ_k^2 (read algorithm description [`uzal_cost`]@ref). If given a vector, the
  minimum result over all `knn ∈ KNN` is returned.
* `L_threshold::Real = 0`: The algorithm breaks, when this threshold is exceeded
  by `ΔL` in an embedding cycle (set as a positive number, i.e. an absolute value
  of `ΔL`).
* `α::Real = 0.05`: The significance level for obtaining the continuity statistic
* `p::Real = 0.5`: The p-parameter for the binomial distribution used for the
  computation of the continuity statistic ⟨ε★⟩.
* `max_cycles = 50`: The algorithm will stop after that many cycles no matter what.
* `econ::Bool = false`: Economy-mode for L-statistic computation. Instead of
  computing L-statistics for time horizons `2:Tw`, here we only compute them for
  `2:2:Tw`, see description for further details.

## Description
The method works iteratively and gradually builds the final embedding vectors
`Y`. Based on the `⟨ε★⟩`-statistic [`pecora`](@ref) the algorithm picks an
optimal delay value `τᵢ` for each embedding cycle `i`.
For achieving that, we take the inpute time series `s`, denoted as the actual
phase space trajectory `Y_actual` and compute the continuity statistic `⟨ε★⟩`.
1. Each local maxima in `⟨ε★⟩` is used for constructing a
candidate embedding trajectory `Y_trial` with a delay corresponding to that
specific peak in `⟨ε★⟩`. 2. We then compute the `L`-statistic [`uzal_cost`](@ref)
for `Y_trial` (`L-trial`) and `Y_actual` (`L_actual`) for increasing prediction
time horizons (free parameter in the `L`-statistic) and save the maximum
difference `max(L-trial - L_actual)` as `ΔL` (Note that this is a
negative number, since the `L`-statistic decreases with better reconstructions).
3. We pick the peak/`τ`-value, for which `ΔL` is minimal (=maximum decrease of
the overall `L`-value) and construct the actual embedding trajectory
`Y_actual` (steps 1.-3. correspond to an embedding cycle). 4. We repeat steps
1.-3. with `Y_actual` as input and stop the algorithm when `ΔL` is > 0,
i.e. when and additional embedding component would not lead to a lower overall
L-value. `Y_actual` -> `Y`.

In case of multivariate embedding, i.e. when embedding a set of M time series
(`s::Dataset`), in each embedding cycle the continuity statistic `⟨ε★⟩` gets
computed for all M time series available. The optimal delay value `τ` in each
embedding cycle is chosen as the peak/`τ`-value for which `ΔL` is
minimal under all available peaks and under all M `⟨ε★⟩`'s. In the first
embedding cycle there will be M² different `⟨ε★⟩`'s to consider, since it is
not clear a priori which time series of the input should consitute the first
component of the embedding vector and form `Y_actual`.

The range of considered delay values is determined in `τs` and for the
nearest neighbor search we respect the Theiler window `w`. The final embedding
vector is stored in `Y` (`Dataset`). The chosen delay values for each embedding
cycle are stored in `τ_vals` and the according time series numbers chosen for
each delay value in `τ_vals` are stored in `ts_vals`. For univariate embedding
(`s::Vector`) `ts_vals` is a vector of ones of length `τ_vals`, because there is
simply just one time series to choose from. The function also returns the
`ΔLs`-values for each embedding cycle and the continuity statistic `⟨ε★⟩`
as an `Array` of `Vector`s.

For distance computations the Euclidean norm is used.

[^Pecora2007]: Pecora, L. M., Moniz, L., Nichols, J., & Carroll, T. L. (2007). [A unified approach to attractor reconstruction. Chaos 17(1)](https://doi.org/10.1063/1.2430294).
[^Uzal2011]: Uzal, L. C., Grinblat, G. L., Verdes, P. F. (2011). [Optimal reconstruction of dynamical systems: A noise amplification approach. Physical Review E 84, 016223](https://doi.org/10.1103/PhysRevE.84.016223).
"""
function pecuzal_embedding_update(s::Vector{T}; τs = 0:50 , w::Int = 1,
    samplesize::Real = 1, K::Int = 13, KNN::Int = 3, L_threshold::Real = 0,
    α::Real = 0.05, p::Real = 0.5, max_cycles::Int = 50, econ::Bool = false
    ) where {T<:Real}

    @assert 0 < samplesize ≤ 1 "Please select a valid `samplesize`, which denotes a fraction of considered fiducial points, i.e. `samplesize` ∈ (0 1]"
    @assert all(x -> x ≥ 0, τs)
    @assert L_threshold ≥ 0
    threshold = -L_threshold # due to the negativity of L-decrease
    metric = Euclidean()

    s_orig = s
    s = regularize(s) # especially important for comparative L-statistics
    # define actual phase space trajectory
    Y_act = Dataset(s)

    # set a flag, in order to tell the while loop when to stop. Each loop
    # stands for encountering a new embedding dimension
    flag, counter = true, 1

    # preallocate output variables
    τ_vals = Int64[0]
    ts_vals = Int64[1]
    Ls = Float64[]
    ε★s = Array{T}(undef, length(τs), max_cycles)

    # loop over increasing embedding dimensions until some break criterion will
    # tell the loop to stop/break
    while flag
        Y_act = pecuzal_embedding_cycle!(
                Y_act, flag, s, τs, w, counter, ε★s, τ_vals, metric,
                Ls, ts_vals, samplesize, K, α, p, KNN, econ)
        flag = pecuzal_break_criterion(Ls, counter, max_cycles, threshold)
        counter += 1
    end
    # construct final reconstruction vector
    NN = (length(s)-sum(τ_vals[1:end-1]))
    Y_final = s_orig
    for i = 2:length(τ_vals[1:end-1])
        Y_final = DelayEmbeddings.hcat_lagged_values(Y_final, s_orig, τ_vals[i])
    end
    return Y_final, τ_vals[1:end-1], ts_vals[1:end-1], Ls[1:end-1], ε★s[:,1:counter-1]

end

function pecuzal_embedding_update(Y::Dataset{D, T}; τs = 0:50 , w::Int = 1,
    samplesize::Real = 1, K::Int = 13, KNN::Int = 3, L_threshold::Real = 0,
    α::Real = 0.05, p::Real = 0.5, max_cycles::Int = 50, econ::Bool = false
    ) where {D, T<:Real}

    @assert 0 < samplesize ≤ 1 "Please select a valid `samplesize`, which denotes a fraction of considered fiducial points, i.e. `samplesize` ∈ (0 1]"
    @assert all(x -> x ≥ 0, τs)
    @assert L_threshold ≥ 0
    threshold = -L_threshold # due to the negativity of L-decrease
    metric = Euclidean()

    Y_orig = Y
    Y = regularize(Y) # especially important for comparative L-statistics

    # define actual phase space trajectory
    Y_act = []

    # set a flag, in order to tell the while loop when to stop. Each loop
    # stands for encountering a new embedding dimension
    flag, counter = true, 1

    # preallocate output variables
    τ_vals = Int64[0]
    ts_vals = Int64[]
    Ls = Float64[]
    ε★s = fill(zeros(T, length(τs), size(Y,2)), 1, max_cycles)

    # loop over increasing embedding dimensions until some break criterion will
    # tell the loop to stop/break
    while flag
        Y_act = pecuzal_multivariate_embedding_cycle!(
                Y_act, flag, Y, τs, w, counter, ε★s, τ_vals, metric,
                Ls, ts_vals, samplesize, K, α, p, KNN, econ)

        flag = pecuzal_break_criterion(Ls, counter, max_cycles, threshold)
        counter += 1
    end
    # construct final reconstruction vector
    Y_final = Y_orig[:,ts_vals[1]]
    for i = 2:length(τ_vals[1:counter-1])
        Y_final = DelayEmbeddings.hcat_lagged_values(Y_final,Y_orig[:,ts_vals[i]],τ_vals[i])
    end

    return Y_final, τ_vals[1:end-1], ts_vals[1:end-1], Ls[1:end-1], ε★s[:,1:counter-1]

end

"""
Perform one univariate embedding cycle on `Y_act`. Return the new `Y_act`
"""
function pecuzal_embedding_cycle!(Y_act::Dataset{D, T}, flag::Bool, s::Vector,
        τs, w::Int, counter::Int, ε★s::AbstractArray, τ_vals::Vector{Int}, metric,
        Ls::Vector{T}, ts_vals::Vector{Int}, samplesize::Real, K::Int, α::Real,
        p::Real, KNN::Int, econ::Bool) where {D, T}

    ε★, _ = pecora(s, Tuple(τ_vals), Tuple(ts_vals); delays = τs, w = w,
                samplesize = samplesize, K = K, metric = metric, α = α,
                p = p, undersampling = false)
    ε★s[:,counter] = ε★

    # zero-padding of ⟨ε★⟩ in order to also cover τ=0 (important for the multivariate case)
    ε★ = vec([0; ε★])
    # get the L-statistic-decrease for each peak in ⟨ε★⟩ and take the maximum one
    L_trials, max_idx = local_L_statistics(ε★, Y_act, s, τs, KNN, w, samplesize, metric, econ)
    L_min, min_idx = findmin(L_trials)

    push!(τ_vals, τs[max_idx[min_idx]-1])
    push!(ts_vals, 1)
    push!(Ls, L_min)

    # create phase space vector for this embedding cycle
    Y_act = DelayEmbeddings.hcat_lagged_values(Y_act, s, τ_vals[counter+1])

    return Y_act
end

"""
Perform one embedding cycle on `Y_act` with a multivariate set Ys
"""
function pecuzal_multivariate_embedding_cycle!(Y_act, flag::Bool,
        Ys::Dataset{DT, T}, τs, w::Int, counter::Int, ε★s::AbstractMatrix,
        τ_vals::Vector{Int}, metric, Ls::Vector{T}, ts_vals::Vector{Int},
        samplesize::Real, K::Int, α::Real, p::Real, KNN::Int, econ::Bool
        ) where {D, DT, T}

    M = size(Ys,2)
    # in the 1st cycle we have to check all (size(Y,2)^2 combinations and pick
    # the tau according to maximum L-statistic decrease)
    if counter == 1
        Y_act = first_embedding_cycle_pecuzal!(Ys, M, τs, w, samplesize, K,
                                metric, α, p, KNN, τ_vals, ts_vals, Ls, ε★s, econ)
    # in all other cycles we just have to check (size(Y,2)) combinations and pick
    # the tau according to maximum L-statistic decrease
    else
        Y_act = embedding_cycle_pecuzal!(Y_act, Ys, counter, M, τs, w, samplesize,
                            K, metric, α, p, KNN, τ_vals, ts_vals, Ls, ε★s, econ)
    end
    return Y_act
end

"""
Perform the first embedding cycle of the multivariate embedding. Return the
actual reconstruction vector `Y_act`.
"""
function first_embedding_cycle_pecuzal!(Ys::Dataset{D, T}, M::Int, τs, w::Int,
            samplesize::Real, K::Int, metric, α::Real, p::Real, KNN::Int,
            τ_vals::Vector{Int}, ts_vals::Vector{Int}, Ls::Vector{T},
            ε★s::AbstractMatrix, econ::Bool) where {D, T}
    counter = 1
    L_min = zeros(M)
    L_act = zeros(M)
    L_min_idx = zeros(Int, M)
    ε★ = zeros(length(τs), M*M)
    idx = zeros(Int, M)
    ξ_min = zeros(M)
    for ts = 1:M
        ε★[:,1+(M*(ts-1)):M*ts], _ = pecora(Ys, (0,), (ts,); delays = τs,
                    w = w, samplesize = samplesize, K = K, metric = metric,
                    α = α, p = p, undersampling = false)
        L_min[ts], L_min_idx[ts], idx[ts]  = choose_right_embedding_params(
                                        ε★[:,1+(M*(ts-1)):M*ts], Ys[:,ts],
                                        Ys, τs, KNN, w, samplesize,
                                        metric, econ)
    end
    L_mini, min_idx = findmin(L_min)
    # update τ_vals, ts_vals, Ls, ε★s
    push!(τ_vals, τs[L_min_idx[min_idx]])
    push!(ts_vals, min_idx)             # time series to start with
    push!(ts_vals, idx[min_idx])        # result of 1st embedding cycle
    push!(Ls, L_mini)                   # L-value of 1st embedding cycle
    ε★s[counter] = ε★[:,1+(M*(ts_vals[1]-1)):M*ts_vals[1]]

    # create phase space vector for this embedding cycle
    Y_act = DelayEmbeddings.hcat_lagged_values(Ys[:,ts_vals[counter]],
                                 Ys[:,ts_vals[counter+1]],τ_vals[counter+1])

    return Y_act
end

"""
Perform an embedding cycle of the multivariate embedding, but the first one.
Return the actual reconstruction vector `Y_act`.
"""
function embedding_cycle_pecuzal!(Y_act::Dataset{D, T}, Ys::Dataset{DT, T},
            counter::Int, M::Int, τs, w::Int, samplesize::Real, K::Int, metric,
            α::Real, p::Real, KNN::Int, τ_vals::Vector{Int}, ts_vals::Vector{Int},
            Ls::Vector{T}, ε★s::AbstractMatrix, econ::Bool) where {D, DT, T}

    ε★, _ = pecora(Ys, Tuple(τ_vals), Tuple(ts_vals); delays = τs, w = w,
            samplesize = samplesize, K = K, metric = metric, α = α,
            p = p, undersampling = false)
    # update τ_vals, ts_vals, Ls, ε★s
    choose_right_embedding_params!(ε★, Y_act, Ys, τ_vals, ts_vals, Ls, ε★s,
                                counter, τs, KNN, w, samplesize, metric, econ)
    # create phase space vector for this embedding cycle
    Y_act = DelayEmbeddings.hcat_lagged_values(Y_act, Ys[:, ts_vals[counter+1]],
                                                        τ_vals[counter+1])
    return Y_act
end

"""
    Choose the maximum L-decrease and corresponding τ for each ε★-statistic,
based on picking the peak in ε★, which corresponds to the minimal `L`-statistic.
"""
function choose_right_embedding_params!(ε★::AbstractMatrix, Y_act,
            Ys::Dataset{D, T}, τ_vals::Vector{Int}, ts_vals::Vector{Int},
            Ls::Vector{T}, ε★s::AbstractMatrix, counter::Int, τs, KNN::Int,
            w::Int, samplesize::Real, metric, econ::Bool) where {D, T}

    L_min_ = zeros(size(Ys,2))
    τ_idx = zeros(Int,size(Ys,2))
    for ts = 1:size(Ys,2)
        # zero-padding of ⟨ε★⟩ in order to also cover τ=0 (important for the multivariate case)
        # get the L-statistic for each peak in ⟨ε★⟩ and take the one according to L_min
        L_trials_, max_idx_ = local_L_statistics(vec([0; ε★[:,ts]]), Y_act, Ys[:,ts],
                                        τs, KNN, w, samplesize, metric, econ)
        L_min_[ts], min_idx_ = findmin(L_trials_)
        τ_idx[ts] = max_idx_[min_idx_]-1
    end
    idx = sortperm(L_min_)
    L_mini, min_idx = findmin(L_min_)
    push!(τ_vals, τs[τ_idx[min_idx]])
    push!(ts_vals, min_idx)
    push!(Ls, L_mini)

    ε★s[counter] = ε★
end

"""
    Choose the right embedding parameters of the ε★-statistic in the first
embedding cycle. Return the `L`-decrease-value, the corresponding index value of
the chosen peak `τ_idx` and the number of the chosen time series to start with
`idx`.
"""
function choose_right_embedding_params(ε★::AbstractMatrix, Y_act,
            Ys::Dataset{D, T}, τs, KNN::Int, w::Int, samplesize::Real, metric,
            econ::Bool) where {D, T}
    L_min_ = zeros(size(Ys,2))
    τ_idx = zeros(Int,size(Ys,2))
    for ts = 1:size(Ys,2)
        # zero-padding of ⟨ε★⟩ in order to also cover τ=0 (important for the multivariate case)
        # get the L-statistic for each peak in ⟨ε★⟩ and take the one according to L_min
        L_trials_, max_idx_ = local_L_statistics(vec([0; ε★[:,ts]]), Dataset(Y_act), Ys[:,ts],
                                        τs, KNN, w, samplesize, metric, econ)
        L_min_[ts], min_idx_ = findmin(L_trials_)
        τ_idx[ts] = max_idx_[min_idx_]-1
    end
    idx = sortperm(L_min_)
    return L_min_[idx[1]], τ_idx[idx[1]], idx[1]
end

"""
    Return the maximum decrease of the L-statistic `L_decrease` and corresponding
delay-indices `max_idx` for all local maxima in ε★
"""
function local_L_statistics(ε★::Vector{T}, Y_act::Dataset{D, T}, s::Vector{T},
        τs, KNN::Int, w::Int, samplesize::Real, metric, econ::Bool
        ) where {D, T}
    _, max_idx = get_maxima(ε★) # determine local maxima in ⟨ε★⟩
    L_decrease = zeros(Float64, length(max_idx))
    for (i,τ_idx) in enumerate(max_idx)
        # create candidate phase space vector for this peak/τ-value
        Y_trial = DelayEmbeddings.hcat_lagged_values(Y_act, s, τs[τ_idx-1])
        # compute L-statistic for Y_act and Y_trial and get the maximum decrease
        L_decrease[i] = uzal_cost_pecuzal(Y_act, Y_trial, τs[end]; K = KNN,
                                w = w, metric = metric, econ = econ)
    end
    return L_decrease, max_idx
end

function pecuzal_break_criterion(Ls::Vector{T}, counter::Int,
        max_num_of_cycles::Int, threshold::Real) where {T}
    flag = true
    if counter == 1 && Ls[end] > threshold
        println("Algorithm stopped due to increasing L-values in the first embedding cycle. "*
                "Valid embedding NOT achieved ⨉.")
        flag = false
    end
    if counter > 1 && Ls[end] > threshold
        println("Algorithm stopped due to increasing L-values. "*
                "VALID embedding achieved ✓.")
        flag = false
    end
    if max_num_of_cycles == counter
        println("Algorithm stopped due to hitting max cycle number. "*
                "Valid embedding NOT achieved ⨉.")
        flag = false
    end
    return flag
end

"""
Return the maxima of the given time series s and its indices
"""
function get_maxima(s::Vector{T}) where {T}
    maximas = T[]
    maximas_idx = Int[]
    N = length(s)
    flag = false
    first_point = 0
    for i = 2:N-1
        if s[i-1] < s[i] && s[i+1] < s[i]
            flag = false
            push!(maximas, s[i])
            push!(maximas_idx, i)
        end
        # handling constant values
        if flag
            if s[i+1] < s[first_point]
                flag = false
                push!(maximas, s[first_point])
                push!(maximas_idx, first_point)
            elseif s[i+1] > s[first_point]
                flag = false
            end
        end
        if s[i-1] < s[i] && s[i+1] == s[i]
            flag = true
            first_point = i
        end
    end
    # make sure there is no empty vector returned
    if isempty(maximas)
        maximas, maximas_idx = findmax(s)
    end
    return maximas, maximas_idx
end

"""
    uzal_cost_pecuzal(Y1::Dataset, Y2::Dataset, Tw; kwargs...) → L_decrease
This function is based on the functionality of [`uzal_cost`](@ref), here
specifically tailored for the needs in the PECUZAL algorithm.
Compute the L-statistics `L1` and `L2` for the input datasets `Y1` and `Y2` for
increasing time horizons `T = 1:Tw`. For each `T`, compute `L1` and `L2` and
decrease `L_decrease = L2 - L1`. If `L_decrease` is a negative value, then `Y2`
can be regarded as a "better" reconstruction that `Y1`. Break, when `L_decrease`
reaches the 1st local minima, since this will typically also be the global
minimum. Return the according minimum `L_decrease`-value.

## Keyword arguments

* `K = 3`: the amount of nearest neighbors considered, in order to compute σ_k^2
  (read algorithm description).
  If given a vector, minimum result over all `k ∈ K` is returned.
* `metric = Euclidean()`: metric used for finding nearest neigbhors in the input
  state space trajectory `Y.
* `w = 1`: Theiler window (neighbors in time with index `w` close to the point,
  that are excluded from being true neighbors). `w=0` means to exclude only the
  point itself, and no temporal neighbors.
* `econ::Bool = false`: Economy-mode for L-statistic computation. Instead of
  computing L-statistics for time horizons `2:Tw`, here we only compute them for
  `2:2:Tw`.
"""
function uzal_cost_pecuzal(Y::Dataset{D, ET}, Y_trial::Dataset{DT, ET}, Tw::Int;
        K::Int = 3, w::Int = 1, econ::Bool = false, metric = Euclidean()
    ) where {D, DT, ET}

    @assert DT == D+1
    @assert Tw ≥ 0

    if econ
        tws = 2:2:Tw # start at 2 will eliminate bad results for noise
    else
        tws = 2:Tw # start at 2 will eliminate bad results for noise
    end

    NNN = length(Y_trial)-1
    # preallocation for 1st dataset
    ϵ² = zeros(NNN)             # neighborhood size
    E² = zeros(NNN, length(tws))         # conditional variance
    ϵ_ball = zeros(ET, K+1, D)  # epsilon neighbourhood
    u_k = zeros(ET, D)          # center of mass
    # preallocation for 2nd dataset
    ϵ²_trial = zeros(NNN)             # neighborhood size
    E²_trial = zeros(NNN, length(tws))         # conditional variance
    ϵ_ball_trial = zeros(ET, K+1, DT) # epsilon neighbourhood
    u_k_trial = zeros(ET, DT)         # center of mass

    dist_former = 9999999 # intial L-decrease

    # loop over each time horizon
    cnt = 1
    for T in tws
        NN = length(Y_trial)-T
        if NN < 1
            error("Time series too short for given possible delays and Theiler window to find enough nearest neighbours")
        end
        ns = 1:NN

        vs = Y[ns] # the fiducial points in the data set
        vs_trial = Y_trial[ns] # the fiducial points in the data set

        vtree = KDTree(Y[1:NN], metric)
        allNNidxs, allNNdist = all_neighbors(vtree, vs, ns, K, w)
        vtree_trial = KDTree(Y_trial[1:NN], metric)
        allNNidxs_trial, allNNdist_trial = all_neighbors(vtree_trial, vs_trial, ns, K, w)

        # compute conditional variances and neighborhood-sizes
        compute_conditional_variances!(ns, vs, vs_trial, allNNidxs,
            allNNidxs_trial, Y, Y_trial, ϵ_ball, ϵ_ball_trial, u_k, u_k_trial,
            T, K, metric, ϵ², ϵ²_trial, E², E²_trial, cnt)

        # compute distance of L-values and check whether that distance can be
        # increased
        dist = compute_L_decrease(E², E²_trial, ϵ², ϵ²_trial, cnt, NN)
        if dist > dist_former && dist_former<0
            break
        else
            dist_former = dist
        end
        cnt += 1
    end
    return dist_former
end

function compute_conditional_variances!(ns, vs, vs_trial, allNNidxs::Vector{Array{Int64,1}},
        allNNidxs_trial::Vector{Array{Int64,1}}, Y::Dataset{D, P},
        Y_trial::Dataset{DT, P}, ϵ_ball::Array{P, 2}, ϵ_ball_trial::Array{P, 2},
        u_k::Vector{P}, u_k_trial::Vector{P}, T::Int, K::Int, metric, ϵ²::Vector,
        ϵ²_trial::Vector, E²::Array{P, 2}, E²_trial::Array{P, 2}, cnt::Int) where {P, D, DT}

    # loop over each point on the trajectories
    for (i,v) in enumerate(vs)
        NNidxs = allNNidxs[i] # indices of k nearest neighbors to v of 1st dataset
        NNidxs_trial = allNNidxs_trial[i] # indices of k nearest neighbors to v of 2nd dataset
        # pairwise distance of fiducial points and all neighbours
        pdsqrd = fiducial_pairwise_dist_sqrd(view(Y.data, NNidxs), v, metric)
        ϵ²[i] = (2/(K*(K+1))) * pdsqrd  # Eq. 16
        pdsqrd_trial = fiducial_pairwise_dist_sqrd(view(Y_trial.data, NNidxs_trial), vs_trial[i], metric)
        ϵ²_trial[i] = (2/(K*(K+1))) * pdsqrd_trial  # Eq. 16

        E²[i, cnt] = comp_Ek2!(ϵ_ball, u_k, Y, ns[i], NNidxs, T, K, metric) # Eqs. 13 & 14
        E²_trial[i, cnt] = comp_Ek2!(ϵ_ball_trial, u_k_trial, Y_trial, ns[i], NNidxs_trial, T, K, metric) # Eqs. 13 & 14
    end
end

function compute_L_decrease(E²::Array{P, 2}, E²_trial::Array{P, 2}, ϵ²::Vector{P},
        ϵ²_trial::Vector{P}, T::Int, NN::Int) where {P}
    # 1st dataset
    # Average E²[T] over all prediction horizons
    E²_avrg = mean(E²[1:NN,1:T], dims=2)                   # Eq. 15
    σ² = E²_avrg ./ ϵ²[1:NN] # noise amplification σ², Eq. 17
    σ²_avrg = mean(σ²) # averaged value of the noise amplification, Eq. 18
    α² = 1 / sum(ϵ²[1:NN].^(-1)) # for normalization, Eq. 21
    L = log10(sqrt(σ²_avrg)*sqrt(α²))
    # 2nd dataset
    # Average E²[T] over all prediction horizons
    E²_avrg_trial = mean(E²_trial[1:NN,1:T], dims=2)                   # Eq. 15
    σ²_trial = E²_avrg_trial ./ ϵ²_trial[1:NN] # noise amplification σ², Eq. 17
    σ²_avrg_trial = mean(σ²_trial) # averaged value of the noise amplification, Eq. 18
    α²_trial = 1 / sum(ϵ²_trial[1:NN].^(-1)) # for normalization, Eq. 21
    L_trial = log10(sqrt(σ²_avrg_trial)*sqrt(α²_trial))

    return L_trial - L
end

function fiducial_pairwise_dist_sqrd(fiducials, v, metric)
    pd = zero(eltype(fiducials[1]))
    pd += evaluate(metric, v, v)^2
    for (i, v1) in enumerate(fiducials)
        pd += evaluate(metric, v1, v)^2
        for j in i+1:length(fiducials)
            @inbounds pd += evaluate(metric, v1, fiducials[j])^2
        end
    end
    return sum(pd)
end

"""
    comp_Ek2!(ϵ_ball,u_k,Y,v,NNidxs,T,K,metric) → E²(T)
Returns the approximated conditional variance for a specific point in state space
`ns` (index value) with its `K`-nearest neighbors, which indices are stored in
`NNidxs`, for a time horizon `T`. This corresponds to Eqs. 13 & 14 in [^Uzal2011].
The norm specified in input parameter `metric` is used for distance computations.
"""
function comp_Ek2!(ϵ_ball::Array{P, 2}, u_k::Vector{P}, Y::Dataset{D, P},
        ns::Int, NNidxs::Vector, T::Int, K::Int, metric
        ) where {D, P}
    # determine neighborhood `T` time steps ahead
    ϵ_ball[1, :] .= Y[ns+T]
    @inbounds for (i, j) in enumerate(NNidxs)
        ϵ_ball[i+1, :] .= Y[j + T]
    end

    # compute center of mass
    @inbounds for i in 1:size(Y)[2]; u_k[i] = sum(view(ϵ_ball, :, i))/(K+1); end # Eq. 14

    E²_sum = 0
    @inbounds for j = 1:K+1
        E²_sum += (evaluate(metric, view(ϵ_ball, j, :), u_k))^2
    end
    E² = E²_sum / (K+1)         # Eq. 13
end

"""
    all_neighbors(vtree, vs, ns, K, w)
Return the `maximum(K)`-th nearest neighbors for all input points `vs`,
with indices `ns` in original data, while respecting the theiler window `w`.

This function is nothing more than a convinience call to `Neighborhood.bulksearch`.
"""
function all_neighbors(vtree, vs, ns, K, w)
    w ≥ length(vtree.data)-1 && error("Theiler window larger than the entire data span!")
    k = maximum(K)
    tw = Theiler(w, ns)
    idxs, dists = bulksearch(vtree, vs, NeighborNumber(k), tw)
end


text:
# io/z-test.jl

using Test

@test ir_dump(:test)
@test fld_read(:test)


text:
using DifferentialEquations
using Plots

function sis_ode(du,u,p,t)
        SH,IH,SL,IL = u
        betaHH,betaHL,betaLH,betaLL,gamma=p
        du[1]=-(betaHH*IH+betaHL*IL)*SH+gamma*IH
        du[2]=+(betaHH*IH+betaHL*IL)*SH-gamma*IH
        du[3]=-(betaLH*IH+betaLL*IL)*SL+gamma*IL
        du[4]=+(betaLH*IH+betaLL*IL)*SL-gamma*IL
    end

parms =[10,0.1,0.1,1,1]
init=[0.19999,0.00001,0.799,0.001]
tspan=tspan = (0.0,15.0)

sis_prob = ODEProblem(sis_ode,init,tspan,parms)
sis_sol = solve(sis_prob,saveat=0.1)

plot(sis_sol,xlabel="Time (Years)",ylabel="Proportion of Population")

max_stars_repo_name: epimodels/epicookbook

max_issues_repo_name: epimodels/epicookbook

max_forks_repo_name: epimodels/epicookbook