text:
# REDUCTION ALGORITHMS

abstract type AbstractPersistenceReduction end
struct StandardReduction <: AbstractPersistenceReduction end
struct TwistReduction <: AbstractPersistenceReduction end

# Boundary matrix reduction algorithms
lastindex(col::AbstractSet) = length(col) == 0 ? -1 : maximum(col)

"""Standart reduction"""
function reduce!(::Type{StandardReduction}, ∂::Vector{<:AbstractSet})
    lowest_one_lookup = fill(-1, length(∂))
    for col in eachindex(∂)
        lowest_one = lastindex(∂[col])
        while lowest_one != -1 && lowest_one_lookup[lowest_one] != -1
            lowest = lowest_one_lookup[lowest_one]
            symdiff!(∂[col], ∂[lowest])
            lowest_one = lastindex(∂[col])
        end
        if lowest_one != -1
            lowest_one_lookup[lowest_one] = col
        end
    end
    return ∂
end

"""Twist reduction"""
function reduce!(::Type{TwistReduction}, ∂::Vector{<:AbstractSet})
    lowest_one_lookup = fill(-1, length(∂))

    for dim in maximum(map(length, ∂)):-1:1
        for col in eachindex(∂)
            if length(∂[col]) == dim
                lowest_one = lastindex(∂[col])
                while lowest_one != -1 && lowest_one_lookup[lowest_one] != -1
                    lowest = lowest_one_lookup[lowest_one]
                    symdiff!(∂[col], ∂[lowest])
                    lowest_one = lastindex(∂[col])
                end
                if lowest_one != -1
                    lowest_one_lookup[lowest_one] = col
                    empty!(∂[lowest_one])
                end
            end
        end
    end
    return ∂
end

reduce(::Type{R}, ∂::Vector{<:AbstractSet}) where {R<:AbstractPersistenceReduction} = reduce!(R, deepcopy(∂))

"""Return persistent diagram (birth-death pairs) per dimension."""
function diagram(flt::Filtration{C,FI}, R::Vector; length0=false,
                 absolute=true, maxoutdim=dim(complex(flt))) where {C <: AbstractComplex, FI <: AbstractFloat}
    # resulting intervals
    intrs = Dict{Int,Vector{Interval{FI}}}()

    # compute intervals
    births = Set{Int}()
    ord = order(flt)
    for (i, (sdim, si, fv)) in enumerate(ord)
        sdim > maxoutdim && continue
        col = R[i]
        if length(col) == 0
            push!(births, i)
        else
            b = lastindex(col)
            d = i
            delete!(births, b)
            delete!(births, d)
            if d > b
                sdim = absolute ? ord[b][1] : ord[d][1]
                bv = ord[b][3]
                dv = ord[d][3]
                if dv > bv || length0
                    !haskey(intrs, sdim) && setindex!(intrs, Interval[], sdim)
                    push!(intrs[sdim], Interval(bv, dv))
                end
            end
        end
    end
    for i in births
        sdim, si, fv = ord[i]
        !haskey(intrs, sdim) && setindex!(intrs, Interval[], sdim)
        push!(intrs[sdim], Interval(fv, FI(Inf)))
    end

    return intrs
end
diagram(::Type{R}, flt::Filtration; kwargs...) where {R <: AbstractPersistenceReduction} =
    diagram(flt, reduce!(R, boundary(flt)); kwargs...)
diagram(flt::Filtration; kwargs...) = diagram(TwistReduction, flt; kwargs...)

"Calculate persistent Betti numbers for a filtration complex of dimension `dim`"
function betti(flt::Filtration, R::Vector, p::Int)
    maxdim =  dim(complex(flt))
    @assert maxdim >= p "Cannot calculate $p-dimensional Betti number for $maxdim-complex"

    # indexes of p-simplices
    p2i = findall(d->d[1]==p, order(flt))

    # the number of zero columns that correspond to p-simplices
    z = sum(map(length, R[p2i]) .== 0)

    # the number of lowest ones in rows that correspond to p-simplices
    l = length( findall(l->l ∈ p2i, map(lastindex, R)) )

    β = z - l
    return β < 0 ? 0 : β
end
betti(::Type{R}, flt::Filtration, p::Int) where {R <: AbstractPersistenceReduction} =
    betti(flt, reduce!(R, boundary(flt)), p)
betti(flt::Filtration, p::Int) = betti(TwistReduction, flt, p)

#
# Persistent Cohomology
#

function coboundarymap(coch::Chain{IX,R}, ∂::Chain{IX,R}) where {IX<:Integer, R}
    sum(coch[be]*bc for (be,bc) in ∂)
end

function cocycles(ccycs, ∂::Chain{IX,R}, w::Float64) where {IX<:Integer, R}
    cᵢ = zeros(R, length(ccycs))
    for (i,intr) in enumerate(ccycs)
        cᵢ[i] = if w ∉ intr
            zero(R)
        else
            coboundarymap(intr.g, ∂)
        end
    end
    cᵢ
end

function update!(cyc::Dict{Int,Vector{AnnotatedInterval{F,Chain{IX,R}}}},
                 cplx, splxs, w::F) where {F<:AbstractFloat, IX<:Integer, R}
    d = first(first(splxs))
    d > maximum(keys(cyc)) && return
    CXₗ₋₁ = cyc[d-1]
    for splx in splxs
        σ = cplx[last(splx), d]
        dαᵢ = cocycles(CXₗ₋₁, boundary(R, σ), w)
        idxs = findall(!iszero, dαᵢ)
        if !isempty(idxs)
            j = idxs[end]
            cⱼ = dαᵢ[j]
            αⱼ = CXₗ₋₁[j].g
            for i in idxs[1:end-1]
                cᵢ = dαᵢ[i]
                CXₗ₋₁[i].g -= (cᵢ/cⱼ)*αⱼ
            end
            CXₗ₋₁[j].d = w
        else
            intr = AnnotatedInterval(w, F(Inf), Chain(d, [splx[2]], [one(R)]))
            push!(cyc[d], intr)
        end
    end
end

function persistentcohomology(::Type{R}, flt::Filtration;
                              length0=false, maxoutdim=dim(complex(flt))) where {R}
    cplx = complex(flt)
    d = dim(cplx)
    itr = Iterators.flatten(
            ((v => filter(e->first(e) == i, splxs))
                for i in (:)(extrema(map(first, splxs))...)) for (v, splxs) in flt )

    cyc = Dict(i=>AnnotatedInterval{R,Chain{UInt64,R}}[] for i in -1:min(d, maxoutdim))
    for (w, splxs) in itr
        @debug "Filtration step" Value=w Simpices=length(splxs)
        update!(cyc, cplx, splxs, w)
    end
    delete!(cyc, -1)

    # remove all intervals of 0-length
    cyc = !length0 ? Dict( d=>filter!(!isempty, c) for (d,c) in cyc) : cyc

    return cyc
end
persistentcohomology(flt::Filtration; kwargs...) = persistentcohomology(Float64, flt; kwargs...)

"""
Persistent Cocycle Reduction

Vin de Silva, Dmitriy Morozov, Mikael Vejdemo-Johansson, Persistent Cohomology and Circular Coordinates,
Discrete Comput Geom (2011) 45: 737–759, DOI 10.1007/s00454-011-9344-x
"""
struct PersistentCocycleReduction{R} <: AbstractPersistenceReduction end

diagram(::Type{PersistentCocycleReduction{R}}, flt::Filtration; kwargs...) where {R} =
    persistentcohomology(R, flt; kwargs...)


text:
using Documenter, MLDatasets

makedocs(
    modules = [MLDatasets],
    clean = false,
    format = :html,
    assets = [
        joinpath("assets", "favicon.ico"),
    ],
    sitename = "MLDatasets.jl",
    authors = "Hiroyuki Shindo, Christof Stocker",
    linkcheck = !("skiplinks" in ARGS),
    pages = Any[
        "Home" => "index.md",
        "Available Datasets" => Any[
            "Image Classification" => Any[
                "MNIST handwritten digits" => "datasets/MNIST.md",
                "Fashion MNIST" => "datasets/FashionMNIST.md",
                "CIFAR-10" => "datasets/CIFAR10.md",
                "CIFAR-100" => "datasets/CIFAR100.md",
                "SVHN format 2" => "datasets/SVHN2.md",
            ],
        ],
        hide("Indices" => "indices.md"),
        "LICENSE.md",
    ],
    html_prettyurls = !("local" in ARGS),
)

deploydocs(
    repo = "github.com/JuliaML/MLDatasets.jl.git",
    target = "build",
    julia = "0.7",
    deps = nothing,
    make = nothing,
)

ext: jl


text:
import Nemo.crt, Nemo.zero, Nemo.iszero, Nemo.isone, Nemo.sub!
export crt_env, crt, crt_inv, modular_init, crt_signed

isone(a::Int) = (a==1)

function zero(a::PolyElem)
  return zero(parent(a))
end

@inline function rem!(a::fmpz, b::fmpz, c::fmpz)
  ccall((:fmpz_mod, :libflint), Nothing, (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), a, b, c)
end

@inline function sub!(a::fmpz, b::fmpz, c::fmpz)
  ccall((:fmpz_sub, :libflint), Nothing, (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), a, b, c)
end

function rem!(a::fmpz_mod_poly, b::fmpz_mod_poly, c::fmpz_mod_poly)
  ccall((:fmpz_mod_poly_rem, :libflint), Nothing, (Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}, Ref{fmpz_mod_poly}), a, b, c)
end

mutable struct crt_env{T}
  pr::Array{T, 1}
  id::Array{T, 1}
  tmp::Array{T, 1}
  t1::T
  t2::T
  n::Int
  M::T #for T=fmpz, holds prod/2

  function crt_env{T}(p::Array{T, 1}) where {T}
    pr = deepcopy(p)
    id = Array{T, 1}()
    i = 1
    while 2*i <= length(pr)
      a = pr[2*i-1]
      b = pr[2*i]
      if false
        g, u, v = gcdx(a, b)
        @assert isone(g)
        push!(id, v*b , u*a )
      else
        # we have 1 = g = u*a + v*b, so v*b = 1-u*a - which saves one mult.
        u = invmod(a, b)
#        @assert  (a*u) % b == 1
        u *= a
        push!(id, 1-u, u)
      end
      push!(pr, a*b)
      i += 1
    end
    r = new{T}()
    r.pr = pr
    r.id = id

    r.tmp = Array{T, 1}()
    n = length(p)
    for i=1:div(n+1, 2)
      push!(r.tmp, zero(p[1]))
    end
    r.t1 = zero(p[1])
    r.t2 = zero(p[1])

    r.n = n
    return r
  end
end

@doc Markdown.doc"""
    crt_env(p::Array{T, 1}) -> crt_env{T}

Given coprime moduli in some euclidean ring (FlintZZ, nmod_poly,
fmpz\_mod_poly), prepare data for fast application of the chinese
remander theorem for those moduli.
"""
function crt_env(p::Array{T, 1}) where T
  return crt_env{T}(p)
end

function show(io::IO, c::crt_env{T}) where T
  print(io, "CRT data for moduli ", c.pr[1:c.n])
end

@doc Markdown.doc"""
    crt{T}(b::Array{T, 1}, a::crt_env{T}) -> T

Given values in b and the environment prepared by crt\_env, return the
unique (modulo the product) solution to $x \equiv b_i \bmod p_i$.
"""
function crt(b::Array{T, 1}, a::crt_env{T}) where T
  res = zero(b[1])
  return crt!(res, b, a)
end

function crt!(res::T, b::Array{T, 1}, a::crt_env{T}) where T
  @assert a.n == length(b)
  bn = div(a.n, 2)
  if isodd(a.n)
    @inbounds zero!(a.tmp[1])
    @inbounds add!(a.tmp[1], a.tmp[1], b[end])
    off = 1
  else
    off = 0
  end

  for i=1:bn
    @inbounds mul!(a.t1, b[2*i-1], a.id[2*i-1])
    @inbounds mul!(a.t2, b[2*i], a.id[2*i])
    @inbounds add!(a.tmp[i+off], a.t1, a.t2)
    @inbounds rem!(a.tmp[i+off], a.tmp[i+off], a.pr[a.n+i])
  end

  if isodd(a.n)
    bn += 1
  end

  id_off = a.n - off
  pr_off = a.n + bn - off
#  println(a.tmp, " id_off=$id_off, pr_off=$pr_off, off=$off, bn=$bn")
  while bn>1
    if isodd(bn)
      off = 1
    else
      off = 0
    end
    bn = div(bn, 2)
    for i=1:bn
      if true  # that means we need only one co-factor!!!
        @inbounds sub!(a.t1, a.tmp[2*i-1], a.tmp[2*i])
        @inbounds mul!(a.t1, a.t1, a.id[id_off + 2*i-1])
        @inbounds add!(a.tmp[i+off], a.t1, a.tmp[2*i])
      else
        @inbounds mul!(a.t1, a.tmp[2*i-1], a.id[id_off + 2*i-1])
        @inbounds mul!(a.t2, a.tmp[2*i], a.id[id_off + 2*i])
        @inbounds add!(a.tmp[i + off], a.t1, a.t2)
      end
      @inbounds rem!(a.tmp[i + off], a.tmp[i + off], a.pr[pr_off+i])
    end
    if off == 1
      @inbounds a.tmp[1], a.tmp[2*bn+1] = a.tmp[2*bn+1], a.tmp[1]
    end
    id_off += 2*bn
    pr_off += bn
    bn += off
#    println(a.tmp, " id_off=$id_off, pr_off=$pr_off, off=$off, bn=$bn")
  end
  zero!(res)
  @inbounds add!(res, res, a.tmp[1])
  return res
end

function crt_signed!(res::fmpz, b::Array{fmpz, 1}, a::crt_env{fmpz})
  crt!(res, b, a)
  if !isdefined(a, :M)
    a.M = div(prod(a.pr[1:a.n]), 2)
  end
  if res>a.M
    sub!(res, res, a.pr[end])
  end
end

function crt_signed(b::Array{fmpz, 1}, a::crt_env{fmpz})
  res = fmpz()
  crt_signed!(res, b, a)
  return res
end

#in .pr we have the products of pairs, ... in the wrong order
# so we traverse this list backwards, while building the remainders...
#.. and then we do it again, efficiently to avoid resorting and re-allocation
#=
function crt_inv{T}(a::T, c::crt_env{T})
  r = Array{T, 1}()
  push!(r, a)
  i = length(c.pr)-1
  j = 1
  while i>1
    push!(r, r[j] % c.pr[i], r[j] %c.pr[i-1])
    i -= 2
    j += 1
  end
  return reverse(r, length(r)-c.n+1:length(r))
end
=#

function crt_inv_iterative!(res::Array{T,1}, a::T, c::crt_env{T}) where T
  for i=1:c.n
    if isassigned(res, i)
      rem!(res[i], a, c.pr[i])
    else
      res[i] = a % c.pr[i]
    end
  end
  return res
end

function crt_inv_tree!(res::Array{T,1}, a::T, c::crt_env{T}) where T
  for i=1:c.n
    if !isassigned(res, i)
      res[i] = zero(a)
    end
  end

  i = length(c.pr)-1
  if i == 0
    rem!(res[1], a, c.pr[1])
    return res
  end
  r = i
  w = r + c.n - 1

  @inbounds zero!(res[r % c.n + 1])
  @inbounds add!(res[r % c.n + 1], res[r % c.n + 1], a)

  while i>1
    @inbounds rem!(res[w % c.n + 1], res[r % c.n + 1], c.pr[i])
    @inbounds rem!(res[(w+c.n - 1) % c.n + 1], res[r % c.n + 1], c.pr[i - 1])
    w += 2*(c.n-1)
    i -= 2
    r += 1*(c.n-1)
  end

  return res
end

@doc Markdown.doc"""
    crt_inv(a::T, crt_env{T}) -> Array{T, 1}

Given a $\code{crt_env}$ and an element a, return
the modular data $a \bmod pr_i$ for all $i$.
This is essentially the inverse to the $\code{crt}$ function.
"""
function crt_inv(a::T, c::crt_env{T}) where T
  res = Array{T}(undef, c.n)
  if c.n < 50
    return crt_inv_iterative!(res, a, c)
  else
    return crt_inv_tree!(res, a, c)
  end
end

function crt_inv!(res::Array{T, 1}, a::T, c::crt_env{T}) where T
  if c.n < 50
    return crt_inv_iterative!(res, a, c)
  else
    return crt_inv_tree!(res, a, c)
  end
end

#explains the idea, but is prone to overflow.
# idea: the tree CRT ..
# given moduli p1 .. pn, we first do (p1, p2), (p2, p3), ...
# then ((p1, p2), (p3, p4)), ...
# until done.
# In every step we need the cofactors, the inverse of pi mod pj
# thus we build a parallel array id for the cofactors
# in id[2i-1], id[2i] are the cofactors for pr[2i-1], pr[2i]
# To recombine, we basically loop through the cofactors:
# use id[1], id[2] to combine b[1], b[2] AND append to b
# The product pr[1]*pr[2] was appended to pr, thus we can walk through the
# growing list till the end
# For the optimized version, we have tmp-array to hold the CRT results
# plus t1, t2 for temporaty products.
function crt(b::Array{Int, 1}, a::crt_env{Int})
  i = a.n+1
  j = 1
  while j <= length(b)
    push!(b, (b[j-1]*a.id[j-1] + b[j]*a.id[j]) % a.pr[i])
    j += 2
    i += 1
  end
  return b[end]
end

function crt_test(a::crt_env{fmpz}, b::Int)
  z = [fmpz(0) for x=1:a.n]
  for i=1:b
    b = rand(0:a.pr[end]-1)
    for j=1:a.n
      rem!(z[j], b, a.pr[j])
    end
    if b != crt(z, a)
      println("problem: $b and $z")
    end
    @assert b == crt(z, a)
  end
end

@doc Markdown.doc"""
  crt(r1::PolyElem, m1::PolyElem, r2::PolyElem, m2::PolyElem) -> PolyElem

Find $r$ such that $r \equiv r_1 \pmod m_1$ and $r \equiv r_2 \pmod m_2$
"""
function crt(r1::PolyElem{T}, m1::PolyElem{T}, r2::PolyElem{T}, m2::PolyElem{T}) where T
  g, u, v = gcdx(m1, m2)
  m = m1*m2
  return (r1*v*m2 + r2*u*m1) % m
end

@doc Markdown.doc"""
  crt_iterative(r::Array{T, 1}, m::Array{T,1}) -> T

Find $r$ such that $r \equiv r_i \pmod m_i$ for all $i$.
A plain iteration is performed.
"""
function crt_iterative(r::Array{T, 1}, m::Array{T, 1}) where T
  p = crt(r[1], m[1], r[2], m[2])
  d = m[1] * m[2]
  for i = 3:length(m)
    p = crt(p, d, r[i], m[i])
    d *= m[i]
  end
  return p
end

@doc Markdown.doc"""
  crt_tree(r::Array{T, 1}, m::Array{T,1}) -> T

Find $r$ such that $r \equiv r_i \pmod m_i$ for all $i$.
A tree based strategy is used that is asymptotically fast.
"""
function crt_tree(r::Array{T, 1}, m::Array{T, 1}) where T
  if isodd(length(m))
    M = [m[end]]
    V = [r[end]]
  else
    M = Array{T, 1}()
    V = Array{T, 1}()
  end

  for i=1:div(length(m), 2)
    push!(V, crt(r[2*i-1], m[2*i-1], r[2*i], m[2*i]))
    push!(M, m[2*i-1]*m[2*i])
  end
  i = 1
  while 2*i <= length(V)
    push!(V, crt(V[2*i-1], M[2*i-1], V[2*i], M[2*i]))
    push!(M, M[2*i-1] * M[2*i])
    i += 1
  end
#  println("M = $M\nV = $V")
  return V[end]
end

@doc Markdown.doc"""
  crt(r::Array{T, 1}, m::Array{T,1}) -> T

Find $r$ such that $r \equiv r_i \pmod m_i$ for all $i$.
"""
function crt(r::Array{T, 1}, m::Array{T, 1}) where T
  length(r) == length(m) || error("Arrays need to be of same size")
  if length(r) < 15
    return crt_iterative(r, m)
  else
    return crt_tree(r, m)
  end
end

function crt_test_time_all(np::Int, n::Int)
  p = next_prime(fmpz(2)^60)
  m = [p]
  x = fmpz(1)
  for i=1:np-1
    push!(m, next_prime(m[end]))
  end
  v = [rand(0:x-1) for x = m]
  println("crt_env...")
  @time ce = crt_env(m)
  @time for i=1:n
    x = crt(v, ce)
  end

  println("iterative...")
  @time for i=1:n
    x = crt_iterative(v, m)
  end

  println("tree...")
  @time for i=1:n
    x = crt_tree(v, m)
  end

  println("inv_tree")
  @time for i=1:n
    crt_inv_tree!(m, x, ce)
  end

  println("inv_iterative")
  @time for i=1:n
    crt_inv_iterative!(m, x, ce)
  end
end

@doc Markdown.doc"""
    induce_crt(a::fmpz_poly, p::fmpz, b::fmpz_poly, q::fmpz, signed::Bool = false) -> fmpz_poly
Given integral polynomials $a$ and $b$ as well as coprime integer moduli
$p$ and $q$, find $f = a \bmod p$ and $f=b \bmod q$.
If signed is set, the symmetric representative is used, the positive one
otherwise.
"""
function induce_crt(a::fmpz_poly, p::fmpz, b::fmpz_poly, q::fmpz, signed::Bool = false)
  c = parent(a)()
  pi = invmod(p, q)
  mul!(pi, pi, p)
  pq = p*q
  if signed
    pq2 = div(pq, 2)
  else
    pq2 = fmpz(0)
  end
  for i=0:max(degree(a), degree(b))
    setcoeff!(c, i, Hecke.inner_crt(coeff(a, i), coeff(b, i), pi, pq, pq2))
  end
  return c, pq
end

@doc Markdown.doc"""
    induce_crt(L::Array{PolyElem, 1}, c::crt_env{fmpz}) -> fmpz_poly
Given fmpz\_poly polynomials $L[i]$ and a {{{crt\_env}}}, apply the
{{{crt}}} function to each coefficient resulting in a polynomial $f = L[i] \bmod p[i]$.
"""
function induce_crt(L::Array{T, 1}, c::crt_env{fmpz}) where {T <: PolyElem}
  Zx, x = FlintZZ["x"]
  res = Zx()
  m = maximum(degree(x) for x = L)

  for i=0:m
    setcoeff!(res, i, crt([lift(coeff(x, i)) for x =L], c))
  end
  return res
end

@doc Markdown.doc"""
  _num_setcoeff!(a::nf_elem, n::Int, c::fmpz)
  _num_setcoeff!(a::nf_elem, n::Int, c::Integer)

Sets the $n$-th coefficient in $a$ to $c$. No checks performed, use
only if you know what you're doing.
"""
function _num_setcoeff!(a::nf_elem, n::Int, c::fmpz)
  K = parent(a)
  @assert n < degree(K) && n >=0
  ra = pointer_from_objref(a)
  if degree(K) == 1
    ccall((:fmpz_set, :libflint), Nothing, (Ref{Nothing}, Ref{fmpz}), ra, c)
    ccall((:fmpq_canonicalise, :libflint), Nothing, (Ref{nf_elem}, ), a)
  elseif degree(K) == 2
     ccall((:fmpz_set, :libflint), Nothing, (Ref{Nothing}, Ref{fmpz}), ra+n*sizeof(Int), c)
  else
    ccall((:fmpq_poly_set_coeff_fmpz, :libflint), Nothing, (Ref{nf_elem}, Int, Ref{fmpz}), a, n, c)
   # includes canonicalisation and treatment of den.
  end
end

function _num_setcoeff!(a::nf_elem, n::Int, c::UInt)
  K = a.parent
  @assert n < degree(K) && n >=0

  ra = pointer_from_objref(a)

  if degree(K) == 1
    ccall((:fmpz_set_ui, :libflint), Nothing, (Ref{Nothing}, UInt), ra, c)
    ccall((:fmpq_canonicalise, :libflint), Nothing, (Ref{nf_elem}, ), a)
  elseif degree(K) == 2
    ccall((:fmpz_set_ui, :libflint), Nothing, (Ref{Nothing}, UInt), ra+n*sizeof(Int), c)
  else
    ccall((:fmpq_poly_set_coeff_ui, :libflint), Nothing, (Ref{nf_elem}, Int, UInt), a, n, c)
   # includes canonicalisation and treatment of den.
  end
end

function _num_setcoeff!(a::nf_elem, n::Int, c::Integer)
  _num_setcoeff!(a, n, fmpz(c))
end

@doc Markdown.doc"""
    induce_crt(L::Array{MatElem, 1}, c::crt_env{fmpz}) -> fmpz_mat
Given matrices $L[i]$ and a {{{crt\_env}}}, apply the
{{{crt}}} function to each coefficient resulting in a matrix $M = L[i] \bmod p[i]$.
"""
function induce_crt(L::Array{T, 1}, c::crt_env{fmpz}, signed::Bool = false) where {T <: MatElem}
  res = zero_matrix(FlintZZ, nrows(L[1]), ncols(L[1]))

  if signed
    cr = crt_signed
  else
    cr = crt
  end

  for i=1:nrows(L[1])
    for j=1:ncols(L[1])
      res[i,j] = cr([lift(x[i,j]) for x =L], c)
    end
  end
  return res
end

mutable struct modular_env
  p::fmpz
  up::UInt
  upinv::UInt

  fld::Array{FqNmodFiniteField, 1}
  fldx::Array{FqNmodPolyRing, 1}
  ce::crt_env{nmod_poly}
  rp::Array{nmod_poly, 1}
  res::Array{fq_nmod, 1}
  Fpx::NmodPolyRing
  K::AnticNumberField
  Rp::Array{fq_nmod_poly, 1}
  Kx::Generic.PolyRing{nf_elem}
  Kxy::Generic.MPolyRing{nf_elem}
  Kpxy::NmodMPolyRing

  function modular_env()
    return new()
  end
end

Base.isempty(me::modular_env) = !isdefined(me, :ce)

function show(io::IO, me::modular_env)
  if isempty(me)
    println("modular environment for p=$(me.p), using $(0) ideals")
  else
    println("modular environment for p=$(me.p), using $(me.ce.n) ideals")
  end
end

@doc Markdown.doc"""
  modular_init(K::AnticNumberField, p::fmpz) -> modular_env
  modular_init(K::AnticNumberField, p::Integer) -> modular_env

Given a number field $K$ and an ``easy'' prime $p$ (ie. fits into an
\code{Int} and is coprime to the polynomial discriminant), compute
the residue class fields of the associated primes ideals above $p$.
Returns data that can be used by \code{modular_proj} and \code{modular_lift}.
"""
function modular_init(K::AnticNumberField, p::fmpz; deg_limit::Int=0, max_split::Int = 0)
  @assert isprime(p)
  UInt(p) # to enforce p being small

  me = modular_env()
  me.Fpx = PolynomialRing(ResidueRing(FlintZZ, Int(p), cached = false), "_x", cached=false)[1]
  fp = me.Fpx(K.pol)
  lp = factor(fp)
  if Set(values(lp.fac)) != Set([1])
    throw(BadPrime(p))
  end
  pols = collect(keys(lp.fac))
  if deg_limit > 0
    pols = pols[findall(x -> degree(x) <= deg_limit, pols)]
  end

  if max_split > 0
    pols = pols[1:min(length(pols), max_split)]
  end
  if length(pols) == 0
    return me
  end
  me.ce = crt_env(pols)
  me.fld = [FqNmodFiniteField(x, :$, false) for x = pols]  #think about F_p!!!
                                   # and chacheing
  me.rp = Array{nmod_poly}(undef, length(pols))
  me.res = Array{fq_nmod}(undef, me.ce.n)

  me.p = p
  me.K = K
  me.up = UInt(p)
  me.upinv = ccall((:n_preinvert_limb, :libflint), UInt, (UInt, ), me.up)
  return me
end

function modular_init(K::AnticNumberField, p::Integer; deg_limit::Int=0, max_split::Int = 0)
  return modular_init(K, fmpz(p), deg_limit = deg_limit, max_split = max_split)
end

@doc Markdown.doc"""
  modular_proj(a::nf_elem, me::modular_env) -> Array{fq_nmod, 1}

Given an algebraic number $a$ and data \code{me} as computed by
\code{modular_init}, project $a$ onto the residue class fields.
"""
function modular_proj(a::nf_elem, me::modular_env)
  ap = me.Fpx(a)
  crt_inv!(me.rp, ap, me.ce)
  for i=1:me.ce.n
    F = me.fld[i]
    if isassigned(me.res, i)
      u = me.res[i]
    else
      u = F()
    end
    ccall((:fq_nmod_set, :libflint), Nothing,
                (Ref{fq_nmod}, Ref{nmod_poly}, Ref{FqNmodFiniteField}),
                u, me.rp[i], F)
    me.res[i] = u
  end
  return me.res
end

@doc Markdown.doc"""
  modular_proj(a::FacElem{nf_elem, AnticNumberField}, me::modular_env) -> Array{fq_nmod, 1}

Given an algebraic number $a$ in factored form and data \code{me} as computed by
\code{modular_init}, project $a$ onto the residue class fields.
"""
function modular_proj(A::FacElem{nf_elem, AnticNumberField}, me::modular_env)
  for i=1:me.ce.n
    me.res[i] = one(me.fld[i])
  end
  for (a, v) = A.fac
    ap = me.Fpx(a)
    crt_inv!(me.rp, ap, me.ce)
    for i=1:me.ce.n
      F = me.fld[i]
      u = F()
      ccall((:fq_nmod_set, :libflint), Nothing,
                  (Ref{fq_nmod}, Ref{nmod_poly}, Ref{FqNmodFiniteField}),
                  u, me.rp[i], F)
      u = u^(v % (size(F)-1))
      me.res[i] *= u
    end
  end
  return me.res
end

@doc Markdown.doc"""
  modular_lift(a::Array{fq_nmod}, me::modular_env) -> nf_elem

Given an array of elements as computed by \code{modular_proj},
compute a global pre-image using some efficient CRT.
"""
function modular_lift(a::Array{fq_nmod, 1}, me::modular_env)
  for i=1:me.ce.n
    ccall((:nmod_poly_set, :libflint), Nothing, (Ref{nmod_poly}, Ref{fq_nmod}), me.rp[i], a[i])
  end
  ap = crt(me.rp, me.ce)
  r = me.K()
  for i=0:ap.length-1
    u = ccall((:nmod_poly_get_coeff_ui, :libflint), UInt, (Ref{nmod_poly}, Int), ap, i)
    _num_setcoeff!(r, i, u)
  end
  return r
end

@doc Markdown.doc"""
  modular_proj(a::Generic.Poly{nf_elem}, me::modular_env) -> Array

Apply the \code{modular_proj} function to each coeficient of $a$.
Computes an array of polynomials over the respective residue class fields.
"""
function modular_proj(a::Generic.Poly{nf_elem}, me::modular_env)

  if !isdefined(me, :fldx)
    me.fldx = [PolynomialRing(x, "_x", cached=false)[1] for x = me.fld]
    me.Rp = Array{fq_nmod_poly}(undef, me.ce.n)
    for i =1:me.ce.n
      me.Rp[i] = me.fldx[i](0)
    end
    me.Kx = parent(a)
  end

  for j=1:me.ce.n
    zero!(me.Rp[j])
  end
  for i=0:length(a)-1
    c = coeff(a, i)
    if iszero(mod(denominator(c), me.p))
      throw(BadPrime(me.p))
    end
    crt_inv!(me.rp, me.Fpx(c), me.ce)
    for j=1:me.ce.n
      u = coeff(me.Rp[j], i)
      ccall((:fq_nmod_set, :libflint), Nothing,
                   (Ref{fq_nmod}, Ref{nmod_poly}, Ref{FqNmodFiniteField}),
                   u, me.rp[j], me.fld[j])
      setcoeff!(me.Rp[j], i, u)
    end
  end
  return me.Rp
end

@doc Markdown.doc"""
  modular_lift(a::Array{fq_nmod_poly}, me::modular_env) -> Generic.Poly{nf_elem}

Apply the \code{modular_lift} function to each coeficient of $a$.
Computes a polynomial over the number field.
"""
function modular_lift(a::Array{fq_nmod_poly, 1}, me::modular_env)
  res = me.Kx()
  d = maximum([x.length for x = a])
  for i=0:d-1
    for j=1:me.ce.n
      ccall((:nmod_poly_set, :libflint), Nothing, (Ref{nmod_poly}, Ref{fq_nmod}), me.rp[j], coeff(a[j], i))
    end
    ap = crt(me.rp, me.ce)
    r = coeff(res, i)
    for j=0:ap.length-1
      u = ccall((:nmod_poly_get_coeff_ui, :libflint), UInt, (Ref{nmod_poly}, Int), ap, j)
      _num_setcoeff!(r, j, u)
    end
    setcoeff!(res, i, r)
  end
  return res
end

@doc Markdown.doc"""
  modular_proj(a::Generic.Mat{nf_elem}, me::modular_env) -> Array{Matrix}
  modular_proj(a::Generic.Mat{NfOrdElem}, me::modular_env) -> Array{Matrix}

Apply the \code{modular_proj} function to each entry of $a$.
Computes an array of matrices over the respective residue class fields.
"""
function modular_proj(a::Generic.Mat{nf_elem}, me::modular_env)
  Mp = fq_nmod_mat[]
  for i=1:me.ce.n
    push!(Mp, zero_matrix(me.fld[i], nrows(a), ncols(a)))
  end
  for i=1:nrows(a)
    for j=1:ncols(a)
      im =modular_proj(a[i,j], me)
      for k=1:me.ce.n
        setindex!(Mp[k], deepcopy(im[k]), i, j)
      end
    end
  end
  return Mp
end

function modular_proj(a::Generic.Mat{NfOrdElem}, me::modular_env)
  Mp = []
  for i=1:me.ce.n
    push!(Mp, zero_matrix(me.fld[i], nrows(a), ncols(a)))
  end
  for i=1:nrows(a)
    for j=1:ncols(a)
      im =modular_proj(me.K(a[i,j]), me)
      for k=1:me.ce.n
        setindex!(Mp[k], deepcopy(im[k]), i, j)
      end
    end
  end
  return Mp
end


text:
module ExponentialFamilies

import Distributions
using VariationalModelIndices

VERSION < v"0.4.0-dev" && using Docile

export sparse_mat_from_tuples
export MatrixTuple
export make_ud_index_matrix, linearize_matrix
export unpack_ud_matrix, unpack_ud_trace_coefficients
export get_mvn_parameters_from_derivs, get_mvn_variational_covariance
export wishart_entropy, wishart_e_log_det, get_wishart_variational_covariance
export get_wishart_sufficient_stats_variational_covariance
export get_wishart_parameters_from_derivs

# Multivariate log gamma related functions.
function multivariate_trigamma{T <: Number}(x::T, p::Int64)
	sum([ trigamma(x + 0.5 * (1 - i)) for i=1:p])
end

function multivariate_digamma{T <: Number}(x::T, p::Int64)
	sum([ digamma(x + 0.5 * (1 - i)) for i=1:p])
end

function multivariate_lgamma{T <: Number}(x::T, p::Int64)
	sum([ lgamma(x + 0.5 * (1 - i)) for i=1:p]) +
	p * (p - 1.0) / 4.0 * log(pi)
end

@doc """
Turn a vector of upper diagonal enries back into a symmetric matrix.
""" ->
function unpack_ud_matrix{T <: Number}(ud_vector::Vector{T}; od_scale=1.0)
	k_tot = linearized_matrix_size(length(ud_vector))
	ud_mat = Array(T, (k_tot, k_tot))
	for k1=1:k_tot, k2=1:k_tot
		ud_mat[k1, k2] =
			(k1 <= k2 ? ud_vector[(k1 + (k2 - 1) * k2 / 2)] :
				        ud_vector[(k2 + (k1 - 1) * k1 / 2)])
		ud_mat[k1, k2] *= k1 != k2 ? od_scale: 1.0
	end
	ud_mat
end

@doc """
Convert a vector of upper diagonal entries into a
matrix with halved off-diagonal entries.
This is what's needed to convert the coefficients
of the derivative wrt mu2 into a matrix V such that
tr(V * mu2) = coeffs' * mu2 """ ->
function unpack_ud_trace_coefficients(ud_vector)
	unpack_ud_matrix(ud_vector, od_scale=0.5)
end

# A tuple representing a matrix element
# [row, column, value]
typealias MatrixTuple Tuple{Int64, Int64, Float64}

function sparse_mat_from_tuples(tup_array::Array{MatrixTuple})
	sparse(Int64[x[1] for x=tup_array],
		   Int64[x[2] for x=tup_array],
		   Float64[x[3] for x=tup_array])
end

###################################################
# Normal functions

@doc """
beta ~ MVN(beta_mean, beta_cov)
...with elements beta_k

Returns: Cov(beta_k1 beta_k2, beta_k3 beta_k4)
""" ->
function get_mvn_fourth_order_cov(beta_mean::Array{Float64, 1},
	                                beta_cov::Array{Float64, 2},
	                                k1::Int64, k2::Int64, k3::Int64, k4::Int64)
  cov_12 = beta_cov[k1, k2];
  cov_13 = beta_cov[k1, k3];
  cov_14 = beta_cov[k1, k4];
  cov_23 = beta_cov[k2, k3];
  cov_24 = beta_cov[k2, k4];
  cov_34 = beta_cov[k3, k4];

	m_1 = beta_mean[k1];
	m_2 = beta_mean[k2];
	m_3 = beta_mean[k3];
	m_4 = beta_mean[k4];

  return (cov_13 * cov_24 +
  	      cov_14 * cov_23 +
  	      cov_13 * m_2 * m_4 +
  	      cov_14 * m_2 * m_3 +
	        cov_23 * m_1 * m_4 +
	        cov_24 * m_1 * m_3);
end

@doc """
Get the covariance matrix of a multivariate multivariate normal
sufficient statistics given the mean parameters.

E(beta) = beta_mean
Cov(beta) = beta_cov
beta_ind_model = A vector of the model indices of the E(beta)
beta2_ind_model = A matrix containing the model indices of E(beta beta')

Returns:
	An array of MatrixTuple.
""" ->
function get_mvn_variational_covariance(
	beta_mean::Array{Float64, 1}, beta_cov::Array{Float64, 2},
	beta_ind_model::Array{Int64}, beta2_ind_model::Matrix{Int64})

	k_tot = length(beta_ind_model)

	@assert k_tot == length(beta_mean) ==
					size(beta_cov, 1) == size(beta_cov, 2) ==
					size(beta2_ind_model, 1) == size(beta2_ind_model, 2)

	beta_cov_tuples = MatrixTuple[]

	# Get the linear covariances.
	for k1=1:k_tot, k2=1:k_tot
		i1 = beta_ind_model[k1]
		i2 = beta_ind_model[k2]
		push!(beta_cov_tuples, (i1, i2, beta_cov[k1, k2]))
	end

	# Get the covariance between the linear and quadratic terms.
	# This will be cov(mu_k1 mu_k2, mu_i3) := cov(mu2_i12, mu_i3).
	# Avoid double counting since only one mu_i1 mu_i2 is recorded.
	for k1=1:k_tot, k2=1:k1, k3=1:k_tot
		i12 = beta2_ind_model[k1, k2]
		i3 = beta_ind_model[k3]
		this_cov = (beta_mean[k1] * beta_cov[k2, k3] +
		            beta_mean[k2] * beta_cov[k1, k3])
		push!(beta_cov_tuples, (i3, i12, this_cov))
		push!(beta_cov_tuples, (i12, i3, this_cov))
	end

	# Get the covariance between the quadratic terms.
	# This will be cov(mu_k1 mu_k2, mu_k3 mu_k4) := cov(mu2_i12, mu_i34).
	# Avoid double counting since only one mu_k1 mu_k2
	# and mu_k3 mu_k4 is recorded.
	for k1=1:k_tot, k2=1:k1, k3=1:k_tot, k4=1:k3
		i12 = beta2_ind_model[k1, k2]
		i34 = beta2_ind_model[k3, k4]
		this_cov = get_mvn_fourth_order_cov(beta_mean, beta_cov,
			k1, k2, k3, k4)

		push!(beta_cov_tuples, (i12, i34, this_cov))
	end

	beta_cov_tuples
end

@doc """
Given a log likelihood's derivatives with respect to a MVN's mean
parameters, return its mean and covariance.  The log likelihood is

beta_deriv' beta + beta2_deriv' beta2

TODO: document beta2 a little more clearly.

Args:
	- beta_deriv: The derivative of the log likelihood with respect to beta
	- beta2_deriv: The derivative of the log likelihood with respect to
		           the upper diagonal of beta beta'.

Returns:
	- E(beta), Cov(beta)
""" ->
function get_mvn_parameters_from_derivs(
		beta_deriv::Array{Float64}, beta2_deriv::Array{Float64})
	k_tot = length(beta_deriv)
	@assert length(beta2_deriv) == k_tot * (k_tot + 1) / 2
	beta_dist =
		Distributions.MvNormalCanon(
			beta_deriv, -2 * unpack_ud_trace_coefficients(beta2_deriv))
	mean(beta_dist), cov(beta_dist)
end

@doc """
Get the normal covariance for a scalar normal with expectation
beta_mean, expecation of the square beta2_mean, and in columns
beta_ind_model and beta2_ind_model respectively.
""" ->
function get_normal_variational_covariance(
		beta_mean, beta2_mean, beta_ind_model, beta2_ind_model)

	norm_cov = MatrixTuple[]

	norm_var = beta2_mean - beta_mean ^ 2
	@assert norm_var >= 0.0

	# Get the linear term variance
	push!(norm_cov, (beta_ind_model, beta_ind_model, norm_var))

	# Get the covariance between the linear and quadratic terms.
	this_cov = 2 * beta_mean * norm_var
	push!(norm_cov, (beta_ind_model, beta2_ind_model, this_cov))
	push!(norm_cov, (beta2_ind_model, beta_ind_model, this_cov))

	# Get the covariance between the quadratic terms.
	this_cov = 2 * norm_var ^ 2 + 4 * norm_var * (beta_mean ^ 2)
	push!(norm_cov, (beta2_ind_model, beta2_ind_model, this_cov))

	norm_cov
end

###################################################
# Wishart functions

@doc """
The covariance of wishart distributed random variables.
""" ->
function get_wishart_variational_covariance(
		v0::Matrix{Float64}, wn::Float64, ud_ind::Matrix{Int64})
	@assert size(v0, 1) == size(v0, 2)
	k_tot = size(v0, 1)
	k_ud = round(Int, k_tot * (k_tot + 1) / 2)
	w_cov = Array(Float64, k_ud, k_ud)

	for j1=1:k_tot, i1=1:j1, j2=1:k_tot, i2=1:j2
		ind_1 = ud_ind[j1, i1]
		ind_2 = ud_ind[j2, i2]
		w_cov[ind_1, ind_2] = wn * (v0[i1, j2] * v0[i2, j1] + v0[i1, i2] * v0[j1, j2])
		if ind_1 != ind_2
			w_cov[ind_2, ind_1] = w_cov[ind_1, ind_2]
		end
	end

	w_cov
end

@doc """
The variance of the log determininant of a wishart matrix.
""" ->
function get_wishart_log_det_variance(wn::Float64, k_tot::Int64)
	multivariate_trigamma(float(wn) / 2, k_tot)
end

@doc """
The covariance of the log determininant of a wishart matrix and wishart variables.
""" ->
function get_wishart_cross_variance(v0::Matrix{Float64}, ud_ind::Matrix{Int64})
	2.0 * linearize_matrix(v0, ud_ind)
end

@doc """
The covariance matrix of the sufficient statistics of a wishart distribution.
""" ->
function get_wishart_sufficient_stats_variational_covariance(
		v0::Matrix{Float64}, wn::Float64, lambda_i, log_det_lambda_i, ud_ind)
	@assert size(v0, 1) == size(v0, 2) == size(ud_ind, 1) == size(ud_ind, 2)
	k_tot = size(v0, 1)
	k_ud = k_tot * (k_tot + 1) / 2
	@assert length(lambda_i) == k_ud

	cov_triplets = MatrixTuple[]

	log_det_cov = get_wishart_log_det_variance(wn, k_tot)
	cross_cov = get_wishart_cross_variance(v0, ud_ind)
	lambda_cov = get_wishart_variational_covariance(v0, wn, ud_ind)

	push!(cov_triplets,
		  (log_det_lambda_i, log_det_lambda_i, log_det_cov))

	for i=1:length(lambda_i)
		push!(cov_triplets,
			  (log_det_lambda_i, lambda_i[i], cross_cov[i]))
  		push!(cov_triplets,
			  (lambda_i[i], log_det_lambda_i, cross_cov[i]))
		for j=1:length(lambda_i)
			push!(cov_triplets,
				  (lambda_i[i], lambda_i[j], lambda_cov[i, j]))
		end
	end

	cov_triplets
end

@doc """
For a Wishart distribution with mean wn v0_inv^{-1},
evaluate the expected log determinant.
""" ->
function wishart_e_log_det(wn::Float64, v0_inv::Matrix{Float64})
	@assert size(v0_inv, 1) == size(v0_inv, 2)
	p = size(v0_inv, 1)
	multivariate_digamma(0.5 * wn, p) + p * log(2) - logdet(v0_inv)
end

@doc """
For a Wishart distribution with mean wn v0_inv^{-1},
evaluate the entropy.
""" ->
function wishart_entropy(wn::Float64, v0_inv::Matrix{Float64}, k_tot::Int64)
	-0.5 * logdet(v0_inv) * (k_tot + 1.0) +
	0.5 * k_tot * (k_tot + 1) * log(2) +
	multivariate_lgamma(0.5 * wn, k_tot) -
	0.5 * (wn - k_tot - 1.0) * multivariate_digamma(0.5 * wn, k_tot) +
	0.5 * wn * k_tot
end

function get_wishart_parameters_from_derivs(
		lambda_deriv::Matrix{Float64}, log_det_lambda_deriv::Float64)

	@assert size(lambda_deriv, 1) == size(lambda_deriv, 2)
	k_tot = size(lambda_deriv, 1)
	wn = 2. * log_det_lambda_deriv + 1. + k_tot
	v0 = -0.5 * inv(lambda_deriv)
	v0_inv = -2.0 * lambda_deriv

	wn, v0, v0_inv
end

###################################################
# Gamma distribution functions

function get_gamma_parameters_from_derivs(
		tau_deriv::Float64, log_tau_deriv::Float64)

	tau_alpha = log_tau_deriv + 1
	tau_beta = -tau_deriv

	@assert tau_alpha >= 0
	@assert tau_beta >= 0

	tau_alpha, tau_beta
end

@doc """
Covariance of sufficient statistics for a gamma distribution with
mean tau_alpha / tau_beta.
""" ->
function get_gamma_variational_covariance(tau_alpha::Float64, tau_beta::Float64,
	                                        e_tau_col::Int64, e_log_tau_col::Int64)
	tau_cov = MatrixTuple[]
	push!(tau_cov, (e_tau_col,     e_tau_col,     tau_alpha / (tau_beta ^ 2)))
	push!(tau_cov, (e_log_tau_col, e_log_tau_col, trigamma(tau_alpha)))
	push!(tau_cov, (e_tau_col,     e_log_tau_col, 1 / tau_beta))
	push!(tau_cov, (e_log_tau_col, e_tau_col,     1 / tau_beta))
	tau_cov
end

function gamma_entropy(tau_alpha::Float64, tau_beta::Float64)
	tau_alpha - log(tau_beta) + lgamma(tau_alpha) +
		(1 - tau_alpha) * digamma(tau_alpha)
end

end # module


text:
"""
    addto!(accumulator, newval)

Add newval to accumulator and return the result. Used in outgrad calculations.  The outgrad
values start as `nothing` representing a 0 gradient. Then they are incremented using values
of type Number, Tuple, AbstractDict, AbstractArray, Nothing and AutoGrad.Sparse. The
accumulator and the newval types must match: Nothing matches all types, Sparse matches types
that match its container, other types must match themselves. `addto!` handles repeated
indices in newval by adding all corresponding values to the `accumulator`.
"""
function addto! end

# addto! needs to be a primitive for higher order gradients:
addto!(a::Value,b::Value)=forw(addto!,a,b)
addto!(a::Value,b)=forw(addto!,a,b)
addto!(a,b::Value)=forw(addto!,a,b)
back(::typeof(addto!),::Type{Arg{N}},dy,y,x1,x2) where N = dy

## Types with exact match
addto!(a::Number, b::Number)=a+b
addto!(a::Tuple, b::Tuple)=tuple([addto!(x,y) for (x,y) in zip(a,b)]...)
addto!(a::AbstractDict, b::AbstractDict) = (z=empty(a); for d in (a,b), (k,v) in d; z[k]=addto!(v,get(z,k,nothing)); end; z)
addto!(a::AbstractArray{T},b::AbstractArray) where T = (if isbitstype(T); (a+b); else; T[addto!(x,y) for (x,y) in zip(a,b)]; end)
# We could have Array{Array} and Array{Any} added, so no parametric type on b.

## Nothing indicates zero gradient and matches any type
addto!(::Nothing,::Nothing)=nothing
addto!(a::Value,::Nothing)=a   # to avoid ambiguity
addto!(::Nothing,a::Value)=a   # to avoid ambiguity
addto!(a,::Nothing)=a
addto!(::Nothing,a)=a

## Sparse matches any type that matches its container
matches(::Any,::Any)=false
matches(::Number,::Number)=true
matches(::AbstractDict,::AbstractDict)=true
matches(a::Tuple,b::Tuple)=(length(a)===length(b))
matches(a::AbstractArray,b::AbstractArray)=(size(a)==size(b))

## If both accumulator and newval are sparse, merge, modifying first arg. Use + if you do not want to modify:
function addto!(a::Sparse, b::Sparse)
    @assert matches(a.container, b.container) "$(summary.((a.container, b.container)))"
    append!(a.values, b.values)
    append!(a.indices, b.indices)
    return a
end

## If sparse is the accumulator, reverse:
addto!(a::Sparse,b::Number)=addto!(b,a)
addto!(a::Sparse,b::Tuple)=addto!(b,a)
addto!(a::Sparse,b::AbstractDict)=addto!(b,a)
addto!(a::Sparse,b::AbstractArray)=addto!(b,a)

## Other types with Sparse:

# This comes up if people use getindex on a number:
function addto!(a::Number, b::Sparse)
    @assert matches(a, b.container)
    for (idx, val) in zip(b.indices, b.values)
        if !(idx == (1,) && isa(val,Number))
            throw(ArgumentError("addto!($a,$b)"))
        end
        a += val
    end
    return a
end

# Convert tuples to Any[] and use the array code
function addto!(a::Tuple,b::Sparse)
    @assert matches(a, b.container)
    ca = collect(Any,a)
    for (idx, val) in zip(b.indices, b.values)
        @assert length(idx) == 1
        cb = (idx[1] isa Real ? val : collect(Any,val))
        addtoindex!(ca, cb, to_indices(ca,idx)...)
    end
    tuple(ca...)
end

# Dictionaries just increment for each index,value pair
function addto!(a::AbstractDict,b::Sparse)
    @assert matches(a, b.container)
    if recording(); a = copy(a); end  # do not overwrite array if in highorder context
    for (idx, val) in zip(b.indices, b.values)
        setindex!(a,addto!(get(a,idx...,nothing),val),idx...)
    end
    return a
end

# Need to be careful with arrays because of possible repeated indices
function addto!(a::AbstractArray,b::Sparse)
    @assert matches(a, b.container)
    if recording(); a = copy(a); end  # do not overwrite array if in highorder context
    for (idx, val) in zip(b.indices, b.values)
        addtoindex!(a, val, to_indices(a,idx)...)
    end
    return a
end

# We need the following function to deal with repeated indices.
# Based on base/multidimensional.jl:634 _unsafe_setindex!
# Instead of last value overriding in case of repeated indices, we must sum.

using Base: unalias, index_lengths, setindex_shape_check
using Base.Cartesian # for @nexprs etc.

@generated function addtoindex!(A::AbstractArray, x, I::Union{Real,AbstractArray}...)
    N = length(I)
    quote
        x′ = unalias(A, x)
        @nexprs $N d->(I_d = unalias(A, I[d]))
        idxlens = @ncall $N index_lengths I
        @ncall $N setindex_shape_check x′ (d->idxlens[d])
        Xy = iterate(x′)
        @inbounds @nloops $N i d->I_d begin
            # This is never reached, but serves as an assumption for
            # the optimizer that it does not need to emit error paths
            Xy === nothing && break
            (val, state) = Xy

            ai = @ncall $N getindex A i # <-- different from _unsafe_setindex!
            val = addto!(ai, val) # <-- different from _unsafe_setindex!

            @ncall $N setindex! A val i
            Xy = iterate(x′, state)
        end
        A
    end
end

# The following methods can assume there are no repeated indices:
# Only Array{Int} allows for repeated indices.
addtoindex!(A::AbstractArray, X, I::Union{Real,Colon,AbstractRange,AbstractArray{Bool}}...) =
    setindex!(A, addto!(getindex(A,I...), X), I...)
addtoindex!(A::AbstractArray, X, I...) =
    setindex!(A, addto!(getindex(A,I...), X), I...)


text:
module Segments

import LinearAlgebra: norm, dot
import Statistics: mean, std

using RealNeuralNetworks.Utils.BoundingBoxes
include("Synapses.jl"); using .Synapses

const Node = NTuple{4,Float32}

# TO-DO: make the data type more general
# the AbstractSynapseList is used to define function input type
# the SynapseList is used to create empty synapse list
const AbstractSynapseList = Vector{<:Union{Missing, Vector{Synapse{Float32}}}}
const SynapseList = Vector{Union{Missing, Vector{Synapse{Float32}}}}

# classes following SWC format
const AXON_CLASS = UInt8(2)
const DENDRITE_CLASS = UInt8(3)
const SOMA_CLASS = one(UInt8)
const UNDEFINED_CLASS = zero(UInt8)

export Segment
mutable struct Segment{T}
    # list of tuple (x,y,z,r)
    nodeList        ::Vector{NTuple{4,T}}
    class           ::UInt8
    #boundingBox     ::BoundingBox
    preSynapseList  ::AbstractSynapseList
    postSynapseList ::AbstractSynapseList
end

function Segment()
    nodeList = Vector{Node}()
    Segment(nodeList)
end

function Segment(nodeList::Vector{Node};
                 class::UInt8=UNDEFINED_CLASS,
                 preSynapseList::AbstractSynapseList = SynapseList(missing, length(nodeList)),
                 postSynapseList::AbstractSynapseList = SynapseList(missing, length(nodeList)))
    @assert length(nodeList) == length(preSynapseList)
    @assert length(nodeList) == length(postSynapseList)
    Segment(nodeList, class, preSynapseList, postSynapseList)
end

###################### IO ###########################

###################### properties ###################
"""
    get_features(self::Segment)
aggregate features to a named tuple
"""
function get_features(self::Segment)
    (pathLength=get_path_length(self),
     surfaceArea=get_surface_area(self),
     volume=get_volume(self),
     meanRadius=mean(get_radius_list(self)),
     stdRadius=std(get_radius_list(self)),
     tortuosity=get_tortuosity(self),
     numPreSynapses=get_num_pre_synapses(self),
     numPostSynapses=get_num_post_synapses(self),
     likeFirstAxon=like_first_axon(self))
end

"""
    get_nodes_distance(self::Node, other::Node)
compute the euclidean distance between two nodes
"""
@inline function get_nodes_distance(self::Union{Vector,Tuple}, other::Union{Vector,Tuple})
    norm( [map((x,y)->x-y, self[1:3], other[1:3]) ...])
end

@inline function get_node_list(self::Segment) self.nodeList end
@inline function get_connectivity_matrix( self::Segment ) self.connectivityMatrix end
@inline function get_bounding_box( self::Segment ) BoundingBox(get_node_list(self)) end
@inline function get_class( self::Segment ) self.class end
@inline function get_pre_synapse_list(self::Segment) self.preSynapseList end
@inline function get_post_synapse_list(self::Segment) self.postSynapseList end

function get_synapse_list(self::AbstractSynapseList)
    ret = Synapse{Float32}[]
    for synapses in self
        if !ismissing(synapses)
            for synapse in synapses
                push!(ret, synapse)
            end
        end
    end
    ret

end

@inline function concatenate_pre_synapses( self::Segment{T} ) where T
    get_synapse_list(self.preSynapseList)
end

@inline function concatenate_post_synapses( self::Segment{T} ) where T
    get_synapse_list(self.postSynapseList)
end

@inline function get_pre_synapse( self::Segment, index::Int ) self.preSynapseList[index] end
@inline function get_post_synapse( self::Segment, index::Int ) self.postSynapseList[index] end

@inline function get_num_pre_synapses(self::Segment)
    concatenate_pre_synapses(self) |> length
end

@inline function get_num_post_synapses(self::Segment)
    concatenate_post_synapses(self) |> length
end

"""
    get_pre_synapse_density(self::Segment)
note that the unit is # / micron
"""
@inline function get_pre_synapse_density(self::Segment)
    numPreSynapses = get_num_pre_synapses(self)
    pathLength = get_path_length(self)
    numPreSynapses / pathLength * 1000
end

"""
    get_post_synapse_density(self::Segment)
note that the unit is # / micron
"""
@inline function get_post_synapse_density(self::Segment)
    numPostSynapses = get_num_post_synapses(self)
    pathLength = get_path_length(self)
    numPostSynapses / pathLength * 1000
end

"""
    get_central_pre_synapse_num(self::Segment{T}) where T

some axons in zebrafish start with some post synapses which is a great confusion of
segment type classification. only considering the central part should benefit the classification.
"""
function get_central_pre_synapse_num(self::Segment{T};
                                        start::Float64=0.2, stop::Float64=0.8) where T
    @assert start < stop
    @assert start>0 && stop<1
    start = max(1, round(Int, start * length(self)))
    stop = min(length(self), round(Int, length(self)*stop))
    self.preSynapseList[start:stop] |> get_synapse_list |> length
end

function get_central_post_synapse_num(self::Segment{T};
                                        start::Float64=0.2, stop::Float64=0.8) where T
    @assert start < stop
    @assert start>=0. && stop<=1.
    start = max(1, round(Int, start * length(self)))
    stop = min(length(self), round(Int, length(self)*stop))
    self.postSynapseList[start:stop] |> get_synapse_list |> length
end

"""
    like_first_axon(self::Segment{T}; ratio::AbstractFloat=0.3) where T

first axon might have a few postsynapses at the starting part.
"""
function like_first_axon(self::Segment{T}; ratio::AbstractFloat=0.3) where T
    @assert ratio>=0 && ratio<=1

    synNum11 = get_central_post_synapse_num(self; start=0., stop=ratio)
    synNum12 = get_central_post_synapse_num(self; start=ratio, stop=1.)
    if synNum11 > 1 && synNum12==0
        return true
    end

    synNum21 = get_central_post_synapse_num(self; start=0.0, stop=1.0-ratio)
    synNum22 = get_central_post_synapse_num(self; start=1.0-ratio, stop=1.0)
    if synNum21==0 && synNum22>1
        return true
    end
    false
end

@inline function set_class(self::Segment, class::UInt8) self.class=class end

@inline function get_bounding_box_distance(self::Segment, point::Union{Tuple, Vector})
    @assert length(point) >= 3
    boundingBox = get_bounding_box(self)
    BoundingBoxes.distance_from(boundingBox, point)
end

@inline function euclidean_distance( n1::NTuple{3,T}, n2::NTuple{3,T}) where T
    norm([map((x,y)->x-y, n1, n2)...])
end

"""
    get_path_length(self::Segment; nodeId::Int=length(self))
accumulate the euclidean distance between neighboring nodes
"""
@inline function get_path_length(self::Segment; nodeId::Int=length(self))
    ret = 0.0
    for i in 2:nodeId
        ret += euclidean_distance(self[i][1:3], self[i-1][1:3])
    end
    ret
end

@inline function get_radius_list( self::Segment ) map(n->n[4], self) end

"""
    get_surface_area(self::Segment)
frustum-based:
http://www.analyzemath.com/Geometry_calculators/surface_volume_frustum.html
"""
function get_surface_area(self::Segment{T}) where T
    ret = zero(T)
    for i in 2:length(self)
        # average diameter
        r1 = self[i][4]
        r2 = self[i-1][4]
        h = euclidean_distance(self[i][1:3], self[i-1][1:3])
        ret += pi* (r1+r2) * sqrt(h*h + (r1-r2)*(r1-r2))
    end
    ret::T
end

"""
    get_volume(self::Segment)
compute frustum-based volume
http://jwilson.coe.uga.edu/emt725/Frustum/Frustum.cone.html
"""
function get_volume(self::Segment{T}) where T
    ret = zero(T)
    for i in 2:length(self)
        r1 = self[i-1][4]
        r2 = self[i][4]
        h = euclidean_distance(self[i-1][1:3], self[i][1:3])
        ret += pi * h * (r1*r1 + r1*r2 + r2*r2) / T(3)
    end
    ret
end

"""
    get_tail_head_radius_ratio( self::Segment )
the spine is normally thick in tail, and thin in the head.
ratio = max_tail / mean_head
The head should point to dendrite. This is a very good feature to identify spine.
"""
@inline function get_tail_head_radius_ratio( self::Segment )
    radiusList = get_radius_list( self )
    N = length(self)
    headRadiusList = radiusList[1:cld(N,2)]
    tailRadiusList = radiusList[cld(N,2):N]
    maximum(tailRadiusList) / mean(headRadiusList)
end

"""
    get_tortuosity( self::Segment )
the ratio of the actual path length to the euclidean distance between head and tail node
"""
@inline function get_tortuosity(self::Segment)
    if length(self) == 1
        return 1.0
    end
    pathLength = get_path_length(self)
    euclideanLength = get_nodes_distance( self[1], self[end] )
    @assert self[1]!=self[end] "segment start is the same with the end: $(self)"
    @assert euclideanLength != 0.0
    pathLength / euclideanLength
end

@inline function get_center(nodeList::Vector{Node})
    (map(i->mean(y->y[i], nodeList), 1:4)...,)
end

@inline function get_center(self::Segment, range::UnitRange)
    center = get_center( self[range] )
end

###################### Base functions ################

function Base.show(self::Segment)
    if get_class(self) == AXON_CLASS
        class = "axon"
    elseif get_class(self) == DENDRITE_CLASS
        class = "dendrite"
    else
        class = "unknown"
    end
    println("segment is a ", class, " and have ", length(self), " nodes, ",
            get_num_pre_synapses(self), " presynapses, ",
            get_num_post_synapses(self), " postsynapses.")
    nothing
end

function Base.iterate(self::Segment, state::Int=1)
    if state > length(self)
        return nothing
    else
        return get_node_list(self)[state], state+1
    end
end

@inline function Base.lastindex(self::Segment) length(self) end
@inline function Base.isempty(self::Segment) isempty(self.nodeList) end

"""
    Base.length(self::Segment)

the number of nodes contained in this segment
"""
function Base.length(self::Segment)
    length(self.nodeList)
end

"""
    Base.merge(self::Segment, other::Segment)
merge two segmentes
"""
function Base.merge(self::Segment, other::Segment)
    nodeList1 = get_node_list(self)
    nodeList2 = get_node_list(other)
    mergedNodeList = vcat( nodeList1, nodeList2 )
    # winner taks all!
    class = length(nodeList1)>length(nodeList2) ? get_class(self) : get_class(other)
    # merge synapse list
    preSynapseList1 = get_pre_synapse_list(self)
    preSynapseList2 = get_pre_synapse_list(other)
    postSynapseList1 = get_post_synapse_list(self)
    postSynapseList2 = get_post_synapse_list(other)

    mergedPreSynapseList = vcat( preSynapseList1, preSynapseList2 )
    mergedPostSynapseList = vcat( postSynapseList1, postSynapseList2 )

    Segment(mergedNodeList; class=class, preSynapseList=mergedPreSynapseList,
            postSynapseList=mergedPostSynapseList)
end

"""
split the segment from the node list index to two segmentes
the indexed node will be included in the second segment
"""
function Base.split(self::Segment{T}, index::Integer) where T
    @assert index >=1 && index<=length(self)
    idx1 = index - 1
    idx2 = index

    if index==1
        if length(self)==1
            idx1 = 1
            idx2 = 1
        else
            idx1 = index
            idx2 = index + 1
        end
    end
    nodeList1 = self.nodeList[1:idx1]
    nodeList2 = self.nodeList[idx2:end]
    preSynapseList1 = self.preSynapseList[1:idx1]
    preSynapseList2 = self.preSynapseList[idx2:end]
    postSynapseList1 = self.postSynapseList[1:idx1]
    postSynapseList2 = self.postSynapseList[idx2:end]

    segment1 = Segment(nodeList1; class=self.class, preSynapseList=preSynapseList1, postSynapseList=postSynapseList1)
    segment2 = Segment(nodeList2; class=self.class, preSynapseList=preSynapseList2, postSynapseList=postSynapseList2)
    @assert length(nodeList1) > 0
    @assert length(nodeList2) > 0
    return segment1, segment2
end

"""
    Base.getindex(self::Segment, index::Integer)
"""
@inline function Base.getindex(self::Segment, index::Integer)
    get_node_list(self)[ index ]
end

@inline function Base.getindex(self::Segment, range::UnitRange)
    get_node_list(self)[range]
end

"""
distance from a point
"""
function distance_from(self::Segment, point::Tuple)
    distance_from(self, [point[1:3]...])
end
function distance_from(self::Segment{T}, point::Vector) where T
    @assert !isempty(self)
    ret = (zero(T), zero(Int))
    nodeList = get_node_list(self)
    @assert !isempty(nodeList)
    @assert length(point) == 3 || length(point) == 4
    distance = typemax(T)
    for (index, node) in enumerate(nodeList)
        d = norm( [node[1:3]...] .- [point[1:3]...] )
        if d < distance
            distance = d
            ret = (d, index)
        end
    end
    @assert ret[2] <= length(self)
    ret
end

################## manipulation ###############################

"""
    attach_pre_synapse!(self::Segment, nodeId::Int, synapse::Synapse{T})
"""
@inline function attach_pre_synapse!(self::Segment, nodeId::Int, synapse::Synapse{T}) where T
    if ismissing(self.preSynapseList[ nodeId ])
        self.preSynapseList[ nodeId ] = [synapse]
    else
        push!(self.preSynapseList[nodeId], synapse)
    end
    nothing
end

"""
    attach_post_synapse!(self::Segment, nodeId::Int, synapse::Synapse)
"""
@inline function attach_post_synapse!(self::Segment, nodeId::Int, synapse::Synapse{T}) where T
    if ismissing(self.postSynapseList[nodeId])
        self.postSynapseList[ nodeId ] = [synapse]
    else
        push!(self.postSynapseList[nodeId], synapse)
    end
    nothing
end

function adjust_class!(self::Segment)
    numPreSynapses = get_num_pre_synapses(self)
    numPostSynapses = get_num_post_synapses(self)
    if self.class == UNDEFINED_CLASS
        if numPreSynapses > numPostSynapses
            # mostly presynapses, so this is an axon
            # Note that this simple metric fails in the axonal hillock,
            # where there exist a lot of post synapses.
            self.class = AXON_CLASS
        elseif numPreSynapses < numPostSynapses
            # mostly postsynapses, so this is a dendrite
            self.class = DENDRITE_CLASS
        elseif numPostSynapses == 0 && get_path_length(self)>5000
            # a long segment without both pre and post synapses was considered axon
            self.class = AXON_CLASS
        else
            self.class = DENDRITE_CLASS
        end
    end
end

function add_offset(self::Segment{T}, offset::Union{Tuple, Vector}) where T
    @assert length(offset) == 3
    nodeList = Vector{NTuple{4,T}}()
    for node in self.nodeList
        newNode = map(+, node, [offset..., zero(T)])
        push!(nodeList, newNode)
    end
    Segment(nodeList, self.class)
end

@inline function remove_node(self::Segment, removeId::Int)
    remove_nodes(self, removeId:removeId)
end

"""
    remove_nodes(self::Segment, removeIdRange::UnitRange{Int})
remove nodes from a segment
"""
function remove_nodes(self::Segment{T}, removeIdRange::UnitRange{Int}) where T
    newLength = length(self) - length(removeIdRange)
    @assert newLength >= 0
    if newLength == 0
        return Segment()
    end
    newNodeList = Vector{NTuple{4, T}}()
    sizehint!(newNodeList, newLength)

    for (index,node) in enumerate(get_node_list(self))
        if !(index in removeIdRange)
            push!(newNodeList, node)
        end
    end

    preSynapseList = get_pre_synapse_list(self)
    postSynapseList = get_post_synapse_list(self)
    newPreSynapseList = SynapseList(missing, newLength)
    newPostSynapseList = SynapseList(missing, newLength)

    newPreSynapseList[1:removeIdRange.start-1]  = preSynapseList[1:removeIdRange.start-1]
    newPostSynapseList[1:removeIdRange.start-1] = postSynapseList[1:removeIdRange.start-1]

    newPreSynapseList[removeIdRange.start:end]  = preSynapseList[removeIdRange.stop+1:end]
    newPostSynapseList[removeIdRange.start:end] = postSynapseList[removeIdRange.stop+1:end]

    Segment(newNodeList; class = get_class(self),
            preSynapseList=newPreSynapseList, postSynapseList=newPostSynapseList)
end

"""
    remove_redundent_nodes!(self::Segment)
remove neighboring nodes that is the same.
"""
function remove_redundent_nodes!(self::Segment{T}) where T
    nodeList = get_node_list(self)
    newNodeList = Vector{NTuple{4, T}}()
    newPreSynapseList = SynapseList()
    newPostSynapseList = SynapseList()
    for index in 1:length(nodeList)-1
        if nodeList[index] != nodeList[index+1]
            push!(newNodeList, nodeList[index])
            push!(newPreSynapseList, self.preSynapseList[index])
            push!(newPostSynapseList, self.postSynapseList[index])
        else
            # merge synapses
            newPreSynapseList[index] = (newPreSynapseList[index]..., self.preSynapseList[index+1]...)
            newPostSynapseList[index] = (newPostSynapseList[index]..., self.postSynapseList[index+1]...)
        end
    end
    # get the final node
    push!(newNodeList, nodeList[end])
    self.nodeList = newNodeList
end

end # module


text:
abstract HTMLNode

# TODO immutable?
type HTMLText <: HTMLNode
    parent::HTMLNode
    text::String
end

# convenience method for defining without parent
HTMLText(text::String) = HTMLText(NullNode(), text)

type NullNode <: HTMLNode end

type HTMLElement{T} <: HTMLNode
    children::Vector{HTMLNode}
    parent::HTMLNode
    attributes::Dict{String, String}
end

# convenience method for defining an empty element
HTMLElement(T::Symbol) = HTMLElement{T}(HTMLNode[],NullNode(),Dict{String,String}())

type HTMLDocument
    doctype::String
    root::HTMLElement
end

type InvalidHTMLException <: Exception
    msg::String
end


text:
using PowerDynamics: AbstractLine
###### knoten ######

function CIGRE_static()

    LoadP = [
        19839E3,
        0.0,
        501.7E3,
        431.65E3,
        727.5E3,
        548.05E3,
        76.5E3,
        586.85E3,
        573.75E3,
        543.3E3,
        329.8E3,
    ]

    cosϕ = [
        0.9737554,
        0.0,
        0.9231825,
        0.9700009,
        0.9699996,
        0.9700017,
        0.8500022,
        0.9699995,
        0.8499989,
        0.9586623,
        0.969997,
    ]

    LoadQ = LoadP .* sin.(acos.(cosϕ)) ./ cosϕ
    LoadQ[2] = 0.0

    begin
        busses_static = Array{PowerDynamics.AbstractNode,1}([])
        push!(busses_static, SlackAlgebraic(U = 110E3 / base_voltage_HV))

        for (P, Q) in zip(LoadP, LoadQ)
            try
                push!(
                    busses_static,
                    PQAlgebraic(S = - complex(P, Q) ./ base_power),
                )
            catch
                push!(
                    busses_static,
                    PQAlgebraic(P = - P / base_power, Q = - Q / base_power),
                )
            end
        end
    end

    ###### Kanten ######

    T = OLTC(
        from = 1,
        to = 2,
        Uos = 110E3 / base_voltage_HV, # Bemessungsspannung Oberspannungsseite in kV
        Uus = 20E3 / base_voltage,# Bemessungsspannung Unterspannungsseite in kV
        k = 0, # Kennzahl der Schaltgruppe, Bsp: YD5-> 5
        ssp = 10.0, # Stufenschalterposition
        stufe_us = 0.625, # Stufung pro Stufenschalterstellung in %
        Sr = 25E6 / base_power, #Bemessungsscheinleistung in MVA
        uk = 12.0, # Kurzschlussspannung in %
        Pvk = 25E3 / base_power, # Kupferverluste in MW
        Pvl = 0.0, # Eisenverluste in kW
        iLeer = 0.0, # Leerlaufstrom in %
    )

    #ldata = CSV.read("$dir/lines.csv"; header = true)
    ldata = [ #L_km(nicht ändern)
        2.82,
        4.42,
        0.61,
        1.3,
        0.56,
        1.54,
        1.67,
        0.32,
        0.77,
        0.33,
        0.24,
    ]

    elist = [
        (2, 3),
        (3, 4),
        (4, 5),
        (4, 9),
        (5, 6),
        (6, 7),
        (8, 9),
        (9, 10),
        (10, 11),
        (11, 12),
        #(7, 8), # Schalter
        #(12, 5), # Schalter
    ]

    R1 = 0.501 # Ω/km
    X1 = 0.716 # Ω/km
    C1 = 0.1511749E-6 # F/km

    Zs = complex(R1, X1) .* ldata #[!, Symbol("L_km(nicht ändern)")]
    Yshs = 1im .* ω .* C1 .* ldata #[!, Symbol("L_km(nicht ändern)")]

    begin
        lines = Array{AbstractLine,1}([])
        push!(lines, T)
        for (e, Z, Ysh) in zip(elist, Zs, Yshs)
            push!(
                lines,
                PiModelLine(
                    from = first(e),
                    to = last(e),
                    y = inv(Z) / base_admittance,
                    y_shunt_km = Ysh / 2.0 / base_admittance,
                    y_shunt_mk = Ysh / 2.0 / base_admittance,
                ),
            )
        end
    end

    busses_static, lines, T, elist, Zs, Yshs
end

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
using TimeSeries
using Dates
using Random
Random.seed!(123)
using PowerSystems
const PSY = PowerSystems

DayAhead = collect(
    DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S"):Hour(1):DateTime(
        "1/1/2024  23:00:00",
        "d/m/y  H:M:S",
    ),
)
#Dispatch_11am =  collect(DateTime("1/1/2024  0:11:00", "d/m/y  H:M:S"):Minute(15):DateTime("1/1/2024  12::00", "d/m/y  H:M:S"))

nodes5() = [
    Bus(1, "nodeA", "PV", 0, 1.0, (min = 0.9, max = 1.05), 230, nothing, nothing),
    Bus(2, "nodeB", "PQ", 0, 1.0, (min = 0.9, max = 1.05), 230, nothing, nothing),
    Bus(3, "nodeC", "PV", 0, 1.0, (min = 0.9, max = 1.05), 230, nothing, nothing),
    Bus(4, "nodeD", "REF", 0, 1.0, (min = 0.9, max = 1.05), 230, nothing, nothing),
    Bus(5, "nodeE", "PV", 0, 1.0, (min = 0.9, max = 1.05), 230, nothing, nothing),
];

branches5_dc(nodes5) = [
    Line(
        "1",
        true,
        0.0,
        0.0,
        Arc(from = nodes5[1], to = nodes5[2]),
        0.00281,
        0.0281,
        (from = 0.00356, to = 0.00356),
        2.0,
        (min = -0.7, max = 0.7),
    ),
    HVDCLine(
        "DCL2",
        true,
        0.0,
        Arc(from = nodes5[1], to = nodes5[4]),
        (min = -3000.0, max = 3000.0),
        (min = -3000, max = 3000),
        (min = -3000.0, max = 3000.0),
        (min = -3000.0, max = 3000.0),
        (l0 = 0.0, l1 = 0.01),
    ),
    Line(
        "3",
        true,
        0.0,
        0.0,
        Arc(from = nodes5[1], to = nodes5[5]),
        0.00064,
        0.0064,
        (from = 0.01563, to = 0.01563),
        18.8120,
        (min = -0.7, max = 0.7),
    ),
    Line(
        "4",
        true,
        0.0,
        0.0,
        Arc(from = nodes5[2], to = nodes5[3]),
        0.00108,
        0.0108,
        (from = 0.00926, to = 0.00926),
        11.1480,
        (min = -0.7, max = 0.7),
    ),
    Line(
        "5",
        true,
        0.0,
        0.0,
        Arc(from = nodes5[3], to = nodes5[4]),
        0.00297,
        0.0297,
        (from = 0.00337, to = 0.00337),
        40.530,
        (min = -0.7, max = 0.7),
    ),
    Line(
        "6",
        true,
        0.0,
        0.0,
        Arc(from = nodes5[4], to = nodes5[5]),
        0.00297,
        0.0297,
        (from = 0.00337, to = 00.00337),
        2.00,
        (min = -0.7, max = 0.7),
    ),
];

branches5(nodes5) = [
    Line(
        "1",
        true,
        0.0,
        0.0,
        Arc(from = nodes5[1], to = nodes5[2]),
        0.00281,
        0.0281,
        (from = 0.00356, to = 0.00356),
        2.0,
        (min = -0.7, max = 0.7),
    ),
    Line(
        "2",
        true,
        0.0,
        0.0,
        Arc(from = nodes5[1], to = nodes5[4]),
        0.00304,
        0.0304,
        (from = 0.00329, to = 0.00329),
        2.0,
        (min = -0.7, max = 0.7),
    ),
    Line(
        "3",
        true,
        0.0,
        0.0,
        Arc(from = nodes5[1], to = nodes5[5]),
        0.00064,
        0.0064,
        (from = 0.01563, to = 0.01563),
        18.8120,
        (min = -0.7, max = 0.7),
    ),
    Line(
        "4",
        true,
        0.0,
        0.0,
        Arc(from = nodes5[2], to = nodes5[3]),
        0.00108,
        0.0108,
        (from = 0.00926, to = 0.00926),
        11.1480,
        (min = -0.7, max = 0.7),
    ),
    Line(
        "5",
        true,
        0.0,
        0.0,
        Arc(from = nodes5[3], to = nodes5[4]),
        0.00297,
        0.0297,
        (from = 0.00337, to = 0.00337),
        40.530,
        (min = -0.7, max = 0.7),
    ),
    Line(
        "6",
        true,
        0.0,
        0.0,
        Arc(from = nodes5[4], to = nodes5[5]),
        0.00297,
        0.0297,
        (from = 0.00337, to = 00.00337),
        2.00,
        (min = -0.7, max = 0.7),
    ),
];

branches5_ml(nodes5) = [
    MonitoredLine(
        "1",
        true,
        0.0,
        0.0,
        Arc(from = nodes5[1], to = nodes5[2]),
        0.00281,
        0.0281,
        (from = 0.00356, to = 0.00356),
        (from_to = 1.0, to_from = 1.0),
        2.0,
        (min = -0.7, max = 0.7),
    ),
    Line(
        "2",
        true,
        0.0,
        0.0,
        Arc(from = nodes5[1], to = nodes5[4]),
        0.00304,
        0.0304,
        (from = 0.00329, to = 0.00329),
        2.0,
        (min = -0.7, max = 0.7),
    ),
    Line(
        "3",
        true,
        0.0,
        0.0,
        Arc(from = nodes5[1], to = nodes5[5]),
        0.00064,
        0.0064,
        (from = 0.01563, to = 0.01563),
        18.8120,
        (min = -0.7, max = 0.7),
    ),
    Line(
        "4",
        true,
        0.0,
        0.0,
        Arc(from = nodes5[2], to = nodes5[3]),
        0.00108,
        0.0108,
        (from = 0.00926, to = 0.00926),
        11.1480,
        (min = -0.7, max = 0.7),
    ),
    Line(
        "5",
        true,
        0.0,
        0.0,
        Arc(from = nodes5[3], to = nodes5[4]),
        0.00297,
        0.0297,
        (from = 0.00337, to = 0.00337),
        40.530,
        (min = -0.7, max = 0.7),
    ),
    Line(
        "6",
        true,
        0.0,
        0.0,
        Arc(from = nodes5[4], to = nodes5[5]),
        0.00297,
        0.0297,
        (from = 0.00337, to = 00.00337),
        2.00,
        (min = -0.7, max = 0.7),
    ),
];

solar_ts_DA = [
    0
    0
    0
    0
    0
    0
    0
    0
    0
    0.351105684
    0.632536266
    0.99463925
    1
    0.944237283
    0.396681234
    0.366511428
    0.155125829
    0.040872694
    0
    0
    0
    0
    0
    0
]

wind_ts_DA = [
    0.985205412
    0.991791369
    0.997654144
    1
    0.998663733
    0.995497149
    0.992414567
    0.98252418
    0.957203427
    0.927650911
    0.907181989
    0.889095913
    0.848186718
    0.766813846
    0.654052531
    0.525336131
    0.396098004
    0.281771509
    0.197790004
    0.153241012
    0.131355854
    0.113688144
    0.099302656
    0.069569628
]

hydro_inflow_ts_DA = [
    0.314300
    0.386684
    0.228582
    0.226677
    0.222867
    0.129530
    0.144768
    0.365731
    0.207628
    0.622885
    0.670507
    0.676221
    0.668602
    0.407638
    0.321919
    0.369541
    0.287632
    0.449544
    0.630505
    0.731462
    0.777178
    0.712413
    0.780988
    0.190485
];

thermal_generators5(nodes5) = [
    ThermalStandard(;
        name = "Alta",
        available = true,
        status = true,
        bus = nodes5[1],
        active_power = 0.40,
        reactive_power = 0.010,
        rating = 0.5,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 0.0, max = 0.40),
        reactive_power_limits = (min = -0.30, max = 0.30),
        ramp_limits = nothing,
        time_limits = nothing,
        operation_cost = ThreePartCost((0.0, 14.0), 0.0, 4.0, 2.0),
        base_power = 100.0,
    ),
    ThermalStandard(
        name = "Park City",
        available = true,
        status = true,
        bus = nodes5[1],
        active_power = 1.70,
        reactive_power = 0.20,
        rating = 2.2125,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 0.0, max = 1.70),
        reactive_power_limits = (min = -1.275, max = 1.275),
        ramp_limits = (up = 0.02 * 2.2125, down = 0.02 * 2.2125),
        time_limits = (up = 2.0, down = 1.0),
        operation_cost = ThreePartCost((0.0, 15.0), 0.0, 1.5, 0.75),
        base_power = 100.0,
    ),
    ThermalStandard(
        name = "Solitude",
        available = true,
        status = true,
        bus = nodes5[3],
        active_power = 5.2,
        reactive_power = 1.00,
        rating = 5.2,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 0.0, max = 5.20),
        reactive_power_limits = (min = -3.90, max = 3.90),
        ramp_limits = (up = 0.012 * 5.2, down = 0.012 * 5.2),
        time_limits = (up = 3.0, down = 2.0),
        operation_cost = ThreePartCost((0.0, 30.0), 0.0, 3.0, 1.5),
        base_power = 100.0,
    ),
    ThermalStandard(
        name = "Sundance",
        available = true,
        status = true,
        bus = nodes5[4],
        active_power = 2.0,
        reactive_power = 0.40,
        rating = 2.5,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 0.0, max = 2.0),
        reactive_power_limits = (min = -1.5, max = 1.5),
        ramp_limits = (up = 0.015 * 2.5, down = 0.015 * 2.5),
        time_limits = (up = 2.0, down = 1.0),
        operation_cost = ThreePartCost((0.0, 40.0), 0.0, 4.0, 2.0),
        base_power = 100.0,
    ),
    ThermalStandard(
        name = "Brighton",
        available = true,
        status = true,
        bus = nodes5[5],
        active_power = 6.0,
        reactive_power = 1.50,
        rating = 0.75,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 0.0, max = 6.0),
        reactive_power_limits = (min = -4.50, max = 4.50),
        ramp_limits = (up = 0.015 * 7.5, down = 0.015 * 7.5),
        time_limits = (up = 5.0, down = 3.0),
        operation_cost = ThreePartCost((0.0, 10.0), 0.0, 1.5, 0.75),
        base_power = 100.0,
    ),
];

thermal_generators5_pwl(nodes5) = [
    ThermalStandard(
        name = "Test PWL",
        available = true,
        status = true,
        bus = nodes5[1],
        active_power = 1.70,
        reactive_power = 0.20,
        rating = 2.2125,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 0.0, max = 1.70),
        reactive_power_limits = (min = -1.275, max = 1.275),
        ramp_limits = (up = 0.02 * 2.2125, down = 0.02 * 2.2125),
        time_limits = (up = 2.0, down = 1.0),
        operation_cost = ThreePartCost(
            VariableCost([(0.0, 50.0), (190.1, 80.0), (582.72, 120.0), (1094.1, 170.0)]),
            0.0,
            1.5,
            0.75,
        ),
        base_power = 100.0,
    ),
];

thermal_generators5_pwl_nonconvex(nodes5) = [
    ThermalStandard(
        name = "Test PWL Nonconvex",
        available = true,
        status = true,
        bus = nodes5[1],
        active_power = 1.70,
        reactive_power = 0.20,
        rating = 2.2125,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 0.0, max = 1.70),
        reactive_power_limits = (min = -1.275, max = 1.275),
        ramp_limits = (up = 0.02 * 2.2125, down = 0.02 * 2.2125),
        time_limits = (up = 2.0, down = 1.0),
        operation_cost = ThreePartCost(
            VariableCost([(0.0, 50.0), (190.1, 80.0), (582.72, 120.0), (825.1, 170.0)]),
            0.0,
            1.5,
            0.75,
        ),
        base_power = 100.0,
    ),
];

thermal_pglib_generators5(nodes5) = [
    ThermalMultiStart(
        "115_STEAM_1",
        true,
        true,
        nodes5[1],
        0.05,
        0.010,
        0.12,
        PrimeMovers.ST,
        ThermalFuels.COAL,
        (min = 0.05, max = 0.12),
        (min = -0.30, max = 0.30),
        (up = 0.002 * 0.12, down = 0.002 * 0.12),
        (startup = 0.05, shutdown = 0.05),
        (up = 4.0, down = 2.0),
        (hot = 2.0, warm = 4.0, cold = 12.0),
        3,
        MultiStartCost(
            VariableCost([(0.0, 5.0), (290.1, 7.33), (582.72, 9.67), (894.1, 12.0)]),
            897.29,
            0.0,
            (hot = 393.28, warm = 455.37, cold = 703.76),
            0.0,
        ),
        100.0,
    ),
    ThermalMultiStart(
        "101_CT_1",
        true,
        true,
        nodes5[1],
        0.08,
        0.020,
        0.12,
        PrimeMovers.ST,
        ThermalFuels.COAL,
        (min = 0.08, max = 0.20),
        (min = -0.30, max = 0.30),
        (up = 0.002 * 0.2, down = 0.002 * 0.2),
        (startup = 0.08, shutdown = 0.08),
        (up = 1.0, down = 1.0),
        (hot = 1.0, warm = 999.0, cold = 999.0),
        1,
        MultiStartCost(
            VariableCost([(0.0, 8.0), (391.45, 12.0), (783.74, 16.0), (1212.28, 20.0)]),
            1085.78,
            0.0,
            (hot = 51.75, warm = PSY.START_COST, cold = PSY.START_COST),
            0.0,
        ),
        100.0,
    ),
];

thermal_generators5_uc_testing(nodes) = [
    ThermalStandard(
        name = "Alta",
        available = true,
        status = false,
        bus = nodes[1],
        active_power = 0.0,
        reactive_power = 0.0,
        rating = 0.5,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 0.2, max = 0.40),
        reactive_power_limits = (min = -0.30, max = 0.30),
        ramp_limits = (up = 0.40, down = 0.40),
        time_limits = (up = 0.0, down = 0.0),
        operation_cost = ThreePartCost((0.0, 14.0), 0.0, 4.0, 2.0),
        base_power = 100.0,
    ),
    ThermalStandard(
        name = "Park City",
        available = true,
        status = false,
        bus = nodes[1],
        active_power = 0.0,
        reactive_power = 0.0,
        rating = 2.2125,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 0.65, max = 1.70),
        reactive_power_limits = (min = -1.275, max = 1.275),
        ramp_limits = (up = 0.02 * 2.2125, down = 0.02 * 2.2125),
        time_limits = (up = 0.0, down = 0.0),
        operation_cost = ThreePartCost((0.0, 15.0), 0.0, 1.5, 0.75),
        base_power = 100.0,
    ),
    ThermalStandard(
        name = "Solitude",
        available = true,
        status = true,
        bus = nodes[3],
        active_power = 2.7,
        reactive_power = 0.00,
        rating = 5.20,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 1.0, max = 5.20),
        reactive_power_limits = (min = -3.90, max = 3.90),
        ramp_limits = (up = 0.0012 * 5.2, down = 0.0012 * 5.2),
        time_limits = (up = 5.0, down = 3.0),
        operation_cost = ThreePartCost((0.0, 30.0), 0.0, 3.0, 1.5),
        base_power = 100.0,
    ),
    ThermalStandard(
        name = "Sundance",
        available = true,
        status = false,
        bus = nodes[4],
        active_power = 0.0,
        reactive_power = 0.00,
        rating = 2.5,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 1.0, max = 2.0),
        reactive_power_limits = (min = -1.5, max = 1.5),
        ramp_limits = (up = 0.015 * 2.5, down = 0.015 * 2.5),
        time_limits = (up = 2.0, down = 1.0),
        operation_cost = ThreePartCost((0.0, 40.0), 0.0, 4.0, 2.0),
        base_power = 100.0,
    ),
    ThermalStandard(
        name = "Brighton",
        available = true,
        status = true,
        bus = nodes[5],
        active_power = 6.0,
        reactive_power = 0.0,
        rating = 7.5,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 3.0, max = 6.0),
        reactive_power_limits = (min = -4.50, max = 4.50),
        ramp_limits = (up = 0.0015 * 7.5, down = 0.0015 * 7.5),
        time_limits = (up = 5.0, down = 3.0),
        operation_cost = ThreePartCost((0.0, 10.0), 0.0, 1.5, 0.75),
        base_power = 100.0,
    ),
];

renewable_generators5(nodes5) = [
    RenewableDispatch(
        "WindBusA",
        true,
        nodes5[5],
        0.0,
        0.0,
        1.200,
        PrimeMovers.WT,
        (min = 0.0, max = 0.0),
        1.0,
        TwoPartCost(0.220, 0.0),
        100.0,
    ),
    RenewableDispatch(
        "WindBusB",
        true,
        nodes5[4],
        0.0,
        0.0,
        1.200,
        PrimeMovers.WT,
        (min = 0.0, max = 0.0),
        1.0,
        TwoPartCost(0.220, 0.0),
        100.0,
    ),
    RenewableDispatch(
        "WindBusC",
        true,
        nodes5[3],
        0.0,
        0.0,
        1.20,
        PrimeMovers.WT,
        (min = -0.800, max = 0.800),
        1.0,
        TwoPartCost(0.220, 0.0),
        100.0,
    ),
];

hydro_generators5(nodes5) = [
    HydroDispatch(
        name = "HydroDispatch",
        available = true,
        bus = nodes5[2],
        active_power = 0.0,
        reactive_power = 0.0,
        rating = 6.0,
        prime_mover = PrimeMovers.HY,
        active_power_limits = (min = 0.0, max = 6.0),
        reactive_power_limits = (min = 0.0, max = 6.0),
        ramp_limits = nothing,
        time_limits = nothing,
        base_power = 100.0,
    ),
    HydroEnergyReservoir(
        name = "HydroEnergyReservoir",
        available = true,
        bus = nodes5[3],
        active_power = 0.0,
        reactive_power = 0.0,
        rating = 7.0,
        prime_mover = PrimeMovers.HY,
        active_power_limits = (min = 0.0, max = 7.0),
        reactive_power_limits = (min = 0.0, max = 7.0),
        ramp_limits = (up = 7.0, down = 7.0),
        time_limits = nothing,
        operation_cost = TwoPartCost(0.15, 0.0),
        base_power = 100.0,
        storage_capacity = 50.0,
        inflow = 4.0,
        conversion_factor = 1.0,
        initial_storage = 0.5,
    ),
];

hydro_generators5_ems(nodes5) = [
    HydroDispatch(
        name = "HydroDispatch",
        available = true,
        bus = nodes5[2],
        active_power = 0.0,
        reactive_power = 0.0,
        rating = 6.0,
        prime_mover = PrimeMovers.HY,
        active_power_limits = (min = 0.0, max = 6.0),
        reactive_power_limits = (min = 0.0, max = 6.0),
        ramp_limits = nothing,
        time_limits = nothing,
        base_power = 100.0,
    ),
    HydroEnergyReservoir(
        name = "HydroEnergyReservoir",
        available = true,
        bus = nodes5[3],
        active_power = 0.0,
        reactive_power = 0.0,
        rating = 7.0,
        prime_mover = PrimeMovers.HY,
        active_power_limits = (min = 0.0, max = 7.0),
        reactive_power_limits = (min = 0.0, max = 7.0),
        ramp_limits = (up = 7.0, down = 7.0),
        time_limits = nothing,
        operation_cost = PSY.StorageManagementCost(
            variable = VariableCost(0.15),
            fixed = 0.0,
            start_up = 0.0,
            shut_down = 0.0,
            energy_shortage_cost = 50.0,
            energy_surplus_cost = 0.0,
        ),
        base_power = 100.0,
        storage_capacity = 50.0,
        inflow = 4.0,
        conversion_factor = 1.0,
        initial_storage = 0.5,
    ),
];

phes5(nodes5) = [
    HydroPumpedStorage(
        name = "HydroPumpedStorage",
        available = true,
        bus = nodes5[3],
        active_power = 0.0,
        reactive_power = 0.0,
        rating = 5.0,
        base_power = 100.0,
        prime_mover = PrimeMovers.HY,
        active_power_limits = (min = 0.0, max = 5.0),
        reactive_power_limits = (min = 0.0, max = 5.0),
        ramp_limits = (up = 10.0 * 0.5, down = 10.0 * 0.5),
        time_limits = nothing,
        operation_cost = TwoPartCost(0.15, 0.0),
        rating_pump = 0.2,
        active_power_limits_pump = (min = 0.0, max = 10.0),
        reactive_power_limits_pump = (min = 0.0, max = 10.0),
        ramp_limits_pump = (up = 10.0 * 0.6, down = 10.0 * 0.6),
        time_limits_pump = nothing,
        storage_capacity = (up = 25.0, down = 25.0), # 50 pu * hr (i.e. 5 GWh)
        inflow = 3.0,
        outflow = 1.0,
        initial_storage = (up = 0.5, down = 0.5),
        storage_target = (up = 0.5, down = 0.75),
        conversion_factor = 1.0,
        pump_efficiency = 1.0,
    ),
];

battery5(nodes5) = [GenericBattery(
    name = "Bat",
    prime_mover = PrimeMovers.BA,
    available = true,
    bus = nodes5[1],
    initial_energy = 2.0,
    state_of_charge_limits = (min = 0.05, max = 4.0),
    rating = 4.0,
    active_power = 4.0,
    input_active_power_limits = (min = 0.0, max = 2.0),
    output_active_power_limits = (min = 0.0, max = 2.0),
    efficiency = (in = 0.80, out = 0.90),
    reactive_power = 0.0,
    reactive_power_limits = (min = -2.0, max = 2.0),
    base_power = 100.0,
)];

batteryems5(nodes5) = [
     PSY.BatteryEMS(;
         name = "Bat2",
         prime_mover = PrimeMovers.BA,
         available = true,
         bus = nodes5[1],
         initial_energy = 5.0,
         state_of_charge_limits = (min = .10, max = 7.0),
         rating = 7.0,
         active_power = 2.0,
         input_active_power_limits = (min = 0.0, max = 2.0),
         output_active_power_limits = (min = 0.0, max = 2.0),
         efficiency = (in = 0.80, out = 0.90),
         reactive_power = 0.0,
         reactive_power_limits = (min = -2.0, max = 2.0),
         base_power = 100.0,
         storage_target=0.2,
         operation_cost = PSY.StorageManagementCost(
            variable = VariableCost(0.0),
            fixed = 0.0,
            start_up = 0.0,
            shut_down = 0.0,
            energy_shortage_cost = 50.0,
            energy_surplus_cost = 40.0,
         ),
     )
 ];

loadbus2_ts_DA = [
    0.792729978
    0.723201574
    0.710952098
    0.677672816
    0.668249175
    0.67166919
    0.687608809
    0.711821241
    0.756320618
    0.7984057
    0.827836527
    0.840362459
    0.84511032
    0.834592803
    0.822949221
    0.816941743
    0.824079963
    0.905735139
    0.989967048
    1
    0.991227765
    0.960842114
    0.921465115
    0.837001437
]

loadbus3_ts_DA = [
    0.831093782
    0.689863228
    0.666058513
    0.627033103
    0.624901388
    0.62858924
    0.650734211
    0.683424321
    0.750876413
    0.828347191
    0.884248576
    0.888523615
    0.87752169
    0.847534405
    0.8227661
    0.803809323
    0.813282799
    0.907575962
    0.98679848
    1
    0.990489904
    0.952520972
    0.906611479
    0.824307054
]

loadbus4_ts_DA = [
    0.871297342
    0.670489749
    0.642812243
    0.630092987
    0.652991383
    0.671971681
    0.716278493
    0.770885833
    0.810075243
    0.85562361
    0.892440566
    0.910660449
    0.922135467
    0.898416969
    0.879816542
    0.896390855
    0.978598576
    0.96523761
    1
    0.969626503
    0.901212601
    0.81894251
    0.771004923
    0.717847996
]

loads5(nodes5) = [
    PowerLoad(
        "Bus2",
        true,
        nodes5[2],
        LoadModels.ConstantPower,
        3.0,
        0.9861,
        100.0,
        3.0,
        0.9861,
    ),
    PowerLoad(
        "Bus3",
        true,
        nodes5[3],
        LoadModels.ConstantPower,
        3.0,
        0.9861,
        100.0,
        3.0,
        0.9861,
    ),
    PowerLoad(
        "Bus4",
        true,
        nodes5[4],
        LoadModels.ConstantPower,
        4.0,
        1.3147,
        100.0,
        4.0,
        1.3147,
    ),
];

interruptible(nodes5) = [InterruptibleLoad(
    "IloadBus4",
    true,
    nodes5[4],
    LoadModels.ConstantPower,
    1.00,
    0.0,
    1.00,
    0.0,
    100.0,
    TwoPartCost(1.50, 24.0),
)]

ORDC_cost = [(9000.0, 0.0), (6000.0, 0.2), (500.0, 0.4), (10.0, 0.6), (0.0, 0.8)]

reserve5(thermal_generators5) = [
    VariableReserve{ReserveUp}(
        "Reserve1",
        true,
        0.6,
        maximum([gen.active_power_limits[:max] for gen in thermal_generators5]) .* 0.001,
    ),
    VariableReserve{ReserveDown}(
        "Reserve2",
        true,
        0.3,
        maximum([gen.active_power_limits[:max] for gen in thermal_generators5]) .* 0.005,
    ),
    VariableReserve{ReserveUp}(
        "Reserve11",
        true,
        0.8,
        maximum([gen.active_power_limits[:max] for gen in thermal_generators5]) .* 0.001,
    ),
    ReserveDemandCurve{ReserveUp}(nothing, "ORDC1", true, 0.6),
    VariableReserveNonSpinning("NonSpinningReserve", true, 0.5, maximum([gen.active_power_limits[:max] for gen in thermal_generators5]) .* 0.001),
]

reserve5_re(renewable_generators5) = [
    VariableReserve{ReserveUp}("Reserve3", true, 30, 100),
    VariableReserve{ReserveDown}("Reserve4", true, 5, 50),
    ReserveDemandCurve{ReserveUp}(nothing, "ORDC1", true, 0.6),
]
reserve5_hy(hydro_generators5) = [
    VariableReserve{ReserveUp}("Reserve5", true, 30, 100),
    VariableReserve{ReserveDown}("Reserve6", true, 5, 50),
    ReserveDemandCurve{ReserveUp}(nothing, "ORDC1", true, 0.6),
]

reserve5_il(interruptible_loads) = [
    VariableReserve{ReserveUp}("Reserve7", true, 30, 100),
    VariableReserve{ReserveDown}("Reserve8", true, 5, 50),
    ReserveDemandCurve{ReserveUp}(nothing, "ORDC1", true, 0.6),
]

ORDC_cost_ts = [
    TimeSeries.TimeArray(DayAhead, repeat([ORDC_cost], 24)),
    TimeSeries.TimeArray(DayAhead + Day(1), repeat([ORDC_cost], 24)),
]

# TODO: add a sensible cost for hybrid devices
hybrid_cost_ts = [
    TimeSeries.TimeArray(DayAhead, repeat([25.0], 24)),
    TimeSeries.TimeArray(DayAhead + Day(1), repeat([25.0], 24)),
]

Reserve_ts = [TimeSeries.TimeArray(DayAhead, rand(24)), TimeSeries.TimeArray(DayAhead + Day(1), rand(24))]

hydro_timeseries_DA = [
    [TimeSeries.TimeArray(DayAhead, hydro_inflow_ts_DA)],
    [TimeSeries.TimeArray(DayAhead + Day(1), ones(24) * 0.1 + hydro_inflow_ts_DA)],
];

storage_target = zeros(24)
storage_target[end] = 0.1
storage_target_DA = [
   [TimeSeries.TimeArray(DayAhead, storage_target)],
   [TimeSeries.TimeArray(DayAhead + Day(1), storage_target)],
];

hydro_budget_DA = [
    [TimeSeries.TimeArray(DayAhead, hydro_inflow_ts_DA * 0.8)],
    [TimeSeries.TimeArray(DayAhead + Day(1), hydro_inflow_ts_DA * 0.8)],
];

RealTime = collect(
    DateTime("1/1/2024 0:00:00", "d/m/y H:M:S"):Minute(5):DateTime(
        "1/1/2024 23:55:00",
        "d/m/y H:M:S",
    ),
)

hydro_timeseries_RT = [
    [TimeSeries.TimeArray(RealTime, repeat(hydro_inflow_ts_DA, inner = 12))],
    [TimeSeries.TimeArray(RealTime + Day(1), ones(288) * 0.1 + repeat(hydro_inflow_ts_DA, inner = 12))],
];

storage_target_RT = [
    [TimeSeries.TimeArray(RealTime, repeat(storage_target, inner = 12))],
    [TimeSeries.TimeArray(RealTime + Day(1), repeat(storage_target, inner = 12))],
];

hydro_budget_RT = [
    [TimeSeries.TimeArray(RealTime, repeat(hydro_inflow_ts_DA  * 0.8, inner = 12))],
    [TimeSeries.TimeArray(RealTime + Day(1), repeat(hydro_inflow_ts_DA  * 0.8, inner = 12))],
];

hybrid_cost_ts_RT = [
    [TimeSeries.TimeArray(RealTime, repeat([25.0], 288))],
    [TimeSeries.TimeArray(RealTime + Day(1), ones(288) * 0.1 + repeat([25.0], 288))],
];

load_timeseries_RT = [
    [
        TimeSeries.TimeArray(RealTime, repeat(loadbus2_ts_DA, inner = 12)),
        TimeSeries.TimeArray(RealTime, repeat(loadbus3_ts_DA, inner = 12)),
        TimeSeries.TimeArray(RealTime, repeat(loadbus4_ts_DA, inner = 12)),
    ],
    [
        TimeSeries.TimeArray(RealTime + Day(1), rand(288) * 0.01 + repeat(loadbus2_ts_DA, inner = 12)),
        TimeSeries.TimeArray(RealTime + Day(1), rand(288) * 0.01 + repeat(loadbus3_ts_DA, inner = 12)),
        TimeSeries.TimeArray(RealTime + Day(1), rand(288) * 0.01 + repeat(loadbus4_ts_DA, inner = 12)),
    ],
]

ren_timeseries_RT = [
    [
        TimeSeries.TimeArray(RealTime, repeat(solar_ts_DA, inner = 12)),
        TimeSeries.TimeArray(RealTime, repeat(wind_ts_DA, inner = 12)),
        TimeSeries.TimeArray(RealTime, repeat(wind_ts_DA, inner = 12)),
    ],
    [
        TimeSeries.TimeArray(RealTime + Day(1), rand(288) * 0.1 + repeat(solar_ts_DA, inner = 12)),
        TimeSeries.TimeArray(RealTime + Day(1), rand(288) * 0.1 + repeat(wind_ts_DA, inner = 12)),
        TimeSeries.TimeArray(RealTime + Day(1), rand(288) * 0.1 + repeat(wind_ts_DA, inner = 12)),
    ],
]

Iload_timeseries_RT = [
    [TimeSeries.TimeArray(RealTime, repeat(loadbus4_ts_DA, inner = 12))],
    [TimeSeries.TimeArray(RealTime + Day(1), rand(288) * 0.1 + repeat(loadbus4_ts_DA, inner = 12))],
]

load_timeseries_DA = [
    [
        TimeSeries.TimeArray(DayAhead, loadbus2_ts_DA),
        TimeSeries.TimeArray(DayAhead, loadbus3_ts_DA),
        TimeSeries.TimeArray(DayAhead, loadbus4_ts_DA),
    ],
    [
        TimeSeries.TimeArray(DayAhead + Day(1), rand(24) * 0.1 + loadbus2_ts_DA),
        TimeSeries.TimeArray(DayAhead + Day(1), rand(24) * 0.1 + loadbus3_ts_DA),
        TimeSeries.TimeArray(DayAhead + Day(1), rand(24) * 0.1 + loadbus4_ts_DA),
    ],
];

ren_timeseries_DA = [
    [
        TimeSeries.TimeArray(DayAhead, solar_ts_DA),
        TimeSeries.TimeArray(DayAhead, wind_ts_DA),
        TimeSeries.TimeArray(DayAhead, wind_ts_DA),
    ],
    [
        TimeSeries.TimeArray(DayAhead + Day(1), rand(24) * 0.1 + solar_ts_DA),
        TimeSeries.TimeArray(DayAhead + Day(1), rand(24) * 0.1 + wind_ts_DA),
        TimeSeries.TimeArray(DayAhead + Day(1), rand(24) * 0.1 + wind_ts_DA),
    ],
];

Iload_timeseries_DA = [
    [TimeSeries.TimeArray(DayAhead, loadbus4_ts_DA)],
    [TimeSeries.TimeArray(DayAhead + Day(1), loadbus4_ts_DA + 0.1 * rand(24))],
]

max_stars_repo_name: NREL-SIIP/PowerSystemBuilder

max_issues_repo_name: NREL-SIIP/PowerSystemBuilder

max_forks_repo_name: NREL-SIIP/PowerSystemBuilder


text:
using LinearAlgebra
using Combinatorics
using BitBasis
using Primes

L = 4

A = with_replacement_combinations([true,false], L) |> collect
A = filter(x->iseven(sum(x)), A)
display(A)


text:
using SFML

window = RenderWindow("Shapes", 800, 600)
event = Event()

texture = Texture(joinpath(dirname(@__FILE__),"wood-texture.jpg"))

triangle = CircleShape()
set_radius(triangle, 100)
set_pointcount(triangle, 3)
set_position(triangle, Vector2f(400, 300))
set_texture(triangle, texture)

convex_shape = ConvexShape()
set_pointcount(convex_shape, 5)
set_point(convex_shape, 0, Vector2f(0, 0))
set_point(convex_shape, 1, Vector2f(150, 10))
set_point(convex_shape, 2, Vector2f(120, 90))
set_point(convex_shape, 3, Vector2f(30, 100))
set_point(convex_shape, 4, Vector2f(0, 50))

set_position(convex_shape, Vector2f(100.0, 400.0))
set_fillcolor(convex_shape, SFML.green)
set_outlinecolor(convex_shape, SFML.red)
set_outline_thickness(convex_shape, 5)

more_points = true

clock = Clock()
restart(clock)

while isopen(window)
    while pollevent(window, event)
        if get_type(event) == EventType.CLOSED
            close(window)
        end
    end

    if as_seconds(get_elapsed_time(clock)) >= 0.75
        triangle_pointcount = get_pointcount(triangle)
        if triangle_pointcount > 8
            more_points = false
        elseif triangle_pointcount <= 3
            more_points = true
        end

        num = more_points ? 1 : -1
        set_pointcount(triangle, get_pointcount(triangle) + num)

        restart(clock)
    end

    clear(window, SFML.black)
    draw(window, triangle)
    draw(window, convex_shape)
    display(window)
end


text:
"""
Implements the `@required_keywords` macro which allows function definitions to have keyword
arguments without a default value. These keywords must always be specified. If they are not,
a `UnassignedKeyword` exception is thrown.
"""
module RequiredKeywords
    import Base: showerror
    export UnassignedKeyword, @required_keywords, @showexp

    """
    Indicated that a funtion has been called without specifying a keyword which has no default.
    """
    struct UnassignedKeyword <:Exception
        msg::String
    end

    Base.showerror(io::IO, e::UnassignedKeyword) = print(io, "Unassigned Keyword:  ", e.msg)

    """
        @required_keywords f(...; x::Int=0, y::Int) = ...

    Allows function definitions to have keyword arguments without a default value. These
    keywords must always be specified. If they are not, a `UnassignedKeyword` exception is thrown.
    """
    macro required_keywords(exp)
    end

    macro showexp(exp)
        show(exp)
    end

end # module


text:
using AutoGrad: Param, @diff, value, grad
using Knet: relu, sigm, tanh, softmax

in_size       = 12
hidden_size   = 20
out_size      = 10

state_size    = 10
memory_size   = 20
location_size = 64

mutable struct Model

    intermediate_state::Param

    mem_attend_value::Param
    mem_attend_func::Param

    hidden_state::Param

    output::Param
    final_state::Param

    out_value::Param
    out_func::Param

    alloc_val::Param
    dealloc_val::Param
    alloc_fn::Param
    dealloc_fn::Param

Model(in_size, hidden_size, out_size, state_size, location_size) = new(

    Param(randn(in_size+state_size, state_size)),

    Param(randn(state_size+location_size, 1)),
    Param(randn(state_size+location_size, 1)),

    Param(randn(in_size+location_size, hidden_size)),

    Param(randn(hidden_size, out_size)),
    Param(randn(hidden_size, state_size)),

    Param(randn(hidden_size, location_size)),
    Param(randn(hidden_size, location_size)),

    Param(randn(location_size*2, 1)),
    Param(randn(location_size*2, 1)),
    Param(randn(location_size*2, 1)),
    Param(randn(location_size*2, 1)),

)
end

(model::Model)(in, state, memory) =
begin

    values, funcs = memory

    intermediate_state = tanh.(hcat(in, state) * model.intermediate_state)
    val_attentions = softmax([(hcat(intermediate_state, location) * model.mem_attend_value)[end] for location in values])
    fn_attentions  = softmax([(hcat(intermediate_state, location) * model.mem_attend_func)[end] for location in funcs])
    attended_val = sum([location .* attention for (location, attention) in zip(values, val_attentions)])
    attended_fn = sum([location .* attention for (location, attention) in zip(funcs, val_attentions)])

    hidden_input = hcat(in, attended_val .* attended_fn)
    hidden_state = tanh.(hidden_input * model.hidden_state)

    output = tanh.(hidden_state * model.output)
    final_state = tanh.(hidden_state * model.final_state)

    out_value = tanh.(hidden_state * model.out_value)
    out_func = tanh.(hidden_state * model.out_func)

    val_attentions = softmax([(hcat(final_state, location) * model.mem_attend_value)[end] for location in values])
    fn_attentions  = softmax([(hcat(final_state, location) * model.mem_attend_func)[end] for location in funcs])

    val_allocs = [sigm.(hcat(out_value, location) * model.alloc_val) for location in values]
    val_deallocs = [sigm.(hcat(out_value, location) * model.dealloc_val) for location in values]
    fn_allocs = [sigm.(hcat(out_value, location) * model.alloc_fn) for location in funcs]
    fn_deallocs = [sigm.(hcat(out_value, location) * model.dealloc_fn) for location in funcs]

    new_values = [(1 .- val_attention) * location + val_attention * (val_dealloc .* location + val_alloc .* out_value) for (location, val_attention, val_dealloc, val_alloc) in zip(values, val_attentions, val_deallocs, val_deallocs)]
    new_funcs = [(1 .- fn_attention) * location + fn_attention * (fn_dealloc .* location + fn_alloc .* out_func) for (location, fn_attention, fn_dealloc, fn_alloc) in zip(funcs, fn_attentions, fn_deallocs, fn_deallocs)]
    new_memory = (new_values, new_funcs)

output, final_state, new_memory
end

in = randn(1, in_size)
state = randn(1, state_size)
memory = [[randn(1, location_size) for _ in 1:memory_size] for _ in 1:2]

model = Model(in_size, hidden_size, out_size, state_size, location_size)

@show state

in, state, memory = model(in, state, memory)

@show state