text:
#=
This code is loosely based on a routine originally copyright Federal Reserve Bank of Atlanta
and written by Iskander Karibzhanov.
=#

"""
```
kalman_filter(data, lead, a, F, b, H, var, z0, vz0, Ny0; allout=false)
kalman_filter(data, lead, a, F, b, H, var, Ny0=0; allout=false)
```

Inputs
------

- `data`: a `Ny x T` matrix containing data `y(1), ... , y(T)`.
- `lead`: the number of steps to forecast after the end of the data.
- `a`: an `Nz x 1` vector for a time-invariant input vector in the transition equation.
- `F`: an `Nz x Nz` matrix for a time-invariant transition matrix in the transition
  equation.
- `b`: an `Ny x 1` vector for a time-invariant input vector in the measurement equation.
- `H`: an `Ny x Nz` matrix for a time-invariant measurement matrix in the measurement
  equation.
- `var`: an `Ny + Nz` x `Ny + Nz` matrix for a time-invariant variance matrix for the
  error in the transition equation and the error in the measurement equation, that is,
  `[η(t)', ϵ(t)']'`.
- `z0`: an optional `Nz x 1` initial state vector.
- `vz0`: an optional `Nz x Nz` covariance matrix of an initial state vector.
- `Ny0`: an optional scalar indicating the number of periods of presample (i.e. the number
  of periods which we don't add to the likelihood)
- `allout`: an optional keyword argument indicating whether we want optional output
  variables returned as well

Outputs
-------

- `logl`: value of the average log likelihood function of the SSM under assumption that
  observation noise ϵ(t) is normally distributed
- `pred`: a `Nz` x `T+lead` matrix containing one-step predicted state vectors.
- `vpred`: a `Nz` x `Nz` x `T+lead` matrix containing mean square errors of predicted
  state vectors.
- `filt`: an optional `Nz` x `T` matrix containing filtered state vectors.
- `vfilt`: an optional `Nz` x `Nz` x `T` matrix containing mean square errors of filtered
  state vectors.

Notes
-----
The state space model is defined as follows:
```
z(t+1) = a+F*z(t)+η(t)     (state or transition equation)
y(t) = b+H*z(t)+ϵ(t)       (observation or measurement equation)
```

When `z0` and `Vz0` are omitted, the initial state vector and its covariance matrix of the
time invariant Kalman filters are computed under the stationarity condition:
```
z0 = (I-F)\a
vz0 = (I-kron(F,F))\(V(:),Nz,Nz)
```
where `F` and `V` are the time invariant transition matrix and the covariance matrix of
transition equation noise, and `vec(V)` is an `Nz^2` x `1` column vector that is constructed
by stacking the `Nz` columns of `V`.  Note that all eigenvalues of `F` are inside the unit
circle when the state space model is stationary.  When the preceding formula cannot be
applied, the initial state vector estimate is set to `a` and its covariance matrix is given
by `1E6I`.  Optionally, you can specify initial values.
"""
function kalman_filter{S<:AbstractFloat}(data::Matrix{S},
                                      lead::Int64,
                                      a::Matrix{S},
                                      F::Matrix{S},
                                      b::Matrix{S},
                                      H::Matrix{S},
                                      var::Matrix{S},
                                      z0::Matrix{S},
                                      vz0::Matrix{S},
                                      Ny0::Int = 0;
                                      allout::Bool = false)
    T = size(data, 2)
    Nz = size(a, 1)
    Ny = size(b, 1)

    z = z0
    P = vz0

    # Check input matrix dimensions
    @assert size(data, 1) == Ny
    @assert size(a, 2) == 1
    @assert size(F) == (Nz, Nz)
    @assert size(b, 2) == 1
    @assert size(H) == (Ny, Nz)
    @assert size(var) == (Ny + Nz, Ny + Nz)
    @assert size(z) == (Nz, 1)
    @assert size(P) == (Nz, Nz)

    # V(t) and R(t) are variances of η(t) and ϵ(t), respectively, and G(t) is a covariance
    # of η(t) and ϵ(t)
    # In dsgelh :
    # --- V is same as QQ
    # --- R is same as EE
    # --- G is same as VV = QQ*MM
    V = var[1:Nz, 1:Nz]
    R = var[(Nz+1):end, (Nz+1):end]
    G = var[1:Nz, (Nz+1):end]

    if allout
        pred          = zeros(S, Nz, T)
        vpred         = zeros(S, Nz, Nz, T)
        yprederror    = NaN*zeros(S, Ny, T)
        ystdprederror = NaN*zeros(S, Ny, T)
        filt          = zeros(Nz, T)
        vfilt         = zeros(Nz, Nz, T)
    end

    L = zero(S)

    for t = 1:T
        # If an element of the vector y(t) is missing (NaN) for the observation t, the
        #   corresponding row is ditched from the measurement equation.
        nonmissing = !isnan(data[:, t])
        data_t = data[nonmissing, t]       # data_t = Y_T = [y1, y2, ..., yT] is matrix of observable data time-series
        H_t = H[nonmissing, :]             # H_t = DD is matrix mapping states to observables
        G_t = G[:, nonmissing]             # G_t = Cov(η_t, ϵ_t)
        R_t = R[nonmissing, nonmissing]    # R_t = Var(ϵ_t)
        Ny_t = length(data_t)              # Ny_t = T is length of time
        b_t = b[nonmissing, :]             # b_t = DD

        ## forecasting
        z = a + F*z                        # z_{t|t-1} = a + F(Θ)*z_{t-1|t-1}
        P = F*P*F' + V                     # P_{t|t-1} = F(Θ)*P_{t-1|t-1}*F(Θ)' + F(Θ)*Var(η_t)*F(Θ)'
        dy = data_t - H_t*z - b_t          # dy = y_t - H(Θ)*z_{t|t-1} - DD is prediction error or innovation
        HG = H_t*G_t                       # HG is ZZ*Cov(η_t, ϵ_t)
        D = H_t*P*H_t' + HG + HG' + R_t    # D = ZZ*P_{t+t-1}*ZZ' + HG + HG' + R_t
        D = (D+D')/2

        if allout
            pred[:, t]                   = z
            vpred[:, :, t]               = P
            yprederror[nonmissing, t]    = dy
            ystdprederror[nonmissing, t] = dy./sqrt(diag(D))
        end

        ddy = D\dy
        # We evaluate the log likelihood function by adding values of L at every iteration
        #   step (for each t = 1,2,...T)
        if t > Ny0
            L += -log(det(D))/2 - first(dy'*ddy/2) - Ny_t*log(2*pi)/2
        end

        ## updating
        PHG = P*H_t' + G_t
        z = z + PHG*ddy                    # z_{t|t} = z_{t|t-1} + P_{t|t-1}*H(Θ)' + ...
        P = P - PHG/D*PHG'                 # P_{t|t} = P_{t|t-1} - PHG*(1/D)*PHG

        if allout
            PH = P*H_t'
            filt[:, t]     = z
            vfilt[:, :, t] = P
        end
    end

    zend = z
    Pend = P

    if allout && lead > 1
        for t = (T+2):(T+lead)
            z = F*z + a
            P = F*P*F' + V
            pred[:, t]     = z
            vpred[:, :, t] = P
        end
    end

    if allout
        rmse = sqrt(mean((yprederror.^2)', 1))
        rmsd = sqrt(mean((ystdprederror.^2)', 1))
        return Kalman(L, zend, Pend, pred, vpred, yprederror, ystdprederror, rmse, rmsd, filt, vfilt)
    else
        return Kalman(L, zend, Pend)
    end

end

function kalman_filter{S<:AbstractFloat}(data::Matrix{S},
                                      lead::Int64,
                                      a::Matrix{S},
                                      F::Matrix{S},
                                      b::Matrix{S},
                                      H::Matrix{S},
                                      var::Matrix{S},
                                      Ny0::Int = 0;
                                      allout::Bool = false)
    Nz = size(a, 1)
    V = var[1:Nz, 1:Nz]

    e, _ = eig(F)
    if countnz(e*e' - eye(Nz)) == Nz^2
        z0 = (eye(Nz) - F)\a
        vz0 = reshape((eye(Nz^2)-kron(F,F))\V, Nz, Nz)
    else
        z0 = a
        vz0 = eye(Nz)*1e6
    end

    return kalman_filter(data, lead, a, F, b, H, var, z0, vz0, Ny0; allout=allout)
end

immutable Kalman{S<:AbstractFloat}
    L::S
    zend::Matrix{S}
    Pend::Matrix{S}
    pred::Matrix{S}
    vpred::Array{S,3}
    yprederror::Matrix{S}
    ystdprederror::Matrix{S}
    rmse::Matrix{S}
    rmsd::Matrix{S}
    filt::Matrix{S}
    vfilt::Array{S,3}
end
function Kalman{S<:AbstractFloat}(L::S,
                                  zend::Matrix{S},
                                  Pend::Matrix{S},
                                  pred::Matrix{S}          = Matrix{S}(),
                                  vpred::Array{S,3}        = Array{S}(0,0,0),
                                  yprederror::Matrix{S}    = Matrix{S}(),
                                  ystdprederror::Matrix{S} = Matrix{S}(),
                                  rmse::Matrix{S}          = Matrix{S}(),
                                  rmsd::Matrix{S}          = Matrix{S}(),
                                  filt::Matrix{S}          = Matrix{S}(),
                                  vfilt::Array{S,3}        = Array{S}(0,0,0))
    return Kalman{S}(L,zend,Pend,pred,vpred,yprederror,ystdprederror,rmse,rmsd,filt,vfilt)
end
function Base.getindex(K::Kalman, d::Symbol)
    if d in (:L, :zend, :Pend, :pred, :vpred, :yprederror, :ystdprederror, :rmse, :rmsd,
             :filt, :vfilt)
        return getfield(K, d)
    else
        throw(KeyError(d))
    end
end


text:
#= Neural Network Variations
Mainly tuning the Multi-Layer Perceptron
A big limitation is that this doesn't use GPU at all
since ScikitLearn doesn't support it
So deeper models will take forever even on a good CPU
And I don't even want to try the 13-feature data set
-Matei
=#

using DataFrames, CSV, ScikitLearn, PyPlot, Statistics

# Data Loading
################################################################################

data = CSV.File("audio_recoded_combined_binary.csv") |> DataFrame

X = convert(Array, data[!,Not(:Class)])
y = convert(Array, data[!,:Class]) # :Class is our target variable

@sk_import model_selection: train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=42)

@sk_import neural_network: MLPClassifier

#=
mlp_1layer10  = MLPClassifier(hidden_layer_sizes=(10))
mlp_1layer25  = MLPClassifier(hidden_layer_sizes=(25))
mlp_1layer40  = MLPClassifier(hidden_layer_sizes=(40))
mlp_1layer50  = MLPClassifier(hidden_layer_sizes=(50))
mlp_1layer100 = MLPClassifier(hidden_layer_sizes=(100))

mlp_6layer    = MLPClassifier(hidden_layer_sizes=(30, 50, 60, 10, 10, 10))
mlp_5layer    = MLPClassifier(hidden_layer_sizes=(50, 10, 10, 10, 10, 10))
=#

@sk_import metrics: classification_report
@sk_import metrics: plot_confusion_matrix
@sk_import model_selection: cross_val_score

function trainModel(model)
    fit!(model,X_train,y_train)
    # Training performance
    y_pred = predict(model,X_train)
    print(classification_report(y_train,y_pred))
    # Confusion matrices
    plot_confusion_matrix(model,X_test,y_test)
    PyPlot.gcf()
end

function crossValidateModel(model,folds=2)
    #fit!(model,X_train,y_train)
    @time cva = cross_val_score(model, X_train, y_train, cv=folds)
    println("mean accuracy: \n" *  string(mean(cva)))
    @time cvp = cross_val_score(model, X_train, y_train, cv=folds, scoring="precision")
    println("mean precision: \n" * string(mean(cvp)))
    @time cvr = cross_val_score(model, X_train, y_train, cv=folds, scoring="recall")
    println("mean recall: \n" *    string(mean(cvr)))
end

function testModel(model)
    # Testing performance
    y_pred = predict(model,X_test)
    print(classification_report(y_test,y_pred))
    # Confusion matrices
    plot_confusion_matrix(model,X_test,y_test)
    PyPlot.gcf()
end

mlp_2layer    = MLPClassifier(hidden_layer_sizes=(50, 10))
trainModel(         mlp_2layer)
crossValidateModel( mlp_2layer)
testModel(          mlp_2layer)

mlp_2layerlarge    = MLPClassifier(hidden_layer_sizes=(100, 50))
trainModel(         mlp_2layerlarge)
crossValidateModel( mlp_2layerlarge)
testModel(          mlp_2layerlarge)

mlp_3layer    = MLPClassifier(hidden_layer_sizes=(50, 10, 10))
trainModel(         mlp_3layer)
crossValidateModel( mlp_3layer)
testModel(          mlp_3layer)

mlp_4layer    = MLPClassifier(hidden_layer_sizes=(50, 10, 10, 10))
trainModel(         mlp_4layer)
crossValidateModel( mlp_4layer)
testModel(          mlp_4layer)

mlp_5layer    = MLPClassifier(hidden_layer_sizes=(50, 10, 10, 10, 10, 10))
trainModel(         mlp_5layer)
crossValidateModel( mlp_5layer)
testModel(          mlp_5layer)

testModel(mlp_1layer50)
@sk_import model_selection: cross_val_score
@time cva = cross_val_score(mlp_1layer50, X_train, y_train, cv=10)
println("mean accuracy:" * mean(cva))
mlp_1layer50_p  = MLPClassifier(hidden_layer_sizes=(50))
fit!(model,X_train,y_train)
@time cvp = cross_val_score(mlp_1layer50_p, X_train, y_train, cv=10, scoring="precision")
println("mean precision:" * mean(cvp))
@time cvr = cross_val_score(model, X_train, y_train, cv=10, scoring="recall")
println("mean recall:" * mean(cvr))

# Cross validation
################################################################################
#=
@sk_import model_selection: cross_val_score

function crossValidateModel(model)

    #fit!(model,X_train,y_train)
    @time cva = cross_val_score(model, X_train, y_train, cv=10)
    println("mean accuracy:" * mean(cva))
    @time cvp = cross_val_score(model, X_train, y_train, cv=10, scoring="precision")
    println("mean precision:" * mean(cvp))
    @time cvr = cross_val_score(model, X_train, y_train, cv=10, scoring="recall")
    println("mean recall:" * mean(cvr))

end

@time crossValidateModel(mlp_1layer50)
=#
# Validation and learning curves
# Will get this to work eventually -Matei
################################################################################

#=
@sk_import model_selection: validation_curve
@sk_import model_selection: learning_curve

using Plots, StatsPlots

validation_curve(mlp_1layer50, X, y,
    param_name="gamma", param_range=np.logspace(-6, -1, 5), scoring="accuracy")

fit!(mlp_1layer50,X_train,y_train)
train_sizes, train_scores, valid_scores = learning_curve(
    mlp_1layer50, X, y, train_sizes=[100, 500, 1000], cv=2)

plot(train_scores, valid_scores)
=#

# Random Forest Comparison
################################################################################

@sk_import ensemble: RandomForestClassifier
rfc    = RandomForestClassifier()
trainModel(         rfc)
crossValidateModel( rfc)
testModel(          rfc)

ext: jl


text:
#### Teacup GIF
## Inspired by an example of Lorez attractor: http://docs.juliaplots.org/latest/
## FFmpeg has to be installed prior to running the code

## Preparation

# Activate project environment
import Pkg; Pkg.activate(".")

# Load required packages
using DifferentialEquations
using Plots

# Load Hastings & Powell's model
include("HPmodel.jl")

## Produce the GIF

# define the model
mutable struct TeaCup
    dt; a1; a2; b1; b2; d1; d2; x; y; z
end

function step!(t::TeaCup)
    dx = t.x*(1-t.x) - t.a1*t.x*t.y/(1 + t.b1*t.x)     ; t.x += t.dt * dx
    dy = t.a1*t.x*t.y/(1 + t.b1*t.x) - t.a2*t.y*t.z/(1 + t.b2*t.y) - t.d1*t.y ; t.y += t.dt * dy
    dz = t.a2*t.y*t.z/(1 + t.b2*t.y) - t.d2*t.z   ; t.z += t.dt * dz
end

# solution from initial values
TeaCup_sol = TeaCup((dt = 0.05, a1 = 5.0, a2 = 0.1, b1 = 3.0, b2 = 2.0, d1 = 0.4, d2 = 0.01, x = 0.76, y = 0.16, z = 9.9)...)

# initialize a 3D plot with 1 empty series
# requires same trick to flip axis Y as for figure2D
plt = plot3d(1, xlabel = "x", xlim = (0.0, 1.0), xticks = 0.0:0.5:1.0, xrotation = 30.0,
                ylabel = "y", ylim = (0.0, 0.5), yticks = (0.0:0.25:0.5, string.(0.5:-0.25:0.0)), yrotation = -10.0,
                zlabel = "z", zlim = (7.5, 10.5), zticks = 7.5:1.5:10.5, zrotation = 0.0,
                linewidth = 0.1,
                marker = :none, legend = :none, size = (600, 600), c = :black,
                right_margin = 5Plots.mm, left_margin = 2Plots.mm,
                )

# build an animated gif by pushing new points to the plot, saving every 20th frame
anim = @animate for i = 1:30000
    step!(TeaCup_sol)
    push!(plt, TeaCup_sol.x, (0.5 - TeaCup_sol.y), TeaCup_sol.z)
end every 20
gif(anim, joinpath("..", "article", "figures", "figure2D.gif"), fps = 40)


text:
@with_kw struct Cell
    boundary::BitVector
    interior::BitVector
    size::Int64 = sum(boundary .| interior)
end

"""
    separator!(g, subgraph_nodes, max_imbalance=0.6, nsample=20)
Call flowcutter (with s and t random) on g and pick nodes at random in the cut matching the 60% imbalance
heuristic to form a separator. Then split the graph g with respect to the separator.
Parameter max_imbalance define the heurestic to choose a cut amoung all generated by
flowcutter. Parameter max_nsample define the number of call to flowcutter.

# Arguments
- `g::SimpleGraph{Int64}`: the subgraph to run flowcutter.
- `subgraph_nodes::Vector{Int64}`: a table where element at index i is corresponding label in the original root graph.
- `max_imbalance::Float64` max imbalance to select cut from flowcutter result.
- `nsample::Int64` max number of try to get a separator.
- `seed::Int64` base seed for RNG choosing input for flowcutter.
# Return
- `sep::Vector{Int64}`: the set of indices of the nodes in separator.
- `split_parts::Array{Array{Int64,1},1}` several parts of the graph resulted from separation.
"""
function separator!(
    g::SimpleGraph{Int64},
    subgraph_nodes::Vector{Int64},
    max_imbalance::Float64,
    nsample::Int64,
    seed::Int64,
)::Tuple{Vector{Int64},Vector{Vector{Int}}}

    # run flowcutter many times and collect all of these cuts
    n = length(subgraph_nodes)
    local_cuts = [Cut[] for _ = 1:Threads.nthreads()]
    rngs = [MersenneTwister(seed + i) for i = 1:Threads.nthreads()]
    Threads.@threads for i = 1:nsample
        s, t = sample(rngs[Threads.threadid()], 1:n, 2, replace=false)
        @inbounds append!(local_cuts[Threads.threadid()], flowcutter(g, s, t))
    end
    cuts = collect(Iterators.flatten(local_cuts))
    filter!(c -> c.imbalance < max_imbalance, cuts)
    # remove dominated cuts
    candidates = Dict{Int64, Cut}()
    for c in cuts
        size = c.size
        # update candidates if they are not previous candidates with its size
        # or in case the candidate is better than its predecessor i.e same size but lower imbalance
        if !(size in keys(candidates)) || c.imbalance < candidates[size].imbalance
            candidates[size] = c
        end
    end

    # select cut with min expansion
    cuts = [c for c in values(candidates)]
    cut = cuts[findmin(map(c->c.expansion, cuts))[2]]
    # select the nodes from the larger side of the cut
    sep = unique(map(a -> a.first, cut.arcs))

    # split the subgraph in several parts
    # be careful with index while removing vertices from sep
    labels = collect(1:n)
    @inbounds for rm_origin_id in sep
        rm = findfirst(id -> id == rm_origin_id, labels)
        rem_vertex!(g, rm)
        labels[rm] = labels[end]
        pop!(labels)
    end
    # catch the multiple parts created from the split
    split_parts = connected_components(g)
    # be careful to return indices which make sense in the original root graph
    return  @inbounds (
        map(node -> subgraph_nodes[node], sep),
        map(
            vector -> map(node -> subgraph_nodes[labels[node]], vector),
            split_parts,
        ),
    )
end

"""
    iterative_dissection(g, best_tw = typemax(Int64), max_imbalance = 0.6, nsample = 20, seed = 42,)
Computes an approximation of the upper bound of the treewidth of the graph g and
an order for elimination using the iterative dissection and the flow cutter algorithm.

# Arguments
- `g::SimpleGraph{Int64}`: the graph to analyse.
- `best_tw::Int64 = typemax(Int64)` stop when hitting a tw bigger than this.
- `max_imbalance::Float64 = 0.6` criteria for selecting cuts that will build the separator.
- `nsample::Int64 = 20` the number of calls to flowcutter with random inputs.
- `seed::Int64 = 4242` the base seed for the RNG sampling the inputs of flowcutter and the choosen one cut for separator.

# Return
- `order::Vector{Int64}` an array of vertices index.
- `treewidth::Int64` the approximation of treewidth.
"""
function iterative_dissection(
        g::SimpleGraph{Int64},
        best_tw::Int64 = typemax(Int64),
        max_imbalance::Float64 = 0.6,
        nsample::Int64 = 20,
        seed::Int64 = 4242,
    )::Pair{Vector{Int64}, Int64}
    n = nv(g)
    order = zeros(Int64, n)
    treewidth = 0
    q = PriorityQueue{Cell, Int64}(Base.Order.Reverse)
    enqueue!(q, Cell(boundary = falses(n), interior = trues(n)), n)
    i = n
    while (!isempty(q)) && (treewidth < best_tw)
        cell = dequeue!(q)
        subgraph_nodes = findall(it -> it > 0, cell.interior)
        n = length(subgraph_nodes)
        if cell.size == treewidth + 1
            while !isempty(q)
                @inbounds order[(i - n + 1):i] .= subgraph_nodes
                i -= n
                subgraph_nodes = findall(it -> it > 0, dequeue!(q).interior)
                n = length(subgraph_nodes)
            end
            @inbounds order[(i - n + 1):i] .= subgraph_nodes
            i -= n
            break
        end

        graph = induced_subgraph(g, subgraph_nodes)[1]
        # if graph is a tree or complete we can stop

        nedges = ne(graph)
        if nedges == n * (n - 1) / 2
            @inbounds order[(i - n + 1):i] .= subgraph_nodes
            i -= n
            treewidth = max(cell.size - 1, treewidth)
            # @debug "bagsize = $(cell.size) tw = $treewidth"
            continue
        elseif nedges == n - 1
            @inbounds order[(i - n + 1):i] .= tree_order!(graph, subgraph_nodes)
            i -= n
            treewidth = max(1 + sum(cell.boundary), treewidth)
            continue
        end

        # compute separator and cut graph in several parts (graph and indices)
        sep, toqueue = separator!(graph, subgraph_nodes, max_imbalance, nsample, seed)

        # update order and treewidth
        k = length(sep)
        @inbounds order[(i - k + 1):i] = sep
        i -= k
        treewidth = max(k + sum(cell.boundary) - 1, treewidth)
        # @debug "bagsize = $(k + sum(cell.boundary)) tw = $treewidth"
        n = nv(g)
        # add next subgraphs to the queue
        @inbounds for new_interior in toqueue
            Ic = falses(n)
            Ic[new_interior] .= true
            Bc = falses(n)
            Bc[sep] .= true
            Bound_union_Sep = findall(Bc .| cell.boundary)
            Bc[sep] .= false
            tmp_Bc = Bound_union_Sep[findall(
                node -> any(
                    node2 -> has_edge(g, node, node2),
                    new_interior
                    ),
                Bound_union_Sep
            )]
            Bc[tmp_Bc] .= true
            new_cell = Cell(boundary = Bc, interior = Ic)
            enqueue!(q, new_cell, new_cell.size)
        end
    end
    return order => treewidth
end


text:
using Documenter, MolecularGraph

makedocs(
    sitename="MolecularGraph.jl",
    pages = [
        "Home" => "index.md",
        "MolecularGraph" => [
            "I/O" => "moleculargraph/io.md",
            "Structure drawing" => "moleculargraph/draw.md",
            "Chemical properties" => "moleculargraph/properties.md",
            "Preprocessing" => "moleculargraph/preprocess.md",
            "Substructure match" => "moleculargraph/substructure.md",
            "MCS" => "moleculargraph/mcs.md",
            "Functional group detection" => "moleculargraph/funcgroup.md"
        ],
        "GraphModel" => [
            "Interface" => "graph/interface.md",
            "Generator" => "graph/generator.md",
            "Traversal" => "graph/traversal.md",
            "Topology" => "graph/topology.md",
            "Operations" => "graph/operation.md",
            "Clique" => "graph/clique.md",
            "Isomorphism" => "graph/isomorphism.md"
        ],
        "Python interface" => "python.md",
        "Design of molecular graph models" => "design.md"
    ]
)

deploydocs(repo="github.com/mojaie/MolecularGraph.jl.git")


text:
##
#
# Draw various histograms
#

using ChristensenSet

DATA_FILE = "ex1i_Float64.dat"

# set up images
images = (
  (RootsImage(-2.0-2.0im, 2.0+2.0im, 2000, 2000), "img1i_1.png", "img1i_1.tex"),
  (RootsImage(-1.0-1.0im, 0.0+0.0im, 2000, 2000), "img1i_2.png", "img1i_2.tex"),
  (RootsImage(-1.5-1.5im, -0.5-0.5im, 2000, 2000), "img1i_3.png", "img1i_3.tex"),
)

# fill images with roots
fill_images!(collect(map(first, images)), DATA_FILE)

# save images to disk
for img in images
  save_image(img[1], img[2], latex_filename=img[3])
  run(`pdflatex $(img[3])`)
end

@info "Done."


text:
# Initially adapted from PiecewiseAffineTransforms.jl
# https://github.com/dfdx/PiecewiseAffineTransforms.jl

global MESHWARP_POLY2SOURCE_MASK = zeros(Bool, 0, 0);

"""
`MESHWARP` - Apply piecewise affine transform to image using bilinear interpolation

    warped_img, [bb.i, bb.j] = meshwarp(img, src, dst, trigs, offset, interp)

* `img`: 2D array, image (todo: extend to Image type)
* `src`: 2xN array, positions of mesh nodes in `img` (global space)
* `dst`: 2xN array, positions of mesh nodes in `warped_img` (global space)
* `trigs`: 3xN array, each column contains indices of nodes of a triangle
* `offset`: 2-element array, position of img[1,1] in global space (optional - default is [0,0])
* `interp`: bool determining whether to use bilinear interpolation or not
    (optional - default is true)
* `warped_img`: with pixel values the same type as original image (for Int type, pixel values are rounded)
* `warped_offset`: 2-element array, position of warped_img[1,1] in global space

The bounding box of the warped image is defined as the smallest
integer-valued rectangle that contains the `dst` mesh.

This means that `warped_offset` is constrained to be integer-valued,
though `offset` is allowed to have floating point values.  The integer
constraint removes the need for further interpolation in any
subsequent fusing of multiple warped tiles.

See definitions in IMWARP documentation for further help.

"""
function meshwarp{T}(img::SharedArray{T},
                    src::Matrix{Float64}, dst::Matrix{Float64},
                    trigs::Matrix{Int64}, offset=[0,0], interp=true)

  bb = snap_bb(find_mesh_bb(dst))
  warped_img = SharedArray{T}(bb.h+1, bb.w+1)
  warped_offset = [bb.i, bb.j];

  Us = Array{Array{Float64, 1}}(size(trigs, 2));
  Vs = Array{Array{Float64, 1}}(size(trigs, 2));
  Ms = Array{Array{Float64, 2}}(size(trigs, 2));

  @everywhere gc();

  src_tri = ones(Float64, 3, 3)
  dst_tri = ones(Float64, 3, 3)

  @fastmath @inbounds for t in 1:size(trigs, 2);
    #println("$t / $(size(trigs, 1))")
    #tr = squeeze(trigs[t, :], 1) = necessary for 0.4.6
    tr = view(trigs, :, t)
    #= slow version
    # coordinates of the source triangle
    X = src[tr, 1]
    Y = src[tr, 2]
    # coordinates of the destination triangle
    U = dst[tr, 1]
    V = dst[tr, 2]
    # Create matrix to transform from destination triangle to source image
    src_tri = [X Y ones(3,1)]
    dst_tri = [U V ones(3,1)]
    =#
    U = dst[1, tr]
    V = dst[2, tr]
    # Create matrix to transform from destination triangle to source image
    src_tri[:, 1] = view(src, 1, tr)
    src_tri[:, 2] = view(src, 2, tr)

    dst_tri[:, 1] = U
    dst_tri[:, 2] = V

    M = dst_tri \ src_tri # dst_tri * M = src_tri
    # Create list of coordinates in warped image that represent this
    # triangle (coordinates in the global space).
    Us[t] = U;
    Vs[t] = V;
    Ms[t] = M;
  end

  function proc_range(idx, arr::Array)
    worker_procs = setdiff(procs(), myid());
    nchunks = length(worker_procs);
    if nchunks == 0 return 1:length(arr); end
    if idx == myid() return 1:0; end
    splits = [round(Int64, s) for s in linspace(0, length(arr), nchunks + 1)];
    return splits[findfirst(worker_procs, idx)]+1:splits[findfirst(worker_procs, idx) + 1]
  end

  #=@sync @fastmath @inbounds for t in 1:size(trigs, 1);
      @async remotecall_wait(procs()[rem(t, nprocs()-1)+2], calculate_pixels_in_trig!, Us[t], Vs[t], Ms[t], img, offset, warped_img, warped_offset)
  end=#

  @sync for p in procs()
          t = proc_range(p, Us);
          @async @fastmath @inbounds remotecall_wait(calculate_pixels_in_trig_chunk!, p, Us[t], Vs[t], Ms[t], img, offset, warped_img, warped_offset, interp ? computepixel_interp : computepixel)
  end

  return warped_img, [bb.i, bb.j]
end

function calculate_pixels_in_trig_chunk!(Us, Vs, Ms, img, offset, warped_img, warped_offset, interp)
 @simd for i in 1:length(Us)
	calculate_pixels_in_trig!(Us[i], Vs[i], Ms[i], img, offset, warped_img, warped_offset, interp)
  end
  global MESHWARP_POLY2SOURCE_MASK = zeros(Bool, 0, 0);
end

function calculate_pixels_in_trig!(U, V, M, img, offset, warped_img, warped_offset, computepixel_func)
  poly::Array{Bool,2}, i_oset::Int64, j_oset::Int64 = poly2source_new!(U, V)
    for j_poly in 1:size(poly,2)
    	for i_poly in 1:size(poly,1)
    @fastmath @inbounds begin
	#check if poly is true
      	if poly[i_poly, j_poly] == false continue end
        # Use warped coordinate in global space for transform - this changes back from the test polygon to global space
	u, v = i_poly + i_oset, j_poly + j_oset
        # Convert warped coordinate to pixel space
        i, j = u-warped_offset[1]+1, v-warped_offset[2]+1
        # x, y = M * [u, v, 1]
        @fastmath @inbounds x, y = M[1,1]*u + M[2,1]*v + M[3,1], M[1,2]*u + M[2,2]*v + M[3,2]
        # Convert original image coordinate to pixel space
        x, y = x-offset[1]+1, y-offset[2]+1
        fx, fy = floor(Int64, x), floor(Int64, y)
        wx, wy = x-fx, y-fy
        if 1 <= fx <= size(img, 1)-1 && 1 <= fy <= size(img, 2)-1
	computepixel_func(warped_img, img, i, j, wx, wy, fx, fy)
        end
      end
      end
    end
#=
   us, vs = poly2source!(U, V)
#    computepixel_func = interp ? computepixel_interp : computepixel
    @simd for ind in 1:length(us)
        # Convert warped coordinate to pixel space
    @fastmath @inbounds begin
        i, j = us[ind]-warped_offset[1]+1, vs[ind]-warped_offset[2]+1
        # Use warped coordinate in global space for transform
        u, v = us[ind], vs[ind]
        # x, y = M * [u, v, 1]
        @fastmath @inbounds x, y = M[1,1]*u + M[2,1]*v + M[3,1], M[1,2]*u + M[2,2]*v + M[3,2]
        # Convert original image coordinate to pixel space
        x, y = x-offset[1]+1, y-offset[2]+1
        fx, fy = floor(Int64, x), floor(Int64, y)
        wx, wy = x-fx, y-fy
        if 1 <= fx <= size(img, 1)-1 && 1 <= fy <= size(img, 2)-1
	computepixel_func(warped_img, img, i, j, wx, wy, fx, fy)
        end
      end
    end
    =#
end

    # Takes weights wx, wy that denotes how close the value is to [fx+1, fy+1] from [fx, fy] and writes the weighted average to [i, j]
    @inline function computepixel_interp(warped_img, img, i, j, wx, wy, fx, fy)
          # Expansion of p = [1-wy wy] * img[fy:fy+1, fx:fx+1] * [1-wx; wx]
          @fastmath @inbounds p1 = ((1-wy)*img[fx,fy] + wy*img[fx,fy+1])
	  @fastmath @inbounds p2 = ((1-wy)*img[fx+1,fy] + wy*img[fx+1,fy+1])
	  @fastmath p1 = p1 * (1-wx);
	  @fastmath p2 = p2 * (wx);
          writepixel(warped_img,i,j,p1+p2)
    end

    function computepixel(warped_img, img, i, j, wx, wy, fx, fy)
      	  @fastmath y = wy < 0.5 ? fy : fy + 1
      	  @fastmath x = wx < 0.5 ? fx : fx + 1
	  @inbounds v = img[x,y]
          writepixel(warped_img,i,j,v)
    end

"""
`POLY2SOURCE` - Run fillpoly and findn over the basic bounding box of a given triangle
"""
function poly2source(pts_i, pts_j)
  # Find bb of vertices (vertices in global space)
  top, bottom = floor(Int64,minimum(pts_i)), ceil(Int64,maximum(pts_i))
  left, right = floor(Int64,minimum(pts_j)), ceil(Int64,maximum(pts_j))
  # Create image based on number of pixels in bb that will identify triangle
  mask = zeros(Bool, bottom-top+1, right-left+1)
  # Convert vertices into pixel space and fill the mask to identify triangle
  fillpoly!(mask, pts_j-left+1, pts_i-top+1, true; convex = false)
  # Create list of pixel coordinates that are contained by the triangle
  us, vs = findn(mask)
  # Convert that list of pixel coordinates back into global space
  us += top-1
  vs += left-1
  return us, vs
end

@inline function poly2source_new!(pts_i, pts_j)
  # Find bb of vertices (vertices in global space)
  top, bottom = floor(Int64,minimum(pts_i)), ceil(Int64,maximum(pts_i))
  left, right = floor(Int64,minimum(pts_j)), ceil(Int64,maximum(pts_j))
  # Create image based on number of pixels in bb that will identify triangle
  if size(MESHWARP_POLY2SOURCE_MASK::Array{Bool, 2})[1] < bottom - top + 1 || size(MESHWARP_POLY2SOURCE_MASK)[2] < right - left + 1
    global MESHWARP_POLY2SOURCE_MASK = zeros(Bool, bottom-top+1, right-left+1)
  end
  (MESHWARP_POLY2SOURCE_MASK::Array{Bool, 2})[:] = false;

  #mask = zeros(Bool, bottom-top+1, right-left+1)
  # Convert vertices into pixel space and fill the mask to identify triangle
  for i in 1:length(pts_i)
    @fastmath @inbounds pts_i[i] = pts_i[i] + 1 - top;
  end
  for i in 1:length(pts_j)
    @fastmath @inbounds pts_j[i] = pts_j[i] + 1 - left;
  end
  #fillpoly!(mask, pts_j-left+1, pts_i-top+1, true)
   fillpoly_convex!(MESHWARP_POLY2SOURCE_MASK::Array{Bool, 2}, pts_j, pts_i, true)
   # return the offsets that have to be added
   return MESHWARP_POLY2SOURCE_MASK, top - 1, left - 1
   #=
  # Create list of pixel coordinates that are contained by the triangle
  us, vs = findn(MESHWARP_POLY2SOURCE_MASK::Array{Bool, 2})
  # Convert that list of pixel coordinates back into global space
  for i in 1:length(us)
    @fastmath @inbounds us[i] = us[i] - 1 + top;
  end
  for i in 1:length(vs)
    @fastmath @inbounds vs[i] = vs[i] - 1 + left;
  end
#  us += top-1
#  vs += left-1
  return us, vs
  =#
end

function poly2source!(pts_i, pts_j)
  # Find bb of vertices (vertices in global space)
  top, bottom = floor(Int64,minimum(pts_i)), ceil(Int64,maximum(pts_i))
  left, right = floor(Int64,minimum(pts_j)), ceil(Int64,maximum(pts_j))
  # Create image based on number of pixels in bb that will identify triangle
  if size(MESHWARP_POLY2SOURCE_MASK::Array{Bool, 2})[1] < bottom - top + 1 || size(MESHWARP_POLY2SOURCE_MASK)[2] < right - left + 1
    global MESHWARP_POLY2SOURCE_MASK = zeros(Bool, bottom-top+1, right-left+1)
  end
  (MESHWARP_POLY2SOURCE_MASK::Array{Bool, 2})[:] = false;
  #mask = zeros(Bool, bottom-top+1, right-left+1)
  # Convert vertices into pixel space and fill the mask to identify triangle
  for i in 1:length(pts_i)
    @fastmath @inbounds pts_i[i] = pts_i[i] + 1 - top;
  end
  for i in 1:length(pts_j)
    @fastmath @inbounds pts_j[i] = pts_j[i] + 1 - left;
  end
  #fillpoly!(mask, pts_j-left+1, pts_i-top+1, true)
   fillpoly_convex!(MESHWARP_POLY2SOURCE_MASK::Array{Bool, 2}, pts_j, pts_i, true)
  # Create list of pixel coordinates that are contained by the triangle
  us, vs = findn(MESHWARP_POLY2SOURCE_MASK::Array{Bool, 2})
  # Convert that list of pixel coordinates back into global space
  for i in 1:length(us)
    @fastmath @inbounds us[i] = us[i] - 1 + top;
  end
  for i in 1:length(vs)
    @fastmath @inbounds vs[i] = vs[i] - 1 + left;
  end
#  us += top-1
#  vs += left-1
  return us, vs
end

"""
`FILLPOLY!` - Fill pixels contained inside a polygon

    M = fillpoly!(M, px, py, value; reverse)

* `M`: 2D array, image
* `px`: 1D array, x-components of polygon vertices
* `py`: 1D array, y-components of polygon vertices
* `value`: fill value
* `convex`: Boolean kwarg, false by default; if true then assumes the polygon is convex and well-formed without duplicates
* `reverse`: Boolean kwarg, false by default; if true then fills in everything outside the polygon

Features:

* Intended to work for nonconvex polygons as well as convex polygons
* The vertices should be listed sequentially in px and py. Clockwise/counterclockwise ordering doesn't matter.
* Duplicate points in sequence are handled by removing until only one remains.
* Polygons need not be closed (i.e. last point needs not be the same as the first point) - if closed, they will be handled as duplicates initially.
* The code treats the "corner case," where the vertex of the polygon intersects the grid line at one point.
* The code treats the "edge case," where an edge of the polygon lies exactly on a grid line.
* The boundary of the polygon is filled whether or not reverse kwarg is set - i.e. the results with `reverse` set to true and false are not perfect complements due to the boundary.

Bugs:

* The original code (https://github.com/dfdx/PiecewiseAffineTransforms.jl/blob/master/src/polyline.jl) used implicit conversion to Int64 (presumably rounding).  Tommy/Shang replaced this by ceil.  This might produce inconsistent results, with the same pixel belonging to more than one triangle.
"""
function fillpoly!{T,P<:Number}(M::Matrix{T}, px::Vector{P}, py::Vector{P}, value::T; convex::Bool=false, reverse::Bool=false)
  convex ? (return fillpoly_convex!(M, px, py, value; reverse = reverse)) : (return fillpoly_nonconvex!(M, px, py, value; reverse = reverse))
end

function fillpoly!{T,P<:Number}(M::SharedArray{T}, px::Vector{P}, py::Vector{P}, value::T; convex::Bool=false, reverse::Bool=false)
  M = convert(Array, M)
  fillpoly!(M, px, py, value; convex=convex, reverse=reverse)
  M = convert(SharedArray, M)
end

function fillpoly_convex!{T,P<:Number}(M::Matrix{T}, px::Vector{P}, py::Vector{P}, value::T; reverse::Bool=false)
  left, right = floor(Int64,minimum(px)), ceil(Int64,maximum(px))
  if reverse xrange = 1:size(M, 2)
  else xrange = left:right
  end
  l = length(px)
  # Scan poly from left to right
  for x=xrange     # loop over grid lines
    top = 0
    bot = 0
    m = length(px)
    xdir_last = sign(px[m] - px[m-1]) # direction of the line segment
    for n=1:length(px)  # loop over edges (m,1),(1,2),(2,3),...,(m-1,m)
      xdir_cur = sign(px[n] - px[m])
      # grid line intersects edge in one of two ways
      if (px[n] <= x <= px[m]) || (px[m] <= x <= px[n])
	if px[n] == px[m]  # intersection is entire edge, do nothing
	else# intersection is point
	  # do not add duplicate points (endpoint of the last segment), unless the direction is reversing in x or the last segment is vertical
	  if px[m] != x || xdir_last + xdir_cur == 0 || xdir_last == 0
            y = py[n] + (x-px[n]) * (py[m]-py[n])/(px[m]-px[n])
	    # deal with rounding error
	    if ceil(Int64, y) > size(M, 1)
	      val = floor(Int64, y)
	    else
              val = ceil(Int64, y)
	    end
	    if top == 0
	      top = val
	    elseif top < val
	      bot = val
	    else
	      bot = top
	      top = val
	      break
	    end
       	  end
      	end
      end
      xdir_last = xdir_cur
      m = n
    end

    if top == 0 || bot == 0 continue end

    if reverse
      @simd for y in 1:top
      @inbounds M[y, x] = value
      end
      @simd for y in bot:size(M,1)
      @inbounds M[y, x] = value
      end
    else
    # Place value in matrix at all the y's between min and max for given x
      @simd for y in top:bot
      @inbounds M[y, x] = value
      end
    end
  end
  return M
end

function fillpoly_nonconvex!{T,P<:Number}(M::Matrix{T}, px::Vector{P}, py::Vector{P}, value::T; reverse::Bool=false)
  @assert length(py) == length(px)
  left, right = floor(Int64,minimum(px)), ceil(Int64,maximum(px))
  if reverse xrange = 1:size(M, 2)
  else xrange = left:right
  end
  l = length(px)
  #force no duplicates, including at the end
    for i in 1:length(px)
	if px[i] == px[l] && py[i] == py[l]
	  px[l] = 0
	  py[l] = 0
	end
	  l = i
    end
    px = px[px .!= 0]
    py = py[py .!= 0]
  # Scan poly from left to right
  for x=xrange     # loop over grid lines
    ys = Array{Int64, 1}()
    tops = Set{Int64}() # tops of vertical edges
    signs = Array{Int64, 1}()
    ys_clean = Array{Int64, 1}()
    m = length(px)
    xdir_last = sign(px[m] - px[m-1]) # direction of the line segment
    for n=1:length(px)  # loop over edges (m,1),(1,2),(2,3),...,(m-1,m)
      xdir_cur = sign(px[n] - px[m])
      # grid line intersects edge in one of two ways
      if (px[n] <= x <= px[m]) || (px[m] <= x <= px[n])
	if px[n] == px[m]  # intersection is entire edge
#	  if xdir_last + xdir_cur == 0 # if two vertical edges in a row remove the last one
#	    pop!(ys)
#	  end
	  push!(tops, ceil(Int64, minimum((py[n], py[m]))))
	  #push!(signs, xdir_cur)
          #push!(ys, ceil(Int64, py[n]))
	else# intersection is point
	  # do not add duplicate points (endpoint of the last segment), unless the direction is reversing in x or the last segment is vertical
	  if px[m] != x || xdir_last + xdir_cur == 0 || xdir_last == 0
            y = py[n] + (x-px[n]) * (py[m]-py[n])/(px[m]-px[n])
	    push!(signs, xdir_cur)
	    # deal with rounding error
	    if ceil(Int64, y) > size(M, 1)
	      push!(ys, floor(Int64, y))
	    else
              push!(ys, ceil(Int64, y))
	    end
       	  end
      	end
      end
      xdir_last = xdir_cur
      m = n
    end

    if length(ys) != 2
    # generically, two intersections for a convex polygon
    perm = sortperm(ys)  # sort the intersection points
    ys = ys[perm]  # sort the intersection points
    signs = signs[perm]  # sort the intersection points' signs
    i = 1
    while (i < length(ys))
	y_entry = ys[i]
	y_exit = ys[i]
	entry_sign = signs[i]
	while (i < length(ys))
	  i += 1
	  y_exit = ys[i]
	  exit_sign = signs[i]
	  if in(y_exit, tops) continue;
	  elseif exit_sign != entry_sign break; end
	end
	push!(ys_clean, y_entry)
	push!(ys_clean, y_exit)
	i += 1
    end
  else
    #if ys length is two then we can simply sort instead of going through the whole algorithm
    ys_clean = sort(ys)
    end

    if reverse
      unshift!(ys_clean, 1)
      push!(ys_clean, size(M,1))
    end
    # Place value in matrix at all the y's between min and max for given x
    for n=1:2:length(ys_clean)
      @simd for y in ys_clean[n]:ys_clean[n+1]
      @inbounds M[y, x] = value
      end
    end
  end
  return M
end

function meshwarp{T}(img::Array{T},
                    src::Matrix{Float64}, dst::Matrix{Float64},
                    trigs::Matrix{Int64}, offset=[0,0], interp=true)
  img = convert(SharedArray, img)
  meshwarp(img, src, dst, trigs, offset, interp)
end


text:
module MLP

export PickModel
export chain, prepare, train!, show_accuracy
export save, load

using Flux
using Flux: onehotbatch, argmax, crossentropy, throttle, @epochs
using BSON: @save, @load
using Base.Iterators: repeated
using MLDatasets # FashionMNIST
using ColorTypes: N0f8, Gray

const Img = Matrix{Gray{N0f8}}

struct PickModel
    title::Symbol
    pick::Vector{Int}
    m::Chain
end

function prepare(data_x, data_y, pick::Vector{Int})
    # 18_000
    inds = find(x->x in pick, data_y)

    imgs = Img.([data_x[:,:,i] for i in inds])
    labels = data_y[inds]

    X = hcat(float.(reshape.(imgs, :))...) |> gpu
    Y = onehotbatch(labels, pick) |> gpu
    X, Y
end

function chain(output::Int)
    Chain(
        Dense(28^2, 32, relu),
        Dense(32, output),
        softmax) |> gpu
end

function train!(title, m::Chain, pick::Vector{Int}, X, Y; n=5)
    loss(x, y) = crossentropy(m(x), y)
    dataset = repeated((X, Y), 20) # 200
    evalcb = () -> @show(loss(X, Y))
    opt = ADAM(params(m))

    @epochs n Flux.train!(loss, dataset, opt, cb = throttle(evalcb, 3))
    PickModel(title, pick, m)
end

show_accuracy(m::Chain, pick::Vector{Int}, x, y) = @show pick mean(argmax(m(x)) .== argmax(y))

model_name(title::Symbol, pick::Vector{Int}, target) =
    string(title, "_mlp_", join(pick), "_", target, ".bson")

function save(p::PickModel)
    name = model_name(p.title, p.pick, "model")
    m = p.m
    @save name m
    weights = Tracker.data.(params(m)) ;
    name = model_name(p.title, p.pick, "weights")
    @save name weights
end

function load(title::Symbol, pick::Vector{Int})
    name = model_name(title, pick, "model")
    @load name m
    name = model_name(title, pick, "weights")
    @load name weights
    Flux.loadparams!(m, weights)
    PickModel(title, pick, m)
end

end # module MLP


text:
using DataFrames, QuasiCopula, LinearAlgebra, GLM, RDatasets
using Test

Mmmec = dataset("mlmRev", "Mmmec");
df = Mmmec
y = :Deaths
grouping = :Region
covariates = [:UVB]

d = Poisson()
link = LogLink()

#### POISSON ####
# forming AR(1) model with Poisson base
Poisson_AR_model = AR_model(df, y, grouping, covariates, d, link)

@test typeof(Poisson_AR_model) == GLMCopulaARModel{Float64, Poisson{Float64}, LogLink}
@test Poisson_AR_model.penalized == false

# forming AR(1) model with Poisson base with L2 ridge penalty
Poisson_AR_model_penalized = AR_model(df, y, grouping, covariates, d, link; penalized = true)

@test typeof(Poisson_AR_model_penalized) == GLMCopulaARModel{Float64, Poisson{Float64}, LogLink}
@test Poisson_AR_model_penalized.penalized == true

# intercept only model
Poisson_AR_model_intercept_only = AR_model(df, y, grouping, d, link)
# check if default to intercept only
@test sum(Poisson_AR_model_intercept_only.data[1].X) == Poisson_AR_model_intercept_only.data[1].n

#### nb ####
d = NegativeBinomial()
link = LogLink()
nb_AR_model = AR_model(df, y, grouping, covariates, d, link)

@test typeof(nb_AR_model) == NBCopulaARModel{Float64, NegativeBinomial{Float64}, LogLink}
@test Poisson_AR_model.penalized == false

# forming AR(1) model with nb base with L2 ridge penalty
nb_AR_model_penalized = AR_model(df, y, grouping, covariates, d, link; penalized = true)

@test typeof(nb_AR_model_penalized) == NBCopulaARModel{Float64, NegativeBinomial{Float64}, LogLink}
@test nb_AR_model_penalized.penalized == true

# intercept only model
nb_AR_model_intercept_only = AR_model(df, y, grouping, d, link)
# check if default to intercept only
@test sum(nb_AR_model_intercept_only.data[1].X) == nb_AR_model_intercept_only.data[1].n

#### gaussian ####
d = Normal()
link = IdentityLink()
gaussian_AR_model = AR_model(df, y, grouping, covariates, d, link)

@test typeof(gaussian_AR_model) == GaussianCopulaARModel{Float64}
@test gaussian_AR_model.penalized == false

# forming AR(1) model with gaussian base with L2 ridge penalty
gaussian_AR_model_penalized = AR_model(df, y, grouping, covariates, d, link; penalized = true)

@test typeof(gaussian_AR_model_penalized) == GaussianCopulaARModel{Float64}
@test gaussian_AR_model_penalized.penalized == true

# intercept only model
gaussian_AR_model_intercept_only = AR_model(df, y, grouping, d, link)
# check if default to intercept only
@test sum(gaussian_AR_model_intercept_only.data[1].X) == gaussian_AR_model_intercept_only.data[1].n


text:
###########################################
# States

"""
    AbstractState{N, W}
Abstract type representing state of the (index) automaton over an alphabet of `N` letters.

The subtypes of `AbstractState{N, W}` should implement the following methods which
constitute `AbstractState` interface:
 * `name`: external name of the state
 * `isterminal`: predicate to determine if state is terminal (i.e. whether the
    state represents left-hand side of some rewriting rule)
 * `rightrule`: if state is terminal, returns the right-hand side of the rewriting rule.
 Undefined in other cases.
 * `inedges`: container (indexed by letters of `AbstractWord`) of source states of in-edges.
 * `outedges`: container (indexed by letters of `AbstractWord`) of destination states of out-edges.
 * `isfailstate`: predicate to determine if state represents "fail" state.
 * `Base.length`: the length of the signature of the shortest path to the
    state, which starts at initial state (length of "simple path")
"""
abstract type AbstractState{N, W} end

"""
    State{N, W<:AbstractWord} <: AbstractState{N, W}
State of an index automaton over an alphabet of `N` letters.
"""
mutable struct State{N, W<:AbstractWord} <: AbstractState{N, W}
    name::W
    terminal::Bool
    failstate::Bool
    rrule::W
    ined::Vector{State{N, W}}
    outed::NTuple{N, State{N, W}}

    function State{N, W}() where {N, W<:AbstractWord}
        x = new{N, W}()
        x.failstate = true
        return x
    end

    function State(name::W, failstate::State{N, W}) where {N, W}
        s = new{N, W}(
            name,
            false,
            false,
            one(W),
            State{N,W}[],
            ntuple(i -> failstate, N),
        )
        return s
    end
end

name(s::State) = s.name
isterminal(s::State) = s.terminal
rightrule(s::State) = s.rrule
inedges(s::State) = s.ined
outedges(s::State) = s.outed
isfailstate(s::State) = s.failstate
Base.length(s::State) = length(name(s))

"""
    declarerightrule!(s::State, w::AbstractWord)
Decalres given word as a right-hand side rule of a given state (and makes this)
state terminal.
"""
function declarerightrule!(s::State, w::AbstractWord)
    s.terminal = true
    s.rrule = w
end

function Base.show(io::IO, s::State)
    if isfailstate(s)
        print(io, "Fail State")
    elseif isterminal(s)
        print(io, "TState($(name(s)) → $(rightrule(s)))")
    else
        print(io, "NTState($(length(inedges(s))), $(length(outedges(s))))")
    end
end

function Base.show(io::IO, ::MIME"text/plain", s::State)
    if isfailstate(s)
        print(io, "Fail state")
    elseif isterminal(s)
        println(io, "Terminal state $(name(s))")
        println(io, "Right rule: $(rightrule(s))")
        println(io, " Edges entering:")
        for (i, e) in enumerate(inedges(s))
            !isfailstate(e) && println(io, "   ", " - with label $(i) from state $(name(e))")
        end
    else
        println(io, "State $(name(s))")
        println(io, " Edges entering:")
        for (i, fromst) in enumerate(inedges(s))
            !isfailstate(fromst) && println(io, "   ", " - with label $(i) from state $(name(fromst))")
        end
        println(io, " Edges leaving:")
        for (i, tost) in enumerate(outedges(s))
            !isfailstate(tost) && println(io, "   ", " - with label $(i) from state $(name(tost))")
        end
    end
end

###########################################
# Automata

"""
    AbstractAutomaton{N, W}
Abstract type representing (index) automaton.

The subtypes of `AbstractAutomaton` should implement the following methods which
constitute `AbstractAutomaton` interface:
 * `alphabet`: return the alphabet of the automaton
 * `states`: return the list of states in the automaton
 * `initialstate`: return the initial state of the automaton
 * `failstate`: return the unique fail state (representing "no edge from source with a given label" situation)
 * `Base.push!`: append a single state to the automaton
 * `addedge!`: add edge between two states
 * `removeedge!`: remove edge between two states
 * `Base.deleteat!`: delete state at given position (together with edges)
 * `walk`: travel through the automaton according to given path and initial state
"""
abstract type AbstractAutomaton{N, W} end

"""
    Automaton{N, W<:AbstractWord} <: AbstractAutomaton{N, W}
Automaton as a vector of states together with alphabet.
"""
struct Automaton{N, W<:AbstractWord} <: AbstractAutomaton{N, W}
    states::Vector{State{N, W}}
    abt::Alphabet
    failstate::State{N, W}
    stateslengths::Vector{Int}
    _past_states::Vector{State{N, W}}
end

# Default type is set to UInt16
function Automaton(abt::Alphabet, W::Type{<:AbstractWord}=Word{UInt16})
    S = State{length(abt), W}
    failstate = S()
    Automaton{length(abt), W}([State(one(W), failstate)], abt, failstate, [0], S[])
end

alphabet(a::Automaton) = a.abt
states(a::Automaton) = a.states
initialstate(a::Automaton) = first(states(a))
failstate(a::Automaton) = a.failstate
stateslengths(a::Automaton) = a.stateslengths

function Base.isempty(a::Automaton)
    σ = initialstate(a)
    length(states(a)) == 1 && isone(name(σ)) && return true
    return false
end

Base.push!(a::Automaton{N, W}, name::W) where {N, W} =
    (push!(states(a), State(name, failstate(a))); push!(stateslengths(a), length(name)); a)
Base.empty!(a::Automaton{N, W}) where{N, W} =
    (empty!(states(a)); empty!(stateslengths(a)); push!(a, one(W)); a)

"""
    replace(k::NTuple{N, T}, val::T, idx::Integer) where {N, T}
Returns a copy of `NTuple` `k` with value at index `idx` changed to `val`.
"""
Base.@propagate_inbounds function replace(k::NTuple{N, T}, val::T, idx::Integer) where {N, T}
    @boundscheck 1 ≤ idx ≤ N || throw(BoundsError(k, idx))
    return @inbounds ntuple( i -> i==idx ? val : k[i], Val(N))
end

"""
    updateoutedges!(src::State, dst::State, label::Integer)
Updates outedges of the state `src` with state `dst` connected through the edge `label`.
"""
Base.@propagate_inbounds function updateoutedges!(src::State{N}, dst::State{N}, label::Integer) where {N}
    @boundscheck 1 ≤ label ≤ N || throw(BoundsError(outedges(src), label))
    @inbounds src.outed = replace(outedges(src), dst, label)
end

"""
    addinedge!(dst::State, src::State)
Add in-edge to `dst` from state `src` to state `dst` storing it in `inedges` of the state `dst`.
"""
addinedge!(dst::State, src::State) = push!(inedges(dst), src)

"""
    addedge!(a::Automaton, label::Integer, src::Integer, dst::Integer)
Adds the edge with a given `label` directed from `src` state to `dst` state.
"""
Base.@propagate_inbounds function addedge!(a::Automaton{N}, label::Integer, src::Integer, dst::Integer) where {N}
    @boundscheck 1 ≤ label ≤ N || throw(BoundsError(outedges(states(a)[src]), label))
    @boundscheck checkbounds(states(a), src)
    @boundscheck checkbounds(states(a), dst)
    @inbounds srcstate = states(a)[src]
    @inbounds dststate = states(a)[dst]
    @inbounds updateoutedges!(srcstate, dststate, label)
    addinedge!(dststate, srcstate)
end

Base.@propagate_inbounds function addedge!(a::Automaton{N}, label::Integer, src::State, dst::State) where {N}
    @boundscheck 1 ≤ label ≤ N || throw(BoundsError(outedges(src), label))
    @assert src in states(a) "Source state is not in automaton!"
    @assert dst in states(a) "Destination state to is not in automaton"
    @inbounds updateoutedges!(src, dst, label)
    addinedge!(dst, src)
end
# Should we also check if the edge corresponding to a given letter does not exist?

"""
    removeedge!(a::Automaton, label::Integer, from::Integer, to::Integer)
Removes the edge with a given `label` `from` one state `to` another.
"""
Base.@propagate_inbounds function removeedge!(a::Automaton{N}, label::Integer, from::Integer, to::Integer) where {N}
    @boundscheck 1 ≤ label ≤ N || throw(BoundsError(outedges(states(a)[from]), label))
    @boundscheck checkbounds(states(a), from)
    @boundscheck checkbounds(states(a), to)
    @inbounds fromstate = states(a)[from]
    @inbounds tostate = states(a)[to]
    @inbounds updateoutedges!(fromstate, failstate(a), label)
    deleteat!(inedges(tostate), findfirst(isequal(fromstate), inedges(tostate)))
end

Base.@propagate_inbounds function Base.deleteat!(a::Automaton, idx::Integer)
    @boundscheck checkbounds(states(a), idx)
    @boundscheck checkbounds(stateslengths(a), idx)
    @inbounds state = states(a)[idx]
    for σ in outedges(state)
        isfailstate(σ) || deleteat!(inedges(σ), findfirst(isequal(state), inedges(σ)))
    end
    for σ in inedges(state)
        updateoutedges!(σ, failstate(a), findfirst(isequal(state), outedges(σ)))
    end
    @inbounds deleteat!(states(a), idx)
    @inbounds deleteat!(stateslengths(a), idx)
end

"""
    walk(a::AbstractAutomaton, signature::AbstractWord[, state=initialstate(a)])
Walk through index automaton according to the path given by the `signature`, starting from `state`.

Return a tuple `(idx, state)` where `idx` is the length of prefix of signature which was successfully traced and
`state` is the final state.
Note that if `idx ≠ length(signature)` there is no path in the automaton corresponding to the full signature.
"""
function walk(a::AbstractAutomaton, signature::AbstractWord, state=initialstate(a))
    idx = 0
    for (i, k) in enumerate(signature)
        next = outedges(state)[k]
        isfailstate(next) ? break : (idx, state) = (i, next)
    end
    return (idx, state)
end

function Base.show(io::IO, a::AbstractAutomaton)
    println(io, "Automaton with $(length(states(a))) states")
    for (i, state) in enumerate(states(a))
        print(io, " $i. Edges leaving state (")
        print_repr(io, name(state), alphabet(a))
        println(io, "):")

        for (i, tostate) in enumerate(outedges(state))
            if !isfailstate(tostate)
                print(io,
                "   ", " - with label ", alphabet(a)[i], " to state (")
                print_repr(io, name(tostate), alphabet(a))
                println(io, ")")
            end
        end
    end
end

###########################################
# Ad index automata

function makeindexautomaton!(a::Automaton, rws::RewritingSystem)
    empty!(a)
    # Determining simple paths
    for (lhs, rhs) in rules(rws)
        σ = initialstate(a)
        for (i, letter) in enumerate(lhs)
            if isfailstate(outedges(σ)[letter])
                push!(a, lhs[1:i])
                addedge!(a, letter, σ, last(states(a)))
            end
            σ = outedges(σ)[letter]
        end
        terminal = last(states(a))
        declarerightrule!(terminal, rhs)
        # Add loops for terminal state
        for i in 1:length(outedges(terminal))
            addedge!(a, i, terminal, terminal)
        end
    end
    # Determining cross paths
    for (i, state) in enumerate(outedges(initialstate(a)))
        isfailstate(state) && addedge!(a, i, 1, 1)
    end
    i = 1
    indcs = findall(isequal(i), stateslengths(a))
    while !isempty(indcs)
        for idx in indcs
            σ = states(a)[idx]
            τ = last(walk(a, @view name(σ)[2:end]))
            @assert !isfailstate(τ) "Tracing $σ in an automaton failed!"
            for (letter, state) in enumerate(outedges(σ))
                isfailstate(state) && addedge!(a, letter, σ, outedges(τ)[letter])
            end
        end
        i += 1
        resize!(indcs, 0)
        for (k,len) in enumerate(stateslengths(a))
            len == i && push!(indcs, k)
        end
    end
    return a
end

"""
    makeindexautomaton(rws::RewritingSystem)
Creates index automaton corresponding to a given rewriting system.
"""
makeindexautomaton(rws::RewritingSystem) =
    (a = Automaton(alphabet(rws)); makeindexautomaton!(a, rws))

"""
    updateautomaton!(a::Automaton, rws::RewritingSystem)
Updates given automaton so it corresponds to a given rewriting system.
"""
updateautomaton!(a::Automaton, rws::RewritingSystem) = makeindexautomaton!(a, rws)

"""
    rewrite_from_left!(v::AbstractWord, w::AbstractWord, a::Automaton)
Rewrites word `w` from left using index automaton `a` and appends the result
to `v`. For standard rewriting `v` should be empty.
"""
function rewrite_from_left!(v::AbstractWord, w::AbstractWord, a::Automaton)
    tape = a._past_states
    resize!(tape, length(w) + 1)
    tape[1] = initialstate(a)
    initial_length = length(v)
    @debug "Rewriting $w with automaton"
    while !isone(w)
        x = popfirst!(w)
        k = length(v) - initial_length + 1
        @debug "" v x w
        s = tape[k]
        state = outedges(s)[x]

        if isterminal(state)
            lhs, rhs = name(state), rightrule(state)
            @debug "applying rewriting rule $lhs → $rhs"
            w = prepend!(w, rhs)
            resize!(v, length(v) - length(lhs) + 1)
        else
            length(tape) <= k ? push!(tape, state) : tape[k+1] = state
            push!(v, x)
            @debug "state is not terminal, pushed $x, v=$v"
        end
    end
    return v
end


text:
# Can we approximate power-2 limb-darkening with polynomial?

using PyPlot
include("/Users/ericagol/Computer/Julia/regress.jl")

alpha = [0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,0.99,1.01,1.05,1.1,1.15,1.2,1.25,1.3,1.35,1.4,1.45,1.5,1.55,1.6,1.65,1.7,1.75,1.8]
#nmu = 1000; mu = logspace(-big(2.5),big(0.0),nmu)
nalpha = length(alpha)
nmu = 1000; mu = linspace(1e-3,1.0,nmu)
nord = 30
coeff = zeros(nalpha,nord,nord-1)
scat = zeros(nalpha,nord)
for ialpha=1:nalpha
  fnon = mu.^alpha[ialpha]-mu
  for i=2:nord
    fn = zeros(i-1,nmu)
#    fn[1,:] .= big(1.0)
#    fn[1,:] .= 1.0
#    fn[2,:] .= mu
    for j=2:i
      fn[j-1,:] .= (j+2)*mu.^j-j*mu.^(j-2)-2mu
    end
    ci,cov = regress(fn,fnon,ones(Float64,nmu))
    fmod = zeros(Float64,nmu)
    for j=1:i-1
      fmod += ci[j]*fn[j,:]
    end
    clf()
    plot(mu,convert(Array{Float64,1},fnon))
    plot(mu,convert(Array{Float64,1},fmod))
    plot(mu,convert(Array{Float64,1},fmod-fnon))
    scat[ialpha,i]=convert(Float64,std(fmod-fnon))
    coeff[ialpha,i,1:i-1]=convert(Array{Float64,1},ci)
    println("alpha: ",convert(Float64,alpha[ialpha])," i: ",i," scatter: ",scat[ialpha,i]," max: ",maximum(abs,convert(Array{Float64,1},fmod-fnon)))
  end
end


text:
module TestDualSimplicialSets
using Test

using LinearAlgebra: Diagonal
using SparseArrays, StaticArrays

using Catlab.CategoricalAlgebra.CSets
using CombinatorialSpaces

const Point2D = SVector{2,Float64}
const Point3D = SVector{3,Float64}

# 1D dual complex
#################

primal_s = DeltaSet1D()
add_vertices!(primal_s, 5)
add_edges!(primal_s, 1:4, repeat([5], 4))
s = DeltaDualComplex1D(primal_s)
@test nparts(s, :DualV) == nv(primal_s) + ne(primal_s)
@test nparts(s, :DualE) == 2 * ne(primal_s)

dual_v = elementary_duals(1,s,4)
@test dual_v == [edge_center(s, 4)]
@test elementary_duals(s, E(4)) == DualV(dual_v)

dual_es = elementary_duals(0,s,5)
@test length(dual_es) == 4
@test s[dual_es, :D_∂v0] == edge_center(s, 1:4)
@test elementary_duals(s, V(5)) == DualE(dual_es)

# 1D oriented dual complex
#-------------------------

primal_s = OrientedDeltaSet1D{Bool}()
add_vertices!(primal_s, 3)
add_edges!(primal_s, [1,2], [2,3], edge_orientation=[true,false])
s = OrientedDeltaDualComplex1D{Bool}(primal_s)
@test s[only(elementary_duals(0,s,1)), :D_edge_orientation] == true
@test s[only(elementary_duals(0,s,3)), :D_edge_orientation] == true

@test ∂(s, DualChain{1}([1,0,1])) isa DualChain{0}
@test d(s, DualForm{0}([1,1])) isa DualForm{1}
@test dual_boundary(1,s) == ∂(1,s)'
@test dual_derivative(0,s) == -d(0,s)'

# 1D embedded dual complex
#-------------------------

# Path graph on 3 vertices with irregular lengths.
primal_s = EmbeddedDeltaSet1D{Bool,Point2D}()
add_vertices!(primal_s, 3, point=[Point2D(1,0), Point2D(0,0), Point2D(0,2)])
add_edges!(primal_s, [1,2], [2,3], edge_orientation=true)
s = EmbeddedDeltaDualComplex1D{Bool,Float64,Point2D}(primal_s)
subdivide_duals!(s, Barycenter())
@test dual_point(s, edge_center(s, [1,2])) ≈ [Point2D(0.5,0), Point2D(0,1)]
@test volume(s, E(1:2)) ≈ [1.0, 2.0]
@test volume(s, elementary_duals(s, V(2))) ≈ [0.5, 1.0]
@test ⋆(0,s) ≈ Diagonal([0.5, 1.5, 1.0])
@test ⋆(1,s) ≈ Diagonal([1, 0.5])
@test ⋆(s, VForm([0,2,0]))::DualForm{1} ≈ DualForm{1}([0,3,0])

@test ∧(0,0,s, [1,2,3], [3,4,7]) ≈ [3,8,21]
@test ∧(s, VForm([1,2,3]), VForm([3,4,7]))::VForm ≈ VForm([3,8,21])
@test ∧(s, VForm([1,1,1]), EForm([2.5, 5.0]))::EForm ≈ EForm([2.5, 5.0])
@test ∧(s, VForm([1,1,0]), EForm([2.5, 5.0])) ≈ EForm([2.5, 2.5])
vform, eform = VForm([1.5, 2, 2.5]), EForm([13, 7])
@test ∧(s, vform, eform) ≈ ∧(s, eform, vform)

# Path graph on 5 vertices with regular lengths.
#
# Equals the graph Laplacian of the underlying graph, except at the boundary.
# Note that the DEC Laplace-Beltrami operator is *not* symmetric.
primal_s = EmbeddedDeltaSet1D{Bool,Point2D}()
add_vertices!(primal_s, 5, point=[Point2D(i,0) for i in -2:2])
add_edges!(primal_s, 1:4, 2:5, edge_orientation=true)
s = EmbeddedDeltaDualComplex1D{Bool,Float64,Point2D}(primal_s)
subdivide_duals!(s, Barycenter())
@test Δ(s, VForm([0,0,1,0,0])) ≈ VForm([0,1,-2,1,0])
@test Δ(0,s) ≈ [-2  2  0  0  0;
                 1 -2  1  0  0;
                 0  1 -2  1  0;
                 0  0  1 -2  1;
                 0  0  0  2 -2]

# 2D dual complex
#################

# Triangulated square.
primal_s = DeltaSet2D()
add_vertices!(primal_s, 4)
glue_triangle!(primal_s, 1, 2, 3)
glue_triangle!(primal_s, 1, 3, 4)
s = DeltaDualComplex2D(primal_s)
@test nparts(s, :DualV) == nv(primal_s) + ne(primal_s) + ntriangles(primal_s)
@test nparts(s, :DualE) == 2*ne(primal_s) + 6*ntriangles(primal_s)
@test nparts(s, :DualTri) == 6*ntriangles(primal_s)
@test primal_vertex(s, subsimplices(s, Tri(1)))::V == V([1,1,2,2,3,3])

dual_vs = elementary_duals(2,s,2)
@test dual_vs == [triangle_center(s,2)]
@test elementary_duals(s, Tri(2)) == DualV(dual_vs)
@test s[elementary_duals(1,s,2), :D_∂v1] == [edge_center(s,2)]
@test s[elementary_duals(1,s,3), :D_∂v1] == repeat([edge_center(s,3)], 2)
@test [length(elementary_duals(s, V(i))) for i in 1:4] == [4,2,4,2]
@test dual_triangle_vertices(s, 1) == [1,7,10]

# 2D oriented dual complex
#-------------------------

# Triangulated square with consistent orientation.
primal_s = OrientedDeltaSet2D{Bool}()
add_vertices!(primal_s, 4)
glue_triangle!(primal_s, 1, 2, 3, tri_orientation=true)
glue_triangle!(primal_s, 1, 3, 4, tri_orientation=true)
primal_s[:edge_orientation] = true
s = OrientedDeltaDualComplex2D{Bool}(primal_s)
@test sum(s[:D_tri_orientation]) == nparts(s, :DualTri) ÷ 2
@test [sum(s[elementary_duals(0,s,i), :D_tri_orientation])
       for i in 1:4] == [2,1,2,1]
@test sum(s[elementary_duals(1,s,3), :D_edge_orientation]) == 1

for k in 0:1
  @test dual_boundary(2-k,s) == (-1)^k * ∂(k+1,s)'
end
for k in 1:2
  # Desbrun, Kanso, Tong 2008, Equation 4.2.
  @test dual_derivative(2-k,s) == (-1)^k * d(k-1,s)'
end

# 2D embedded dual complex
#-------------------------

# Single triangle: numerical example from Gillette's notes on DEC, §2.13.
#
# Compared with Gillette, edges #2 and #3 are swapped in the ordering, which
# changes the discrete exterior derivative and other operators. The numerical
# values remain the same, as we verify.
primal_s = EmbeddedDeltaSet2D{Bool,Point2D}()
add_vertices!(primal_s, 3, point=[Point2D(0,0), Point2D(1,0), Point2D(0,1)])
glue_triangle!(primal_s, 1, 2, 3, tri_orientation=true)
primal_s[:edge_orientation] = true
s = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(primal_s)

subdivide_duals!(s, Barycenter())
@test dual_point(s, triangle_center(s, 1)) ≈ Point2D(1/3, 1/3)
@test volume(s, Tri(1)) ≈ 1/2
@test volume(s, elementary_duals(s, V(1))) ≈ [1/12, 1/12]
@test [sum(volume(s, elementary_duals(s, V(i)))) for i in 1:3] ≈ [1/6, 1/6, 1/6]
@test ⋆(0,s) ≈ Diagonal([1/6, 1/6, 1/6])
@test ⋆(1,s) ≈ Diagonal([√5/6, 1/6, √5/6])
@test ⋆(s, VForm([1,2,3]))::DualForm{2} ≈ DualForm{2}([1/6, 1/3, 1/2])
@test isapprox(δ(1,s), [ 2.236  0  2.236;
                        -2.236  1  0;
                         0     -1 -2.236], atol=1e-3)
@test δ(s, EForm([0.5,1.5,0.5])) isa VForm

@test ∧(s, VForm([2,2,2]), TriForm([2.5]))::TriForm ≈ TriForm([2.5])
vform, triform = VForm([1.5, 2, 2.5]), TriForm([7.5])
@test ∧(s, vform, triform) ≈ ∧(s, triform, vform)
eform1, eform2 = EForm([1.5, 2, 2.5]), EForm([3, 7, 10])
@test ∧(s, eform1, eform1)::TriForm ≈ TriForm([0])
@test ∧(s, eform1, eform2) ≈ -∧(s, eform2, eform1)

X♭, α = EForm([1.5, 2, 2.5]), DualForm{1}([3, 7, 10])
@test interior_product(s, X♭, α) isa DualForm{0}
@test length(interior_product_flat(1,s, X♭.data, α.data)) == 1
@test lie_derivative(s, X♭, α) isa DualForm{1}
@test length(lie_derivative_flat(1,s, X♭.data, α.data)) == 3

subdivide_duals!(s, Circumcenter())
@test dual_point(s, triangle_center(s, 1)) ≈ Point2D(1/2, 1/2)
@test ⋆(0,s) ≈ Diagonal([1/4, 1/8, 1/8])
@test ⋆(1,s) ≈ Diagonal([0.5, 0.0, 0.5])
@test δ(1,s) ≈ [ 2  0  2;
                -4  0  0;
                 0  0 -4]

subdivide_duals!(s, Incenter())
@test dual_point(s, triangle_center(s, 1)) ≈ Point2D(1/(2+√2), 1/(2+√2))
@test isapprox(⋆(0,s), Diagonal([0.146, 0.177, 0.177]), atol=1e-3)
@test isapprox(⋆(1,s), Diagonal([0.359, 0.207, 0.359]), atol=1e-3)
@test isapprox(δ(1,s), [ 2.449  0      2.449;
                        -2.029  1.172  0;
                         0     -1.172 -2.029], atol=1e-3)

end


text:
# Utilities for taking gradients through simulations

export
    extract_parameters,
    inject_gradients

"""
    extract_parameters(system, force_field)

Form a `Dict` of all parameters in a `System`, allowing gradients
to be tracked.
"""
function extract_parameters(sys, ff)
    params_dic = Dict()

    for at_data in sys.atoms_data
        key_prefix = "atom_$(at_data.atom_type)_"
        if !haskey(params_dic, key_prefix * "mass")
            at = ff.atom_types[at_data.atom_type]
            params_dic[key_prefix * "mass"] = at.mass
            params_dic[key_prefix * "σ"   ] = at.σ
            params_dic[key_prefix * "ϵ"   ] = at.ϵ
        end
    end

    for inter in values(sys.pairwise_inters)
        if inter isa LennardJones
            key_prefix = "inter_LJ_"
            params_dic[key_prefix * "weight_14"] = inter.weight_14
            params_dic[key_prefix * "weight_solute_solvent"] = inter.weight_solute_solvent
        elseif inter isa Coulomb
            key_prefix = "inter_CO_"
            params_dic[key_prefix * "weight_14"] = inter.weight_14
            params_dic[key_prefix * "coulomb_const"] = inter.coulomb_const
        elseif inter isa CoulombReactionField
            key_prefix = "inter_CRF_"
            params_dic[key_prefix * "dist_cutoff"] = inter.dist_cutoff
            params_dic[key_prefix * "solvent_dielectric"] = inter.solvent_dielectric
            params_dic[key_prefix * "weight_14"] = inter.weight_14
            params_dic[key_prefix * "coulomb_const"] = inter.coulomb_const
        end
    end

    for inter in values(sys.specific_inter_lists)
        if Molly.interaction_type(inter) <: HarmonicBond
            for bond_type in inter.types
                key_prefix = "inter_HB_$(bond_type)_"
                if !haskey(params_dic, key_prefix * "b0")
                    bond = ff.bond_types[atom_types_to_tuple(bond_type)]
                    params_dic[key_prefix * "b0"] = bond.b0
                    params_dic[key_prefix * "kb"] = bond.kb
                end
            end
        elseif Molly.interaction_type(inter) <: HarmonicAngle
            for angle_type in inter.types
                key_prefix = "inter_HA_$(angle_type)_"
                if !haskey(params_dic, key_prefix * "th0")
                    angle = ff.angle_types[atom_types_to_tuple(angle_type)]
                    params_dic[key_prefix * "th0"] = angle.th0
                    params_dic[key_prefix * "cth"] = angle.cth
                end
            end
        elseif Molly.interaction_type(inter) <: PeriodicTorsion
            for (torsion_type, torsion_inter) in zip(inter.types, Array(inter.inters))
                if torsion_inter.proper
                    key_prefix = "inter_PT_$(torsion_type)_"
                else
                    key_prefix = "inter_IT_$(torsion_type)_"
                end
                if !haskey(params_dic, key_prefix * "phase_1")
                    torsion = ff.torsion_types[atom_types_to_tuple(torsion_type)]
                    for i in 1:length(torsion.phases)
                        params_dic[key_prefix * "phase_$i"] = torsion.phases[i]
                        params_dic[key_prefix * "k_$i"    ] = torsion.ks[i]
                    end
                end
            end
        end
    end

    return params_dic
end

"""
    inject_gradients(sys, params_dic)

Add parameters from a dictionary to a `System`.
Allows gradients for individual parameters to be tracked.
Returns atoms, pairwise interactions and specific interaction lists.
"""
function inject_gradients(sys, params_dic, gpu::Bool=isa(sys.coords, CuArray))
    if gpu
        atoms_grad = cu(inject_atom.(Array(sys.atoms), sys.atoms_data, (params_dic,)))
    else
        atoms_grad = inject_atom.(sys.atoms, sys.atoms_data, (params_dic,))
    end
    pis_grad = inject_interaction.(sys.pairwise_inters, (params_dic,))
    sis_grad = inject_interaction_list.(sys.specific_inter_lists, (params_dic,), gpu)
    return atoms_grad, pis_grad, sis_grad
end

# get function errors with AD
dict_get(dic, key, default) = haskey(dic, key) ? dic[key] : default

function inject_atom(at, at_data, params_dic)
    key_prefix = "atom_$(at_data.atom_type)_"
    Atom(
        at.index,
        at.charge, # Residue-specific
        dict_get(params_dic, key_prefix * "mass"  , at.mass  ),
        dict_get(params_dic, key_prefix * "σ"     , at.σ     ),
        dict_get(params_dic, key_prefix * "ϵ"     , at.ϵ     ),
        at.solute,
    )
end

function inject_interaction_list(inter::InteractionList2Atoms, params_dic, gpu)
    if gpu
        inters_grad = cu(inject_interaction.(Array(inter.inters), inter.types, (params_dic,)))
    else
        inters_grad = inject_interaction.(inter.inters, inter.types, (params_dic,))
    end
    InteractionList2Atoms(inter.is, inter.js, inter.types, inters_grad)
end

function inject_interaction_list(inter::InteractionList3Atoms, params_dic, gpu)
    if gpu
        inters_grad = cu(inject_interaction.(Array(inter.inters), inter.types, (params_dic,)))
    else
        inters_grad = inject_interaction.(inter.inters, inter.types, (params_dic,))
    end
    InteractionList3Atoms(inter.is, inter.js, inter.ks, inter.types, inters_grad)
end

function inject_interaction_list(inter::InteractionList4Atoms, params_dic, gpu)
    if gpu
        inters_grad = cu(inject_interaction.(Array(inter.inters), inter.types, (params_dic,)))
    else
        inters_grad = inject_interaction.(inter.inters, inter.types, (params_dic,))
    end
    InteractionList4Atoms(inter.is, inter.js, inter.ks, inter.ls,
                          inter.types, inters_grad)
end

function inject_interaction(inter::LennardJones{S, C, W, WS, F, E}, params_dic) where {S, C, W, WS, F, E}
    key_prefix = "inter_LJ_"
    LennardJones{S, C, W, WS, F, E}(
        inter.cutoff,
        inter.nl_only,
        inter.lorentz_mixing,
        dict_get(params_dic, key_prefix * "weight_14", inter.weight_14),
        dict_get(params_dic, key_prefix * "weight_solute_solvent", inter.weight_solute_solvent),
        inter.force_units,
        inter.energy_units,
    )
end

function inject_interaction(inter::Coulomb, params_dic)
    key_prefix = "inter_CO_"
    Coulomb(
        inter.cutoff,
        inter.nl_only,
        dict_get(params_dic, key_prefix * "weight_14", inter.weight_14),
        dict_get(params_dic, key_prefix * "coulomb_const", inter.coulomb_const),
        inter.force_units,
        inter.energy_units,
    )
end

function inject_interaction(inter::CoulombReactionField, params_dic)
    key_prefix = "inter_CRF_"
    CoulombReactionField(
        dict_get(params_dic, key_prefix * "dist_cutoff", inter.dist_cutoff),
        dict_get(params_dic, key_prefix * "solvent_dielectric", inter.solvent_dielectric),
        inter.nl_only,
        dict_get(params_dic, key_prefix * "weight_14", inter.weight_14),
        dict_get(params_dic, key_prefix * "coulomb_const", inter.coulomb_const),
        inter.force_units,
        inter.energy_units,
    )
end

function inject_interaction(inter::HarmonicBond, inter_type, params_dic)
    key_prefix = "inter_HB_$(inter_type)_"
    HarmonicBond(
        dict_get(params_dic, key_prefix * "b0", inter.b0),
        dict_get(params_dic, key_prefix * "kb", inter.kb),
    )
end

function inject_interaction(inter::HarmonicAngle, inter_type, params_dic)
    key_prefix = "inter_HA_$(inter_type)_"
    HarmonicAngle(
        dict_get(params_dic, key_prefix * "th0", inter.th0),
        dict_get(params_dic, key_prefix * "cth", inter.cth),
    )
end

function inject_interaction(inter::PeriodicTorsion{N, T, E}, inter_type, params_dic) where {N, T, E}
    if inter.proper
        key_prefix = "inter_PT_$(inter_type)_"
    else
        key_prefix = "inter_IT_$(inter_type)_"
    end
    PeriodicTorsion{N, T, E}(
        inter.periodicities,
        ntuple(i -> dict_get(params_dic, key_prefix * "phase_$i", inter.phases[i]), N),
        ntuple(i -> dict_get(params_dic, key_prefix * "k_$i"    , inter.ks[i]    ), N),
        inter.proper,
    )
end


text:
using LatBo: Simulation
using LatBo.LB: D2Q9, D3Q19, Lattice
using LatBo.Indices: Indexing

facts("Lattice direction inversion") do
    for lattice_name in [:D2Q9, :D3Q19]
        context("for $(lattice_name)") do
            lattice = eval(lattice_name)
            ncomponents = size(lattice.celerities, 2)

            @fact lattice.inversion .>= 1 --> all
            @fact lattice.inversion .<= ncomponents --> all
            @fact length(unique(lattice.inversion)) --> ncomponents
            @fact lattice.celerities[:, lattice.inversion] --> -lattice.celerities
        end
    end
end

type FakeSimulation <: Simulation
    indexing::Indexing
    lattice::Lattice
end

facts("Compute neighbor index from sim/lattice object") do
    sim = FakeSimulation(Cartesian([10, 7]), D2Q9)
    for from = ([2, 3], [5, 6]), direction = 1:size(sim.lattice.celerities, 2)
        site = index(sim.indexing, from)
        to = neighbor_index(sim, site, direction)
        @fact gridcoords(sim.indexing, to) - from --> sim.lattice.celerities[:, direction]
    end
    @fact_throws AssertionError neighbor_index(sim, 0, 1)
    @fact_throws AssertionError neighbor_index(sim, length(sim.indexing) + 1, 1)
end


text:
using ParserCombinator

mutable struct Direction
    current::Char
    next::Array{Direction}
end

function unwrap_any(a)
    if typeof(a) == Array{Any,1}
        return unwrap_any(a[1])
    end
    return a
end

function Direction(d::Array)::Direction
    return Direction(d[1][1], [])
end

function Joiner(a)
    if length(a) == 1
        return a
    elseif length(a) == 2
        a2 = collect(Iterators.flatten(a[2]))
        for n in a2
            push!(a[1].next, n[1][1])
        end
        return a2
    else
        println("LEngth more than 2")
        println(typeof(a))
        println(unwrap_any(a[2]))
    end
end

Dir = p"[NSEW]"|> Direction
All = Delayed()
Normal = Delayed()
Branch = Delayed()

Normal.matcher = Dir + Star(All) |> Joiner
Branch.matcher = E"(" + All + Star(E"|" + All) + E")" |> identity
All.matcher = Normal | Branch |> identity
Parser = All + Eos()

parse_one("N(S)W", Parser)


text:
# Copyright 2019, Patricia Levi and contributors

###################
# install & load packages
###################
## NB: this setup may take some time to run - on the order of a handful of minutes.
## You only need to run "Pkg.add" functions once, though
## The first time you call 'using' on a package, it will precomile, which will also take a few minutes
using Pkg
Pkg.add(["CapacityExpansion","Clp", "StatsPlots", "DataFrames", "CSV"])

# load packages - you must do this every time you restart julia
using CapacityExpansion, Clp, StatsPlots, DataFrames, CSV

### Add utils functions ###
## IF YOU ARE ABLE TO UPLOAD 'homework_utils.jl' TO FARMSHARE/
## THE DIRECTORY WHERE YOU ARE RUNNING THIS SCRIPT:
include(joinpath(@__DIR__,"homework_utils.jl"))
## OTHERWISE: COPY-PASTE THE CONTENTS OF plot_scripts.jl
## INTO THE COMMAND LINE

###################
### Set up data ###
###################
state="GER_18" # data is also provided for "GER_1", "CA_1", and "CA_14" - https://youngfaithful.github.io/CapacityExpansion.jl/dev/csv_structure/#Provided-Data-1
years=[2015]

# load non-timeseries data
cep_data = load_cep_data_provided(state)

# load timeseries data
ts_input_data = load_timeseries_data_provided(state; #input the states-name
        T=24, #input the number of time-steps per original period -> 24 hours per period
        years=years #input the years of data used
        )

# cluster timeseries data - run aggregation with a hierarchical clustering algorithm
# An aggregation of the temporal data is useful to reduce the computational
# complexity of a capacity expansion problem - further information: https://holgerteichgraeber.github.io/TimeSeriesClustering.jl/stable/
ts_clust_data = run_clust(ts_input_data;
        method="hierarchical", # use a hierarchical clustering method to groupt the periods
        representation="centroid", # represent the cluster with a centroid representation (the mean of the group)
        n_init=1, # choose the best result from a single run (partitional clustering methods instead of the hierarchical clusterin method need multiple initializations)
        n_clust=10 # choose the number of representative periods that represent the full time-series (e.g. we use 10 representative days to represent the original 365 days of a year to reduce the temporal complexity of the model)
        )
##############################
## PART 2: Initial CO2 Runs ##
##############################
# set up an optimizer
# An optimizer is used to solve the linear optimization problem that is constructed
# by CapacityExpansion.jl. https://youngfaithful.github.io/CapacityExpansion.jl/dev/opt_cep/#Solver-1
optimizer=Clp.Optimizer

# tweak the CO2 level
co2_low = run_opt(ts_clust_data.clust_data,cep_data,optimizer;
    descriptor = "co2_low", limit_emission=Dict{String,Number}("CO2/electricity"=>#=your choice=#))

co2_mid = run_opt(ts_clust_data.clust_data,cep_data,optimizer;
    descriptor = "co2_mid",limit_emission=Dict{String,Number}("CO2/electricity"=>#=your choice=#))

co2_high = run_opt(ts_clust_data.clust_data,cep_data,optimizer;
    descriptor = "co2_high",limit_emission=Dict{String,Number}("CO2/electricity"=>#=your choice=#))

results = [co2_low,co2_mid,co2_high]

# plot & inspect the results
plotcapacity(results,"baseruncapacity")
plotgen(results,"baserungeneration")
plotcost(results,"baseruncost")
# these functions also save a CSV of the displayed data, with the same name
# as the plot image
###################################
## PART 3: testing assumptions   ##
###################################
co2_limit_mid=co2_mid.config["limit_emission"]["CO2"]["electricity"]
##----------------
## Testing greenfield/brownfield
brownfield = run_opt(ts_clust_data.clust_data,cep_data,optimizer;
    descriptor = "brownf",
    limit_emission=Dict{String,Number}("CO2/electricity"=>co2_limit_mid),
    infrastructure=Dict{String,Array}("existing"=>["all"]))

# plot & inspect the results
results = [co2_mid, brownfield]
plotcapacity(results,"brownfieldcapacity")
plotgen(results,"brownfieldgeneration")
plotcost(results,"brownfieldcost")

##----------------
## Testing transmission constraints
trans_constrain = run_opt(ts_clust_data.clust_data,cep_data,optimizer;
    descriptor = "trans_constrain",
    limit_emission=Dict{String,Number}("CO2/electricity"=>co2_limit_mid),
    transmission = true
    )

# plot & inspect the results
results = [co2_mid, trans_constrain]
plotcapacity(results,"trans_constraincapacity")
plotgen(results,"trans_constraingeneration")
plotcost(results,"trans_constraincost")

###################################
## PART 4: storage sensitivities ##
###################################

##----------------
## how do the results change with storage?

simplestor = run_opt(ts_clust_data.clust_data,cep_data,optimizer;
    descriptor = "seasonal_stor",
    limit_emission=Dict{String,Number}("CO2/electricity"=>co2_limit_mid),
    storage_type="simple",
    conversion=true
    )
# In order to store hydrogen it needs to be converted from electriciy to hydrogen: need to activate conversion and select a storage_type
results = [co2_mid,seasonalstor]

# plot & inspect the results
plotcapacity(results,"storage_capacity")
plotgen(results,"storage_generation")
plotcost(results,"storage_cost")

##----------------
## change the 'simple' battery costs ##

# see current battery cost:
first(cep_data.costs["bat_e",:,:,"cap_fix","EUR"])

# make new cep_data with new battery costs
batt_cost_change = #your choice - fraction of previous cost#
cep_data2 = deepcopy(cep_data)
cep_data2.costs["bat_e",:,:,"cap_fix","EUR"] = cep_data.costs["bat_e",:,:,"cap_fix","EUR"].*batt_cost_change;
#= copy paste the above 3 lines, changing the name of "cep_data"
to create different datasets with different battery cost assumptions =#

# run model with your new cep_data versions
batt1 = run_opt(ts_clust_data.clust_data,#=new cep data=#,optimizer;
    descriptor = #=your choice =#,
    storage = "simple",
    conversion = true,
    limit_emission=Dict{String,Number}("CO2/electricity"=>co2_limit_mid))
batt2 = #...

# gather all results and plot
results =[simplestor,batt1, batt2 #=other batt versions =#]
plotcapacity(results,"battery_capacity")
plotgen(results,"battery_generation")
plotcost(results,"battery_cost")

##----------------
## change the seasonal storage costs
## by changing hydrogen costs

# see current hydrogen cost:
first(cep_data.costs["h2_in",:,:,"cap_fix","EUR"])

# make new cep_data with new hydrogen costs
h2_cost_change = #your choice#
cep_data2 = deepcopy(cep_data)
cep_data2.costs["h2_in",:,:,"cap_fix","EUR"] = cep_data.costs["h2_in",:,:,"cap_fix","EUR"].*h2_cost_change;
#= copy paste the above 3 lines, changing the name of "cep_data"
to create different datasets with different battery cost assumptions =#

run model with your new cep_data versions
hy1 = run_opt(ts_clust_data.clust_data,#=new cep data =#,optimizer;
    descriptor =#=your choice =#, co2_limit=#=your choice =#)
hy2 = #...

# gather all results and plot
results =[hy1, hy2 #=other hydrogen sensitivities =#]

plotcapacity(results,"hydrogen_capacity")
plotgen(results,"hydrogen_generation")
plotcost(results,"hydrogen_cost")

###################################
## PART 5: final analysis        ##
###################################
#=
if you want to run additional scenarios to inform your
understanding, do so here, following the workflow of
(1) modify costs
(2) run model
(3) combine model results that youd like to compare
(4) call plotting functions
You should be able to repurpose code from above to do so.
=#

#General alternative formulation for multiple calculations:
results = Array{OptResult,1}() #Setting up a new results Array that will be filled with results
for co2_limit in collect(100:100:500) #Loop through different configurations (e.g. a co2_limit)
    push!(results, run_opt(ts_clust_data.clust_data,cep_data,optimizer;
        descriptor = string(co2_limit),
        limit_emission=Dict{String,Number}("CO2/electricity"=>co2_limit))) # Push the new result from run_opt(your configuration) to the results-Array
end


text:
module Kronecker

export Kronecker

"""
    Kronecker
"""
struct Kronecker{MT <: AbstractMatrix} <: AbstractMatrix
    # data
    A::Vector{MT}
end

# constructor
Kronecker(args...) = Kronecker([args...])

end


text:
#### Hyperdiffusion Model Functions
using DocStringExtensions
using LinearAlgebra
export HyperDiffusion, NoHyperDiffusion, StandardHyperDiffusion

abstract type HyperDiffusion end
vars_state_conservative(::HyperDiffusion, FT) = @vars()
vars_state_auxiliary(::HyperDiffusion, FT) = @vars()
vars_state_gradient(::HyperDiffusion, FT) = @vars()
vars_gradient_laplacian(::HyperDiffusion, FT) = @vars()
vars_state_gradient_flux(::HyperDiffusion, FT) = @vars()
vars_hyperdiffusive(::HyperDiffusion, FT) = @vars()
function atmos_init_aux!(
    ::HyperDiffusion,
    ::AtmosModel,
    aux::Vars,
    geom::LocalGeometry,
) end
function atmos_nodal_update_auxiliary_state!(
    ::HyperDiffusion,
    ::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function compute_gradient_argument!(
    ::HyperDiffusion,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function transform_post_gradient_laplacian!(
    h::HyperDiffusion,
    hyperdiffusive::Vars,
    gradvars::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function flux_second_order!(
    h::HyperDiffusion,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    aux::Vars,
    t::Real,
) end
function compute_gradient_flux!(
    h::HyperDiffusion,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
) end

"""
  NoHyperDiffusion <: HyperDiffusion
Defines a default hyperdiffusion model with zero diffusive fluxes.
"""
struct NoHyperDiffusion <: HyperDiffusion end

"""
  StandardHyperDiffusion{FT} <: HyperDiffusion
Horizontal hyperdiffusion methods for application in GCM and LES settings
Timescales are prescribed by the user while the diffusion coefficient is
computed as a function of the grid lengthscale.
"""
struct StandardHyperDiffusion{FT} <: HyperDiffusion
    τ_timescale::FT
end

vars_state_auxiliary(::StandardHyperDiffusion, FT) = @vars(Δ::FT)
vars_state_gradient(::StandardHyperDiffusion, FT) =
    @vars(u::SVector{3, FT}, h_tot::FT)
vars_gradient_laplacian(::StandardHyperDiffusion, FT) =
    @vars(u::SVector{3, FT}, h_tot::FT)
vars_hyperdiffusive(::StandardHyperDiffusion, FT) =
    @vars(ν∇³u::SMatrix{3, 3, FT, 9}, ν∇³h_tot::SVector{3, FT})

function atmos_init_aux!(
    ::StandardHyperDiffusion,
    ::AtmosModel,
    aux::Vars,
    geom::LocalGeometry,
)
    aux.hyperdiffusion.Δ = lengthscale(geom)
end

function compute_gradient_argument!(
    h::StandardHyperDiffusion,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρinv = 1 / state.ρ
    transform.hyperdiffusion.u = ρinv * state.ρu
    transform.hyperdiffusion.h_tot = transform.h_tot
end

function transform_post_gradient_laplacian!(
    h::StandardHyperDiffusion,
    hyperdiffusive::Vars,
    hypertransform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ∇Δu = hypertransform.hyperdiffusion.u
    ∇Δh_tot = hypertransform.hyperdiffusion.h_tot
    # Unpack
    τ_timescale = h.τ_timescale
    # Compute hyperviscosity coefficient
    ν₄ = (aux.hyperdiffusion.Δ / 2)^4 / 2 / τ_timescale
    hyperdiffusive.hyperdiffusion.ν∇³u = ν₄ * ∇Δu
    hyperdiffusive.hyperdiffusion.ν∇³h_tot = ν₄ * ∇Δh_tot
end

function flux_second_order!(
    h::StandardHyperDiffusion,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    aux::Vars,
    t::Real,
)
    flux.ρu += state.ρ * hyperdiffusive.hyperdiffusion.ν∇³u
    flux.ρe += state.ρ * hyperdiffusive.hyperdiffusion.ν∇³h_tot
end


text:
"""
Type representing a third-order tensor
"""
struct ThirdOrderTensorValue{D1,D2,D3,T,L} <: MultiValue{Tuple{D1,D2,D3},T,3,L}
    data::NTuple{L,T}
    function ThirdOrderTensorValue{D1,D2,D3,T}(data::NTuple{L,T}) where {D1,D2,D3,T,L}
        @assert L == D1*D2*D3
        new{D1,D2,D3,T,L}(data)
    end
end

# Empty ThirdOrderTensorValue constructor

ThirdOrderTensorValue()                       = ThirdOrderTensorValue{0,0,0,Int}(NTuple{0,Int}())
ThirdOrderTensorValue{0,0,0}()                = ThirdOrderTensorValue{0,0,0,Int}(NTuple{0,Int}())
ThirdOrderTensorValue{0,0,0,T}() where{T}     = ThirdOrderTensorValue{0,0,0,T}(NTuple{0,T}())
ThirdOrderTensorValue(data::NTuple{0})        = ThirdOrderTensorValue{0,0,0,Int}(data)
ThirdOrderTensorValue{0,0,0}(data::NTuple{0}) = ThirdOrderTensorValue{0,0,0,Int}(data)

# ThirdOrderTensorValue single NTuple argument constructor

@generated function ThirdOrderTensorValue(data::NTuple{L,T}) where {L,T}
  D=Int(cbrt(L))
  quote
    ThirdOrderTensorValue{$D,$D,$D,T}(data)
  end
end
ThirdOrderTensorValue{D}(data::NTuple{L,T}) where {D,L,T}                         = ThirdOrderTensorValue{D,D,D,T}(data)
ThirdOrderTensorValue{D1,D2,D3}(data::NTuple{L,T}) where {D1,D2,D3,L,T}           = ThirdOrderTensorValue{D1,D2,D3,T}(data)
ThirdOrderTensorValue{D1,D2,D3,T1}(data::NTuple{L,T2}) where {D1,D2,D3,L,T1,T2}   = ThirdOrderTensorValue{D1,D2,D3,T1}(NTuple{L,T1}(data))
ThirdOrderTensorValue{D1,D2,D3,T1,L}(data::NTuple{L,T2}) where {D1,D2,D3,L,T1,T2} = ThirdOrderTensorValue{D1,D2,D3,T1}(NTuple{L,T1}(data))

# ThirdOrderTensorValue single Tuple argument constructor

ThirdOrderTensorValue(data::Tuple) = ThirdOrderTensorValue(promote(data...))
ThirdOrderTensorValue{D}(data::Tuple) where {D} = ThirdOrderTensorValue{D,D,D}(promote(data...))
ThirdOrderTensorValue{D1,D2,D3}(data::Tuple) where {D1,D2,D3} = ThirdOrderTensorValue{D1,D2,D3}(promote(data...))
ThirdOrderTensorValue{D1,D2,D3,T1}(data::Tuple) where {D1,D2,D3,T1}  = ThirdOrderTensorValue{D1,D2,D3,T1}(NTuple{length(data),T1}(data))
ThirdOrderTensorValue{D1,D2,D3,T1,L}(data::Tuple) where {D1,D2,D3,T1,L}  = ThirdOrderTensorValue{D1,D2,D3,T1}(NTuple{L,T1}(data))

# ThirdOrderTensorValue Vararg constructor

ThirdOrderTensorValue(data...) = ThirdOrderTensorValue(data)
ThirdOrderTensorValue{D}(data...) where {D} = ThirdOrderTensorValue{D}(data)
ThirdOrderTensorValue{D1,D2,D3}(data...) where {D1,D2,D3} = ThirdOrderTensorValue{D1,D2,D3}(data)
ThirdOrderTensorValue{D1,D2,D3,T1}(data...) where {D1,D2,D3,T1} =  ThirdOrderTensorValue{D1,D2,D3,T1}(data)
ThirdOrderTensorValue{D1,D2,D3,T1,L}(data...) where {D1,D2,D3,T1,L} =  ThirdOrderTensorValue{D1,D2,D3,T1}(data)

# ThirdOrderTensorValue single AbstractArray{3,T} argument constructor

ThirdOrderTensorValue(data::AbstractArray{T,3}) where {T} = ((D1,D2,D3)=size(data);L=length(data);ThirdOrderTensorValue{D1,D2,D3,T}(NTuple{L,T}(data)))
ThirdOrderTensorValue{D}(data::AbstractArray{T,3}) where {D,T} = (L=length(data);ThirdOrderTensorValue{D,D,D,T}(NTuple{L,T}(data)))
ThirdOrderTensorValue{D1,D2,D3}(data::AbstractArray{T,3}) where {D1,D2,D3,T} = (L=length(data);ThirdOrderTensorValue{D1,D2,D3,T}(NTuple{L,T}(data)))
ThirdOrderTensorValue{D1,D2,D3,T1}(data::AbstractArray{T2,3}) where {D1,D2,D3,T1,T2} = (L=length(data);ThirdOrderTensorValue{D1,D2,D3,T1}(NTuple{L,T1}(data)))
ThirdOrderTensorValue{D1,D2,D3,T1,L}(data::AbstractArray{T2,3}) where {D1,D2,D3,T1,T2,L} = ThirdOrderTensorValue{D1,D2,D3,T1}(NTuple{L,T1}(data))

###############################################################
# Conversions (ThirdOrderTensorValue)
###############################################################

# Direct conversion
convert(::Type{<:ThirdOrderTensorValue{D1,D2,D3,T}}, arg::AbstractArray) where {D1,D2,D3,T} = ThirdOrderTensorValue{D1,D2,D3,T}(arg)
convert(::Type{<:ThirdOrderTensorValue{D1,D2,D3,T}}, arg::Tuple) where {D1,D2,D3,T} = ThirdOrderTensorValue{D1,D2,D3,T}(arg)

# Inverse conversion
convert(::Type{<:SMatrix{D1,D2,D3,T}}, arg::ThirdOrderTensorValue) where {D1,D2,D3,T} = SMatrix{D1,D2,D3,T}(Tuple(arg))
convert(::Type{<:MMatrix{D1,D2,D3,T}}, arg::ThirdOrderTensorValue) where {D1,D2,D3,T} = MMatrix{D1,D2,D3,T}(Tuple(arg))
convert(::Type{<:NTuple{L,T1}}, arg::ThirdOrderTensorValue) where {L,T1} = NTuple{L,T1}(Tuple(arg))

# Internal conversion
convert(::Type{<:ThirdOrderTensorValue{D1,D2,D3,T}}, arg::ThirdOrderTensorValue{D1,D2,D3}) where {D1,D2,D3,T} = ThirdOrderTensorValue{D1,D2,D3,T}(Tuple(arg))
convert(::Type{<:ThirdOrderTensorValue{D1,D2,D3,T}}, arg::ThirdOrderTensorValue{D1,D2,D3,T}) where {D1,D2,D3,T} = arg

# other

change_eltype(::Type{ThirdOrderTensorValue{D1,D2,D3,T1,L}},::Type{T2}) where {D1,D2,D3,T1,T2,L} = ThirdOrderTensorValue{D1,D2,D3,T2,L}
change_eltype(::T,::Type{T2}) where {T<:ThirdOrderTensorValue,T2} = change_eltype(T,T2)

zero(::Type{<:ThirdOrderTensorValue{D1,D2,D3,T}}) where {D1,D2,D3,T} = ThirdOrderTensorValue{D1,D2,D3,T}(tfill(zero(T),Val{D1*D2*D3}()))
zero(::ThirdOrderTensorValue{D1,D2,D3,T}) where {D1,D2,D3,T} = zero(ThirdOrderTensorValue{D1,D2,D3,T})

Mutable(::Type{<:ThirdOrderTensorValue{D1,D2,D3,T}}) where {D1,D2,D3,T} = MArray{Tuple{D1,D2,D3},T}
Mutable(::ThirdOrderTensorValue{D1,D2,D3,T}) where {D1,D2,D3,T} = Mutable(ThirdOrderTensorValue{D1,D2,D3,T})
mutable(a::ThirdOrderTensorValue{D1,D2,D3}) where {D1,D2,D3} = MArray{Tuple{D1,D2,D3}}(a.data)

###############################################################
# Introspection (ThirdOrderTensorValue)
###############################################################

eltype(::Type{<:ThirdOrderTensorValue{D1,D2,D3,T}}) where {D1,D2,D3,T} = T
eltype(::ThirdOrderTensorValue{D1,D2,D3,T}) where {D1,D2,D3,T} = eltype(ThirdOrderTensorValue{D1,D2,D3,T})

size(::Type{<:ThirdOrderTensorValue{D1,D2,D3}}) where {D1,D2,D3} = (D1,D2,D3)
size(::ThirdOrderTensorValue{D1,D2,D3}) where {D1,D2,D3} = size(ThirdOrderTensorValue{D1,D2,D3})

length(::Type{<:ThirdOrderTensorValue{D1,D2,D3}}) where {D1,D2,D3} = D1*D2*D3
length(::ThirdOrderTensorValue{D1,D2,D3}) where {D1,D2,D3} = length(ThirdOrderTensorValue{D1,D2,D3})

num_components(::Type{<:ThirdOrderTensorValue{D1,D2,D3}}) where {D1,D2,D3} = length(ThirdOrderTensorValue{D1,D2,D3})
num_components(::ThirdOrderTensorValue{D1,D2,D3}) where {D1,D2,D3} = num_components(ThirdOrderTensorValue{D1,D2,D3})


text:
## Parameter types used by Junet for multiple dispatch ##

"""
    DirParam

Abstract parent for parameter types encoding directedness.

In Junet, all data is stored in data structures that are agnostic of graphs'
directionality. The directionality is encoded in type parameters, for which
`DirParam` provides an abstract parent type.

The interpretation of those parameter types depends on each particular type that
uses them, but there are some commonalities:
  * `Forward` means interpreting data right as it's stored, e.g., look in the
    `forward` field to find info about edges originating from a node.
  * `Reverse` is the opposite of `Forward`.
  * `Both` means to look in both directions. Typically leads to processing two
    more data than with the previous ones.

For good presentation, they have two parents — `Directed` and `Undirected`,
by the type of graphs in which they are used.
"""
abstract type DirParam end
abstract type Directed <: DirParam end
abstract type Forward  <: Directed end
abstract type Reverse  <: Directed end
abstract type Undirected <: DirParam end
abstract type Both       <: Undirected end

dir_xor(::Type{D}, ::Type{R}) where {D<:Forward, R<:Forward} = Forward
dir_xor(::Type{D}, ::Type{R}) where {D<:Forward, R<:Reverse} = Reverse
dir_xor(::Type{D}, ::Type{R}) where {D<:Reverse, R<:Forward} = Reverse
dir_xor(::Type{D}, ::Type{R}) where {D<:Reverse, R<:Reverse} = Forward
dir_xor(::Type{D}, ::Type{R}) where {D<:DirParam, R<:DirParam} = Both

"""
    MultiParam

Similar to `DirParam`, this abstract type is for encoding whether graph may
have multiple edges. There are two child types: `Multi` amd `Simple`.
"""
abstract type MultiParam end
abstract type Multi  <: MultiParam end
abstract type Simple <: MultiParam end