text:
#=
In this file the core computations for creating a recurrence matrix
are defined (via multiple dispatch).

The low level interface is contained in the function
`recurrence_matrix`, and this is where any specialization should happen.
=#
################################################################################
# AbstractRecurrenceMatrix type hierarchy and extensions of methods
################################################################################
const FAN = NeighborNumber

abstract type AbstractRecurrenceMatrix{SearchType} end
const ARM = AbstractRecurrenceMatrix

struct RecurrenceMatrix{SearchType} <: AbstractRecurrenceMatrix{SearchType}
    data::SparseMatrixCSC{Bool,Int}
end
struct CrossRecurrenceMatrix{SearchType} <: AbstractRecurrenceMatrix{SearchType}
    data::SparseMatrixCSC{Bool,Int}
end
struct JointRecurrenceMatrix{SearchType} <: AbstractRecurrenceMatrix{SearchType}
    data::SparseMatrixCSC{Bool,Int}
end

function Base.summary(R::AbstractRecurrenceMatrix)
    N = nnz(R.data)
    return "$(nameof(typeof(R))) of size $(size(R.data)) with $N entries"
end
Base.show(io::IO, R::AbstractRecurrenceMatrix) = println(io, summary(R))

# Propagate used functions:
begin
    extentions = [
        (:Base, (:getindex, :size, :length, :view, :iterate,
            :eachindex, :axes, :CartesianIndices)),
        (:LinearAlgebra, (:diag, :triu, :tril, :issymmetric)),
        (:SparseArrays, (:nnz, :rowvals, :nzrange, :nonzeros))
    ]
    for (M, fs) in extentions
        for f in fs
            @eval $M.$(f)(x::ARM, args...) = $(f)(x.data, args...)
        end
    end
end

for operator in [:(==), :(!=)]
    @eval Base.$operator(x::ARM, y::ARM) = $operator(x.data, y.data)
end

LinearAlgebra.issymmetric(::RecurrenceMatrix{WithinRange}) = true
LinearAlgebra.issymmetric(::JointRecurrenceMatrix{WithinRange}) = true
LinearAlgebra.issymmetric(::ARM) = false

# column values in sparse matrix (parallel to rowvals)
function colvals(x::SparseMatrixCSC)
    cv = zeros(Int, nnz(x))
    @inbounds for c in axes(x,2)
        cv[nzrange(x,c)] .= c
    end
    cv
end
colvals(x::ARM) = colvals(x.data)

# Convert to matrices
Base.Array(R::ARM) = Matrix(R.data)
Base.Matrix(R::ARM) = Matrix(R.data)
SparseArrays.SparseMatrixCSC(R::ARM) = SparseMatrixCSC(R.data)

################################################################################
# Concrete Implementations & Documentation
################################################################################
"""
    RecurrenceMatrix(x, ε::Real; kwargs...)

Create a recurrence matrix from trajectory `x` (either a `Dataset` or a `Vector`)
and with distance threshold `ε`.
Objects of type `<:AbstractRecurrenceMatrix` are displayed as a [`recurrenceplot`](@ref).

See also: [`CrossRecurrenceMatrix`](@ref), [`JointRecurrenceMatrix`](@ref) and
use [`recurrenceplot`](@ref) to turn the result of these functions into a plottable format.

## Keyword Arguments
* `metric = "euclidean"` : metric of the distances, either `Metric` or a string,
   as in [`distancematrix`](@ref).
* `scale = 1` : a function of the distance matrix (see [`distancematrix`](@ref)),
  or a fixed number, used to scale the value of `ε`. Typical choices are
  `maximum` or `mean`, such that the threshold `ε` is defined as a ratio of the
  maximum or the mean distance between data points, respectively (using
  `mean` or `maximum` calls specialized versions that are faster than the naive
  approach).  Use `1` to keep the distances unscaled (default).
* `fixedrate::Bool = false` : a flag that indicates if `ε` should be
  taken as a target fixed recurrence rate (see [`recurrencerate`](@ref)).
  If `fixedrate` is set to `true`, `ε` must be a value between 0 and 1,
  and `scale` is ignored.
* `parallel::Bool` : whether to parallelize the computation of the recurrence
   matrix using threads. Default depends on available threads and length of `x`.

## Description

The recurrence matrix is a numeric representation of a recurrence plot,
described in detail in [^Marwan2007] and [^Marwan2015]. It represents a square
a sparse square matrix of Boolean values that quantifies recurrences in the trajectory,
i.e., points where the trajectory returns close to itself.
If `d(x[i], x[j]) ≤ ε` (with `d` the distance function),
then the entry `(i, j)` of the matrix will have a `true`
value. The criteria to evaluate distances between data points are defined
based on the keyword arguments.

[^Marwan2007]:
    N. Marwan *et al.*, "Recurrence plots for the analysis of complex systems",
    [Phys. Reports 438*(5-6), 237-329 (2007)](https://doi.org/10.1016/j.physrep.2006.11.001)

[^Marwan2015]:
    N. Marwan & C.L. Webber, *Recurrence Quantification Analysis. Theory and Best Practices*
    [Springer (2015)](https://link.springer.com/book/10.1007/978-3-319-07155-8)
"""
function RecurrenceMatrix{WithinRange}(x, ε; metric = DEFAULT_METRIC,
    parallel::Bool = length(x) > 500 && Threads.nthreads() > 1,
    kwargs...)
    m = getmetric(metric)
    s = resolve_scale(x, m, ε; kwargs...)
    m = recurrence_matrix(x, m, s, Val(parallel))
    return RecurrenceMatrix{WithinRange}(m)
end

RecurrenceMatrix(args...; kwargs...) = RecurrenceMatrix{WithinRange}(args...; kwargs...)

"""
    RecurrenceMatrix{FAN}(x, r::Real; metric = Euclidean(), parallel::Bool)

The parametrized constructor `RecurrenceMatrix{FAN}` creates the recurrence matrix
with a fixed number of neighbors for each point in the phase space, i.e. the number
of recurrences is the same for all columns of the recurrence matrix.
In such case, `r` is taken as the target fixed local recurrence rate,
defined as a value between 0 and 1 which means that `k = round(Int, r*N)` recurrences
for each point are identified with `N = length(x)`.
This is often referred to in the literature as the method of
"Fixed Amount of Nearest Neighbors" (or FAN for short).

Note that while recurrence matrices
with neighbors defined within a given range are always symmetric, those defined
by a fixed amount of neighbors can be non-symmetric.
"""
function RecurrenceMatrix{NeighborNumber}(x, ε; metric = DEFAULT_METRIC,
        parallel::Bool = length(x) > 500 && Threads.nthreads() > 1
    )
    m = getmetric(metric)
    s = get_fan_threshold(x, m, ε)
    m = recurrence_matrix(x, x, m, s, Val(parallel)) # to allow for non-symmetry
    return RecurrenceMatrix{NeighborNumber}(m)
end

"""
    CrossRecurrenceMatrix(x, y, ε::Real; kwargs...)
    CrossRecurrenceMatrix{FAN}(x, y, k::Int; kwargs...)

Create a cross recurrence matrix from trajectories `x` and `y`.

The cross recurrence matrix is a bivariate extension of the recurrence matrix.
For the time series `x`, `y`, of length `n` and `m`, respectively, it is a
sparse `n×m` matrix of Boolean values, such that if `d(x[i], y[j]) ≤ ε`,
then the cell `(i, j)` of the matrix will have a `true` value.

Note that, unlike univariate recurrence matrices, cross recurrence matrices
are not generally symmetric, regardless of the method used to make them.

See [`RecurrenceMatrix`](@ref) for details, references and keywords.
See also: [`JointRecurrenceMatrix`](@ref).
"""
function CrossRecurrenceMatrix{WithinRange}(x, y, ε; metric = DEFAULT_METRIC,
    parallel::Bool = length(x) > 500 && Threads.nthreads() > 1,
    kwargs...)
    m = getmetric(metric)
    s = resolve_scale(x, y, m, ε; kwargs...)
    m = recurrence_matrix(x, y, m, s, Val(parallel))
    return CrossRecurrenceMatrix{WithinRange}(m)
end

CrossRecurrenceMatrix(args...; kwargs...) = CrossRecurrenceMatrix{WithinRange}(args...; kwargs...)

function CrossRecurrenceMatrix{NeighborNumber}(x, y, ε; metric = DEFAULT_METRIC,
    parallel::Bool = length(x) > 500 && Threads.nthreads() > 1,
    kwargs...)
    m = getmetric(metric)
    s = get_fan_threshold(x, y, m, ε)
    m = recurrence_matrix(x, y, m, s, Val(parallel))
    return CrossRecurrenceMatrix{NeighborNumber}(m)
end

"""
    JointRecurrenceMatrix(x, y, ε; kwargs...)
    JointRecurrenceMatrix{FAN}(x, y, ε; kwargs...)

Create a joint recurrence matrix from `x` and `y`.

The joint recurrence matrix considers the recurrences of the trajectories
of `x` and `y` separately, and looks for points where both recur
simultaneously. It is calculated by the element-wise multiplication
of the recurrence matrices of `x` and `y`. If `x` and `y` are of different
length, the recurrences are only calculated until the length of the shortest one.

See [`RecurrenceMatrix`](@ref) for details, references and keywords.
See also: [`CrossRecurrenceMatrix`](@ref).
"""
function JointRecurrenceMatrix{SearchType}(x, y, ε; kwargs...) where SearchType
    n = min(size(x,1), size(y,1))
    if n == size(x,1) && n == size(y,1)
        rm1 = RecurrenceMatrix{SearchType}(x, ε; kwargs...)
        rm2 = RecurrenceMatrix{SearchType}(y, ε; kwargs...)
    else
        rm1 = RecurrenceMatrix{SearchType}(x[1:n,:], ε; kwargs...)
        rm2 = RecurrenceMatrix{SearchType}(y[1:n,:], ε; kwargs...)
    end
    return JointRecurrenceMatrix{SearchType}(rm1.data .* rm2.data)
end

JointRecurrenceMatrix(args...; kwargs...) = JointRecurrenceMatrix{WithinRange}(args...; kwargs...)

"""
    JointRecurrenceMatrix(R1, R2; kwargs...)

Create a joint recurrence matrix from given recurrence matrices `R1, R2`.
"""
function JointRecurrenceMatrix(
    R1::AbstractRecurrenceMatrix{SearchType}, R2::AbstractRecurrenceMatrix{SearchType}; kwargs...
    ) where SearchType
    R3 = R1.data .* R2.data
    return JointRecurrenceMatrix{SearchType}(R3)
end

################################################################################
# Scaling / fixed rate / fixed amount of neighbors
################################################################################
# here args... is (x, y, metric, ε) or just (x, metric, ε)
function resolve_scale(x, y, metric, ε; scale=1, fixedrate=false)
    if fixedrate
        @assert 0 < ε < 1 "ε must be within (0, 1) for fixed rate"
        sfun = (m) -> quantile(vec(m), ε)
        return resolve_scale(x, y, metric, 1.0; scale=sfun, fixedrate=false)
    else
        scale_value = _computescale(scale, x, y, metric)
        return ε*scale_value
    end
end
resolve_scale(x, metric, ε) = resolve_scale(x, x, metric, ε)

# If `scale` is a function, compute the numeric value of the scale based on the
# distance matrix; otherwise return the value of `scale` itself
_computescale(scale::Real, args...) = scale

function _computescale(scale::Function, x, y, metric)
    if x===y
        distances = zeros(Int(length(x)*(length(x)-1)/2), 1)
        c = 0
        @inbounds for i in 1:length(x)-1, j=(i+1):length(x)
            distances[c+=1] = evaluate(metric, x[i], y[j])
        end
    else
        distances = distancematrix(x, y, metric)
    end
    return scale(distances)
end
# specific methods to avoid `distancematrix`
function _computescale(scale::typeof(maximum), x, y, metric::Metric)
    maxvalue = zero(eltype(x))
    if x===y
        @inbounds for i in 1:length(x)-1, j=(i+1):length(x)
            newvalue = evaluate(metric, x[i], y[j])
            (newvalue > maxvalue) && (maxvalue = newvalue)
        end
    else
        @inbounds for xi in x, yj in y
            newvalue = evaluate(metric, xi, yj)
            (newvalue > maxvalue) && (maxvalue = newvalue)
        end
    end
    return maxvalue
end
function _computescale(scale::typeof(mean), x, y, metric::Metric)
    meanvalue = 0.0
    if x===y
        @inbounds for i in 1:length(x)-1, j=(i+1):length(x)
            meanvalue += evaluate(metric, x[i], y[j])
        end
        denominator = length(x)*(length(x)-1)/2
    else
        @inbounds for xi in x, yj in y
            meanvalue += evaluate(metric, xi, yj)
        end
        denominator = Float64(length(x)*length(y))
    end
    return meanvalue/denominator
end

"""
    get_fan_threshold(args) → ε_fan
Compute the adaptive Fixed Amount of Neibours (FAN) `ε_fan`
Here args... is (x, y, metric, ε) or just (x, metric, ε)
"""
function get_fan_threshold(x, y, metric, ε)
    @assert 0 < ε < 1 "Global recurrence rate must be ∈ (0, 1)"
    fan_threshold = zeros(length(y))
    d = distancematrix(x, y, metric)
    if x === y
        ε += 1/length(x)
    end
    for i in axes(d, 2)
        fan_threshold[i] = quantile(view(d, : ,i), ε)
    end
    return fan_threshold
end

get_fan_threshold(x, metric, ε) = get_fan_threshold(x, x, metric, ε)

################################################################################
# recurrence_matrix - Low level function
################################################################################
# TODO: increase the efficiency here by not computing everything!

# Core function

# First, we define the non-parallel versions.

# For two datasets

"""
    recurrence_matrix(x, y, metric::Metric, ε, parallel::Val)

Return a sparse matrix which encodes recurrence points.

Note that `parallel` may be either `Val(true)` or `Val(false)`.
"""
function recurrence_matrix(xx::AbstractDataset, yy::AbstractDataset, metric::Metric, ε, ::Val{false})
    x = xx.data
    y = yy.data
    @assert ε isa Real || length(ε) == length(y)
    rowvals = Vector{Int}()
    colvals = Vector{Int}()
    for j in eachindex(y)
        nzcol = 0
        for i in eachindex(x)
            @inbounds if evaluate(metric, x[i], y[j]) ≤ ( (ε isa Real) ? ε : ε[j] )
                push!(rowvals, i)
                nzcol += 1
            end
        end
        append!(colvals, fill(j, (nzcol,)))
    end
    nzvals = fill(true, (length(rowvals),))
    return sparse(rowvals, colvals, nzvals, length(x), length(y))
end

# Vector version can be more specialized (and metric is irrelevant)
function recurrence_matrix(x::AbstractVector, y::AbstractVector, metric::Metric, ε, ::Val{false})
    @assert ε isa Real || length(ε) == length(y)
    rowvals = Vector{Int}()
    colvals = Vector{Int}()
    for j in eachindex(y)
        nzcol = 0
        for i in eachindex(x)
            if @inbounds abs(x[i] - y[j]) ≤ ( (ε isa Real) ? ε : ε[j] )
                push!(rowvals, i)
                nzcol += 1
            end
        end
        append!(colvals, fill(j, (nzcol,)))
    end
    nzvals = fill(true, (length(rowvals),))
    return sparse(rowvals, colvals, nzvals, length(x), length(y))
end

# For one dataset
function recurrence_matrix(x::AbstractVector, metric::Metric, ε, ::Val{false})
    @assert ε isa Real || length(ε) == length(y)
    rowvals = Vector{Int}()
    colvals = Vector{Int}()
    for j in eachindex(x)
        nzcol = 0
        for i in 1:j
            if @inbounds abs(x[i] - x[j]) ≤ ( (ε isa Real) ? ε : ε[j] )
                push!(rowvals, i)
                nzcol += 1
            end
        end
        append!(colvals, fill(j, (nzcol,)))
    end
    nzvals = fill(true, (length(rowvals),))
    return Symmetric(sparse(rowvals, colvals, nzvals, length(x), length(x)), :U)
end

function recurrence_matrix(xx::AbstractDataset, metric::Metric, ε, ::Val{false})
    x = xx.data
    @assert ε isa Real || length(ε) == length(y)
    rowvals = Vector{Int}()
    colvals = Vector{Int}()
    for j in eachindex(x)
        nzcol = 0
        for i in 1:j
            @inbounds if evaluate(metric, x[i], x[j]) ≤ ( (ε isa Real) ? ε : ε[j] )
                push!(rowvals, i)
                nzcol += 1
            end
        end
        append!(colvals, fill(j, (nzcol,)))
    end
    nzvals = fill(true, (length(rowvals),))
    return Symmetric(sparse(rowvals, colvals, nzvals, length(x), length(x)), :U)
end

# Now, we define the parallel versions of these functions.

# Core function
function recurrence_matrix(xx::AbstractDataset, yy::AbstractDataset, metric::Metric, ε, ::Val{true})
    x = xx.data
    y = yy.data
    @assert ε isa Real || length(ε) == length(y)
    # We create an `Array` of `Array`s, for each thread to have its
    # own array to push to.  This avoids race conditions with
    # multiple threads pushing to the same `Array` (`Array`s are not atomic).
    rowvals = [Vector{Int}() for _ in 1:Threads.nthreads()]
    colvals = [Vector{Int}() for _ in 1:Threads.nthreads()]

    # This is the same logic as the serial function, but parallelized.
    Threads.@threads for j in eachindex(y)
        threadn = Threads.threadid()
        nzcol = 0
        for i in eachindex(x)
            @inbounds if evaluate(metric, x[i], y[j]) ≤ ( (ε isa Real) ? ε : ε[j] )
                push!(rowvals[threadn], i) # push to the thread-specific row array
                nzcol += 1
            end
        end
        append!(colvals[threadn], fill(j, (nzcol,)))
    end
    finalrows = vcat(rowvals...) # merge into one array
    finalcols = vcat(colvals...) # merge into one array
    nzvals = fill(true, (length(finalrows),))
    return sparse(finalrows, finalcols, nzvals, length(x), length(y))
end

# Vector version can be more specialized (and metric is irrelevant)
function recurrence_matrix(x::AbstractVector, y::AbstractVector, metric::Metric, ε, ::Val{true})
    ε isa Real || length(ε) == length(x)
    # We create an `Array` of `Array`s, for each thread to have its
    # own array to push to.  This avoids race conditions with
    # multiple threads pushing to the same `Array` (`Array`s are not atomic).
    rowvals = [Vector{Int}() for _ in 1:Threads.nthreads()]
    colvals = [Vector{Int}() for _ in 1:Threads.nthreads()]

    # This is the same logic as the serial function, but parallelized.
    Threads.@threads for j in eachindex(y)
        threadn = Threads.threadid()
        nzcol = 0
        for i in eachindex(x)
            @inbounds if abs(x[i] - y[j]) ≤ ( (ε isa Real) ? ε : ε[j] )
                push!(rowvals[threadn], i) # push to the thread-specific row array
                nzcol += 1
            end
        end
        append!(colvals[threadn], fill(j, (nzcol,)))
    end
    finalrows = vcat(rowvals...) # merge into one array
    finalcols = vcat(colvals...) # merge into one array
    nzvals = fill(true, (length(finalrows),))
    return sparse(finalrows, finalcols, nzvals, length(x), length(y))
end

function recurrence_matrix(x::AbstractVector, metric::Metric, ε, ::Val{true})
    @assert ε isa Real || length(ε) == length(y)
    # We create an `Array` of `Array`s, for each thread to have its
    # own array to push to.  This avoids race conditions with
    # multiple threads pushing to the same `Array` (`Array`s are not atomic).
    rowvals = [Vector{Int}() for _ in 1:Threads.nthreads()]
    colvals = [Vector{Int}() for _ in 1:Threads.nthreads()]

    # This is the same logic as the serial function, but parallelized.
    Threads.@threads for k in partition_indices(length(x))
        threadn = Threads.threadid()
        for j in k
            nzcol = 0
            for i in 1:j
                @inbounds if abs(x[i] - x[j]) ≤ ( (ε isa Real) ? ε : ε[j] )
                    push!(rowvals[threadn], i) # push to the thread-specific row array
                    nzcol += 1
                end
            end
            append!(colvals[threadn], fill(j, (nzcol,)))
        end
    end
    finalrows = vcat(rowvals...) # merge into one array
    finalcols = vcat(colvals...) # merge into one array
    nzvals = fill(true, (length(finalrows),))
    return Symmetric(sparse(finalrows, finalcols, nzvals, length(x), length(x)), :U)
end

function recurrence_matrix(xx::AbstractDataset, metric::Metric, ε, ::Val{true})
    x = xx.data
    @assert ε isa Real || length(ε) == length(y)
    # We create an `Array` of `Array`s, for each thread to have its
    # own array to push to.  This avoids race conditions with
    # multiple threads pushing to the same `Array` (`Array`s are not atomic).
    rowvals = [Vector{Int}() for _ in 1:Threads.nthreads()]
    colvals = [Vector{Int}() for _ in 1:Threads.nthreads()]

    # This is the same logic as the serial function, but parallelized.
    Threads.@threads for k in partition_indices(length(x))
        threadn = Threads.threadid()
        for j in k
            nzcol = 0
            for i in 1:j
                @inbounds if evaluate(metric, x[i], x[j]) ≤ ( (ε isa Real) ? ε : ε[j] )
                    push!(rowvals[threadn], i) # push to the thread-specific row array
                    nzcol += 1
                end
            end
            append!(colvals[threadn], fill(j, (nzcol,)))
        end
    end
    finalrows = vcat(rowvals...) # merge into one array
    finalcols = vcat(colvals...) # merge into one array
    nzvals = fill(true, (length(finalrows),))
    return Symmetric(sparse(finalrows, finalcols, nzvals, length(x), length(x)), :U)
end


text:
export
    xf2eul,
    xf2rav,
    xfmsta,
    xpose6,
    xpose,
    xposeg

"""
    xf2eul(xform, axisa, axisb, axisc)

Convert a state transformation matrix to Euler angles and their derivatives
with respect to a specified set of axes.

### Arguments ###

- `xform`: A state transformation matrix
- `axisa`: Axis A of the Euler angle factorization
- `axisb`: Axis B of the Euler angle factorization
- `axisc`: Axis C of the Euler angle factorization

### Output ###

Returns a tuple of an array of Euler angles and their derivatives and a boolean
that indicates whether these are a unique representation.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/xf2eul_c.html)
"""
function xf2eul(xform, axisa, axisb, axisc)
    @checkdims 6 6 xform
    eulang = Array{SpiceDouble}(undef, 6)
    unique = Ref{SpiceBoolean}()
    ccall((:xf2eul_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, SpiceInt, SpiceInt, SpiceInt, Ref{SpiceDouble}, Ref{SpiceBoolean}),
          xform, axisa, axisb, axisc, eulang, unique)
    handleerror()
    eulang, Bool(unique[])
end

"""
    xf2rav(xform)

Determines the rotation matrix and angular velocity of the rotation from a
state transformation matrix.

### Arguments ###

- `xform`: State transformation matrix

### Output ###

Returns a tuple of the rotation matrix and the angular velocity vector
associated with `xform`.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/xf2rav_c.html)
"""
function xf2rav(xform)
    @checkdims 6 6 xform
    rot = Array{SpiceDouble}(undef, 3, 3)
    av = Array{SpiceDouble}(undef, 3)
    ccall((:xf2rav_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Ref{SpiceDouble}, Ref{SpiceDouble}),
          permutedims(xform), rot, av)
    permutedims(rot), av
end

"""
    xfmsta(input_state, input_coord_sys, output_coord_sys, body)

Transform a state between coordinate systems.

### Arguments ###

- `input_state`: Input state
- `input_coord_sys`: Current (input) coordinate system
- `output_coord_sys: Desired (output) coordinate system
- `body`: Name or NAIF ID of body with which coordinates are associated (if applicable)

### Output ###

Returns the converted output state.

### References ###

- [NAIF Documentation](https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/xfmsta_c.html)
"""
function xfmsta(input_state, input_coord_sys, output_coord_sys, body)
    @checkdims 6 input_state
    output_state = Array{SpiceDouble}(undef, 6)
    ccall((:xfmsta_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, Cstring, Cstring, Cstring, Ref{SpiceDouble}),
          input_state, input_coord_sys, output_coord_sys, body, output_state)
    handleerror()
    output_state
end

function _xpose6(m)
    mout = Array{SpiceDouble}(undef, 6, 6)
    ccall((:xpose6_c, libcspice), Cvoid, (Ref{SpiceDouble}, Ref{SpiceDouble}), m, mout)
    mout
end

"""
    xpose6(matrix)

!!! warning "Deprecated"
    Use `transpose(matrix)` instead.
"""
xpose6

@deprecate xpose6 transpose

function _xpose(m)
    mout = Array{SpiceDouble}(undef, 3, 3)
    ccall((:xpose_c, libcspice), Cvoid, (Ref{SpiceDouble}, Ref{SpiceDouble}), m, mout)
    mout
end

"""
    xpose(matrix)

!!! warning "Deprecated"
    Use `transpose(matrix)` instead.
"""
xpose

@deprecate xpose transpose

function _xposeg(matrix)
    m, n = size(matrix)
    mout = Array{SpiceDouble}(undef, n, m)
    ccall((:xposeg_c, libcspice), Cvoid,
          (Ref{SpiceDouble}, SpiceInt, SpiceInt, Ref{SpiceDouble}), matrix, n, m, mout)
    mout
end

"""
    xposeg(matrix)

!!! warning "Deprecated"
    Use `transpose(matrix)` instead.
"""
xposeg

@deprecate xposeg transpose

ext: jl


text:
using ArtifactUtils, Artifacts
add_artifact!(
    "../Artifacts.toml",
    "refractiveindex.info",
    "https://github.com/polyanskiy/refractiveindex.info-database/archive/v2020-01-19.tar.gz",
    # lazy=true, # need to use LazyArtifacts for this
    force=true,
    clear=false,
)


text:
"""
    $TYPEDEF

Generic state machine functor type.

Example
```julia
bar!(usrdata) = IncrementalInference.exitStateMachine
foo!(usrdata) = bar!

sm = StateMachine(next=foo!)
usrdata = nothing
while st(usrdata); end
```

Notes
- Also see IncrementalInference/test/testStateMachine.jl
"""
mutable struct StateMachine{T}
  next::Function
  iter::Int
  history::Vector{Tuple{DateTime, Int, Function, T}}
  name::String
  StateMachine{T}(;next=emptyState, iter::Int=0, name::AbstractString="") where T = new{T}(next, iter, Vector{Tuple{DateTime, Int, Function, T}}(), name)
end

"""
    $SIGNATURES

Run state machine function (as functor).

Notes
- `timeout::Union{Real,Nothing}` is optional with default `=nothing`.
  - this code is skipped in lowered and llvm code if not used
  - subroutine will use `pollinterval::Real` [seconds] to interrogate during `timeout::Real` [seconds] period.
- can stop FSM early by using any of the following:
  - `breakafter`, `iterlimit`.
- Can `injectDelayBefore` a function `st.next` to help with debugging.
- can print FSM steps with `verbose=true`.
  - `verbosefid::IOStream` is used as destination for verbose output, default is `stdout`.
- FSM steps and `userdata` can be recorded in standard `history` format using `recordhistory=true`.
- `housekeeping_cb` is callback to give user access to `StateMachine` internals and opportunity to insert bespoke operations.

Example
```julia
bar!(usrdata) = IncrementalInference.exitStateMachine
foo!(usrdata) = bar!

sm = StateMachine(next=foo!)
usrdata = nothing
while st(usrdata); end
```
"""
function (st::StateMachine{T})( userdata::T=nothing,
                                timeout::Union{Nothing,<:Real}=nothing;
                                pollinterval::Real=0.05,
                                breakafter::Function=exitStateMachine,
                                verbose::Bool=false,
                                verbosefid=stdout,
                                verboseXtra="",
                                iterlimit::Int=-1,
                                injectDelayBefore::Union{Nothing,Pair{<:Function, <:Real}}=nothing,
                                recordhistory::Bool=false,
                                housekeeping_cb::Function=(st)->()  ) where {T}
  #
  st.iter += 1
  # verbose print to help debugging
  !verbose ? nothing : println(verbosefid, "FSM $(st.name), iter=$(st.iter) -- $(st.next), $(verboseXtra)")
  # early exit plumbing
  retval = st.next != breakafter && (iterlimit == -1 || st.iter < iterlimit)
  # record steps for later
  T0 = Dates.now()
  recordhistory ? push!(st.history, (T0, st.iter, deepcopy(st.next), deepcopy(userdata))) : nothing
  # user has some special situation going on.
  housekeeping_cb(st)
  (injectDelayBefore !== nothing && injectDelayBefore[1] == st.next) ? sleep(injectDelayBefore[2]) : nothing
  if timeout === nothing
    # no watchdog, just go and optimize llvm lowered code
    st.next = st.next(userdata)
  else
    # add the watchdog into the llvm lowered code
    currtsk = current_task()
    # small amount of memory usage, but must guarantee InterruptException is not accidently fired during next step.
    doneWatchdog = Base.RefValue{Int}(0)
    wdt = @async begin
      # wait for watchdog timeperiod in a seperate co-routine
      res = timedwait(()->doneWatchdog[]==1, timeout, pollint=pollinterval)
      # Two requirements needed to interrupt FSM step
      res == :timed_out && doneWatchdog[] == 0 ? schedule(currtsk, InterruptException(), error=true) : nothing
    end
    st.next = st.next(userdata)
    doneWatchdog[] = 1
  end
  return retval
end

"""
    $SIGNATURES

Dummy function in case a `statemachine.next` is not initialized properly.
"""
function emptyState(dummy)
  @warn "Empty state machine, assign `next` to entry function -- i.e. StateMachine(next=foo)"
  return exitStateMachine
end

"""
    $SIGNATURES

Default function used for exiting any state machine.
"""
function exitStateMachine(dummy)
  return emptyState
end

"""
    $SIGNATURES

How many iterations has this `::StateMachine` stepped through.
"""
getIterCount(st::StateMachine) = st.iter

"""
  $SIGNATURES

Repeat a state machine step without changing history or primary values.
"""
function sandboxStateMachineStep(hist::Vector{Tuple{DateTime, Int, <:Function, T}},
                                 step::Int  ) where T
  #
  usrdata = deepcopy(hist[step][4])
  @time nextfnc = hist[step][3](usrdata)
  return (hist[step][1], step+1, nextfnc, usrdata)
end

getStateLabel(state) = Symbol(split(split(string(state), '_')[1],'.')[end])

function histStateMachineTransitions(hist::T;
                                     allStates=Vector{Symbol}(),
                                     stateVisits = Dict{Symbol, Vector{Symbol}}() ) where {T <: Array}
  # local memory

  # find all states and transitions and add as outgoing vector of edges from state lbl
  for i in 1:(length(hist)-1)
    sta  = string(hist[i][3])
    lbl = getStateLabel(sta)
    nsta = string(hist[i+1][3])
    nlbl = getStateLabel(nsta)
    if !haskey(stateVisits, lbl)
      stateVisits[lbl] = Symbol[nlbl;]
    else
      push!(stateVisits[lbl], nlbl)
    end
    union!(allStates, [lbl; nlbl])
  end

  return stateVisits, allStates
end

#


text:
## S2MktInt -----------------------------------------------------
function S2MktInt(ds2_mkt_int::Dict{Symbol, Any})
  shock_object = :CapMkt
  shock = ds2_mkt_int[:shock]
  spot_up_abs_min = ds2_mkt_int[:spot_up_abs_min]
  balance = DataFrame()
  scr = zeros(Float64, 2)
  scen_up = false
  return S2MktInt(shock_object,
                  shock,
                  spot_up_abs_min, balance, scr, scen_up)
end

function S2MktInt(param::ProjParam,
                  s2_balance::DataFrame,
                  ds2_mkt_int::Dict{Symbol, Any})
  p = deepcopy(param)
  mkt_int = S2MktInt(ds2_mkt_int)
  mkt_int.balance = deepcopy(s2_balance)
  for 𝑖𝑛𝑡_𝑡𝑦𝑝𝑒_𝑠𝑦𝑚𝑏 ∈ collect(keys(mkt_int.shock))
    append!(mkt_int.balance,
            s2bal(p, mkt_int,
                  (inv, s2_int) ->
                  mktintshock!(inv, s2_int, 𝑖𝑛𝑡_𝑡𝑦𝑝𝑒_𝑠𝑦𝑚𝑏),
                  𝑖𝑛𝑡_𝑡𝑦𝑝𝑒_𝑠𝑦𝑚𝑏))
  end
  scr!(mkt_int)
  return mkt_int
end

## S2MktEq ------------------------------------------------------
function S2MktEq(ds2_mkt_eq::Dict{Symbol, Any})
  shock_object = :CapMkt_AdjVal0
  eq_type = Dict{Symbol, Symbol}()
  balance = DataFrame()
  corr = ds2_mkt_eq[:corr]
  shock = ds2_mkt_eq[:shock]
  scr = zeros(Float64, 2)
  return S2MktEq(shock_object, eq_type, shock,
                 balance, corr, scr)
end

function S2MktEq(param::ProjParam,
                 s2_balance::DataFrame,
                 ds2_mkt_eq::Dict{Symbol, Any},
                 eq2type)
  p = deepcopy(param)
  mkt_eq = S2MktEq(ds2_mkt_eq)
  merge!(mkt_eq.eq2type, eq2type)
  mkt_eq.balance = deepcopy(s2_balance)
  for 𝑠𝑦𝑚𝑏 ∈ collect(keys(mkt_eq.shock))
    append!(mkt_eq.balance,
            s2bal(p, mkt_eq,
                  (invs, s2_eq) ->
                  mkteqshock!(invs, s2_eq, 𝑠𝑦𝑚𝑏),
                  𝑠𝑦𝑚𝑏))
  end
  scr!(mkt_eq)
  return mkt_eq
end

## S2Mkt --------------------------------------------------------
function S2Mkt(ds2_mkt::Dict{Symbol, Any})
  mds = Array{S2Module}(undef, 0)
  corr_up = ds2_mkt[:corr](ds2_mkt[:raw], ds2_mkt[:adj], :up)
  corr_down = ds2_mkt[:corr](ds2_mkt[:raw], ds2_mkt[:adj], :down)
  scr = zeros(Float64, 2)
  return S2Mkt(mds, corr_up, corr_down, scr)
end

function S2Mkt(param::ProjParam,
               s2_balance::DataFrame,
               eq2type::Dict,
               ds2_mkt_all::Dict)
  p = deepcopy(param)
  mkt = S2Mkt(ds2_mkt_all[:mkt])
  push!(mkt.mds, S2MktInt(p, s2_balance,
                 ds2_mkt_all[:mkt_int]))
  push!(mkt.mds, S2MktEq(p, s2_balance,
                 ds2_mkt_all[:mkt_eq], eq2type))
  push!(mkt.mds, S2MktProp(p, s2_balance))
  push!(mkt.mds, S2MktSpread(p, s2_balance))
  push!(mkt.mds, S2MktFx(p, s2_balance))
  push!(mkt.mds, S2MktConc(p, s2_balance))
  scen_up = false
  for 𝑖 ∈ 1:length(mkt.mds)
    if :scen_up ∈ fieldnames(typeof(mkt.mds[𝑖]))
      scen_up = mkt.mds[𝑖].scen_up
    end
  end
  corr = (scen_up ? mkt.corr_up : mkt.corr_down)
  mkt.scr = scr(mkt, corr)
  return mkt
end

## S2Def1 -------------------------------------------------------
function S2Def1(ds2_def_1)
  tlgd = Array{Float64}(undef, 0)
  slgd = Array{Float64}(undef, 0)
  u = Array{Float64}(undef, 0,0)
  v = Array{Float64}(undef, 0)
  scr_par = Dict{Symbol, Vector{Float64}}()
  for 𝑖 ∈ 1:nrow(ds2_def_1[:scr_par])
    merge!(scr_par,
           Dict(ds2_def_1[:scr_par][𝑖, :range] =>
            [ds2_def_1[:scr_par][𝑖, :threshold_upper],
             ds2_def_1[:scr_par][𝑖, :multiplier]]))
  end
  scr = zeros(Float64, 2)
  return S2Def1(tlgd, slgd, u, v, scr_par, scr)
end

function S2Def1(param::ProjParam,
                ds2_def_1::Dict{Symbol, Any})
  p = deepcopy(param)
  def = S2Def1(ds2_def_1)
  cqs_vec = filter(x -> occursin("cqs", string(x)),
                   names(ds2_def_1[:prob]))
  prob = [ds2_def_1[:prob][1, 𝑐𝑞𝑠] for 𝑐𝑞𝑠 ∈ cqs_vec]
  def.tlgd = zeros(Float64, length(cqs_vec))
  def.slgd = zeros(Float64, length(cqs_vec))
  def.u = Array{Float64}(undef, length(cqs_vec), length(cqs_vec))
  def.v = Array{Float64}(undef, length(cqs_vec))

  def.v = 1.5 * prob .* (1 .- prob) ./ (2.5 .- prob)
  for 𝑖 ∈ 1:size(def.u,1), 𝑗 ∈ 1:1:size(def.u,2)
    def.u[𝑖,𝑗] =
      (1-prob[𝑖]) * prob[𝑖] * (1-prob[𝑗]) * prob[𝑗] /
      (1.25 * (prob[𝑖] + prob[𝑗]) - prob[𝑖] * prob[𝑗])
  end
  invs = InvPort(p.t_0, p.T, p.cap_mkt, p.invs_par...)
  for 𝑖 = 1:length(invs.igs[:IGCash].investments)
    j = indexin([invs.igs[:IGCash].investments[𝑖].cqs],
          Symbol.(cqs_vec))[1]
    lgd =
      invs.igs[:IGCash].investments[𝑖].lgd *
      invs.igs[:IGCash].investments[𝑖].mv_0
    def.tlgd[j] += lgd
    def.slgd[j] += lgd * lgd
  end
  scr!(def)
  return def
end

## S2Def --------------------------------------------------------
function S2Def(ds2_def)
  mds = Array{S2Module}(undef, 0)
  corr = ds2_def[:corr]
  scr = zeros(Float64, 2)
  return S2Def(mds, corr, scr)
end

function S2Def(param::ProjParam,
               s2_balance::DataFrame,
               ds2_def_all::Dict)
  p = deepcopy(param)
  def = S2Def(ds2_def_all[:def])
  push!(def.mds, S2Def1(p, ds2_def_all[:def_1]))
  push!(def.mds, S2Def2(p, s2_balance))
  def.scr = scr(def, def.corr)
  return def
end

## S2LifeBio (mortality risk, longevity risk, lapse risk---------
function S2LifeBio(ds2_bio::Dict{Symbol, Any})
  shock_object = :LiabIns
  if :shock_param in collect(keys(ds2_bio))
    shock_param = ds2_bio[:shock_param]
  else
    shock_param = Dict()
  end
  shock = ds2_bio[:shock]
  balance = DataFrame()
  mp_select = Dict{Symbol, Vector{Bool}}()
  scr = zeros(Float64, 2)
  return S2LifeBio(shock_object, shock,
                   shock_param,
                   balance, mp_select, scr)
end

function S2LifeBio(param::ProjParam,
                   s2_balance::DataFrame,
                   ds2_bio::Dict{Symbol, Any})
  p = deepcopy(param)
  bio = S2LifeBio(ds2_bio)
  bio.balance = deepcopy(s2_balance)
  selectmort!(p, bio)
  for 𝑠𝑦𝑚𝑏 ∈ collect(keys(bio.shock))
    append!(bio.balance,
            s2bal(p, bio,
                  (l_ins, wx) -> bioshock!(l_ins, wx, 𝑠𝑦𝑚𝑏),
                  𝑠𝑦𝑚𝑏))
  end
  scr!(bio)
  return bio
end

## S2LifeCost ---------------------------------------------------
function S2LifeCost(ds2_cost::Dict{Symbol, Any})
  shock_object = :InvPort_LiabIns
  shock = ds2_cost[:shock]
  if :shock_param in collect(keys(ds2_cost))
    shock_param = ds2_cost[:shock_param]
  else
    shock_param = Dict()
  end
  balance = DataFrame()
  scr = zeros(Float64, 2)
  return S2LifeCost(shock_object, shock, shock_param,
                    balance, scr)
end

function S2LifeCost(param::ProjParam,
                    s2_balance::DataFrame,
                    ds2_cost::Dict{Symbol, Any})
  p = deepcopy(param)
  cost = S2LifeCost(ds2_cost)
  cost.balance = deepcopy(s2_balance)
  for 𝑠𝑦𝑚𝑏 ∈ collect(keys(cost.shock))
    append!(cost.balance,
            s2bal(p, cost, (invs, l_ins, cst) ->
                  costshock!(invs, l_ins, cst), 𝑠𝑦𝑚𝑏))
  end
  scr!(cost)
  return cost
end

## S2Life -------------------------------------------------------
function S2Life(ds2_life::Dict{Symbol, Any})
  mds = Array{S2Module}(undef, 0)
  corr = ds2_life[:corr]
  scr = zeros(Float64, 2)
  return S2Life(mds, corr, scr)
end

function S2Life(param::ProjParam,
                s2_balance::DataFrame,
                d_life::Dict)
  p = deepcopy(param)
  life = S2Life(d_life[:life])
  push!(life.mds, S2LifeBio(p, s2_balance, d_life[:life_qx]))
  push!(life.mds, S2LifeBio(p, s2_balance, d_life[:life_px]))
  push!(life.mds, S2LifeMorb(p, s2_balance))
  push!(life.mds, S2LifeBio(p, s2_balance, d_life[:life_sx]))
  push!(life.mds, S2LifeCost(p, s2_balance, d_life[:life_cost]))
  push!(life.mds, S2LifeRevision(p, s2_balance))
  push!(life.mds, S2LifeBio(p, s2_balance, d_life[:life_cat]))
  life.scr = scr(life, life.corr)
  return life
end

## S2Op ---------------------------------------------------------
S2Op(ds2_op::Dict{Symbol, Float64}, s2_op::Dict) =
  S2Op(ds2_op,
       s2_op[:prem_earned],
       s2_op[:prem_earned_prev],
       s2_op[:tp], 0, 0, s2_op[:cost_ul], 0.0)

## S2 -----------------------------------------------------------
function   S2(ds2_op, s2_op)
  mds = Array{S2Module}(undef, 0)
  balance = DataFrame()
  corr = zeros(Float64, 5, 5)
  bscr = zeros(Float64, 2)
  adj_tp = 0.0
  adj_dt = 0.0
  op = S2Op(ds2_op, s2_op)
  scr = 0.0
  S2(mds, balance, corr, bscr,
     adj_tp, adj_dt, op, scr, 0, 0, 0, 0, 0)
end

function  S2(param::ProjParam,
             eq2type::Dict{Symbol, Symbol},
             ds2_mkt_all::Dict,
             ds2_def_all::Dict,
             ds2_life_all::Dict,
             s2_op::Dict,
             ds2_op::Dict{Symbol, Float64},
             ds2::Dict{Symbol, Any})
  p = deepcopy(param)
  s2 = S2(ds2_op, s2_op)
  s2.corr = ds2[:corr]
  s2.coc = ds2[:coc]
  s2.balance = s2bal(p)
  s2.op.tp = s2.balance[1, :tpg] + s2.balance[1, :bonus]
  push!(s2.mds, S2Mkt(p, s2.balance, eq2type, ds2_mkt_all))
  push!(s2.mds, S2Def(p, s2.balance, ds2_def_all))
  push!(s2.mds, S2Life(p, s2.balance, ds2_life_all))
  push!(s2.mds, S2Health(p, s2.balance))
  push!(s2.mds, S2NonLife(p, s2.balance))
  scr!(s2, p.tax_credit_0)
  s2.liabs_mod = s2.balance[1,:tpg] +
                 s2.balance[1,:bonus] +
                 s2.balance[1,:cost_prov]
  s2.invest_mod =  s2.balance[1,:invest]
  s2.risk_margin = s2riskmargin(p, s2.scr, s2.coc)
  s2.scr_ratio =
    (s2.invest_mod - s2.liabs_mod - s2.risk_margin)/ s2.scr
  return s2
end


text:
# Combined 2D plots across multiple frames, multi-process version.
#
# To run on a single node,
# julia -p $ncores demo_2d_parallel_pyplot.jl
#
# Hongyang Zhou, hyzhou@umich.edu

using Distributed, ParallelDataTransfer, Glob
@everywhere using Vlasiator, PyPlot, Printf, LaTeXStrings
@everywhere using Vlasiator: set_args, prep2d, set_colorbar, set_plot

@everywhere function init_figure()
   fig, axs = plt.subplots(2, 3; num=myid(),
      figsize=(13.5, 9.5), sharex=true, sharey=true, constrained_layout=true)

   for ax in axs
      ax.set_aspect("equal")
      # Set border line widths
      for loc in ("left", "bottom", "right", "top")
         edge = get(ax.spines, loc, nothing)
         edge.set_linewidth(2.0)
      end
      ax.xaxis.set_tick_params(width=2.0, length=3)
      ax.yaxis.set_tick_params(width=2.0, length=3)
      ax.xaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
      ax.yaxis.set_minor_locator(matplotlib.ticker.AutoMinorLocator())
   end

   axs[2,1].set_xlabel(L"X [$R_E$]"; fontsize)
   axs[2,2].set_xlabel(L"X [$R_E$]"; fontsize)
   axs[2,3].set_xlabel(L"X [$R_E$]"; fontsize)
   axs[1,1].set_ylabel(L"Y [$R_E$]"; fontsize)
   axs[2,1].set_ylabel(L"Y [$R_E$]"; fontsize)

   axs[1,1].set_title("Density"; fontsize)
   axs[1,2].set_title("Velocity x"; fontsize)
   axs[1,3].set_title("Velocity y"; fontsize)
   axs[2,1].set_title("Thermal pressure"; fontsize)
   axs[2,2].set_title("Magnetic field"; fontsize)
   axs[2,3].set_title("Electric field"; fontsize)

   x, y = Vlasiator.get_axis(pArgs1)

   fakedata = zeros(Float32, length(y), length(x))
   c1 = axs[1,1].pcolormesh(x, y, fakedata; norm=cnorm1, cmap=cmap, shading="nearest")
   c2 = axs[1,2].pcolormesh(x, y, fakedata; norm=cnorm2, cmap=cmap, shading="nearest")
   c3 = axs[1,3].pcolormesh(x, y, fakedata;
      norm=cnorm3, cmap=matplotlib.cm.RdBu, shading="nearest")
   c4 = axs[2,1].pcolormesh(x, y, fakedata; norm=cnorm4, cmap=cmap, shading="nearest")
   c5 = axs[2,2].pcolormesh(x, y, fakedata; norm=cnorm5, cmap=cmap, shading="nearest")
   c6 = axs[2,3].pcolormesh(x, y, fakedata; norm=cnorm6, cmap=cmap, shading="nearest")

   format = matplotlib.ticker.FormatStrFormatter("%.1f")
   cb1 = colorbar(c1; ax=axs[1,1], ticks=cticks1, format)
   cb1.ax.set_ylabel("[amu/cc]"; fontsize)
   cb1.outline.set_linewidth(1.0)

   cb2 = colorbar(c2; ax=axs[1,2], ticks=cticks2)
   cb2.ax.set_ylabel("[km/s]"; fontsize)
   cb2.outline.set_linewidth(1.0)

   cb3 = colorbar(c3; ax=axs[1,3], ticks=cticks3)
   cb3.ax.set_ylabel("[km/s]"; fontsize)
   cb3.outline.set_linewidth(1.0)

   cb4 = colorbar(c4; ax=axs[2,1], ticks=cticks4, format)
   cb4.ax.set_ylabel("[nPa]"; fontsize)
   cb4.outline.set_linewidth(1.0)

   cb5 = colorbar(c5; ax=axs[2,2], ticks=cticks5, format)
   cb5.ax.set_ylabel("[nT]"; fontsize)
   cb5.outline.set_linewidth(1.0)

   cb6 = colorbar(c6; ax=axs[2,3], ticks=cticks6)
   cb6.ax.set_ylabel(L"[$\mu V/m$]"; fontsize)
   cb6.outline.set_linewidth(1.0)

   cs = (c1, c2, c3, c4, c5, c6)

   return fig, cs
end

@everywhere function process(fig, cs, file)
   isfile("out/"*file[end-8:end-5]*".png") && return

   println("file = $file")
   meta = load(file)

   data = prep2d(meta, "proton/vg_rho", :mag)'
   cs[1].set_array(data ./ 1e6)

   data = prep2d(meta, "proton/vg_v", :x)'
   cs[2].set_array(data ./ 1e3)

   data = prep2d(meta, "proton/vg_v", :y)'
   cs[3].set_array(data ./ 1e3)

   data = prep2d(meta, "vg_pressure", :mag)'
   cs[4].set_array(data .* 1e9)

   data = prep2d(meta, "vg_b_vol", :z)'
   cs[5].set_array(data .* 1e9)

   data = prep2d(meta, "vg_e_vol", :mag)'
   cs[6].set_array(data .* 1e6)

   str_title = @sprintf "Density pulse run, t= %4.1fs" meta.time
   fig.suptitle(str_title; fontsize="xx-large")

   savefig("out/"*file[end-8:end-5]*".png", bbox_inches="tight")
end

function make_jobs(files)
   for f in files
       put!(jobs, f)
   end
end

@everywhere function do_work(jobs, status)
   fig, cs = init_figure()
   while true
      file = take!(jobs)
      process(fig, cs, file)
      put!(status, true)
   end
   close(fig)
end

############################################################################################
files = glob("bulk*.vlsv", ".")

nfile = length(files)

const jobs   = RemoteChannel(()->Channel{String}(nfile))
const status = RemoteChannel(()->Channel{Bool}(nworkers()))

@passobj 1 workers() files

@broadcast begin # on all workers
   # Set contour plots' axes and colorbars
   const cmap = matplotlib.cm.turbo
   colorscale = Linear
   axisunit = RE

   # Upper/lower limits for each variable
   const ρmin, ρmax   = 0.0, 9.0      # [amu/cc]
   const vxmin, vxmax = -640.0, 100.0 # [km/s]
   const vymin, vymax = -300.0, 300.0 # [km/s]
   const pmin, pmax   = 0.0, 0.7      # [nPa]
   const bmin, bmax   = -5.0, 50.     # [nT]
   const emin, emax   = 5.0, 3e4      # [muV/m]

   meta = load(files[1])

   const pArgs1 = set_args(meta, "proton/vg_rho", axisunit; normal=:none)
   const cnorm1, cticks1 = set_colorbar(colorscale, ρmin, ρmax)
   const cnorm2, cticks2 = set_colorbar(colorscale, vxmin, vxmax)
   const cnorm3, cticks3 = set_colorbar(colorscale, vymin, vymax)
   const cnorm4, cticks4 = set_colorbar(colorscale, pmin, pmax)
   const cnorm5, cticks5 = set_colorbar(Linear, bmin, bmax)
   const cnorm6, cticks6 = set_colorbar(Log, emin, emax)

   const fontsize = "x-large"
end

println("Total number of files: $nfile")
println("Running with $(nworkers()) workers...")

@async make_jobs(files) # Feed the jobs channel with all files to process.

@sync for p in workers()
   @async remote_do(do_work, p, jobs, status)
end

let n = nfile
   t = @elapsed while n > 0 # wait for all jobs to complete
      take!(status)
      n -= 1
   end
   println("Finished in $(round(t, digits=2))s.")
end


text:
module SharpArithmetic

export CompSum

export
    comp_sorted_add,
    comp_add,
    comp_sorted_sum,
    comp_sum,
    naive_sum,
    naive_sum2,
    value

struct CompSum{T}
    main::T
    err::T
end

value(x::CompSum) = x.main + x.err

@inline function comp_sorted_add(x::Number, y::Number)
    # @assert abs(x) >= abs(y)
    s = x + y
    e = (x - s) + y
    CompSum(s, e)
end

@inline function comp_add(x::Number, y::Number)
    if abs(x) < abs(y)
        x, y = y, x
    end
    comp_sorted_add(x, y)
end

@inline function comp_sorted_add(x::CompSum{<:Number}, y::Number)
    cs = comp_sorted_add(x.main, y)
    CompSum(cs.main, x.err + cs.err)
end

@inline function comp_add(x::CompSum{<:Number}, y::Number)
    cs = comp_add(x.main, y)
    CompSum(cs.main, x.err + cs.err)
end

function comp_sorted_sum(xs)
    cs = CompSum(zero(eltype(xs)), zero(eltype(xs)))
    for x in xs
        cs = comp_sorted_add(cs, x)
    end
    cs
end

function comp_sum(xs)
    cs = CompSum(zero(eltype(xs)), zero(eltype(xs)))
    for x in xs
        cs = comp_add(cs, x)
    end
    cs
end

function naive_sum(xs)
    s = zero(eltype(xs))
    for x in xs
        s += x
    end
    s
end

function naive_sum2(xs)
    sum(x for x in xs)
end

end


text:
using Documenter, KernelAbstractions

makedocs(
    modules = [KernelAbstractions],
    sitename = "KernelAbstractions",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages = [
        "Home" => "index.md",
        "Writing kernels" => "kernels.md",
        "Examples" => [
            "examples/memcopy.md"
            "examples/memcopy_static.md"
            "examples/naive_transpose.md"
            "examples/performance.md"
            "examples/matmul.md"
        ],
        "API"          => "api.md",
        "Extras" => [
            "extras/unrolling.md"
        ]
    ],
    doctest = true
)


text:
field = Mmap.mmap("extract_contour2",Array{Float64,3},(8192,8192,1024))

short = zeros(Float64,(1024,1024,128))

#I = N*i - N + 1
for k=1:128
  for j=1:1024
    for i=1:1024
      short[i,j,k] = field[8*i-7,8*j-7,8*k-7]
    end
  end
end

write("short_extract_contour2",short)


text:
module ReactiveRuntime

using WebSockets, MacroTools, Observables, JSServe
using JSServe.DOM
using JSServe: @js_str
export @variable

function Base.show(io::IO, m::MIME"application/vnd.webio.application+html", x::Observable)
    s = JSServe.with_session() do session
        JSServe.jsrender(session, x)
    end
    show(io, m, s)
end

const cell_outputs = Dict{Symbol, Any}()
const obs_func_names = Dict{Symbol, Symbol}()
const channels = Dict{Symbol, Channel}()

function collect_globals!(globals, expr, current)
    vars = Set{Symbol}()
    MacroTools.postwalk(expr) do var
        var isa Symbol && var != current && haskey(globals, var) && push!(vars, var)
    end
    return vars
end

struct OnUpdater{F}
    f::F
    name::Symbol
    observable::Observable{Any}
end

function (f::OnUpdater)(val)
    newval = f.f()
    if haskey(channels, f.name)
        close(pop!(channels, f.name))
    end
    if newval isa Channel
        channels[f.name] = newval
    end
    f.observable[] = newval
end

function lift_cell_reactive(expr)
    var = nothing
    body = if @capture(expr, var_ = begin; body__; end)
        Expr(:block, body...)
    elseif @capture(expr, var_ = simple_)
        simple
    else
        # Don't modify
        expr
    end
    obs_vars = collect_globals!(cell_outputs, body, var)
    varname = var === nothing ? gensym("tmp") : var
    funcname = get!(obs_func_names, varname, gensym(varname))
    var_quote = QuoteNode(varname)
    func_quote = QuoteNode(funcname)
    lifted = quote
        $(varname) = let
            function $(funcname)($(obs_vars...))
                $(body)
            end
            new_parents = ($(obs_vars...),)
            observable, old_parents = get!(ReactiveRuntime.cell_outputs, $(var_quote)) do
                ReactiveRuntime.Observable{Any}(nothing), new_parents
            end
            observable, old_parents = ReactiveRuntime.cell_outputs[$(var_quote)]
            new_updater = ReactiveRuntime.OnUpdater($(func_quote), observable) do
                 Base.invokelatest($(funcname), getindex.(new_parents)...)
            end
            # remove old listeners
            for old_parent in old_parents
                filter!(old_parent.listeners) do callback
                    !(callback isa ReactiveRuntime.OnUpdater && callback.name == $(func_quote))
                end
            end
            # add new listeners
            for new_parent in new_parents
                ReactiveRuntime.on(new_updater, new_parent)
            end
            # call the updater, to update value with the new function body
            new_updater(nothing)
            observable
        end
    end

    return lifted
end

macro variable(expr)
    expr = macroexpand(__module__, expr)
    return esc(lift_cell_reactive(expr))
end

end # module


text:
"""
`SolutionODE`: Solution of an ordinary differential equation

Contains all fields necessary to store the solution of an ODE.

### Fields

* `nd`: dimension of the dynamical variable ``q``
* `nt`: number of time steps to store
* `n0`: number of initial conditions
* `t`:  time steps
* `x`:  solution `x[nd, nt+1, n0]` with `x[:,0,:]` the initial conditions
* `ntime`: number of time steps to compute
* `nsave`: save every nsave'th time step

"""
immutable SolutionODE{T} <: Solution{T,3}
    nd::Int
    nt::Int
    n0::Int
    t::Timeseries{T}
    x::Array{T,3}
    ntime::Int
    nsave::Int

    function SolutionODE(nd, n0, ntime, nsave, Δt)
        @assert T <: Real
        @assert nd > 0
        @assert n0 > 0
        @assert nsave > 0
        @assert ntime ≥ nsave
        @assert mod(ntime, nsave) == 0

        nt = div(ntime, nsave)
        t = Timeseries{T}(nt, Δt, nsave)
        x = zeros(T, nd, nt+1, n0)
        new(nd, nt, n0, t, x, ntime, nsave)
    end
end

function SolutionODE{T}(equation::ODE{T}, Δt::T, ntime::Int, nsave::Int=1)
    s = SolutionODE{T}(equation.d, equation.n, ntime, nsave, Δt)
    set_initial_conditions!(s, equation)
    return s
end

function set_initial_conditions!(solution::SolutionODE, equation::ODE)
    set_initial_conditions!(solution, equation.t₀, equation.q₀)
end

function set_initial_conditions!(solution::SolutionODE, t₀, q₀)
    simd_copy_yx_first_last!(q₀, solution, 0)
    solution.t[0] = t₀
    compute_timeseries!(solution.t)
end

function reset(s::SolutionODE)
    for k in 1:solution.n0
        for i in 1:solution.nd
            solution[i, 0, k] = solution[i, end, k]
        end
    end
    solution.t[0] = solution.t[end]
    compute_timeseries!(solution.t)
end

Base.indices(s::SolutionODE) = (1:s.nd, 0:s.nt, 1:s.n0)
Base.strides(s::SolutionODE) = (1, s.nd, s.nd*s.nt)
# Base.linearindexing{T<:SolutionODE}(::Type{T}) = LinearFast()

@inline function Base.getindex(s::SolutionODE, i::Int, j::Int, k::Int)
    @boundscheck checkbounds(s.x, i, j+1, k)
    @inbounds r = getindex(s.x, i, j+1, k)
    return r
end

@inline function Base.getindex(s::SolutionODE, i::Int, j::Int)
    if s.n0 == 1
        @boundscheck checkbounds(s.x, i, j+1, 1)
        @inbounds r = getindex(s.x, i, j+1, 1)
    else
        @boundscheck checkbounds(s.x, 1:s.nd, i+1, j)
        @inbounds r = getindex(s.x, 1:s.nd, i+1, j)
    end
    return r
end

@inline function Base.getindex(s::SolutionODE, k::Int)
    if s.n0 == 1
        @boundscheck checkbounds(s.x, 1:s.nd, k+1, 1)
        @inbounds r = getindex(s.x, 1:s.nd, k+1, 1)
    else
        @boundscheck checkbounds(s.x, 1:s.nd, 1:s.nt, k)
        @inbounds r = getindex(s.x, 1:s.nd, 1:s.nt, k)
    end
    return r
end

@inline function Base.setindex!(s::SolutionODE, x, i::Int, j::Int, k::Int)
    @boundscheck checkbounds(s.x, i, j+1, k)
    @inbounds setindex!(s.x, x, i, j+1, k)
end

@inline function Base.setindex!(s::SolutionODE, x, j::Int, k::Int)
    if s.n0 == 1
        @assert length(x) == 1
        @boundscheck checkbounds(s.x, j, k+1, 1)
        @inbounds setindex!(s.x, x, j, k+1, 1)
    else
        @assert length(x) == s.nd
        @boundscheck checkbounds(s.x, :, j+1, k)
        @inbounds setindex!(s.x, x, :, j+1, k)
    end
end

@inline function Base.setindex!(s::SolutionODE, x, k::Int)
    @assert length(x) == s.nd*s.n0
    @boundscheck checkbounds(s.x, :, :, k)
    @inbounds setindex!(s.x, x, :, :, k)
end

"Creates HDF5 file and initialises datasets for ODE solution object."
function createHDF5{T}(solution::SolutionODE{T}, file::AbstractString, ntime::Int=1)
    @assert ntime ≥ 1

    info("Creating HDF5 file ", file)
    # TODO Put warning if file exists.
    h5 = h5open(file, "w")

    # create dataset and copy initial conditions
    # ntime can be used to set the expected total number of timesteps
    # so that the size of the array does not need to be adapted dynamically.
    # Right now, it has to be set as dynamical size adaptation is not yet
    # working.
    x = d_create(h5, "x", datatype(T), dataspace(solution.nd, ntime))
    x[1:solution.nd, 1] = solution[1:solution.nd, 0]

    return h5
end

"Append solution to HDF5 file."
function writeSolutionToHDF5(solution::SolutionODE, h5::HDF5.HDF5File, offset=0)
    # aquire dataset from HDF5 file
    x = h5["x"]

    # set convenience variables and compute ranges
    d  = solution.nd
    n  = solution.nt
    j1 = offset+2
    j2 = offset+1+n

    # # extend dataset if necessary
    # if size(x, 2) < j2
    #     set_dims!(x, (d, j2))
    # end

    # copy data from solution to HDF5 dataset
    x[1:d, j1:j2] = solution[1:d, 1:n]

    return nothing
end


text:
abstract type Colonization{R,W} <: SetCellRule{R,W} end

struct RandomColonization{R,W,C} <: Colonization{R,W}
    probability::C
end
function RandomColonization{R,W}(; probability = 0.1) where {R,W}
    RandomColonization{R,W}(probability)
end
function DynamicGrids.applyrule!(data, rule::RandomColonization{R,W}, N, I) where {R,W}
    # potentially type asset N is occupancy, later
    N > zero(N) || return zero(N)
    rand() < rule.probability && return one(N)
    return zero(N)
end

"""
    TODO
"""
struct SpatiallyExplicitLevinsColonization{R,W,C} <: Colonization{R,W}
    probability::C
end
function SpatiallyExplicitLevinsColonization{R,W}(; probability = 0.1) where {R,W}
    SpatiallyExplicitLevinsColonization{R,W}(probability)
end
function DynamicGrids.applyrule!(
    data,
    rule::SpatiallyExplicitLevinsColonization{R,W},
    N,
    I,
) where {R,W}
    # potentially type assert N is occupancy, later
    N > zero(N) || return zero(N)
    rand() < rule.probability && return one(N)
    return zero(N)
end

struct LevinsColonization{R,W,C} <: Colonization{R,W}
    probability::C
end
function LevinsColonization{R,W}(; probability = 0.1) where {R,W}
    LevinsColonization{R,W}(probability)
end
function DynamicGrids.applyrule!(data, rule::LevinsColonization{R,W}, O, I) where {R,W}
    # potentially type asset N is occupancy, later
    O > zero(O) || return zero(O)

    p = sum(data[:O]) / size(data[:O])
    c = rule.probability

    rand() < c * p * (1 - p) && return one(O)
    return zero(O)
end


text:
# the show function for general graphs

function show(io::IO, v::ExVertex)
    if isempty(v.label)
        print(io, "vertex [$(v.index)]")
    else
        print(io, "vertex [$(v.index)] \"$(v.label)\"")
    end
end

function show(io::IO, e::Union(Edge, ExEdge))
    print(io, "edge [$(e.index)]: $(e.source) -- $(e.target)")
end

function show(io::IO, graph::AbstractGraph)
    title = is_directed(graph) ? "Directed Graph" : "Undirected Graph"
    print(io, "$title ($(num_vertices(graph)) vertices, $(num_edges(graph)) edges)")
end

function show_details(io::IO, graph::AbstractGraph)
    title = is_directed(graph) ? "Directed Graph" : "Undirected Graph"
    println(io, "$title with $(num_vertices(graph)) vertices and $(num_edges(graph)) edges:")

    if !implements_vertex_list(graph)
        return
    end

    if implements_incidence_list(graph)
        println(io, "Incidence List:")
        for v in vertices(graph)
            println(io, "$(v): ")
            for e in out_edges(v, graph)
                print(io, "    ")
                println(io, e)
            end
        end

    elseif implements_adjacency_list(graph)
        println(io, "Adjacency List:")
        for v in vertices(graph)
            println(io, "$(v): ")
            for u in out_neighbors(v, graph)
                print(io, "    ")
                println(io, "$v -- $u")
            end
        end

    else
        println(io, "Vertices:")
        for v in vertices(graph)
            println(io, "$(v)")
        end

        if implements_edge_list(graph)
            println(io, "Edges:")
            for e in edges(graph)
                println(io, "$(e)")
            end
        end
    end
end


text:
using BinDeps

@BinDeps.setup

@unix_only begin
    libclp = library_dependency("libclp",aliases=["libClp"])
    libcoinmp = library_dependency("libcoinmp",aliases=["libCoinMP"])
end
@windows_only begin
    if Int != Int32
        error("Win64 platform is not yet supported by this package")
    end
    libclp = library_dependency("libclp",aliases=["CoinMP"])
    libcoinmp = library_dependency("libcoinmp",aliases=["CoinMP"])
end

coinmpname = "CoinMP-1.7.6"

provides(Sources, URI("http://www.coin-or.org/download/source/CoinMP/$coinmpname.tgz"),
    [libclp,libcoinmp], os = :Unix)

provides(Binaries, URI("http://www.mit.edu/~mlubin/CoinMP_julia_20130903.tar.gz"),
    [libclp,libcoinmp], os = :Windows)

@osx_only begin
    using Homebrew
    provides( Homebrew.HB, "coinmp", [libclp, libcoinmp], os = :Darwin )
end

prefix=joinpath(BinDeps.depsdir(libclp),"usr")
patchdir=BinDeps.depsdir(libclp)
srcdir = joinpath(BinDeps.depsdir(libclp),"src",coinmpname)

provides(SimpleBuild,
    (@build_steps begin
        GetSources(libclp)
        @build_steps begin
            ChangeDirectory(srcdir)
            `cat $patchdir/CoinMP-emptyproblem.patch` |> `patch -N -p1`
            `./configure --prefix=$prefix --enable-dependency-linking`
            `make install`
        end
    end),[libclp,libcoinmp], os = :Unix)

@BinDeps.install [:libclp => :libclp, :libcoinmp => :libcoinmp]


text:
using Test, Joseki
import HTTP.IOExtras.bytes

@testset "Middleware" begin
    req = sample_request()
    # Make sure this doesn't throw an error
    Joseki.require_json_body!(req)

    # Now make a bad request and make sure it does throw
    bad_body = "{definitely not valid json}"
    req = HTTP.Request("GET", "", [], bytes(bad_body))
    @test_throws Exception Joseki.require_json_body!(req)
end


text:
using Gnuplot, LinearAlgebra
using ParameterizedFunctions, DifferentialEquations
g = @ode_def begin
  dx = σ*(y - x)
  dy = x*(ρ - z) - y
  dz = x*y - β*z
end σ ρ β
u0 = [1.0; 0.0; 0.0]
tspan = (0.0, 80.0)
p = [10.0,28.0,8/3]
prob = ODEProblem(g, u0, tspan, p)
sol = solve(prob, Tsit5(), saveat = 0.01)
x, y, z = sol[1,:], sol[2,:], sol[3,:]
tempo = sol.t

@gsp "set xyplane at -3" "set auto fix" "set grid" :-
@gsp :- "set title 'Lorenz attractor'" "set cblabel 'time'"
@gsp :- xlab = "x" ylab = "y"
@gsp :- xr=extrema(x) yr=extrema(y) zr=extrema(z)
# animation
@gsp :- "i=2"
@gsp :- x y z tempo "every ::1::i w l notit lc palette" :-
for frame in 2:40:length(x)
  @gsp :- frame "i=$frame" "replot" :-
end
@gsp # display animation
save(term="gif animate size 1000,800 delay 1", output="lorenzAnim.gif")

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
struct Bounds
    lb
    ub

    function Bounds(lb, ub)
        lb = isnothing(lb) ? () : lb
        ub = isnothing(ub) ? () : ub

        size(lb) == () && size(ub) == () && !isnothing(lb) && !isnothing(ub) && return new(lb, ub)

        sizes = map(x -> size(x), [lb, ub])
        nonempty_bounds = findall(sizes .!= Ref(()))
        @assert !isempty(nonempty_bounds)
        length = sizes[nonempty_bounds[1]]

        lb = 1 in nonempty_bounds ? lb : repeat([-Inf], length)
        ub = 2 in nonempty_bounds ? ub : repeat([-Inf], length)

        new(lb, ub)
    end
end

abstract type PSProblem end

@with_kw mutable struct InitializedProblem{T} <: PSProblem
    dynamics
    path_constraints
    cost
    x₀::AbstractArray{T}
    u₀::AbstractArray{T}
    δ
    state_bounds
    control_bounds
end

@with_kw mutable struct UninitializedProblem <: PSProblem
    dynamics
    path_constraints
    cost
    n_states
    n_controls
    n
    δ
    state_bounds::Bounds
    control_bounds::Bounds
    end_points
end

UninitializedProblem(prob::InitializedProblem) = UninitializedProblem(prob.dynamics, prob.path_constraints, prob.cost, size(prob.x₀)[1], size(prob.u₀)..., prob.δ, prob.end_points)

function symbolic_to_str(symbolic_expr::Basic; extend_dims=true, interpolate_variables=false)
    str_sym = string(symbolic_expr)

    # TODO: make it generic over variable names
    regex_dims = r"(\w+)(\d)\[(\d+)\]"
    regex_slice = r"([X|U]\[.*?\])"

    if extend_dims
        str_sym = replace(str_sym, regex_dims => s"\1[\2,\3]")
    end
    if interpolate_variables
        str_sym = replace(str_sym, regex_slice => s"$(\1)")
    end

    str_sym
end

substitute_args(expr, vars) = expr
substitute_args(expr::Symbol, vars) = get(vars, expr, expr)
function substitute_args(expr::Expr, vars)
    for (i, arg) in enumerate(expr.args)
        if arg in keys(vars)
            expr.args[i] = vars[arg]
        else
            expr.args[i] = substitute_args(arg, vars)
        end
    end
    return expr
end

function ode_constraint_exprs(M, dynamics, n, δ, ns, nc)

    function ode_error(dynamics, X, U)
        X * M' .- dynamics(X, U)
    end

    X_symbolic = zeros(Basic, (ns, n + 1))
    U_symbolic = zeros(Basic, (nc, n + 1))
    for (is, ic) in zip(1:ns, 1:nc), i in 1:n + 1
        X_symbolic[is, i] = symbols("X$is[$i]")[1]
        U_symbolic[ic, i] = symbols("U$ic[$i]")[1]
    end

    err = ode_error(dynamics, X_symbolic, U_symbolic)
    strs = symbolic_to_str.(err; interpolate_variables=false)
    strs .*= "≤ $δ"

    Meta.parse.(strs)
end

arr2d_to_arrlist(arr) = getindex.(Ref(arr), :, 1:size(arr)[2])
arrlist_to_arr2d(arr) = hcat(arr...)

function build_problem(prob::UninitializedProblem)
    n = prob.n
    τ, w = gausslobatto(n + 1)

    M = zeros(Float64, (n + 1, n + 1))
    diff_matrix!(M, n)

    xb = prob.state_bounds
    ub = prob.control_bounds
    broadcast_dynamics(X, U) = arrlist_to_arr2d(prob.dynamics.(arr2d_to_arrlist.([X, U])...))
    exprs = ode_constraint_exprs(M, broadcast_dynamics, n, prob.δ, prob.n_states, prob.n_controls)

    model = Model(Ipopt.Optimizer)

    @variables model begin
        xb.lb[i] <= X[i=1:prob.n_states, 1:n + 1] <= xb.ub[i], (base_name = "X$i")
        ub.lb[i] <= U[i=1:prob.n_controls, 1:n + 1] <= ub.ub[i], (base_name = "U$i")
    end

    @objective model Min prob.cost(X, U)
    @constraint model path_constraints 0 .<= prob.path_constraints(X, U)
    @constraint model initial_point X[:, 1] .== prob.end_points[1]
    @constraint model end_point X[:, end] .== prob.end_points[2]

    vars = Dict()
    for (is, ic) in zip(1:prob.n_states, 1:prob.n_controls), j in 1:n + 1
        vars[Meta.parse("X[$is, $j]")] = X[is, j]
        vars[Meta.parse("U[$ic, $j]")] = U[ic, j]
    end

    spliced_exprs = substitute_args.(exprs, Ref(vars))

    for expr in spliced_exprs
        add_NL_constraint(model, expr)
    end

    optimize!(model)

    X_sol = value.(X)
    U_sol = value.(U)

    cost_sol = prob.cost(X_sol, U_sol)

    X_sol, U_sol, cost_sol, model, τ
end


text:
function crr(O::Option; N::Int64 = 1000)
    Δt = O.t / N
    U = exp(O.σ * √Δt)
    D = 1 / U
    R = exp(O.r * Δt)
    p = (R - D) / (U - D)
    q = (U - R) / (U - D)

    if O.call_or_put == "Put"
        Z = [max(0, O.k - O.s * exp((2 * i - N) * O.σ * √Δt)) for i = 0:N]
    else
        Z = [max(0, O.s * exp((2 * i - N) * O.σ * √Δt) - O.k) for i = 0:N]
    end

    for n = N-1:-1:0
        for i = 0:n
            if O.call_or_put == "Put"
                x = O.k - O.s * exp((2 * i - n) * O.σ * √Δt)
            else
                x = O.s * exp((2 * i - n) * O.σ * √Δt) - O.k
            end
            y = (q * Z[i+1] + p * Z[i+2]) / R
            Z[i+1] = max(x, y)
        end
    end

    return Z[1]
end


text:
using MLDatasets:Titanic

function get_processed_data(args)
    labels = Titanic.targets()
    features = Titanic.features()

    # Split into training and test sets, 2/3 for training, 1/3 for test.
    train_indices = [1:3:891; 2:3:891]

    X_train = features[:, train_indices]
    y_train = labels[:, train_indices]

    X_test = normed_features[:, 3:3:891]
    y_test = onehot_labels[:, 3:3:891]

    #repeat the data `args.repeat` times
    train_data = Iterators.repeated((X_train, y_train), args.repeat)
    test_data = (X_test, y_test)

    return train_data, test_data
end


text:
using Clang
using Test

# only test this on Linux since "time.h" is missing on other platforms by default.
if get(ENV, "TRAVIS_OS_NAME", "") == "linux"

# LIBCLANG_HEADERS are those headers to be wrapped.
const LIBCLANG_INCLUDE = joinpath(@__DIR__, "..", "deps", "usr", "include", "clang-c") |> normpath
const LIBCLANG_HEADERS = [joinpath(LIBCLANG_INCLUDE, header) for header in readdir(LIBCLANG_INCLUDE) if endswith(header, ".h")]

@testset "generator" begin
    api_path = joinpath(@__DIR__, "..", "gen", "libclang_api.jl") |> normpath
    common_path = joinpath(@__DIR__, "..", "gen", "libclang_common.jl") |> normpath
    test_api_path = joinpath(@__DIR__, "test_api.jl")
    test_common_path = joinpath(@__DIR__, "test_common.jl")
    wc = init(; headers = LIBCLANG_HEADERS,
                output_file = test_api_path,
                common_file = test_common_path,
                clang_includes = vcat(LIBCLANG_INCLUDE, CLANG_INCLUDE),
                clang_args = ["-I", joinpath(LIBCLANG_INCLUDE, "..")],
                header_wrapped = (root, current)->root == current,
                header_library = x->"libclang",
                clang_diagnostics = true)
    run(wc)
    @test read(test_api_path, String) == read(api_path, String)
    @test read(test_common_path, String) == read(common_path, String)
end

end


text:
@Generated 1 immutable Cell{Name, T}
    if !isa(Name, Symbol)
        str = "Cell parameter 1 (Name) is expected to be a symbol, got $Name"
        return :(error(str))
    end

    if !isa(T, DataType)
        str = "Cell parameter 2 (T) is expected to be a DataType, got $T"
        return :(error(str))
    end

    :(Name::T)
end

function Base.call{Name, T}(::Type{Cell{Name}}, data::T)
    Cell{Name,T}(data)
end

Base.call{Name, T, Tnew}(::Type{Cell{Name,Tnew}}, cell::Cell{Name,T})    = Cell{Name,Tnew}(convert(Tnew, cell.(1)))
Base.convert{Name, Name_new, T, T_new}(::Type{Cell{Name_new,T_new}}, cell::Cell{Name,T}) = Cell{Name_new,T_new}(convert(T_new, cell.(1)))

@inline colname{Name}(::Cell{Name}) = Name
@inline colname{Name, T}(::Type{Cell{Name,T}}) = Name
@inline colname{C <: Cell}(::Type{C}) = colname(super(C))
@inline Base.eltype{Name, T}(::Cell{Name,T}) = T
@inline Base.eltype{Name, T}(::Type{Cell{Name,T}}) = T
@inline Base.eltype{C <: Cell}(::Type{C}) = eltype(super(C))

@inline nrow(::Cell) = 1
@inline ncol(::Cell) = 1

Base.start(cell::Cell) = false
Base.next(cell::Cell, state) = (cell.(1), true)
Base.done(cell::Cell, state) = state
Base.endof(cell::Cell) = 1
Base.length(cell::Cell) = 1

@inline getindex(c::Cell) = c.(1)
@inline getindex(c::Cell, i) = i == 1 ? c.(1) : error("Cannot index Cell at $i")

# copy
Base.copy{Name}(cell::Cell{Name}) = Cell{Name}(copy(cell.(1)))


text:
# ! FIXME: It should be safe to return internal structures as long as
# noone is going to modify them. I see no such scenario.

@doc raw"""
    Directional.l2(x :: CorrelationTracker, phase)

Return $L_2^{\text{phase}}$ function for an underlying system of the
tracker `x`.
"""
Directional.l2(tracker :: CorrelationTracker{T}, phase) where T =
    tracker.corrdata[L2Tracker{T}(phase)]

@doc raw"""
    Directional.s2(x :: CorrelationTracker, phase)

Return $S_2^{\text{phase}}$ function for an underlying system of the
tracker `x`.
"""
Directional.s2(tracker :: CorrelationTracker{T}, phase) where T =
    tracker.corrdata[S2Tracker{T}(phase)]

@doc raw"""
    Directional.surfsurf(x :: CorrelationTracker, phase)

Return $F_{ss}_2^{\text{phase}}$ function for an underlying system of the
tracker `x`.
"""
Directional.surfsurf(tracker :: CorrelationTracker{T}, phase) where T =
    tracker.corrdata[SSTracker{T}(phase)]

@doc raw"""
    Directional.surfvoid(x :: CorrelationTracker, phase)

Return $F_{sv}_2^{\text{phase}}$ function for an underlying system of the
tracker `x`.
"""
Directional.surfvoid(tracker :: CorrelationTracker{T}, phase) where T =
    tracker.corrdata[SVTracker{T}(phase)]

const TrackerAlike{T, N} = Union{CorrelationTracker{T, N}, ExtrapolatedData{T, N}}

"""
    CorrelationTracker(system, like)

Create a `CorrelationTracker` on top of array `system` with parameters
matching those in the tracker `like`.
"""
function CorrelationTracker(system :: AbstractArray{T, N},
                            like   :: TrackerAlike{T, N}) where {T, N}
    return CorrelationTracker(copy(system);
                              tracking   = like |> tracked_data |> collect,
                              periodic   = like.periodic,
                              directions = tracked_directions(like),
                              len        = tracked_length(like))
end

# Make AbstractTracker callable for convenience
(tracked :: L2Tracker{T})(tracker :: TrackerAlike{T}) where T =
    tracker.corrdata[tracked]
(tracked :: S2Tracker{T})(tracker :: TrackerAlike{T}) where T =
    tracker.corrdata[tracked]
(tracked :: SSTracker{T})(tracker :: TrackerAlike{T}) where T =
    tracker.corrdata[tracked]
(tracked :: SVTracker{T})(tracker :: TrackerAlike{T}) where T =
    tracker.corrdata[tracked]
(tracked :: CCTracker{T})(tracker :: TrackerAlike{T}) where T =
    tracker.corrdata[tracked]

"""
    cross_correlation(tracker :: CorrelationTracker, phase1, phase2)

Return precalculated values of cross correlation function for phases
`phase1` and `phase2`.
"""
cross_correlation(tracker :: CorrelationTracker{T}, phase1, phase2) where T =
    tracker.corrdata[CCTracker{T}(phase1, phase2)]


text:
# A cache for results from line search methods (to avoid recomputation)
type LineSearchResults{T}
    alpha::Vector{T}
    value::Vector{T}
    slope::Vector{T}
    nfailures::Int
end

LineSearchResults{T}(::Type{T}) = LineSearchResults(T[], T[], T[], 0)

Base.length(lsr::LineSearchResults) = length(lsr.alpha)

function Base.push!{T}(lsr::LineSearchResults{T}, a::T, v::T, d::T)
    push!(lsr.alpha, a)
    push!(lsr.value, v)
    push!(lsr.slope, d)
    return
end

function clear!(lsr::LineSearchResults)
    empty!(lsr.alpha)
    empty!(lsr.value)
    empty!(lsr.slope)
    return
    # nfailures is deliberately not set to 0
end

type LinesearchException{T<:Real} <: Exception
    message::AbstractString
    alpha::T
    f_update::Int
    g_update::Int
    lsr::LineSearchResults
end

immutable DifferentiableFunction
    f::Function
    g!::Function
    fg!::Function
end

function DifferentiableFunction(f::Function, g!::Function)
    function fg!(x::Array, storage::Array)
        g!(x, storage)
        return f(x)
    end
    return DifferentiableFunction(f, g!, fg!)
end


text:
@export function scr_test(sim; λ=DEFAULT_λ, indices=[37], weights=[1])
    ρ0 = sim.γL
    ρ0 = sum(zip(indices, weights)) do (i, w)
        v0 = eigen(sim.H).vectors[:,i]
        w*v0*v0'
    end
    Γ = sim.ΓR

    J = sim.H - im*Γ/2
    j, U = eigen(J)
    N = size(J, 1)

    _ρ_ = inv(U)*ρ0*inv(U')
    _Γ_ = U'Γ*U
    _L_ = [λ * U'L*inv(U') for L in sim.L]

    # function β(t)
    #     As = Matrix{Complex{Float64}}[]
    #     for L in _L_
    #         A = Matrix{Complex{Float64}}(undef, N, N)

    #         for idx in CartesianIndices(A)
    #             l, m = idx[1], idx[2]
    #             if j[l] == j[m]
    #                 A[l,m] = -im * t * exp(-im*j[l]*t) * L[l,m]
    #             else
    #                 A[l,m] = inv(j[l]- j[m]) * (exp(-im*j[l]*t)*L[l,m] - L[l,m]*exp(-im*j[m]*t))
    #             end
    #         end

    #         push!(As, U*A*inv(U))
    #     end
    #     return As
    # end

    function f(t)
        integrated_T = sum(CartesianIndices((N, N))) do idx
            k, l = idx[1], idx[2]
            x = j[k] - j[l]'
            if iszero(x)
                return _ρ_[k,l] * _Γ_[l,k] * t
            else
                return _ρ_[k,l] * _Γ_[l,k] * im*(exp(-im*x*t) - 1) / x
            end
        end

        integrated_S = sum(CartesianIndices((N, N, N))) do idx
            l, m, n = idx[1], idx[2], idx[3]
            x = j[l]' - j[m]'
            if iszero(x)
                a = j[n] - j[m]'
                return _ρ_[n,l] * _L_[3][l,m] * _Γ_[m,n] * im*(
                    im*t*exp(-im*a*t) / a +
                    (exp(-im*a*t) - 1) / a^2
                )
            else
                a = j[n] - j[l]'
                b = j[n] - j[m]'
                return _ρ_[n,l] * _L_[3][l,m] * _Γ_[m,n] * inv(b - a) * im*(
                    (exp(-im*a*t) - 1) / a - (exp(-im*b*t) - 1) / b
                )
                # return _ρ_[n,l] * _L_[3][l,m] * _Γ_[m,n] * inv(j[l]' - j[m]') * im*(
                #     (exp(-im*(j[n] - j[l]')*t) - 1) / (j[n] - j[l]') -
                #     (exp(-im*(j[n] - j[m]')*t) - 1) / (j[n] - j[m]')
                # )
            end
        end

        return integrated_T, integrated_S
        u = LinearAlgebra.exp!(-im*J*t)
        A = ρ0*u'Γ
        return tr(A*u), integrated_T, [tr(A*b) for b in β(t)], integrated_S
    end

    return f
end

@export function scr_test_full(sim; λ=DEFAULT_λ, indices=[Int(length(sim.eig.values) / 2 + 1)], weights=[1], α=1.0)
    ρ0 = sum(zip(indices, weights)) do (i, w)
        v0 = eigen(sim.H).vectors[:,i]
        w*v0*v0'
    end
    Γ = sim.ΓL
    Σ = - im * (α*sim.ΓR + sim.ΓL) / 2

    J = (sim.H + Σ) ⊗ σ0 + sum(λ * map(⊗, sim.L, σ))
    j, U = eigen(J)

    _ρ_ = inv(U) * (ρ0 ⊗ σ0) * inv(U')
    _Γ_ = U' * (Γ ⊗ σ0) * U
    _Γi_ = [U' * (Γ ⊗ σi) * U for σi in σ]

    N = size(J, 1)

    function f(t)
        integrated_T = sum(CartesianIndices((N, N))) do idx
            k, l = idx[1], idx[2]
            x_kl = im*(j[k] - j[l]')
            if iszero(x_kl)
                return _ρ_[k,l] * _Γ_[l,k] * t
            else
                return _ρ_[k,l] * _Γ_[l,k] * (1 - exp(-x_kl * t)) / x_kl
            end
        end

        integrated_S = sum(CartesianIndices((N, N))) do idx
            k, l = idx[1], idx[2]
            x_kl = im*(j[k] - j[l]')
            if iszero(x_kl)
                return _ρ_[k,l] * _Γi_[3][l,k] * t
            else
                return _ρ_[k,l] * _Γi_[3][l,k] * (1 - exp(-x_kl * t)) / x_kl
            end
        end

        # u = LinearAlgebra.exp!(-im*J*t)
        # A = (ρ0 ⊗ σ0) * u' * (Γ ⊗ σ0)
        return integrated_T, integrated_S#tr(A*u)
    end

    return f
end

# Magnet
@export function scr_test_magnet(sim; λ=DEFAULT_λ, indices=[Int(length(sim.eig.values) / 2 + 1)], weights=[1], α=1.0, m=0.5)
    ρ0 = sum(zip(indices, weights)) do (i, w)
        v0 = eigen(sim.H).vectors[:,i]
        w*v0*v0'
    end
    Γ = α*sim.ΓR
    Σ = - im * (α*sim.ΓR + sim.ΓL) / 2
    ΔΣ = - im * (m*sim.γL) / 2

    J = ((sim.H + Σ) ⊗ σ0 + ΔΣ ⊗ σ[:z] + sum(λ * map(⊗, sim.L, σ)),
         (sim.H + Σ) ⊗ σ0 - ΔΣ ⊗ σ[:z] + sum(λ * map(⊗, sim.L, σ)))
    j, U = begin
        tmp = eigen.(J)
        map(x -> x.values, tmp), map(x -> x.vectors, tmp)
    end

    ρ = map(U -> inv(U) * (ρ0 ⊗ σ0) * inv(adjoint(U)), U)
    Γ = map(U -> adjoint(U) * (Γ ⊗ σ0) * U, U)

    N = size(J[1], 1)

    function f(t)
        integrated_T = map(j, ρ, Γ) do j, ρ, Γ
            T = sum(CartesianIndices((N, N))) do idx
                k, l = idx[1], idx[2]
                x_kl = im*(j[k] - j[l]')
                real(x_kl) < 0 && error("negative `real(x_kl)` found: $(real(x_kl))")
                return ρ[k,l] * Γ[l,k] * (1 - exp(-x_kl * t)) / x_kl
            end
            real(T) < 0 && error("negative `real(T)` found: $(real(T))")
            real(T)
        end

        return integrated_T
    end

    return f
end

# Calculate the photon emission
@export function scr_photon_emission_full(
    sim;
    λ=DEFAULT_λ,
    indices=[37],
    weights=[1],
    seed = 824,
    W = 1e-1,
    rng = MersenneTwister(seed)
)
    ρ0 = sum(zip(indices, weights)) do (i, w)
        v0 = sim.eig.vectors[:,i]
        w*v0*v0'
    end
    Γ = sim.ΓR

    J = (sim.H - im*Γ/2) ⊗ σ0 + sum(λ * Li ⊗ σi for (Li, σi) in zip(sim.L, σ))
    j, U = eigen(J)
    N = size(J, 1)

    V = let
        local i, j = rand(rng, 1:size(sim.H, 1), 2)
        while i == j i, j = rand(rng, 1:N, 2) end
        W*sim.eig.vectors[:,i]*sim.eig.vectors[:,j]'
    end

    _ρ_ = inv(U) * (ρ0 ⊗ σ0) * inv(U')
    _Γ_ = U' * (Γ ⊗ σ0) * U
    _V_ = inv(U) * (V ⊗ σ0) * U

    function f()
        K(a, b) = inv(j[a]' - j[b])
        K_(a, b) = inv(j[a]' - j[b]')

        output = sum(CartesianIndices((N, N, N, N))) do idx
            k, l, m, n = idx[1], idx[2], idx[3], idx[4]

            X = 2π*im * _ρ_[k,l] * _V_'[l,m] * _Γ_[m,n] * _V_[m,n]

            if iszero(j[l]' - j[m]')
                return -X * K(l,n) * K(m,k) * (K(l,n) + K(m,k))
            else
                return X * K_(l,m) * (K(l,n)*K(l,k) - K(m,n)*K(m,k))
            end
        end

        return output
    end

    return f
end

##################################################

@export function scr_magnet(sim; E=sim.μ, λ=DEFAULT_λ, m=0.5, n=1, α=0.0, ϕ=0)
    ρ0 = let
        v0 = eigenstate(sim, E = E, ϕ=ϕ)
        0.5 * v0*v0'
    end

    Γ = sim.ΓL
    Σ = - im * (α*sim.ΓL + sim.ΓR) / 2
    ΔΣ = - im * (m*sim.ΔΓR) / 2

    J = ((sim.H + Σ) ⊗ σ0 + ΔΣ ⊗ σ[:z] + sum(λ * map(⊗, sim.L, σ)),
         (sim.H + Σ) ⊗ σ0 - ΔΣ ⊗ σ[:z] + sum(λ * map(⊗, sim.L, σ)))
    j, U = begin
        tmp = eigen.(J)
        map(x -> x.values, tmp), map(x -> x.vectors, tmp)
    end

    ρ = map(U -> inv(U) * (ρ0 ⊗ σ0) * inv(adjoint(U)), U)
    Γ = map(U -> adjoint(U) * (Γ ⊗ σ0) * U, U)

    N = size(J[1], 1)

    Ts = map(j, ρ, Γ) do j, ρ, Γ
        T = sum(CartesianIndices((N, N))) do idx
            k, l = idx[1], idx[2]
            x_kl = im*(j[k] - j[l]')
            real(x_kl) < 0 && error("negative `real(x_kl)` found: $(real(x_kl))")
            return ρ[k,l] * Γ[l,k] / x_kl
        end
        real(T) < 0 && error("negative `real(T)` found: $(real(T))")
        real(T)
    end

    return (Ts = Ts,)
end

@export function time_dependent_escape(
    sim;
    E = sim.μ,
    λ = DEFAULT_λ,
    m = 1.0,
    α = DEFAULT_α,
    ϕ = 0,
    hbar = HBAR
)
    ρ0 = let
        v0 = eigenstate(sim, E = E, ϕ=ϕ)
        v0*v0'
    end

    Γ = α*sim.ΓR
    Σ = - im * α*sim.ΓR / 2

    J = (sim.H + Σ) ⊗ σ0 + sum(λ * map(⊗, sim.L, σ))
    j, U = eigen(J)
    N = length(j)

    ρ = inv(U) * (ρ0 ⊗ ((σ0 + m * σ[:z]) / 2)) * inv(adjoint(U))
    Γ = adjoint(U) * (Γ ⊗ σ0) * U

    ξs = map(CartesianIndices((N, N))) do idx
        k, l = idx[1], idx[2]
        x_kl = im*(j[k] - j[l]')
        real(x_kl) < 0 && error("negative `real(x_kl)` found: $(real(x_kl))")
        return x_kl
    end
    Ts = map(CartesianIndices((N, N))) do idx
        k, l = idx[1], idx[2]
        x_kl = im*(j[k] - j[l]')
        real(x_kl) < 0 && error("negative `real(x_kl)` found: $(real(x_kl))")
        return ρ[k,l] * Γ[l,k] / x_kl
    end
    fs = map(CartesianIndices((N, N))) do idx
        k, l = idx[1], idx[2]
        x_kl = im*(j[k] - j[l]')
        real(x_kl) < 0 && error("negative `real(x_kl)` found: $(real(x_kl))")
        return f(t) = 1 - exp(-x_kl * t / hbar)
    end

    return (ξs = ξs, Ts = Ts, fs = fs)
end

@export function scr_f(
    δx = 0.0,
    δy = 0.1;
    t0 = 0,
    t1 = 1,
    nsamples = 500,
    ts = range(t0, t1, length=nsamples),
    seed = 42,
    α = 1,
)
    x0, y0, y1 = let
        xsymbol = :δz
        x0 = 0.8
        y0 = -4.2
        mol = Helicene(N=7)

        x0, y0, y1 = TheoryOfCISS.find_crossing(
            x0, y0, mol,
            xsymbol = xsymbol,
            threshold=1e-12,
            attenuation=1e-2,
            rounding = :both
        )
    end

    x = x0
    y = (y0 + y1) / 2

    cols_up = gen_bands(Helicene(N=7), xsymbol = :δz, bounds=(x + δx, x + δx), nsamples=1, seed=seed);
    cols_down = gen_bands(Helicene(N=7), xsymbol = :δz, bounds=(x + δx, x + δx), nsamples=1, seed=seed);

    gen_on_bands!(cols_up, ysymbol=:E, f = time_dependent_escape, m=1.0, bounds=(y - δy, y + δy), α=α);
    gen_on_bands!(cols_down, ysymbol=:E, f = time_dependent_escape, m=-1.0, bounds=(y - δy, y + δy), α=α);

    _f(t, x) = sum(map((T, f) -> T * f(t), x.Ts, x.fs))
    f_x(t, i) = real(_f(t, cols_up[1].data[i])), real(_f(t, cols_down[1].data[i]))
    f_xs = map(i -> mapreduce(t -> [f_x(t, i)...]', vcat, ts), 1:2)
    pol(a, b) = 100 * !iszero(a - b) * ((a - b) / (a + b))
    dat = map(x -> map(i -> pol(i...), eachrow(x)), f_xs)
    return ts, f_xs, dat
end


text:
# Worked with Kasra #
cd("D:\\Semester 3_Fall 2020\\ECON 6453 - Econometrics III\\ProblemSets\\PS6-ccp")
using Random, Distributions, JLD2, DataFrames, CSV, LinearAlgebra, Statistics, FreqTables, TexTables, Optim, HTTP, GLM, ForwardDiff, DataFramesMeta

include("create_grids.jl")

function problemset6()

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 1
#:::::::::::::::::::::::::::::::::::::::::::::::::::
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS5-ddc/busdata.csv"
df = CSV.read(HTTP.get(url).body)

df = @transform(df, bus_id = 1:size(df,1))

#---------------------------------------------------
# reshape from wide to long (must do this twice be-
# cause DataFrames.stack() requires doing it one
# variable at a time)
#---------------------------------------------------
# first reshape the decision variable
dfy = @select(df, :bus_id,:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20,:RouteUsage,:Branded)
dfy_long = DataFrames.stack(dfy, Not([:bus_id,:RouteUsage,:Branded]))
rename!(dfy_long, :value => :Y)
dfy_long = @transform(dfy_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfy_long, Not(:variable))

# next reshape the odometer variable
dfx = @select(df, :bus_id,:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20)
dfx_long = DataFrames.stack(dfx, Not([:bus_id]))
rename!(dfx_long, :value => :Odometer)
dfx_long = @transform(dfx_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfx_long, Not(:variable))

# join reshaped df's back together
df_long = leftjoin(dfy_long, dfx_long, on = [:bus_id,:time])
sort!(df_long,[:bus_id,:time])

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 2
#:::::::::::::::::::::::::::::::::::::::::::::::::::
theta=glm(@formula(Y ~ Odometer * Odometer^2 * RouteUsage * RouteUsage^2 * Branded * time * time^2), df_long, Binomial(), LogitLink())

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 3
#:::::::::::::::::::::::::::::::::::::::::::::::::::
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS5-ddc/busdata.csv"
df2 = CSV.read(HTTP.get(url).body)
Y = Matrix(df2[:,[:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20]])
X = Matrix(df2[:,[:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20]])
Z = Vector(df2[:,:RouteUsage])
B = Vector(df2[:,:Branded])
N = size(Y,1)
T = size(Y,2)
Xstate = Matrix(df2[:,[:Xst1,:Xst2,:Xst3,:Xst4,:Xst5,:Xst6,:Xst7,:Xst8,:Xst9,:Xst10,:Xst11,:Xst12,:Xst13,:Xst14,:Xst15,:Xst16,:Xst17,:Xst18,:Xst19,:Xst20]])
Zstate = Vector(df2[:,:Zst])
zval,zbin,xval,xbin,xtran = create_grids()

log_df=DataFrame()
log_df.Odometer=kron(ones(zbin),xval) #20301
log_df.RouteUsage =kron(ones(xbin),zval) #20301
log_df.Branded=zeros((size(log_df,1)))
log_df.time=zeros(size(log_df,1))

function values(din,Zstate,Xstate,xtran,zbin,xbin,xval)
FV=zeros(size(xtran,1),2,T+1)
for t=2:T
for b=0:1
din.time[t]=t
din.Branded[t]=b
p0 = predict(θ, din)

FV[:, b+1, t] = - .9 .*log.(p0)
end
end
FVT1 = zeros(size(df,1), T)
for i=1:size(df,1)
row0 = (Zstate[i]-1)*xbin+1
for t=1:T
row1  = Xstate[i,t] + (Zstate[i]-1)*xbin
FVT1[i,t] = (xtran[row1,:].-xtran[row0,:])'*FV[row0:row0+xbin-1,df.Branded[i]+1,t+1]
end
end

return FVT1'[:]
end

fvt1=values(log_df,Zstate,Xstate,xtran,zbin,xbin,xval)
df_long = @transform(df_long,fv = fvt1)

theta_hat_ccp_glm = glm(@formula(Y ~ Odometer + Branded),df_long, Binomial(), LogitLink(),offset=df_long.fv)

end

@time problemset6()

max_stars_repo_path: ProblemSets/PS6_Ahmed Chaudhry.jl

max_issues_repo_path: ProblemSets/PS6_Ahmed Chaudhry.jl

max_forks_repo_path: ProblemSets/PS6_Ahmed Chaudhry.jl


text:
export FileHeader

struct FileHeader
    th::String
    bfh::BinaryFileHeader
end

function show(io::IO, fh::FileHeader)
    show(fh.bfh)
end

FileHeader() = FileHeader(" "^3200, BinaryFileHeader())


text:
if VERSION >= v"0.7.0-beta.0"
    using Random
end
using Compat
import NCDatasets

varname = "varname"
filename = tempname()

samples = [
    # Ints
    [1,2,3,4],

    # Floats
    [1. 2. 3.; 4. 5. 6.],

    # chars
    ['a','b','c'],
    ['a' 'b' 'c'; 'd' 'e' 'f'],

    # strings
    ["wieso","weshalb","warum"],
    ["wieso" "weshalb" "warum"; "why" "why" "whyyy"],
    map(x -> randstring(rand(3:10)), zeros(2,3,4)),
]

for sampledata in samples
    rm(filename;force=true)

    # write data
    ncid = NCDatasets.nc_create(filename,NCDatasets.NC_CLOBBER | NCDatasets.NC_NETCDF4)

    dimids = zeros(Cint,ndims(sampledata))
    for i = 1:ndims(sampledata)
        dimids[i] = NCDatasets.nc_def_dim(ncid, "dim-$(i)", size(sampledata,i))
    end

    T = eltype(sampledata)
    xtype = NCDatasets.ncType[T]
    # reverse order
    varid = NCDatasets.nc_def_var(ncid, varname, xtype, reverse(dimids))
    NCDatasets.nc_put_att(ncid, varid, "attr-string-list",["one","two"])
    NCDatasets.nc_put_var(ncid, varid, sampledata)
    NCDatasets.nc_close(ncid)

    # load data

    ncid = NCDatasets.nc_open(filename,NCDatasets.NC_NOWRITE)
    varid = NCDatasets.nc_inq_varid(ncid,varname)
    xtype2 = NCDatasets.nc_inq_vartype(ncid,varid)
    @test xtype == xtype

    attrval = NCDatasets.nc_get_att(ncid, varid, "attr-string-list")
    @test attrval == ["one","two"]

    sampledata2 = Array{T,ndims(sampledata)}(undef,size(sampledata))
    NCDatasets.nc_get_var!(ncid,varid,sampledata2)
    @test sampledata == sampledata2

    # start = [1,1]
    # count = [size(

    # NCDatasets.nc_get_vara!(ncid,varid,start,count,sampledata2)
    # @test sampledata == sampledata2

    NCDatasets.nc_close(ncid)
end

# accept SubString as file name argument

ncid = NCDatasets.nc_create(filename,NCDatasets.NC_CLOBBER | NCDatasets.NC_NETCDF4)
NCDatasets.nc_close(ncid)

ncid = NCDatasets.nc_open(split("$(filename)#foo",'#')[1],NCDatasets.NC_NOWRITE)
NCDatasets.nc_close(ncid)


text:
# https://adventofcode.com/2020/day/15

example = "0,3,6"

function read_array(string::AbstractString)::Vector{Int}
    return parse.(Int, split(string, ',', keepempty=false))
end

"""
Use array as a lookup table since keys are integers and we know the
upper bound. Access to array's elements is O(1).
"""
function memoize_sequence(starting_sequence, turns)
    n = max(turns + 1, maximum(starting_sequence))
    memory = zeros(Int32, n) .- (2 * n)

    for (i, prev) in enumerate(starting_sequence[1:(end - 1)])
        memory[prev + 1] = i
    end

    prev = starting_sequence[end]
    pos = 0

    for i = length(starting_sequence):(turns - 1)
        pos = memory[prev + 1]
        memory[prev + 1] = i
        prev = pos > 0 ? i - pos : 0
    end
    return prev
end

function part1(string, turns=2020)
    starting_sequence = read_array(string)
    return memoize_sequence(starting_sequence, turns)
end

@assert part1(example, 10) == 0
@assert part1(example) == 436
@assert part1("1,3,2") == 1
@assert part1("2,1,3") == 10
@assert part1("1,2,3") == 27
@assert part1("2,3,1") == 78
@assert part1("3,2,1") == 438
@assert part1("3,1,2") == 1836

part2(string) = part1(string, 30000000)

@assert part2("0,3,6") == 175594
@assert part2("1,3,2") == 2578
@assert part2("2,1,3") == 3544142
@assert part2("1,2,3") == 261214
@assert part2("2,3,1") == 6895259
@assert part2("3,2,1") == 18
@assert part2("3,1,2") == 362

test = "0,13,16,17,1,10,6"
println("Part 1: $(result1 = part1(test))")
println("Part 2: $(result2 = part2(test))")

@assert result1 == 276
@assert result2 == 31916

max_stars_repo_name: twolodzko/advent-of-code

max_issues_repo_name: twolodzko/advent-of-code

max_forks_repo_name: twolodzko/advent-of-code


text:
function apply!(::NOP_0, organism::TierrianOrganism, model::TierraModel) end
function apply!(::NOP_1, organism::TierrianOrganism, model::TierraModel) end

function apply!(::OR_1, organism::TierrianOrganism, model::TierraModel) organism.c = organism.c | 1 end
function apply!(::SH_L, organism::TierrianOrganism, model::TierraModel) organism.c = organism.c << 1 end
function apply!(::ZERO, organism::TierrianOrganism, model::TierraModel) organism.c = 0 end

function apply!(::IF_0, organism::TierrianOrganism, model::TierraModel)
    if (organism.c != 0) organism.ip = _address_within(organism, organism.ip + one(UInt16)) end
end

function apply!(::SUB_AB, organism::TierrianOrganism, model::TierraModel) organism.c = organism.a - organism.b end
function apply!(::SUB_AC, organism::TierrianOrganism, model::TierraModel) organism.a = organism.a - organism.c end

function apply!(::INC_A, organism::TierrianOrganism, model::TierraModel) organism.a += one(UInt16) end
function apply!(::INC_B, organism::TierrianOrganism, model::TierraModel) organism.b += one(UInt16) end
function apply!(::INC_C, organism::TierrianOrganism, model::TierraModel) organism.c += one(UInt16) end
function apply!(::DEC_C, organism::TierrianOrganism, model::TierraModel) organism.c -= one(UInt16) end

function apply!(::PUSH_A, organism::TierrianOrganism, model::TierraModel) push!(organism.stack, organism.a) end
function apply!(::PUSH_B, organism::TierrianOrganism, model::TierraModel) push!(organism.stack, organism.b) end
function apply!(::PUSH_C, organism::TierrianOrganism, model::TierraModel) push!(organism.stack, organism.c) end
function apply!(::PUSH_D, organism::TierrianOrganism, model::TierraModel) push!(organism.stack, organism.d) end

function apply!(::POP_A, organism::TierrianOrganism, model::TierraModel) organism.a = _pop_stack!(organism) end
function apply!(::POP_B, organism::TierrianOrganism, model::TierraModel) organism.b = _pop_stack!(organism) end
function apply!(::POP_C, organism::TierrianOrganism, model::TierraModel) organism.c = _pop_stack!(organism) end
function apply!(::POP_D, organism::TierrianOrganism, model::TierraModel) organism.d = _pop_stack!(organism) end

function apply!(::JMP, organism::TierrianOrganism, model::TierraModel)
    organism.ip, _ = search_template(model, organism, organism.ip, SearchBoth())
end
function apply!(::JMP_BACK, organism::TierrianOrganism, model::TierraModel)
    organism.ip, _ = search_template(model, organism, organism.ip, SearchBackward())
end

function apply!(::CALL, organism::TierrianOrganism, model::TierraModel)
    callee_address, _ = search_template(model, organism, organism.ip, SearchBoth())
    push!(organism.stack, organism.ip)
    organism.ip = callee_address
end
function apply!(::RET, organism::TierrianOrganism, model::TierraModel) organism.ip = _pop_stack!(organism) end

function apply!(::MOV_AB, organism::TierrianOrganism, model::TierraModel) organism.b = organism.a end
function apply!(::MOV_CD, organism::TierrianOrganism, model::TierraModel) organism.d = organism.c end
function apply!(::MOVI_BA, organism::TierrianOrganism, model::TierraModel)
    instruction = read_memory(model, organism.b)
    write_memory(model, organism, organism.a, instruction)
end

function apply!(::ADR_A, organism::TierrianOrganism, model::TierraModel)
    organism.a, organism.c = search_template(model, organism, organism.ip, SearchBoth())
end
function apply!(::ADR_BACK_A, organism::TierrianOrganism, model::TierraModel)
    organism.a, organism.c = search_template(model, organism, organism.ip, SearchBackward())
end
function apply!(::ADR_FORW_A, organism::TierrianOrganism, model::TierraModel)
    organism.a, organism.c = search_template(model, organism, organism.ip, SearchForward())
end

function apply!(::MALLOC, organism::TierrianOrganism, model::TierraModel) allocate_daughter(model, organism) end
function apply!(::DIVIDE, organism::TierrianOrganism, model::TierraModel) divide(model, organism) end

function _pop_stack!(organism::TierrianOrganism)
    if isempty(organism.stack)
        organism.error_flag = true
        return 0
    else
        return pop!(organism.stack)
    end
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
#######################
#### bitwise shift left
#######################

function shift_left!(a::ArbUInt, shift::T) where T <: Integer
    shift < zero(T) && return shift_right!(a, abs(shift))
    iszero(shift) && return a
    bits = T(BITS)
    digits = shift ÷ bits
    digits > typemax(UInt) && throw(DomainError(shift, "requested shift too large to fit in memory"))
    digits = UInt(digits)
    shift = UInt8(shift % bits)
    _shift_left!(a, digits, shift)
end
const shl! = shift_left!

function _shift_left!(a::ArbUInt, digits::UInt, shift::UInt8)
    data = if iszero(digits)
        a.data
    else
        # we may crash here due to lack of memory, nothing we can do about that
        Base._growbeg0!(a.data, signed(digits)) # REMOVEME: once we can grow with guaranteed zeroed memory, since this is internal
        a.data
    end

    digits += one(digits) # one-based adjustment
    if shift > 0
        carry = zero(eltype(data))
        carry_shift = UInt8(BITS) - shift
        for i in digits:lastindex(data)
            elem = data[i]
            new_carry = elem >> carry_shift
            data[i] = (elem << shift) | carry
            carry = new_carry
        end
        !iszero(carry) && push!(data, carry)
    end

    normalize!(a) # unlike rust, we modify the original - not sure if that's a good idea
end

Base.:(<<)(a::ArbUInt, b::T) where T <: Integer = shift_left!(deepcopy(a), b)
Base.:(<<)(a::ArbUInt, b::T) where T <: Unsigned = shift_left!(deepcopy(a), b) # ambiguity
Base.:(<<)(_, b::ArbUInt) = throw(ArgumentError("can't leftshift by ArbUInt"))
Base.:(<<)(_::ArbUInt, _::ArbUInt) = throw(ArgumentError("can't leftshift ArbUInt by another ArbUInt"))

########################
#### bitwise shift right
########################

function shift_right!(a::ArbUInt, shift::T) where T <: Integer
    shift < zero(T) && return shift_left!(a, abs(shift))
    iszero(shift) && return a
    bits = T(BITS)
    more_digits = shift ÷ bits
    more_digits = more_digits > typemax(UInt) ? typemax(UInt) : UInt(more_digits)
    shift = UInt8(shift % bits)
    _shift_right!(a, more_digits, shift)
end
const shr! = shift_right!

function _shift_right!(a::ArbUInt, digits::UInt, shift::UInt8)
    digits >= sizeof(ArbDigit)*length(a.data) && return set_zero!(a)
    data = @view a.data[digits+1:end]

    if shift > 0
        borrow = 0
        borrow_shift = UInt8(BITS) - shift
        for i in lastindex(data):-1:1
            elem = data[i]
            new_borrow = elem << borrow_shift
            data[i] = (elem >> shift) | borrow
            borrow = new_borrow
        end
    end

    # maybe dangerous when aliased?
    a.data[begin:length(data)] .= data
    resize!(a.data, length(data))
    normalize!(a) # unlike rust, we modify the original
end

# this makes a surprising amount of sense
function Base.:(>>)(a::T, b::ArbUInt) where T <: Integer
    iszero(b) && return a
    if length(b.data) == 1
        return a >> b.data[begin]
    else
        return zero(T)
    end
end
Base.:(>>)(a::ArbUInt, b::T) where T <: Integer = shift_right!(deepcopy(a), b)
Base.:(>>)(_::ArbUInt, _::ArbUInt) = throw(ArgumentError("can't rightshift ArbUInt by another ArbUInt"))