text:
@doc raw"""
    differential_geodesic_startpoint(M, p, q, t, X)
computes $D_p g(t;p,q)[\eta]$.

# See also
 [`differential_geodesic_endpoint`](@ref), [`jacobi_field`](@ref)
"""
differential_geodesic_startpoint(M::mT,x,y,t,η) where {mT <: Manifold} = jacobi_field(M,x,y,t,η,βdifferential_geodesic_startpoint)
@doc raw"""
    differential_geodesic_endpoint(M,x,y,t,η)
computes $D_qg(t;p,q)[\eta]$.

# See also
 [`differential_geodesic_startpoint`](@ref), [`jacobi_field`](@ref)
"""
differential_geodesic_endpoint(M::mT, x, y, t, η) where {mT <: Manifold} = jacobi_field(M,y,x,1-t,η,βdifferential_geodesic_startpoint)
@doc raw"""
    differential_exp_basepoint(M, p, X, Y)

Compute $D_p\exp_p X[Y]$.

# See also
[`differential_exp_argument`](@ref), [`jacobi_field`](@ref)
"""
function differential_exp_basepoint(M::Manifold,p,X,Y)
    return jacobi_field(M, p, exp(M,p,X), 1.0, Y, βdifferential_exp_basepoint)
end
@doc raw"""
    differential_exp_argument(M, p, X, Y)
computes $D_X\exp_pX[Y]$.
Note that $X ∈  T_X(T_p\mathcal M) = T_p\mathcal M$ is still a tangent vector.

# See also
 [`differential_exp_basepoint`](@ref), [`jacobi_field`](@ref)
"""
function differential_exp_argument(M::Manifold, p, X, Y)
    return jacobi_field(M,p,exp(M,p,X),1.0,Y,βdifferential_exp_argument)
end

@doc raw"""
    differential_log_basepoint(M, p, q, X)
computes $D_p\log_pq[X]$.

# See also
 [`differential_log_argument`](@ref), [`jacobi_field`](@ref)
"""
function differential_log_basepoint(M::Manifold, p, q, X)
    return jacobi_field(M,p,q,0.0,X,βdifferential_log_basepoint)
end

@doc raw"""
    differential_log_argument(M,p,q,X)
computes $D_q\log_p,q[X]$.

# See also
 [`differential_log_argument`](@ref), [`jacobi_field`](@ref)
"""
function differential_log_argument(M::Manifold, p, q, X)
    return jacobi_field(M,q,p,1.0,X,βdifferential_log_argument)
end

@doc raw"""
    Y = differential_forward_logs(M, p, X)

compute the differenital of [`forward_logs`](@ref) $F$ on the `PowerManifold` manifold
`M` at `p` and direction `X` ,
in the power manifold array, the differential of the function

```math
F_i(x) = \sum_{j ∈ \mathcal I_i} \log_{p_i} p_j$, \quad i  ∈  \mathcal G,
```

where $\mathcal G$ is the set of indices of the `PowerManifold` manifold `M`
and $\mathcal I_i$ denotes the forward neighbors of $i$.

# Input
* `M`     – a `PowerManifold` manifold
* `p`     – a point.
* `X`     – a tangent vector.

# Ouput
* `Y` – resulting tangent vector in $T_x\mathcal N$ representing the differentials of the logs, where
  $\mathcal N$ is thw power manifold with the number of dimensions added to `size(x)`.
"""
function differential_forward_logs(M::PowerManifold, p, X)
    power_size = power_dimensions(M)
    R = CartesianIndices(Tuple(power_size))
    d = length(power_size)
    maxInd = [last(R).I...]
    d2 = (d>1) ? ones(Int,d+1) + (d-1)*(1:(d+1) .== d+1 ) : 1
    if d > 1
        N = PowerManifold(M.manifold, NestedPowerRepresentation(), power_size...,d)
    else
        N = PowerManifold(M.manifold, NestedPowerRepresentation(), power_size...)
    end
    Y = zero_tangent_vector(N, repeat(p,inner=d2) )
    for i in R # iterate over all pixel
        for k in 1:d # for all direction combinations
            I = [i.I...] # array of index
            J = I .+ 1 .* (1:d .== k) #i + e_k is j
            if all( J .<= maxInd )
                # this is neighbor in range,
                j = CartesianIndex{d}(J...) # neigbbor index as Cartesian Index
                # collects two, namely in kth direction since xi appears as base and arg
                # Y[i,k] = differential_log_basepoint(M.manifold,p[i],p[j],X[i])
                #            + differential_log_argument(M.manifold,p[i],p[j],X[j])
                set_component!(
                    N,
                    Y,
                    differential_log_basepoint(
                        M.manifold,
                        get_component(M,p,i),
                        get_component(M,p,j),
                        get_component(M,X,i),
                    ) + differential_log_argument(
                        M.manifold,
                        get_component(M,p,i),
                        get_component(M,p,j),
                        get_component(M,X,j),
                    ),
                    i,
                    k,
                )
            end
        end # directions
    end # i in R
    return Y
end


text:
module Gherkin

import Base: ==, hash

export Scenario, ScenarioOutline, Feature, FeatureHeader, Given, When, Then
export parsefeature, hastag

"A good or bad result when parsing Gherkin."
abstract type ParseResult{T} end

"A successful parse that results in the expected value."
struct OKParseResult{T} <: ParseResult{T}
    value::T
end

"An unsuccessful parse that results in an error."
struct BadParseResult{T} <: ParseResult{T}
    reason::Symbol
    expected::Symbol
    actual::Symbol
end

"Was the parsing successful?"
issuccessful(::OKParseResult{T}) where {T} = true
issuccessful(::BadParseResult{T}) where {T} = false

"A step in a Gherkin Scenario."
abstract type ScenarioStep end

"Equality for scenario steps is their text and their block text."
function ==(a::T, b::T) where {T <: ScenarioStep}
    a.text == b.text && a.block_text == b.block_text
end
"Hash scenario steps by their text and block text."
hash(a::T, h::UInt) where {T <: ScenarioStep} = hash((a.text, a.block_text), h)

"A Given scenario step."
struct Given <: ScenarioStep
    text::String
    block_text::String

    Given(text::AbstractString; block_text = "") = new(text, block_text)
end
"A When scenario step."
struct When <: ScenarioStep
    text::String
    block_text::String

    When(text::AbstractString; block_text="") = new(text, block_text)
end
"A Then scenario step."
struct Then <: ScenarioStep
    text::String
    block_text::String

    Then(text::AbstractString; block_text="") = new(text, block_text)
end

"An AbstractScenario is a Scenario or a Scenario Outline in Gherkin."
abstract type AbstractScenario end

"""
A Gherkin Scenario.

# Example
```
@tag1 @tag2
Scenario: Some description
    Given some precondition
     When some action is taken
     Then some postcondition holds
```

becomes a `Scenario` struct with description "Some description", tags `["@tag1", "@tag2"]`, and
steps
```[
    Given("some precondition"),
    When("some action is taken"),
    Then("some postcondition holds")
]```.
"""
struct Scenario <: AbstractScenario
    description::String
    tags::Vector{String}
    steps::Vector{ScenarioStep}
end

"""
A `ScenarioOutline` is a scenario with multiple examples.

# Example
```
@tag1 @tag2
Scenario: Some description
    Given some precondition
     When some action is taken
     Then some postcondition holds
```
"""
struct ScenarioOutline <: AbstractScenario
    description::String
    tags::Vector{String}
    steps::Vector{ScenarioStep}
    placeholders::Vector{String}
    examples::AbstractArray
end

"""
A FeatureHeader has a (short) description, a longer description, and a list of applicable tags.

# Example
```
@tag1 @tag2
Feature: This is a description
    This is
    a longer description.
```
"""
struct FeatureHeader
    description::String
    long_description::Vector{String}
    tags::Vector{String}
end

"""
A Feature has a feature header and a list of Scenarios and Scenario Outlines.
"""
struct Feature
    header::FeatureHeader
    scenarios::Vector{AbstractScenario}
end

"""
ByLineParser takes a long text and lets the Gherkin parser work line by line.
"""
mutable struct ByLineParser
    current::String
    rest::Vector{String}
    isempty::Bool

    function ByLineParser(text::String)
        lines = split(text, "\n")
        current = lines[1]
        rest = lines[2:end]
        new(current, rest, false)
    end
end

"""
    consume!(p::ByLineParser)

The current line has been consumed by the Gherkin parser.
"""
function consume!(p::ByLineParser)
    if isempty(p.rest)
        p.current = ""
        p.isempty = true
    else
        p.current = p.rest[1]
        p.rest = p.rest[2:end]
    end
end

"""
    @untilemptyline(byline::Symbol,)

Execute the function for each line, until an empty line is encountered, or no further lines are
available.
"""
macro untilemptyline(ex::Expr)
    esc(quote
        while !isempty(byline)
            if iscurrentlineempty(byline)
                consume!(byline)
                break
            end
            $ex
            consume!(byline)
        end
    end)
end

"""
    ignoringemptylines!(f::Function, byline::ByLineParser)

Execute function `f` for all non-empty lines, until the end of the file.
"""
macro ignoringemptylines(ex::Expr)
    esc(quote
        while !isempty(byline)
            if iscurrentlineempty(byline)
                consume!(byline)
                continue
            end

            $ex
        end
    end)
end

"""
Override the normal `isempty` method for `ByLineParser`.
"""
Base.isempty(p::ByLineParser) = p.isempty

"""
    iscurrentlineempty(p::ByLineParser)

Check if the current line in the parser is empty.
"""
iscurrentlineempty(p::ByLineParser) = strip(p.current) == ""

"""
    parsetags!(byline::ByLineParser)

Parse all tags on form @tagname until a non-tag is encountered. Return a possibly empty list of
tags.
"""
function parsetags!(byline::ByLineParser)
    tags = []
    while !isempty(byline)
        # Use a regular expression to find all occurrences of @tagname on a line.
        tag_match = collect((m.match for m = eachmatch(r"(@[^\s]+)", byline.current)))

        # Break if something that isn't a list of tags is encountered.
        if isempty(tag_match)
            break
        end

        consume!(byline)
        append!(tags, tag_match)
    end

    tags
end

"""
    parsefeatureheader!(byline::ByLineParser) :: ParseResult{FeatureHeader}

Parse a feature header and return the feature header on success, or a bad parse result on failure.

# Example of a feature header
```
@tag1
@tag2
Feature: Some description
    A longer description
    on multiple lines.
```
"""
function parsefeatureheader!(byline::ByLineParser) :: ParseResult{FeatureHeader}
    feature_tags = parsetags!(byline)

    description_match = match(r"Feature: (?<description>.+)", byline.current)
    if description_match == nothing
        return BadParseResult{FeatureHeader}(:unexpected_construct, :feature, :scenario)
    end
    consume!(byline)

    # Consume all lines after the `Feature:` row as a long description, until an empty line is
    # encountered.
    long_description_lines = []
    # untilemptyline!(byline) do
    #     push!(long_description_lines, strip(byline.current))
    # end
    @untilemptyline begin
        push!(long_description_lines, strip(byline.current))
    end

    feature_header = FeatureHeader(description_match[:description],
                                   long_description_lines,
                                   feature_tags)
    return OKParseResult{FeatureHeader}(feature_header)
end

"""
    parseblocktext!(byline::ByLineParser)

Parse a block of text that may occur as part of a scenario step.
A block of text is multiple lines of text surrounded by three double quotes.
For example, this comment is a valid block text.

Precondition: This function assumes that the current line is the starting line of three double
quotes.
"""
function parseblocktext!(byline::ByLineParser)
    # Consume the current line with three double quotes """
    consume!(byline)
    block_text_lines = []

    # Consume and save all lines, until we reach one that is the ending line of three double quotes.
    @untilemptyline begin
        line = byline.current
        if occursin(r"\"\"\"", line)
            consume!(byline)
            break
        end
        push!(block_text_lines, strip(line))
    end
    return OKParseResult{String}(join(block_text_lines, "\n"))
end

"""
    parsescenariosteps!(byline::ByLineParser)

Parse all scenario steps following a Scenario or Scenario Outline.
"""
function parsescenariosteps!(byline::ByLineParser)
    steps = []
    allowed_step_types = Set([Given, When, Then])

    @untilemptyline begin
        # Match Given, When, or Then on the line, or match a block text.
        # Note: This is a place where English Gherkin is hard coded.
        step_match = match(r"(?<step_type>Given|When|Then|And) (?<step_definition>.+)", byline.current)
        block_text_start_match = match(r"\"\"\"", byline.current)
        # A line must either be a new scenario step, or a block text following the previous scenario
        # step.
        if step_match == nothing && block_text_start_match == nothing
            return BadParseResult{Vector{ScenarioStep}}(:invalid_step, :step_definition, :invalid_step_definition)
        end

        # The current line starts a block text. Parse the rest of the block text and continue with
        # the following line.
        # Note: The method parseblocktext!(byline) consumes the current line, so it doesn't need to
        # be done here.
        if block_text_start_match != nothing
            block_text_result = parseblocktext!(byline)
            prev_step_type = typeof(steps[end])
            steps[end] = prev_step_type(steps[end].text; block_text=block_text_result.value)
            continue
        end

        # The current line is a scenario step.
        # Check that the scenario step is allowed, and create a Given/When/Then object.
        # Note that scenario steps must occur in the order Given, When, Then. A Given may not follow
        # once a When has been seen, and a When must not follow when a Then has been seen.
        # This is what `allowed_step_types` keeps track of.
        step_type = step_match[:step_type]
        step_definition = step_match[:step_definition]
        if step_type == "Given"
            if !(Given in allowed_step_types)
                return BadParseResult{Vector{ScenarioStep}}(:bad_step_order, :NotGiven, :Given)
            end
            step = Given(step_definition)
        elseif step_type == "When"
            if !(When in allowed_step_types)
                return BadParseResult{Vector{ScenarioStep}}(:bad_step_order, :NotWhen, :When)
            end
            step = When(step_definition)
            delete!(allowed_step_types, Given)
        elseif step_type == "Then"
            step = Then(step_definition)
            delete!(allowed_step_types, Given)
            delete!(allowed_step_types, When)
        elseif step_type == "And"
            # A scenario step may be And, in which case it's the same type as the previous step.
            # This means that an And may not be the first scenario step.
            if isempty(steps)
                return BadParseResult{Vector{ScenarioStep}}(:leading_and, :specific_step, :and_step)
            end
            last_specific_type = typeof(steps[end])
            step = last_specific_type(step_definition)
        end
        push!(steps, step)
    end
    return OKParseResult{Vector{ScenarioStep}}(steps)
end

"""
    parsescenario!(byline::ByLineParser)

Parse a Scenario or a Scenario Outline.

# Example of a scenario
```
@tag1 @tag2
Scenario: Some description
    Given some precondition
     When some action is taken
     Then some postcondition holds
```
"""
function parsescenario!(byline::ByLineParser)
    tags = parsetags!(byline)

    # The scenario is either a Scenario or a Scenario Outline. Check for Scenario Outline first.
    scenario_outline_match = match(r"Scenario Outline: (?<description>.+)", byline.current)
    if scenario_outline_match != nothing
        description = scenario_outline_match[:description]
        consume!(byline)

        # Parse scenario outline steps.
        steps_result = parsescenariosteps!(byline)
        if !issuccessful(steps_result)
            return steps_result
        end
        steps = steps_result.value

        # Consume the Example: line.
        consume!(byline)

        # Get the name of each placeholder variable.
        placeholders = collect((m.match for m = eachmatch(r"(\w+)", byline.current)))
        consume!(byline)

        # Parse the examples, until we hit an empty line.
        # TODO: This needs to be updated to allow for multiple Examples sections.
        examples = Array{String,2}(undef, length(placeholders), 0)
        @untilemptyline begin
            # Each variable is in a column, separated by |
            example = split(strip(byline.current), "|")
            filter!(x -> !isempty(x), example)
            # Remove surrounding whitespace around each value.
            example = map(strip, example)
            examples = [examples example]
        end
        return OKParseResult{ScenarioOutline}(
            ScenarioOutline(description, tags, steps, placeholders, examples))
    end

    # Here we parse a normal Scenario instead.
    scenario_match = match(r"Scenario: (?<description>.+)", byline.current)
    description = scenario_match[:description]
    consume!(byline)

    steps_result = parsescenariosteps!(byline)
    if !issuccessful(steps_result)
        return steps_result
    end
    steps = steps_result.value

    OKParseResult{Scenario}(Scenario(description, tags, steps))
end

"""
    parsefeature(::String)

Parse an entire feature file.

# Example of a feature file
```
@tag1 @tag2
Feature: Some feature description
    This is some block text,
    and it may be multiple lines.

    Scenario: Some scenario description
        Given some precondition
         When some action is taken
         Then some postcondition holds
```
"""
function parsefeature(text::String) :: ParseResult{Feature}
    byline = ByLineParser(text)

    # The feature header includes all feature level tags, the description, and the long multiline
    # description.
    feature_header_result = parsefeatureheader!(byline)
    if !issuccessful(feature_header_result)
        return BadParseResult{Feature}(feature_header_result.reason,
                                       feature_header_result.expected,
                                       feature_header_result.actual)
    end

    # Each `parsescenario!`
    scenarios = []
    @ignoringemptylines begin
        scenario_parse_result = parsescenario!(byline)
        if issuccessful(scenario_parse_result)
            push!(scenarios, scenario_parse_result.value)
        end
    end

    OKParseResult{Feature}(
        Feature(feature_header_result.value, scenarios))
end

"""
    hastag(feature::Feature, tag::String)

Check if a feature has a given tag.

# Example
```
feature = parsefeature(featuretext)
hassometag = hastag(feature, "@sometag")
```
"""
hastag(feature::Feature, tag::String) = tag in feature.header.tags
hastag(scenario::Scenario, tag::String) = tag in scenario.tags

end

ext: jl


text:
using JuliaPlot
using Test

@testset "JuliaPlot.jl" begin
    # Write your tests here.
end


text:
local a::int

a = 1


text:
module Bar
using ExportPublic

function foo()
  1
end

function bar()
  2
end

function barBar()
  3
end

function fooFoo()
    _secret_implementation_detail()
end

function _secret_implementation_detail()
    4
end

@exportPublic()
end


text:
"""
    run_gpregression(samplemetafile, betafile, m=50; verbose=true, usepriors=true)

Run analysis for model selection between constant, linear or non-linear (Mat52) kernels

Input:
- `samplefmetafile` - File with sample meta
- `betafile` - `RData` for with probe beta.
- `outdir` - Directory to store results, with following contents
  - `gpstats.tsv` - Tab separated file with marginal likelihoods and optimised kernel parameters for all kernels
  - `gpmeanvar.rdata` - RData file containing GP mean and variance for each probe
- `m = 50` - Number of points to calculate GP mean and variance at between min and max Age.
- usepriors=true - Use default priors on hyperparameters of kernels
"""
function run_gpregression(samplemetafile, betafile, outdir, m=50; verbose=true, usepriors=true)
    nt = Threads.nthreads()

    verbose && println("[GPM]\tRunning GP regression on methylation time series with ", nt, " threads.")
    meta, probes, beta = loadall(samplemetafile, betafile, verbose)
    ap = Float64.(meta.Age)
    st = range(minimum(ap), maximum(ap), length=m)
    gpf = ifelse(usepriors, gpmodels_prior, gpmodels)

    if usepriors
        verbose && println("[GPM]\tUsing priors on kernel parameters.")
    end
    starttime = time()
    dfgp, μ, v = gpreg_all_threads(ap, beta, st, gpf, nt = 16)
    verbose && println("[GPM]\tComplete in ", time() - starttime, " seconds")

    gps = gpstats(probes, dfgp, beta, μ)
    mkpath(outdir)

    f = CSV.write(joinpath(outdir, "gpstats.tsv"), gps, delim='\t')
    verbose && println("[GPM]\tWritten ", f)
    saverdata(probes, st, μ, v, outdir; verbose=verbose)
end

"""
    saverdata(probes, samples, μ, v, outdir; verbose=true)

Save gp mean and var in RData format.
"""
function saverdata(probes, samples, μ, v, outdir; verbose=true)

    file = joinpath(outdir, "gpmeanvar.rdata")

    gp_mean = μ'
    gp_var  = v'

    @rput gp_mean gp_var probes samples

    R"""
        colnames(gp_mean) <- samples
        rownames(gp_mean) <- probes
        colnames(gp_var) <- samples
        rownames(gp_var) <- probes

        save(gp_mean, gp_var, file=$file)
    """

    verbose && println("[GPM]\tWritten ", file)
    file

end

"""
    loadresults_rdata(file)

Loads gp_mean and gp_var results from `RData` file back into Julia.
"""
function loadresults_rdata(file)

    R"""
        load($file)
        st <- colnames(gp_mean)
        probes <- rownames(gp_mean)
    """
    @rget gp_mean gp_var st probes

    (μ=Matrix(gp_mean'), v=Matrix(gp_var'), st=parse.(Float64, st), probes=probes)
end

"""
    loadresults(outdir)

Loads gpstats and GP mean and var
"""
function loadresults(outdir)
    gps = CSV.read(joinpath(outdir, "gpstats.tsv"), DataFrame)
    gpres = loadresults_rdata(joinpath(outdir, "gpmeanvar.rdata"))

    gps, gpres
end

"""
    gpstats(probes, dfgp, beta, gpμ)

Writes tab-separated file with gp stats. Can accept Union{T, Missing} beta, however, will fail if all values of a column are missing.
"""
function gpstats(probes, dfgp, beta, gpμ)

    @assert length(probes) == size(dfgp, 1)

    probe_min = minimum.(skipmissing.(eachcol(beta)))
    probe_max = maximum.(skipmissing.(eachcol(beta)))
    probe_env  = probe_max .- probe_min

    gp_min = minimum.(skipmissing.(eachcol(gpμ)))
    gp_max = maximum.(skipmissing.(eachcol(gpμ)))
    gp_env  = gp_max .- gp_min

    evdf = DataFrame(probe_min=probe_min, probe_max=probe_max, probe_env=probe_env,
                     gp_min=gp_min, gp_max=gp_max, gp_env=gp_env)

    df = [DataFrame(Probe=probes) dfgp evdf]
     rename!(df,
        :const_σn    => :const_sign,
        :const_σf    => :const_sigf,
        :linear_mll  => :linear_mll,
        :linear_σn   => :linear_sign,
        :linear_σf   => :linear_sigf,
        :linear_ℓ    => :linear_ell,
        :mat52_mll   => :mat52_mll,
        :mat52_σn    => :mat52_sign,
        :mat52_σf    => :mat52_sigf,
        :mat52_ℓ     => :mat52_ell)

    df[!, :LR_mat52_const]  = df.mat52_mll .- df.const_mll
    df[!, :LR_mat52_linear] = df.mat52_mll .- df.linear_mll
    df[!, :LR_linear_const] = df.linear_mll .- df.const_mll

    bic_term = log.(vec(sum(.!ismissing.(beta), dims=1)))./2
    df[!, :BIC_mat52_const]  = df.LR_mat52_const - bic_term
    #df[!, :BIC_mat52_linear] = df.LR_mat52_linear ## mat52 and linear models have identical number of params
    df[!, :BIC_linear_const] = df.LR_linear_const - bic_term

    df[!, :PreferedModelLR]   = pref_model.(df.mat52_mll, df.linear_mll, df.const_mll)

     df
end

function pref_model(mat52_mll, linear_mll, const_mll, labels=("Mat52", "Linear", "Const"))
    i = argmax((mat52_mll, linear_mll, const_mll))
    labels[i]
end


text:
export preprocess_describe

using mlpack._Internal.io

import mlpack_jll
const preprocess_describeLibrary = mlpack_jll.libmlpack_julia_preprocess_describe

# Call the C binding of the mlpack preprocess_describe binding.
function preprocess_describe_mlpackMain()
  success = ccall((:preprocess_describe, preprocess_describeLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module preprocess_describe_internal
  import ..preprocess_describeLibrary

end # module

"""
    preprocess_describe(input; [dimension, population, precision, row_major, verbose, width])

This utility takes a dataset and prints out the descriptive statistics of the
data. Descriptive statistics is the discipline of quantitatively describing the
main features of a collection of information, or the quantitative description
itself. The program does not modify the original file, but instead prints out
the statistics to the console. The printed result will look like a table.

Optionally, width and precision of the output can be adjusted by a user using
the `width` and `precision` parameters. A user can also select a specific
dimension to analyze if there are too many dimensions. The `population`
parameter can be specified when the dataset should be considered as a
population.  Otherwise, the dataset will be considered as a sample.

So, a simple example where we want to print out statistical facts about the
dataset `X` using the default settings, we could run

```julia
julia> using CSV
julia> X = CSV.read("X.csv")
julia> preprocess_describe(X; verbose=1)
```

If we want to customize the width to 10 and precision to 5 and consider the
dataset as a population, we could run

```julia
julia> using CSV
julia> X = CSV.read("X.csv")
julia> preprocess_describe(X; precision=5, verbose=1, width=10)
```

# Arguments

 - `input::Array{Float64, 2}`: Matrix containing data,
 - `dimension::Int`: Dimension of the data. Use this to specify a
      dimension  Default value `0`.

 - `population::Bool`: If specified, the program will calculate statistics
      assuming the dataset is the population. By default, the program will
      assume the dataset as a sample.  Default value `false`.

 - `precision::Int`: Precision of the output statistics.  Default value
      `4`.

 - `row_major::Bool`: If specified, the program will calculate statistics
      across rows, not across columns.  (Remember that in mlpack, a column
      represents a point, so this option is generally not necessary.)  Default
      value `false`.

 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

 - `width::Int`: Width of the output table.  Default value `8`.

# Return values

"""
function preprocess_describe(input;
                             dimension::Union{Int, Missing} = missing,
                             population::Union{Bool, Missing} = missing,
                             precision::Union{Int, Missing} = missing,
                             row_major::Union{Bool, Missing} = missing,
                             verbose::Union{Bool, Missing} = missing,
                             width::Union{Int, Missing} = missing,
                             points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, preprocess_describeLibrary), Nothing, ());

  IORestoreSettings("Descriptive Statistics")

  # Process each input argument before calling mlpackMain().
  IOSetParamMat("input", input, points_are_rows)
  if !ismissing(dimension)
    IOSetParam("dimension", convert(Int, dimension))
  end
  if !ismissing(population)
    IOSetParam("population", convert(Bool, population))
  end
  if !ismissing(precision)
    IOSetParam("precision", convert(Int, precision))
  end
  if !ismissing(row_major)
    IOSetParam("row_major", convert(Bool, row_major))
  end
  if !ismissing(width)
    IOSetParam("width", convert(Int, width))
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  # Call the program.
  preprocess_describe_mlpackMain()

  return
end


text:
using Documenter, FixedPointAcceleration

makedocs(
    format = Documenter.HTML(),
    sitename = "FixedPointAcceleration",
    modules = [FixedPointAcceleration],
    pages = ["index.md",
             "1_FixedPoints.md",
             "2_Algorithms.md",
             "3_UsingAdvice.md",
             "4_Applications.md",
             "5_TerminationConditions.md",
             "99_refs.md"]
)

deploydocs(
    repo   = "github.com/s-baumann/FixedPointAcceleration.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing
)


text:
module SampleProblems
	using LinearAlgebra
	using ..RocketlandDefns
	using ..Aerodynamics
	function normalize_problem(dp::DescentProblem)::DescentProblem
		Ul = maximum(dp.rIi)
		Ut = dp.tf_guess
		Um = dp.mwet
		return DescentProblem(
			g=dp.g/(Ul/Ut^2), mdry=dp.mdry/Um, mwet=dp.mwet/Um,
			Tmin = dp.Tmin/(Um * Ul/Ut^2), Tmax = dp.Tmax/(Um * Ul/Ut^2),
			omMax = dp.omMax/Ut, jB = broadcast(*, dp.jB, 1/(Um * Ul^2)),
			rTB = broadcast(*, dp.rTB, 1/Ul), rIi = broadcast(*, dp.rIi, 1/Ul),
			rIf = broadcast(*, dp.rIf, 1/Ul), vIi = broadcast(*, dp.vIi, 1/(Ul/Ut)),
			vIf = broadcast(*, dp.vIi, 1/(Ul/Ut)), qBIf = dp.qBIf, qBIi = dp.qBIi,
			wBi = dp.wBi, wBf = dp.wBf, rFB = broadcast(*, dp.rFB, 1/Ut),
			deltaMax = dp.deltaMax, thetaMax = dp.thetaMax, gammaGs = dp.gammaGs,
			alpha = dp.alpha/(Ut^2/Ul), K=dp.K, imax = dp.imax, wNu = dp.wNu, wID = dp.wID,
			wDS = dp.wDS, wCst = dp.wCst, wTviol = dp.wTviol, delTol = dp.delTol,
			tf_guess = dp.tf_guess/Ut, ri = dp.ri, rh0 = dp.rh0, rh1 = dp.rh1,
			rh2 = dp.rh2, alph = dp.alph, bet = dp.bet, dpMax=dp.dpMax/(Um/(Ul*Ut^2)), rho = dp.rho/(Um/Ul^3), sos=dp.sos/(Ul/Ut),
			aero = Aerodynamics.rescale_aerodata(dp.aero, Ul, Ut, Um))
	end

	aero_info = Aerodynamics.load_aerodata("aero/lift_drag.csv")
	base_prob = DescentProblem(g=9.82, mwet=66018, mdry=65947, Tmin=0.1*4.686588e6, Tmax=4.686588e6, jB=diagm(0=>[72487.03125,2.0734175e6,2.0734175e6]),
										  alpha=0.000345, rTB=[-4.26114,0,0], rFB = [2.0,0,0], rIi = [1000.0,1000.0,100.0], rIf=[0.0,0.0,0.0], vIi = [-100.0,-200.0,0], sos=352.0, wNu = 1e4)
	base_prob_scaled = normalize_problem(base_prob)

	base_prob_aero = DescentProblem(g=9.82, mwet=66018, mdry=65947, Tmin=0.1*4.686588e6, Tmax=4.686588e6, jB=diagm(0=>[72487.03125,2.0734175e6,2.0734175e6]),
										  alpha=0.000345, rTB=[-4.26114,0,0], rFB = [2.0,0,0], rIi = [1000.0,1000.0,100.0], rIf=[0.0,0.0,0.0], vIi = [-100.0,-200.0,0], sos=352.0, aero = aero_info, wNu = 1e4)
	base_prob_aero_scaled = normalize_problem(base_prob_aero)
end

max_stars_repo_name: BenChung/SuccessiveConvexification

max_issues_repo_name: BenChung/SuccessiveConvexification

max_forks_repo_name: BenChung/SuccessiveConvexification


text:
@testset "laplacian" begin
    h = 1.0
    rho = 1000.0
    G = 9.8
    len_z = 16
    len_x = 32
    nz = Int(len_z/h + 1)
    nx = Int(len_x/h + 1)
    tf_h=constant(1.0)
    coef = rand(nz, nx)
    func = rand(nz, nx)

    tf_coef = constant(coef)
    tf_func = constant(func)

    # gradient check -- v
    function scalar_function(m)
        # return sum(tanh(laplacian(m, tf_func, tf_h, constant(rho*G))))
        return sum(tanh(laplacian(tf_coef, m, tf_h, constant(rho*G))))
    end

    # m_  = tf_coef
    m_ = tf_func
    v_ = 0.01*rand(nz, nx)
    y_ = scalar_function(m_)
    dy_ = gradients(y_, m_)
    ms_ = Array{Any}(undef, 5)
    ys_ = Array{Any}(undef, 5)
    s_ = Array{Any}(undef, 5)
    w_ = Array{Any}(undef, 5)
    gs_ =  @. 1 / 10^(1:5)

    for i = 1:5
        g_ = gs_[i]
        ms_[i] = m_ + g_*v_
        ys_[i] = scalar_function(ms_[i])
        s_[i] = ys_[i] - y_
        w_[i] = s_[i] - g_*sum(v_.*dy_)
    end

    sess = Session()
    init(sess)
    sval_ = run(sess, s_)
    wval_ = run(sess, w_)
    close("all")
    loglog(gs_, abs.(sval_), "*-", label="finite difference")
    loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
    loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
    loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

    plt[:gca]()[:invert_xaxis]()
    legend()
    xlabel("\$\\gamma\$")
    ylabel("Error")
    savefig("laplacian.png"); close("all")
end

@testset "poisson" begin
    # gradient check -- v
    h = 1.0
    rho = 1000.0
    G = 9.8
    len_z = 16
    len_x = 32
    nz = Int(len_z/h + 1)
    nx = Int(len_x/h + 1)
    tf_h=constant(1.0)
    coef = zeros(nz, nx)
    rhs  = zeros(nz, nx)
    for i = 1:nz
        for j = 1:nx
            rhs[i,j] = -sin(2*pi/len_z*(i-1)*h) * sin(2*pi/len_x*(j-1)*h)
            coef[i,j] = 1.0 - cos(2*pi/len_z*(i-1)*h) * sin(2*pi/len_x*(j-1)*h) * len_z / (2*pi*rho*G)

            # rhs[i,j] = 2.0*(i-1)*h*exp(-(((i-1)*h)^2) -(((j-1)*h)^2)) * rho * G
            # coef[i,j] = 1.0 + exp(-(((i-1)*h)^2) -(((j-1)*h)^2))
        end
    end

    tf_coef = constant(coef)
    tf_rhs = constant(rhs)
    function scalar_function(m)
        return sum(tanh(poisson_op(tf_coef,m,tf_h,constant(rho*G), constant(0))))
        # return sum(tanh(poisson_op(m,tf_rhs,tf_h, constant(rho*G), constant(0))))
    end

    m_ = tf_rhs
    # m_  = tf_coef
    v_ = 0.01*rand(nz, nx)
    y_ = scalar_function(m_)
    dy_ = gradients(y_, m_)
    ms_ = Array{Any}(undef, 5)
    ys_ = Array{Any}(undef, 5)
    s_ = Array{Any}(undef, 5)
    w_ = Array{Any}(undef, 5)
    gs_ =  @. 1 / 20^(1:5)

    for i = 1:5
        g_ = gs_[i]
        ms_[i] = m_ + g_*v_
        ys_[i] = scalar_function(ms_[i])
        s_[i] = ys_[i] - y_
        w_[i] = s_[i] - g_*sum(v_.*dy_)
    end

    sess = Session()
    init(sess)
    sval_ = run(sess, s_)
    wval_ = run(sess, w_)
    close("all")
    loglog(gs_, abs.(sval_), "*-", label="finite difference")
    loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
    loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
    loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

    plt.gca().invert_xaxis()
    legend()
    xlabel("\$\\gamma\$")
    ylabel("Error")
    savefig("poisson.png"); close("all")
end

@testset "sat_op" begin
    function ave_normal(quantity, m, n)
      aa = sum(quantity)
      return aa/(m*n)
    end

    # TODO:
    # const ALPHA = 0.006323996017182
    const ALPHA = 1.0
    const SRC_CONST = 86400.0
    const K_CONST =  9.869232667160130e-16 * 86400
    nz=20
    nx=30
    sw = constant(zeros(nz, nx))
    swref = constant(zeros(nz,nx))
    μw = constant(0.001)
    μo = constant(0.003)
    K = constant(100.0 .* ones(nz, nx))
    ϕ = constant(0.25 .* ones(nz, nx))
    dt = constant(30.0)
    h = constant(100.0 * 0.3048)
    q1 = zeros(nz,nx)
    q2 = zeros(nz,nx)
    q1[10,5] = 0.002 * (1/(100.0 * 0.3048)^2)/20.0/0.3048 * SRC_CONST
    q2[10,25] = -0.002 * (1/(100.0 * 0.3048)^2)/20.0/0.3048 * SRC_CONST
    qw = constant(q1)
    qo = constant(q2)

    λw = sw.*sw/μw
    λo = (1-sw).*(1-sw)/μo
    λ = λw + λo
    f = λw/λ
    q = qw + qo + λw/(λo+1e-16).*qo

    # Θ = laplacian_op(K.*λo, potential_c, h, constant(0.0))
    Θ = upwlap_op(K*K_CONST, λo, constant(zeros(nz,nx)), h, constant(0.0))

    load_normal = (Θ+q/ALPHA) - ave_normal(Θ+q/ALPHA, nz, nx)

    tf_comp_p0 = upwps_op(K*K_CONST, λ, load_normal, constant(zeros(nz,nx)), h, constant(0.0), constant(2))
    sess = Session()
    init(sess)
    p0 = run(sess, tf_comp_p0)
    tf_p0 = constant(p0)

    # s = sat_op(sw,p0,K,ϕ,qw,qo,sw,dt,h)

    # function step(sw)
    #     λw = sw.*sw
    #     λo = (1-sw).*(1-sw)
    #     λ = λw + λo
    #     f = λw/λ
    #     q = qw + qo + λw/(λo+1e-16).*qo

    #     # Θ = laplacian_op(K.*λo, constant(zeros(nz,nx)), h, constant(0.0))
    #     Θ = upwlap_op(K, λo, constant(zeros(nz,nx)), h, constant(0.0))
    #     # Θ = constant(zeros(nz,nx))

    #     load_normal = (Θ+q/ALPHA) - ave_normal(Θ+q/ALPHA, nz, nx)

    #     p = poisson_op(λ.*K, load_normal, h, constant(0.0), constant(0)) # potential p = pw - ρw*g*h
    #     # p = upwps_op(K, λ, load_normal, constant(zeros(nz,nx)), h, constant(0.0), constant(0))
        # sw = sat_op(sw,p,K,ϕ,qw,qo,μw,μo,sw,dt,h)
    #     return sw
    # end

    # NT=100
    # function evolve(sw, NT, qw, qo)
    #     # qw_arr = constant(qw) # qw: NT x m x n array
    #     # qo_arr = constant(qo)
    #     tf_sw = TensorArray(NT+1)
    #     function condition(i, ta)
    #         tf.less(i, NT+1)
    #     end
    #     function body(i, tf_sw)
    #         sw_local = step(read(tf_sw, i))
    #         i+1, write(tf_sw, i+1, sw_local)
    #     end
    #     tf_sw = write(tf_sw, 1, sw)
    #     i = constant(1, dtype=Int32)
    #     _, out = while_loop(condition, body, [i;tf_sw])
    #     read(out, NT+1)
    # end

    # s = evolve(sw, NT, qw, qo)

    # J = tf.nn.l2_loss(s)
    # tf_grad_K = gradients(J, K)
    # sess = Session()
    # init(sess)
    # # P = run(sess,p0)

    # # error("")
    # S=run(sess, s)
    # imshow(S);colorbar();

    # error("")

    # grad_K = run(sess, tf_grad_K)
    # imshow(grad_K);colorbar();
    # error("")
    # TODO:

    # gradient check -- v
    function scalar_function(m)
        # return sum(tanh(sat_op(m,tf_p0,K*K_CONST,ϕ,qw,qo,μw,μo,constant(zeros(nz,nx)),dt,h)))
        return sum(tanh(sat_op(sw,m,K*K_CONST,ϕ,qw,qo,μw,μo,constant(zeros(nz,nx)),dt,h)))
        # return sum(tanh(sat_op(sw,tf_p0,m,ϕ,qw,qo,μw,μo,constant(zeros(nz,nx)),dt,h)))
        # return sum(tanh(sat_op(sw,tf_p0,K*K_CONST,m,qw,qo,μw,μo,constant(zeros(nz,nx)),dt,h)))
    end

    # m_ = sw
    # v_ = 0.1 * rand(nz,nx)

    m_ = tf_p0
    v_ = 5e5 .* rand(nz,nx)

    # m_ = K*K_CONST
    # v_ = 10 .* rand(nz,nx) *K_CONST

    # m_ = ϕ
    # v_ = 0.1 * rand(nz,nx)

    y_ = scalar_function(m_)
    dy_ = gradients(y_, m_)
    ms_ = Array{Any}(undef, 5)
    ys_ = Array{Any}(undef, 5)
    s_ = Array{Any}(undef, 5)
    w_ = Array{Any}(undef, 5)
    gs_ =  @. 1 / 10^(1:5)

    for i = 1:5
        g_ = gs_[i]
        ms_[i] = m_ + g_*v_
        ys_[i] = scalar_function(ms_[i])
        s_[i] = ys_[i] - y_
        w_[i] = s_[i] - g_*sum(v_.*dy_)
    end

    sess = Session()
    init(sess)
    sval_ = run(sess, s_)
    wval_ = run(sess, w_)
    close("all")
    loglog(gs_, abs.(sval_), "*-", label="finite difference")
    loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
    loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
    loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

    plt.gca().invert_xaxis()
    legend()
    xlabel("\$\\gamma\$")
    ylabel("Error")
    savefig("sat_op.png"); close("all")
end

@testset "upwlap_op" begin
    h = 1.0
    rho = 1000.0
    G = 9.8
    len_z = 16
    len_x = 32
    nz = Int(len_z/h + 1)
    nx = Int(len_x/h + 1)
    tf_h=constant(1.0)

    perm = rand(nz, nx)
    mobi = rand(nz, nx)
    func = rand(nz, nx)

    tf_perm = constant(perm)
    tf_mobi = constant(mobi)
    tf_func = constant(func)

    # gradient check -- v
    function scalar_function(m)
        # return sum(tanh(upwlap_op(m, tf_mobi, tf_func, tf_h, constant(rho*G))))
        # return sum(tanh(upwlap_op(tf_perm, m, tf_func, tf_h, constant(rho*G))))
        return sum(tanh(upwlap_op(tf_perm, tf_mobi, m, tf_h, constant(rho*G))))
    end

    # m_ = constant(rand(10,20))
    # m_ = tf_perm
    # m_ = tf_mobi
    m_ = tf_func
    v_ = rand(nz, nx)
    y_ = scalar_function(m_)
    dy_ = gradients(y_, m_)
    ms_ = Array{Any}(undef, 5)
    ys_ = Array{Any}(undef, 5)
    s_ = Array{Any}(undef, 5)
    w_ = Array{Any}(undef, 5)
    gs_ =  @. 1 / 20^(1:5)

    for i = 1:5
        g_ = gs_[i]
        ms_[i] = m_ + g_*v_
        ys_[i] = scalar_function(ms_[i])
        s_[i] = ys_[i] - y_
        w_[i] = s_[i] - g_*sum(v_.*dy_)
    end

    sess = Session()
    init(sess)
    sval_ = run(sess, s_)
    wval_ = run(sess, w_)
    close("all")
    loglog(gs_, abs.(sval_), "*-", label="finite difference")
    loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
    loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
    loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

    plt.gca().invert_xaxis()
    legend()
    xlabel("\$\\gamma\$")
    ylabel("Error")
    savefig("upwlap_op.png"); close("all")
end

@testset "upwps_op" begin
    h = 20.0
    rho = 1000.0
    G = 9.8
    len_z = 16*20
    len_x = 32*20
    nz = Int(len_z/h + 1)
    nx = Int(len_x/h + 1)
    tf_h=constant(1.0)
    coef = zeros(nz, nx)
    ones_cons = ones(nz, nx)
    rhs  = zeros(nz, nx)
    for i = 1:nz
        for j = 1:nx
            rhs[i,j] = -sin(2*pi/len_z*(i-1)*h) * sin(2*pi/len_x*(j-1)*h)
            coef[i,j] = 1.0 - cos(2*pi/len_z*(i-1)*h) * sin(2*pi/len_x*(j-1)*h) * len_z / (2*pi*rho*G)

            # rhs[i,j] = 2.0*(i-1)*h*exp(-(((i-1)*h)^2) -(((j-1)*h)^2)) * rho * G
            # coef[i,j] = 1.0 + exp(-(((i-1)*h)^2) -(((j-1)*h)^2))
        end
    end
    # coef = 1.0 .+ rand(nz, nx)
    # rhs = rand(nz, nx)

    tf_coef = constant(coef)
    tf_rhs = constant(rhs)
    tf_funcref = constant(rand(nz, nx))
    tf_ones = constant(ones_cons)
    function scalar_function(m)
        # return sum(tanh(upwps_op(tf_coef, tf_ones, m, tf_funcref, tf_h, constant(rho*G), constant(0))))
        # return sum(tanh(upwps_op(m, tf_ones, tf_rhs, tf_ones, tf_h, constant(rho*G), constant(0))))
        return sum(tanh(upwps_op(tf_ones, m, tf_rhs, tf_ones, tf_h, constant(rho*G), constant(0))))
    end

    # m_ = tf_rhs
    m_  = tf_coef
    v_ = 0.01*rand(nz, nx)
    y_ = scalar_function(m_)
    dy_ = gradients(y_, m_)
    ms_ = Array{Any}(undef, 5)
    ys_ = Array{Any}(undef, 5)
    s_ = Array{Any}(undef, 5)
    w_ = Array{Any}(undef, 5)
    gs_ =  @. 1 / 10^(1:5)

    for i = 1:5
        g_ = gs_[i]
        ms_[i] = m_ + g_*v_
        ys_[i] = scalar_function(ms_[i])
        s_[i] = ys_[i] - y_
        w_[i] = s_[i] - g_*sum(v_.*dy_)
    end

    sess = Session()
    init(sess)
    sval_ = run(sess, s_)
    wval_ = run(sess, w_)
    close("all")
    loglog(gs_, abs.(sval_), "*-", label="finite difference")
    loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
    loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
    loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

    plt.gca().invert_xaxis()
    legend()
    xlabel("\$\\gamma\$")
    ylabel("Error")
    savefig("upwps_op.png"); close("all")
end


text:
module CustomTypes

export Parameters, ClassicSS, DualSS, AdjDual, AdjClassic, RegUncSS

mutable struct Parameters
    η
    r
    s
    λ
    α
    m
    σ
    γ
    F
    δ
    b
    ϵ
    Parameters() = new()
end

mutable struct ClassicSS
    zp
    zc
    n
    nᶜ
    n⁰
    u
    θ
    μ
    z_bar
    W
    ξ
    rU
    ClassicSS() = new()
end

mutable struct DualSS
    zp
    zc
    zs
    zf
    npᶜ
    np⁰
    np
    nf
    θ
    rU
    ξ
    μᵖ
    μᶠ
    z_bar_p
    z_bar_f
    z_bar
    W
    u
    DualSS() = new()
end

mutable struct AdjDual
    τ_p
    τ_f
    AdjDual() = new()
end

mutable struct AdjClassic
    τ
    AdjClassic() = new()
end

mutable struct RegUncSS
    F
    zp
    zc
    zs
    zf
    npᶜ
    np⁰
    np
    nf
    u
    q
    θ
    ξ
    μᵖ
    μᶠ
    z_bar_p
    z_bar_f
    z_bar
    W
    Wpᶜ
    Wp⁰
    Wp
    flag
    RegUncSS() = new()
end

end


text:
using BinaryBuilder

name = "NLopt"
version = v"2.6.2"

# Collection of sources required to build NLopt
sources = [
    GitSource("https://github.com/stevengj/nlopt.git",
              "41967f1981d82d8495c0b27151a126bc35ae0d96"), # v2.6.2
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/nlopt
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DNLOPT_PYTHON=Off -DNLOPT_OCTAVE=Off -DNLOPT_MATLAB=Off -DNLOPT_GUILE=Off -DNLOPT_SWIG=Off -DNLOPT_LINK_PYTHON=Off ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms()) # build on all supported platforms

# The products that we will ensure are always built
products = [
    LibraryProduct("libnlopt", :libnlopt),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct UpgradeOperationHistoricalStatusInfoProperties <: SwaggerModel
    runningStatus::Any # spec type: Union{ Nothing, UpgradeOperationHistoryStatus } # spec name: runningStatus
    progress::Any # spec type: Union{ Nothing, RollingUpgradeProgressInfo } # spec name: progress
    error::Any # spec type: Union{ Nothing, ApiError } # spec name: error
    startedBy::Any # spec type: Union{ Nothing, String } # spec name: startedBy
    targetImageReference::Any # spec type: Union{ Nothing, ImageReference } # spec name: targetImageReference
    rollbackInfo::Any # spec type: Union{ Nothing, RollbackStatusInfo } # spec name: rollbackInfo

    function UpgradeOperationHistoricalStatusInfoProperties(;runningStatus=nothing, progress=nothing, error=nothing, startedBy=nothing, targetImageReference=nothing, rollbackInfo=nothing)
        o = new()
        validate_property(UpgradeOperationHistoricalStatusInfoProperties, Symbol("runningStatus"), runningStatus)
        setfield!(o, Symbol("runningStatus"), runningStatus)
        validate_property(UpgradeOperationHistoricalStatusInfoProperties, Symbol("progress"), progress)
        setfield!(o, Symbol("progress"), progress)
        validate_property(UpgradeOperationHistoricalStatusInfoProperties, Symbol("error"), error)
        setfield!(o, Symbol("error"), error)
        validate_property(UpgradeOperationHistoricalStatusInfoProperties, Symbol("startedBy"), startedBy)
        setfield!(o, Symbol("startedBy"), startedBy)
        validate_property(UpgradeOperationHistoricalStatusInfoProperties, Symbol("targetImageReference"), targetImageReference)
        setfield!(o, Symbol("targetImageReference"), targetImageReference)
        validate_property(UpgradeOperationHistoricalStatusInfoProperties, Symbol("rollbackInfo"), rollbackInfo)
        setfield!(o, Symbol("rollbackInfo"), rollbackInfo)
        o
    end
end # type UpgradeOperationHistoricalStatusInfoProperties

const _property_map_UpgradeOperationHistoricalStatusInfoProperties = Dict{Symbol,Symbol}(Symbol("runningStatus")=>Symbol("runningStatus"), Symbol("progress")=>Symbol("progress"), Symbol("error")=>Symbol("error"), Symbol("startedBy")=>Symbol("startedBy"), Symbol("targetImageReference")=>Symbol("targetImageReference"), Symbol("rollbackInfo")=>Symbol("rollbackInfo"))
const _property_types_UpgradeOperationHistoricalStatusInfoProperties = Dict{Symbol,String}(Symbol("runningStatus")=>"UpgradeOperationHistoryStatus", Symbol("progress")=>"RollingUpgradeProgressInfo", Symbol("error")=>"ApiError", Symbol("startedBy")=>"String", Symbol("targetImageReference")=>"ImageReference", Symbol("rollbackInfo")=>"RollbackStatusInfo")
Base.propertynames(::Type{ UpgradeOperationHistoricalStatusInfoProperties }) = collect(keys(_property_map_UpgradeOperationHistoricalStatusInfoProperties))
Swagger.property_type(::Type{ UpgradeOperationHistoricalStatusInfoProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_UpgradeOperationHistoricalStatusInfoProperties[name]))}
Swagger.field_name(::Type{ UpgradeOperationHistoricalStatusInfoProperties }, property_name::Symbol) =  _property_map_UpgradeOperationHistoricalStatusInfoProperties[property_name]

const _allowed_UpgradeOperationHistoricalStatusInfoProperties_startedBy = ["Unknown", "User", "Platform"]

function check_required(o::UpgradeOperationHistoricalStatusInfoProperties)
    true
end

function validate_property(::Type{ UpgradeOperationHistoricalStatusInfoProperties }, name::Symbol, val)
    if name === Symbol("startedBy")
        Swagger.validate_param(name, "UpgradeOperationHistoricalStatusInfoProperties", :enum, val, _allowed_UpgradeOperationHistoricalStatusInfoProperties_startedBy)
    end
end


text:
import Pkg: TOML
import Pkg.Types: VersionSpec

include("Compress.jl")

const registry_path = joinpath(DEPOT_PATH[1], "registries", "General")
const registry_file = joinpath(registry_path, "Registry.toml")
const packages = TOML.parsefile(registry_file)["packages"]

const version_map = Dict{String,Vector{VersionNumber}}()

for (uuid, info) in packages
    path = joinpath(registry_path, info["path"])
    versions_file = joinpath(path, "Versions.toml")
    versions = Compress.load(versions_file)
    version_map[uuid] = sort!(collect(keys(versions)))
end

for (_, info) in packages
    path = joinpath(registry_path, info["path"])
    # load and normalize Deps.toml
    deps_file = joinpath(path, "Deps.toml")
    isfile(deps_file) || continue
    deps = Compress.load(deps_file)
    Compress.save(deps_file, deps)
    # load and normalize Compat.toml
    compat_file = joinpath(path, "Compat.toml")
    isfile(compat_file) || continue
    compat = Compress.load(compat_file)
    for (ver, data) in compat
        for (dep, spec) in data
            ranges = VersionSpec(spec).ranges
            compat[ver][dep] =
                length(ranges) == 1 ? string(ranges[1]) : map(string, ranges)
        end
    end
    Compress.save(compat_file, compat)
end


text:
using Test
using AsyncFinalizers

# Since the world age inside the task is determined by the time when the task is started in
# old Julia, wait a bit for it to start.
sleep(0.01)

include("utils.jl")

CHAOS_COUNTER = Threads.Atomic{Int}(0)

ASYNCFINALIZERSTESTS_CHAOS = get(ENV, "ASYNCFINALIZERSTESTS_CHAOS", "takemanyto")
if ASYNCFINALIZERSTESTS_CHAOS == "takemanyto"
    function AsyncFinalizers.Internal.chaos_takemanyto()
        if Threads.atomic_add!(CHAOS_COUNTER, 1) == 1
            error("!!! CHAOS !!!")
        end
    end
    wait_recovery() = nothing
elseif ASYNCFINALIZERSTESTS_CHAOS == "run_finalizers"
    function AsyncFinalizers.Internal.chaos_run_finalizers()
        # FIXME: a yield point here breaks the test
        # yield()
        if Threads.atomic_add!(CHAOS_COUNTER, 1) == 1
            error("!!! CHAOS !!!")
        end
    end
    CHAOS_RECOVERED = Channel{Nothing}(1)
    function AsyncFinalizers.Internal.queue_reset_done()
        put!(CHAOS_RECOVERED, nothing)
    end
    function wait_recovery()
        if CHAOS_COUNTER[] == 2
            take!(CHAOS_RECOVERED)
        end
    end
else
    error("Unknown target: ASYNCFINALIZERSTESTS_CHAOS = ", ASYNCFINALIZERSTESTS_CHAOS)
end
@info "Testing: ASYNCFINALIZERSTESTS_CHAOS = $ASYNCFINALIZERSTESTS_CHAOS"

AsyncFinalizers.Internal.reinit()

nfinalizers = 5
FINALIZER_CALLED = zeros(Int, nfinalizers)

for i in 1:nfinalizers
    @info "Testing $i-th finalizer"
    function async_finalizer(_)
        FINALIZER_CALLED[i] += 1
    end
    Utils.collect_garbage(; async_finalizer = async_finalizer)
    wait_recovery()
end

@test FINALIZER_CALLED == ones(Int, nfinalizers)


text:
module TestPaths

using ..Utils: runtests_process, ⊜

using Test

const SCRIPT = joinpath(@__DIR__, "../../samples/paths/runtests.jl")

test() = @test runtests_process(SCRIPT) ⊜ true

end  # module


text:
testpath(f) = joinpath(dirname(@__FILE__), f)

for test_file in [
  "base.jl",
  "events.jl",
  "operators.jl",
  "simulations.jl",
  "processes.jl",
  "resources/containers.jl",
  "resources/stores.jl",
  "utils/time.jl",
  ]
  include(testpath(test_file))
end


text:
using ElectrodermalActivity

# Specify the folder containing the raw Empatica data.
input_directory = "/home/zygmunt/Downloads/edadata/EDA LAB DATA/107/experiment"
output_directory = "/home/zygmunt/Downloads/edadata/EDA LAB DATA/107/experiment/processed"
#input_directory = "C:/Users/Spock/Documents/Empatica Data/"
#output_directory = "C:/Users/Spock/Documents/Empatica Data/processed"

# Create the output_directory if it doesn't already exist.
mkpath(output_directory)

#=
 The marked interval starts a specified number of seconds prior to the
 first marker and ends at the first marker.
=#
interval_name_1 = "Interval 1"
seconds_prior_marker_1 = 60

#=
  The second marked interval starts at the first marker and continues for a specified
  number of seconds.
=#
interval_name_2 = "Interval 2"
seconds_after_marker_1 = 60

#=
 The third marked interval starts a specified number of seconds prior to the
 second marker and ends at the second marker.
=#
interval_name_3 = "Interval 3"
seconds_prior_marker_2 = 60

#=
  The fourth marked interval starts at the second marker and continues for a specified
  number of seconds.
=#
interval_name_4 = "Interval 4"
seconds_after_marker_2 = 60

interval_names = (name_1 = interval_name_1,
                  name_2 = interval_name_2,
                  name_3 = interval_name_3,
                  name_4 = interval_name_4)
offsets = (offset_1 = seconds_prior_marker_1,
           offset_2 = seconds_after_marker_1,
           offset_3 = seconds_prior_marker_2,
           offset_4 = seconds_after_marker_2)

a = prepare_empatica(input_directory, output_directory, interval_names, offsets)


text:
u, v, w = ℝ⁴(rand(4)), ℝ⁴(rand(4)), ℝ⁴(rand(4))
a, b = ℝ⁴(rand(4)), ℝ⁴(rand(4))
zero = ℝ⁴(0, 0, 0, 0)
α, β = rand(2)

## Abstract Vector Space Tests ##

@test isapprox(u + (v + w), (u + v) + w) # 1. associativity of addition
@test isapprox(u + v, v + u) # 2. commutativity of addition
@test isapprox(u + zero, u) # 3. the zero vector
@test isapprox(u - u, zero) # 4. the inverse element
@test isapprox(α * (u + v), α * u + α * v) # 5. distributivity Ι
@test isapprox((α + β) * u, α * u + β * u) # 6. distributivity ΙΙ
@test isapprox(α * (β * u), (α * β) * u) # 7. associativity of scalar multiplication
@test isapprox(1u, u) # 8. the unit scalar 1

## Product Space Tests ##

@test isapprox(dot(u, v), dot(v, u)) # symmetric
@test isapprox(dot(u, α * a + β * b), α * dot(u, a) + β * dot(u, b)) # linear
@test dot(u, u) ≥ 0 && isapprox(dot(u, zero), 0) # positive semidefinite
@test size(outer(u, v)) == (length(u), length(v))

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
using HDF5

function loadHDF5(filename="MHD_256_Mag.h5")
    """Load magnetic field from HDF5 data"""
    f = h5open(filename, "r")
    mm = read(f, "maxmin")'
    Bx = read(f, "BX") * ((mm[1, 1] - mm[2, 1]) / 256 + mm[2, 1])
    By = read(f, "BX") * ((mm[1, 2] - mm[2, 2]) / 256 + mm[2, 2])
    Bz = read(f, "BX") * ((mm[1, 3] - mm[2, 3]) / 256 + mm[2, 3])
    close(f)

    return(Bx, By, Bz)
end


text:
## Compatibility layer to work with Julia 1.0

# https://github.com/JuliaLang/julia/pull/29679
@static if VERSION < v"1.1.0-DEV.472"
    isnothing(::Any) = false
    isnothing(::Nothing) = true
    export isnothing
end

# https://github.com/JuliaLang/julia/pull/32148
@static if VERSION < v"1.3.0"
    import Base: print
    print(io::IO, ::Nothing) = print(io, "nothing")
end

# https://github.com/JuliaLang/julia/pull/29173
@static if VERSION < v"1.1.0"
    import Future: copy!
end


text:
module Codegen
using Reexport
@reexport using Cxx

import CUDAnativelib.Version: __sm__, __triple__, CUDA_VERSION, LLVM_VERSION, CUDA_HOME

function new_cuda_compiler(;target = __triple__, sm = __sm__)
  C = Cxx.setup_instance(makeCCompiler = false, target = target, CPU = sm)

  if CUDA_VERSION >= 8000 && LLVM_VERSION < v"4.0.0"
    warn("CUDA 8.0+ requires at least LLVM 4.0. Using the compatibility headers from 4.0")
    warn("Forefully disabling intrinsics from SM_60")
    addHeaderDir(C,
                 joinpath(@__DIR__, "..", "deps", "include", "clang40"),
                 kind = C_ExternCSystem)
  elseif LLVM_VERSION < v"3.8.1"
    error("LLVM is to old. At least version 3.8.1 is required.")
  end

  Cxx.initialize_instance!(C; register_boot=false)
  push!(Cxx.active_instances, C)
  return Cxx.CxxInstance{length(Cxx.active_instances)}()
end

const __current_compiler__ = new_cuda_compiler()
addHeaderDir(__current_compiler__, joinpath(CUDA_HOME, "include"), kind = C_System)

cxxinclude(__current_compiler__, "__clang_cuda_runtime_wrapper.h")

# Cxx.jl needs to have this function defined and this is fine as long as we don't call it.
cxx"""
extern "C" {
extern int __cxxjl_personality_v0();
}
"""
end # module


text:
using Catlab.CategoricalAlgebra
using Catlab.Present
using Catlab.Theories
using Catlab.Graphs.BasicGraphs: TheoryGraph
using Catlab.Graphs
using DataStructures: OrderedDict

"""
Reference: CT for computing science:
https://www.math.mcgill.ca/triples/Barr-Wells-ctcs.pdf

We are concerned with "Regular" sketches, where no node is the vertex of more
than one cone.

Here are also interesting examples + explanation of connection to Essentially
Algebraic Theories: https://www.math.mcgill.ca/barr/papers/sketch.pdf

Section 7.7 describes how this is essentially the same as the
syntactic version. Maybe useful to help convert between.

FD theories add cocones and model sum types. This goes beyond what we want to
model with a GAT, though may be interesting to consider for the model
enumeration project.

Also check Fiore: categorical semantics of dependent type theory
"""

"""
Section 7.3+10.1.3 Shorthand Caveats:
2. nodes a×b×c IMPLICITLY have a cone with projection legs.
3. likewise nodes labeled 1 are implicitly terminal (empty base limit)
4. arrows labeled ⟨f₁,f₂,...⟩:a->b₁×b₂×... , where edges fₙ all share a codomain, are implicitly assumed to have diagrams:
                fᵢ
              a ⟶ bᵢ
 ⟨f₁,f₂,...⟩  ↓   ↗ pᵢ
          b₁×b₂×...

5. arrows labeled f₁×f₂×...: a₁×a₂×...->b₁×b₂×...
are implicitly assumed to have diagrams:
                 pᵢ
      a₁×a₂×... --> aᵢ
f₁×f₂×...↓          ↓ fᵢ
      b₁×b₂×... --> bᵢ
                pᵢ

6. diagrams of the form of two paths with a common start and end can be specified as s₁;s₂;... = t₁;t₂;...
7. If one leg of the above paths is empty, we require an identity arrow and set the other path to be equal to it.
8. Nodes labeled a×ᵧb implies the existence of a cone
        a×ᵧb
       ↙   ↘
      a→ γ ←b
9. an arrow s: a↣b (i.e. s is monic) implies the existence of a cone
        a
   id ╱ | ╲ id
    ↙   ↓s  ↘
   a ⟶ b ⟵ a
     s     s
"""

"""
Combinatorial representation of data of a finite limit sketch

It is not yet clear whether this captures everything we want to capture, e.g.
we expect that the diagrams are trees without cycles so that there are only a
finite number of equations to check, but there could be alternative styles of
diagrams that would work better in a CSet.
"""
@present TheoryFLSketch(FreeSchema) begin
  # Main Graph
  (V, E)::Ob
  (src, tgt)::Hom(E,V)

  # Diagrams
  (Dv, De)::Ob
  root::Hom(V, Dv)
  (dSrc, dTgt)::Hom(De,Dv) # graph data of all diagrams (disjoint union)
  dV::Hom(Dv, V) # Partition union of graphs
  dE::Hom(De, V) # via map into V

  # Cones
  (Cone, Leg, Cv, Ce)::Ob
  (cSrc, cTgt)::Hom(Ce,Cv) # graph data of all cone bases (disjoint union)
  apex::Hom(Cone, V)    # Which object is the cone vertex?
  legCone::Hom(Leg, Cone) # Which Cone does this leg belong to?
  legTgt::Hom(Leg, Cv)
  legEdge::Hom(Leg, E)
  cV::Hom(Cv, Cone)    # Partition cone graph
  cE::Hom(Ce, Cone)    # via map into Cone

  # Homorphisms data
  cvMap::Hom(Cv, V)
  ceMap::Hom(Ce, E)
  dvMap::Hom(Dv, V)
  deMap::Hom(De, E)

  # Diagrams/cones don't touch each other
  compose(cSrc, cV) == cE    # EQUATION ON Ce
  compose(cTgt, cV) == cE    # EQUATION ON Ce
  compose(dSrc, dV) == dE    # EQUATION ON De
  compose(dTgt, dV) == dE    # EQUATION ON De

  # Root of each diagram is the right vertex
  compose(root, dV) == id(V)  # EQUATION ON V
  compose(root, dvMap) == id(V)  # EQUATION ON V

  # Homomorphism properties
  compose(legEdge, src) == compose(legCone, apex) # EQUATION ON CONE LEGS
  compose(legEdge, tgt) == compose(legTgt, cvMap) # EQUATION ON CONE LEGS
  compose(deMap, src) == compose(dSrc, dvMap) # EQUATION ON De
  compose(deMap, tgt) == compose(dTgt, dvMap) # EQUATION ON De
  compose(ceMap, src) == compose(cSrc, cvMap) # EQUATION ON Ce
  compose(ceMap, tgt) == compose(cTgt, cvMap) # EQUATION ON Ce
end;

"""Append either _src or _tgt to a symbol"""
function add_srctgt(fk::Symbol, src::Bool)::Symbol
  return Symbol(string(fk) * "_" * (src ? "src" : "tgt"))
end
"""Append both _src and _tgt to a symbol"""
function add_srctgt(fk::Symbol)::Pair{Symbol, Symbol}
  return add_srctgt(fk, true) => add_srctgt(fk, false)
end

if !isdefined(Main, :Fls) # because we are currently running this as a script
const Fls = CSetType(TheoryFLSketch, index=[
  :src, :tgt, :dSrc, :dV, :dE, :cV, :cE, :apex, :legCone]);
end

"""Edges and vertices labeled by symbols"""
@present TheoryLabeledGraph <: TheoryGraph begin
  Label::Data
  vlabel::AttrType(V,Label)
  elabel::AttrType(E,Label)
end;
# we don't want unique_index=[:vlabel, :elabel] because we want to use
# labels to specify homomorphisms which may not be injective
@acset_type LabeledGraph(TheoryLabeledGraph, index=[:src,:tgt]){Symbol};

"""Annotate a FLS CSet with labels"""
struct LabeledFLS
  fls::Fls
  labels::Pair{Vector{Symbol}, Vector{Symbol}} # label nodes and edges
  function LabeledFLS(f::Fls, l::Pair{Vector{Symbol}, Vector{Symbol}})
    for labs in l
        length(labs)==length(Set(labs)) || error("labels must be unique: $labs")
    end
    length(l[1]) == nparts(f, :V) || error("incorrect # of vlabels $(l[1])")
    length(l[2]) == nparts(f, :E) || error("incorrect # of elabels $(l[2])")
    return new(f, l)
  end
end

"""Access vertex labels of graph underlying FLS"""
function vi(fls::LabeledFLS, v::Symbol)::Int
  return findfirst(==(v), fls.labels[1])
end

"""Access edge labels of graph underlying FLS"""
function ei(fls::LabeledFLS, v::Symbol)::Int
  return findfirst(==(v), fls.labels[2])
end

"""No cones, trivial diagrams"""
function FLSinit(g::LabeledGraph)::LabeledFLS
  fls = Fls()
  add_parts!(fls, :V, nv(g))
  add_parts!(fls, :Dv, nv(g), dV=1:nv(g), dvMap=1:nv(g))
  add_parts!(fls, :E, ne(g))
  set_subpart!(fls, :src, g[:src])
  set_subpart!(fls, :tgt, g[:tgt])
  set_subpart!(fls, :root, 1:nv(g))
  return LabeledFLS(fls, g[:vlabel] => g[:elabel])
end

"""Get the commutivity diagram of a FLS"""
function get_diagram(fls::LabeledFLS, root::Symbol)::LabeledGraph
  res = LabeledGraph()
  f = fls.fls
  rooti = vi(fls, root)
  nodes = f.indices[:dV][rooti]
  edges = f.indices[:dE][rooti]
  nodedict = Dict([i=>e for (e, i) in enumerate(nodes)])
  vlab = [fls.labels[1][f[:dvMap][i]] for i in nodes]
  add_parts!(res, :V, length(nodes), vlabel=vlab)
  edata = [(e, f[:dSrc][e], f[:dTgt][e]) for e in edges]
  for (e, s, t) in edata
    elab = fls.labels[2][f[:deMap][e]]
    add_part!(res, :E, src=nodedict[s], tgt=nodedict[t], elabel=elab)
  end
  return res
end

"""Get whole cone (apex last index), or just the base of of the cone"""
function get_cone(fls::LabeledFLS, apex::Symbol, base::Bool=false
         )::Union{Nothing, LabeledGraph}
  f = fls.fls
  apexes = fls.labels[1][f[:apex]]
  cone_ind = findfirst(==(apex), apexes)
  if cone_ind === nothing
    return nothing
  end
  res = LabeledGraph()
  cv, ce = [f.indices[x][cone_ind] for x in [:cV, :cE]]
  reind = [findfirst(==(i), cv) for i in 1:(isempty(cv) ? 0 : maximum(cv))]
  vlab = fls.labels[1][f[:cvMap][cv]]
  add_parts!(res, :V, length(cv), vlabel=vlab)

  elab = fls.labels[2][f[:ceMap][ce]]
  es = reind[f[:cSrc][ce]]
  et = reind[f[:cTgt][ce]]
  add_parts!(res, :E, length(ce), elabel=elab, src=es, tgt=et)
  if !base
    pex = add_part!(res, :V, vlabel=fls.labels[1][f[:apex][cone_ind]])
    legs = f.indices[:legCone][cone_ind]
    el = fls.labels[2][f[:legEdge][legs]]
    tg = reind[f[:legTgt][legs]]
    add_parts!(res, :E, length(legs), elabel=el, src=pex, tgt=tg)
  end
  return res
end

"""Extract the underlying graph (no diagrams/cones)"""
function get_schema(fls::LabeledFLS)::LabeledGraph
  res = LabeledGraph()
  add_parts!(res, :V, length(fls.labels[1]), vlabel=fls.labels[1])
  add_parts!(res, :E, length(fls.labels[2]), src=fls.fls[:src], tgt=fls.fls[:tgt], elabel=fls.labels[2])
  return res
end

"""Get all commuting paths starting at a particular node"""
function all_paths(fls, root::Symbol)::OrderedDict{Vector{Symbol}, Int}
  res = OrderedDict{Vector{Symbol}, Int}()
  rooti = vi(fls, root)
  f = fls.fls
  rootind = f[:root][rooti]
  stack = Tuple{Int,Vector{Symbol}}[(rootind, [])]
  seen = Set()
  while !isempty(stack)
    currnode, currpath = pop!(stack)
    res[currpath] = currnode
    if !(currnode in seen)
      push!(seen, currnode)
      for e in f.indices[:dSrc][currnode]
        orig_e = f[:deMap][e]
        e_symbol = fls.labels[2][orig_e]
        nextnode = f[:dTgt][e]
        push!(stack, (nextnode, vcat(currpath, [e_symbol])))
      end
    end
  end

  return res
end

"""
Add a diagram where the first node is the root. Labels indicate homomorphism.
"""
function add_diagram!(fls::LabeledFLS, d::LabeledGraph)::Nothing
  rootless = LabeledGraph()
  add_parts!(rootless, :V, nv(d), vlabel=d[:vlabel])
  for (e, (ss, tt)) in enumerate(zip(d[:src], d[:tgt]))
    if tt != 1
      add_part!(rootless, :E, elabel=d[:elabel][e], src=ss, tgt=tt)
    end
  end
  topological_sort(rootless)  # must be acyclic, after removing loops to root
  root = d[:vlabel][1]
  rooti = vi(fls, root)
  # paths that already exist in the FLS
  location = all_paths(fls, root)
  f = fls.fls
  seen = Dict{Int, Int}() # index in d => index in FLS diagram
  # (index in d, index in FLS diagram, path in d)
  stack = Tuple{Int,Int,Vector{Symbol}}[(1, f[:root][rooti], Symbol[])]
  while !isempty(stack)
    dnode, currnode, currpath = pop!(stack)
    seen[dnode] = currnode
    vlabel = d[:vlabel][dnode]
    nextedges = d.indices[:src][dnode]
    for e in nextedges
      elab = d[:elabel][e]
      orig_edge = ei(fls, elab)
      nextpath = vcat(currpath, [elab])
      seenflag = false
      if haskey(location, nextpath)
        nextnode = location[nextpath]
      else
        next_d = d[:tgt][e]
        next_orig = vi(fls, d[:vlabel][next_d])
        if haskey(seen, next_d)
          nextnode = seen[next_d]
          seenflag=true
        else
          nextnode = add_part!(f, :Dv, dV=rooti, dvMap=next_orig)
        end
        add_part!(f, :De, dSrc=currnode, dTgt=nextnode,
        dE=rooti, deMap=orig_edge)
      end
      if !seenflag
        push!(stack, (d[:tgt][e], nextnode, nextpath))
      end
    end
  end
end
"""
Add a cone, assuming the last vertex of the graph is the apex
"""
function add_cone!(fls::LabeledFLS, c::LabeledGraph)::Nothing
  f = fls.fls
  conetgt = vi(fls, c[:vlabel][end])
  prevcones = vcat(f.indices[:apex]...)
  @assert !(conetgt in prevcones) "No more than one cone on a given vertex"
  cone_id = add_part!(f, :Cone)
  vert_ids = add_parts!(f, :Cv, nv(c)-1, cV=cone_id,
              cvMap=[vi(fls, l) for l in c[:vlabel][1:end-1]])
  nonlegs = [i for (i, t) in enumerate(c[:src]) if t!=nv(c)]
  legs  = [i for (i, t) in enumerate(c[:src]) if t==nv(c)]
  emap, lemap = [[ei(fls, l) for l in c[:elabel][x]]
           for x in [nonlegs, legs]]
  esrc = vert_ids[c[:src][nonlegs]]
  etgt, letgt = [vert_ids[c[:tgt][x]] for x in [nonlegs, legs]]
  add_parts!(f, :Ce, length(nonlegs), cE=cone_id, ceMap=emap, cSrc=esrc, cTgt=etgt)
  set_subpart!(f, cone_id, :apex, vi(fls, c[:vlabel][nv(c)]))
  add_parts!(f, :Leg, length(legs), legCone=cone_id,
         legEdge=lemap, legTgt=letgt)
  return nothing
end

"""
Create a linear sketch (no cones) based on a CSet presentation
"""
function catpres_to_linear(pres::Presentation)::LabeledFLS
  obs = pres.generators[:Ob]
  homs = pres.generators[:Hom]
  vlab = [x.args[1] for x in obs]
  odict = Dict([o => i for (i, o) in enumerate(obs)])
  idsymbs = [Symbol("_id_"*string(v)) for v in vlab]
  elab = vcat([x.args[1] for x in homs], idsymbs)
  # Create initial graph
  g = LabeledGraph()
  n, nh = length(vlab), length(homs)
  add_parts!(g, :V, n, vlabel=vlab)
  srcs = vcat([odict[h.type_args[1]] for h in homs], 1:n)
  tgts = vcat([odict[h.type_args[2]] for h in homs], 1:n)
  add_parts!(g, :E, nh+n, src=srcs, tgt=tgts, elabel=elab)
  fls = FLSinit(g)
  # Add diagrams
  for (p1, p2) in pres.equations
    d = paths_to_diagram(p1, p2)
    add_diagram!(fls, d)
  end
  return fls
end

"""Create a diagram from a path equality"""
function paths_to_diagram(p1, p2)::LabeledGraph
  d = LabeledGraph()
  root = p1.type_args[1].args[1]
  add_part!(d, :V, vlabel=root)
  function add_path!(start::Bool, p)::Nothing
    typ = typeof(p).parameters[1]
    if typ == :id
      if start
        tgt = add_part!(d, :V, vlabel=root)
      else
        tgt = nparts(d, :V)
      end
      elab = Symbol("_id_"*string(root))
      add_part!(d, :E, src=1, tgt=tgt, elabel=elab)
    elseif typ == :compose
      homs = p.args
      curr = 1
      if start
        last = length(homs) + 1
      else
        last = nparts(d, :V)
      end
      for (i, hom) in enumerate(homs)
        add_part!(d, :V, vlabel=hom.type_args[2].args[1])
        if i == length(homs)
          tgt = last
        else
          tgt = curr+1
        end
        add_part!(d, :E, src=curr, tgt=tgt, elabel=hom.args[1])
        curr += 1
      end
    elseif typ == :generator
      if start
        last = add_part!(d, :V, vlabel=p.type_args[2].args[1])
      else
        last = nparts(d, :V)
      end
      add_part!(d, :E, src=1, tgt = last, elabel=p.args[1])
    else
      @assert false "Not prepared for type $typ"
    end
    return nothing
  end
  add_path!(true, p1)
  add_path!(false, p2)
  return d
end


text:
using Documenter, MixedModels

makedocs()

deploydocs(
    deps = Deps.pip("pygments", "mkdocs", "mkdocs-material"),
    repo = "github.com/dmbates/MixedModels.jl.git",
    julia = "release",
    osname = "linux"
)

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
@testset "Reader" begin

function single_field_doc_(type::UInt8, value, endian_convert = true)
    io = IOBuffer()
    write(io, htol(Int32(4 + 1 + 2 + sizeof(value) + 1)))
    write(io, type)
    write(io, UInt8('x'))
    write(io, 0x0)
    if isbits(value)
        unsafe_write(io, Ref(endian_convert ? htol(value) : value), sizeof(value))
    else
        write(io, value)
    end
    write(io, 0x0)
    take!(io)
end

@testset "double" begin
    reader = BSONReader(single_field_doc_(BSON_TYPE_DOUBLE, 1.25))
    @test reader["x"][Float64] == 1.25
    @test reader["x"][Number] == 1.25
    @test reader["x"][AbstractFloat] == 1.25
    @test reader["x"][Any] == 1.25
end

@testset "decimal" begin
    reader = BSONReader(single_field_doc_(BSON_TYPE_DECIMAL128, d128"1.25"))
    @test reader["x"][Dec128] == d128"1.25"
    @test reader["x"][Number] == d128"1.25"
    @test reader["x"][AbstractFloat] == d128"1.25"
    @test reader["x"][Any] == d128"1.25"
end

@testset "int32" begin
    reader = BSONReader(single_field_doc_(BSON_TYPE_INT32, Int32(123)))
    @test reader["x"][Int32] == 123
    @test reader["x"][Number] == 123
    @test reader["x"][Integer] == 123
    @test reader["x"][Any] == 123
end

@testset "int64" begin
    reader = BSONReader(single_field_doc_(BSON_TYPE_INT64, Int64(123)))
    @test reader["x"][Int64] == 123
    @test reader["x"][Number] == 123
    @test reader["x"][Integer] == 123
    @test reader["x"][Any] == 123
end

@testset "bool" begin
    reader = BSONReader(single_field_doc_(BSON_TYPE_BOOL, 0x1))
    @test reader["x"][Bool] == true
    @test reader["x"][Any] == true
    reader = BSONReader(single_field_doc_(BSON_TYPE_BOOL, 0x0))
    @test reader["x"][Bool] == false
    @test reader["x"][Any] == false
end

@testset "datetime" begin
    t = DateTime(2021, 1, 2, 9, 30)
    v = trunc(Int64, datetime2unix(t) * 1000)
    reader = BSONReader(single_field_doc_(BSON_TYPE_DATETIME, v))
    @test reader["x"][DateTime] == t
    @test reader["x"][Any] == t
end

@testset "timestamp" begin
    x = BSONTimestamp(1, 2)
    reader = BSONReader(single_field_doc_(BSON_TYPE_TIMESTAMP, x))
    @test reader["x"][BSONTimestamp] == x
    @test reader["x"][Any] == x
end

@testset "ObjectId" begin
    x = BSONObjectId((
        0x1, 0x2, 0x3, 0x4,
        0x5, 0x6, 0x7, 0x8,
        0x9, 0xA, 0xB, 0xC,
    ))
    reader = BSONReader(single_field_doc_(BSON_TYPE_OBJECTID, x, false))
    @test reader["x"][BSONObjectId] == x
    @test reader["x"][Any] == x
end

@testset "string" begin
    x = "test"
    io = IOBuffer()
    write(io, htol(Int32(length(x) + 1)))
    write(io, x)
    write(io, 0x0)
    buf = take!(io)
    reader = BSONReader(single_field_doc_(BSON_TYPE_STRING, buf))
    @test reader["x"][String] == x
    @test reader["x"][Any] == x
    reader = BSONReader(single_field_doc_(BSON_TYPE_CODE, buf))
    @test reader["x"][String] == x
    @test reader["x"][Any] == BSONCode(x)
end

@testset "binary" begin
    x = rand(UInt8, 10)
    io = IOBuffer()
    write(io, htol(Int32(length(x))))
    write(io, BSON_SUBTYPE_GENERIC_BINARY)
    write(io, x)
    reader = BSONReader(single_field_doc_(BSON_TYPE_BINARY, take!(io)))
    x2 = reader["x"][BSONBinary]
    @test x2.data == x
    @test x2.subtype == BSON_SUBTYPE_GENERIC_BINARY
    x2 = reader["x"][UnsafeBSONBinary]
    @test x2.data == x
    @test x2.subtype == BSON_SUBTYPE_GENERIC_BINARY
    x2 = reader["x"][Any]
    @test x2 isa BSONBinary
    @test x2.data == x
    @test x2.subtype == BSON_SUBTYPE_GENERIC_BINARY
end

@testset "uuid" begin
    x = uuid4()
    io = IOBuffer()
    write(io, htol(Int32(16)))
    write(io, BSON_SUBTYPE_UUID)
    unsafe_write(io, Ref(x), 16)
    reader = BSONReader(single_field_doc_(BSON_TYPE_BINARY, take!(io)))
    @test reader["x"][UUID] == x
    @test reader["x"][Any] == x
end

@testset "heterogenous array" begin
    io = IOBuffer()
    len = Int32(4 + 11 + 11 + 4 + 1)
    write(io, htol(len))
    write(io, BSON_TYPE_DOUBLE)
    write(io, UInt8('1'))
    write(io, 0x0)
    write(io, 1.25)
    write(io, BSON_TYPE_INT64)
    write(io, UInt8('2'))
    write(io, 0x0)
    write(io, Int64(123))
    write(io, BSON_TYPE_BOOL)
    write(io, UInt8('3'))
    write(io, 0x0)
    write(io, 0x1)
    write(io, 0x0)
    buf = take!(io)
    @assert len == length(buf)
    reader = BSONReader(single_field_doc_(BSON_TYPE_ARRAY, buf))
    @test reader["x"][Vector{Any}] == [1.25, Int64(123), true]
    @test reader["x"][Any] == [1.25, Int64(123), true]
end

@testset "homogenous array" begin
    io = IOBuffer()
    len = Int32(4 + 11 * 3 + 1)
    write(io, htol(len))
    write(io, BSON_TYPE_INT64)
    write(io, UInt8('1'))
    write(io, 0x0)
    write(io, Int64(1))
    write(io, BSON_TYPE_INT64)
    write(io, UInt8('2'))
    write(io, 0x0)
    write(io, Int64(2))
    write(io, BSON_TYPE_INT64)
    write(io, UInt8('3'))
    write(io, 0x0)
    write(io, Int64(3))
    write(io, 0x0)
    buf = take!(io)
    @assert len == length(buf)
    reader = BSONReader(single_field_doc_(BSON_TYPE_ARRAY, buf))
    @test reader["x"][Vector{Int64}] == Int64[1, 2, 3]
    @test reader["x"][1][Int64] == 1
    @test reader["x"][3][Int64] == 3
end

@testset "regex" begin
    io = IOBuffer()
    write(io, "test")
    write(io, 0x0)
    write(io, "abc")
    write(io, 0x0)
    reader = BSONReader(single_field_doc_(BSON_TYPE_REGEX, take!(io)))
    x = reader["x"][BSONRegex]
    @test x.pattern == "test"
    @test x.options == "abc"
    x = reader["x"][Any]
    @test x.pattern == "test"
    @test x.options == "abc"
end

@testset "conversion error $T" for T in [
    Float64,
    Int64,
    Int32,
    Bool,
    DateTime,
    Dec128,
    UUID,
    String,
    Nothing,
    BSONTimestamp,
    BSONObjectId,
    BSONBinary,
    UnsafeBSONBinary,
    BSONRegex,
    BSONCode,
    BSONSymbol,
    BSONMinKey,
    BSONMaxKey,
    BSONUndefined,
    BSONDBPointer,
    BSONUUIDOld
]
    buf = UInt8[]
    writer = BSONWriter(buf)
    writer["x"] = T == Int64 ? 1.25 : Int64(123)
    close(writer)
    reader = BSONReader(buf)
    @test_throws BSONConversionError(T == Int64 ? BSON_TYPE_DOUBLE : BSON_TYPE_INT64, T) reader["x"][T]
end

@testset "long field name skip" begin
    buf = UInt8[]
    writer = BSONWriter(buf)
    writer["some_very_long_field_name"] = 123
    writer["x"] = 456
    close(writer)
    reader = BSONReader(buf)
    @test reader["x"][Int] == 456
end

@testset "missing field" begin
    buf = UInt8[]
    writer = BSONWriter(buf)
    writer["x"] = 123
    close(writer)
    reader = BSONReader(buf)
    @test_throws KeyError reader["y"]
end

@testset "foreach" begin
    buf = UInt8[]
    writer = BSONWriter(buf)
    writer["x"] = 1
    writer["y"] = 2
    writer["z"] = 3
    close(writer)
    reader = BSONReader(buf)
    values = Pair{String, Int}[]
    foreach(x -> push!(values, x.first => x.second[Int]), reader)
    @test values == ["x" => 1, "y" => 2, "z" => 3]
end

@testset "optional" begin
    buf = UInt8[]
    writer = BSONWriter(buf)
    writer["x"] = nothing
    writer["y"] = 123
    close(writer)
    reader = BSONReader(buf)
    @test reader["x"][Union{Nothing, Int}] === nothing
    @test reader["y"][Union{Nothing, Int}] === 123
end

end


text:
@testset "RandomVariable" begin
    dist = Normal(0, 1)
    name = :x
    x = RandomVariable(dist, name)

    @test isa(x, RandomVariable)
    @test x.dist == dist
    @test x.name == name
    @test dimensions(x) == 1

    @test size(sample(x)) == (1, 1)
    @test size(sample(x, 10)) == (10, 1)
end

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
using Distributions
using JLD
using ArgParse
using NCDatasets
using LinearAlgebra
#using DelimitedFiles
using CSV
using DataFrames
# Import EnsembleKalmanProcesses modules
using EnsembleKalmanProcesses.EnsembleKalmanProcessModule
using EnsembleKalmanProcesses.Observations
using EnsembleKalmanProcesses.ParameterDistributionStorage

include(joinpath(@__DIR__, "helper_funcs.jl"))

"""
ek_update(iteration_::Int64)

Update MIMAParameters ensemble using Ensemble Kalman Inversion,
return the ensemble and the parameter names.
"""
function ek_update(iteration_::Int64)
    # Recover versions from last iteration
    versions = readlines("versions_$(iteration_).txt")

    # Use known u_names and n_params, could eventually set this up by reading from output file (as done in ClimateMachine example)
    u_names = ["cwtropics"]
    n_params = length(u_names)

    # Recover ensemble from last iteration, [N_ens, N_params]
    N = length(versions)
    u = zeros(N, n_params)

    # Get ensemble parameter inputs from txt file
    basedir = "/scratch/users/lauraman/MiMA/ekp_runs/"
    filename = "QBO_metrics.csv"

    for (ens_index, version_) in enumerate(versions)
        input_filename = "$(basedir)/$(version_)"
        open(input_filename, "r") do io
            line = readline(io)
            for (u_index, u_name) in enumerate(u_names)
                u[ens_index, u_index] = parse(Float64, (strip(string(split(line, "$(u_name) = ")[2]), [' '])) )
            end
        end
    end
    u = Array(u')
    println(u)

    println("Get truth")
    # Get observations / truth output
    y_names = ["std20", "std77"]
    yt = zeros(0)
    yt_var_list = []

    truthdir = "/scratch/users/lauraman/MiMA/runs/033/"
    runname = "$(truthdir)/$(filename)"
    df = DataFrame(CSV.File(runname))
    yt_ = zeros(length(y_names))
    for (y_index, y_name) in enumerate(y_names)
        yt_[y_index] = df[1, y_name]
    end
    println(yt_)

    # Also get variance for truth
    y_names_var = ["std20_var", "std77_var"]
    yt_var_ = zeros(length(y_names_var))
    for (y_index, y_name) in enumerate(y_names_var)
        yt_var_[y_index] = df[1, y_name]
    end
    # Convert to diagonal covariance matrix
    yt_var_ = diagm(yt_var_)
    println(yt_var_)

    # Add nugget to variance (regularization)
    println(Matrix(0.1I, size(yt_var_)[1], size(yt_var_)[2]))
    yt_var_ = yt_var_ + Matrix(0.1I, size(yt_var_)[1], size(yt_var_)[2])
    println(yt_var_)
    println(yt_var_list)
    append!(yt, yt_)
    push!(yt_var_list, yt_var_)

    # Get outputs from .csv file
    println("Get outputs")
    g_names = y_names
    g_ens = zeros(N, length(g_names))
    basedir = "/scratch/users/lauraman/MiMA/ekp_runs/"
    filename = "QBO_metrics.csv"
    for (ens_index, version_) in enumerate(versions)
        runname = "$(basedir)/$(iteration_)_$(ens_index)/$(filename)"
        df = DataFrame(CSV.File(runname))
        colnames = names(df)
        for (g_index, g_name) in enumerate(g_names)
            g_ens[ens_index, g_index] = df[1, g_name]
        end

    end
    g_ens = Array(g_ens')
    print(g_ens)

    # Construct EKP
    ekobj = EnsembleKalmanProcess(u, yt_, yt_var_, Inversion())
    println(ekobj)
    # Advance EKP
    update_ensemble!(ekobj, g_ens)
    # Get new step
    u_new = get_u_final(ekobj)
    println(u_new)
    return u_new, u_names
end

# Read iteration number of ensemble to be recovered
s = ArgParseSettings()
@add_arg_table s begin
    "--iteration"
    help = "Calibration iteration number"
    arg_type = Int
    default = 1
end
parsed_args = parse_args(ARGS, s)
iteration_ = parsed_args["iteration"]

println(iteration_)
# Perform update
ens_new, param_names = ek_update(iteration_)
params_arr = [col[:] for col in eachcol(ens_new)]

# Generate new identifiers
versions = map(param -> generate_mima_params(param, param_names), params_arr)
open("versions_$(iteration_+1).txt", "w") do io
    for version in versions
        write(io, "mima_param_defs_$(version).txt\n")
    end
end


text:
#!/usr/bin/env julia

# test polygon simplification

using Luxor

using Test

using Random
Random.seed!(42)

function sinecurves()
    gsave()
    translate(100, 100)
    setline(0.25)
    sethue(0, 0, 0)
    x_vals = collect(0:pi/180: 40pi)
    # generate array of Points
    plist = [Point(d * 2pi , -sin(d) * cos(20 * d) * 8 * sin(d/10)) for d in x_vals]
    prettypoly(plist, :stroke, () -> (gsave(); randomhue(); circle(O, 0.5, :stroke); grestore()))
    text("original " * string(length(plist)) * " vertices" , 0, -15)
    for detail in [0.01, 0.05, 0.1, 0.2, 0.5, 0.75, 1.0, 2.0]
        translate(0, 50)
        simplified = simplify(plist, detail)
        prettypoly(simplified, :stroke, () -> (gsave(); randomhue(); circle(O, 0.5, :stroke); grestore()))
        text("detail  $(detail), " * string(length(simplified)) * " vertices" , 0, -15)
    end
    grestore()
end

function test(pagewidth, pageheight)
    translate(100, pageheight/2)
    g = GridRect(O, 0, 60)
    setline(0.25)
    sethue(0, 0, 0)
    x_vals = collect(0:pi/100: 4pi)
    polyline = [Point(d * 20pi, 5 * -sin(d) * cos(12 * d) * 8 * sin(d/10)) for d in x_vals]
    poly(polyline, :stroke)
    prettypoly(polyline)
    text("original " * string(length(polyline)) * " vertices" , 0, -30)
    for detail in [0.01, 0.05, 0.075, 0.1, 0.2, 0.5, 0.75, 1.0, 2.0]
        gsave()
        translate(nextgridpoint(g))
        polysimple = simplify(polyline, detail)
        prettypoly(polysimple, :stroke, () -> (gsave(); randomhue(); circle(O, 2, :stroke); grestore()))
        text(" detail ($detail), " * string(length(polysimple)) * " vertices" , 0, -20)
        grestore()
    end
end

function simplify_poly(fname)
    pagewidth  = 1190.0 # points
    pageheight = 1684.0 # points
    Drawing(pagewidth, pageheight, fname)
    sinecurves()
    test(pagewidth, pageheight)
    @test finish() == true
    println("...finished test: output in $(fname)")
end

function simplify_square()
    poly = Point[Point(-1.0, 0.0), Point(-1.0, -1.0), Point(0.0, -1.0), Point(1.0, -1.0), Point(1.0, 0.0), Point(1.0, 1.0), Point(0.0, 1.0), Point(-1.0, 1.0)]

    poly = simplify(poly, 0.01)
    @test poly == Point[Point(-1.0, 0.0), Point(-1.0, -1.0), Point(1.0, -1.0), Point(1.0, 1.0), Point(-1.0, 1.0)]
end

simplify_poly("simplify-poly.pdf")

simplify_square()


text:
@model function simple_model(x::Vector{Float64})
    loc = sample(Normal())
    log_scale = sample(Normal())
    data = plate(observe, Normal(loc, exp(log_scale)), x)
end

@testset "model expansion 1" begin
    results = inference(
        simple_model, LikelihoodWeighting();
        params = (3.0 .+ randn(10),), inference_params = Dict("num_iterations" => 1000)
    )
    post_loc = mean(results, "loc")
    @info "Posterior location estimate = $post_loc"
end

@model function g(x::Float64)
    x = track(x)
end

@testset "model expansion 2" begin
    t = trace()
    g(t, 3.0)
    @info t
end

function string_track(t::Trace, s::String)
    track(t, "test_track_string", s)
end

@testset "track strings" begin
    t = trace()
    string_track(t, "hi")
    @info t
end


text:
using AdventOfCodeSolutions
using Test

function input(puzzle::Puzzle{2020, 22, n}) where n
    io = openInput(puzzle)
    split(read(io, String), "\n", keepempty=false)
end

function solve(::Puzzle{2020, 22, 1}, input)
    nothing
end

function solve(::Puzzle{2020, 22, 2}, input)
    nothing
end

@testset "Day 22" begin
    TEST_INPUT_1 = """
"""

    @test_skip(solve(Puzzle(22, 1), TEST_INPUT_1) == nothing)
    @test_skip(solve(Puzzle(22, 1), input(Puzzle(22, 1))) == nothing)

    TEST_INPUT_2 = """
"""

    @test_skip(solve(Puzzle(22, 2), TEST_INPUT_2) == nothing)
    @test_skip(solve(Puzzle(22, 2), input(Puzzle(22, 2))) == nothing)
end


text:
module TestCore

using Test
using Try

function test_ok_nothing()
    @test Ok() === Ok(nothing)::Ok{Nothing}
end

end  # module


text:
module PenaltyFunctions

using LinearAlgebra, InteractiveUtils, Reexport, RecipesBase
@reexport using LearnBase
import LearnBase: prox, prox!, deriv, value, grad, grad!, addgrad!, scaled

export
    Penalty,
        ElementPenalty,
            ProxableElementPenalty,
                NoPenalty,
                L1Penalty,
                L2Penalty,
                ElasticNetPenalty,
            SCADPenalty,
            MCPPenalty,
            LogPenalty,
        ArrayPenalty,
            NuclearNormPenalty,
            GroupLassoPenalty,
            MahalanobisPenalty,
    addgrad

const AA{T, N} = AbstractArray{T, N}

# common functions
soft_thresh(x::Number, λ::Number) = sign(x) * max(zero(x), abs(x) - λ)

function soft_thresh!(x::AA{<:Number}, λ::Number)
    for i in eachindex(x)
        @inbounds x[i] = soft_thresh(x[i], λ)
    end
    x
end

function name(p::Penalty)
    s = replace(string(typeof(p)), "PenaltyFunctions." => "")
    # s = replace(s, r"\{(.*)", "")
    f = fieldnames(typeof(p))
    flength = length(f)
    if flength > 0
        s *= "("
        for (i, field) in enumerate(f)
            s *= "$field = $(getfield(p, field))"
            if i < flength
                s *= ", "
            end
        end
        s *= ")"
    end
    s
end
Base.show(io::IO, p::Penalty) = print(io, name(p))

include("elementpenalty.jl")
include("arraypenalty.jl")

# Make Penalties Callable
for T in filter(!isabstracttype, union(subtypes(ElementPenalty),
                                      subtypes(ProxableElementPenalty),
                                      subtypes(ArrayPenalty)))
    @eval (pen::$T)(args...) = value(pen, args...)
end

end


text:
using CausalityTools
using TransferEntropy, PerronFrobenius, StateSpaceReconstruction, Simplices
using CrossMappings
using TimeseriesSurrogates
using Documenter, DocumenterMarkdown
using PyCall, Conda
Conda.add("scipy")
using Plots
using DynamicalSystems
using Distributions
using StaticArrays
using Statistics
using StatsBase

ENV["GKSwstype"] = "100"

PAGES = [
    "index.md",
    "Discretization" => "discretization.md",

    #"StateSpaceReconstruction" => [
        #"Delay embeddings" => "embed.md",
        #"Rectangular partitioning" => "partitioning_rectangular.md",
        #"Triangulation partitioning" => "partitioning_triangulation.md"
    #],

    "PerronFrobenius" => [
        "Estimator documentation" => [
            "Transfer operator grid docs" => "transferoperator/transferoperator_grid_docs.md",
        #    "Transfer operator triangulation docs" => "transferoperator/transferoperator_triang_docs.md"
        ],
        #"Examples" => [
        #    "Transfer operator grid example" => "transferoperator/transferoperator_grid_example.md",
        #    "Transfer operator triangulation" => "transferoperator/transferoperator_triang_example.md"
        #],

        "Invariant measure" => "invariantmeasure/invariantmeasure_docs.md",#[
            #"Documentation" => ,
            #"Examples" => "invariantmeasure/invariantmeasure_example.md"
        #]
    ],

    "Causality algorithms" => [
        "TransferEntropy" => [
            #"Common interface" => "transferentropy/commoninterface_TE.md",
            "kNN estimator" => "transferentropy/estimator_TE_kNN.md",
            "Transfer operator grid estimator" => "transferentropy/estimator_TE_transferoperator_grid.md",
            "Visitation frequency estimator" => "transferentropy/estimator_TE_visitfreq.md",
            "Wrapper" => "transferentropy/wrapper_TE.md"
        ],
        "Cross mappings" => [
            "Convergent cross mapping" => "crossmappings/crossmapping.md"
        ]
    ],

    "Surrogate data" => [
        "The method of surrogate data" => "surrogates/surrogates_overview.md",
        "Docs" => [
            "randomshuffle" => "surrogates/randomshuffle_docs.md",
            "randomamplitudes" => "surrogates/randomamplitudes_docs.md",
            "randomphases" => "surrogates/randomphases_docs.md",
            "aaft" => "surrogates/aaft_docs.md",
            "iaaft" => "surrogates/iaaft_docs.md"
        ],
        "Examples" => [
            "randomshuffle" => "surrogates/randomshuffle_example.md",
            "randomamplitudes" => "surrogates/randomamplitudes_example.md",
            "randomphases" => "surrogates/randomphases_example.md",
            "aaft" => "surrogates/aaft_example.md",
            "iaaft" => "surrogates/iaaft_example.md"
        ]
    ],

    #"Workflow" => "workflow.md",

    #"Interop with DynamicalSystems.jl" => [
    #    "Example 1" => "interop_dynamicalsystems_infoflow1.md"
    #],

    #"Examples" => [
    #    "TE partitioning schemes"=> "transferentropy/examples_TE_different_partitionings.md"
    #],

    "Example systems" => [
		"ar1" => "example_systems/ar1.md",
        "henon2" => "example_systems/henon2.md",
        "logistic2" => "example_systems/logistic2.md",
        "logistic3" => "example_systems/logistic3.md",
		"linearmap3d_nonlinearcoupling" => "example_systems/linearmap3d_nonlinearcoupling.md",
		"verdes" => "example_systems/verdes.md"
    ]
]

makedocs(
    sitename = "CausalityTools.jl documentation",
    modules = [CausalityTools,
                TransferEntropy,
                PerronFrobenius,
                StateSpaceReconstruction,
                TimeseriesSurrogates,
                CrossMappings],
    format = :markdown,
    pages = PAGES
)

if !Sys.iswindows()
    deploydocs(
        deps   = Deps.pip("mkdocs==0.17.5", "mkdocs-material==2.9.4",
        "python-markdown-math", "pygments", "pymdown-extensions"),
        repo   = "github.com/kahaaga/CausalityTools.jl.git",
        target = "site",
        make = () -> run(`mkdocs build`)
    )
end


text:
struct TimeDomainData{Q,V}
    ts::AbstractArray{Q,1}
    vals::Vector{V}
    function TimeDomainData(t::AbstractArray{Q,1}, v::Vector{V}) where {Q<:Real,V<:Real}
        if length(t) !== length(v)
            error("Number of time points much match number of sampled values.")
        else
            new{Q,V}(t,v)
        end
    end
end

struct FrequencyDomainData{R,C}
    fs::AbstractArray{R,1}
    vals::Vector{C}
    function
        FrequencyDomainData(f::AbstractArray{R,1}, v::Vector{C}) where {R<:Real, C<:Complex}
        if length(f) !== length(v)
            error("Number of time points much match number of sampled values.")
        else
            new{R,C}(f,v)
        end
    end
end


text:
import REPL: ends_with_semicolon

Base.push!(x::Set{Cell}) = x

"Run given cells and all the cells that depend on them."
function run_reactive!(notebook::Notebook, cells::Array{Cell,1}; deletion_hook::Function=WorkspaceManager.delete_vars)::CellTopology
	# make sure that we're the only `run_reactive!` being executed - like a semaphor
	token = take!(notebook.executetoken)

	# save the old topology - we'll delete variables assigned from it and re-evalutate its cells
	old_topology = topological_order(notebook, cells)

	old_runnable = old_topology.runnable
	to_delete_vars = union!(Set{Symbol}(), (runnable.symstate.assignments for runnable in old_runnable)...)
	to_delete_funcs = union!(Set{Vector{Symbol}}(), (keys(runnable.symstate.funcdefs) for runnable in old_runnable)...)

	# update the cache using the new code and compute the new topology
	update_caches!(notebook, cells)

	new_topology = topological_order(notebook, union(cells, keys(old_topology.errable)))
	to_run = setdiff(union(new_topology.runnable, old_topology.runnable), keys(new_topology.errable))::Array{Cell,1} # TODO: think if old error cell order matters

	# change the bar on the sides of cells to "running"
	for cell in to_run
		cell.running = true
		putnotebookupdates!(notebook, clientupdate_cell_running(notebook, cell))
	end
	for (cell, error) in new_topology.errable
		cell.running = false
		relay_reactivity_error!(cell, error)
		putnotebookupdates!(notebook, clientupdate_cell_output(notebook, cell))
	end

	# delete new variables that will be defined by a cell
	new_runnable = new_topology.runnable
	to_delete_vars = union!(to_delete_vars, (runnable.symstate.assignments for runnable in new_runnable)...)
	to_delete_funcs = union!(to_delete_funcs, (keys(runnable.symstate.funcdefs) for runnable in new_runnable)...)

	# delete new variables in case a cell errors (then the later cells show an UndefVarError)
	new_errable = keys(new_topology.errable)
	to_delete_vars = union!(to_delete_vars, (errable.symstate.assignments for errable in new_errable)...)
	to_delete_funcs = union!(to_delete_funcs, (keys(errable.symstate.funcdefs) for errable in new_errable)...)

	to_reimport = union(Set{Expr}(), map(c->c.module_usings, setdiff(notebook.cells, to_run))...)
	deletion_hook(notebook, to_delete_vars, to_delete_funcs, to_reimport; to_run=to_run) # `deletion_hook` defaults to `WorkspaceManager.delete_vars`

	local any_interrupted = false
	for (i, cell) in enumerate(to_run)
		if any_interrupted
			relay_reactivity_error!(cell, InterruptException())
		else
			run = run_single!(notebook, cell)
			any_interrupted |= run.interrupted
		end
		cell.running = false
		putnotebookupdates!(notebook, clientupdate_cell_output(notebook, cell))
	end

	# allow other `run_reactive!` calls to be executed
	put!(notebook.executetoken, token)
	return new_topology
end

"See `run_reactive`."
function run_reactive_async!(notebook::Notebook, cells::Array{Cell,1}; kwargs...)::Task
	@async begin
		# because this is being run asynchronously, we need to catch exceptions manually
		try
			run_reactive!(notebook, cells; kwargs...)
		catch ex
			bt = stacktrace(catch_backtrace())
			showerror(stderr, ex, bt)
		end
	end
end

run_reactive!(notebook::Notebook, cell::Cell; kwargs...) = run_reactive!(notebook, [cell]; kwargs...)::CellTopology
run_reactive_async!(notebook::Notebook, cell::Cell; kwargs...) = run_reactive_async!(notebook, [cell]; kwargs...)::Task

"Run a single cell non-reactively, return run information."
function run_single!(notebook::Notebook, cell::Cell)
	run = WorkspaceManager.eval_fetch_in_workspace(notebook, cell.parsedcode, cell.cell_id, ends_with_semicolon(cell.code))
	cell.runtime = run.runtime

	cell.output_repr = run.output_formatted[1]
	cell.repr_mime = run.output_formatted[2]
	cell.errored = run.errored

	return run
	# TODO: capture stdout and display it somehwere, but let's keep using the actual terminal for now
end

"Parse and analyze code for the given `cells`; analyze (indirect) function calls."
function update_caches!(notebook, cells)
	for cell in cells
		start_cache!(notebook, cell)
	end
	update_funcdefs!(notebook)
	for cell in cells
		finish_cache!(notebook, cell)
	end
end

"Update a single cell's cache - parsed code etc"
function start_cache!(notebook::Notebook, cell::Cell)
	cell.parsedcode = parse_custom(notebook, cell)
	cell.module_usings = ExploreExpression.compute_usings(cell.parsedcode)
	cell.rootassignee = ExploreExpression.get_rootassignee(cell.parsedcode)
	cell.symstate = try
		ExploreExpression.compute_symbolreferences(cell.parsedcode)
	catch ex
		@error "Expression explorer failed on: " cell.code
		showerror(stderr, ex, stacktrace(backtrace()))
		ExploreExpression.SymbolsState()
	end
end

"Account for globals referenced in function calls by including `SymbolsState`s from called functions in the cell itself."
function finish_cache!(notebook::Notebook, cell::Cell)
	calls = all_indirect_calls(notebook, cell.symstate)
	calls = union!(calls, keys(cell.symstate.funcdefs)) # _assume_ that all defined functions are called inside the cell to trigger eager reactivity.
	filter!(in(keys(notebook.combined_funcdefs)), calls)

	union!(cell.symstate.references, (notebook.combined_funcdefs[func].references for func in calls)...)
	union!(cell.symstate.assignments, (notebook.combined_funcdefs[func].assignments for func in calls)...)

	add_funcnames!(notebook, cell, calls)
end

"""Add method calls and definitions as symbol references and definition, resp.

Will add `Module.func` (stored as `Symbol[:Module, :func]`) as Symbol("Module.func") (which is not the same as the expression `:(Module.func)`)."""
function add_funcnames!(notebook::Notebook, cell::Cell, calls::Set{Vector{Symbol}})
	push!(cell.symstate.references, (cell.symstate.funccalls .|> join_funcname_parts)...)
	push!(cell.symstate.assignments, (keys(cell.symstate.funcdefs) .|> join_funcname_parts)...)

	union!(cell.symstate.references, (notebook.combined_funcdefs[func].funccalls .|> join_funcname_parts for func in calls)...)
	union!(cell.symstate.assignments, (keys(notebook.combined_funcdefs[func].funcdefs) .|> join_funcname_parts for func in calls)...)
end