text:
# Specify contrasts for coding categorical data in model matrix. Contrasts types
# are a subtype of AbstractContrasts. ContrastsMatrix types hold a contrast
# matrix, levels, and term names and provide the interface for creating model
# matrix columns and coefficient names.
#
# Contrasts types themselves can be instantiated to provide containers for
# contrast settings (currently, just the base level).
#
# ModelFrame will hold a Dict{Symbol, ContrastsMatrix} that maps column
# names to contrasts.
#
# ModelMatrix will check this dict when evaluating terms, falling back to a
# default for any categorical data without a specified contrast.

"""
Interface to describe contrast coding systems for categorical variables.

Concrete subtypes of `AbstractContrasts` describe a particular way of converting a
categorical data vector into numeric columns in a `ModelMatrix`. Each
instantiation optionally includes the levels to generate columns for and the base
level. If not specified these will be taken from the data when a `ContrastsMatrix` is
generated (during `ModelFrame` construction).

# Constructors

For `C <: AbstractContrast`:

```julia
C()                                     # levels are inferred later
C(levels = ::Vector{Any})               # levels checked against data later
C(base = ::Any)                         # specify base level
C(levels = ::Vector{Any}, base = ::Any) # specify levels and base
```

# Arguments

* `levels`: Optionally, the data levels can be specified here.  This allows you
  to specify the order of the levels.  If specified, the levels will be checked
  against the levels actually present in the data when the `ContrastsMatrix` is
  constructed. Any mismatch will result in an error, because levels missing in
  the data would lead to empty columns in the model matrix, and levels missing
  from the contrasts would lead to empty or undefined rows.
* `base`: The base level may also be specified.  The actual interpretation
  of this depends on the particular contrast type, but in general it can be
  thought of as a "reference" level.  It defaults to the first level.

# Contrast coding systems

* [`DummyCoding`](@ref) - Code each non-base level as a 0-1 indicator column.
* [`EffectsCoding`](@ref) - Code each non-base level as 1, and base as -1.
* [`HelmertCoding`](@ref) - Code each non-base level as the difference from the
  mean of the lower levels
* [`SeqDiffCoding`](@ref) - Code for differences between sequential levels of
  the variable.
* [`HypothesisCoding`](@ref) - Manually specify contrasts via a hypothesis
  matrix, which gives the weighting for the average response for each level
* [`StatsModels.ContrastsCoding`](@ref) - Manually specify contrasts matrix,
  which is directly copied into the model matrix.

The last two coding types, `HypothesisCoding` and `StatsModels.ContrastsCoding`,
provide a way to manually specify a contrasts matrix. For a variable `x` with
`k` levels, a contrasts matrix `M` is a `k×k-1` matrix, that maps the `k` levels
onto `k-1` model matrix columns.  Specifically, let `X` be the full-rank
indicator matrix for `x`, where `X[i,j] = 1` if `x[i] == levels(x)[j]`, and 0
otherwise. Then the model matrix columns generated by the contrasts matrix `M`
are `Y = X * M`.

The *hypothesis matrix* is the `k-1×k` matrix that gives the weighted
combinations of group mean responses that are represented by regression
coefficients for the generated contrasts.  The contrasts matrix is the
generalized pseudo-inverse (e.g. `LinearAlgebra.pinv`) of the hypothesis matrix.
See [`HypothesisCoding`](@ref) or Schad et al. (2020) for more information.

# Extending

The easiest way to specify custom contrasts is with `HypothesisCoding` or
`StatsModels.ContrastsCoding`.  But if you want to actually implement a custom
contrast coding system, you can subtype `AbstractContrasts`.  This requires a
constructor, a `contrasts_matrix` method for constructing the actual contrasts
matrix that maps from levels to `ModelMatrix` column values, and (optionally) a
`termnames` method:

```julia
mutable struct MyCoding <: AbstractContrasts
    ...
end

contrasts_matrix(C::MyCoding, baseind, n) = ...
termnames(C::MyCoding, levels, baseind) = ...
```

# References

Schad, D. J., Vasishth, S., Hohenstein, S., & Kliegl, R. (2020). How to
capitalize on a priori contrasts in linear (mixed) models: A tutorial. _Journal
of Memory and Language, 110_, 104038. https://doi.org/10.1016/j.jml.2019.104038

"""
abstract type AbstractContrasts end

# Contrasts + Levels (usually from data) = ContrastsMatrix
struct ContrastsMatrix{C <: AbstractContrasts, T, U}
    matrix::Matrix{Float64}
    termnames::Vector{U}
    levels::Vector{T}
    contrasts::C
    invindex::Dict{T,Int}
    function ContrastsMatrix(matrix::AbstractMatrix,
                             termnames::Vector{U},
                             levels::Vector{T},
                             contrasts::C) where {U,T,C <: AbstractContrasts}
        invindex = Dict{T,Int}(x=>i for (i,x) in enumerate(levels))
        new{C,T,U}(matrix, termnames, levels, contrasts, invindex)
    end
end

# only check equality of matrix, termnames, and levels, and that the type is the
# same for the contrasts (values are irrelevant).  This ensures that the two
# will behave identically in creating modelmatrix columns
Base.:(==)(a::ContrastsMatrix{C,T}, b::ContrastsMatrix{C,T}) where {C<:AbstractContrasts,T} =
    a.matrix == b.matrix &&
    a.termnames == b.termnames &&
    a.levels == b.levels

Base.hash(a::ContrastsMatrix{C}, h::UInt) where {C} =
    hash(C, hash(a.matrix, hash(a.termnames, hash(a.levels, h))))

"""
An instantiation of a contrast coding system for particular levels

This type is used internally for generating model matrices based on categorical
data, and **most users will not need to deal with it directly**.  Conceptually,
a `ContrastsMatrix` object stands for an instantiation of a contrast coding
*system* for a particular set of categorical *data levels*.

If levels are specified in the `AbstractContrasts`, those will be used, and likewise
for the base level (which defaults to the first level).

# Constructors

```julia
ContrastsMatrix(contrasts::AbstractContrasts, levels::AbstractVector)
ContrastsMatrix(contrasts_matrix::ContrastsMatrix, levels::AbstractVector)
```

# Arguments

* `contrasts::AbstractContrasts`: The contrast coding system to use.
* `levels::AbstractVector`: The levels to generate contrasts for.
* `contrasts_matrix::ContrastsMatrix`: Constructing a `ContrastsMatrix` from
  another will check that the levels match.  This is used, for example, in
  constructing a model matrix from a `ModelFrame` using different data.

"""
function ContrastsMatrix(contrasts::C, levels::AbstractVector{T}) where {C<:AbstractContrasts, T}

    # if levels are defined on contrasts, use those, validating that they line up.
    # what does that mean? either:
    #
    # 1. DataAPI.levels(contrasts) == levels (best case)
    # 2. data levels missing from contrast: would generate empty/undefined rows.
    #    better to filter data frame first
    # 3. contrast levels missing from data: would have empty columns, generate a
    #    rank-deficient model matrix.
    c_levels = something(DataAPI.levels(contrasts), levels)
    if eltype(c_levels) != eltype(levels)
        throw(ArgumentError("mismatching levels types: got $(eltype(levels)), expected " *
                            "$(eltype(c_levels)) based on contrasts levels."))
    end
    mismatched_levels = symdiff(c_levels, levels)
    if !isempty(mismatched_levels)
        throw(ArgumentError("contrasts levels not found in data or vice-versa: " *
                            "$mismatched_levels." *
                            "\n  Data levels: $levels." *
                            "\n  Contrast levels: $c_levels"))
    end

    n = length(c_levels)
    if n == 0
        throw(ArgumentError("empty set of levels found (need at least two to compute " *
                            "contrasts)."))
    elseif n == 1
        throw(ArgumentError("only one level found: $(c_levels[1]) (need at least two to " *
                            "compute contrasts)."))
    end

    # find index of base level. use contrasts.base, then default (1).
    base_level = baselevel(contrasts)
    baseind = base_level === nothing ?
              1 :
              findfirst(isequal(base_level), c_levels)
    if baseind === nothing
        throw(ArgumentError("base level $(base_level) not found in levels " *
                            "$c_levels."))
    end

    tnames = termnames(contrasts, c_levels, baseind)

    mat = contrasts_matrix(contrasts, baseind, n)

    ContrastsMatrix(mat, tnames, c_levels, contrasts)
end

ContrastsMatrix(c::Type{<:AbstractContrasts}, levels::AbstractVector) =
    throw(ArgumentError("contrast types must be instantiated (use $c() instead of $c)"))

# given an existing ContrastsMatrix, check that all passed levels are present
# in the contrasts. Note that this behavior is different from the
# ContrastsMatrix constructor, which requires that the levels be exactly the same.
# This method exists to support things like `predict` that can operate on new data
# which may contain only a subset of the original data's levels. Checking here
# (instead of in `modelmat_cols`) allows an informative error message.
function ContrastsMatrix(c::ContrastsMatrix, levels::AbstractVector)
    if !isempty(setdiff(levels, c.levels))
         throw(ArgumentError("there are levels in data that are not in ContrastsMatrix: " *
                             "$(setdiff(levels, c.levels))" *
                             "\n  Data levels: $(levels)" *
                             "\n  Contrast levels: $(c.levels)"))
    end
    return c
end

function termnames(C::AbstractContrasts, levels::AbstractVector, baseind::Integer)
    not_base = [1:(baseind-1); (baseind+1):length(levels)]
    levels[not_base]
end

Base.getindex(contrasts::ContrastsMatrix{C,T}, rowinds, colinds) where {C,T} =
    getindex(contrasts.matrix, getindex.(Ref(contrasts.invindex), rowinds), colinds)

# Making a contrast type T only requires that there be a method for
# contrasts_matrix(T,  baseind, n) and optionally termnames(T, levels, baseind)
# The rest is boilerplate.
for contrastType in [:DummyCoding, :EffectsCoding, :HelmertCoding, :SeqDiffCoding]
    @eval begin
        mutable struct $contrastType <: AbstractContrasts
            base::Any
            levels::Union{AbstractVector,Nothing}
        end
        ## constructor with optional keyword arguments, defaulting to nothing
        $contrastType(; base=nothing, levels::Union{AbstractVector,Nothing}=nothing) = $contrastType(base, levels)
        baselevel(c::$contrastType) = c.base
        DataAPI.levels(c::$contrastType) = c.levels
    end
end

# fallback method for other types that might not have base or level fields
baselevel(c::AbstractContrasts) = nothing
DataAPI.levels(c::AbstractContrasts) = nothing

"""
    FullDummyCoding()

Full-rank dummy coding generates one indicator (1 or 0) column for each level,
**including** the base level. This is sometimes known as
[one-hot encoding](https://en.wikipedia.org/wiki/One-hot).

Not exported but included here for the sake of completeness.
Needed internally for some situations where a categorical variable with ``k``
levels needs to be converted into ``k`` model matrix columns instead of the
standard ``k-1``.  This occurs when there are missing lower-order terms, as in
discussed below in [Categorical variables in Formulas](@ref).

# Examples

```julia
julia> StatsModels.ContrastsMatrix(StatsModels.FullDummyCoding(), ["a", "b", "c", "d"]).matrix
4×4 Array{Float64,2}:
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0
```
"""
mutable struct FullDummyCoding <: AbstractContrasts
# Dummy contrasts have no base level (since all levels produce a column)
end

ContrastsMatrix(C::FullDummyCoding, levels::AbstractVector{T}) where {T} =
    ContrastsMatrix(Matrix(1.0I, length(levels), length(levels)), levels, levels, C)

"Promote contrasts matrix to full rank version"
Base.convert(::Type{ContrastsMatrix{FullDummyCoding}}, C::ContrastsMatrix) =
    ContrastsMatrix(FullDummyCoding(), C.levels)

"""
    DummyCoding([base[, levels]])
    DummyCoding(; base=nothing, levels=nothing)

Dummy coding generates one indicator column (1 or 0) for each non-base level.

If `levels` are omitted or `nothing`, they are determined from the data
by calling the `levels` function on the data when constructing `ContrastsMatrix`.
If `base` is omitted or `nothing`, the first level is used as the base.

Columns have non-zero mean and are collinear with an intercept column (and
lower-order columns for interactions) but are orthogonal to each other. In a
regression model, dummy coding leads to an intercept that is the mean of the
dependent variable for base level.

Also known as "treatment coding" or "one-hot encoding".

# Examples

```julia
julia> StatsModels.ContrastsMatrix(DummyCoding(), ["a", "b", "c", "d"]).matrix
4×3 Array{Float64,2}:
 0.0  0.0  0.0
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0
```
"""
DummyCoding

contrasts_matrix(C::DummyCoding, baseind, n) =
    Matrix(1.0I, n, n)[:, [1:(baseind-1); (baseind+1):n]]

"""
    EffectsCoding([base[, levels]])
    EffectsCoding(; base=nothing, levels=nothing)

Effects coding generates columns that code each non-base level as the
deviation from the base level.  For each non-base level `x` of `variable`, a
column is generated with 1 where `variable .== x` and -1 where `variable .== base`.

`EffectsCoding` is like `DummyCoding`, but using -1 for the base level instead
of 0.

If `levels` are omitted or `nothing`, they are determined from the data
by calling the `levels` function when constructing `ContrastsMatrix`.
If `base` is omitted or `nothing`, the first level is used as the base.

When all levels are equally frequent, effects coding generates model matrix
columns that are mean centered (have mean 0).  For more than two levels the
generated columns are not orthogonal.  In a regression model with an
effects-coded variable, the intercept corresponds to the grand mean.

Also known as "sum coding" or "simple coding". Note
though that the default in R and SPSS is to use the *last* level as the base.
Here we use the *first* level as the base, for consistency with other coding
systems.

# Examples

```julia
julia> StatsModels.ContrastsMatrix(EffectsCoding(), ["a", "b", "c", "d"]).matrix
4×3 Array{Float64,2}:
 -1.0  -1.0  -1.0
  1.0   0.0   0.0
  0.0   1.0   0.0
  0.0   0.0   1.0
```

"""
EffectsCoding

function contrasts_matrix(C::EffectsCoding, baseind, n)
    not_base = [1:(baseind-1); (baseind+1):n]
    mat = Matrix(1.0I, n, n)[:, not_base]
    mat[baseind, :] .= -1
    return mat
end

"""
    HelmertCoding([base[, levels]])
    HelmertCoding(; base=nothing, levels=nothing)

Helmert coding codes each level as the difference from the average of the lower
levels.

If `levels` are omitted or `nothing`, they are determined from the data
by calling the `levels` function when constructing `Contrastsmatrix`.
If `base` is omitted or `nothing`, the first level is used as the base.
For each non-base level, Helmert coding generates a columns with -1 for each of
n levels below, n for that level, and 0 above.

When all levels are equally frequent, Helmert coding generates columns that are
mean-centered (mean 0) and orthogonal.

# Examples

```julia
julia> StatsModels.ContrastsMatrix(HelmertCoding(), ["a", "b", "c", "d"]).matrix
4×3 Array{Float64,2}:
 -1.0  -1.0  -1.0
  1.0  -1.0  -1.0
  0.0   2.0  -1.0
  0.0   0.0   3.0
```
"""
HelmertCoding

function contrasts_matrix(C::HelmertCoding, baseind, n)
    mat = zeros(n, n-1)
    for i in 1:n-1
        mat[1:i, i] .= -1
        mat[i+1, i] = i
    end

    # re-shuffle the rows such that base is the all -1.0 row (currently first)
    mat = mat[[baseind; 1:(baseind-1); (baseind+1):end], :]
    return mat
end

"""
    SeqDiffCoding([base[, levels]])

Code each level in order to test "sequential difference" hypotheses, which
compares each level to the level below it (starting with the second level).
Specifically, the ``n``th predictor tests the hypothesis that the difference
between levels ``n`` and ``n+1`` is zero.

Differences are computed in order of `levels`.  If `levels` are omitted or
`nothing`, they are determined from the data by calling the `levels` function
when constructing `ContrastsMatrix`.
If `base` is omitted or `nothing`, the first level is used as the base.

# Examples

```jldoctest seqdiff
julia> seqdiff = StatsModels.ContrastsMatrix(SeqDiffCoding(), ["a", "b", "c", "d"]).matrix
4×3 Array{Float64,2}:
 -0.75  -0.5  -0.25
  0.25  -0.5  -0.25
  0.25   0.5  -0.25
  0.25   0.5   0.75
```

The interpretation of sequential difference coding may be hard to see from the
contrasts matrix itself.  The corresponding hypothesis matrix shows a clearer
picture.  From the rows of the hypothesis matrix, we can see that these
contrasts test the difference between the first and second levels, the second
and third, and the third and fourth, respectively:

```jldoctest seqdiff
julia> StatsModels.hypothesis_matrix(seqdiff)
3×4 Array{Int64,2}:
 -1   1   0  0
  0  -1   1  0
  0   0  -1  1
```

"""
SeqDiffCoding

function contrasts_matrix(C::SeqDiffCoding, baseind, n)
    mat = zeros(n, n-1)
    for col in 1:n-1
        mat[1:col, col] .= col-n
        mat[col+1:end, col] .= col
    end
    return mat ./ n
end

# TODO: consider customizing term names:
# termnames(C::SeqDiffCoding, levels::AbstractVector, baseind::Integer) =
#     ["$(levels[i])-$(levels[i-1])" for i in 2:length(levels)]

"""
    HypothesisCoding(hypotheses::AbstractMatrix; levels=nothing, labels=nothing)

Specify how to code a categorical variable in terms of a *hypothesis matrix*.
For a variable with ``k`` levels, this should be a ``k-1 \times k`` matrix.
Each row of the matrix corresponds to a hypothesis about the mean
outcomes under each of the ``k`` levels of the predictor.  The entries in the
row give the weights assigned to each of these ``k`` means, and the
corresponding predictor in a regression model estimates the weighted sum of
these cell means.

For instance, if we have a variable which has four levels A, B, C, and D, and we
want to test the hypothesis that the difference between the average outcomes for
levels A and B is different from zero, the corresponding row of the hypothesis
matrix would be `[-1, 1, 0, 0]`.  Likewise, to test whether the difference
between B and C is different from zero, the hypothesis vector would be `[0, -1,
1, 0]`.  To test each "successive difference" hypothesis, the full hypothesis
matrix would be

```jldoctest hyp
julia> sdiff_hypothesis = [-1  1  0  0
                            0 -1  1  0
                            0  0 -1  1];
```

Contrasts are derived the hypothesis matrix by taking the pseudoinverse:

```jldoctest hyp
julia> using LinearAlgebra

julia> sdiff_contrasts = pinv(sdiff_hypothesis)
4×3 Array{Float64,2}:
 -0.75  -0.5  -0.25
  0.25  -0.5  -0.25
  0.25   0.5  -0.25
  0.25   0.5   0.75
```

The above matrix is what is produced by constructing a [`ContrastsMatrix`](@ref) from a
`HypothesisCoding` instance:

```jldoctest hyp
julia> StatsModels.ContrastsMatrix(HypothesisCoding(sdiff_hypothesis), ["a", "b", "c", "d"]).matrix
4×3 Array{Float64,2}:
 -0.75  -0.5  -0.25
  0.25  -0.5  -0.25
  0.25   0.5  -0.25
  0.25   0.5   0.75
```

The interpretation of the such "sequential difference" contrasts are clear when
expressed as a hypothesis matrix, but it is not obvious just from looking at the
contrasts matrix.  For this reason `HypothesisCoding` is preferred for
specifying custom contrast coding schemes over `ContrastsCoding`.

Optional arguments `levels` and `labels` give the names (in order) of the
hypothesis matrix columns (corresponding to levels of the data) and rows
(corresponding to the tested hypothesis).  The `labels` also determine the names
of the model matrix columns generated by these contrasts.

# References

Schad, D. J., Vasishth, S., Hohenstein, S., & Kliegl, R. (2020). How to
capitalize on a priori contrasts in linear (mixed) models: A tutorial. _Journal
of Memory and Language, 110_, 104038. https://doi.org/10.1016/j.jml.2019.104038

"""
mutable struct HypothesisCoding{T<:AbstractMatrix, S<:AbstractMatrix} <: AbstractContrasts
    hypotheses::T
    contrasts::S
    levels::Union{AbstractVector,Nothing}
    labels::Union{AbstractVector,Nothing}

    function HypothesisCoding(hypotheses::T;
                              levels::Union{AbstractVector,Nothing}=nothing,
                              labels::Union{AbstractVector,Nothing}=nothing) where {T<:AbstractMatrix}
        labels == nothing &&
            Base.depwarn(
                "HypothesisCoding without specified contrast labels is deprecated.  " *
                "Specify contrast labels with `HypothesisCoding(; labels=[...])` " *
                "or HypothesisCoding(Dict(label1=>hyp1, label2=>hyp2, ...))",
                :HypothesisCoding)
        contrasts = pinv(hypotheses)
        S = typeof(contrasts)
        check_contrasts_size(contrasts, levels)
        new{T,S}(hypotheses, contrasts, levels, labels)
    end
end

"""
    HypothesisCoding(hypotheses::Vector{Pair}; levels=nothing)

Specify hypotheses as a vector of `label=>hypothesis_vector` pairs.  Labels are
generated from the keys of the pairs.  The `levels` argument provides the names
of the levels for each entry in the hypothesis vectors.  If omitted, the
`levels` function will be called on the data when constructing a
`ContrastsMatrix`

# Example

```jldoctest
julia> hc = HypothesisCoding(["a_and_b" => [1, 1, 0], "b_and_c" => [0, 1, 1]]);

julia> hc.hypotheses
2×3 Array{Int64,2}:
 1  1  0
 0  1  1

julia> hc.labels
2-element Array{String,1}:
 "a_and_b"
 "b_and_c"
```
"""
function HypothesisCoding(hypotheses::Vector{<:Pair{<:Any,<:AbstractVector}};
                          levels::Union{AbstractVector,Nothing}=nothing)
    labels = first.(hypotheses)
    mat = reduce(vcat, adjoint.(last.(hypotheses)))
    HypothesisCoding(mat; labels=labels, levels=levels)
end

function contrasts_matrix(C::HypothesisCoding, baseind, n)
    check_contrasts_size(C.contrasts, n)
    C.contrasts
end

termnames(C::HypothesisCoding, levels::AbstractVector, baseind::Int) =
    something(C.labels, levels[1:length(levels) .!= baseind])

DataAPI.levels(c::HypothesisCoding) = c.levels

"""
    StatsModels.ContrastsCoding(mat::AbstractMatrix[, levels]])
    StatsModels.ContrastsCoding(mat::AbstractMatrix[; levels=nothing])

Coding by manual specification of contrasts matrix. For k levels, the contrasts
must be a k by k-1 Matrix.  The contrasts in this matrix will be copied directly
into the model matrix; if you want to specify your contrasts as hypotheses (i.e.,
weights assigned to each level's cell mean), you should use
[`HypothesisCoding`](@ref) instead.
"""
mutable struct ContrastsCoding{T<:AbstractMatrix} <: AbstractContrasts
    mat::T
    levels::Union{Vector,Nothing}

    function ContrastsCoding(mat::T, levels::Union{AbstractVector,Nothing}) where {T<:AbstractMatrix}
        Base.depwarn("`ContrastsCoding(contrasts)` is deprecated and will not be" *
                     " exported in the future.  Future versions will require" *
                     " `StatsModels.ContrastsCoding` or `using StatsModels: " *
                     "ContrastsCoding`.  For general users we recommend " *
                     "`HypothesisCoding(pinv(contrasts))` instead.",
                     :ContrastsCoding)
        check_contrasts_size(mat, levels)
        new{T}(mat, levels)
    end
end

check_contrasts_size(mat::Matrix, ::Nothing) = check_contrasts_size(mat, size(mat,1))
check_contrasts_size(mat::Matrix, levels::Vector) = check_contrasts_size(mat, length(levels))
check_contrasts_size(mat::Matrix, n_lev::Int) =
    size(mat) == (n_lev, n_lev-1) ||
    throw(ArgumentError("contrasts matrix wrong size for $n_lev levels. " *
                        "Expected $((n_lev, n_lev-1)), got $(size(mat))"))

## constructor with optional keyword arguments, defaulting to nothing
ContrastsCoding(mat::AbstractMatrix; levels::Union{AbstractVector,Nothing}=nothing) =
    ContrastsCoding(mat, levels)

function contrasts_matrix(C::ContrastsCoding, baseind, n)
    check_contrasts_size(C.mat, n)
    C.mat
end

DataAPI.levels(c::ContrastsCoding) = c.levels

## hypothesis matrix
"""
    needs_intercept(mat::AbstractMatrix)

Returns `true` if a contrasts matrix `mat` needs to consider the intercept to
determine the corresponding hypothesis matrix.  This is determined by whether
`rank(mat)` is less than the number of rows and any column is not orthogonal to
the intercept column.
"""
needs_intercept(mat::AbstractMatrix) =
    (rank(mat) < size(mat, 1)) &&
    !all(isapprox.(sum(mat, dims=1), 0, rtol=1e-5))

"""
    hypothesis_matrix(cmat::AbstractMatrix; intercept=needs_intercept(cm), tolerance=1e-5)
    hypothesis_matrix(contrasts::AbstractContrasts, n; baseind=1, kwargs...)
    hypothesis_matrix(cmat::ContrastsMatrix; kwargs...)

Compute the hypothesis matrix for a contrasts matrix using the generalized
pseudo-inverse (`LinearAlgebra.pinv`).  `intercept` determines whether a column
of ones is included before taking the pseudoinverse, which is needed for
contrasts where the columns are not orthogonal to the intercept (e.g., have
non-zero mean).  If `tolerance != 0` (the default), the hypotheses are rounded
to `Int`s if possible and `Rational`s if not, using the given tolerance.  If
`tolerance == 0`, then the hypothesis matrix is returned as-is.

The orientation of the hypothesis matrix is _opposite_ that of the contrast
matrix: each row of the contrasts matrix is a data level and each column is a
predictor, whereas each row of the hypothesis matrix is the interpretation of a
predictor with weights for each level given in the columns.

Note that this assumes a *balanced design* where there are the same number of
observations in every cell.  This is only important for non-orthgonal contrasts
(including contrasts that are not orthogonal with the intercept).

# Examples

```jldoctest hypmat
julia> cmat = StatsModels.contrasts_matrix(DummyCoding(), 1, 4)
4×3 Array{Float64,2}:
 0.0  0.0  0.0
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

julia> StatsModels.hypothesis_matrix(cmat)
4×4 Array{Int64,2}:
  1  0  0  0
 -1  1  0  0
 -1  0  1  0
 -1  0  0  1
```

For non-centered contrasts like `DummyCoding`, without including the intercept
the hypothesis matrix is incorrect.  So while `intercept=true` is the default for
non-centered contrasts, you can see the (wrong) hypothesis matrix when ignoring
it by forcing `intercept=false`:

```jldoctest hypmat
julia> StatsModels.hypothesis_matrix(cmat, intercept=false)
3×4 Array{Int64,2}:
 0  1  0  0
 0  0  1  0
 0  0  0  1
```

The default behavior is to coerce to the nearest integer or rational value, with
a tolerance of the `tolerance` kwarg (defaults to `1e-5`).  The raw
pseudo-inverse matrix can be obtained as `Float64` by setting `tolerance=0`:

```julia-repl
julia> StatsModels.hypothesis_matrix(cmat, tolerance=0) # ugly
4×4 Array{Float64,2}:
  1.0  -2.23753e-16   6.91749e-18  -1.31485e-16
 -1.0   1.0          -2.42066e-16   9.93754e-17
 -1.0   4.94472e-17   1.0           9.93754e-17
 -1.0   1.04958e-16  -1.31044e-16   1.0
```

Finally, the hypothesis matrix for a constructed `ContrastsMatrix` (as stored by
`CategoricalTerm`s) can also be extracted:

```jldoctest hypmat
julia> StatsModels.hypothesis_matrix(StatsModels.ContrastsMatrix(DummyCoding(), ["a", "b", "c", "d"]))
4×4 Array{Int64,2}:
  1  0  0  0
 -1  1  0  0
 -1  0  1  0
 -1  0  0  1

```
"""
function hypothesis_matrix(cm::AbstractMatrix;
                           intercept::Bool=needs_intercept(cm), tolerance::Real=1e-5)
    if intercept
        cm = hcat(ones(size(cm, 1)), cm)
    end
    hypotheses = pinv(cm)
    iszero(tolerance) ? hypotheses : pretty_mat(hypotheses, tol=tolerance)
end

hypothesis_matrix(contrasts::AbstractContrasts, n::Integer; baseind::Integer=1, kwargs...) =
    hypothesis_matrix(contrasts_matrix(contrasts, baseind, n); kwargs...)

hypothesis_matrix(cmat::ContrastsMatrix; kwargs...) =
    hypothesis_matrix(cmat.matrix; kwargs...)

function pretty_mat(mat::AbstractMatrix; tol::Real=10*eps(eltype(mat)))
    fracs = rationalize.(mat, tol=tol)
    if all(x -> denominator(x) == 1, fracs)
        return Int.(fracs)
    else
        return fracs
    end
end


text:
module Alchemy

using GLMakie
using GLMakie: GLFW, to_native
using Makie
using Molecules

include("Constants.jl")
include("Molecules.jl")

end # module

ext: jl


text:
# test 1: minimal options/arguments, auto-generated help/version

using ArgParse
using Base.Test

function ap_test1(args)

    s = ArgParseSettings()

    @add_arg_table s begin
        "--opt1"               # an option (will take an argument)
        "--opt2", "-o"         # another option, with short form
        "arg1"                 # a positional argument
    end

    s.exc_handler = (settings, err)->error(err.text)

    parsed_args = parse_args(args, s) # the result is a Dict{String,Any}
end

@test ap_test1([]) == (String=>Any)["opt1"=>nothing, "opt2"=>nothing, "arg1"=>nothing]
@test ap_test1(["arg"]) == (String=>Any)["opt1"=>nothing, "opt2"=>nothing, "arg1"=>"arg"]
@test ap_test1(["--opt1", "X", "-o=5", "--", "-arg"]) == (String=>Any)["opt1"=>"X", "opt2"=>"5", "arg1"=>"-arg"]
@test_throws ap_test1(["--opt1", "X", "-o=5", "-arg"])
@test ap_test1(["--opt1=", "--opt2=5"]) == (String=>Any)["opt1"=>"", "opt2"=>"5", "arg1"=>nothing]


text:
export GaussianKernel

"Gaussian smoothing kernel."
struct GaussianKernel{N} <: KernelFormula{N}
    support::Float
    σ::Float

    GaussianKernel{1}(; support::Number = 3.0) =
        @check(support > 0) && new(support, 1 / √π)
    GaussianKernel{2}(; support::Number = 3.0) =
        @check(support > 0) && new(support, 1 / π)
    GaussianKernel{3}(; support::Number = 3.0) =
        @check(support > 0) && new(support, 1 / (π * √π))
end

Base.show(io::IO, ::Type{GaussianKernel{N}}) where {N} =
    print(io, "GaussianKernel")

"""
    (W::GaussianKernel)(q::Number) -> Float

Evaluate the kernel with unit width at relative distance `q`.
"""
(W::GaussianKernel)(q::Number) = @dbgassert(q ≥ 0) && W.σ * exp(-q^2)

"""
    ddr(W::GaussianKernel, q::Number) -> Float

Evaluate the derivative of the kernel with unit width at relative distance `q`.
The derivative is with respect to `r`.
"""
ddr(W::GaussianKernel, q::Number) = -2q * W(q)

"""
    ddh(W::GaussianKernel, q::Number) -> Float

Evaluate the derivative of the kernel with unit width at relative distance `q`.
The derivative is with respect to `h`.
"""
ddh(W::GaussianKernel{N}, q::Number) where {N} = (2 * q^2 - N) * W(q)

get_support(W::GaussianKernel) = W.support

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
# implementation of a sorted dict (not optimized for speed) for storing
# the factorization of an integer

struct Factorization{T<:Integer} <: AbstractDict{T, Int}
    pe::Vector{Pair{T, Int}} # Prime-Exponent

    Factorization{T}() where {T<:Integer} = new{T}(Vector{Pair{T, Int}}())
end

function Factorization{T}(d::AbstractDict) where T<:Integer
    f = Factorization{T}()
    append!(f.pe, sort!(collect(d)))
    f
end

Factorization(d::AbstractDict{T}) where {T<:Integer} = Factorization{T}(d)
Base.convert(::Type{Factorization}, d::AbstractDict) = Factorization(d)

Base.iterate(f::Factorization, state...) = iterate(f.pe, state...)

function Base.get(f::Factorization, p, default)
    found = searchsorted(f.pe, p, by=first)
    isempty(found) ?
        default :
        last(f.pe[first(found)])
end

Base.getindex(f::Factorization, p::Integer) = get(f, p, 0)

function Base.setindex!(f::Factorization{T}, e::Int, p::Integer) where T
    found = searchsorted(f.pe, p, by=first)
    if isempty(found)
        insert!(f.pe, first(found), T(p)=>e)
    else
        f.pe[first(found)] = T(p)=>e
    end
    f
end

Base.length(f::Factorization) = length(f.pe)

Base.show(io::IO, ::MIME{Symbol("text/plain")}, f::Factorization) =
    join(io, isempty(f) ? "1" : [(e == 1 ? "$p" : "$p^$e") for (p,e) in f.pe], " * ")


text:
#OLS QR-Decomposition function using in GlobalSearchRegression.jl

function qrreg(data)
    qrf = qr(expvars_data)
    b = qrf \ depvar_data                   # estimate
    ŷ = expvars_data * b                    # predicted values
    er = depvar_data - ŷ                    # in-sample residuals
    er2 = er .^ 2                           # squared errors
    sse = sum(er2)                          # residual sum of squares
    df_e = nobs - ncoef                     # degrees of freedom
    rmse = sqrt(sse / nobs)                 # root mean squared error
    r2 = 1 - var(er) / var(depvar_data)     # model R-squared
    bstd = sqrt.( sum( (UpperTriangular(qrf.R) \ Matrix(1.0LinearAlgebra.I, ncoef, ncoef) ) .^ 2, dims=2) * (sse / df_e) ) # std deviation of coefficients

end


text:
function WeightedObservedPoint(arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return WeightedObservedPoint((jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function get_weight(obj::WeightedObservedPoint)
    return jcall(obj, "getWeight", jdouble, ())
end

function get_x(obj::WeightedObservedPoint)
    return jcall(obj, "getX", jdouble, ())
end

function get_y(obj::WeightedObservedPoint)
    return jcall(obj, "getY", jdouble, ())
end


text:
include("setup_1d.jl")

println("dwt by filtering (N=",N,"), ", L, " levels")
f(x0, L, wf) = for i = 1:tn; dwt(x0, wf, L); end
f(x0, L, wf);
@time f(x0, L, wf);


text:
using DiffEqBase: CallbackSet

using NQCDynamics: InitialConditions, DynamicsOutputs
using NQCDistributions: DynamicalDistribution, ProductDistribution

abstract type AbstractSelection end

"""
Select the initial conditions from the distribution in order.
"""
struct OrderedSelection{D,I} <: AbstractSelection
    "Distribution that is sampled."
    distribution::D
    indices::I
end

"""
Obtain initial conditions by randomly sampling the distribution.
"""
struct RandomSelection{D} <: AbstractSelection
    "Distribution that is sampled."
    distribution::D
end

function Selection(distribution, selection::AbstractVector)
    @info "Sampling the provided distribution in the range $selection."
    OrderedSelection(distribution, selection)
end

function Selection(distribution, ::Any)
    @info "Sampling randomly from provided distribution."
    RandomSelection(distribution)
end

function (select::OrderedSelection)(prob, i, repeat)
    j = select.indices[i]
    u0 = sample_distribution(prob.p, select.distribution, j)
    DynamicsMethods.create_problem(u0, prob.tspan, prob.p)
end

function sample_distribution(sim::AbstractSimulation, distribution::DynamicalDistribution, i)
    u = distribution[i]
    DynamicsMethods.DynamicsVariables(sim, u.v, u.r)
end

function sample_distribution(sim::AbstractSimulation, distribution::ProductDistribution, i)
    u = distribution.nuclear[i]
    DynamicsMethods.DynamicsVariables(sim, u.v, u.r, distribution.electronic)
end

function (select::RandomSelection)(prob, i, repeat)
    u0 = sample_distribution(prob.p, select.distribution)
    DynamicsMethods.create_problem(u0, prob.tspan, prob.p)
end

function sample_distribution(sim::AbstractSimulation, distribution::DynamicalDistribution)
    u = rand(distribution)
    DynamicsMethods.DynamicsVariables(sim, u.v, u.r)
end

function sample_distribution(sim::AbstractSimulation, distribution::ProductDistribution)
    u = rand(distribution.nuclear)
    DynamicsMethods.DynamicsVariables(sim, u.v, u.r, distribution.electronic)
end


text:
function eval_TIFₜ(xₜ::Vector{Float64}, coreₜ::MaximinObjectiveCore) # x = I
    return minimum(LinearAlgebra.BLAS.gemv('N', coreₜ.IIF, xₜ))::Float64
end

function eval_TIFₜ(xₜ::Vector{Float64}, coreₜ::CCMaximinObjectiveCore) # x = I
    IIF = coreₜ.IIF
    α = coreₜ.alpha
    K, I, R = size(IIF)
    alphaR = Int(ceil(R * (α)))
    αQle = Inf
    if K > 1
        for k = 1:K
            αQle = min(αQle, sort(LinearAlgebra.BLAS.gemv('T', IIF[k, :, :], xₜ))[alphaR])
        end
    else
        αQle = sort(_gemvblasT(IIF[1, :, :], xₜ, R))[alphaR]
    end
    return αQle::Float64
end

function eval_TIFₜ(xₜ::Vector{Float64}, coreₜ::MinimaxObjectiveCore) # x = I
    IIF = coreₜ.IIF
    targets = coreₜ.targets
    return -maximum(abs.(LinearAlgebra.BLAS.gemv('N', IIF, xₜ) - targets))::Float64
end


text:
"""
    struct Equation{T, TL, G<:AbstractFloat}

The equation to be solved `∂u/∂t = L*u + N(u)`. Array `L` includes the coefficients
of the linear term `L*u` and `calcN!` is a function which computes the nonlinear
term `N(u)`. The struct also includes the problem's `grid` and the float type of the
state vector (and consequently of `N(u)`).

$(TYPEDFIELDS)
"""
struct Equation{T, TL, G<:AbstractFloat}
    "array with coefficient for the linear part of the equation"
       L :: TL
    "function that computes the nonlinear part of the equation"
  calcN! :: Function
    "the grid"
    grid :: AbstractGrid{G}
    "the dimensions of `L`"
    dims :: Tuple
    "the float type for the state vector"
       T :: T # eltype or tuple of eltypes of sol and N
end

"""
    Equation(L, calcN!, grid; dims=supersize(L), T=nothing)

The equation constructor from the array `L` of the coefficients of the linear term, the function
`calcN!` that computes the nonlinear term and the `grid` for the problem.
"""
function Equation(L, calcN!, grid::AbstractGrid{G}; dims=supersize(L), T=nothing) where G
  T = T == nothing ? T = cxtype(G) : T

  return Equation(L, calcN!, grid, dims, T)
end

"""
    mutable struct Clock{T<:AbstractFloat}

Represents the clock of a problem.

$(TYPEDFIELDS)
"""
mutable struct Clock{T<:AbstractFloat}
    "the time-step"
    dt :: T
    "the time"
     t :: T
    "the step number"
  step :: Int
end

"""
    struct Problem{T, A<:AbstractArray, Tg<:AbstractFloat, TL}

A problem that represents a partial differential equation.

$(TYPEDFIELDS)
"""
struct Problem{T, A<:AbstractArray, Tg<:AbstractFloat, TL}
    "the state vector"
          sol :: A
    "the problem's"
        clock :: Clock{Tg}
    "the equation"
          eqn :: Equation{T, TL, Tg}
    "the grid"
         grid :: AbstractGrid{Tg}
    "the variables"
         vars :: AbstractVars
    "the parameters"
       params :: AbstractParams
    "the timestepper"
  timestepper :: AbstractTimeStepper{A}
end

"""
    EmptyParams <: AbstractParams

A placeholder struct for parameters.
"""
struct EmptyParams <: AbstractParams end

"""
    EmptyVars <: AbstractVars

A placeholder struct for variables.
"""
struct EmptyVars <: AbstractVars end

"""
    Problem(eqn::Equation, stepper, dt, grid::AbstractGrid{T},
            vars=EmptyVars, params=EmptyParams, dev::Device=CPU(); stepperkwargs...) where T

Construct a `Problem` for equation `eqn` using the time`stepper` with timestep
`dt`, on `grid` and on `dev`ice. Optionally, use the keyword arguments to provide
variables with `vars` and parameters with `params`. The `stepperkwargs` are passed
to the time-stepper constructor.
"""
function Problem(eqn::Equation, stepper, dt, grid::AbstractGrid{T},
                 vars=EmptyVars, params=EmptyParams, dev::Device=CPU(); stepperkwargs...) where T
  clock = Clock{T}(dt, 0, 0)

  timestepper = TimeStepper(stepper, eqn, dt, dev; stepperkwargs...)

  sol = devzeros(dev, eqn.T, eqn.dims)

  return Problem(sol, clock, eqn, grid, vars, params, timestepper)
end

show(io::IO, clock::FourierFlows.Clock) =
     print(io, "Clock\n",
               "  ├─── timestep dt: ", clock.dt, '\n',
               "  ├────────── step: ", clock.step, '\n',
               "  └──────── time t: ", clock.t)

show(io::IO, eqn::FourierFlows.Equation) =
     print(io, "Equation\n",
               "  ├──────── linear coefficients: L", '\n',
               "  │                              ├───type: ", eltype(eqn.L), '\n',
               "  │                              └───size: ", size(eqn.L), '\n',
               "  ├───────────── nonlinear term: calcN!()", '\n',
               "  └─── type of state vector sol: ", eqn.T)

show(io::IO, problem::FourierFlows.Problem) =
    print(io, "Problem\n",
              "  ├─────────── grid: grid (on " * FourierFlows.griddevice(problem.grid) * ")", '\n',
              "  ├───── parameters: params", '\n',
              "  ├────── variables: vars", '\n',
              "  ├─── state vector: sol", '\n',
              "  ├─────── equation: eqn", '\n',
              "  ├────────── clock: clock", '\n',
              "  └──── timestepper: ", string(nameof(typeof(problem.timestepper))))


text:
# Test Parameter Learning routines

# @testset "SPN structure learning: LearnSPN" begin

#     println()

#     @info " Generating true model..."
#     # First create SPN to generate data from
#     A = Vector([
#         SumNode(),     # 1
#         ProductNode(), # 2
#         ProductNode(), # 3
#         ProductNode(), # 4
#         CategoricalDistribution(1, [0.6,0.4]), # 5
#         CategoricalDistribution(1, [0.1,0.9]), # 6
#         CategoricalDistribution(2, [0.3,0.7]), # 7
#         CategoricalDistribution(2, [0.8,0.2])  # 8
#     ])
#     B = [
#         [2, 3, 4], # 1 Sum P(A, B) 0.2,0.5,0.3
#         [5, 7],    # 2 Prod P(A, B)=P(A)*P(B)
#         [5, 8],    # 3 Prod P(A, B)=P(A)*P(B)
#         [6, 8],    # 4 Prod P(A, B)=P(A)*P(B)
#     ]
#     C = sparse([2,   3,     4],
#                [1,   1,     1],
#                [0.2, 0.5, 0.3])

#     SPN = SumProductNetwork{Int64}(A, B, C)
#     @info " Generating data set..."
#     N = 3000
#     println("N = $N")
#     data = rand(SPN,N)

#     @info " Learning SPN structure from data..."
#     @time SPN2 = learnspn(data,[2,2])
#     println(SPN2)

#     # @info "## Printing graph structure"
#     # io = IOBuffer()
#     # todot(io,SPN2)
#     # open("spn.dot", "w") do f
#     #     write(f, String(take!(io)))
#     # end
#     # println("Call 'dot -Tpdf spn.dot -o spn.pdf' to generate image")

#     #@test MAE(SPN,SPN2,data) < 0.05

#     println()
#     Z = 0.0
#     for a in 1:2, b in 1:2
#         ref = SPN([a,b])
#         est = SPN2([a,b])
#         println("S($a,$b) = $ref ≈ ", est)
#         Z += est
#     end
#     @test Z ≈ 1.0
#     println()

# end

@testset "LearnSPN: NLTCS" begin

    using DataFrames, CSV

    @info "Reading and processing dataset..."
    # Load data set
    @time dataset = CSV.read("../data/nltcs.train.csv",
                             header=collect(1:16) # columns names
    );

    println("Dataset has $(size(dataset,1)) instances, $(size(dataset,2)) columns.")

    data = convert(Matrix, dataset) .+ 1
    # @info "done"

    # build Dense SPN
    @info "Learning structure..."
    @time S = learnspn(data,repeat([2],16,),100)
    @info "done." S
    #save(S,"../examples/nltcs.learnspn.notuning.spn")

    @info "Evaluating learned model..."
    @time testset = CSV.read("../data/nltcs.test.csv",
                             header=collect(1:16) # columns names
    );
    println("Testset has $(size(testset,1)) instances, $(size(testset,2)) columns.")
    test = convert(Matrix, testset) .+ 1
    println("Avg Training Negative LL: ", NLL(S,data))
    println("Avg Testset Negative LL: ", NLL(S,test))

    @info "Tuning parameters..."
    learner = EMParamLearner()
    learner.tolerance = 1e-6
    #println("Initializing weights...")
    #initialize(S) # randomly initialize weights
    # Now learn network parameters from data
    while !converged(learner) && learner.steps < 500
        SumProductNetworks.step(learner,S,data)
        println("It: $(learner.steps) \t ANLL: $(-learner.score) \t ANLL: $(NLL(S,test))")
        if !isfinite(learner.score)
            break
        end
    end
    println("Avg Training Negative LL: ", NLL(S,data))
    println("Avg Testset Negative LL: ", NLL(S,test))

    #save(S,"../examples/nltcs.learnspn.em.spn")

end

# @testset "LearnSPN: ZIP dataset" begin

#     using DataFrames, CSV

#     width,height = 16,16

#     @info "Reading and processing dataset..."
#     # Load data set
#     @time dataset = CSV.read("../data/zip.test",
#          header=["digit"; collect(map(string, 1:width*height))], # columns names
#          copycols=true # so that we can manipulate columns
#     );
#     nine = dataset[dataset.digit .== 9, 2:257];
#     # count instances of digit 9 in training set
#     println("Dataset has $(size(nine,1)) instances, $(size(nine,2)) columns.")
#     # binarize data
#     bin_nine = copy(nine)
#     for i=1:(width*height)
#         bin_nine[nine[i] .<= -0.5, i] .= 1
#         idx = (nine[i] .> -0.5) .& (nine[i] .<= 0.5)
#         bin_nine[idx, i] .= 2
#         bin_nine[nine[i] .> 0.5, i] .= 3
#     end
#     # threshold = 0.0
#     # for i=1:256
#     #     bin_nine[nine[i] .<= threshold, i] .= 1
#     #     bin_nine[nine[i] .> threshold, i] .= 2
#     # end
#     #println(first(bin_nine,5))
#     data = convert(Matrix, bin_nine);
#     # @info "done"

#     # build Dense SPN
#     @info "Learning structure..."
#     @time S = learnspn(data,repeat([3],width*height))
#     @info "done." S
#     print("Cardinalities: ")
#     for node in leaves(S)
#         print("$(length(node.values)) ")
#     end
#     println()
#     computed = falses(length(S))
#     for i in Iterators.Reverse(1:length(S))
#         computed[i] = true
#         if isa(S[i],SumNode) || isa(S[i],ProductNode)
#             if i > size(S._backward,1)
#                 @error "Inner node is not in backward" i size(S._backward,1)
#             end
#             for j in children(S,i)
#                 if computed[j] == false
#                     @error "Parent computed before child" i j Type{S[i]} Type{S[j]}
#                     error("Invalid spn")
#                 end
#             end
#         end
#     end

#     # @info "Printing graph structure to file..."
#     # io = IOBuffer()
#     # todot(io,S)
#     # open("../examples/densespn.dot", "w") do f
#     #     write(f, String(take!(io)))
#     # end
#     # println("done! Call dot -Tpdf examples/densespn.dot -o examples/densespn.pdf to generate image")

#     @info "Learning parameters..."
#     learner = EMParamLearner()
#     learner.tolerance = 1e-3
#     learner.minimumvariance = 0.1
#     #println("Initializing weights...")
#     #@time initialize(S) # randomly initialize weights
#     # Now learn network parameters from data
#     while !converged(learner) && learner.steps < 10
#         SumProductNetworks.step(learner,S,data)
#         println("It: $(learner.steps) \t NLL: $(-learner.score)")
#         if !isfinite(learner.score)
#             break
#         end
#     end
#     # @info "done"

#     # run inference
#     @info "Testing if log evaluation is consistent..."
#     @time logpdf(S,view(data,2,:))
#     @test logpdf(S,data[1,:]) == logpdf(S,view(data,1,:))
#     # @info "done."

#     @info "Sampling from network..."
#     @time x = rand(S)
#     # @info "done."

# end # END of test set


text:
struct ∇SemFIML{
        X <: AbstractArray,
        Y <: AbstractArray,
        Z <: AbstractArray,
        W,
        F,
        F2,
        I,
        I2,
        T,
        F3,
        TM,
        F4,
        TBm,
        I3,
        INV <: AbstractArray,
        C <: AbstractArray,
        L <: AbstractArray,
        M <: AbstractArray,
        IM <: AbstractArray,
        I4 <: AbstractArray,
        C2} <: DiffFunction
    B::X
    B!::F
    E::Y
    E!::F2
    F::Z
    C::W
    S_ind_vec::I
    A_ind_vec::I2
    matsize::T
    M!::F3
    M::TM
    Bm!::F4
    Bm::TBm
    M_ind_vec::I3
    inverses::INV #preallocated inverses of imp_cov
    choleskys::C #preallocated choleskys
    logdets::L #logdets of implied covmats
    imp_mean::IM
    meandiff::M
    imp_inv::I4
    counter::C2
end

function ∇SemFIML(
    observed::O where {O <: SemObs},
    imply::I where {I <: Imply},
    A::Spa1,
    S::Spa2,
    F::Spa3,
    parameters,
    start_val;
    M = nothing
        ) where {
        Spa1 <: SparseMatrixCSC,
        Spa2 <: SparseMatrixCSC,
        Spa3 <: SparseMatrixCSC
        }

    A = copy(A)
    S = copy(S)
    F = copy(F)

    invia = neumann_series(A)

    #imp_cov_sym = F*invia*S*permutedims(invia)*permutedims(F)
    #imp_cov_sym = Array(imp_cov_sym)
    invia = ModelingToolkit.simplify.(invia)
    B = invia
    E = B*S*B'
    E = E*permutedims(F)
    E = Array(E)
    E = ModelingToolkit.simplify.(E)
    B = F*B
    B = Array(B)
    B = ModelingToolkit.simplify.(B)

    B! =
        eval(ModelingToolkit.build_function(
            B,
            parameters
            )[2])

    E! =
        eval(ModelingToolkit.build_function(
            E,
            parameters
            )[2])

    B_pre = zeros(size(F)...)
    E_pre = zeros(size(F, 2), size(F, 1))

    grad = similar(start_val)
    matsize = size(A)
    C_pre = zeros(size(F, 1), size(F, 1))

    #S_ind_vec = Vector{Tuple{Array{Int64,1},Array{Int64,1},Array{Float64,1}}}()
    #A_ind_vec = Vector{Tuple{Array{Int64,1},Array{Int64,1},Array{Float64,1}}}()
    S_ind_vec = Vector{Vector{CartesianIndex{2}}}()
    A_ind_vec = Vector{Vector{CartesianIndex{2}}}()

    for i = 1:size(parameters, 1)
        S_ind = findall(var -> isequal(var, parameters[i]), S)
        A_ind = findall(var -> isequal(var, parameters[i]), A)

        push!(S_ind_vec, S_ind)
        push!(A_ind_vec, A_ind)
    end

    if !isnothing(M)

        fun_mean =
            eval(ModelingToolkit.build_function(
                M,
                parameters
            )[2])
        M_pre = zeros(size(F)[2])

        Bm = B*M
        Bm = ModelingToolkit.simplify.(Bm)
        fun_bm =
            eval(ModelingToolkit.build_function(
                Bm,
                parameters
            )[2])
        Bm_pre = zeros(size(F)[2])

        M_ind_vec = Vector{Vector{Int64}}()
        for i = 1:size(parameters, 1)
            M_ind = findall(var -> isequal(var, parameters[i]), M)
            push!(M_ind_vec, M_ind)
        end
    else
        M_ind_vec = nothing
        M_pre = nothing
        Bm_pre = nothing
        fun_mean = nothing
        fun_bm = nothing
    end

    inverses = broadcast(x -> zeros(x, x), Int64.(observed.pattern_nvar_obs))
    choleskys = Array{Cholesky{Float64,Array{Float64,2}},1}(undef, length(inverses))

    n_patterns = size(observed.rows, 1)
    logdets = zeros(n_patterns)

    imp_mean = zeros.(Int64.(observed.pattern_nvar_obs))
    meandiff = zeros.(Int64.(observed.pattern_nvar_obs))
    imp_inv = zeros(size(observed.data, 2), size(observed.data, 2))

    return ∇SemFIML(
        B_pre,
        B!,
        E_pre,
        E!,
        F,
        C_pre,
        S_ind_vec,
        A_ind_vec,
        matsize,
        fun_mean,
        M_pre,
        fun_bm,
        Bm_pre,
        M_ind_vec,
        inverses,
        choleskys,
        logdets,
        imp_mean,
        meandiff,
        imp_inv,
        [0.0])
end

function compute_C(inv_cov, obs_cov, N)
    if N > one(N)
        C = LinearAlgebra.I - inv_cov*obs_cov
        return C
    end
    return reshape([0.0],1,1)
end

function grad_per_pattern(pattern, Σ_inv_j, Σ_der, µ_der, N, C, b)
    Σ_der_j = Σ_der[pattern, pattern]
    µ_der_j = µ_der[pattern]
    if N > one(N)
        d1 = dot(Σ_inv_j,Σ_der_j,C)
    else
        d1 = dot(Σ_inv_j,Σ_der_j)
    end
    d2 = (b'*Σ_inv_j*Σ_der_j + 2*µ_der_j')*Σ_inv_j*b
    d = N*(d1 - d2)
    return d
end

function (diff::∇SemFIML)(par, grad, model::Sem{O, I, L, D}) where
            {O <: SemObs, L <: Loss, I <: Imply, D <: SemAnalyticDiff}
    #ld = logdet(a)
    if check_fiml_2(diff, model)
        copy_per_pattern!(diff, model)
        #batch_cholesky!(diff, model)
        #diff.logdets .= logdet.(diff.choleskys)
        #batch_inv!(diff, model)
        copy_per_pattern!(diff, model)
        batch_sym_inv_update!(diff, model)
        diff.logdets .= -logdet.(diff.inverses)

        diff.B!(diff.B, par) # B = inv(I-A)
        diff.E!(diff.E, par) # E = B*S*B'
        for i in 1:size(diff.meandiff, 1)
            diff.meandiff[i] .=
                model.observed.obs_mean[i] -
                model.imply.imp_mean[model.observed.patterns[i]]
        end

        C = [compute_C(diff.inverses[i], model.observed.obs_cov[i], model.observed.pattern_n_obs[i])
                for i in 1:size(diff.meandiff, 1)]

        let B = diff.B, E = diff.E, F = diff.F, Bm = diff.Bm
            Threads.@threads for i = 1:size(par, 1)
                term = similar(diff.C)
                term2 = similar(diff.C)
                sparse_outer_mul!(term, B, E, diff.A_ind_vec[i])
                sparse_outer_mul!(term2, B, B', diff.S_ind_vec[i])
                Σ_der = term2 + term + term'

                term3 = Vector{Float64}(undef, size(B, 1))
                term4 = similar(term3)
                sparse_outer_mul!(term3, B, Bm, diff.A_ind_vec[i])
                sparse_outer_mul!(term4, B, diff.M_ind_vec[i])
                µ_der = term3 + term4

                for j = 1:size(model.observed.patterns, 1)
                    grad[i] +=
                        grad_per_pattern(
                            model.observed.patterns[j],
                            diff.inverses[j],
                            Σ_der,
                            µ_der,
                            model.observed.pattern_n_obs[j],
                            C[j],
                            diff.meandiff[j])
                end
            end
        end
    else
        grad .= 0
    end
end


text:
using Base:Float64

# concrete result type -
struct VLResult{T}
    value::T
end

abstract type VLAbstractAsset end
abstract type VLAbstractLattice end

mutable struct VLEquityAsset <: VLAbstractAsset

    # data -
    assetSymbol::String
    purchasePricePerShare::Float64
    numberOfShares::Int64
    purchaseDate::Date

    function VLEquityAsset(assetSymbol::String, purchasePricePerShare::Float64, numberOfShares::Int64, purchaseDate::Date)
        this = new(assetSymbol, purchasePricePerShare, numberOfShares, purchaseDate)
    end
end

mutable struct VLCallOptionContract <: VLAbstractAsset

    # data -
    assetSymbol::String
    strikePrice::Float64
    expirationDate::Date
    premimumValue::Float64
    numberOfContracts::Int64
    sense::Symbol
    contractMultiplier::Float64

    function VLCallOptionContract(assetSymbol::String, expirationDate::Union{Nothing,Date}, strikePrice::Float64, premimumValue::Float64,
        numberOfContracts::Int64; sense::Symbol=:buy, contractMultiplier::Float64=100.0)

        this = new(assetSymbol, strikePrice, expirationDate, premimumValue, numberOfContracts, sense, contractMultiplier)
    end
end

mutable struct VLPutOptionContract <: VLAbstractAsset

    # data -
    assetSymbol::String
    strikePrice::Float64
    expirationDate::Union{Nothing,Date}
    premimumValue::Float64
    numberOfContracts::Int64
    sense::Symbol
    contractMultiplier::Float64

    function VLPutOptionContract(assetSymbol::String, strikePrice::Float64, premimumValue::Float64;
        expirationDate::Union{Nothing,Date}=nothing, numberOfContracts::Int64=1, sense::Symbol=:buy, contractMultiplier::Float64=100.0)

        _ = new(assetSymbol, strikePrice, expirationDate, premimumValue, numberOfContracts, sense, contractMultiplier)
    end
end

mutable struct VLBinomialLattice <: VLAbstractLattice

    # data -
    μ::Float64
    σ::Float64
    𝝙t::Float64
    numberOfLevels::Int64

    function VLBinomialLattice(μ::Float64, σ::Float64, 𝝙t::Float64, numberOfLevels::Int64)
        _ = new(μ, σ, 𝝙t, numberOfLevels)
    end
end


text:
#=

Want to process `@constraint f(f(x)) ∈ [0.3, 0.4]`
where `f(x) = 4x * (1-x)`

Given code `f(f(x))`, we need `f_forward` and `f_backward`.

Each "copy" of `f` uses the same actual forward and back functions,
`f.forward` and `f.backward`.

```
@function f(x) = 4x * (1-x)
```
should generate these forward and backward functions, and register the function
`f`.

"""
=#

type FunctionObject
    input_args::Vector{Symbol}
    output_args::Vector{Symbol}
    local_vars::Vector{Symbol}
    code::Expr
end

doc"""
A `ConstraintFunction` contains the created forward and backward
code
"""
type ConstraintFunction{F <: Function, G <: Function}
    input::Vector{Symbol}  # input arguments for forward function
    output::Vector{Symbol} # output arguments for forward function
    forward::F
    backward::G
end

type FunctionArguments
    # input::Vector{Symbol}  # input arguments for forward function
    # output::Vector{Symbol}  # output arguments
    # generated::Vector{Symbol} # local variables generated

    input
    generated
    return_arguments
end

const registered_functions = Dict{Symbol, FunctionArguments}()

@doc """
`@function` registers a function to be used in forwards and backwards mode.

Example: `@function f(x, y) = x^2 + y^2`
"""  # this docstring does not work!

@eval macro ($(:function))(ex)   # workaround to define macro @function

    # (f, args, code) = @match ex begin
    #     ( f_(args__) = code_ ) => (f, args, code)
    # end

    (f, args, code) = match_function(ex)

    @show f, args, code

    #root, all_vars, generated, code2 = IntervalConstraintProgramming.flatten!(code)
    return_arguments, flatAST = IntervalConstraintProgramming.flatten!(code)

    @show return_arguments

    #@show root, all_vars, generated, code2

    # rearrange so actual return arguments come first:

    @show flatAST
    @show return_arguments

    # make into an array:
    if !(isa(return_arguments, Array))
        return_arguments = [return_arguments]
    end

    flatAST.intermediate = setdiff(flatAST.intermediate, return_arguments)

    println("HERE")
    flatAST.intermediate = [return_arguments; flatAST.intermediate]

    forward_code = forward_pass(flatAST) #root, all_vars, generated, code2)
    backward_code = backward_pass(flatAST) #root, all_vars, generated, code2)

    @show forward_code, backward_code

    registered_functions[f] = FunctionArguments(flatAST.variables, flatAST.intermediate, return_arguments)

    return quote
        #$(esc(Meta.quot(f))) = ConstraintFunction($(all_vars), $(generated), $(forward_code), $(backward_code))
        $(esc(f)) = ConstraintFunction($(flatAST.variables), $(flatAST.intermediate), $(forward_code), $(backward_code))
        #registered_functions[$(Meta.quot(f))] =  ConstraintFunction($(all_vars), $(generated), $(forward_code), $(backward_code))
        #$(Meta.quot(f)) =  ConstraintFunction($(all_vars), $(generated), $(forward_code), $(backward_code))
    end
end

function match_function(ex)

    try
        (f, args, body) =
            @match ex begin
             ( (f_(args__) = body_) |
              (function f_(args__) body_ end) ) => (f, args, body)
           end

         return (f, args, rmlines(body))  # rmlines is from MacroTools package

    catch
        throw(ArgumentError("$ex does not have the form of a function"))
    end
end

doc"""
Generate code for an anonymous function with given
input arguments, output arguments, and code block.
"""
function make_function(input_args, output_args, code)

    input = Expr(:tuple, input_args...)  # make a tuple of the variables
    output = Expr(:tuple, output_args...)  # make a tuple of the variables

    new_code = copy(code)
    push!(new_code.args, :(return $output))

    complete_code = :( $input -> $new_code )

    return GeneratedFunction(input_args, output_args, complete_code)
end


text:
facts("Order object methods") do
  context("Constructor checks") do
    @fact_throws TradingLogic.Order("error", 10, 150.0,
                                    :obscuresym,
                                    :market, :complete, "string")
    @fact_throws TradingLogic.Order("error", 10, 150.0,
                                    :buy, :obscuresym,
                                    :complete, "string")
    @fact_throws TradingLogic.Order("error", 10, 150.0,
                                    :buy, :market,
                                    :obscuresym, "string")
  end
  context("Get info and modify methods") do
    emptord = TradingLogic.emptyorder()
    @fact TradingLogic.getorderposchg(emptord) --> 0
    @fact TradingLogic.ispending(emptord) --> false
    emptord.status = :pending
    @fact TradingLogic.ispending(emptord) --> true
    TradingLogic.setcancelled!(emptord)
    @fact TradingLogic.ispending(emptord) --> false

    # targeted position change with proper sign
    emptord.side = :qq
    @fact_throws TradingLogic.getorderposchg(emptord)
    posord = TradingLogic.Order("pos", 10, 150.0,
                                :buy, :limit, :pending, "aux")
    negord = TradingLogic.Order("neg", 10, 150.0,
                                :sell, :limit, :pending, "aux")
    @fact TradingLogic.getorderposchg(posord) > 0 --> true
    @fact TradingLogic.getorderposchg(negord) < 0 --> true
  end
end

facts("Order handling based on target input") do
  context("Limit and marker order from target input") do
    nostoplim = Array(Float64, 0)
    ord = TradingLogic.emptyorder()
    # should not be calling targ2order! with zero position change
    @fact TradingLogic.targ2order!(ord, (0, nostoplim),
                                   "", 0, true) --> false
    # stoplimit is not handled by targ2order!
    @fact TradingLogic.targ2order!(ord, (10, [1.0, 2.0]),
                                   "", 0, true) --> false

    # when no limit price is given: submit market
    ord.ordertype = :limit
    qty = -10
    @fact TradingLogic.targ2order!(ord, (qty, nostoplim),
                                   "", 0, true) --> true
    @fact ord.ordertype --> :market
    @fact isnan(ord.price) --> true
    @fact TradingLogic.getorderposchg(ord) --> qty
    @fact ord.side --> :sell
    @fact TradingLogic.ispending(ord) --> true

    # limit order for submission
    ord = TradingLogic.emptyorder()
    @fact ord.ordertype --> :market
    @fact ord.side --> :buy
    @fact TradingLogic.ispending(ord) --> false
    qty, prc = -50, 100.0
    @fact TradingLogic.targ2order!(ord, (qty, [prc]),
                                   "", 150, true) --> true
    @fact ord.ordertype --> :limit
    @fact isnan(ord.price) --> false
    @fact ord.price --> roughly(prc)
    @fact TradingLogic.getorderposchg(ord) --> qty
    @fact ord.side --> :sell
    @fact TradingLogic.ispending(ord) --> true
  end

  # orderhandling! input
  tnow() = unix2datetime(time())
  tinit = tnow()
  posactual = [0]
  blotter = TradingLogic.emptyblotter()
  backtest = true

  penter, pexit = 150.0, 110.0
  qty = 10
  shortgain = -(pexit - penter) * qty

  context("Order handling: pending order now complete") do
    # any target (irrelevant at this step -> evaluated next step)
    targ = (0, Array(Float64, 0))

    # market order was pending
    ord = TradingLogic.Order("neg", qty, NaN, :sell,
                             :market, :pending, "aux")
    @fact TradingLogic.ispending(ord) --> true
    @fact TradingLogic.orderhandling!(
      targ, penter, tinit, posactual, ord,
      blotter, backtest) --> (true, 0.0)
    # position updated, no new order generated yet
    # (need to re-evaluate target after position update)
    @fact blotter[tinit] --> (-qty, penter)
    @fact TradingLogic.ispending(ord) --> false
    @fact posactual[1] --> -qty

    # same with limit-order that was pending
    #  and is satisfied at current price
    ord = TradingLogic.Order("pos", qty, 120.0, :buy,
                             :limit, :pending, "aux")
    @fact TradingLogic.ispending(ord) --> true
    @fact TradingLogic.orderhandling!(
      targ, pexit, tnow(), posactual, ord,
      blotter, backtest) --> (true, shortgain)
    # position updated, no new order generated yet
    # (need to re-evaluate target after position update)
    @fact TradingLogic.ispending(ord) --> false
    @fact posactual[1] --> 0
    @fact length(keys(blotter)) --> 2
  end

  context("Order handling: no position change targeted") do
    targ = (0, Array(Float64, 0))

    # nothing was pending
    ord = TradingLogic.emptyorder()
    @fact TradingLogic.ispending(ord) --> false
    @fact TradingLogic.orderhandling!(
      targ, 100.0, tnow(), posactual, ord,
      blotter, backtest) --> (true, shortgain)
    # now still nothing pending, actual position is the same
    @fact TradingLogic.ispending(ord) --> false
    @fact posactual[1] --> 0

    # order was pending that is still pending (limit)
    ord = TradingLogic.Order("pos", 10, 50.0, :buy,
                             :limit, :pending, "aux")
    @fact TradingLogic.ispending(ord) --> true
    @fact TradingLogic.orderhandling!(
      targ, 100.0, tnow(), posactual, ord,
      blotter, backtest) --> (true, shortgain)
    # now pending order gets cancelled
    @fact TradingLogic.ispending(ord) --> false
    @fact ord.status --> :cancelled
    # no position updates
    @fact posactual[1] --> 0
  end

  context("Order handling: market buy/sell") do
    targ = (10, Array(Float64, 0))

    # nothing was pending
    ord = TradingLogic.emptyorder()
    @fact TradingLogic.ispending(ord) --> false
    @fact TradingLogic.orderhandling!(
      targ, 100.0, tnow(), posactual, ord,
      blotter, backtest) --> (true, shortgain)
    # now pending, actual position still the same
    @fact TradingLogic.ispending(ord) --> true
    @fact ord.ordertype --> :market
    @fact posactual[1] --> 0

    # order was pending that is still pending (limit)
    ord = TradingLogic.Order("pos", 10, 50.0, :buy,
                             :limit, :pending, "aux")
    @fact TradingLogic.ispending(ord) --> true
    @fact TradingLogic.orderhandling!(
      targ, 100.0, tnow(), posactual, ord,
      blotter, backtest) --> (true, shortgain)
    # now market target so cancel pending order
    @fact TradingLogic.ispending(ord) --> false
    @fact ord.status --> :cancelled
    # no position updates
    @fact posactual[1] --> 0
  end

  context("Order handling: limit buy/sell") do
    targ = (30, [75.0])

    # nothing was pending
    ord = TradingLogic.emptyorder()
    @fact TradingLogic.ispending(ord) --> false
    @fact TradingLogic.orderhandling!(
      targ, 100.0, tnow(), posactual, ord,
      blotter, backtest) --> (true, shortgain)
    # now limit pending, actual position still the same
    @fact TradingLogic.ispending(ord) --> true
    @fact ord.ordertype --> :limit
    @fact ord.side --> :buy
    @fact posactual[1] --> 0

    # limit order was pending that is still pending
    # in line with the target
    ord = TradingLogic.Order("pos", 30, 75.0, :buy,
                             :limit, :pending, "aux")
    @fact TradingLogic.orderhandling!(
      targ, 100.0, tnow(), posactual, ord,
      blotter, backtest) --> (true, shortgain)
    # leave it pending
    @fact TradingLogic.ispending(ord) --> true
    @fact ord.ordertype --> :limit
    @fact ord.side --> :buy
    @fact posactual[1] --> 0

    # limit order was pending that is still pending
    # position change not in line with the target (value)
    ord = TradingLogic.Order("pos", 25, 75.0, :buy,
                             :limit, :pending, "aux")
    @fact TradingLogic.orderhandling!(
      targ, 100.0, tnow(), posactual, ord,
      blotter, backtest) --> (true, shortgain)
    # cancel different order
    @fact TradingLogic.ispending(ord) --> false
    @fact ord.status --> :cancelled
    # no position updates
    @fact posactual[1] --> 0

    # limit order was pending that is still pending
    # position change not in line with the target (side)
    ord = TradingLogic.Order("pos", 30, 75.0, :sell,
                             :limit, :pending, "aux")
    @fact TradingLogic.orderhandling!(
      targ, 70.0, tnow(), posactual, ord,
      blotter, backtest) --> (true, shortgain)
    # cancel different order
    @fact TradingLogic.ispending(ord) --> false
    @fact ord.status --> :cancelled
    # no position updates
    @fact posactual[1] --> 0

    # limit order was pending that is still pending
    # limit price not in line with the target
    ord = TradingLogic.Order("pos", 30, 85.0, :buy,
                             :limit, :pending, "aux")
    @fact TradingLogic.orderhandling!(
      targ, 100.0, tnow(), posactual, ord,
      blotter, backtest) --> (true, shortgain)
    # cancel different order
    @fact TradingLogic.ispending(ord) --> false
    @fact ord.status --> :cancelled
    # no position updates
    @fact posactual[1] --> 0
  end

  context("Order handling: stoplimit to track, buy-side") do
    targ = (50, [80.0, 75.0]) #[limitprice, stopprice]

    # nothing was pending
    # stop price not reached
    ord = TradingLogic.emptyorder()
    @fact TradingLogic.ispending(ord) --> false
    @fact TradingLogic.orderhandling!(
      targ, 70.0, tnow(), posactual, ord,
      blotter, backtest) --> (true, shortgain)
    # no change, keep tracking the price
    @fact TradingLogic.ispending(ord) --> false
    @fact posactual[1] --> 0

    # nothing was pending
    # stop price reached
    ord = TradingLogic.emptyorder()
    @fact TradingLogic.ispending(ord) --> false
    @fact TradingLogic.orderhandling!(
      targ, 75.1, tnow(), posactual, ord,
      blotter, backtest) --> (true, shortgain)
    # now limit pending, actual position still the same
    @fact TradingLogic.ispending(ord) --> true
    @fact ord.ordertype --> :limit
    @fact ord.side --> :buy
    @fact ord.price --> roughly(80.0)
    @fact posactual[1] --> 0

    # limit order was pending that is still pending
    # stop price not reached
    ord = TradingLogic.Order("pos", 30, 65.0, :buy,
                             :limit, :pending, "aux")
    @fact TradingLogic.orderhandling!(
      targ, 70.0, tnow(), posactual, ord,
      blotter, backtest) --> (true, shortgain)
    # cancel previous order, keep tracking the price
    @fact TradingLogic.ispending(ord) --> false
    @fact ord.status --> :cancelled
    # no position updates
    @fact posactual[1] --> 0

    # limit order was pending that is still pending
    # stop price reached
    # limit-target and prev. limit-pending can not match
    # (stop < limit for buy, sell symmetrically)
    # here we have (buy) limit_prev < stop < price_curr
    ord = TradingLogic.Order("pos", 50, 70.0, :buy,
                             :limit, :pending, "aux")
    @fact TradingLogic.orderhandling!(
      targ, 75.1, tnow(), posactual, ord,
      blotter, backtest) --> (true, shortgain)
    # cancel previous order with different limit-price
    @fact TradingLogic.ispending(ord) --> false
    @fact ord.status --> :cancelled
    # no position updates
    @fact posactual[1] --> 0
  end

  context("Order handling: stoplimit to track, sell-side") do
    targ = (-50, [60.0, 75.0]) #[limitprice, stopprice]

    # nothing was pending
    # stop price not reached
    ord = TradingLogic.emptyorder()
    @fact TradingLogic.ispending(ord) --> false
    @fact TradingLogic.orderhandling!(
      targ, 80.0, tnow(), posactual, ord,
      blotter, backtest) --> (true, shortgain)
    # no change, keep tracking the price
    @fact TradingLogic.ispending(ord) --> false
    @fact posactual[1] --> 0

    # nothing was pending
    # stop price reached
    ord = TradingLogic.emptyorder()
    @fact TradingLogic.ispending(ord) --> false
    @fact TradingLogic.orderhandling!(
      targ, 74.9, tnow(), posactual, ord,
      blotter, backtest) --> (true, shortgain)
    # now limit pending, actual position still the same
    @fact TradingLogic.ispending(ord) --> true
    @fact ord.ordertype --> :limit
    @fact ord.side --> :sell
    @fact ord.price --> roughly(60.0)
    @fact posactual[1] --> 0

    # limit order was pending that is still pending
    # stop price not reached
    ord = TradingLogic.Order("neg", 30, 85.0, :sell,
                             :limit, :pending, "aux")
    @fact TradingLogic.orderhandling!(
      targ, 80.0, tnow(), posactual, ord,
      blotter, backtest) --> (true, shortgain)
    # cancel previous order, keep tracking the price
    @fact TradingLogic.ispending(ord) --> false
    @fact ord.status --> :cancelled
    # no position updates
    @fact posactual[1] --> 0

    # limit order was pending that is still pending
    # stop price reached
    # limit-target and prev. limit-pending can not match
    # (stop < limit for buy, sell symmetrically)
    # here we have (sell) limit_prev > stop > price_curr
    ord = TradingLogic.Order("neg", 50, 80.0, :sell,
                             :limit, :pending, "aux")
    @fact TradingLogic.orderhandling!(
      targ, 74.9, tnow(), posactual, ord,
      blotter, backtest) --> (true, shortgain)
    # cancel previous order with different limit-price
    @fact TradingLogic.ispending(ord) --> false
    @fact ord.status --> :cancelled
    # no position updates
    @fact posactual[1] --> 0
  end
  #println(blotter)
end


text:
type SettlementPhysical <: SettlementType end
type SettlementCash <: SettlementType end

type SwaptionResults{S <: AbstractString}
  value::Float64
  additionalResults::Dict{S, Float64}
end

SwaptionResults() = SwaptionResults(0.0, Dict("spreadCorrection" => 0.0, "strike" => 0.0, "atmForward" => 0.0, "annuity" => 0.0, "swapLength" => 0.0, "stdDev" => 0.0, "vega" => 0.0))

function reset!(results::SwaptionResults)
  results.value = 0.0
  for k in keys(results.additionalResults)
    results.additionalResults[k] = 0.0
  end

  return results
end

type Swaption{E <: Exercise, S <: SettlementType, P <: PricingEngine} <: Option
  lazyMixin::LazyMixin
  swap::VanillaSwap
  exercise::E
  delivery::S
  results::SwaptionResults
  pricingEngine::P

  # function call{E, S}(::Type{Swaption}, lz::LazyMixin, swap::VanillaSwap, exercise::E, delivery::S, results::SwaptionResults)
  #   new{E, S, PricingEngine}(lz, swap, exercise, delivery, results)
  # end

  function call{E, S, P}(::Type{Swaption}, lz::LazyMixin, swap::VanillaSwap, exercise::E, delivery::S, results::SwaptionResults, pricingEngine::P = NullSwaptionEngine())
    new{E, S, P}(lz, swap, exercise, delivery, results, pricingEngine)
  end
end

Swaption{E <: Exercise}(swap::VanillaSwap, exercise::E) = Swaption(LazyMixin(), swap, exercise, SettlementPhysical(), SwaptionResults())
Swaption{E <: Exercise, P <: PricingEngine}(swap::VanillaSwap, exercise::E, pe::P) = Swaption(LazyMixin(), swap, exercise, SettlementPhysical(), SwaptionResults(), pe)

type NonstandardSwaption{E <: Exercise, S <: SettlementType, P <: PricingEngine} <: Option
  lazyMixin::LazyMixin
  swap::NonstandardSwap
  exercise::E
  delivery::S
  results::SwaptionResults
  pricingEngine::P

  function call{E, S, P}(::Type{NonstandardSwaption}, lz::LazyMixin, swap::NonstandardSwap, exercise::E, delivery::S, results::SwaptionResults, pricingEngine::P = NullSwaptionEngine())
    new{E, S, P}(lz, swap, exercise, delivery, results, pricingEngine)
  end
end

NonstandardSwaption{E <: Exercise}(swap::NonstandardSwap, exercise::E) = NonstandardSwaption(LazyMixin(), swap, exercise, SettlementPhysical(), SwaptionResults())
NonstandardSwaption{E <: Exercise, P <: PricingEngine}(swap::NonstandardSwap, exercise::E, pe::P) = NonstandardSwaption(LazyMixin(), swap, exercise, SettlementPhysical(), SwaptionResults(), pe)

function perform_calculations!(swaption::Union{Swaption, NonstandardSwaption})
  reset!(swaption.results)
  _calculate!(swaption.pricingEngine, swaption)

  return swaption
end

function npv(swaption::Union{Swaption, NonstandardSwaption})
  calculate!(swaption)

  return swaption.results.value
end

function calibration_basket(swaption::NonstandardSwaption, swaptionEngine::PricingEngine, swapIndex::SwapIndex, swaptionVol::SwaptionVolatilityStructure,
                            basketType::CalibrationBasketType)
  return calibration_basket(swaptionEngine, swaption, swapIndex, swaptionVol, basketType)
end

function clone(swaption::Swaption, pe::PricingEngine = swaption.pricingEngine)
  lazyMixin, res = pe == swaption.pricingEngine ? (swaption.lazyMixin, swaption.results) : (LazyMixin(), SwaptionResults())

  newSwaption = Swaption(lazyMixin, swaption.swap, swaption.exercise, swaption.delivery, res, pe)

  return newSwaption
end

get_pricing_engine_type{E, S, P}(::Swaption{E, S, P}) = P


text:
"""
    heapify(array::Array{Int64,1})

The heapify function merges the two subroutines for:
    1) Put elements of the array in heap order
    2) Repair the heap elements `i` whose root element, if is not at the index `largest`
For more information see: https://en.wikipedia.org/wiki/Heapsort

# Arguments
- `array::Array{Int64,1}`: Unsorted array of integers
- `n::Int64`: Length of the array
- `i::Int64`: Current index of the array

"""
function heapify(array::Array{Int64,1}, n::Int64, i::Int64)
    # Initialize largest as root and left and right subarrays
    largest = i
    # Julia notations, which starts with 1
    left = 2 * i
    right = 2 * i + 1

    # See if left child of root exists and is greater than root
    if left <= n && array[i] < array[left]
        largest = left
    end
    # See if right child of root exists and is greater than root
    if right <= n && array[largest] < array[right]
        largest = right
    end
    # Change root, if needed
    if largest != i
        array[i], array[largest] = array[largest], array[i]
        # Heapify the root again
        array = heapify(array, n, largest)
    end
    return array
end

"""
    heap_sorting(array::Array{Int64,1})

As a comparison-based sorting algorithm, the heapsort algorithm (HSA) divides its input
into a sorted and an unsorted region, and it iteratively shrinks the unsorted region by
extracting the largest element from it and inserting it into the sorted region. A specialty
is that the HSA keeps the unsorted region in a heap data structure to find the largest
element in each step more quickly. In more detail, in the first part of the HSA
(while-loop), the largest value has to be found and set to position one. In the second part
of the HSA (while-loop), the array's first and largest value has to be swap to the last
index of the array, and the swapping-procedure starts again for a new interval n-1. For
more information see: [https://en.wikipedia.org/wiki/Heapsort](https://en.wikipedia.org/wiki/Heapsort)

# Arguments
- `array::Array{Int64,1}`: Unsorted array of integers

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> arr = [64, 34, 25, 12, 22, 11, 90]
julia> ClassicAlgorithmsCollections.heap_sorting(arr)
[11, 12, 22, 25, 34, 64, 90]
```
"""
function heap_sorting(array::Array{Int64,1})
    # Base case. A array of zero or one elements is sorted, by definition.
    n = length(array)

    # Shrink the active space
    i = fld(n, 2)
    # Build a maxheap
    while i > 0
        array = heapify(array, n, i)
        i -= 1

    end

    # One by one extract elements
    i = n
    while i > 1
        array[i], array[1] = array[1], array[i]
        # Julia specific reducing effectiv index
        array = heapify(array, i - 1, 1)
        i -= 1
    end

    return array
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
# GPUArrays.jl interface

#
# Device functionality
#

## execution

struct oneArrayBackend <: AbstractGPUBackend end

struct oneKernelContext <: AbstractKernelContext end

@inline function GPUArrays.launch_heuristic(::oneArrayBackend, f::F, args::Vararg{Any,N};
                                             elements::Int, elements_per_thread::Int) where {F,N}
    kernel = @oneapi launch=false f(oneKernelContext(), args...)

    items = suggest_groupsize(kernel.fun, elements).x
    # XXX: the z dimension of the suggested group size is often non-zero.
    #      preserve this in GPUArrays?
    # XXX: how many groups is a good number? the API doesn't tell us.
    #      measured on a low-end IGP, 32 blocks seems like a good sweet spot.
    #      note that this only matters for grid-stride kernels, like broadcast.
    return (threads=items, blocks=32)
end

function GPUArrays.gpu_call(::oneArrayBackend, f, args, threads::Int, blocks::Int;
                            name::Union{String,Nothing})
    @oneapi items=threads groups=blocks name=name f(oneKernelContext(), args...)
end

## on-device

# indexing

GPUArrays.blockidx(ctx::oneKernelContext) = oneAPI.get_group_id(0)
GPUArrays.blockdim(ctx::oneKernelContext) = oneAPI.get_local_size(0)
GPUArrays.threadidx(ctx::oneKernelContext) = oneAPI.get_local_id(0)
GPUArrays.griddim(ctx::oneKernelContext) = oneAPI.get_num_groups(0)

# math

@inline GPUArrays.cos(ctx::oneKernelContext, x) = oneAPI.cos(x)
@inline GPUArrays.sin(ctx::oneKernelContext, x) = oneAPI.sin(x)
@inline GPUArrays.sqrt(ctx::oneKernelContext, x) = oneAPI.sqrt(x)
@inline GPUArrays.log(ctx::oneKernelContext, x) = oneAPI.log(x)

# memory

@inline function GPUArrays.LocalMemory(::oneKernelContext, ::Type{T}, ::Val{dims}, ::Val{id}
                                      ) where {T, dims, id}
    ptr = oneAPI.emit_localmemory(Val(id), T, Val(prod(dims)))
    oneDeviceArray(dims, LLVMPtr{T, onePI.AS.Local}(ptr))
end

# synchronization

@inline GPUArrays.synchronize_threads(::oneKernelContext) = oneAPI.barrier()

#
# Host abstractions
#

GPUArrays.backend(::Type{<:oneArray}) = oneArrayBackend()

const GLOBAL_RNGs = Dict{ZeDevice,GPUArrays.RNG}()
function GPUArrays.default_rng(::Type{<:oneArray})
    dev = device()
    get!(GLOBAL_RNGs, dev) do
        N = oneL0.compute_properties(dev).maxTotalGroupSize
        state = oneArray{NTuple{4, UInt32}}(undef, N)
        rng = GPUArrays.RNG(state)
        Random.seed!(rng)
        rng
    end
end


text:
using Base.Test
import WaterData

@testset "Phase boundaries" begin
    @testset "have correct long-short code mappings" begin
        iskeyof(dict) = k -> k in keys(dict)
        isvaluein(dict) = v -> v in values(dict)
        phasedict = WaterData.phase_mappings

        @test "liquid" |> iskeyof(phasedict)
        @test "L" |> isvaluein(phasedict)
        @test "L" |> iskeyof(phasedict)
        @test "ice VII" |> iskeyof(phasedict)
        @test "VII" |> iskeyof(phasedict)
        @test "VII" |> isvaluein(phasedict)
        @test phasedict["liquid"] == "L"
        @test phasedict["ice VII"] == "VII"
    end

    @testset "have correct end values on the boundaries" begin
        pb = WaterData.PhaseBoundary

        maxtemp(pb) = maximum(pb.T)
        mintemp(pb) = minimum(pb.T)
        maxpress(pb) = maximum(pb.P)
        minpress(pb) = minimum(pb.P)

        @test_throws BoundsError pb(:L, :II)
        @test_throws BoundsError pb("L", "II")

        @test isapprox(mintemp(pb(:L, :VII)), 355.0, rtol=0.01)
        @test isapprox(maxpress(pb(:L, :VII)), 400000e5, rtol=0.01)
        @test isapprox(mintemp(pb("L", :VII)), 355.0, rtol=0.01)
        @test isapprox(minpress(pb("L", :VII)), 22160e5, rtol=0.01)
        @test isapprox(mintemp(pb(:L, "VII")), 355.0, rtol=0.01)
        @test isapprox(mintemp(pb("L", "VII")), 355.0, rtol=0.01)
        @test isapprox(maxtemp(pb(:VII, :X)), 1500, rtol=0.01)
    end

    @testset "Testing points on and across the boundary" begin
        Ps = [1,2,3,4,5]
        Ts = [5,4,3,2,1]
        pb = WaterData.OtherPhaseBoundary(Ps, Ts)

        line1 = ([1, 1], [5, 5])
        line2 = ([1, 1], [2, 5])
        line3 = ([0.5, 0.5], [0.5, 6])
        line4 = ([4, 4], [6, 6])
        line5 = ([4, 1], [4, 3])
        line6 = ([2, 3], [3, 3])

        @test WaterData.intersects(line1, pb)
        @test WaterData.intersects(line2, pb)
        @test !WaterData.intersects(line3, pb)
        @test !WaterData.intersects(line4, pb)
        @test WaterData.intersects(line5, pb)
        # we treat "just touches the boundary" as "intersects"
        @test WaterData.intersects(line6, pb)
    end
end


text:
"""
    orthophoto(
    	txtpotreedirs::String,
    	outputimage::String,
    	bbin::Union{String,AABB},
    	GSD::Float64,
    	PO::String,
    	quota::Union{Float64,Nothing},
    	thickness::Union{Float64,Nothing},
    	ucs::Union{String,Matrix{Float64}},
    	BGcolor::Array{Float64,1},
    	pc::Bool
    	)

Orthographic projection of point cloud.

Description:
 - txtpotreedirs: path to Potree project or a text file containing a list of Potree projects
 - outputimage: output filename for image
 - bbin: region of interest, aabb ("xmin ymin zmin xmax ymax zmax") or volume (path to json file of described volume)
 - GSD: ground sampling distance
 - PO: projection plane and sight direction. Option: "XY+","XY-","XZ+","XZ-","YZ+","YZ-"
 - ucs: user coordinate system
 - BGcolor: background color of image
 - pc: if true saves the segmented point cloud in region of interest

It is possible to limit the region of interest for a section of point cloud. The slice is described with:
 - quota: distance of plane to origin
 - thickness: thickness of plane

Put both of them to nothing if you are not interested in a slice.
"""
function orthophoto(
    txtpotreedirs::String,
    outputimage::String,
    bbin::Union{String,AABB},
    GSD::Float64,
    PO::String,
    altitude::Union{Float64,Nothing},
    thickness::Union{Float64,Nothing},
    ucs::Union{String,Matrix{Float64}},
    BGcolor::Array{Float64,1},
    pc::Bool,
    epsg::Union{Nothing,Integer},
)

    # initialization
    params = ParametersOrthophoto(
        txtpotreedirs,
        outputimage,
        bbin,
        GSD,
        PO,
        altitude,
        thickness,
        ucs,
        BGcolor,
        pc,
        epsg,
    )

    proj_folder = splitdir(params.outputfile)[1]

    if params.pc
        temp = joinpath(proj_folder, "tmp.las")
        params.stream_tmp = open(temp, "w")
    end

    flushprintln(" ")
    flushprintln("========= PROCESSING =========")

    for potree in params.potreedirs
        traversal(potree, params)
    end

    if params.pc
        close(params.stream_tmp)
    end

    flushprintln(" ")
    flushprintln("========= SAVES =========")

    # saves image
    saveimage(params)

    FileManager.successful(params.numPointsProcessed != 0, proj_folder::String)
    flushprintln("Processed $(params.numPointsProcessed) points")

    # saves point cloud extracted
    if pc
        savepointcloud(params, temp)
    end

end


text:
function get_branch_to_pm(ix::Int64, branch::PowerSystems.Transformer2W)
    PM_branch = Dict{String,Any}(
        "br_r"        => branch.r,
        "rate_a"      => branch.rate,
        "shift"       => 0.0,
        "rate_b"      => branch.rate,
        "br_x"        => branch.r,
        "rate_c"      => branch.rate,
        "g_to"        => 0.0,
        "g_fr"        => 0.0,
        "b_fr"        => branch.primaryshunt/2,
        "f_bus"       => branch.connectionpoints.from.number,
        "br_status"   => Float64(branch.available),
        "t_bus"       => branch.connectionpoints.to.number,
        "b_to"        => branch.primaryshunt/2,
        "index"       => ix,
        "angmin"      => -1.50,
        "angmax"      =>  1.50,
        "transformer" => true,
        "tap"         => 1.0,
    )
    return PM_branch
end

function get_branch_to_pm(ix::Int64, branch::PowerSystems.TapTransformer)
    PM_branch = Dict{String,Any}(
        "br_r"        => branch.r,
        "rate_a"      => branch.rate,
        "shift"       => 0.0,
        "rate_b"      => branch.rate,
        "br_x"        => branch.r,
        "rate_c"      => branch.rate,
        "g_to"        => 0.0,
        "g_fr"        => 0.0,
        "b_fr"        => branch.primaryshunt/2,
        "f_bus"       => branch.connectionpoints.from.number,
        "br_status"   => Float64(branch.available),
        "t_bus"       => branch.connectionpoints.to.number,
        "b_to"        => branch.primaryshunt/2,
        "index"       => ix,
        "angmin"      => -1.50,
        "angmax"      =>  1.50,
        "transformer" => true,
        "tap"         => branch.tap
    )
    return PM_branch
end

function get_branch_to_pm(ix::Int64, branch::PowerSystems.Line)
    PM_branch = Dict{String,Any}(
        "br_r"        => branch.r,
        "rate_a"      => branch.rate.from_to,
        "shift"       => 0.0,
        "rate_b"      => branch.rate.from_to,
        "br_x"        => branch.r,
        "rate_c"      => branch.rate.from_to,
        "g_to"        => 0.0,
        "g_fr"        => 0.0,
        "b_fr"        => branch.b.from,
        "f_bus"       => branch.connectionpoints.from.number,
        "br_status"   => Float64(branch.available),
        "t_bus"       => branch.connectionpoints.to.number,
        "b_to"        => branch.b.to,
        "index"       => ix,
        "angmin"      => branch.anglelimits.min,
        "angmax"      => branch.anglelimits.max,
        "transformer" => false,
        "tap"         => 1.0,
    )
    return PM_branch
end

function get_branches_to_pm(branches::Array{T}) where {T <: PowerSystems.Branch}
        PM_branches = Dict{String,Any}()

        for (ix, branch) in enumerate(branches)
            PM_branches["$(ix)"] = get_branch_to_pm(ix, branch)
        end

    return PM_branches
end

function get_buses_to_pm(buses::Array{PowerSystems.Bus})
    PM_buses = Dict{String,Any}()
    for bus in buses
        PM_bus = Dict{String,Any}(
        "zone"     => 1,
        "bus_i"    => bus.number,
        "bus_type" => 2,
        "vmax"     => bus.voltagelimits.max,
        "area"     => 1,
        "vmin"     => bus.voltagelimits.min,
        "index"    => bus.number,
        "va"       => bus.angle,
        "vm"       => bus.voltage,
        "base_kv"  => bus.basevoltage,
        "pni"      => 0.0,
        "qni"      => 0.0,
        )
        PM_buses["$(bus.number)"] = PM_bus
    end
    PM_buses["$(buses[1].number)"]["bus_type"] = 3
    return PM_buses
end

#=
function expression_to_pm_active(PM_dict::Dict{String,Any}, netinjection::BalanceNamedTuple, sys::PowerSystems.PowerSystem)

    for bus in sys.buses, time in 1:sys.time_periods

        PM_dict["nw"]["$(time)"]["bus"]["$(bus.number)"]["pni"] = netinjection.var_active[bus.number,time]

    end

end

function expression_to_pm_reactive(PM_dict::Dict{String,Any}, netinjection::BalanceNamedTuple, sys::PowerSystems.PowerSystem)

    for bus in sys.buses, time in 1:sys.time_periods

        PM_dict["nw"]["$(time)"]["bus"]["$(bus.number)"]["qni"] = netinjection.var_reactive[bus.number,time]

    end

end
=#

function pass_to_pm(sys::PowerSystems.PowerSystem, netinjection::BalanceNamedTuple)

    PM_translation = Dict{String,Any}(
    "bus" => get_buses_to_pm(sys.buses),
    "branch" => get_branches_to_pm(sys.branches),
    "baseMVA" => sys.basepower,
    "per_unit" => true,
    "dcline"         => Dict{String,Any}(),
    "gen"            => Dict{String,Any}(),
    "shunt"          => Dict{String,Any}(),
    "load"           => Dict{String,Any}(),
    )

    # TODO: this function adds overhead in large number of time_steps
    # We can do better later.

    PM_translation = IM.replicate(PM_translation,sys.time_periods)

    return PM_translation

end


text:
using Optim

function lmfit{I <: Integer}(f::Function, p0::Vector{Float64}, maxIter::I)
	# this is a convenience function for the curve_fit() methods
	# which assume f(p) is the cost functionj i.e. the residual of a
	# model where
	#   model(xpts, params...) = ydata + error (noise)

	# this minimizes f(p) using a least squares sum of squared error:
	#   sse = sum(f(p)^2)
	# This is currently embedded in Optim.levelberg_marquardt()
	# which calls sumabs2
	#
	# returns p, f(p), g(p) where
	#   p    : best fit parameters
	#   f(p) : function evaluated at best fit p, (weighted) residuals
	#   g(p) : estimated Jacobian at p (Jacobian with respect to p)

	# construct Jacobian function, which uses finite difference method
	g = Calculus.jacobian(f)

	results = Optim.levenberg_marquardt(f, g, p0; maxIter=maxIter, tolG=1e-8, lambda=1.0)
	# p = results.minimum
	# resid = f(p)
	# dof = length(resid) - length(p)

  return results
end