text:
using JuAFEM, SparseArrays
using JuAFEM:Vec
using StaticArrays

using FEMRefinements

v1 = Vec((2.0, 0.0))
v2 = Vec((0.0, 2.0))

grid = generate_grid2D(1, 2, v1, v2);
# FEMRefinements.refine!(grid, [2])
# FEMRefinements.refine!(grid, [1])
FEMRefinements.refine!(grid, collect(1:getncells(grid)))
FEMRefinements.refine!(grid, [1,4,8])
# FEMRefinements.refine!(grid, collect(1:getncells(grid)))
FEMRefinements.refine!(grid, [1,4])
# FEMRefinements.refine!(grid, collect(1:getncells(grid)))
# FEMRefinements.refine!(grid, [getncells(grid) - 1,getncells(grid)])
# FEMRefinements.refine!(grid, [getncells(grid)])
# FEMRefinements.refine!(grid, [getncells(grid) - 1,getncells(grid)])
# FEMRefinements.refine!(grid, [15])

dim = 2
ip = Lagrange{dim,RefCube,1}()
qr = QuadratureRule{dim,RefCube}(2)
cellvalues = CellScalarValues(qr, ip);
dh = DofHandler(grid)
push!(dh, :u, 1)
dh, vertexdict, _, _ = JuAFEM.__close!(dh);
K = create_sparsity_pattern(dh);

# using UnicodePlots
# fill!(K.nzval, 1.0)
# spy(K; height = 15)

ch = ConstraintHandler(dh);

∂Ω = union(getfaceset.((grid,), ["left"])...);

dbc = Dirichlet(:u, ∂Ω, (x, t) -> 0)

function doassemble(cellvalues::CellScalarValues{dim}, K::SparseMatrixCSC, dh::DofHandler) where {dim}
    n_basefuncs = getnbasefunctions(cellvalues)
    Ke = zeros(n_basefuncs, n_basefuncs)
    fe = zeros(n_basefuncs)
    f = zeros(ndofs(dh))
    assembler = start_assemble(K, f)

    @inbounds for cell in CellIterator(dh)

        fill!(Ke, 0)
        fill!(fe, 0)

        reinit!(cellvalues, cell)

        for q_point in 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)

            for i in 1:n_basefuncs
                v  = shape_value(cellvalues, q_point, i)
                ∇v = shape_gradient(cellvalues, q_point, i)
                fe[i] += v * dΩ
                for j in 1:n_basefuncs
                    ∇u = shape_gradient(cellvalues, q_point, j)
                    Ke[i, j] += (∇v ⋅ ∇u) * dΩ
                end
            end
        end

        assemble!(assembler, celldofs(cell), fe, Ke)
    end
    return K, f
end

K, f = doassemble(cellvalues, K, dh);
ncinfo = getncinfo(grid)
hangingnodes = getindex.(ncinfo, 1)
mastercells = getindex.(ncinfo, 2)
slavecells = getindex.(ncinfo, 3)
ncdofs = get.((vertexdict[1],), hangingnodes, NaN) ## ist das mapping nicht andersrum?
_ndofs = ndofs(dh)

I, J, V = Int64[], Int64[], Float64[]

ip = Lagrange{2,RefCube,1}()
tdof = 0 # TODO ++ if conforming dof
for dof in 1:_ndofs
    ncidx = findfirst(x -> x == dof, ncdofs)
    if ncidx !== nothing
        occurence = findall(x -> x == mastercells[ncidx], mastercells)
        # celldepth =
        if length(occurence) == 1
            mastercellidx = getindex.(mastercells, 1)[ncidx]
            mastercoords = getcoordinates(grid, mastercellidx)
            masterdofs = celldofs(dh, mastercellidx)
            hangingnode = grid.nodes[dof]
            facecoords = JuAFEM.faces(grid.cells[mastercellidx])
            masterfacedofs, masterdofs = facedofs(dh, mastercellidx, getindex.(mastercells, 2)[ncidx])
            local_coordinate = find_local_coordinate(Lagrange{2,RefCube,1}(), mastercoords, hangingnode.x)
            quad_rule = QuadratureRule{2,RefCube,Float64}([1], [local_coordinate])
            cv = CellScalarValues(quad_rule, ip)
            reinit!(cv, mastercoords)
            for i in 1:length(masterfacedofs)
                Pvalue = JuAFEM.shape_value(cv, 1, findfirst(x -> x == masterfacedofs[i], masterdofs))
                push!(I, dof)
                push!(J, masterfacedofs[i])
                push!(V, Pvalue)
            end
        end
    else
        push!(I, dof)
        push!(J, dof)
        push!(V, 1.0)
    end
end

P = sparse(I, J, V, _ndofs, _ndofs)

JuAFEM.add!(ch, dbc,ncdofs);
JuAFEM.close!(ch,ncdofs)
JuAFEM.update!(ch, ncdofs, 0.0);
P = P[:, 1:end .∉ [ncdofs]]
A = (P' * K) * P
f = P' * f
apply!(A, f, ch)
u = A \ f;

# vtk_grid("heat_equation", dh) do vtk
#    vtk_point_data(vtk, dh, P*u)
# end

# using Test                        #src
# @test norm(u) ≈ 3.307743912641305 #src

FEMRefinements.plot_solution3D(dh,P * u; plotgrid=true)


text:
@testset "Processes" begin

  # need expanding
  process_result = Process("Proc 1")
  @test process_result.Name == "Proc 1"
  @test name(process_result) == "Proc 1"

  proc_2 = JSON.parse("""
  {
    "Name": "Proc 2"
  }
  """)

  @test name(Process(proc_2)) == "Proc 2"
  @test Process(proc_2) != process_result

end

ext: jl


text:
#### Domain && Support

immutable RealInterval
    lb::Float64
    ub::Float64

    RealInterval(lb::Real, ub::Real) = new(Float64(lb), Float64(ub))
end

minimum(r::RealInterval) = r.lb
maximum(r::RealInterval) = r.ub
in(x::Real, r::RealInterval) = (r.lb <= Float64(x) <= r.ub)

isbounded{D<:UnivariateDistribution}(d::Union{D,Type{D}}) = isupperbounded(d) && islowerbounded(d)

islowerbounded{D<:UnivariateDistribution}(d::Union{D,Type{D}}) = minimum(d) > -Inf
isupperbounded{D<:UnivariateDistribution}(d::Union{D,Type{D}}) = maximum(d) < +Inf

hasfinitesupport{D<:DiscreteUnivariateDistribution}(d::Union{D,Type{D}}) = isbounded(d)
hasfinitesupport{D<:ContinuousUnivariateDistribution}(d::Union{D,Type{D}}) = false

function insupport!{D<:UnivariateDistribution}(r::AbstractArray, d::Union{D,Type{D}}, X::AbstractArray)
    length(r) == length(X) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    for i in 1 : length(X)
        @inbounds r[i] = insupport(d, X[i])
    end
    return r
end

insupport{D<:UnivariateDistribution}(d::Union{D,Type{D}}, X::AbstractArray) =
     insupport!(BitArray(size(X)), d, X)

insupport{D<:ContinuousUnivariateDistribution}(d::Union{D,Type{D}},x::Real) = minimum(d) <= x <= maximum(d)
insupport{D<:DiscreteUnivariateDistribution}(d::Union{D,Type{D}},x::Real) = isinteger(x) && minimum(d) <= x <= maximum(d)

support{D<:ContinuousUnivariateDistribution}(d::Union{D,Type{D}}) = RealInterval(minimum(d), maximum(d))
support{D<:DiscreteUnivariateDistribution}(d::Union{D,Type{D}}) = round(Int, minimum(d)):round(Int, maximum(d))

# Type used for dispatch on finite support
# T = true or false
immutable FiniteSupport{T} end

## macros to declare support

macro distr_support(D, lb, ub)
    D_has_constantbounds = (isa(ub, Number) || ub == :Inf) &&
                           (isa(lb, Number) || lb == :(-Inf))

    paramdecl = D_has_constantbounds ? :(d::Union{$D, Type{$D}}) : :(d::$D)

    # overall
    esc(quote
        minimum($(paramdecl)) = $lb
        maximum($(paramdecl)) = $ub
    end)
end

##### generic methods (fallback) #####

## sampling

rand(d::UnivariateDistribution) = quantile(d, rand())

rand!(d::UnivariateDistribution, A::AbstractArray) = _rand!(sampler(d), A)
rand(d::UnivariateDistribution, n::Int) = _rand!(sampler(d), Array(eltype(d), n))
rand(d::UnivariateDistribution, shp::Dims) = _rand!(sampler(d), Array(eltype(d), shp))

## statistics

std(d::UnivariateDistribution) = sqrt(var(d))
median(d::UnivariateDistribution) = quantile(d, 0.5)
modes(d::UnivariateDistribution) = [mode(d)]
entropy(d::UnivariateDistribution, b::Real) = entropy(d) / log(b)

isplatykurtic(d::UnivariateDistribution) = kurtosis(d) > 0.0
isleptokurtic(d::UnivariateDistribution) = kurtosis(d) < 0.0
ismesokurtic(d::UnivariateDistribution) = kurtosis(d) == 0.0

function kurtosis(d::Distribution, correction::Bool)
    if correction
        return kurtosis(d)
    else
        return kurtosis(d) + 3.0
    end
end

excess(d::Distribution) = kurtosis(d)
excess_kurtosis(d::Distribution) = kurtosis(d)
proper_kurtosis(d::Distribution) = kurtosis(d, false)

#### pdf, cdf, and friends

# pdf

pdf(d::DiscreteUnivariateDistribution, x::Int) = throw(MethodError(pdf, (d, x)))
pdf(d::DiscreteUnivariateDistribution, x::Integer) = pdf(d, round(Int, x))
pdf(d::DiscreteUnivariateDistribution, x::Real) = isinteger(x) ? pdf(d, round(Int, x)) : 0.0

pdf(d::ContinuousUnivariateDistribution, x::Real) = throw(MethodError(pdf, (d, x)))

# logpdf

logpdf(d::DiscreteUnivariateDistribution, x::Int) = log(pdf(d, x))
logpdf(d::DiscreteUnivariateDistribution, x::Integer) = logpdf(d, round( Int, x))
logpdf(d::DiscreteUnivariateDistribution, x::Real) = isinteger(x) ? logpdf(d, round(Int, x)) : -Inf

logpdf(d::ContinuousUnivariateDistribution, x::Real) = log(pdf(d, x))

# cdf
cdf(d::DiscreteUnivariateDistribution, x::Int) = cdf(d, x, FiniteSupport{hasfinitesupport(d)})

# Discrete univariate with infinite support
function cdf(d::DiscreteUnivariateDistribution, x::Int, ::Type{FiniteSupport{false}})
    c = 0.0
    for y = minimum(d):x
        c += pdf(d, y)
    end
    return c
end

# Discrete univariate with finite support
function cdf(d::DiscreteUnivariateDistribution, x::Int, ::Type{FiniteSupport{true}})
    # calculate from left if x < (min + max)/2
    # (same as infinite support version)
    x <= div(minimum(d) + maximum(d),2) && return cdf(d, x, FiniteSupport{false})

    # otherwise, calculate from the right
    c = 1.0
    for y = x+1:maximum(d)
        c -= pdf(d, y)
    end
    return c
end

cdf(d::DiscreteUnivariateDistribution, x::Real) = cdf(d, floor(Int,x))

cdf(d::ContinuousUnivariateDistribution, x::Real) = throw(MethodError(cdf, (d, x)))

# ccdf

ccdf(d::DiscreteUnivariateDistribution, x::Int) = 1.0 - cdf(d, x)
ccdf(d::DiscreteUnivariateDistribution, x::Real) = ccdf(d, floor(Int,x))
ccdf(d::ContinuousUnivariateDistribution, x::Real) = 1.0 - cdf(d, x)

# logcdf

logcdf(d::DiscreteUnivariateDistribution, x::Int) = log(cdf(d, x))
logcdf(d::DiscreteUnivariateDistribution, x::Real) = logcdf(d, floor(Int,x))
logcdf(d::ContinuousUnivariateDistribution, x::Real) = log(cdf(d, x))

# logccdf

logccdf(d::DiscreteUnivariateDistribution, x::Int) = log(ccdf(d, x))
logccdf(d::DiscreteUnivariateDistribution, x::Real) = logccdf(d, floor(Int,x))

logccdf(d::ContinuousUnivariateDistribution, x::Real) = log(ccdf(d, x))

# quantile

quantile(d::UnivariateDistribution, p::Real) = throw(MethodError(quantile, (d, p)))

# cquantile

cquantile(d::UnivariateDistribution, p::Real) = quantile(d, 1.0 - p)

# invlogcdf

invlogcdf(d::UnivariateDistribution, lp::Real) = quantile(d, exp(lp))

# invlogccdf

invlogccdf(d::UnivariateDistribution, lp::Real) = quantile(d, -expm1(lp))

# gradlogpdf

gradlogpdf(d::ContinuousUnivariateDistribution, x::Real) = throw(MethodError(gradlogpdf, (d, x)))

# vectorized versions
for fun in [:pdf, :logpdf,
            :cdf, :logcdf,
            :ccdf, :logccdf,
            :invlogcdf, :invlogccdf,
            :quantile, :cquantile]

    _fun! = Symbol('_', fun, '!')
    fun! = Symbol(fun, '!')

    @eval begin
        function ($_fun!)(r::AbstractArray, d::UnivariateDistribution, X::AbstractArray)
            for i in 1 : length(X)
                r[i] = ($fun)(d, X[i])
            end
            return r
        end

        function ($fun!)(r::AbstractArray, d::UnivariateDistribution, X::AbstractArray)
            length(r) == length(X) ||
                throw(ArgumentError("Inconsistent array dimensions."))
            $(_fun!)(r, d, X)
        end

        ($fun)(d::UnivariateDistribution, X::AbstractArray) =
            $(_fun!)(Array(promote_type(partype(d), eltype(X)), size(X)), d, X)
    end
end

function _pdf_fill_outside!(r::AbstractArray, d::DiscreteUnivariateDistribution, X::UnitRange)
    vl = vfirst = first(X)
    vr = vlast = last(X)
    n = vlast - vfirst + 1
    if islowerbounded(d)
        lb = minimum(d)
        if vl < lb
            vl = lb
        end
    end
    if isupperbounded(d)
        ub = maximum(d)
        if vr > ub
            vr = ub
        end
    end

    # fill left part
    if vl > vfirst
        for i = 1:(vl - vfirst)
            r[i] = 0.0
        end
    end

    # fill central part: with non-zero pdf
    fm1 = vfirst - 1
    for v = vl:vr
        r[v - fm1] = pdf(d, v)
    end

    # fill right part
    if vr < vlast
        for i = (vr-vfirst+2):n
            r[i] = 0.0
        end
    end
    return vl, vr, vfirst, vlast
end

function _pdf!(r::AbstractArray, d::DiscreteUnivariateDistribution, X::UnitRange)
    vl,vr, vfirst, vlast = _pdf_fill_outside!(r, d, X)

    # fill central part: with non-zero pdf
    fm1 = vfirst - 1
    for v = vl:vr
        r[v - fm1] = pdf(d, v)
    end
    return r
end

abstract RecursiveProbabilityEvaluator

function _pdf!(r::AbstractArray, d::DiscreteUnivariateDistribution, X::UnitRange, rpe::RecursiveProbabilityEvaluator)
    vl,vr, vfirst, vlast = _pdf_fill_outside!(r, d, X)

    # fill central part: with non-zero pdf
    if vl <= vr
        fm1 = vfirst - 1
        r[vl - fm1] = pv = pdf(d, vl)
        for v = (vl+1):vr
            r[v - fm1] = pv = nextpdf(rpe, pv, v)
        end
    end

    return r
end

pdf(d::DiscreteUnivariateDistribution) = isbounded(d) ? pdf(d, minimum(d):maximum(d)) :
                                                        error("pdf(d) is not allowed when d is unbounded.")

## loglikelihood

function _loglikelihood(d::UnivariateDistribution, X::AbstractArray)
    ll = 0.0
    for i in 1:length(X)
        @inbounds ll += logpdf(d, X[i])
    end
    return ll
end

loglikelihood(d::UnivariateDistribution, X::AbstractArray) =
    _loglikelihood(d, X)

### macros to use StatsFuns for method implementation

macro _delegate_statsfuns(D, fpre, psyms...)
    dt = eval(D)
    T = dt <: DiscreteUnivariateDistribution ? :Int : :Real

    # function names from StatsFuns
    fpdf = Symbol(fpre, "pdf")
    flogpdf = Symbol(fpre, "logpdf")
    fcdf = Symbol(fpre, "cdf")
    fccdf = Symbol(fpre, "ccdf")
    flogcdf = Symbol(fpre, "logcdf")
    flogccdf = Symbol(fpre, "logccdf")
    finvcdf = Symbol(fpre, "invcdf")
    finvccdf = Symbol(fpre, "invccdf")
    finvlogcdf = Symbol(fpre, "invlogcdf")
    finvlogccdf = Symbol(fpre, "invlogccdf")

    # parameter fields
    pargs = [Expr(:(.), :d, Expr(:quote, s)) for s in psyms]

    esc(quote
        pdf(d::$D, x::$T) = $(fpdf)($(pargs...), x)
        logpdf(d::$D, x::$T) = $(flogpdf)($(pargs...), x)

        cdf(d::$D, x::$T) = $(fcdf)($(pargs...), x)
        ccdf(d::$D, x::$T) = $(fccdf)($(pargs...), x)
        logcdf(d::$D, x::$T) = $(flogcdf)($(pargs...), x)
        logccdf(d::$D, x::$T) = $(flogccdf)($(pargs...), x)

        quantile(d::$D, q::Real) = convert($T, $(finvcdf)($(pargs...), q))
        cquantile(d::$D, q::Real) = convert($T, $(finvccdf)($(pargs...), q))
        invlogcdf(d::$D, lq::Real) = convert($T, $(finvlogcdf)($(pargs...), lq))
        invlogccdf(d::$D, lq::Real) = convert($T, $(finvlogccdf)($(pargs...), lq))
    end)
end

##### specific distributions #####

const discrete_distributions = [
    "bernoulli",
    "betabinomial",
    "binomial",
    "categorical",
    "discreteuniform",
    "geometric",
    "hypergeometric",
    "negativebinomial",
    "noncentralhypergeometric",
    "poisson",
    "skellam",
    "poissonbinomial"
]

const continuous_distributions = [
    "arcsine",
    "beta",
    "betaprime",
    "biweight",
    "cauchy",
    "chisq",    # Chi depends on Chisq
    "chi",
    "cosine",
    "epanechnikov",
    "exponential",
    "fdist",
    "frechet",
    "gamma", "erlang",
    "generalizedpareto",
    "generalizedextremevalue",
    "gumbel",
    "inversegamma",
    "inversegaussian",
    "kolmogorov",
    "ksdist",
    "ksonesided",
    "laplace",
    "levy",
    "logistic",
    "noncentralbeta",
    "noncentralchisq",
    "noncentralf",
    "noncentralt",
    "normal",
    "normalcanon",
    "normalinversegaussian",
    "lognormal",    # LogNormal depends on Normal
    "pareto",
    "rayleigh",
    "symtriangular",
    "tdist",
    "triangular",
    "triweight",
    "uniform",
    "vonmises",
    "weibull"
]

for dname in discrete_distributions
    include(joinpath("univariate", "discrete", "$(dname).jl"))
end

for dname in continuous_distributions
    include(joinpath("univariate", "continuous", "$(dname).jl"))
end


text:
using Test
using LogicCircuits
using ProbabilisticCircuits

@testset "Probability of constraint" begin

    # two nodes
    simplevtree = zoo_vtree_file("simple2.vtree")
    simplepsdd = zoo_psdd_file("simple2.4.psdd")
    pc = read((simplepsdd, simplevtree), StructProbCircuit)

    @test pr_constraint(pc, pc) ≈ 1.0

    # Test with two psdds
    pc1files = (zoo_psdd_file("simple2.5.psdd"), simplevtree)
    pc2files = (zoo_psdd_file("simple2.6.psdd"), simplevtree)
    pc1 = read(pc1files, StructProbCircuit)
    pc2 = read(pc2files, StructProbCircuit)

    @test pr_constraint(pc1, pc2) ≈ 1

end


text:
module Constants

export astronomical_unit

const AU = 1.49597870700e8 # km

astronomical_unit(::Type{Float64}) = AU
astronomical_unit() = astronomical_unit(Float64)

end


text:
using Test

# the source code to be tested
include("../src/SBRG.jl")

# @testset "SUVector" begin
@testset "constructor" begin
    vs = (
        [1, 2, 5, 5, 7],
        [1, 2, 5, 7, 7],
        [1, 1, 2, 5, 7],
        [1, 5, 2, 5, 7],
        [1, 5, 2, 7],
        [2, 1, 5, 7],
        [1, 2, 7, 5],
    )
    for v in vs
        @test_throws AssertionError SUVector(v)
    end
end
@testset "indexing, push, append" begin
    v = SUVector([1, 5, 10])
    @test size(v) == (3,)
    @test v[2] == 5
    @test v[end] == 10
    @test (v[2] = 4) == 4
    @test_throws AssertionError v[2] = 11
    push!(v, 11)
    append!(v, [14, 18])
    @test v[4] == 11
    @test v[end] == 18
end
@testset "intersect" begin
    v1 = SUVector([1, 2, 3, 4])
    v2 = SUVector([1, 2])
    v3 = SUVector([2, 4])
    v4 = SUVector{Int}()
    args = ((v1, v2), (v1, v3), (v1, v4), (v4, v4))
    outs = (
        SUVector([1, 2]),
        SUVector([2, 4]),
        v4,
        v4,
    )
    for (arg, out) in zip(args, outs)
        @test all(intersect(arg...) .== out)
    end
end
@testset "symdiff" begin
    v1 = SUVector([2, 5, 7, 8, 10, 12, 15])
    v2 = SUVector([5, 8, 11, 12])
    v3 = SUVector{Int}()
    args = ((v1, v2), (v2, v3), (v3, v3))
    outs = (
        SUVector([2, 7, 10, 11, 15]),
        v2,
        v3,
    )
    for (arg, out) in zip(args, outs)
        @test all(symdiff(arg...) .== out)
    end
end
@testset "union" begin
    v1 = SUVector([1, 4, 6, 7])
    v2 = SUVector([2, 5, 6])
    v3 = SUVector([4, 7, 9])
    v4 = SUVector{Int}()
    args = ((v1, v2), (v1, v3), (v1, v4))
    outs = (
        SUVector([1, 2, 4, 5, 6, 7]),
        SUVector([1, 4, 6, 7, 9]),
        v1,
    )
    # check for correctness
    for (arg, out) in zip(args, outs)
        @test all(union(arg...) .== out)
    end
end
# end

@testset "Pauli & Term" begin
    # set up test arguments and outputs
    X1Y2 = Pauli(1, [1, 2], [2])
    X2Y3 = Pauli(1, [2, 3], [3])
    X1X2 = Pauli(0, [1, 2], Int64[])
    Z1Z2 = Pauli(0, Int64[], [1, 2])
    Y1Y2 = Pauli(0, [1, 2], [1, 2])
    T1 = Term(2, X1Y2)
    T2 = Term(3, X2Y3)
    # set up outputs
    X1Y2X2Y3 = Pauli(0, [1, 3], [2, 3])
    X2Y3X1Y2 = Pauli(2, [1, 3], [2, 3])
    T1T2 = Term(6, X1Y2X2Y3)
    # check for correctness
    @test X1Y2 == X1Y2
    @test X1Y2 != X2Y3
    @test X1Y2 * X2Y3 == X1Y2X2Y3
    @test X2Y3 * X1Y2 == X2Y3X1Y2
    @test !commute(X1Y2, X2Y3)
    @test commute(X1X2, Z1Z2)
    @test commute(X1X2, Y1Y2)
    @test commute(Y1Y2, Z1Z2)
    @test_throws AssertionError Pauli(5, Int64[], Int64[])
    @test_throws AssertionError Pauli(2, [2, 1], [1])
    @test all(sites(X1Y2) .== [1, 2])
    @test center(X1Y2) == 1.5
    @test T1*T2 == T1T2
end

@testset "RGHamiltonian" begin
    # some pauli operators on three sites
    X1X2 = Pauli(0, [1,2], Int[]);
    X2X3 = Pauli(0, [2,3], Int[]);
    Z1 = Pauli(0, Int[], [1])
    Z2 = Pauli(0, Int[], [2])
    Z3 = Pauli(0, Int[], [3])
    Os = (X1X2, X2X3, Z1, Z2, Z3)
    # coefficients
    hs = rand(length(Os))
    # hamiltonian
    H = RGHamiltonian([Term(h, O) for (h, O) in zip(hs, Os)])
end

max_stars_repo_name: aormorningstar/SBRG

max_issues_repo_name: aormorningstar/SBRG

max_forks_repo_name: aormorningstar/SBRG


text:
const Δt_sec = 86400.0
const ldd_mv = 255

# read the staticmaps into memory
nc = NCDataset(staticmaps_rhine_path)
# helper function to get the axis order and directionality right
read_right(nc, var) = reverse(permutedims(nc[var][:]); dims = 2)
ldd_2d = read_right(nc, "ldd")

inds, _ = Wflow.active_indices(ldd_2d, ldd_mv)
n = length(inds)

# take out only the active cells
ldd = ldd_2d[inds]
slope = read_right(nc, "slope")[inds]
N = read_right(nc, "N")[inds]
Qold = read_right(nc, "Qold")[inds]
Bw = read_right(nc, "Bw")[inds]
waterlevel = read_right(nc, "waterlevel")[inds]
DCL = read_right(nc, "DCL")[inds]
close(nc)

# create the directed acyclic graph from the drainage direction array
graph = Wflow.flowgraph(ldd, inds, Wflow.pcr_dir)
# a topological sort is used for visiting nodes in order from upstream to downstream
toposort = topological_sort_by_dfs(graph)
sink = toposort[end]
@test ldd[sink] == 5  # the most downstream node must be a sink

# calculate parameters of kinematic wave
const q = 0.000001
const β = 0.6
const AlpPow = (2.0 / 3.0) * β
AlpTermR = (N ./ sqrt.(slope)) .^ β
P = Bw + (2.0 * waterlevel)
α = AlpTermR .* P .^ AlpPow

Q = zeros(n)
Q = Wflow.kin_wave!(Q, graph, toposort, Qold, q, α, β, DCL, Δt_sec)

@testset "flow rate" begin
    @test sum(Q) ≈ 2.957806043289641e6
    @test Q[toposort[1]] ≈ 0.007260052312634069f0
    @test Q[toposort[n-100]] ≈ 3945.762718338739f0
    @test Q[sink] ≈ 4131.101474418251
end

@testset "kinematic wave subsurface flow" begin
    @test all(
        isapprox.(
            Wflow.kinematic_wave_ssf(
                210128378079.0733,
                215395179156.82645,
                1540.34273559,
                1.238,
                18021.0,
                0.25,
                0.346,
                0.0017669756,
                1800.0,
                1.0,
                1697.05 * 1000.0,
                1200.0 * 1000.0,
                2443723.716252628,
            ),
            (7.410313985168225e10, 1540.1496836278836, -0.0),
        ),
    )
end

@testset "accucapacity" begin
    # test based on a subset of the examples at
    # https://pcraster.geo.uu.nl/pcraster/4.3.0/documentation/pcraster_manual/sphinx/op_accucapacity.html#examples
    # of the node at (row 3, column 2) and upstream nodes
    g = DiGraph(6)
    add_edge!(g, 1, 4)
    add_edge!(g, 2, 5)
    add_edge!(g, 3, 5)
    add_edge!(g, 4, 6)
    add_edge!(g, 5, 6)
    network = (graph = g, order = [1, 2, 3, 4, 5, 6])

    # example 1, accucapacityflux
    material = Float64[0.5, 2, 2, 0.5, 2, 0.5]
    capacity = fill(1.5, 6)
    flux, new_material = Wflow.accucapacityflux(material, network, capacity)
    @test new_material != material
    @test new_material == [0.0, 0.5, 0.5, 0.0, 3.5, 1.5]
    @test flux == Float64[0.5, 1.5, 1.5, 1, 1.5, 1.5]

    # example 2, accucapacityflux
    material = fill(10.0, 6)
    capacity = Float64[2, 30, 30, 2, 30, 2]
    flux, new_material = Wflow.accucapacityflux(material, network, capacity)
    @test new_material != material
    @test new_material == [8.0, 0.0, 0.0, 10.0, 0.0, 40.0]
    @test flux == Float64[2, 10, 10, 2, 30, 2]

    # example 1, accucapacitystate
    material = Float64[0.5, 2, 2, 0.5, 2, 0.5]
    capacity = fill(1.5, 6)
    new_material = Wflow.accucapacitystate(material, network, capacity)
    @test new_material != material
    @test new_material == [0.0, 0.5, 0.5, 0.0, 3.5, 1.5]

    # example 2, accucapacitystate
    material = fill(10.0, 6)
    capacity = Float64[2, 30, 30, 2, 30, 2]
    new_material = Wflow.accucapacitystate(material, network, capacity)
    @test new_material != material
    @test new_material == [8.0, 0.0, 0.0, 10.0, 0.0, 40.0]
end

@testset "local inertial long channel MacDonald (1997)" begin

    g = 9.80665
    L = 1000.0
    dx = 5.0
    n = Int(L / dx)

    # analytical solution MacDonald (1997) for channel with length L of 1000.0 m, Manning's
    # n of 0.03, constant inflow of 20.0 m3/s at upper boundary and channel width of 10.0 m
    # water depth profile h(x)
    h(x) = (4 / g)^(1.0 / 3.0) * (1.0 + 0.5 * exp(-16.0 * (x / L - 0.5)^2.0))
    # spatial derivative of h(x)
    h_acc(x) =
        -(4 / g)^(1.0 / 3.0) * 16.0 / L * (x / L - 0.5) * exp(-16 * (x / L - 0.5)^2.0)
    # solution for channel slope s(x)
    s(x) =
        (1.0 - 4.0 / (g * h(x)^(3.0))) * h_acc(x) +
        0.36 * (2 * h(x) + 10.0)^(4.0 / 3.0) / ((10.0 * h(x))^(10.0 / 3.0))

    h_a = h.([dx:dx:L;]) # water depth profile (analytical solution)
    # integrate slope to get elevation (bed level) z
    x = [dx:dx:L;]
    zb = first.([quadgk(s, xi, L, rtol = 1e-12) for xi in x])

    # initialize ShallowWaterRiver
    graph = DiGraph(n)
    for i = 1:n
        add_edge!(graph, i, i + 1)
    end

    dl = fill(dx, n)
    width = fill(10.0, n)
    n_river = fill(0.03, n)

    # for each link the src and dst node is required
    nodes_at_link = Wflow.adjacent_nodes_at_link(graph)
    _ne = ne(graph)

    # determine z, width, length and manning's n at links
    zb_max = fill(0.0, _ne)
    width_at_link = fill(0.0, _ne)
    length_at_link = fill(0.0, _ne)
    mannings_n = fill(0.0, _ne)
    for i = 1:_ne
        zb_max[i] = max(zb[nodes_at_link.src[i]], zb[nodes_at_link.dst[i]])
        width_at_link[i] = min(width[nodes_at_link.dst[i]], width[nodes_at_link.src[i]])
        length_at_link[i] = 0.5 * (dl[nodes_at_link.dst[i]] + dl[nodes_at_link.src[i]])
        mannings_n[i] =
            (
                n_river[nodes_at_link.dst[i]] * dl[nodes_at_link.dst[i]] +
                n_river[nodes_at_link.src[i]] * dl[nodes_at_link.src[i]]
            ) / (dl[nodes_at_link.dst[i]] + dl[nodes_at_link.src[i]])
    end

    network = (
        nodes_at_link = nodes_at_link,
        links_at_node = Wflow.adjacent_links_at_node(graph, nodes_at_link),
    )

    alpha = 0.7
    Δt = 1.0
    h_thresh = 1.0e-03
    froude_limit = true
    h_init = zeros(n - 1)
    push!(h_init, h_a[n])

    sw_river = Wflow.ShallowWaterRiver(
        n = n - 1,
        ne = _ne,
        g = 9.80665,
        α = alpha,
        h_thresh = h_thresh,
        Δt = Δt,
        q0 = zeros(_ne),
        q = zeros(_ne),
        q_av = zeros(_ne),
        zb_max = zb_max,
        mannings_n = mannings_n,
        h = h_init,
        η_max = zeros(_ne),
        hf = zeros(_ne),
        h_av = zeros(n),
        width = width,
        width_at_link = width_at_link,
        a = zeros(_ne),
        r = zeros(_ne),
        volume = fill(0.0, n),
        error = zeros(n),
        inflow = zeros(n),
        inwater = zeros(n),
        dl = dl,
        dl_at_link = length_at_link,
        bankfull_volume = fill(Wflow.mv, n),
        bankfull_depth = fill(Wflow.mv, n),
        zb = zb,
        froude_limit = froude_limit,
        reservoir_index = zeros(n),
        lake_index = zeros(n),
        reservoir = nothing,
        lake = nothing,
    )

    # run until steady state is reached
    ϵ = 1.0e-12
    while true
        sw_river.inwater[1] = 20.0
        h0 = mean(sw_river.h)
        Δt = Wflow.stable_timestep(sw_river)
        Wflow.shallowwater_river_update(sw_river, network, Δt, nothing, 0.0, true)
        d = abs(h0 - mean(sw_river.h))
        if d <= ϵ
            break
        end
    end

    # test for mean absolute error [cm]
    @test mean(abs.(sw_river.h .- h_a)) * 100.0 ≈ 1.873574206931199

end


text:
#!/usr/local/bin/julia
# LSL.jl: Julia interface for Lab Streaming Layer
# Copyright (C) 2019 Samuel Powell

# lib_send_data.jl: replicate SendDataC.c LSL application example

using LSL

function main(name, uid)

  # Create new stream info
  info = LSL.lsl_create_streaminfo(name, "EEG", 8, 500, LSL.cft_float32, uid);

  # Add meta-data fields to stream info
  desc = LSL.lsl_get_desc(info)
  LSL.lsl_append_child_value(desc, "manufacturer", "LSL")
  channels = ["C3","C4","Cz","FPz","POz","CPz","O1","O2"]
  chns = LSL.lsl_append_child(desc,"channels")
  for c in 1:8
    chn = LSL.lsl_append_child(chns,"channel")
    LSL.lsl_append_child_value(chn,"label",channels[c])
    LSL.lsl_append_child_value(chn,"unit","microvolts")
    LSL.lsl_append_child_value(chn,"type","EEG")
  end

  # Create a new outlet (chunking: default, buffering: 360 seconds)
  outlet = LSL.lsl_create_outlet(info,0,360);

  println("Waiting for consumers")
  while LSL.lsl_wait_for_consumers(outlet, 120) == 0
    sleep(0.1)
  end

  println("Now sending data...")

  # Send data until the last consumer has disconnected
  t = 0
  cursample = zeros(Float32, 8)

  while LSL.lsl_have_consumers(outlet) != 0
    cursample[1] = t;
    for c in 2:8
      cursample[c] = Float32.((rand()%1500)/500.0-1.5);
    end
    LSL.lsl_push_sample_f(outlet,cursample);
    t += 1
  end

  println("Lost the last consumer, shutting down")
  LSL.lsl_destroy_outlet(outlet);

end

println("lib_send_data example program: sends 8 float channels as fast as possible")
println("Usage: lib_send_data.jl [streamname] [streamuid]")
println("Using lsl $(LSL.lsl_library_version()), lsl_library info: $(unsafe_string(LSL.lsl_library_info()))")

const name = length(ARGS) > 0 ? ARGS[1] : "send_data_lib"
const uid = length(ARGS) > 1 ? ARGS[2] : "325wqer4354"

main(name, uid)


text:
module IsDef
export isdef, Out, NotApplicable, ∨, apply

using Compat
import InteractiveUtils

"""
just applies a given function to arguments and keyword arguments

This little helper is crucial if you want to typeinfer
when only knowing the function type instead of the function instance.
"""
apply(f, args...; kwargs...) = f(args...; kwargs...)

# alias for promote_type to deal with types more compact
# we choose the symbol for join `∨` because promote_type is kind of a maximum (in type-hierarchy, with Any being the top)
# see https://en.wikipedia.org/wiki/Join_and_meet
@doc raw"""
    TypeA ∨ TypeB = promote_type(TypeA, TypeB)
    ∨(values...) = ∨(typeof.(values)...)

`∨` (latex `\vee`) is alias for `promote_type`.

When called on values, the values will be cast to types via use of `typeof` for convenience.
"""
function ∨ end
T1::Type ∨ T2::Type = promote_type(T1, T2)
∨(Ts::Type...) = promote_type(Ts...)
val1 ∨ val2 = typeof(val1) ∨ typeof(val2)
∨(values...) = ∨(typeof.(values)...)

# Core Interface
# ==============

struct NotApplicable end
# this will be overwritten below for certain types like Any
isapplicable(type) = !(NotApplicable <: type)
# despite NotApplicable <: Any, we interprete Any always as applicable
isapplicable(::Type{Any}) = true

"""
    isdef(func, ArgType1, ArgType2, ...)::Bool
    isdef(f, args...) = isdef(f, typeof.(args)...)

Checks whether the function is defined for the actual types or not.
This works in compile time and hence can be used to optimize code.

When called on values, the values will be cast to types via use of `typeof` for convenience.

!!! warning "CAUTION"
    If `isdef(...) == true`, still a MethodError might happen at runtime. This is due to incomplete type inference.

    SOLUTION: Overload `IsDef.return_type` if you experience unexpected behaviour for your types

    For instance to say that some call like `myfunc(::Int, ::String)` is not defined, then define the following
    ```julia
    function IsDef.return_type(::Type{Tuple{typeof(myfunc), TypeArg1, TypeArg2}})
      Union{}  # return empty Union to indicate something is not defined
    end
    ```

Semantics
---------
The implementation follows a closed-world semantics in that `isdef` will return true
for a given abstract type, if the function is defined for all currently defined leaf types of that
abstract type.

The only exceptions currently are `Any`, `Function`,  and `Exception`, for which `isdef` returns `true`.
"""
function isdef(f, types::Vararg{<:Type})
  signature_type = Tuple{typeof(f), types...}
  isapplicable(return_type(signature_type))
end
# we need to handle Type constructors specifically, as the typeof function forgets
function isdef(f::Type{T}, types::Vararg{<:Type}) where T
  signature_type = Tuple{Type{T}, types...}
  isapplicable(return_type(signature_type))
end
function isdef(::typeof(apply), types::Vararg{<:Type})
  signature_type = Tuple{types...}
  isapplicable(return_type(signature_type))
end
isdef(f, args...) = isdef(f, typeof.(args)...)

"""
    Out(func, ArgType1, ArgType2, ...)::ReturnType
    Out(f, args...) = isdef(f, typeof.(args)...)

Returns outputtype of function application. Returns `Traits.NotApplicable` if compiler notices that no Method can be
found.

When called on values, the values will be cast to types via use of `typeof` for convenience.

!!! warning "CAUTION"
    If `Out(...) == Any`, still a MethodError might happen at runtime. This is due to incomplete type inference.

    SOLUTION: Overload `IsDef.return_type` if you experience unexpected behaviour for your types

    For instance to say that some call like `myfunc(::Int, ::String)` is not defined, then define the following
    ```julia
    function IsDef.return_type(::Type{Tuple{typeof(myfunc), TypeArg1, TypeArg2}})
      Union{}  # return empty Union to indicate something is not defined
    end
    ```

UnionAll types and abstract types are concretified to the Union of their existing subtypes, in order
to improve type-inference. The only exceptions so far are `Any`, `Function` and `Exception`, as they have way too
many subtypes to be of practical use.
"""
function Out(f, types::Vararg{<:Type})
  signature_type = Tuple{typeof(f), types...}
  return_type(signature_type)
end
# we need to handle Type constructors specifically, as the typeof function forgets
function Out(f::Type{T}, types::Vararg{<:Type}) where T
  signature_type = Tuple{Type{T}, types...}
  return_type(signature_type)
end
function Out(::typeof(apply), types::Vararg{<:Type})
  signature_type = Tuple{types...}
  return_type(signature_type)
end
Out(f, args...) = Out(f, typeof.(args)...)

# Core Background
# ===============

"""
  wrapper arround Julia's type inference

This should be overloaded if you want to fix certain wrong typeinferences for your custom types.
It is used internally by both `isdef` and `Out`.

Returning `Union{}` is interpreted as `MethodError`.
"""
return_type(::Type{Ts}) where {Ts <: Tuple} = core_return_type(abstracttype_to_list_of_concretetypes(Ts))
# TODO we tried making this function the @generated one (and not abstracttype_to_list_of_concretetypes),
# however we immediately get NotApplicable for the very first test
# @test Out(Base.map, typeof(x->2x), Vector{Int}) == Vector{Int}

"""
`core_return_type` can either except
- a Union (which may resolved type-wise)
- or a Tuple-Type (which will be the final call signature)

When something fails, it won't return `Union{}` (or `Array{Union{}, 1}`), but `NotApplicable`.
This also simplifies the Union part, as Union{} would simply disappear as soon as put into Union{args...}.
A `NotApplicable` will stay.
"""
function core_return_type(T::Union)
  Union{core_return_type(T.a), core_return_type(T.b)}
end

# NOTE we cannot dispatch on Type{<:Tuple}, as this is more specific than Union, and we only have Union of Tuple,
# hence the Union clause would never be called
# TODO we need to add another function layer as otherwise typeinference again fails...
# The test which failed in all cases is `@test Out(Out, typeof(Base.map), typeof(x->2x), Vector{Int}) == Type{Vector{Int}}`
# TODO we need to forward Tuple{T} and cannot use Tuple{T} in the same function as Core.Compiler.return_type
# Otherwise the meta test fails
# Expression: Out(Out, typeof(Base.map), typeof((x->2x)), Vector{Int}) == Type{Vector{Int}}
# Evaluated: DataType == Type{Array{Int64,1}}
function core_return_type(T::Type)
  _core_return_type(Tuple{T})
end

# TODO for some unclarified and not yet minimized reason, the following function signature does not infere well
# `_core_return_type(::Type{T}) where T` and hence neither `_core_return_type(::Type{T}) where T <: Tuple`
# apparently the Type{T} looses some information given by the value, by restricting to type-level only.
# TODO restricting type to `Type{<:Tuple}` somehow looses typeinformation similarly
# The test which failed in all cases is `@test Out(Out, typeof(Base.map), typeof(x->2x), Vector{Int}) == Type{Vector{Int}}`
function _core_return_type(T::Type)
  # TODO we need to workaround some surprising type inference issues. See https://github.com/JuliaLang/julia/issues/36626
  apply_internal(full_call) = full_call[1](Base.tail(full_call)...)
  emptyunion_to_notapplicable(Core.Compiler.return_type(apply_internal, T))
end

# check whether a type means that there was an error
# ==================================================

emptyunion_to_notapplicable(::Type{Union{}}) = NotApplicable
@generated function emptyunion_to_notapplicable(::Type{T}) where T
  _emptyunion_to_notapplicable(T)
end
_emptyunion_to_notapplicable(T) = isbottom(T) ? NotApplicable : T

isbottom(T::Type{Union{}}) = true
isbottom(value) = false
@generated function isbottom(::Type{T}) where T
  _isbottom(T)
end

function _isbottom(T::DataType)
  any(T.parameters) do p
    isbottom(p)
  end
end
_isbottom(T::UnionAll) = _isbottom(Base.unwrap_unionall(T))
_isbottom(T::Union) = _isbottom(T.a) || _isbottom(T.b)

# Helpers
# =======

# Some Types need to stay abstract because they have way too many subtypes
const DontTouchTypes = (Any, Function, Exception)

for type in DontTouchTypes
  eval(quote
    Type2Union(T::Type{$type}) = $type
    leaftypes(::Type{$type}) = [$type]
    allsubtypes(::Type{$type}) = [$type]
  end)
end

# Tuples{Union{...}} -> Union{Tuple{...}}
# however we represent it as a tuple instead of a Union because this gives us better type inference
# (we just tried it out to get to know)

# we only generate the tuple version, as this is the entry point we need
# all other Type2Union are just for recursion and nested types
# but also bringing them @generated does not add much, hence we follow recommendations and leave them as functions
"""
takes a type like
```
Tuple{Int, AbstractFloat}
```
and returns Union
```
Union{
  Tuple{Int, BigFloat},
  Tuple{Int, Float16},
  Tuple{Int, Float32},
  Tuple{Int, Float64},
}
```
"""
@generated function abstracttype_to_list_of_concretetypes(::Type{T}) where T <: Tuple
  _abstracttype_to_list_of_concretetypes(T)
end
_abstracttype_to_list_of_concretetypes(T) = Type2Union(T)

function Type2Union(::Type{T}) where T <: Tuple
  alltypes = Base.uniontypes.(Type2Union.(T.parameters))
  combinations = collect(Iterators.product(alltypes...))
  Union{map(c -> Tuple{c...}, combinations)...}
end

# Other Types can be mapped to the Union of their concrete subtypes
Type2Union(T::Union) = Union{Type2Union(T.a), Type2Union(T.b)}
function Type2Union(::Type{T}) where T
  leaftypes_with_abstract_typevariables = leaftypes(T)
  # apparently the current typeinference already works super well if instead of a typevariable-upperbound
  # `Vector{<: Number}`
  # we just use the union
  # `Union{Vector{subtype1}, Vector{subtype2}, ... for subtype in allsubtypes(Number)}`
  plain_leaftypes = unionall_to_union.(leaftypes_with_abstract_typevariables)
  # TODO apparently julia's type-inference stops after a Union of three types with the approximate result `Any`
  # that is of course not optimal...
  # I guess best is to wait for a better solution and just overload the return_type function with your specific needs
  Union{plain_leaftypes...}
end

# TODO is it possible to combine Type2Union and unionall_to_union together? looks like it could be possible
unionall_to_union(type) = unionall_to_union(type, IdDict())
unionall_to_union(any, dict) = any  # bytetypes in parameters
function unionall_to_union(typevar::TypeVar, dict::IdDict)
  # as soon as we visit a typevariable WITHIN the typeparameters (not within the wheres)
  # this should already be mapped in the dict
  dict[typevar]
end
function unionall_to_union(type::Union, dict::IdDict)
  Union{unionall_to_union(type.a, dict), unionall_to_union(type.b, dict)}
end
function unionall_to_union(type::DataType, dict::IdDict)
  if isempty(type.parameters)
    type
  else
    unionall = type.name.wrapper
    newparameters = [unionall_to_union(p, dict) for p in type.parameters]
    unionall{newparameters...}
  end
end
function unionall_to_union(type::UnionAll, dict::IdDict)
  subtypes = allsubtypes(type.var.ub)
  function recurse_with_new_dict(subtype)
    newdict = copy(dict)
    if subtype in DontTouchTypes
      # in case of Any, we do not do anything practically
      newdict[type.var] = type.var
      newbody = unionall_to_union(type.body, newdict)
      _UnionAll(type.var, newbody) # rewrap into UnionAll
    else
      newdict[type.var] = subtype
      unionall_to_union(type.body, newdict)
      # here the type.var was replaced by an actual type, hence the result is already a valid type and does not need to be rewrapped into UnionALl
    end
  end
  Union{recurse_with_new_dict.(subtypes)...}
end

# UnionAll of Unions have buggy behaviour currently
# see https://github.com/JuliaLang/julia/issues/35910 for updates
_UnionAll(typevar, type::Union) = Union{_UnionAll(typevar, type.a), _UnionAll(typevar, type.b)}
_UnionAll(typevar, other) = UnionAll(typevar, other)

function subtypes(T)
  # InteractiveUtils.subtypes is super mighty in that it can already deal with UnionAll types
  sub = InteractiveUtils.subtypes(T)
  if length(sub) == 1 && sub[1] === T
    # InteractiveUtils has a small bug that sometimes the  type is itself returned as a subtype
    # which is in theory true, however leads to infinite recursion
    []
  else
    sub
  end
end

function leaftypes(T)
  sub = subtypes(T)
  if isempty(sub)
    [T]
  else
    vcat((leaftypes(S) for S in sub)...)
  end
end

function allsubtypes(T)
  sub = subtypes(T)
  if isempty(sub)
    [T]
  else
    vcat(T, (allsubtypes(S) for S in sub)...)
  end
end

"""
for some reasons the @generated macro triggers too early and leads Union{} for some types

calling this macro will redefine the @generated functions accordingly so that new TypeDefinitions are properly included
"""
macro redefine_generated()
  quote
    @generated function abstracttype_to_list_of_concretetypes(::Type{T}) where T <: Tuple
      _abstracttype_to_list_of_concretetypes(T)
    end
    @generated function IsDef.emptyunion_to_notapplicable(::Type{T}) where T
      _emptyunion_to_notapplicable(T)
    end
    @generated function IsDef.isbottom(::Type{T}) where T
      _isbottom(T)
    end
    nothing
  end
end

end # module


text:
# ---
# title: 1208. Get Equal Substrings Within Budget
# id: problem1208
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Sliding Window
# link: <https://leetcode.com/problems/get-equal-substrings-within-budget/description/>
# hidden: true
# ---
#
# You are given two strings `s` and `t` of the same length. You want to change
# `s` to `t`. Changing the `i`-th character of `s` to `i`-th character of `t`
# costs `|s[i] - t[i]|` that is, the absolute difference between the ASCII
# values of the characters.
#
# You are also given an integer `maxCost`.
#
# Return the maximum length of a substring of `s` that can be changed to be the
# same as the corresponding substring of `t`with a cost less than or equal to
# `maxCost`.
#
# If there is no substring from `s` that can be changed to its corresponding
# substring from `t`, return `0`.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "abcd", t = "bcdf", maxCost = 3
#     Output: 3
#     Explanation: "abc" of s can change to "bcd". That costs 3, so the maximum length is 3.
#
# **Example 2:**
#
#
#
#     Input: s = "abcd", t = "cdef", maxCost = 3
#     Output: 1
#     Explanation: Each character in s costs 2 to change to charactor in t, so the maximum length is 1.
#
#
# **Example 3:**
#
#
#
#     Input: s = "abcd", t = "acde", maxCost = 0
#     Output: 1
#     Explanation: You can't make any change, so the maximum length is 1.
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length, t.length <= 10^5`
#   * `0 <= maxCost <= 10^6`
#   * `s` and `t` only contain lower case English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
x = parse(Int, readline())
println(x^3)


text:
using Econometrics, Statistics, Flux, Random, LinearAlgebra
using Base.Iterators
using BSON: @load
using BSON: @save

# unsmoothed QR objective
function QRobj(yhat, y, τ)
    ξ = y .- yhat
    sum(max.(τ.*ξ,(τ .-1.0).*ξ))
end
# smoothed QR objective
function SmoothQRobj(yhat, y, τ, α=0.02)
    ξ = y .- yhat
    sum(τ.*ξ .+ α.*log.(1.0 .+ exp.(-1.0 .* ξ./α)))
end

function Train()
    @load "cooked_data.bson" params statistics
    params = Float32.(params)
    statistics = Float32.(statistics)
    S = TrainingTestingSize # number of draws from prior
    trainsize = Int(TrainingProportion*S)
    yin = params[1:trainsize, :]'
    yout = params[trainsize+1:end, :]'
    xin = statistics[1:trainsize, :]'
    xout = statistics[trainsize+1:end, :]'
    # model
    nStats = size(xin,1)
    model = Chain(
        Dense(nStats,5*nParams, tanh),
        Dense(5*nParams, nParams)
    )
    θ = Flux.params(model)
    opt = AdaMax()
    # weight by inverse std. dev. of params, to put equal weight
    s = Float32.(std(params,dims=1)')
    # Select the objective here: ordinary regression or (smoothed) quantile regression
    #loss(x,y) = sqrt.(Flux.mse(model(x)./s,y./s))
    τ = 0.5  # choose the quantile you want
    #loss(x,y) = QRobj(model(x),y, τ)
    loss(x,y) = SmoothQRobj(model(x),y, τ)
    function monitor(e)
        println("epoch $(lpad(e, 4)): (training) loss = $(round(loss(xin,yin); digits=4)) (testing) loss = $(round(loss(xout,yout); digits=4))| ")
    end
    bestsofar = 1.0e10
    pred = 0.0 # define it here to have it outside the for loop
    batches = [(xin[:,ind],yin[:,ind])  for ind in partition(1:size(yin,2), 32)];
    for i = 1:Epochs
        Flux.train!(loss, θ, batches, opt)
        current = loss(xout,yout)
        if current < bestsofar
            bestsofar = current
            @save "best.bson" model
            xx = xout
            yy = yout
            println("________________________________________________________________________________________________")
            monitor(i)
            pred = model(xx) # map pred to param space
            error = yy .- pred
            results = [pred;error]
            rmse = sqrt.(mean(error.^Float32(2.0),dims=2))
            println(" ")
            println("RMSE for model parameters ")
            prettyprint(reshape(round.(rmse,digits=3),1,3))
            println(" ")
            println("dstats prediction of parameters:")
            dstats(pred');
            println(" ")
            println("dstats prediction error:")
            dstats(error');
        end
    end
    return nothing
end


text:
using SeaPearl

"""
Proposition for using SeaPearl CP solver for the all interval series problem
Reference https://www.csplib.org/Problems/prob007/
"""

struct OutputDataAIS
    nb_sols ::Int
    indices ::Matrix{Int}
    diffs ::Matrix{Int}
end

"""
    model_AIS(n::Int; benchmark=false, variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())

return the SeaPearl model for to the All-interval_Series (AIS) problem, using SeaPearl.MinDomainVariableSelection heuristique
and  SeaPearl.AllDifferent  and SeaPearl.Asolute (without solving it)

# Arguments
- `n::Int`: dimension
- 'variableSelection': SeaPearl variable selection. By default: SeaPearl.MinDomainVariableSelection{false}()
- 'valueSelection': SeaPearl value selection. By default: =SeaPearl.BasicHeuristic()
"""
function model_AIS(n::Int; benchmark=false, variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())
    trailer = SeaPearl.Trailer()
    model = SeaPearl.CPModel(trailer)

    series = Vector{SeaPearl.AbstractIntVar}(undef, n)
    series[1] = SeaPearl.IntVar(0, n-1, "s_"*string(1), trailer)
    SeaPearl.addVariable!(model, series[1]; branchable=true)

    interval_vector = Vector{SeaPearl.AbstractIntVar}(undef, n-1) #differences
    interval_vector2 = Vector{SeaPearl.AbstractIntVar}(undef, n-1) #Absolute of the differences

    for i = 1:n-1
        series[i+1] = SeaPearl.IntVar(0, n-1, "s_"*string(i+1), trailer)
        SeaPearl.addVariable!(model, series[i+1]; branchable=true)

        interval_vector[i] = SeaPearl.IntVar(-n+1, n-1, "v1_"*string(i), trailer)

        minus = SeaPearl.IntVarViewOpposite(series[i+1], "-s["*string(i+1)*"]")
        vars = SeaPearl.AbstractIntVar[interval_vector[i], minus, series[i]]
        transition = SeaPearl.SumToZero(vars, trailer)
        push!(model.constraints, transition)

        interval_vector2[i] = SeaPearl.IntVar(1, n-1, "v2_"*string(i), trailer)

        abs_constraint = SeaPearl.Absolute(interval_vector[i], interval_vector2[i], trailer)
        push!(model.constraints, abs_constraint)
    end

    push!(model.constraints, SeaPearl.AllDifferent(series, trailer))
    push!(model.constraints, SeaPearl.AllDifferent(interval_vector2, trailer))

    return model
end

"""
    solve_AIS(n::Int; benchmark=false, variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())

Solve the SeaPearl model for to the All-interval_Series (AIS) problem, using SeaPearl.MinDomainVariableSelection heuristique
and  SeaPearl.AllDifferent and SeaPearl.Absolute, and the function model_AIS

# Arguments
- `n::Int`: dimension
- 'variableSelection': SeaPearl variable selection. By default: SeaPearl.MinDomainVariableSelection{false}()
- 'valueSelection': SeaPearl value selection. By default: =SeaPearl.BasicHeuristic()
"""
function solve_AIS(n::Int; benchmark=false, variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())
    model = model_AIS(n; benchmark=false, variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())
    status = @time SeaPearl.solve!(model; variableHeuristic=variableSelection, valueSelection=valueSelection)
    return model
end

"""
    solve_AIS(model::SeaPearl.CPModel)

Solve the SeaPearl model for to the AIS problem, using an existing model (from model_AIS)

# Arguments
- `model::SeaPearl.CPModel`: model (from model_AIS)

"""
function solve_AIS(model::SeaPearl.CPModel)
    status = @time SeaPearl.solve!(model; variableHeuristic=variableSelection, valueSelection=valueSelection)
    return model
end

"""
    outputFromSeaPearl(model::SeaPearl.CPModel; optimality=false)

Return the results of the AIS problem as type OutputDataAIS, taking the solved model (from solve_AIS) as argument

# Arguments
- `model::SeaPearl.CPModel`: needs the model to be already solved (by solve_AIS)
"""
function outputFromSeaPearl(model::SeaPearl.CPModel; optimality=false)
    solutions = model.solutions
    nb_sols = length(solutions)
    n = length(model.variables)
    indices = Matrix{Int}(undef, n, nb_sols)
    diffs = Matrix{Int}(undef, n-1, nb_sols)
    for (ind,sol) in enumerate(solutions)
        for i in 1:n
            indices[i, ind]= sol["s_"*string(i)]
            if (i<n) diffs[i,ind] = sol["s_"*string(i+1)]- sol["s_"*string(i)] end
        end
    end
    return OutputDataAIS(nb_sols, indices, diffs)
end

"""
outputFromSeaPearl(n::Int; benchmark=false, variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())

Shows the results of the AIS problem as type OutputDataAIS, taking dimension as argument

# Arguments
- `n::Int`: dimension
- 'variableSelection': SeaPearl variable selection. By default: SeaPearl.MinDomainVariableSelection{false}()
- 'valueSelection': SeaPearl value selection. By default: =SeaPearl.BasicHeuristic()
"""
function outputFromSeaPearl(n::Int; benchmark=false, variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())
    model = solve_AIS(n::Int; benchmark=false, variableSelection=SeaPearl.MinDomainVariableSelection{false}(), valueSelection=SeaPearl.BasicHeuristic())
    return outputFromSeaPearl(model)
end

max_stars_repo_name: corail-research/SeaPearlZoo

max_issues_repo_name: corail-research/SeaPearlZoo

max_forks_repo_name: corail-research/SeaPearlZoo


text:
#using Queryverse
using DataFrames
using DataFramesMeta
using Dates
using DelimitedFiles
using JuMP
using Logging
using CSV
using Query
using Statistics

const FILE_ACCESS_OVER = "w+";
const FILE_ACCESS_APPEND = "a+";
#include("import_post_data.jl")

DF_Gen_Sch = CSV.read(".//inputs//UnitSchedule_(DEC-DEP 2019).csv",
                  dateformat="mm-dd-yyyy THH:MM:SS.000Z", header=true, copycols=true, DataFrame);

DF_Gens = CSV.read(".//inputs//data_generators.csv", header=true, DataFrame);
DF_Peakers = CSV.read(".//inputs//data_peakers.csv", header=true, DataFrame);

select!(DF_Gens, Not(:Cogen));

#=
cycles_units_header    = ["Section", "Time", "Note1", "Note2", "Note3", "Note4"]
open(".//outputs//postproc_gens_cycle_units.csv", FILE_ACCESS_OVER) do io
    writedlm(io, permutedims(cycles_units_header), ',')
end; # closes file
=#

function calc_gens_cyclesbyID(df_results::DataFrame, df_gen_sch::DataFrame,
    df_gens_all::DataFrame, ini_date::String, end_date::String)
    #Selects a subset of schedules that fall into the months that users specify.
    # @where could be extended to include other month using the OR operation (||)
    ini_date = DateTime(ini_date, "mm/dd/yyyy HH:MM:SS");
    end_date = DateTime(end_date, "mm/dd/yyyy HH:MM:SS");

    # add UNIT_ID
    insertcols!(df_gen_sch, 4, :UNIT_ID => df_gen_sch.UNIT_NAME .* "_" .* string.(df_gen_sch.CC_KEY));
    #df_gen_sch = df_gen_sch

    df_gens_sch_drange = @from i in df_gen_sch begin
                @where DateTime(i.BEGIN_DATE) >= ini_date && DateTime(i.BEGIN_DATE) <= end_date && i.EDITION_NAME=="Actual"
                @select {i.PORTFOLIO_NAME, i.EDITION_NAME, i.UNIT_NAME, i.CC_KEY, i.UNIT_ID, i.UNIT_TYPE, i.BEGIN_DATE, i.MW}
                @collect DataFrame
           end

    df_convgen_sch = @from i in df_gens_sch_drange begin
                @where  i.UNIT_TYPE=="Steam" || i.UNIT_TYPE=="Steam Dual Boiler" || i.UNIT_TYPE=="Peaker" || i.UNIT_TYPE=="Nuclear" || i.UNIT_TYPE=="Combined-Cycle"
                @select {i.PORTFOLIO_NAME, i.EDITION_NAME, i.UNIT_NAME, i.CC_KEY, i.UNIT_ID, i.UNIT_TYPE, i.BEGIN_DATE, i.MW}
                @collect DataFrame
           end
    insertcols!(df_convgen_sch, size(df_convgen_sch,2), :On_Off => Int64.(df_convgen_sch.MW .> 0));
    insertcols!(df_convgen_sch, size(df_convgen_sch,2), :ShutDown => zeros(Int64));
    insertcols!(df_convgen_sch, size(df_convgen_sch,2), :StartUp => zeros(Int64));

    for unit in eachrow(df_gens_all)
        #unit_name = unit.UNIT_NAME;
        unit_id = unit.UNIT_ID;
        df_unit_sch = @where(df_convgen_sch, (:UNIT_ID .== unit_id))
        insertcols!(df_unit_sch, 1, :UNIT => unit.UNIT);
        #df_unit_sch = @where(df_convgen_sch, (:UNIT_NAME .== "ALLE_UN01"))

        sort!(df_unit_sch, :BEGIN_DATE)

        was_on = zero(Int64);
        for row in eachrow(df_unit_sch)
            if rownumber(row) == 1  #First row, save initial ON/OFF value
                was_on = row.On_Off;
                row.StartUp = row.On_Off;
            elseif ((row.On_Off == 1)  & (was_on == 0))
                was_on = 1;
                row.StartUp = 1;
            elseif ((row.On_Off == 0)  & (was_on == 1))
                was_on = 0;
                row.ShutDown = 1;
            end;
        end;

        if isempty(df_results)
            df_results = copy(df_unit_sch);
        else
            df_results = vcat(df_results, df_unit_sch);
        end;
        #println("Number of cycles for unit:", df_unit_sch.UNIT_NAME, "is: ", n_cycles)
    end; # for loop all units in dataframe df_gens_all

    return df_results;
end; # function

#Combines the slow generators and peakers specs into one dataframe
DF_Gens_All = vcat(DF_Gens, DF_Peakers);

# Filtering the Actual generation outcomes from UNIT_Schedules.csv
# Reformats the date-time data given in column BEGIN_DATE
DF_Gen_Sch.BEGIN_DATE = DateTime.(DF_Gen_Sch.BEGIN_DATE, "mm/dd/yyyy HH:MM:SS p");

Start_date = "01/01/2019 00:00:00"
Final_date = "01/02/2019 02:00:00"
DF_Results = DataFrame();

DF_Results = calc_gens_cycles(DF_Results, DF_Gen_Sch, DF_Gens_All, Start_date, Final_date);

Start_date = replace(Start_date, "/" => "-")
Final_date = replace(Final_date, "/" => "-")
Start_date = Start_date[1:findfirst(isequal(' '), Start_date)-1]
Final_date = Final_date[1:findfirst(isequal(' '), Final_date)-1]

CSV.write(".//outputs//postprocess//postproc_gens_cycles_ByID_$(Start_date)_$(Final_date).csv", DF_Results)

# Get Summary of results by unit and save them
#Reference: https://julia-data-query.readthedocs.io/en/latest/dplyr.html
DF_Results_ByID = groupby(DF_Results, [:UNIT_ID, :UNIT_TYPE]);

DF_ResultsID_Summary = @combine(DF_Results_ByID,
  Sum_On_Off = sum(:On_Off),
  Avg_On_Off = mean(:On_Off),
  Sum_StartUp = sum(:StartUp),
  Sum_ShutDown = sum(:ShutDown))

CSV.write(".//outputs//postprocess//postproc_gens_cycles_SUMMARYbyID_$(Start_date)_$(Final_date).csv", DF_ResultsID_Summary)

#=TODO:
Emissions calculations for 358 days
Shutdown is 1, one hour after the unit has been shutoff.
Compare the number of shutdowns and startup.
startup if at time 0 the units are on.
=#

#include(".//postprocessing//calc_input_gens_cyclesbyID.jl")

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
"""
```julia
StudentTLikelihood(ν::T,σ::Real=one(T))
```
[Student-t likelihood](https://en.wikipedia.org/wiki/Student%27s_t-distribution) for regression:
```math
    p(y|f,ν,σ) = Γ(0.5(ν+1))/(sqrt(νπ) σ Γ(0.5ν)) * (1+(y-f)^2/(σ^2ν))^(-0.5(ν+1))
```
`ν` is the number of degrees of freedom and `σ` is the variance for local scale of the data.

---

For the analytical solution, it is augmented via:
```math
    p(y|f,ω) = N(y|f,σ^2 ω)
```
Where `ω ~ IG(0.5ν,,0.5ν)` where `IG` is the inverse gamma distribution
See paper [Robust Gaussian Process Regression with a Student-t Likelihood](http://www.jmlr.org/papers/volume12/jylanki11a/jylanki11a.pdf)
"""
mutable struct StudentTLikelihood{T<:Real} <: RegressionLikelihood{T}
    ν::T
    α::T
    σ::T
    c::Vector{T}
    θ::Vector{T}
    function StudentTLikelihood{T}(ν::T,σ::T=one(T)) where {T<:Real}
        new{T}(ν,(ν+one(T))/2.0,σ)
    end
    function StudentTLikelihood{T}(ν::T,σ::T,c::AbstractVector{T},θ::AbstractVector{T}) where {T<:Real}
        new{T}(ν,(ν+one(T))/2.0,σ,c,θ)
    end
end

function StudentTLikelihood(ν::T,σ::T=one(T)) where {T<:Real}
    StudentTLikelihood{T}(ν,σ)
end

function init_likelihood(likelihood::StudentTLikelihood{T},inference::Inference{T},nLatent::Int,nSamplesUsed::Int,nFeatures::Int) where T
    if inference isa AnalyticVI || inference isa GibbsSampling
        StudentTLikelihood{T}(
        likelihood.ν,likelihood.σ,
        rand(T,nSamplesUsed),
        zeros(T,nSamplesUsed))
    else
        StudentTLikelihood{T}(likelihood.ν,likelihood.σ)
    end
end

function pdf(l::StudentTLikelihood,y::Real,f::Real)
    tdistpdf(l.ν,(y-f)/l.σ)
end

function Base.show(io::IO,model::StudentTLikelihood{T}) where T
    print(io,"Student-t likelihood")
end

function compute_proba(l::StudentTLikelihood{T},μ::AbstractVector{T},σ²::AbstractVector{T}) where {T<:Real}
    return μ,max.(σ²,zero(σ²)).+0.5*l.ν*l.σ^2/(0.5*l.ν-1)
end

## Local Updates ##

function local_updates!(l::StudentTLikelihood{T},y::AbstractVector,μ::AbstractVector,diag_cov::AbstractVector) where {T}
    l.c .= 0.5*(diag_cov+abs2.(μ-y).+l.σ^2*l.ν)
    l.θ = l.α./l.c
end

function sample_local!(l::StudentTLikelihood{T},y::AbstractVector,f::AbstractVector) where {T}
    l.c .= rand.(InverseGamma.(l.α,0.5*(abs2.(μ-y).+l.σ^2*l.ν)))
    set_ω!(l,inv.(l.c))
    return nothing
end

## Global Gradients ##

@inline ∇E_μ(l::StudentTLikelihood,::AOptimizer,y::AbstractVector) = (l.θ.*y,)
@inline ∇E_Σ(l::StudentTLikelihood,::AOptimizer,y::AbstractVector) = (0.5.*l.θ,)

## ELBO Section ##

function expec_log_likelihood(l::StudentTLikelihood{T},i::AnalyticVI,y::AbstractVector,μ::AbstractVector,diag_cov::AbstractVector) where {T}
    tot = -0.5*length(y)*(log(twoπ*l.σ^2))
    tot += -sum(log.(l.c).-digamma(l.α))
    tot += -0.5*(dot(l.θ,diag_cov)+dot(l.θ,abs2.(μ))-2.0*dot(l.θ,μ.*y)+dot(l.θ,abs2.(y)))
    return tot
end

AugmentedKL(l::StudentTLikelihood,::AbstractVector)= InverseGammaKL(l)

function InverseGammaKL(l::StudentTLikelihood{T}) where {T}
    α_p = l.ν/2; β_p= α_p*l.σ^2
    InverseGammaKL(l.α,l.c,α_p,β_p)
end

## PDF and Log PDF Gradients ## (verified gradients)

function grad_log_pdf(l::StudentTLikelihood{T},y::Real,f::Real) where {T<:Real}
    (one(T)+l.ν) * (y-f) / ((f-y)^2 + l.σ^2*l.ν)
end

function hessian_log_pdf(l::StudentTLikelihood{T},y::Real,f::Real) where {T<:Real}
    v = l.ν * l.σ^2; Δ² = (f-y)^2
    (one(T)+l.ν) * (-v + Δ²) / (v+Δ²)^2
end


text:
# module Variant

export Variant, getChrom, getPos, getId, getRef, getAlt, getQual, getFilter, getInfo, getFormat

#Stores data for each variant
#type Variant
mutable struct  Variant
    CHROM::AbstractString
    POS::Int64
    ID::Array{AbstractString,1}
    REF::AbstractString
    ALT::Array{AbstractString,1}
    #QUAL::Array{Float64, 1}
    QUAL
    FILTER::Array{AbstractString,1}
    INFO::Dict{AbstractString, Any}
    FORMAT::Dict{AbstractString, Dict{AbstractString, Any}}
end

#Get methods for variants
getChrom(v::Variant) = v.CHROM
getPos(v::Variant) = v.POS
getId(v::Variant) = v.ID
getRef(v::Variant) = v.REF
getAlt(v::Variant) = v.ALT
getQual(v::Variant) = v.QUAL
getFilter(v::Variant) = v.FILTER
getInfo(v::Variant) = v.INFO
getFormat(v::Variant) = v.FORMAT

#end


text:
# this is used to test that relative include paths work on other processes

i27429 = 27429


text:
# Helper... Refactoring needed.
function getfieldinterpolation(dh::Ferrite.DofHandler, field_name::Symbol)
    field_idx = indexin(dh.field_names, [:b])
    field_idx == nothing && error("did not find field $field_name")
    dh.field_interpolations[field_idx]
end

function field_offset(dh::Ferrite.DofHandler, field_idx::Int)
    offset = 0
    for i in 1:(field_idx-1)
        offset += Ferrite.getnbasefunctions(dh.field_interpolations[i])::Int * dh.field_dims[i]
    end
    return offset
end

function facedofs(cell, local_face_idx::Int, field_idx::Int)
    celldofs_ = Ferrite.celldofs(cell)
    offset_ = field_offset(cell.dh, field_idx)
    fip = Ferrite.getfieldinterpolation(cell.dh, field_idx)
    faces_ = Ferrite.faces(fip)
    isempty(faces_) && return () # no face dofs :)
    indices = faces_[local_face_idx] .+ offset_
    celldofs_[[indices...]]
end

Ferrite.vertices(cell::Ferrite.Cell{3,3,1}) = cell.nodes

Ferrite.default_interpolation(::Type{Ferrite.Cell{3,3,1}}) = Ferrite.Lagrange{2,Ferrite.RefTetrahedron,1}()

# Note: This extracts the face spanned by the vertices, not the actual face!
linear_face_cell(cell::Ferrite.Cell{3,N,4}, local_face_idx::Int) where N =  Ferrite.Cell{3,3,1}(Ferrite.faces(cell)[local_face_idx])
linear_face_cell(cell::Ferrite.Cell{3,N,6}, local_face_idx::Int) where N = Ferrite.Quadrilateral3D(Ferrite.faces(cell)[local_face_idx])

# Obtain the face interpolation on regular geometries.
getfaceip(ip::Ferrite.Interpolation{dim, shape, order}, local_face_idx::Int) where {dim, shape <: Union{Ferrite.RefTetrahedron, Ferrite.RefCube}, order} = Ferrite.getlowerdim(ip)

struct MakiePlotter{dim,DH<:Ferrite.AbstractDofHandler,T} <: AbstractPlotter
    dh::DH
    u::Makie.Observable{Vector{T}} # original solution on the original mesh (i.e. dh.mesh)

    cells_connectivity::Vector

    gridnodes::Matrix{AbstractFloat} # coordinates of grid nodes in matrix form
    physical_coords::Matrix{AbstractFloat} # coordinates in physical space of a vertex
    triangles::Matrix{Int} # each row carries a triple with the indices into the coords matrix defining the triangle
    reference_coords::Matrix{AbstractFloat} # coordinates on the reference cell for the corresponding vertex
end

function MakiePlotter(dh::Ferrite.AbstractDofHandler, u::Vector)
    cells = Ferrite.getcells(dh.grid)
    dim = Ferrite.getdim(dh.grid)

    # We do not take any assumptions on the mesh, so first we have to loopkup
    num_triangles = 0
    for cell in cells
        num_triangles += ntriangles(cell)
    end

    # Preallocate the matrices carrying the triangulation information
    triangles = Matrix{Int}(undef, num_triangles, 3)
    physical_coords = Matrix{Float64}(undef, num_triangles*3, dim)
    gridnodes = [node[i] for node in Ferrite.getcoordinates.(Ferrite.getnodes(dh.grid)), i in 1:dim]
    reference_coords = Matrix{Float64}(undef, num_triangles*3, 2) # for now no 1D

    # Decompose does the heavy lifting for us
    coord_offset = 1
    triangle_offset = 1
    for cell in cells
        (coord_offset, triangle_offset) = decompose!(coord_offset, physical_coords, reference_coords, triangle_offset, triangles, dh.grid, cell)
    end
    return MakiePlotter{dim,typeof(dh),eltype(u)}(dh,Observable(u),[],gridnodes,physical_coords,triangles,reference_coords);
end

"""
Total number of vertices
"""
num_vertices(p::MakiePlotter) = size(p.physical_coords,1)

#Visualization is just fancy triangle plotting, every element needs to be translatable to a triangle *sadnoises*
to_triangle(cell::Ferrite.AbstractCell{2,N,3}) where N = [Ferrite.vertices(cell)[1:3]]
to_triangle(cell::Ferrite.AbstractCell{3,N,4}) where N = [Ferrite.vertices(cell)[[1,2,3]], Ferrite.vertices(cell)[[1,2,4]], Ferrite.vertices(cell)[[2,3,4]], Ferrite.vertices(cell)[[1,4,3]]]
to_triangle(cell::Ferrite.AbstractCell{2,N,4}) where N = [Ferrite.vertices(cell)[[1,2,3]], Ferrite.vertices(cell)[[3,4,1]]]
to_triangle(cell::Ferrite.AbstractCell{3,N,6}) where N = [Ferrite.vertices(cell)[[1,2,3]], Ferrite.vertices(cell)[[3,4,1]],
                                                          Ferrite.vertices(cell)[[1,5,6]], Ferrite.vertices(cell)[[6,2,1]],
                                                          Ferrite.vertices(cell)[[2,6,7]], Ferrite.vertices(cell)[[7,3,2]],
                                                          Ferrite.vertices(cell)[[3,7,8]], Ferrite.vertices(cell)[[8,4,3]],
                                                          Ferrite.vertices(cell)[[1,4,8]], Ferrite.vertices(cell)[[8,5,1]],
                                                          Ferrite.vertices(cell)[[5,8,7]], Ferrite.vertices(cell)[[7,6,5]]]

"""
TODO this looks faulty...think harder.
"""
# Helper to count triangles e.g. for preallocations.
ntriangles(cell::Ferrite.AbstractCell{2,3,3}) where {N} = 1 # Tris in 2D
ntriangles(cell::Ferrite.AbstractCell{3,3,1}) where {N} = 1 # Tris in 3D
ntriangles(cell::Ferrite.AbstractCell{dim,N,4}) where {dim,N} = 4 # Quads in 2D and 3D
ntriangles(cell::Ferrite.AbstractCell{3,N,1}) where {dim,N} = 4 # Tets as a special case of a Quad, obviously :)
ntriangles(cell::Ferrite.AbstractCell{3,N,6}) where {dim,N} = 6*4 # Hex

"""
Get the vertices represented as a list of coordinates of a cell.

@TODO refactor into Ferrite core.
"""
function vertices(grid::Ferrite.AbstractGrid, cell::Ferrite.AbstractCell{dim,N,M}) where {dim,N,M}
    Ferrite.getnodes(grid)[[Ferrite.vertices(cell)...]]
end

"""
Decompose a triangle into a coordinates and a triangle index list to disconnect it properly. Guarantees to preserve orderings and orientations.
"""
function decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, cell::Union{Ferrite.AbstractCell{2,N,3}, Ferrite.AbstractCell{3,3,1}}) where {N}
    for (i,v) in enumerate(vertices(grid, cell))
        coord_matrix[coord_offset, :] = Ferrite.getcoordinates(v)
        ref_coord_matrix[coord_offset, :] = Ferrite.reference_coordinates(Ferrite.default_interpolation(typeof(cell)))[i]
        triangle_matrix[triangle_offset, i] = coord_offset
        coord_offset+=1
    end
    triangle_offset+=1
    (coord_offset, triangle_offset)
end

"""
Decompose a tetrahedron into a coordinates and a triangle index list to disconnect it properly. Guarantees to preserve orderings and orientations.
"""
function decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, cell::Ferrite.AbstractCell{3,N,4}) where {N}
    # Is just 4 triangles :)
    for ti ∈ Ferrite.faces(cell)
        (coord_offset, triangle_offset) = decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, Ferrite.Cell{3,3,1}(ti))
    end
    (coord_offset, triangle_offset)
end

"""
Decompose a quadrilateral into a coordinates and a triangle index list to disconnect it properly. Guarantees to preserve orderings and orientations.

Assumes a CCW quadrilateral numbering:
4---3
|   |
1---2
and creates the decomposition
4-------3
| \\ C / |
|  \\ /  |
|D  5  B|
|  / \\  |
| / A \\ |
1-------2
where A=(1,2,5),B=(2,3,5),C=(3,4,5),D=(4,1,5) are the generated triangles in this order.
"""
function decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, cell::Union{Ferrite.AbstractCell{2,N,4}, Ferrite.AbstractCell{3,4,1}}) where {N}
    # A bit more complicated. The default diagonal decomposition into 2 triangles is missing a solution mode.
    # To resolve this we make a more expensive decomposition in 4 triangles which correctly displays the solution mode (in linear problems).
    coord_offset_initial = coord_offset
    vts = vertices(grid, cell)
    space_dim = size(coord_matrix,2)

    # Compute coordinate of vertex 5
    center = [ntuple(x->0.0, space_dim)...]
    for v in vts
        center += Ferrite.getcoordinates(v)
    end
    center /= 4.0

    # Generate triangles in order
    for i = 1:length(vts)
        v1 = vts[i]
        # next index on the outer chain CCW
        i2 = i+1
        if i2 > length(vts)
            i2 = 1 # dunno how to modulo this correctly :/
        end
        v2 = vts[i2]

        # current vertex
        coord_matrix[coord_offset, :] = Ferrite.getcoordinates(v1)
        ref_coord_matrix[coord_offset, :] = Ferrite.reference_coordinates(Ferrite.default_interpolation(typeof(cell)))[i]
        coord_offset+=1

        # next vertex in chain
        coord_matrix[coord_offset, :] = Ferrite.getcoordinates(v2)
        ref_coord_matrix[coord_offset, :] = Ferrite.reference_coordinates(Ferrite.default_interpolation(typeof(cell)))[i2]
        coord_offset+=1

        # center vertex (5)
        coord_matrix[coord_offset, :] = center
        ref_coord_matrix[coord_offset, :] = [ntuple(x->0.0, 2)...]
        coord_offset+=1

        # collect indices
        triangle_matrix[triangle_offset, :] = (0:2) .+ (coord_offset_initial+3*(i-1))
        triangle_offset+=1
    end

    (coord_offset, triangle_offset)
end

"""
Decompose a hexahedron into a coordinates and a triangle index list to disconnect it properly. Guarantees to preserve orderings and orientations.
"""
function decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, cell::Ferrite.AbstractCell{3,N,6}) where {N}
    # Just 6 quadrilaterals :)
    for ti ∈ Ferrite.faces(cell)
        (coord_offset, triangle_offset) = decompose!(coord_offset, coord_matrix, ref_coord_matrix, triangle_offset, triangle_matrix, grid, Ferrite.Cell{3,4,1}(ti))
    end
    (coord_offset, triangle_offset)
end

refshape(cell::Ferrite.AbstractCell) = typeof(Ferrite.default_interpolation(typeof(cell))).parameters[2]

x₁(x) = x[1]
x₂(x) = x[2]
x₃(x) = x[3]
l2(x) = LinearAlgebra.norm(x,2)
l1(x) = LinearAlgebra.norm(x,1)

midpoint(cell::Ferrite.AbstractCell{2,N,3}, points) where N = Point2f((1/3) * (points[cell.nodes[1],:] + points[cell.nodes[2],:] + points[cell.nodes[3],:]))
midpoint(cell::Ferrite.AbstractCell{2,N,4}, points) where N = Point2f(0.5 * (points[cell.nodes[1],:] + points[cell.nodes[3],:]))
midpoint(cell::Ferrite.AbstractCell{3,N,4}, points) where N = Point3f((1/4) * (points[cell.nodes[1],:] + points[cell.nodes[2],:] + points[cell.nodes[3],:] + points[cell.nodes[4],:]))
midpoint(cell::Ferrite.AbstractCell{3,N,6}, points) where N = Point3f(0.5 * (points[cell.nodes[1],:] + points[cell.nodes[7],:]))

function postprocess(node_values)
    dim = length(node_values)
    if dim == 1
        return node_values
    else
        return sqrt(sum(node_values.^2))
    end
end

"""
Transfer the solution of a plotter to the new mesh in 2D.

@precondition: Assumes that the number of basis function for each dof is equal.

@TODO Refactor. This is peak inefficiency.
"""
function transfer_solution(plotter::MakiePlotter{2}, u::Vector; field_idx::Int=1, process::Function=FerriteViz.postprocess)
    n_vertices_per_tri = 3 # we have 3 vertices per triangle...

    # select objects from plotter
    dh = plotter.dh
    ref_coords = plotter.reference_coords
    grid = dh.grid

    # field related variables
    field_name = Ferrite.getfieldnames(dh)[field_idx]
    field_dim = Ferrite.getfielddim(dh, field_idx)
    ip = dh.field_interpolations[field_idx]

    # actual data
    local_dof_range = Ferrite.dof_range(dh, field_name)

    data = fill(0.0, num_vertices(plotter), field_dim)
    current_vertex_index = 1
    for (cell_index, cell) in enumerate(Ferrite.getcells(plotter.dh.grid))
        cell_geo = Ferrite.getcells(dh.grid,cell_index)
        _celldofs_field = reshape(Ferrite.celldofs(dh,cell_index)[local_dof_range], (field_dim, Ferrite.getnbasefunctions(ip)))

        # Loop over all local triangle vertices
        for i in 1:(ntriangles(cell_geo)*n_vertices_per_tri)
            ξ = Tensors.Vec(ref_coords[current_vertex_index, :]...)
            for d in 1:field_dim
                for node_idx ∈ 1:Ferrite.getnbasefunctions(ip)
                    data[current_vertex_index, d] += Ferrite.value(ip, node_idx, ξ) ⋅ u[_celldofs_field[d, node_idx]]
                end
            end
            current_vertex_index += 1
        end
    end
    return mapslices(process, data, dims=[2])
end

"""
Transfer the solution of a plotter to the new mesh in 3D. We just evaluate the faces.

@TODO Refactor. This is peak inefficiency.
"""
function transfer_solution(plotter::MakiePlotter{3}, u::Vector; field_idx::Int=1, process::Function=FerriteViz.postprocess) where T
    n_vertices_per_tri = 3 # we have 3 vertices per triangle...

    # select objects from plotter
    dh = plotter.dh
    ref_coords = plotter.reference_coords
    grid = dh.grid

    # field related variables
    field_name = Ferrite.getfieldnames(dh)[field_idx]
    field_dim = Ferrite.getfielddim(dh, field_idx)

    # @FIXME decouple the interpolation from the geometry, because for discontinuous interpolations
    #   the "interpolation faces" (which we need for dof-assignment) do not coincide with the geometric
    #   faces... Alternatively we could try something similar to FaceValues.
    ip_cell = dh.field_interpolations[field_idx]
    _faces = Ferrite.faces(ip_cell) # faces of the cell with local dofs
    @assert !isempty(_faces) # Discontinuous interpolations in 3d not supported yet. See above.

    # actual data
    local_dof_range = Ferrite.dof_range(dh, field_name)

    current_vertex_index = 1
    data = fill(0.0, num_vertices(plotter), field_dim)
    for (cell_index, cell) in enumerate(Ferrite.getcells(plotter.dh.grid))
        cell_geo = Ferrite.getcells(dh.grid,cell_index)

        _local_celldofs = Ferrite.celldofs(dh,cell_index)[local_dof_range]
        _celldofs_field = reshape(_local_celldofs, (field_dim, Ferrite.getnbasefunctions(ip_cell)))

        for (local_face_idx,_) in enumerate(Ferrite.faces(cell_geo))
            ip_face = getfaceip(ip_cell, local_face_idx)
            face_geo = linear_face_cell(cell_geo, local_face_idx)

            # extract face dofs
            _facedofs_field = _celldofs_field[:,[_faces[local_face_idx]...]]

            # Loop over vertices
            for i in 1:(ntriangles(face_geo)*n_vertices_per_tri)
                ξ = Tensors.Vec(ref_coords[current_vertex_index, :]...)
                for d in 1:field_dim
                    for node_idx ∈ 1:Ferrite.getnbasefunctions(ip_face)
                        data[current_vertex_index, d] += Ferrite.value(ip_face, node_idx, ξ) ⋅ u[_facedofs_field[d, node_idx]]
                    end
                end
                current_vertex_index += 1
            end
        end
    end

    return mapslices(process, data, dims=[2])
end

function transfer_scalar_celldata(plotter::MakiePlotter{3}, u::Vector; process::Function=FerriteViz.postprocess)
    n_vertices = 3 # we have 3 vertices per triangle...

    # select objects from plotter
    dh = plotter.dh
    grid = dh.grid

    current_vertex_index = 1
    data = fill(0.0, num_vertices(plotter), 1)
    for (cell_index, cell) in enumerate(Ferrite.getcells(grid))
        cell_geo = grid.cells[cell_index]
        for (local_face_idx,_) in enumerate(Ferrite.faces(cell_geo))
            face_geo = linear_face_cell(cell_geo, local_face_idx)
            # Loop over vertices
            for i in 1:(ntriangles(face_geo)*n_vertices)
                data[current_vertex_index, 1] = u[cell_index]
                current_vertex_index += 1
            end
        end
    end

    return mapslices(process, data, dims=[2])
end

function transfer_scalar_celldata(plotter::MakiePlotter{2}, u::Vector;  process::Function=FerriteViz.postprocess)
    n_vertices = 3 # we have 3 vertices per triangle...

    # select objects from plotter
    dh = plotter.dh
    grid = dh.grid

    current_vertex_index = 1
    data = fill(0.0, num_vertices(plotter), 1)
    for (cell_index, cell) in enumerate(Ferrite.getcells(grid))
        cell_geo = grid.cells[cell_index]
        for i in 1:(ntriangles(cell_geo)*n_vertices)
            data[current_vertex_index, 1] = u[cell_index]
            current_vertex_index += 1
        end
    end

    return mapslices(process, data, dims=[2])
end

function transfer_scalar_celldata(grid::Ferrite.AbstractGrid{3}, num_vertices::Number, u::Vector; process::Function=FerriteViz.postprocess)
    n_vertices = 3 # we have 3 vertices per triangle...
    current_vertex_index = 1
    data = fill(0.0, num_vertices, 1)
    for (cell_index, cell) in enumerate(Ferrite.getcells(grid))
        cell_geo = grid.cells[cell_index]
        for (local_face_idx,_) in enumerate(Ferrite.faces(cell_geo))
            face_geo = linear_face_cell(cell_geo, local_face_idx)
            # Loop over vertices
            for i in 1:(ntriangles(face_geo)*n_vertices)
                data[current_vertex_index, 1] = u[cell_index]
                current_vertex_index += 1
            end
        end
    end
    return mapslices(process, data, dims=[2])
end

function transfer_scalar_celldata(grid::Ferrite.AbstractGrid{2}, num_vertices::Number, u::Vector;  process::Function=FerriteViz.postprocess)
    n_vertices = 3 # we have 3 vertices per triangle...
    current_vertex_index = 1
    data = fill(0.0, num_vertices, 1)
    for (cell_index, cell) in enumerate(Ferrite.getcells(grid))
        cell_geo = grid.cells[cell_index]
        for i in 1:(ntriangles(cell_geo)*n_vertices)
            data[current_vertex_index, 1] = u[cell_index]
            current_vertex_index += 1
        end
    end
    return mapslices(process, data, dims=[2])
end

function dof_to_node(dh::Ferrite.AbstractDofHandler, u::Array{T,1}; field::Int=1, process::Function=postprocess) where T
    fieldnames = Ferrite.getfieldnames(dh)
    field_dim = Ferrite.getfielddim(dh, field)
    data = fill(NaN, Ferrite.getnnodes(dh.grid), field_dim)
    offset = Ferrite.field_offset(dh, fieldnames[field])

    for cell in Ferrite.CellIterator(dh)
        _celldofs = Ferrite.celldofs(cell)
        counter = 1
        for node in cell.nodes
            for d in 1:field_dim
                data[node, d] = u[_celldofs[counter + offset]]
                counter += 1
            end
        end
    end
    return mapslices(process, data, dims=[2])
end


text:
export ObservationGrup

"""
Observation Groups help us iterate through a panel

Examples:
- A vector of production from a well
- The full set of leases before or after 1st well drilled
- The set of actions associated w/ 1 lease
"""
struct ObservationGroup{D<:AbstractDataStructure,I} <: AbstractObservationGroup
    data::D
    i::I
    function ObservationGroup(data::D, i::I) where {D<:AbstractDataStructure,I}
        i in eachindex(data) || throw(BoundsError(data,i))
        return new{D,I}(data,i)
    end
end

const ObservationGroupEmpty = ObservationGroup{EmptyDataSet}

# Functions for these data structures
#-----------------------

# ObservationGroup
length(o::AbstractObservation) = length(_y(o))

# default methods for iteration through an AbstractDataStructure
firstindex(d::AbstractDataStructure) = 1
lastindex( d::AbstractDataStructure) = length(d)
IndexStyle(d::AbstractDataStructure) = IndexLinear()
eachindex( d::AbstractDataStructure) = OneTo(length(d))

# Default iteration method is to create an ObservationGroup
getindex(  d::AbstractDataStructure, i) = ObservationGroup(d,i)

function iterate(d::AbstractDataStructure, i=firstindex(d))
    if i <= lastindex(d)
        return getindex(d,i), i+1
    else
        return nothing
    end
end

# AbstractDataSet iteration utilties
#------------------------------------------

_data(d::AbstractDataSet) = d

group_ptr(d::AbstractDataSet) = throw(error("group_ptr not defined for $(typeof(d))"))
obs_ptr(  d::AbstractDataSet) = throw(error("obs_ptr not defined for $(typeof(d))"))

# default method
length(   d::AbstractDataSet) = length(group_ptr(d))-1
_num_obs( d::AbstractDataSet) = length(obs_ptr(d))-1

groupstart( d::AbstractDataSet, i::Integer) = getindex(group_ptr(d), i)
groupstop(  d::AbstractDataSet, i::Integer) = groupstart(d,i+1)-1
grouplength(d::AbstractDataSet, i::Integer) = groupstop(d,i) - groupstart(d,i) + 1
grouprange( d::AbstractDataSet, i::Integer) = groupstart(d,i) : groupstop(d,i)

obsstart( d::AbstractDataSet, j::Integer) = getindex(obs_ptr(d),j)
obsstop(  d::AbstractDataSet, j::Integer) = obsstart(d,j+1)-1
obslength(d::AbstractDataSet, j::Integer) = obsstop(d,j) - obsstart(d,j) + 1
obsrange( d::AbstractDataSet, j::Integer) = obsstart(d,j) : obsstop(d,j)

# DataSet or Observation
_model(d::DataOrObs) = d.model
_y(    d::DataOrObs) = d.y
_x(    d::DataOrObs) = d.x
_num_x(d::DataOrObs) = size(_x(d), 1)

# ObservationGroup iteration utilties
#------------------------------------------

_data( g::AbstractObservationGroup) = g.data
_i(    g::AbstractObservationGroup) = g.i
_num_x(g::AbstractObservationGroup) = _num_x(_data(g))
_nparm(g::AbstractObservationGroup) = _nparm(_data(g))
_model(g::AbstractObservationGroup) = _model(_data(g))

length(    g::AbstractObservationGroup) = grouplength(_data(g), _i(g))
grouprange(g::AbstractObservationGroup) = grouprange( _data(g), _i(g))

obsstart( g::AbstractObservationGroup, k) = obsstart( _data(g), getindex(grouprange(g), k))
obsrange( g::AbstractObservationGroup, k) = obsrange( _data(g), getindex(grouprange(g), k))
obslength(g::AbstractObservationGroup, k) = obslength(_data(g), getindex(grouprange(g), k))

# Observation iteration utilties
#------------------------------------------

function update_over_obs(f!::Function, data::AbstractDataSet)
    for j in OneTo(_num_obs(data))
        f!( Observation(data, j) )
    end
    return nothing
end

# default method
update!(d::AbstractDataSet, theta) = nothing


text:
import Base: det, logdet, diagm, Diagonal, diag
export diag, diagm, Diagonal

push!(ops, DiffOp(:(Base.diag), :(Tuple{DLA.AM}), [true]))
function ∇(::typeof(diag), ::Arg1, p, y::AV, ȳ::AV, X::AM)
    X̄ = fill!(similar(X), zero(eltype(X)))
    X̄[diagind(X̄)] = ȳ
    return X̄
end
function ∇(X̄::AM, ::typeof(diag), ::Arg1, p, y::AV, ȳ::AV, X::AM)
    X̄_diag = view(X̄, diagind(X̄))
    X̄_diag .+= ȳ
    return X̄
end

push!(ops, DiffOp(:(Base.diag), :(Tuple{DLA.AM, Integer}), [true, false]))
function ∇(::typeof(diag), ::Arg1, p, y::AV, ȳ::AV, X::AM, k::Integer)
    X̄ = fill!(similar(X), zero(eltype(X)))
    X̄[diagind(X̄, k)] = ȳ
    return X̄
end
function ∇(X̄::AM, ::typeof(diag), ::Arg1, p, y::AV, ȳ::AV, X::AM, k::Integer)
    X̄_diag = view(X̄, diagind(X̄, k))
    X̄_diag .+= ȳ
    return X̄
end

push!(ops, DiffOp(:(Base.diagm), :(Tuple{DLA.AV}), [true]))
∇(::typeof(diagm), ::Arg1, p, Y::AM, Ȳ::AM, x::AV) =
    copy!(similar(x), view(Ȳ, diagind(Ȳ)))
∇(x̄::AV, ::typeof(diagm), ::Arg1, p, Y::AM, Ȳ::AM, x::AV) =
    broadcast!(+, x̄, x̄, view(Ȳ, diagind(Ȳ)))

push!(ops, DiffOp(:(Base.diagm), :(Tuple{DLA.AV, Integer}), [true, false]))
∇(::typeof(diagm), ::Arg1, p, Y::AM, Ȳ::AM, x::AV, k::Integer) =
    copy!(similar(x), view(Ȳ, diagind(Ȳ, k)))
∇(x̄::AV, ::typeof(diagm), ::Arg1, p, Y::AM, Ȳ::AM, x::AV, k::Integer) =
    broadcast!(+, x̄, x̄, view(Ȳ, diagind(Ȳ, k)))

push!(ops, DiffOp(:(Base.diagm), :(Tuple{Real}), [true]))
function ∇(::typeof(diagm), ::Arg1, p, Y::AM, Ȳ::AM, x::Real)
    length(Ȳ) != 1 && throw(error("Ȳ isn't a 1x1 matrix."))
    return Ȳ[1]
end

push!(ops, DiffOp(:(Base.Diagonal), :(Tuple{DLA.AV}), [true]))
∇(::Type{Diagonal}, ::Arg1, p, Y::Diagonal{<:Real}, Ȳ::Diagonal{<:Real}, x::AV) =
    copy!(similar(x), Ȳ.diag)
∇(x̄::AV, ::Type{Diagonal}, ::Arg1, p, Y::Diagonal{<:Real}, Ȳ::Diagonal{<:Real}, x::AV) =
    broadcast!(+, x̄, x̄, Ȳ.diag)

push!(ops, DiffOp(:(Base.Diagonal), :(Tuple{DLA.AM}), [true]))
function ∇(::Type{Diagonal}, ::Arg1, p, Y::Diagonal{<:Real}, Ȳ::Diagonal{<:Real}, X::AM)
    X̄ = zeros(X)
    copy!(view(X̄, diagind(X)), Ȳ.diag)
    return X̄
end
function ∇(X̄::AM, ::Type{Diagonal}, ::Arg1, p, Y::Diagonal{<:Real}, Ȳ::Diagonal{<:Real}, X::AM)
    X̄_diag = view(X̄, diagind(X̄))
    broadcast!(+, X̄_diag, X̄_diag, Ȳ.diag)
    return X̄
end

push!(ops, DiffOp(:(Base.det), :(Tuple{Diagonal{<:Real}}), [true]))
∇(::typeof(det), ::Arg1, p, y::Real, ȳ::Real, X::Diagonal{<:Real}) =
    Diagonal(ȳ .* y ./ X.diag)
function ∇(X̄::Diagonal{<:Real}, ::typeof(det), ::Arg1, p, y::Real, ȳ::Real, X::Diagonal{<:Real})
    broadcast!((x̄, x, y, ȳ)->x̄ + ȳ * y / x, X̄.diag, X̄.diag, X.diag, y, ȳ)
    return X̄
end

push!(ops, DiffOp(:(Base.logdet), :(Tuple{Diagonal{<:Real}}), [true]))
∇(::typeof(logdet), ::Arg1, p, y::Real, ȳ::Real, X::Diagonal{<:Real}) = Diagonal(ȳ ./ X.diag)
function ∇(X̄::Diagonal{<:Real}, ::typeof(logdet), ::Arg1, p, y::Real, ȳ::Real, X::Diagonal{<:Real})
    broadcast!((x̄, x, ȳ)->x̄ + ȳ / x, X̄.diag, X̄.diag, X.diag, ȳ)
    return X̄
end


text:
using Smg2s, Test, SparseArrays, LinearAlgebra

@testset "Initialization of Matrix" begin
    @testset "Init 1" begin
        Am=spzeros(ComplexF64, 50, 50)
        initMat!(Am, -20, -10, 50, sparsity=0.1)
        @test nnz(Am) ≈ 38 atol=5
        @test UpperTriangular(Am)-diagm(diag(Am)) == zeros(50,50)
    end
    @testset "Init 2" begin
        Am = spzeros(Float64, 50, 50);
        initMat!(Am, -20, -10, 50, sparsity=0.05)
        @test nnz(Am) ≈ 19 atol=5
        @test UpperTriangular(Am)-diagm(diag(Am)) == zeros(50,50)
        @test diag(Am,-1) == zeros(49)
    end
end


text:
function explicit_euler_step(qᵢ, vᵢ, v̇ᵢ, h)
    qᵢ₊₁ = qᵢ + h * vᵢ
    vᵢ₊₁ = vᵢ + h * v̇ᵢ

    [qᵢ₊₁
     vᵢ₊₁]
end

function semi_implicit_euler_step(qᵢ, vᵢ, v̇ᵢ, h)
    vᵢ₊₁ = vᵢ + h * v̇ᵢ
    qᵢ₊₁ = qᵢ + h * vᵢ₊₁

    [qᵢ₊₁
     vᵢ₊₁]
end

function forward_dynamics_defects!(defects, robot, x::AbstractArray{T}, h) where {T}
    # Indices of the decision variables
    ind_qᵢ   = (1:robot.n_q)
    ind_vᵢ   = (1:robot.n_v) .+ (robot.n_q)
    ind_τᵢ   = (1:robot.n_τ) .+ (robot.n_q + robot.n_v)
    ind_qᵢ₊₁ = (1:robot.n_q) .+ (robot.n_q + robot.n_v + robot.n_τ)
    ind_vᵢ₊₁ = (1:robot.n_v) .+ (robot.n_q + robot.n_v + robot.n_τ + robot.n_q)

    # Retrieve values of the decision variables
    qᵢ  , vᵢ  , τᵢ = x[ind_qᵢ  ], x[ind_vᵢ  ], x[ind_τᵢ]
    qᵢ₊₁, vᵢ₊₁     = x[ind_qᵢ₊₁], x[ind_vᵢ₊₁]

    # Update robot state
    state = robot.statecache[T]
    set_configuration!(state, qᵢ)
    set_velocity!(state, vᵢ)

    # Compute forward dynamics
    result = robot.dynamicsresultcache[T]
    dynamics!(result, state, τᵢ)
    v̇ᵢ = result.v̇

    # Integrate the robot state
    xᵢ₊₁_integrated = semi_implicit_euler_step(qᵢ, vᵢ, v̇ᵢ, h)

    # Next robot state discretised
    xᵢ₊₁_discretised = [qᵢ₊₁
                        vᵢ₊₁]

    # Evaluate state defects (joint positions and joint velocities)
    copyto!(defects, xᵢ₊₁_integrated - xᵢ₊₁_discretised)
end

function inverse_dynamics_defects!(defects, robot, x::AbstractArray{T}, h) where {T}
    # Indices of the decision variables
    ind_qᵢ   = (1:robot.n_q)
    ind_vᵢ   = (1:robot.n_v) .+ (robot.n_q)
    ind_τᵢ   = (1:robot.n_τ) .+ (robot.n_q + robot.n_v)
    ind_qᵢ₊₁ = (1:robot.n_q) .+ (robot.n_q + robot.n_v + robot.n_τ)
    ind_vᵢ₊₁ = (1:robot.n_v) .+ (robot.n_q + robot.n_v + robot.n_τ + robot.n_q)

    # Retrieve values of the decision variables
    qᵢ  , vᵢ  , τᵢ = x[ind_qᵢ  ], x[ind_vᵢ  ], x[ind_τᵢ]
    qᵢ₊₁, vᵢ₊₁     = x[ind_qᵢ₊₁], x[ind_vᵢ₊₁]

    # Update robot state
    state = robot.statecache[T]
    set_configuration!(state, qᵢ)
    set_velocity!(state, vᵢ)

    # Implicit joint accelerations
    v̇ᵢ = similar(velocity(state))
    copyto!(v̇ᵢ, (vᵢ₊₁ - vᵢ) / h)

    # Compute inverse dynamics
    torques = similar(velocity(state))  # create a SegmentedVector instance
    dynamics_result = robot.dynamicsresultcache[T]
    wrenches = dynamics_result.jointwrenches
    accelerations = dynamics_result.accelerations
    inverse_dynamics!(torques, wrenches, accelerations, state, v̇ᵢ)

    # position_defects = (qᵢ + h * vᵢ) - qᵢ₊₁  # Explicit Euler
    position_defects = (qᵢ + h * vᵢ₊₁) - qᵢ₊₁  # Semi-Implicit Euler
    torque_defects = torques - τᵢ

    # Evaluate defects (joint positions and joint torques)
    copyto!(defects, [position_defects
                        torque_defects])
end

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

function make_cb_eval_fc_con_dyn(problem, robot, jacdata, eval_fc!::Function)
    nₓ = robot.n_q + robot.n_v + robot.n_τ  # dimension of each mesh points
    length_c = size(jacdata.jac, 1)  # length of constraint per pair of knots

    function cb_eval_fc_con_dyn(kc, cb, evalRequest, evalResult, userParams)
        x = evalRequest.x

        for i = 0:problem.num_knots - 2
            # Calculate the indices of the appropriate decision variables
            ind_vars = range(1 + i * nₓ, length=size(jacdata.jac, 2))

            # Evaluate constraints for Knitro
            offset_con = i * length_c
            ind_con = (1:length_c) .+ offset_con
            @views eval_fc!(evalResult.c[ind_con], robot, x[ind_vars], problem.dt)
        end

        return 0
    end
end

function make_cb_eval_ga_con_dyn(problem, robot, jacdata)
    nₓ = robot.n_q + robot.n_v + robot.n_τ  # dimension of each mesh points

    function cb_eval_ga_con_dyn(kc, cb, evalRequest, evalResult, userParams)
        x = evalRequest.x

        for i = 0:problem.num_knots - 2
            # Calculate the indices of the appropriate decision variables
            ind_vars = range(1 + i * nₓ, length=size(jacdata.jac, 2))

            # Evaluate the Jacobian at that point
            jacdata(x[ind_vars])

            # Pass the Jacobian to Knitro
            offset_jac = i * jacdata.length_jac
            ind_jac = (1:jacdata.length_jac) .+ offset_jac
            evalResult.jac[ind_jac] = nonzeros(jacdata.jac)
        end

        return 0
    end
end


text:
# Function declarations for documention

"""
    rownumber(e::Entry)
    rownumber(r::Row)

Returns the row of Pascal's triangle that the entry `e` lies within or that the row `r`
represents

# Examples
```jldoctest
julia> rownumber(Entry(6,3))
6

julia> rownumber(Row(6))
6
```
"""
function rownumber end

"""
    rowposition(e::Entry)

Returns the position of the entry `e` within the row of Pascal's triangle.

# Examples
```jldoctest
julia> rowposition(Entry(6,3))
3
```
"""
function rowposition end

"""
    value(e::Entry)
    value(r::Row) or values(r::Row)
    value(c::Column) or values(c::Column)
    value(c::Centre) or values(c::Centre)

Returns the value of the entry `e`, or the values of the row `r`, column `c`
or a sequence of central elements `c`.

# Examples
```jldoctest
julia> value(Entry(6,3))
20

julia> value(Row(4))
5-element Array{Integer,1}:
 1
 4
 6
 4
 1

julia> value(Column(4,3))
3-element Array{Integer,1}:
  1
  5
 15
```
"""
function value end

"""
    isfirst(e::Entry)
    isfirst(r::Row)
    isfirst(c::Column)

Returns whether the given entry, row or column is at the top-left of Pascal's triangle.
That is, if `e` is the top entry, `r` is the top row or `c` is the left-most column.
Note that rows and columns are numbered from zero.

# Examples
```jldoctest
julia> isfirst(Entry(0,0))
true

julia> isfirst(Row(0))
true

julia> isfirst(Column(0,0))
true
```
"""
function isfirst end

"""
    isatleft(e::Entry)
    isatleft(c::Column)

Returns true if the entry or column are at the left edge of Pascal's triangle. That is,
if entry `e` is in the left-most column, or `c` is column number zero. This corresponds
to entries (n,0).

# Examples
```jldoctest
julia> isatleft(Entry(4,0))
true

julia> isatleft(Column(0,3))
true
```
"""
function isatleft end

"""
    isatright(e::Entry)

Returns true if the entry is at the right edge of Pascal's triangle. That is, if entry
`e` corresponds to an entry (n,n).
"""
function isatright end

"""
    isadjacent(a::Entry, b::Entry) or areadjacent(a::Entry, b::Entry)

Returns true if the two entries are directly adjacent on the same row of Pascal's
triangle.

# Examples
```jldoctest
julia> isadjacent(Entry(6,3), Entry(6,4))
true

julia> areadjacent(Entry(6,2), Entry(6,6))
false
```
"""
function isadjacent end
const areadjacent = isadjacent

"""
    isinterior(a::Entry)

Returns true if the entry is in the interior of Pascal's triangle (that is, the
value isn't 1)
"""
function isinterior end

"""
    issubtractable(a::Entry, b::Entry) or aresubtractable(a::Entry, b::Entry)

Returns true if `a` and `b` are arranged so that they can be subtracted. This is
equivalent to `a` being in the row beneath `b` and either directly beneath or
one place to the right.

# Examples
```jldoctest
julia> issubtractable(Entry(6,3), Entry(5,3))
true

julia> aresubtractable(Entry(6,3), Entry(5,2))
true
```
"""
function issubtractable end
const aresubtractable = issubtractable

"""
    isvalid(e::Entry)
    isvalid(r::Row)
    isvalid(c::Column)
    isvalid(c::Centre)
    isvalid(l::LazyCentre)

Returns true if the data contained in `e`, `r`, `c` or `l` represents a
valid entry, row, column or central elements. This is useful for checking
if the internal state has been corrupted but can be slow and is not normally
necessary when using the methods of this package.

See [`Row`](@ref) for the internal structure of a `Row`.
"""
function isvalid end

"""
    toarray(r::Row)
    toarray(c::Column)
    toarray(c::Centre)
    toarray(l::LazyCentre)

Returns an array of `Entry`s equivalent to the given row, column or central
elements. For a `LazyCentre` this will only include the elements that have
actually been calculated (which may be more than those that have been
explicitly requested).

```jldoctest
julia> toarray(Row(2))
3-element Array{Entry{Int64,Int64},1}:
 Entry{Int64,Int64}(2, 0, 1)
 Entry{Int64,Int64}(2, 1, 2)
 Entry{Int64,Int64}(2, 2, 1)

julia> toarray(Column(2,3))
3-element Array{Entry{Int64,Int64},1}:
 Entry{Int64,Int64}(2, 2, 1)
 Entry{Int64,Int64}(3, 2, 3)
 Entry{Int64,Int64}(4, 2, 6)
```
"""
function toarray end

"""
    up(e::Entry)
    up(r::Row)
    up!(r::Row)

Returns the entry or row above the given one in Pascal's triangle. Throws
an [`OutOfBoundsError`](@ref) if there is no entry or row above.

# Examples
```jldoctest
julia> up(Entry(4,3))
Entry{Int64,Int64}(3, 3, 1)

julia> up(Row(2))
2-element Row{Integer} with indices PascalsTriangle.ZeroRange(1):
 1
 1

julia> up(Entry(4,4))
ERROR: OutOfBoundsError: no entry above
[...]
```
"""
function up end
function up! end

"""
    down(e::Entry)
    down(r::Row)
    down!(r::Row)

Returns the entry or row directly below the given one in Pascal's triangle.

# Examples
```jldoctest
julia> down(Entry(3,2))
Entry{Int64,Int64}(4, 2, 6)

julia> down(Row(2))
4-element Row{Integer} with indices PascalsTriangle.ZeroRange(3):
 1
 3
 3
 1
```
"""
function down end
function down! end

"""
    left(e::Entry)
    left(c::Column)
    left!(c::Column)

Returns the entry or column to the left of the given one in Pascal's triangle.
Throws an [`OutOfBoundsError`](@ref) if the entry or column is at the left edge
(ie column number 0, or an entry (n,0)).

# Examples
```jldoctest
julia> left(Entry(4,2))
Entry{Int64,Int64}(4, 1, 4)

julia> left(Column(2,3))
3-element Column{Integer}:
 1
 2
 3

julia> left!(Column(0,0))
ERROR: OutOfBoundsError: no previous column
[...]
```
"""
function left end
function left! end

"""
    right(e::Entry)
    right(c::Column)
    right!(c::Column)

Returns the entry or column to the right of the given one in Pascal's triangle.
Throws an [`OutOfBoundsError`](@ref) if the entry is at the right edge (ie an
entry (n,n))

# Examples
```jldoctest
julia> right(Entry(4,2))
Entry{Int64,Int64}(4, 3, 4)

julia> right(Column(4,0))
0-element Column{Integer}

julia> right(Entry(4,4))
ERROR: OutOfBoundsError: no entry to the right
[...]
```
"""
function right end
function right! end

"""
    prev(e::Entry)
    prev(r::Row)
    prev(c::Column)
    prev!(r::Row)
    prev!(c::Column)

Returns the previous entry, row or column. The previous row is the one
above the given `r` and the previous column is the one to the left of
the given `c`. The previous entry is to the left of the given `e` unless
its at the start of the row, in which case it will wrap to the right of
the previous row.

Throws an [`OutOfBoundsError`](@ref) if there is no previous entry, row
or column.

# Examples
```jldoctest
julia> prev(Entry(4,0))
Entry{Int64,Int64}(3, 3, 1)

julia> prev(Row(3))
3-element Row{Integer} with indices PascalsTriangle.ZeroRange(2):
 1
 2
 1

julia> prev(Column(0,0))
ERROR: OutOfBoundsError: no previous column
[...]
```
"""
function prev end
function prev! end

"""
    next(e::Entry)
    next(r::Row)
    next(c::Column)
    next!(r::Row)
    next!(c::Column)

Returns the next entry, row or column in Pascal's triangle. The next row
is the one beneath the given `r` and the next column is the one to the
right of the given `c`. The next entry is to the right unless `e` is at the
end of the row, in which case it will wrap to the beginning of the next
row.

# Examples
```jldoctest
julia> next(Entry(4,4))
Entry{Int64,Int64}(5, 0, 1)

julia> next(Row(3))
5-element Row{Integer} with indices PascalsTriangle.ZeroRange(4):
 1
 4
 9
 4
 1

julia> next(Column(4,0))
0-element Column{Integer}
```
"""
function next end
function next! end