text:
using Catlab
using Catlab.Doctrines
using Catlab.Graphics
using Catlab.WiringDiagrams
using Catlab.Programs
import Base.Multimedia: display
import Catlab.Graphics: to_graphviz, LeftToRight
println("Done Importing Catlab")
import Base: (==), length, show
using Test
using Petri
using SemanticModels.CategoryTheory
import SemanticModels.CategoryTheory: undecorate, ⊔
using SemanticModels.PetriModels
using SemanticModels.PetriCospans
import SemanticModels.PetriCospans: otimes_ipm, compose_pushout
println("Done Importing SemanticModels")

import Catlab.Doctrines:
  Ob, Hom, dom, codom, compose, ⋅, ∘, id, oplus, otimes, ⊗, ⊕, munit, mzero, braid,
  dagger, dunit, dcounit, mcopy, Δ, delete, ◊, mmerge, ∇, create, □,
  plus, zero, coplus, cozero, meet, top, join, bottom

wd = to_wiring_diagram
draw(d::WiringDiagram) = to_graphviz(add_junctions(d), orientation=LeftToRight, labels=true)
draw(d::HomExpr) = draw(wd(d))

println("Making COVID Model")

@present Disease(FreeEpidemiology) begin
    S::Ob
    E::Ob
    I::Ob
    R::Ob
end

S,E,I,R = generators(Disease)
# draw(spontaneous(E,I)⋅spontaneous(I,R)⊗(spontaneous(E,I)))
sei = compose(exposure(S,I,E), otimes(spontaneous(E,I), id(I)), mmerge(I))
# draw(sei)
seir = sei⋅Δ(I)⋅(id(I)⊗spontaneous(I, R))
# draw(seir)
seir2 = compose(mcopy(S)⊗id(I), id(S)⊗seir)
# draw(seir2)

# draw(death(S))

d = @program Disease (s::S, e::E, i::I) begin
    e1, i1 = exposure{S,I,E}(s,i)
    i2 = spontaneous{E,I}(e1)
    e = [e, e1]
    e_out = spontaneous{E,E}(e)
    i1 = [i1, i2]
    r = spontaneous{I,R}(i1)
    s_out = spontaneous{S,S}(s)
    return s_out, e_out, spontaneous{I,I}(i1)
end
# draw(d)

# draw(d⋅d)

seirdef = to_hom_expr(FreeEpidemiology, d)
try
    add_definition!(Disease, :seir, seirdef)
catch
    println("INFO: definition already added.")
end

# if the disease is fatal, we need to add a death component
seird = @program Disease (s::S, e::E, i::I) begin
    e1, i1 = exposure{S,I,E}(s,i)
    i2 = spontaneous{E,I}(e1)
    e = [e, e1]
    e_out = spontaneous{E,E}(e)
    i1 = [i1, i2]
    r = spontaneous{I,R}(i1)
    s_out = spontaneous{S,S}(s)
    death{I}(i1)
    return s_out, e_out, spontaneous{I,I}(i1)
end
#TODO: This does not get translated correctly, bug?
seirddef = to_hom_expr(FreeEpidemiology, seird)
try
    add_definition!(Disease, :seird, seirddef)
catch
    println("INFO: definition already added.")
end

seirgen = generator(Disease, :seir)
seirdgen = generator(Disease, :seird)

ncities(city,n::Int) = compose([city for i in 1:n]...)
city³ = ncities(seirgen, 3)
# draw(city³)

dcity³ = wd(city³)
dc3 = substitute(dcity³, box_ids(dcity³), [d,d,d])
@show dc3 == ncities(d, 3)
# draw(dc3)

import Base: repeat
repeat(d::WiringDiagram, n::Int) = compose([d for i in 1:n]...)
repeat(d::FreeEpidemiology.Hom, n::Int) = compose([d for i in 1:n]...)

# draw(ncities(seirdgen, 3))
# draw(repeat(seird, 3))

# draw(seirddef)
#
# using TikzPictures
# using Catlab.Graphics.TikZWiringDiagrams
# using Convex
# using SCS
#
# to_tikz(seirddef, labels=true)
X = FinSet(1)
F(ex) = functor((FinSet, FinSetMorph),
                ex,
                generators=Dict(S=>X, I=>X, E=>X, R=>X))

Fseir = compose(exposure(X,X,X),otimes(spontaneous(X,X),id(X)),mmerge(X),mcopy(X),otimes(id(X),spontaneous(X,X)))
# states are [S, I, E, R]
seir_petri = left(Fseir.f).d[1]
f = FinSetMorph(1:4, [1, 2, 3])
g = FinSetMorph(1:4, [1, 2, 3])
Fseir′ = PetriCospan(Cospan(Decorated(f,seir_petri), Decorated(g, seir_petri)))
s = spontaneous(X,X)
Fflow = otimes(s, s, s)
Fcity = Fseir′ ⋅ Fflow
Fcity³ = Fcity⋅Fcity⋅Fcity
@test left(Fcity³.f).d[1].model.S == 1:15
@test left(Fcity³.f).d[1].model.Δ == [([1, 2], [3, 2]),
                                      ([3], [2]),
                                      ([2], [4]),
                                      #outflow 1
                                      ([1], [5]),
                                      ([2], [6]),
                                      ([3], [7]),
                                      # SEIR 2
                                      ([5, 6], [7, 6]),
                                      ([7], [6]),
                                      ([6], [8]),
                                      #outflow 2
                                      ([5], [9]),
                                      ([6], [10]),
                                      ([7], [11]),
                                      # SEIR 3
                                      ([9, 10], [11, 10]),
                                      ([11], [10]),
                                      ([10], [12]),
                                      #outflow 3
                                      ([9], [13]),
                                      ([10], [14]),
                                      ([11], [15])
                                      ]
Fcity₀ = Fseir′ ⋅ Fflow
Fcity₁ = Fflow ⋅ Fseir′ ⋅ Fflow
Fcityₑ = Fflow ⋅ Fseir′
Fcity³ = Fcity₀⋅Fcity₁⋅Fcityₑ

@test left(Fcity³.f).d[1].model.S == 1:18
@test left(Fcity³.f).d[1].model.Δ == [
    # SEIR 1
     ([1, 2], [3, 2]),
     ([3], [2]),
     ([2], [4]),
     # outflow 1→2
     ([1], [5]),
     ([2], [6]),
     ([3], [7]),
     # inflow 1→2
     ([5], [8]),
     ([6], [9]),
     ([7], [10]),
    # SEIR 2
     ([8, 9], [10, 9]),
     ([10], [9]),
     ([9], [11]),
     # outflow 2→3
     ([8], [12]),
     ([9], [13]),
     ([10], [14]),
     # inflow 2→3
     ([12], [15]),
     ([13], [16]),
     ([14], [17]),
    # SEIR 3
     ([15, 16], [17, 16]),
     ([17], [16]),
     ([16], [18])
]

Fcity₀ = Fseir′ ⋅ Fflow
Fcity₁ = Fseir′ ⋅ Fflow
Fcityₑ = Fseir′
Fcity³ = Fcity₀⋅Fcity₁⋅Fcityₑ

@test left(Fcity³.f).d[1].model.S == 1:12
@test left(Fcity³.f).d[1].model.Δ == [
    # City 1 SEIR
    ([1, 2], [3, 2]),
    ([3], [2]), # E→I
    ([2], [4]), # I→R
    # outflow 1→2
    ([1], [5]), # S→S′
    ([2], [6]), # I→I′
    ([3], [7]), # E→E′
    # City 2 SEIR
    ([5, 6], [7, 6]),
    ([7], [6]), # E→I
    ([6], [8]), # I→R
    # outflow 2→3
    ([5], [9]), # S→S′
    ([6], [10]),# I→I′
    ([7], [11]),# E→E′
    # City 3 SEIR
    ([9, 10], [11, 10]),
    ([11], [10]), # E→I
    ([10], [12]) # I→R
    ]

Pseird = PetriModel(
          Petri.Model(1:5,[
            ([1,2],[3,2]), # exposure
            ([3],[2]),     # onset
            ([2],[4]),     # recovery
            ([2],[5]),     # death
            ], missing, missing))
inputs = FinSetMorph(1:5, [1,2,3])
outputs = FinSetMorph(1:5, [1,2,3])
Fcityd = PetriCospan(Cospan(Decorated(inputs, Pseird),
                            Decorated(outputs, Pseird)))
Fcity₀ = Fcityd ⋅ Fflow
Fcity₁ = Fcityd ⋅ Fflow
Fcityₑ = Fcityd
Fcity³ = Fcity₀⋅Fcity₁⋅Fcityₑ

left(Fcity³.f).d[1].model.S == 1:15
length(left(Fcity³.f).d[1].model.Δ) == 18


text:
module Library

using BioSequences
using ..RNASplicing: OpaqueTransElement, Interaction, Spliceosome, only
# The core splicing machinery can be expressed in this framework

# # Core Machinery

"5’ splice site - 9-nucleotide motif that generally conforms to the sequence"
donorsite = biore"(C|A)AG|GU(A|G)AGU"rna
# donorsite(seq) = occursin(biore"(C|A)AG|GU(A|G)AGU"rna, seq) ? 1.0 : 0.0
u1 = OpaqueTransElement(:u1)
u2 = OpaqueTransElement(:u2)
u2af = OpaqueTransElement(:u2af)
e1 = Interaction(donorsite, only(u1), 0, 10.0)
coremachinery = [e1]

acceptsite = biore"A(C|U)AG"rna

# # SREs

# `seq` have a sequence of at least three Gs?
grun = biore"GGG+"rna

# Major splicesosome
majorsos = Spliceosome(Dict(u1 => 1.0))
end

ext: jl


text:
using BSR
using Images, ImageView

function main()
    # load image from data directory
    pkgdir = joinpath(Pkg.dir(), "BSR")
    datadir = joinpath(pkgdir, "deps/downloads/BSR/grouping/data")
    imagefile = joinpath(datadir, "101087_small.jpg")
    img = permutedims(data(reinterpret(Float64, float64(imread(imagefile)))), [3, 2, 1])

    # display image
    view(img, pixelspacing = [1, 1])
    println("original image (press enter to continue)")
    readline()

    # compute gPb and display result
    println("computing gpb (this may take some time...)")
    gpbresponse = gpb(img)
    nrows = 2
    ncols = 4
    grid = canvasgrid(nrows, ncols)
    ops = [:pixelspacing => [1, 1]]
    for d = 1:size(gpbresponse, 3)
        row = div(d - 1, ncols) + 1
        col = mod1(d, ncols)
        view(grid[row, col], gpbresponse[:, :, d]; ops...)
    end
    println("gpb response (press enter to continue)")
    readline()
end

main()


text:
using NeuralQuantum
using NeuralQuantum: state_uview
using Test
using StatsBase: Histogram, fit
using HypothesisTests

function pdf(net, hilb)
    cnet = cached(net)
    ψ = ComplexF64.(Vector(cnet, hilb, false))
    probs = (abs2.(ψ)./sum(abs2.(ψ)))
    return probs
end

function sample_freqs(sampler, net, hilb)
    sampl = NeuralQuantum.SimpleIterativeSampler(net, sampler, hilb, batch_sz=16)
    NeuralQuantum.sample!(sampl)
    # extract indices
    ids = index(hilb, sampl.samples)

    hist = fit(Histogram, vec(ids), 1:spacedimension(hilb)+1)
    nhist = normalize(hist, mode=:pdf)

    return hist.weights
end

function test_sampler(sampler, net, hilb)
    exact_dist    = pdf(net, hilb)
    sampled_freqs = sample_freqs(sampler, net, hilb)
    tst = ChisqTest(sampled_freqs, exact_dist)
    @test pvalue(tst) >= 0.01
end

N = 4
n_samples = 1000
SEED = 123

lat  = HyperCube([N], periodic=true)
hilb = HomogeneousSpin(lat, 1//2)
H    = quantum_ising_hamiltonian(ComplexF32, lat, hilb, g=-1.0, V=2.0)

net = RBM(N, 1)
init_random_pars!(net, sigma=0.2^2)

n_samples = max(20000, 50*spacedimension(hilb))

@testset "Hamiltonian" begin
    @testset "ExactSampler" begin
        sampler = ExactSampler(n_samples, seed=SEED)
        test_sampler(sampler, net, hilb)
    end

    @testset "LocalRule" begin
        sampler = MetropolisSampler(LocalRule(), N*n_samples, 1, seed=SEED)
        test_sampler(sampler, net, hilb)

        sampler = MetropolisSampler(LocalRule(), n_samples, N+1, seed=SEED)
        test_sampler(sampler, net, hilb)
    end

    @testset "OperatorRule" begin
        sampler = MetropolisSampler(OperatorRule(H), n_samples, N+1, seed=SEED)
        test_sampler(sampler, net, hilb)
    end
end

@testset "Liouvillian" begin
    hilb2 = SuperOpSpace(hilb)
    net = NDM(N, 1, 1, af_logcosh)
    n_samples = max(N*10000, 40*spacedimension(hilb2))

    @testset "ExactSampler" begin
        sampler = ExactSampler(n_samples)
        test_sampler(sampler, net, hilb2)
    end

    @testset "LocalRule" begin
        sampler = MetropolisSampler(LocalRule(), n_samples, N+1, seed=SEED)
        test_sampler(sampler, net, hilb2)
    end

    @testset "NagyRule" begin
        sampler = MetropolisSampler(NagyRule(H), n_samples, N+1, seed=SEED)
        test_sampler(sampler, net, hilb2)
    end
end


text:
using Pkg
Pkg.activate(".")
println("Loading SpineBasedRecordLinkage.jl")
using SpineBasedRecordLinkage
if length(ARGS) != 3
    error("The number of arguments is $(length(ARGS)), should be 3.")
end
directory1 = ARGS[1]
directory2 = ARGS[2]
outfile    = ARGS[3]
compare_linkage_runs(directory1, directory2, outfile)


text:
"""

```
mddf_linkedcells_parallel(trajectory::Trajectory, options::Options,samples::Samples, RNG, mddf_compute!)
```

Computes the MDDF using linked cells, in parallel.

"""
function mddf_linkedcells_parallel(trajectory::Trajectory, options::Options,
                                   samples::Samples, RNG, mddf_compute!)

  # Simplify code by assigning some shortened names
  solute = trajectory.solute
  solvent = trajectory.solvent
  x_solute = trajectory.x_solute
  x_solvent = trajectory.x_solvent

  # Number of threads
  if options.nthreads < 0
    nspawn = Threads.nthreads() - 1
  else
    nspawn = options.nthreads - 1
  end
  if nspawn == 0
    error(" Parallel version must be executed only with more than 1 thread. ")
  end

  # Initializing the structure that carries the result per thread
  R0 = Result(trajectory,options)
  R = [ Result(trajectory,options,irefatom=R0.irefatom) for i in 1:nspawn ]

  # Safe passing of frame counter to the threads
  tframe = zeros(Int,nspawn)

  # Create data structures required for multithreading
  framedata = Vector{FrameData}(undef,nspawn)
  for ispawn in 1:nspawn
    framedata[ispawn] = FrameData(deepcopy(trajectory),R[1])
  end

  # Tasks
  t = Vector{Task}(undef,nspawn)
  free = ones(Bool,nspawn)

  # Skip initial frames if desired
  iframe = 0 # Counter for all frames of the input file
  while iframe < options.firstframe-1
    nextframe!(trajectory)
    iframe += 1
  end

  # Print some information about this run
  options.silent || title(R[1],solute,solvent,nspawn)

  ndone = 0
  iframe_read = 0 # Counter for the frames that are actually being considered
  if ! options.silent
    progress = Progress(R[1].nframes_read,1)
  end
  while ndone < R[1].nframes_read

    # Launch for each free thread the computation of a frame
    while iframe_read < R[1].nframes_read && count(free) > 0
      ifree = findfirst(x->x==true,free)
      # Skip frames that are not going to be considered
      while (iframe+1)%options.stride != 0
        nextframe!(trajectory)
        iframe += 1
      end
      # Actually read the data of the frame that will be considered
      nextframe!(trajectory)
      iframe += 1
      iframe_read += 1
      tframe[ifree] = iframe
      @. framedata[ifree].trajectory.x_solute = trajectory.x_solute
      @. framedata[ifree].trajectory.x_solvent = trajectory.x_solvent
      @. framedata[ifree].trajectory.sides = trajectory.sides
      # Spawn the calculations for this frame
      t[ifree] = ThreadPools.@tspawnat ifree+1 mddf_compute!(tframe[ifree],framedata[ifree],options,RNG,R[ifree])
      free[ifree] = false
    end

    # Wait a little bit before checking
    sleep(options.sleep)

    # Check thread status
    for ispawn in 1:nspawn
      if ! free[ispawn]
        if istaskfailed(t[ispawn])
          error(" Computation of MDDF failed in thread: $ispawn", fetch(t[ispawn]))
        end
        if istaskdone(t[ispawn])
          ndone += 1
          free[ispawn] = true
          options.silent || next!(progress)
          if options.GC && (Sys.free_memory() / Sys.total_memory() < options.GC_threshold)
            GC.gc() # why we need this anyway??? There should not be so much garbage.
          end
        end
      end
    end

  end
  closetraj(trajectory)

  # Sum up the results of all threads into the data of thread one (R1<-R1+R2)
  for ispawn in 2:nspawn
    sum!(R[1],R[ispawn])
  end

  # Setup the final data structure with final values averaged over the number of frames,
  # sampling, etc, and computes final distributions and integrals
  finalresults!(R[1],options,trajectory,samples)
  options.silent || println(bars)

  return R[1]

end


text:
TSAD_UNIVARIATE = ["ACSF1", "Adiac", "AllGestureWiimoteX", "AllGestureWiimoteY", "AllGestureWiimoteZ", "ArrowHead", "BME", "Beef", "BeetleFly", "BirdChicken", "CBF", "Car", "Chinatown", "ChlorineConcentration", "CinCECGTorso", "Coffee", "Computers", "CricketX", "CricketY", "CricketZ", "Crop", "DiatomSizeReduction", "DistalPhalanxOutlineAgeGroup", "DistalPhalanxOutlineCorrect", "DistalPhalanxTW", "DodgerLoopDay", "DodgerLoopGame", "DodgerLoopWeekend", "ECG200", "ECG5000", "ECGFiveDays", "EOGHorizontalSignal", "EOGVerticalSignal", "Earthquakes", "ElectricDevices", "EthanolLevel", "FaceAll", "FaceFour", "FacesUCR", "FiftyWords", "Fish", "FordA", "FordB", "FreezerRegularTrain", "FreezerSmallTrain", "Fungi", "GestureMidAirD1", "GestureMidAirD2", "GestureMidAirD3", "GesturePebbleZ1", "GesturePebbleZ2", "GunPoint", "GunPointAgeSpan", "GunPointMaleVersusFemale", "GunPointOldVersusYoung", "Ham", "HandOutlines", "Haptics", "Herring", "HouseTwenty", "InlineSkate", "InsectEPGRegularTrain", "InsectEPGSmallTrain", "InsectWingbeatSound", "ItalyPowerDemand", "LargeKitchenAppliances", "Lightning2", "Lightning7", "Mallat", "Meat", "MelbournePedestrian", "MiddlePhalanxOutlineAgeGroup", "MiddlePhalanxOutlineCorrect", "MiddlePhalanxTW", "MixedShapesRegularTrain", "MixedShapesSmallTrain", "MoteStrain", "NonInvasiveFetalECGThorax1", "NonInvasiveFetalECGThorax2", "OSULeaf", "OliveOil", "PLAID", "PhalangesOutlinesCorrect", "Phoneme", "PickupGestureWiimoteZ", "PigAirwayPressure", "PigArtPressure", "PigCVP", "Plane", "PowerCons", "ProximalPhalanxOutlineAgeGroup", "ProximalPhalanxOutlineCorrect", "ProximalPhalanxTW", "RefrigerationDevices", "Rock", "ScreenType", "SemgHandGenderCh2", "SemgHandMovementCh2", "SemgHandSubjectCh2", "ShakeGestureWiimoteZ", "ShapeletSim", "ShapesAll", "SmallKitchenAppliances", "SmoothSubspace", "SonyAIBORobotSurface1", "SonyAIBORobotSurface2", "StarLightCurves", "Strawberry", "SwedishLeaf", "Symbols", "SyntheticControl", "ToeSegmentation1", "ToeSegmentation2", "Trace", "TwoLeadECG", "TwoPatterns", "UMD", "UWaveGestureLibraryAll", "UWaveGestureLibraryX", "UWaveGestureLibraryY", "UWaveGestureLibraryZ", "Wafer", "Wine", "WordSynonyms", "Worms", "WormsTwoClass", "Yoga"]

# the default anomaly class is defined as: the lexically first class with the least occurences
# Example: Dict([("class2", 10), ("class1", 10), ("class3": 20))]) -> "class1"
TSAD_UNIVARIATE_DEFAULTS = Dict("ECG5000" => "5", "FordA" => "1", "DodgerLoopGame" => "2", "Trace" => "1", "DodgerLoopDay" => "3", "CricketZ" => "1", "GunPointOldVersusYoung" => "1", "HandOutlines" => "0", "Ham" => "1", "CricketY" => "1", "DiatomSizeReduction" => "1", "SmoothSubspace" => "1", "TwoPatterns" => "4", "Phoneme" => "39", "GunPointMaleVersusFemale" => "2", "SmallKitchenAppliances" => "1", "UWaveGestureLibraryZ" => "1", "Strawberry" => "1", "Mallat" => "1", "Fish" => "5", "GestureMidAirD3" => "1", "ShakeGestureWiimoteZ" => "1", "FreezerRegularTrain" => "1", "FaceFour" => "1", "WordSynonyms" => "25", "TwoLeadECG" => "1", "Meat" => "1", "BirdChicken" => "1", "EOGVerticalSignal" => "12", "ProximalPhalanxOutlineAgeGroup" => "1", "SwedishLeaf" => "1", "UMD" => "1", "MiddlePhalanxTW" => "7", "EOGHorizontalSignal" => "12", "GestureMidAirD1" => "1", "GesturePebbleZ2" => "2", "NonInvasiveFetalECGThorax2" => "11", "ScreenType" => "1", "BeetleFly" => "1", "SyntheticControl" => "1", "WormsTwoClass" => "1", "Yoga" => "1", "AllGestureWiimoteZ" => "1", "DistalPhalanxOutlineCorrect" => "0", "Adiac" => "1", "FacesUCR" => "11", "BME" => "1", "InsectEPGSmallTrain" => "3", "Coffee" => "1", "StarLightCurves" => "1", "SemgHandGenderCh2" => "2", "PigCVP" => "1", "Wine" => "2", "PickupGestureWiimoteZ" => "1", "InlineSkate" => "1", "MixedShapesRegularTrain" => "5", "MixedShapesSmallTrain" => "5", "Wafer" => "-1", "UWaveGestureLibraryX" => "1", "AllGestureWiimoteY" => "1", "CBF" => "1", "Crop" => "1", "Earthquakes" => "1", "GunPointAgeSpan" => "2", "SonyAIBORobotSurface2" => "1", "Worms" => "5", "ECGFiveDays" => "1", "PhalangesOutlinesCorrect" => "0", "Haptics" => "1", "InsectWingbeatSound" => "1", "ECG200" => "-1", "OliveOil" => "3", "ShapesAll" => "1", "Chinatown" => "1", "Beef" => "2", "SemgHandMovementCh2" => "1", "FreezerSmallTrain" => "1", "InsectEPGRegularTrain" => "3", "Rock" => "3", "GunPoint" => "1", "Symbols" => "2", "Lightning7" => "2", "PowerCons" => "1", "Computers" => "1", "GestureMidAirD2" => "1", "HouseTwenty" => "2", "ArrowHead" => "1", "OSULeaf" => "6", "ACSF1" => "0", "CricketX" => "1", "Plane" => "1", "FordB" => "1", "FaceAll" => "11", "ToeSegmentation2" => "1", "CinCECGTorso" => "1", "MiddlePhalanxOutlineAgeGroup" => "1", "ProximalPhalanxOutlineCorrect" => "0", "AllGestureWiimoteX" => "1", "SemgHandSubjectCh2" => "1", "UWaveGestureLibraryAll" => "1", "Fungi" => "1", "Herring" => "2", "DodgerLoopWeekend" => "2", "PigAirwayPressure" => "1", "ProximalPhalanxTW" => "3", "PLAID" => "10", "ChlorineConcentration" => "1", "PigArtPressure" => "1", "FiftyWords" => "49", "EthanolLevel" => "2", "MiddlePhalanxOutlineCorrect" => "0", "ElectricDevices" => "6", "UWaveGestureLibraryY" => "1", "SonyAIBORobotSurface1" => "2", "DistalPhalanxOutlineAgeGroup" => "1", "GesturePebbleZ1" => "2", "ToeSegmentation1" => "1", "DistalPhalanxTW" => "7", "LargeKitchenAppliances" => "1", "NonInvasiveFetalECGThorax1" => "11", "Car" => "1", "RefrigerationDevices" => "1", "ItalyPowerDemand" => "1", "MelbournePedestrian" => "8", "ShapeletSim" => "0", "MoteStrain" => "2", "Lightning2" => "-1")

# a dictionary mapping the dataset names to their respective class counts
TSAD_UNIVARIATE_CLASSES = Dict{String, Dict}("ECG5000" => Dict("4" => 194, "1" => 2919, "5" => 24, "2" => 1767, "3" => 96), "FordA" => Dict("1" => 2394, "-1" => 2527), "DodgerLoopGame" => Dict("1" => 82, "2" => 76), "Trace" => Dict("4" => 50, "1" => 50, "2" => 50, "3" => 50), "DodgerLoopDay" => Dict("4" => 22, "1" => 24, "5" => 24, "2" => 21, "6" => 25, "7" => 22, "3" => 20), "CricketZ" => Dict("4" => 65, "1" => 65, "12" => 65, "2" => 65, "6" => 65, "7" => 65, "11" => 65, "8" => 65, "10" => 65, "5" => 65, "9" => 65, "3" => 65), "GunPointOldVersusYoung" => Dict("1" => 215, "2" => 236), "HandOutlines" => Dict("1" => 874, "0" => 494), "Ham" => Dict("1" => 103, "2" => 111), "CricketY" => Dict("4" => 65, "1" => 65, "12" => 65, "2" => 65, "6" => 65, "7" => 65, "11" => 65, "8" => 65, "10" => 65, "5" => 65,
"9" => 65, "3" => 65), "DiatomSizeReduction" => Dict("4" => 91, "1" => 34, "2" => 98, "3" => 99), "SmoothSubspace" => Dict("1" => 100, "2" => 100, "3" => 100), "TwoPatterns" => Dict("4" => 1201, "1" => 1306, "2" => 1248, "3" => 1245), "Phoneme" => Dict("24" => 22, "29" => 123, "1" => 42, "12" => 68, "32" => 9, "20" => 100, "4" => 26, "2" => 54, "25" => 32, "6" => 33, "23" => 149, "22" => 71, "11" => 66, "35" => 28, "13" => 41, "15" => 26, "27" => 64, "31" => 145, "5" => 9, "33" => 6, "28" => 113, "16" => 17, "14" => 38, "39" => 2, "21" => 125, "38" => 24, "7" => 47, "34" => 23, "8" => 12, "36" => 20, "26" => 4, "17" => 109, "10" => 3, "19" => 18, "37" => 14, "18" => 83, "9" => 78, "30" => 28, "3" => 238), "GunPointMaleVersusFemale" => Dict("1" => 237, "2" => 214), "SmallKitchenAppliances" => Dict("1" => 250, "2" => 250, "3" => 250), "UWaveGestureLibraryZ" => Dict("8" => 560, "4" => 560, "1" => 559, "5" => 559, "2" => 560, "6" => 559, "7" => 559, "3" => 560), "Strawberry" => Dict("1" => 351, "2" => 632), "Mallat" => Dict("8" => 300, "4" => 300, "1" => 300, "5" => 300, "2" => 300, "6" => 300, "7" => 300, "3" => 300), "Fish" => Dict("4" => 50, "1" => 50, "5" => 49, "2" => 50, "6" => 50, "7" => 49, "3" => 50), "GestureMidAirD3" => Dict("24" => 13, "4" => 13, "1" => 13, "12" => 13, "20" => 13, "2" => 13, "25" => 13, "6" => 13, "23" => 13, "22" => 13, "11" => 13, "13" => 13, "15" => 13, "5" => 13, "16" => 13, "14" => 13, "21" => 13, "7" => 13, "8" => 13, "26" => 13, "17" => 13, "10" => 13, "19" => 13, "18" => 13, "9" => 13, "3" => 13), "ShakeGestureWiimoteZ" => Dict("8" => 10, "4" => 10, "1" => 10, "10" => 10, "2" => 10, "5" => 10, "6" => 10, "7" => 10, "9" => 10, "3" => 10), "FreezerRegularTrain" => Dict("1" => 1498, "2" => 1500), "FaceFour" => Dict("4" => 27, "1" => 22, "2" => 34, "3" => 29), "WordSynonyms" => Dict("24" => 20, "4" => 115, "1" => 20, "12" => 44, "20" => 26, "2" => 200, "25" => 12, "6" => 72, "23" => 14, "22" => 23, "11" => 16, "13" => 16, "15" => 16, "5" => 18, "16" => 30, "14" => 39, "21" => 14, "7" => 18, "8" => 56, "17" => 14, "10" => 46, "19" => 14, "18" => 27, "9" => 16, "3" => 19), "TwoLeadECG" => Dict("1" => 581, "2" => 581), "Meat" => Dict("1" => 40, "2" => 40, "3" => 40), "BirdChicken" => Dict("1" => 20, "2" => 20), "EOGVerticalSignal" => Dict("4" => 60, "1" => 60, "12" => 59, "2" => 60, "6" => 61, "7" => 61, "11" => 61, "8" => 60, "10" => 60, "5" => 61, "9" => 61, "3" => 60), "ProximalPhalanxOutlineAgeGroup" => Dict("1" => 89, "2" => 226, "3" =>
288), "SwedishLeaf" => Dict("3" => 75, "4" => 75, "1" => 75, "12" => 75, "2" => 75, "6" => 75, "7" => 75, "11" => 75, "13" => 75, "8" => 75, "10" => 75, "15" => 75, "5" => 75, "9" => 75, "14" => 75), "UMD" => Dict("1" => 60, "2" => 60, "3" => 60), "MiddlePhalanxTW" => Dict("8" => 200, "4" => 128, "5" => 66, "6" => 41, "7" => 34, "3" => 82), "EOGHorizontalSignal" => Dict("4" => 60, "1" => 60, "12" => 59, "2" => 60, "6" => 61, "7" => 61, "11" => 61, "8" => 60, "10" =>
60, "5" => 61, "9" => 61, "3" => 60), "GestureMidAirD1" => Dict("24" => 13, "4" => 13, "1" => 13, "12" => 13, "20" => 13, "2" => 13, "25" => 13, "6" => 13, "23" => 13, "22" => 13, "11" => 13, "13" => 13, "15" => 13, "5" => 13, "16" => 13, "14" => 13, "21" => 13, "7" => 13, "8" => 13, "26" => 13, "17" => 13, "10"
=> 13, "19" => 13, "18" => 13, "9" => 13, "3" => 13), "GesturePebbleZ2" => Dict("4" => 53, "1" => 49, "5" => 52, "2" => 45, "6" => 51, "3" => 54), "NonInvasiveFetalECGThorax2" => Dict("24" => 87, "29" => 96, "1" => 84, "12" => 90, "32" => 90, "20" => 90, "4" => 85, "2" => 91, "25" => 90, "6" => 94, "23" => 87, "22" => 86, "11" => 75, "35" => 94, "13" => 97, "41" => 93, "15" => 98, "27" => 89, "31" => 94, "42" => 93, "33" => 85, "28" => 85, "16" => 95, "5" => 96, "14" => 89, "39" => 90, "40" => 91, "21" => 89, "38" => 91, "7" => 93, "34" => 94, "8" => 91, "36" => 91, "26" => 78, "17" => 85, "10" => 92, "19" => 88, "37" => 87, "18" => 91, "9" => 88, "30" => 86, "3" => 87), "ScreenType" => Dict("1" => 250, "2" => 250, "3" => 250), "BeetleFly" => Dict("1" => 20, "2" => 20), "SyntheticControl" => Dict("4" => 100, "1" => 100, "5" => 100, "2" => 100, "6" => 100, "3" => 100), "WormsTwoClass" => Dict("1" => 109, "2" => 149), "Yoga" => Dict("1" => 1530, "2" => 1770), "AllGestureWiimoteZ" => Dict("8" => 100, "4" => 100, "1" => 100, "10" => 100, "2" => 100, "5" => 100, "6" => 100, "7" => 100, "9" => 100, "3" => 100), "DistalPhalanxOutlineCorrect" => Dict("1" => 537, "0" => 337), "Adiac" => Dict("24" => 23, "29" => 20, "1" => 20, "12" => 26, "32" => 20, "20" => 21, "4" => 20, "2" => 23, "25" => 29, "6" => 20, "23" => 20, "22" => 20, "11" => 20, "35" => 20, "13" => 20, "15" => 21, "27" => 20, "31" => 20, "5" => 20, "33" => 20, "28" => 20, "16" => 20, "14" => 20, "21" => 20, "7" => 20, "34" => 20, "8" => 25, "36" => 25, "26" => 20, "17" => 20, "10" => 20, "19" => 26, "37" => 20, "18" => 22, "9" => 20, "30" => 20, "3" => 20), "FacesUCR" => Dict("3" => 176, "4" => 171, "1" => 112, "12" => 106, "2" => 178, "6" => 187, "7" => 149, "11" => 48, "13" => 327, "8" => 273, "10" => 135, "5" => 176, "9" => 140, "14" => 72), "BME" => Dict("1" => 60, "2" => 60, "3" => 60), "InsectEPGSmallTrain" => Dict("1" => 95, "2" => 126, "3" => 45), "Coffee" => Dict("1" => 27, "0" => 29), "StarLightCurves" => Dict("1" => 1329, "2" => 2580, "3" =>
5327), "SemgHandGenderCh2" => Dict("1" => 540, "2" => 360), "PigCVP" => Dict("32" => 6, "29" => 6, "1" => 6, "12" => 6, "9" => 6, "20" => 6, "2" => 6, "6" => 6, "41" => 6, "51" => 6, "27" => 6, "42" => 6, "33" => 6, "28" => 6, "50" => 6, "52" => 6, "14" => 6, "7" => 6, "8" => 6, "26" => 6, "10" => 6, "19" => 6, "24" => 6, "4" => 6, "25" => 6, "23" => 6, "22" => 6, "11" => 6, "35" => 6, "13" => 6, "49" => 6, "15" => 6, "31" => 6, "5" => 6, "43" => 6, "16" => 6, "39" => 6, "40" => 6, "21" => 6, "38" => 6, "46" => 6, "34" => 6, "44" => 6, "45" => 6, "36" => 6, "17" => 6, "37" => 6, "47" => 6, "18" => 6, "48" => 6, "30" => 6, "3" => 6), "Wine" => Dict("1" => 57, "2" => 54), "PickupGestureWiimoteZ" => Dict("8" => 10, "4" => 10, "1" => 10, "10" => 10, "2" => 10, "5" => 10, "6" =>
10, "7" => 10, "9" => 10, "3" => 10), "InlineSkate" => Dict("4" => 108, "1" => 62, "5" => 117, "2" => 100, "6" => 98, "7" => 62, "3" => 103), "MixedShapesRegularTrain" => Dict("4" => 681, "1" => 558, "5" => 414, "2" => 754, "3" => 518), "MixedShapesSmallTrain" => Dict("4" => 601, "1" => 478, "5" => 334, "2" => 674, "3" => 438), "Wafer" => Dict("1" => 6402, "-1" => 762), "UWaveGestureLibraryX" => Dict("8" => 560, "4" => 560, "1" => 559, "5" => 559, "2" => 560, "6" =>
559, "7" => 559, "3" => 560), "AllGestureWiimoteY" => Dict("8" => 100, "4" => 100, "1" => 100, "10" => 100, "2" => 100, "5" => 100, "6" => 100, "7" => 100, "9" => 100, "3" => 100), "CBF" => Dict("1" => 310, "2" => 310, "3" => 310), "Crop" => Dict("24" => 1000, "4" => 1000, "1" => 1000, "12" => 1000, "20" => 1000, "2" => 1000, "6" => 1000, "23" => 1000, "22" => 1000, "11" => 1000, "13" => 1000, "15" => 1000, "5" => 1000, "16" => 1000, "14" => 1000, "21" => 1000, "7" => 1000, "8" => 1000, "17" => 1000, "10" => 1000, "19" => 1000, "18" => 1000, "9" => 1000, "3" => 1000), "Earthquakes" => Dict("1" => 93, "0" => 368), "GunPointAgeSpan" => Dict("1" => 228, "2" => 223), "SonyAIBORobotSurface2" => Dict("1" => 376, "2" => 604), "Worms" => Dict("4" => 45, "1" => 109, "5" => 25, "2" => 44, "3" => 35), "ECGFiveDays" => Dict("1" => 442, "2" => 442), "PhalangesOutlinesCorrect" => Dict("1" => 1696, "0" => 960), "Haptics" => Dict("4" => 100, "1" => 78, "5" => 100, "2" => 92, "3" => 93), "InsectWingbeatSound" => Dict("4" => 200, "1" => 200, "2" => 200, "6" => 200, "7" => 200, "11" => 200, "8" => 200, "10" => 200, "5" => 200, "9" => 200, "3" => 200), "ECG200" => Dict("1" => 133, "-1" => 67), "OliveOil" => Dict("4" => 25, "1" => 10, "2" => 17, "3" => 8),
"ShapesAll" => Dict("32" => 20, "29" => 20, "1" => 20, "12" => 20, "54" => 20, "20" => 20, "9" => 20, "2" => 20, "6" => 20, "41" => 20, "51" => 20, "53" => 20, "27" => 20, "42" => 20, "33" => 20, "28" => 20, "50" => 20, "52" => 20, "14" => 20, "7" => 20, "57" => 20, "8" => 20, "26" => 20, "10" => 20, "19" => 20,
"60" => 20, "24" => 20, "4" => 20, "58" => 20, "25" => 20, "23" => 20, "22" => 20, "11" => 20, "35" => 20, "13" => 20, "49" => 20, "59" => 20, "15" => 20, "31" => 20, "5" => 20, "43" => 20, "16" => 20, "39" => 20, "40" => 20, "21" => 20, "38" => 20, "46" => 20, "34" => 20, "44" => 20, "45" => 20, "56" => 20, "55" => 20, "36" => 20, "17" => 20, "37" => 20, "47" => 20, "18" => 20, "48" => 20, "30" => 20, "3" => 20), "Chinatown" => Dict("1" => 104, "2" => 259), "Beef" => Dict("4" => 12, "5" => 12, "2" => 2, "3" => 12), "SemgHandMovementCh2" => Dict("4" => 150, "1" => 150, "5" => 150, "2" => 150, "6" => 150, "3" => 150), "FreezerSmallTrain" => Dict("1" => 1437, "2" => 1439), "InsectEPGRegularTrain" => Dict("1" => 111, "2" => 148, "3" => 52), "Rock" => Dict("4" => 19, "1" => 26,
"2" => 14, "3" => 11), "GunPoint" => Dict("1" => 100, "2" => 100), "Symbols" => Dict("4" => 181, "1" => 181, "5" => 162, "2" => 162, "6" => 167, "3" => 167), "Lightning7" => Dict("4" => 15, "1" => 17, "5" => 38, "0" => 18, "2" => 14, "6" => 22, "3" => 19), "PowerCons" => Dict("1" => 180, "2" => 180), "Computers"
=> Dict("1" => 250, "2" => 250), "GestureMidAirD2" => Dict("24" => 13, "4" => 13, "1" => 13, "12" => 13, "20" => 13, "2" => 13, "25" => 13, "6" => 13, "23" => 13, "22" => 13, "11" => 13, "13" => 13, "15" => 13, "5" => 13, "16" => 13, "14" => 13, "21" => 13, "7" => 13, "8" => 13, "26" => 13, "17" => 13, "10" => 13, "19" => 13, "18" => 13, "9" => 13, "3" => 13), "HouseTwenty" => Dict("1" => 87, "2" => 70), "ArrowHead" => Dict("1" => 65, "0" => 81, "2" => 65), "OSULeaf" => Dict("4" => 97, "1" => 66, "5" => 81, "2" => 84, "6" => 37, "3" => 75), "ACSF1" => Dict("8" => 20, "4" => 20, "1" => 20, "5" => 20, "0" => 20, "2" => 20, "6" => 20, "7" => 20, "9" => 20, "3" => 20), "CricketX" => Dict("4" => 65, "1" => 65, "12" => 65, "2" => 65, "6" => 65, "7" => 65, "11" => 65, "8" => 65, "10" => 65, "5" => 65, "9" => 65, "3" => 65), "Plane" => Dict("4" => 30, "1" => 30, "5" => 30, "2" => 30, "6" => 30, "7" => 30, "3" => 30), "FordB" => Dict("1" => 2185, "-1" => 2261), "FaceAll" => Dict("3" => 176, "4" => 171, "1" => 112, "12" => 106, "2" => 178, "6" => 187, "7" => 149, "11" => 48, "13" => 327, "8"
=> 273, "10" => 135, "5" => 176, "9" => 140, "14" => 72), "ToeSegmentation2" => Dict("1" => 42, "0" => 124), "CinCECGTorso" => Dict("4" => 355, "1" => 355, "2" => 355, "3" => 355), "MiddlePhalanxOutlineAgeGroup" => Dict("1" => 92, "2" => 194, "3" => 266), "ProximalPhalanxOutlineCorrect" => Dict("1" => 603, "0" => 286), "AllGestureWiimoteX" => Dict("8" => 100, "4" => 100, "1" => 100, "10" => 100, "2" => 100, "5" => 100, "6" => 100, "7" => 100, "9" => 100, "3" => 100), "SemgHandSubjectCh2" => Dict("4" => 180, "1" => 180, "5" => 180, "2" => 180, "3" => 180), "UWaveGestureLibraryAll" => Dict("8" => 560, "4" => 560, "1" => 559, "5" => 559, "2" => 560, "6" => 559, "7" => 559, "3" => 560), "Fungi" => Dict("3" => 8, "4" => 20, "1" => 8, "12" => 16, "2" => 12, "6" => 12, "7" => 12, "11" => 8, "13" => 8, "8" => 8, "17" => 8, "10" => 8, "15" => 12, "5" => 8, "16" => 12, "18" => 16, "9" => 20, "14" => 8), "Herring" => Dict("1" => 77, "2" => 51), "DodgerLoopWeekend" => Dict("1" => 112, "2" => 46), "PigAirwayPressure" => Dict("32" => 6, "29" => 6, "1" => 6, "12" => 6, "9" => 6, "20" => 6, "2" =>
6, "6" => 6, "41" => 6, "51" => 6, "27" => 6, "42" => 6, "33" => 6, "28" => 6, "50" => 6, "52" => 6, "14" => 6, "7" => 6, "8" => 6, "26" => 6, "10" => 6, "19" => 6, "24" => 6, "4" => 6, "25" => 6, "23" => 6, "22" => 6, "11" => 6, "35" => 6, "13" => 6, "49" => 6, "15" => 6, "31" => 6, "5" => 6, "43" => 6, "16" =>
6, "39" => 6, "40" => 6, "21" => 6, "38" => 6, "46" => 6, "34" => 6, "44" => 6, "45" => 6, "36" => 6, "17" => 6, "37" => 6, "47" => 6, "18" => 6, "48" => 6,
"30" => 6, "3" => 6), "ProximalPhalanxTW" => Dict("8" => 251, "4" => 176, "5" => 98, "6" => 30, "7" => 31, "3" => 17), "PLAID" => Dict("4" => 156, "1" => 175, "2" => 115, "6" => 114, "7" => 172, "8" => 139, "10" => 26, "0" => 66, "5" => 35, "9" => 38, "3" => 38), "ChlorineConcentration" => Dict("1" => 1000, "2" => 1000, "3" => 2307), "PigArtPressure" => Dict("32" => 6, "29" => 6, "1" => 6, "12" => 6, "9" => 6, "20" => 6, "2" => 6, "6" => 6, "41" => 6, "51" => 6, "27" => 6, "42" => 6, "33" => 6, "28" => 6, "50" => 6, "52" => 6, "14" => 6, "7" => 6, "8" => 6, "26" => 6, "10" => 6, "19" => 6, "24" => 6, "4" => 6, "25" => 6, "23" => 6, "22" => 6, "11" => 6, "35" => 6, "13" => 6, "49" => 6, "15" => 6, "31" => 6, "5" => 6, "43" => 6, "16" => 6, "39" => 6, "40" => 6, "21" => 6, "38" => 6, "46" => 6, "34" => 6, "44" => 6, "45" => 6, "36" => 6, "17" => 6, "37" => 6, "47" => 6, "18" => 6, "48" => 6, "30" => 6, "3" => 6), "FiftyWords" => Dict("32" => 9, "29" => 9, "1" => 109, "12" => 22, "9" => 23, "20" => 13, "2" => 91, "6" => 34, "41" => 7, "27" => 10, "42" => 7, "33" => 8, "28" => 9, "50" => 6, "14" => 18, "7" => 32, "8" => 24, "26" => 10, "10" => 22, "19" => 13, "24" => 10, "4" => 53, "25" => 10, "23" => 10, "22" => 11, "11" => 22, "35" => 8,
"13" => 21, "49" => 6, "15" => 16, "31" => 9, "5" => 38, "43" => 7, "16" => 14, "39" => 8, "40" => 8, "21" => 12, "38" => 8, "46" => 7, "34" => 8, "44" => 7, "45" => 7, "36" => 8, "17" => 14, "37" => 8, "47" => 7, "18" => 13, "48" => 7, "30" => 9, "3" => 61), "EthanolLevel" => Dict("4" => 250, "1" => 252, "2" =>
250, "3" => 252), "MiddlePhalanxOutlineCorrect" => Dict("1" => 553, "0" => 336), "ElectricDevices" => Dict("4" => 2639, "1" => 1394, "5" => 4275, "2" => 4187, "6" => 1252, "7" => 1284, "3" => 1606), "UWaveGestureLibraryY" => Dict("8" => 560, "4" => 560, "1" => 559, "5" => 559, "2" => 560, "6" => 559, "7" => 559,
"3" => 560), "SonyAIBORobotSurface1" => Dict("1" => 349, "2" => 272), "DistalPhalanxOutlineAgeGroup" => Dict("1" => 45, "2" => 170, "3" => 322), "GesturePebbleZ1" => Dict("4" => 53, "1" => 49, "5" => 52, "2" => 45, "6" => 51, "3" => 54), "ToeSegmentation1" => Dict("1" => 128, "0" => 140), "DistalPhalanxTW" => Dict("8" => 254, "4" => 46, "5" => 121, "6" => 43, "7" => 26, "3" => 47), "LargeKitchenAppliances" => Dict("1" => 250, "2" => 250, "3" => 250), "NonInvasiveFetalECGThorax1" => Dict("24" => 87, "29" => 96, "1" => 84, "12" => 90, "32" => 90, "20" => 90, "4" => 85, "2" => 91, "25" => 90, "6" => 94, "23" => 87, "22" =>
86, "11" => 75, "35" => 94, "13" => 97, "41" => 93, "15" => 98, "27" => 89, "31" => 94, "42" => 93, "33" => 85, "28" => 85, "16" => 95, "5" => 96, "14" => 89, "39" => 90, "40" => 91, "21" => 89, "38" => 91, "7" => 93, "34" => 94, "8" => 91, "36" => 91, "26" => 78, "17" => 85, "10" => 92, "19" => 88, "37" => 87, "18" => 91, "9" => 88, "30" => 86, "3" => 87), "Car" => Dict("4" => 30, "1" => 30, "2" => 30, "3" => 30), "RefrigerationDevices" => Dict("1" => 250, "2" => 250, "3" => 250), "ItalyPowerDemand" => Dict("1" => 545, "2" => 547), "MelbournePedestrian" => Dict("8" => 349, "4" => 365, "1" => 365, "10" => 365, "2" => 365, "5" => 365, "6" => 365, "7" => 365, "9" => 364, "3" => 365), "ShapeletSim" => Dict("1" => 100, "0" => 100), "MoteStrain" => Dict("1" => 685, "2" => 587), "Lightning2" => Dict("1" => 73, "-1" => 48))

# multivariate datasets are not yet supported
TSAD_MULTIVARIATE = String[]


text:
lstm(; n=0, fbias=0) =
quote
    x  = input()
    i  = add2(x,h;n=$n,f=sigm)
    f  = add2(x,h;n=$n,f=sigm,b=$fbias)
    o  = add2(x,h;n=$n,f=sigm)
    g  = add2(x,h;n=$n,f=tanh)
    ig = mul(i,g)
    fc = mul(f,c)
    c  = add(ig,fc)
    tc = tanh(c)
    h  = mul(tc,o)
end

add2(; n=1,f=sigm,b=0) =
quote
    x1 = input()
    w1 = par($n,0)
    y1 = dot(w1,x1)
    x2 = input()
    w2 = par($n,0)
    y2 = dot(w2,x2)
    x3 = add(y1,y2)
    w3 = par(0; init=Constant($b))
    y3 = add(w3,x3)
    ou = $f(y3)
end

### DEAD CODE

# using Knet

# function LSTM(n; fbias=0, dropout=0)
#     Net(
#         Drop(dropout),	# 1. drop

#         (Mmul(n),1),
#         (Mmul(n),26),
#         Add2(),
#         Bias(),
#         Sigm(),         # 6. input

#         (Mmul(n),1),
#         (Mmul(n),26),
#         Add2(),
#         Bias(init=fill!, initp=fbias),
#         Sigm(),         # 11. forget

#         (Mmul(n),1),
#         (Mmul(n),26),
#         Add2(),
#         Bias(),
#         Sigm(),         # 16. output

#         (Mmul(n),1),
#         (Mmul(n),26),
#         Add2(),
#         Bias(),
#         Tanh(),         # 21. g

#         (Mul2(),6,21),
#         (Mul2(),11,24),
#         Add2(),         # 24. c = i*g + f*c'

#         Tanh(),
#         (Mul2(),16,25)  # 26. h = o*tanh(c)
#     )
# end

# # More compact definition without dropout:
# function LSTM1(n; fbias=0)
#     m = Net((add2(n),0,13), Sigm(),     	# 1-2. input  (op 1-5)
#             (add2(n),0,13), Sigm(),             # 3-4. forget (op 6-10)
#             (add2(n),0,13), Sigm(),             # 5-6. output (op 11-15)
#             (add2(n),0,13), Tanh(),             # 7-8. g      (op 16-20)
#             (Mul2(),2,8),(Mul2(),4,11),Add2(),  # 9-11. c = i*g + f*c[t-1] (op 21-23)
#             Tanh(), (Mul2(),6,12))              # 12-13. h = o * tanh(c) (op 24-25)
#     fbias != 0 && setparam!(m.op[9]; init=fill!, initp=fbias)
#     return m
# end

# add2(n)=Net(Mmul(n), (Mmul(n),-1), Add2(), Bias())


text:
"""
    RowVector(vector)

A lazy-view wrapper of an `AbstractVector`, which turns a length-`n` vector into a `1×n`
shaped row vector and represents the transpose of a vector (the elements are also transposed
recursively). This type is usually constructed (and unwrapped) via the [`transpose`](@ref)
function or `.'` operator (or related [`ctranspose`](@ref) or `'` operator).

By convention, a vector can be multiplied by a matrix on its left (`A * v`) whereas a row
vector can be multiplied by a matrix on its right (such that `v.' * A = (A.' * v).'`). It
differs from a `1×n`-sized matrix by the facts that its transpose returns a vector and the
inner product `v1.' * v2` returns a scalar, but will otherwise behave similarly.
"""
struct RowVector{T,V<:AbstractVector} <: AbstractMatrix{T}
    vec::V
    function RowVector{T,V}(v::V) where V<:AbstractVector where T
        check_types(T,v)
        new(v)
    end
end

@inline check_types{T1,T2}(::Type{T1},::AbstractVector{T2}) = check_types(T1, T2)
@pure check_types{T1,T2}(::Type{T1},::Type{T2}) = T1 === transpose_type(T2) ? nothing :
    error("Element type mismatch. Tried to create a `RowVector{$T1}` from an `AbstractVector{$T2}`")

ConjRowVector{T, CV <: ConjVector} = RowVector{T, CV}

# The element type may be transformed as transpose is recursive
@inline transpose_type{T}(::Type{T}) = promote_op(transpose, T)

# Constructors that take a vector
@inline RowVector(vec::AbstractVector{T}) where {T} = RowVector{transpose_type(T),typeof(vec)}(vec)
@inline RowVector{T}(vec::AbstractVector{T}) where {T} = RowVector{T,typeof(vec)}(vec)

# Constructors that take a size and default to Array
@inline RowVector{T}(n::Int) where {T} = RowVector{T}(Vector{transpose_type(T)}(n))
@inline RowVector{T}(n1::Int, n2::Int) where {T} = n1 == 1 ?
    RowVector{T}(Vector{transpose_type(T)}(n2)) :
    error("RowVector expects 1×N size, got ($n1,$n2)")
@inline RowVector{T}(n::Tuple{Int}) where {T} = RowVector{T}(Vector{transpose_type(T)}(n[1]))
@inline RowVector{T}(n::Tuple{Int,Int}) where {T} = n[1] == 1 ?
    RowVector{T}(Vector{transpose_type(T)}(n[2])) :
    error("RowVector expects 1×N size, got $n")

# Conversion of underlying storage
convert{T,V<:AbstractVector}(::Type{RowVector{T,V}}, rowvec::RowVector) =
    RowVector{T,V}(convert(V,rowvec.vec))

# similar tries to maintain the RowVector wrapper and the parent type
@inline similar(rowvec::RowVector) = RowVector(similar(parent(rowvec)))
@inline similar{T}(rowvec::RowVector, ::Type{T}) = RowVector(similar(parent(rowvec), transpose_type(T)))

# Resizing similar currently loses its RowVector property.
@inline similar{T,N}(rowvec::RowVector, ::Type{T}, dims::Dims{N}) = similar(parent(rowvec), T, dims)

# Basic methods
"""
    transpose(v::AbstractVector)

The transposition operator (`.'`).

# Example

```jldoctest
julia> v = [1,2,3]
3-element Array{Int64,1}:
 1
 2
 3

julia> transpose(v)
1×3 RowVector{Int64,Array{Int64,1}}:
 1  2  3
```
"""
@inline transpose(vec::AbstractVector) = RowVector(vec)
@inline ctranspose(vec::AbstractVector) = RowVector(_conj(vec))
@inline ctranspose(vec::AbstractVector{<:Real}) = RowVector(vec)

@inline transpose(rowvec::RowVector) = rowvec.vec
@inline transpose(rowvec::ConjRowVector) = copy(rowvec.vec) # remove the ConjArray wrapper from any raw vector
@inline ctranspose(rowvec::RowVector) = conj(rowvec.vec)
@inline ctranspose(rowvec::RowVector{<:Real}) = rowvec.vec

parent(rowvec::RowVector) = rowvec.vec

"""
    conj(v::RowVector)

Returns a [`ConjArray`](@ref) lazy view of the input, where each element is conjugated.

### Example

```jldoctest
julia> v = [1+im, 1-im].'
1×2 RowVector{Complex{Int64},Array{Complex{Int64},1}}:
 1+1im  1-1im

julia> conj(v)
1×2 RowVector{Complex{Int64},ConjArray{Complex{Int64},1,Array{Complex{Int64},1}}}:
 1-1im  1+1im
```
"""
@inline conj(rowvec::RowVector) = RowVector(_conj(rowvec.vec))
@inline conj(rowvec::RowVector{<:Real}) = rowvec

# AbstractArray interface
@inline length(rowvec::RowVector) =  length(rowvec.vec)
@inline size(rowvec::RowVector) = (1, length(rowvec.vec))
@inline size(rowvec::RowVector, d) = ifelse(d==2, length(rowvec.vec), 1)
@inline indices(rowvec::RowVector) = (Base.OneTo(1), indices(rowvec.vec)[1])
@inline indices(rowvec::RowVector, d) = ifelse(d == 2, indices(rowvec.vec)[1], Base.OneTo(1))
IndexStyle(::RowVector) = IndexLinear()
IndexStyle(::Type{<:RowVector}) = IndexLinear()

@propagate_inbounds getindex(rowvec::RowVector, i) = transpose(rowvec.vec[i])
@propagate_inbounds setindex!(rowvec::RowVector, v, i) = setindex!(rowvec.vec, transpose(v), i)

# Cartesian indexing is distorted by getindex
# Furthermore, Cartesian indexes don't have to match shape, apparently!
@inline function getindex(rowvec::RowVector, i::CartesianIndex)
    @boundscheck if !(i.I[1] == 1 && i.I[2] ∈ indices(rowvec.vec)[1] && check_tail_indices(i.I...))
        throw(BoundsError(rowvec, i.I))
    end
    @inbounds return transpose(rowvec.vec[i.I[2]])
end
@inline function setindex!(rowvec::RowVector, v, i::CartesianIndex)
    @boundscheck if !(i.I[1] == 1 && i.I[2] ∈ indices(rowvec.vec)[1] && check_tail_indices(i.I...))
        throw(BoundsError(rowvec, i.I))
    end
    @inbounds rowvec.vec[i.I[2]] = transpose(v)
end

@propagate_inbounds getindex(rowvec::RowVector, ::CartesianIndex{0}) = getindex(rowvec)
@propagate_inbounds getindex(rowvec::RowVector, i::CartesianIndex{1}) = getindex(rowvec, i.I[1])

@propagate_inbounds setindex!(rowvec::RowVector, v, ::CartesianIndex{0}) = setindex!(rowvec, v)
@propagate_inbounds setindex!(rowvec::RowVector, v, i::CartesianIndex{1}) = setindex!(rowvec, v, i.I[1])

@inline check_tail_indices(i1, i2) = true
@inline check_tail_indices(i1, i2, i3, is...) = i3 == 1 ? check_tail_indices(i1, i2, is...) : false

# helper function for below
@inline to_vec(rowvec::RowVector) = transpose(rowvec)
@inline to_vec(x::Number) = x
@inline to_vecs(rowvecs...) = (map(to_vec, rowvecs)...)

# map
@inline map(f, rowvecs::RowVector...) = RowVector(map(f, to_vecs(rowvecs...)...))

# broacast (other combinations default to higher-dimensional array)
@inline broadcast(f, rowvecs::Union{Number,RowVector}...) =
    RowVector(broadcast(f, to_vecs(rowvecs...)...))

# Horizontal concatenation #

@inline hcat(X::RowVector...) = transpose(vcat(map(transpose, X)...))
@inline hcat(X::Union{RowVector,Number}...) = transpose(vcat(map(transpose, X)...))

@inline typed_hcat{T}(::Type{T}, X::RowVector...) =
    transpose(typed_vcat(T, map(transpose, X)...))
@inline typed_hcat{T}(::Type{T}, X::Union{RowVector,Number}...) =
    transpose(typed_vcat(T, map(transpose, X)...))

# Multiplication #

# inner product -> dot product specializations
@inline *{T<:Real}(rowvec::RowVector{T}, vec::AbstractVector{T}) = dot(parent(rowvec), vec)
@inline *(rowvec::ConjRowVector, vec::AbstractVector) = dot(rowvec', vec)

# Generic behavior
@inline function *(rowvec::RowVector, vec::AbstractVector)
    if length(rowvec) != length(vec)
        throw(DimensionMismatch("A has dimensions $(size(rowvec)) but B has dimensions $(size(vec))"))
    end
    sum(@inbounds(return rowvec[i]*vec[i]) for i = 1:length(vec))
end
@inline *(rowvec::RowVector, mat::AbstractMatrix) = transpose(mat.' * transpose(rowvec))
*(::RowVector, ::RowVector) = throw(DimensionMismatch("Cannot multiply two transposed vectors"))
@inline *(vec::AbstractVector, rowvec::RowVector) = vec .* rowvec
*(vec::AbstractVector, rowvec::AbstractVector) = throw(DimensionMismatch("Cannot multiply two vectors"))

# Transposed forms
A_mul_Bt(::RowVector, ::AbstractVector) = throw(DimensionMismatch("Cannot multiply two transposed vectors"))
@inline A_mul_Bt(rowvec::RowVector, mat::AbstractMatrix) = transpose(mat * transpose(rowvec))
@inline A_mul_Bt(rowvec1::RowVector, rowvec2::RowVector) = rowvec1*transpose(rowvec2)
A_mul_Bt(vec::AbstractVector, rowvec::RowVector) = throw(DimensionMismatch("Cannot multiply two vectors"))
@inline A_mul_Bt(vec1::AbstractVector, vec2::AbstractVector) = vec1 * transpose(vec2)
@inline A_mul_Bt(mat::AbstractMatrix, rowvec::RowVector) = mat * transpose(rowvec)

@inline At_mul_Bt(rowvec::RowVector, vec::AbstractVector) = transpose(rowvec) * transpose(vec)
@inline At_mul_Bt(vec::AbstractVector, mat::AbstractMatrix) = transpose(mat * vec)
At_mul_Bt(rowvec1::RowVector, rowvec2::RowVector) = throw(DimensionMismatch("Cannot multiply two vectors"))
@inline At_mul_Bt(vec::AbstractVector, rowvec::RowVector) = transpose(vec)*transpose(rowvec)
At_mul_Bt(vec::AbstractVector, rowvec::AbstractVector) = throw(DimensionMismatch(
    "Cannot multiply two transposed vectors"))
@inline At_mul_Bt(mat::AbstractMatrix, rowvec::RowVector) = mat.' * transpose(rowvec)

At_mul_B(::RowVector, ::AbstractVector) = throw(DimensionMismatch("Cannot multiply two vectors"))
@inline At_mul_B(vec::AbstractVector, mat::AbstractMatrix) = transpose(At_mul_B(mat,vec))
@inline At_mul_B(rowvec1::RowVector, rowvec2::RowVector) = transpose(rowvec1) * rowvec2
At_mul_B(vec::AbstractVector, rowvec::RowVector) = throw(DimensionMismatch(
    "Cannot multiply two transposed vectors"))
@inline At_mul_B{T<:Real}(vec1::AbstractVector{T}, vec2::AbstractVector{T}) =
    reduce(+, map(At_mul_B, vec1, vec2)) # Seems to be overloaded...
@inline At_mul_B(vec1::AbstractVector, vec2::AbstractVector) = transpose(vec1) * vec2

# Conjugated forms
A_mul_Bc(::RowVector, ::AbstractVector) = throw(DimensionMismatch("Cannot multiply two transposed vectors"))
@inline A_mul_Bc(rowvec::RowVector, mat::AbstractMatrix) = ctranspose(mat * ctranspose(rowvec))
@inline A_mul_Bc(rowvec1::RowVector, rowvec2::RowVector) = rowvec1 * ctranspose(rowvec2)
A_mul_Bc(vec::AbstractVector, rowvec::RowVector) = throw(DimensionMismatch("Cannot multiply two vectors"))
@inline A_mul_Bc(vec1::AbstractVector, vec2::AbstractVector) = vec1 * ctranspose(vec2)
@inline A_mul_Bc(mat::AbstractMatrix, rowvec::RowVector) = mat * ctranspose(rowvec)

@inline Ac_mul_Bc(rowvec::RowVector, vec::AbstractVector) = ctranspose(rowvec) * ctranspose(vec)
@inline Ac_mul_Bc(vec::AbstractVector, mat::AbstractMatrix) = ctranspose(mat * vec)
Ac_mul_Bc(rowvec1::RowVector, rowvec2::RowVector) = throw(DimensionMismatch("Cannot multiply two vectors"))
@inline Ac_mul_Bc(vec::AbstractVector, rowvec::RowVector) = ctranspose(vec)*ctranspose(rowvec)
Ac_mul_Bc(vec::AbstractVector, rowvec::AbstractVector) = throw(DimensionMismatch("Cannot multiply two transposed vectors"))
@inline Ac_mul_Bc(mat::AbstractMatrix, rowvec::RowVector) = mat' * ctranspose(rowvec)

Ac_mul_B(::RowVector, ::AbstractVector) = throw(DimensionMismatch("Cannot multiply two vectors"))
@inline Ac_mul_B(vec::AbstractVector, mat::AbstractMatrix) = ctranspose(Ac_mul_B(mat,vec))
@inline Ac_mul_B(rowvec1::RowVector, rowvec2::RowVector) = ctranspose(rowvec1) * rowvec2
Ac_mul_B(vec::AbstractVector, rowvec::RowVector) = throw(DimensionMismatch("Cannot multiply two transposed vectors"))
@inline Ac_mul_B(vec1::AbstractVector, vec2::AbstractVector) = ctranspose(vec1)*vec2

# Left Division #

\(mat::AbstractMatrix, rowvec::RowVector) = throw(DimensionMismatch("Cannot left-divide transposed vector by matrix"))
At_ldiv_B(mat::AbstractMatrix, rowvec::RowVector) = throw(DimensionMismatch("Cannot left-divide transposed vector by matrix"))
Ac_ldiv_B(mat::AbstractMatrix, rowvec::RowVector) = throw(DimensionMismatch("Cannot left-divide transposed vector by matrix"))

# Right Division #

@inline /(rowvec::RowVector, mat::AbstractMatrix) = transpose(transpose(mat) \ transpose(rowvec))
@inline A_rdiv_Bt(rowvec::RowVector, mat::AbstractMatrix) = transpose(mat \ transpose(rowvec))
@inline A_rdiv_Bc(rowvec::RowVector, mat::AbstractMatrix) = ctranspose(mat  \ ctranspose(rowvec))

max_stars_repo_name: felixrehren/julia

max_issues_repo_name: felixrehren/julia

max_forks_repo_name: felixrehren/julia


text:
@testset "corenums" begin
    A = load_matrix_network("cores_example")
    A_ref = [1 0 1 1; 0 1 1 1; 1 1 0 1; 1 1 1 1]
    d_ref = [3;3;3;3]
    rt_ref = [1;2;3;4]
    @testset "MatrixNetwork" begin
        corenums(MatrixNetwork(A))
        (d,rt) = corenums(MatrixNetwork(sparse(A_ref)))
        @test d_ref == d
        @test rt_ref == rt
    end
    @testset "SparseMatrixCSC" begin
        corenums(sparse(A))
        (d,rt) = corenums(sparse(A_ref))
        @test d_ref == d
        @test rt_ref == rt
    end
    @testset "CSR" begin
        corenums(sparse_to_csr(sparse(A))...)
        (d,rt) = corenums(sparse_to_csr(sparse(A_ref))...)
        @test d_ref == d
        @test rt_ref == rt
    end
    @testset "triplet" begin
        corenums(findnz(sparse(A_ref))[1], findnz(sparse(A_ref))[2])
        (d,rt) = corenums(findnz(sparse(A_ref))[1], findnz(sparse(A_ref))[2])
        @test d_ref == d
        @test rt_ref == rt
    end
end


text:
#Create constes like Mat4f0, Point2, Point2f0
for i=1:4
    for T=[:Point, :Vec]
        name = Symbol("$T$i")
        namef0 = Symbol("$T$(i)f0")
        @eval begin
            const $name = $T{$i}
            const $namef0 = $T{$i, Float32}
            export $name
            export $namef0
        end
    end
    name   = Symbol("Mat$i")
    namef0 = Symbol("Mat$(i)f0")
    @eval begin
        const $name{T} = $Mat{$i,$i, T, $(i*i)}
        const $namef0 = $name{Float32}
        export $name
        export $namef0
    end
end

const Vecf0{N} = Vec{N, Float32}
const Pointf0{N} = Point{N, Float32}
export Vecf0, Pointf0

#Type aliases
"""
An alias for a one-simplex, corresponding to `LineSegment{T}` -> `Simplex{2,T}`.
"""
const LineSegment{T} = Simplex{2,T}

const ZeroIndex{T <: Integer} = OffsetInteger{-1, T}
const OneIndex{T <: Integer} = OffsetInteger{0, T}
const GLIndex = ZeroIndex{Cuint}

const Triangle{T} = Face{3, T}
const GLFace{Dim} = Face{Dim, GLIndex} #offset is relative to julia, so -1 is 0-indexed
const GLTriangle = Face{3, GLIndex}
const GLQuad = Face{4, GLIndex}

const Cube{T} = HyperCube{3, T}

"""
An alias for a HyperSphere of dimension 2. i.e. `Circle{T}` -> `HyperSphere{2, T}`
"""
const Circle{T} = HyperSphere{2, T}

"""
An alias for a HyperSphere of dimension 3. i.e. `Sphere{T}` -> `HyperSphere{3, T}`
"""
const Sphere{T} = HyperSphere{3, T}

const Rect{N, T} = HyperRectangle{N, T}
const Rect2D{T} = Rect{2, T}
const Rect3D{T} = Rect{3, T}
const TRect{T} = Rect{N, T} where N

const FRect{N} = Rect{N, Float32}
const FRect2D = Rect2D{Float32}
const FRect3D = Rect3D{Float32}

const IRect{N} = HyperRectangle{N, Int}
const IRect2D = Rect2D{Int}
const IRect3D = Rect3D{Int}

const VecTypes{N, T} = Union{StaticVector{N, T}, NTuple{N, T}}

"""
AABB, or Axis Aligned Bounding Box, is an alias for a 3D `HyperRectangle`.
"""
const AABB{T} = HyperRectangle{3, T}

const HMesh = HomogenousMesh

"""
A `Cylinder2` or `Cylinder3` is a 2D/3D cylinder defined by its origin point,
its extremity and a radius. `origin`, `extremity` and `r`, must be specified.
"""
const Cylinder2{T} = Cylinder{2, T}
const Cylinder3{T} = Cylinder{3, T}

"""
A `Capsule2` or `Capsule3` is a 2D/3D capsule defined by its origin point,
its extremity and a radius. `origin`, `extremity` and `r`, must be specified.
"""
const Capsule2{T} = Capsule{2, T}
const Capsule3{T} = Capsule{3, T}

const UV{T} = TextureCoordinate{2, T}
const UVW{T} = TextureCoordinate{3, T}

"""
A `SimpleMesh` is an alias for a `HomogenousMesh` parameterized only by
vertex and face types.
"""
const SimpleMesh{VT, FT} = HMesh{VT, FT, Nothing, Nothing, Nothing, Nothing, Nothing}
const PlainMesh{VT, FT} = HMesh{Point{3, VT}, FT, Nothing, Nothing, Nothing, Nothing, Nothing}
const GLPlainMesh = PlainMesh{Float32, GLTriangle}

const Mesh2D{VT, FT} = HMesh{Point{2, VT}, FT, Nothing, Nothing, Nothing, Nothing, Nothing}
const GLMesh2D = Mesh2D{Float32, GLTriangle}

const UVMesh{VT, FT, UVT} = HMesh{Point{3, VT}, FT, Nothing, UV{UVT}, Nothing, Nothing, Nothing}
const GLUVMesh = UVMesh{Float32, GLTriangle, Float32}

const UVWMesh{VT, FT, UVT} = HMesh{Point{3, VT}, FT, Nothing, UVW{UVT}, Nothing, Nothing, Nothing}
const GLUVWMesh = UVWMesh{Float32, GLTriangle, Float32}

const NormalMesh{VT, FT, NT} = HMesh{Point{3, VT}, FT, Normal{3, NT}, Nothing, Nothing, Nothing, Nothing}
const GLNormalMesh = NormalMesh{Float32, GLTriangle, Float32}

const UVMesh2D{VT, FT, UVT} = HMesh{Point{2, VT}, FT, Nothing, UV{UVT}, Nothing, Nothing, Nothing}
const GLUVMesh2D = UVMesh2D{Float32, GLTriangle, Float32}

const NormalColorMesh{VT, FT, NT, CT} = HMesh{Point{3, VT}, FT, Normal{3, NT}, Nothing, CT, Nothing, Nothing}
const GLNormalColorMesh = NormalColorMesh{Float32, GLTriangle, Float32, RGBA{Float32}}

const NormalVertexcolorMesh{VT, FT, NT, CT} = HMesh{Point{3, VT}, FT, Normal{3, NT}, Nothing, Vector{CT}, Nothing, Nothing}
const GLNormalVertexcolorMesh = NormalVertexcolorMesh{Float32, GLTriangle, Float32, RGBA{Float32}}

const NormalAttributeMesh{VT, FT, NT, AT, A_ID_T} = HMesh{Point{3, VT}, FT, Normal{3, NT}, Nothing, Nothing, AT, A_ID_T}
const GLNormalAttributeMesh = NormalAttributeMesh{Float32, GLTriangle, Float32, Vector{RGBA{N0f8}}, Float32}

const NormalUVWMesh{VT, FT, NT, UVT} = HMesh{Point{3, VT}, FT, Normal{3, NT}, UVW{UVT}, Nothing, Nothing, Nothing}
const GLNormalUVWMesh = NormalUVWMesh{Float32, GLTriangle, Float32, Float32}

const NormalUVMesh{VT, FT, NT, UVT} = HMesh{Point{3, VT}, FT, Normal{3, NT}, UV{UVT}, Nothing, Nothing, Nothing}
const GLNormalUVMesh = NormalUVMesh{Float32, GLTriangle, Float32, Float32}

for T in (
        :SimpleMesh,
        :PlainMesh,
        :GLPlainMesh,
        :Mesh2D,
        :GLMesh2D,
        :UVMesh,
        :GLUVMesh,
        :UVWMesh,
        :GLUVWMesh,
        :NormalMesh,
        :GLNormalMesh,
        :UVMesh2D,
        :GLUVMesh2D,
        :NormalColorMesh,
        :GLNormalColorMesh,
        :NormalVertexcolorMesh,
        :GLNormalVertexcolorMesh,
        :NormalAttributeMesh,
        :GLNormalAttributeMesh,
        :NormalUVWMesh,
        :GLNormalUVWMesh,
        :NormalUVMesh,
        :GLNormalUVMesh,
    )
    @eval Base.show(io::IO, ::Type{<: $T}) = print(io, $(string(T)))
end


text:
 using GPUifyLoops: @loop

#####
##### Useful kernels for doing diagnostics
#####

function velocity_div!(grid, u, v, w, div)
    @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
        @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
            @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
                @inbounds div[i, j, k] = divᶜᶜᶜ(i, j, k, grid, u, v, w)
            end
        end
    end
end


text:
module ErrorStateKalmanFilter

using LinearAlgebra
using Quaternions

export
    NominalState,
    ErrorState,
    Measurement,
    Noise,
    skew_matrix,
    predict

include("types.jl")
include("ESKF.jl")

end


text:
@testset "P-n19-k2 by coordinates" begin
    cvrp, _, _ = CVRPLIB.readCVRPLIB("P-n19-k2")
    ap = AlgorithmParameters(timeLimit=1.8)
    result = solve_cvrp(cvrp, ap; verbose=true)
    @test result.cost <= 212 * 1.01
end

@testset "P-n19-k2 by dist_mtx with coordinates" begin
    cvrp, _, _ = CVRPLIB.readCVRPLIB("P-n19-k2")
    ap = AlgorithmParameters(timeLimit=2.1)
    result = solve_cvrp(cvrp, ap; verbose=true, use_dist_mtx=true)
    @test result.cost <= 212 * 1.01
end

@testset "P-n19-k2 by dist_mtx without coordinates" begin
    cvrp, _, _ = CVRPLIB.readCVRPLIB("P-n19-k2")
    cvrp.coordinates = zeros(size(cvrp.coordinates))
    ap = AlgorithmParameters(timeLimit=1.23)
    result = solve_cvrp(cvrp, ap; verbose=true, use_dist_mtx=true)
    @test result.cost <= 212 * 1.01
end

@testset "X-n101-k25" begin
    cvrp, _, _ = CVRPLIB.readCVRPLIB("X-n101-k25")

    ap = AlgorithmParameters(timeLimit=3)
    result = solve_cvrp(cvrp, ap; verbose=true)

    @show result.cost
    @show result.time
    visited_customers = result.routes

    @show typeof(visited_customers), size(visited_customers)

    report_visited_customers = reporting(visited_customers)
    for i in eachindex(report_visited_customers)
        @test report_visited_customers[i] == visited_customers[i] .- 1
    end
end

@testset "A-n32-k5.vrp file read" begin
    cvrp = CVRPLIB.readCVRP("A-n32-k5.vrp")
    ap = AlgorithmParameters(timeLimit=2.3)
    result = solve_cvrp(cvrp, ap; verbose=true)
    @test result.cost <= 784 * 1.01

    result = solve_cvrp("A-n32-k5.vrp", ap)
    @test result.cost <= 784 * 1.01
end

@testset "x, y, dist_mtx CVRP" begin
    cvrp, _, _ = CVRPLIB.readCVRPLIB("A-n32-k5")
    x = cvrp.coordinates[:, 1]
    y = cvrp.coordinates[:, 2]
    dist_mtx = cvrp.weights
    service_time = zeros(cvrp.dimension)
    demand = cvrp.demand
    capacity = cvrp.capacity
    n_vehicles = 5

    ap = AlgorithmParameters(timeLimit=1.8)

    result1 = solve_cvrp(x, y, service_time, demand, cvrp.capacity, n_vehicles, ap; verbose=true)
    result2 = solve_cvrp(dist_mtx, service_time, demand, cvrp.capacity, n_vehicles, ap; verbose=true)
    result3 = solve_cvrp(dist_mtx, service_time, demand, cvrp.capacity, n_vehicles, ap; x_coords=x, y_coords=y, verbose=true)

    @test result1.cost == result2.cost == result3.cost
end


text:
export admm

mutable struct ADMM{rT,matT,opT,ropT,vecT,rvecT,preconT} <: AbstractLinearSolver where {vecT <: AbstractVector{Union{rT, Complex{rT}}}, rvecT <: AbstractVector{rT}}
  # oerators and regularization
  A::matT
  reg::Vector{Regularization}
  regTrafo::Vector{ropT}
  # fields and operators for x update
  op::opT
  β::vecT
  β_y::vecT
  # fields for primal & dual variables
  x::vecT
  z::Vector{vecT}
  zᵒˡᵈ::Vector{vecT}
  u::Vector{vecT}
  # other parameters
  precon::preconT
  ρ::rvecT # TODO: Switch all these vectors to Tuple
  iterations::Int64
  iterationsInner::Int64
  # state variables for CG
  cgStateVars::CGStateVariables
  # convergence parameters
  rᵏ::rvecT
  sᵏ::vecT
  ɛᵖʳⁱ::rvecT
  ɛ_dt::vecT
  σᵃᵇˢ::rT
  absTol::rT
  relTol::rT
  tolInner::rT
  normalizeReg::Bool
  regFac::rT
end

"""
    ADMM(A, x::vecT=zeros(eltype(A),size(A,2))
          ; reg=nothing, regName=["L1"], λ=[0.0], kargs...)

creates an `ADMM` object for the system matrix `A`.

# Arguments
* `A`                           - system matrix
* `x::vecT`                     - Array with the same type and size as the solution
* (`reg=nothing`)               - Regularization object
* (`regName=["L1"]`)            - name of the Regularization to use (if reg==nothing)
* (`λ=[0.0]`)                   - Regularization paramter
* (`regTrafo=nothing`)          - transformations applied inside each regularizer
* (`precon=Identity()`)         - preconditionner for the internal CG algorithm
* (`ρ::Real=1.e-2`)          - penalty of the augmented lagrangian
* (`adaptRho::Bool=false`)      - adapt rho to balance primal and dual feasibility
* (`iterations::Int64=50`)      - max number of ADMM iterations
* (`iterationsInner::Int64=10`) - max number of internal CG iterations
* (`absTol::Real=eps()`)     - abs tolerance for stopping criterion
* (`relTol::Real=eps()`)     - rel tolerance for stopping criterion
* (`tolInner::Real=1.e-5`)   - rel tolerance for CG stopping criterion
"""
function ADMM(A::matT, x::Vector{T}=zeros(eltype(A),size(A,2)); reg=nothing, regName=["L1"]
            , λ=[0.0]
            , regTrafo=nothing
            , AHA::opT=nothing
            , precon=Identity()
            , ρ=1e-1
            , iterations::Integer=50
            , iterationsInner::Integer=10
            , absTol::Real=eps(real(T))
            , relTol::Real=eps(real(T))
            , tolInner::Real=1e-5
            , normalizeReg::Bool=false
            , kargs...) where {T,matT,opT}
# TODO: The constructor is not type stable

  # unify Floating types
  if typeof(ρ) <: Number
    ρ_vec = [real(T).(ρ)]
  else
    ρ_vec = real(T).(ρ)
  end
  λ = real(T).(λ)
  absTol = real(T)(absTol)
  relTol = real(T)(relTol)
  tolInner = real(T)(tolInner)

  if reg == nothing
    reg = Regularization(regName, λ, kargs...)
  end

  if regTrafo == nothing
    regTrafo = [LinearOperator{T}(length(x), length(x), true, true, x -> identity(x), x -> identity(x), x -> identity(x)) for i=1:length(vec(reg))]
  end

  # fields for primal & dual variables
  z = [similar(x, size(regTrafo[i],1)) for i=1:length(vec(reg))]
  zᵒˡᵈ = [similar(z[i]) for i=1:length(vec(reg))]
  u = [similar(z[i]) for i=1:length(vec(reg))]

  # operator and fields for the update of x
  op = ( AHA!=nothing ? AHA : A'*A )
  for i=1:length(vec(reg))
    op += ρ[i]*adjoint(regTrafo[i])*regTrafo[i]
  end
  β = similar(x)
  β_y = similar(x)

  # statevariables for CG
  # we store them here to prevent CG from allocating new fields at each call
  statevars = CGStateVariables(zero(x),similar(x),similar(x))

  # convergence parameters
  rk = similar( real.(x), length(vec(reg)) ) #[0.0 for i=1:length(vec(reg))]
  sk = similar(x)
  eps_pri = similar( real.(x), length(vec(reg)) ) # [0.0 for i=1:length(vec(reg))]
  eps_dt = similar(x)

  return ADMM(A,vec(reg),regTrafo,op,β,β_y,x,z,zᵒˡᵈ,u,precon,ρ_vec,iterations
              ,iterationsInner,statevars, rk,sk,eps_pri,eps_dt,zero(real(T)),absTol,relTol,tolInner
              ,normalizeReg,one(real(T)))
end

"""
  init!(solver::ADMM{matT,opT,vecT,rvecT,preconT}, b::vecT
              ; A::matT=solver.A
              , AHA::opT=solver.op
              , x::vecT=similar(b,0)
              , kargs...) where {matT,opT,vecT,rvecT,preconT}

(re-) initializes the ADMM iterator
"""
function init!(solver::ADMM{rT,matT,opT,ropT,vecT,rvecT,preconT}, b::vecT
              ; A::matT=solver.A
              , AHA::opT=solver.op
              , x::vecT=similar(b,0)
              , kargs...) where {rT,matT,opT,ropT,vecT,rvecT,preconT}

  # operators
  if A != solver.A
    solver.A = A
    solver.op = ( AHA!=nothing ? AHA : A'*A )
    for i=1:length(vec(reg))
      solver.op += ρ[i]*adjoint(regTrafo[i])*regTrafo[i]
    end
  end

  # start vector
  if isempty(x)
    if !isempty(b)
      solver.x .= adjoint(A) * b
    else
      solver.x .= 0
   end
 else
   solver.x[:] .= x
 end

  # primal and dual variables
  for i=1:length(solver.reg)
    solver.z[i][:] .= solver.regTrafo[i]*solver.x
    solver.zᵒˡᵈ[i][:] .= 0
    solver.u[i] .= 0
  end

  # right hand side for the x-update
  solver.β_y[:] .= adjoint(A) * b

  # convergence parameter
  solver.rᵏ .= 0
  solver.sᵏ .= 0
  solver.ɛᵖʳⁱ .= 0
  solver.ɛ_dt .= 0
  solver.σᵃᵇˢ = sqrt(length(b))*solver.absTol

  # normalization of regularization parameters
  if solver.normalizeReg
    solver.regFac = norm(b,1)/length(b)
  else
    solver.regFac = 1
  end
end

"""
    solve(solver::ADMM, b::vecT
          ; A::matT=solver.A
          , startVector::vecT=similar(b,0)
          , solverInfo=nothing
          , kargs...) where {matT,vecT}

solves an inverse problem using ADMM.

# Arguments
* `solver::ADMM`                  - the solver containing both system matrix and regularizer
* `b::Vector`                     - data vector
* (`A::matT=solver.A`)            - operator for the data-term of the problem
* (`startVector::Vector{T}=T[]`)  - initial guess for the solution
* (`solverInfo=nothing`)          - solverInfo for logging

when a `SolverInfo` objects is passed, the primal residuals `solver.rk`
and the dual residual `norm(solver.sk)` are stored in `solverInfo.convMeas`.
"""
function solve(solver::ADMM{rT,matT,opT,ropT,vecT,rvecT,preconT}, b::vecT; A=solver.A, startVector::vecT=similar(b,0), solverInfo=nothing, kargs...) where {rT,matT,opT,ropT,vecT,rvecT,preconT}
  # initialize solver parameters
  init!(solver, b; A=A, x=startVector)

  # log solver information
  solverInfo != nothing && storeInfo(solverInfo,solver.z,solver.rk...,norm(solver.sk))

  # perform ADMM iterations
  for (iteration, item) = enumerate(solver)
    solverInfo != nothing && storeInfo(solverInfo,solver.z,solver.rk...,norm(solver.sk))
  end

  return solver.x
end

"""
  iterate(it::ADMM, iteration::Int=0)

performs one ADMM iteration.
"""
function iterate(solver::ADMM, iteration::Integer=0)
  if done(solver, iteration) return nothing end

  # 1. solve arg min_x 1/2|| Ax-b ||² + ρ/2 Σ_i||Φi*x+ui-zi||²
  # <=> (A'A+ρ Σ_i Φi'Φi)*x = A'b+ρΣ_i Φi'(zi-ui)
  copyto!(solver.β, solver.β_y)
  for i=1:length(solver.reg)
    solver.β[:] .+= solver.ρ[i]*adjoint(solver.regTrafo[i])*(solver.z[i].-solver.u[i])
  end
  cg!(solver.x, solver.op, solver.β, Pl=solver.precon
      , maxiter=solver.iterationsInner, reltol=solver.tolInner, statevars=solver.cgStateVars)

  # 2. update z using the proximal map of 1/ρ*g(x)
  for i=1:length(solver.reg)
    copyto!(solver.zᵒˡᵈ[i], solver.z[i])
    solver.z[i][:] .= solver.regTrafo[i]*solver.x .+ solver.u[i]
    if solver.ρ[i] != 0
      solver.reg[i].prox!(solver.z[i], solver.regFac*solver.reg[i].λ/solver.ρ[i]; solver.reg[i].params...)
    end
  end

  # 3. update u
  for i=1:length(solver.reg)
    solver.u[i][:] .+= solver.regTrafo[i]*solver.x .- solver.z[i]
  end

  # update convergence measures
  # primal residuals norms (one for each constraint)
  for i=1:length(solver.reg)
    solver.rᵏ[i] = norm(solver.regTrafo[i]*solver.x-solver.z[i])  # primal residual (x-z)
    solver.ɛᵖʳⁱ[i] = solver.σᵃᵇˢ + solver.relTol*max( norm(solver.regTrafo[i]*solver.x), norm(solver.z[i]) )
  end
  # accumulated dual residual
  # effectively this corresponds to combining all constraints into one larger constraint.
  solver.sᵏ[:] .= 0.0
  solver.ɛ_dt[:] .= 0.0
  for i=1:length(solver.reg)
    solver.sᵏ[:] .+= solver.ρ[i]*adjoint(solver.regTrafo[i])*(solver.z[i] .- solver.zᵒˡᵈ[i]) # dual residual (concerning f(x))
    solver.ɛ_dt[:] .+= solver.ρ[i]*adjoint(solver.regTrafo[i])*solver.u[i]
  end

  # return the primal feasibilty measure as item and iteration number as state
  return solver.rᵏ, iteration+1
end

function converged(solver::ADMM)
  if norm(solver.sᵏ) >= solver.σᵃᵇˢ+solver.relTol*norm(solver.ɛ_dt)
    return false
  else
    for i=1:length(solver.reg)
      (solver.rᵏ[i] >= solver.ɛᵖʳⁱ[i]) && return false
    end
  end

  return true
end

@inline done(solver::ADMM,iteration::Int) = converged(solver) || iteration>=solver.iterations


text:
using PowerSystems
using PowerSimulations
using PowerSystemCaseBuilder

using Dates
using DataFrames

using Cbc # mip solver
solver = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 1, "ratioGap" => 0.5)
using Ipopt # solver that supports duals
ipopt_solver = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)

sys_DA = build_system(SIIPExampleSystems, "5_bus_matpower_DA")
sys_RT = build_system(SIIPExampleSystems, "5_bus_matpower_RT")

template_uc = template_unit_commitment()
template_ed = template_economic_dispatch()

problems = SimulationProblems(
    UC = OperationsProblem(
        template_uc,
        sys_DA,
        optimizer = solver,
        balance_slack_variables = true,
    ),
    ED = OperationsProblem(
        template_ed,
        sys_RT,
        optimizer = ipopt_solver,
        constraint_duals = [:CopperPlateBalance],
    ),
)

feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24))

feedforward = Dict(
    ("ED", :devices, :ThermalStandard) => SemiContinuousFF(
        binary_source_problem = ON,
        affected_variables = [ACTIVE_POWER],
    ),
)

#cache = Dict("UC" => [TimeStatusChange(ThermalStandard, PSI.ON)])
intervals = Dict("UC" => (Hour(24), Consecutive()), "ED" => (Hour(1), Consecutive()))

DA_RT_sequence = SimulationSequence(
    problems = problems,
    intervals = intervals,
    ini_cond_chronology = InterProblemChronology(),
    feedforward_chronologies = feedforward_chronologies,
    feedforward = feedforward,
)

file_path = mkpath(joinpath(".", "5-bus-simulation"))
sim = Simulation(
    name = "5bus-test",
    steps = 1,
    problems = problems,
    sequence = DA_RT_sequence,
    simulation_folder = file_path,
)

build!(sim)

execute!(sim, enable_progress_bar = false)

# Results

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
for N in [Float64, Rational{Int}, Float32]
    # random zero set
    rand(ZeroSet)

    Z = ZeroSet{N}(2)
    B = BallInf(ones(N, 2), N(1))

    # support vector
    d = N[1, 0]
    @test σ(d, Z) == zeros(N, 2)

    # testing that the zero set is neutral element for the Minkowski sum
    @test B ⊕ Z == B && Z ⊕ B == B

    cpa = MinkowskiSumArray([B, N(2) * B, N(3) * B])
    @test cpa ⊕ Z == cpa && Z ⊕ cpa == cpa

    # test M-sum of zero set with itself
    @test ZeroSet{N}(2) ⊕ ZeroSet{N}(2) == ZeroSet{N}(2)

    # boundedness
    @test isbounded(Z)

    # element & an_element function
    @test element(Z) ∈ Z
    @test element(Z, 1) == 0
    @test an_element(Z) ∈ Z

    # high and low
    @test high(Z) == zeros(N, 2)
    @test low(Z) == zeros(N, 2)

    # isempty
    @test !isempty(Z)

    # isuniversal
    answer, w = isuniversal(Z, true)
    @test !isuniversal(Z) && !answer && w ∉ Z

    # subset
    z = ZeroSet{N}(1)
    s1 = Singleton(N[0])
    s2 = Singleton(N[2])
    @test z ⊆ s1 && ⊆(z, s1, true)[1]
    subset, point = ⊆(z, s2, true)
    @test z ⊈ s2 && !subset && point ∈ z && point ∉ s2
    @test z ⊆ z && z ⊈ ZeroSet{N}(2)

    # linear map (concrete)
    M = hcat(N[1])
    Mz = linear_map(M, z)
    @test Mz isa ZeroSet && dim(Mz) == 1

    M = N[-1 -2;]
    MZ = linear_map(M, Z)
    @test MZ isa ZeroSet && dim(MZ) == 1

    # translation
    @test translate(Z, N[1, 2]) == Singleton(N[1, 2])

    # base type and element type
    @test basetype(Z) == basetype(typeof(Z)) == ZeroSet
    @test eltype(Z) == eltype(typeof(Z)) == N
end


text:
"""
    index_data(indices, args...)

Index into each arg with `indices`, placing `NaN`s where `indices` are zero.
"""
function index_data(indices, pts::AbstractVector{<:NTuple{N}}) where {N}
    return map(indices) do i
        i == 0 ? ntuple(_ -> NaN, N) : Float64.(pts[i])
    end
end
function index_data(indices, x::AbstractVector)
    return map(indices) do i
        i == 0 ? NaN : Float64(x[i])
    end
end
# images
function index_data(indices, x::AbstractMatrix)
    cart = CartesianIndices(x)
    return map(indices) do i
        i == 0 ? (NaN, NaN) : (cart[i][2], cart[i][1])
    end
end
function index_data(indices, x, y)
    return index_data(indices, x), index_data(indices, y)
end
function index_data(indices, x, y, z)
    return index_data(indices, x), index_data(indices, y), index_data(indices, z)
end

"""
    plottable(sx::AbstractCell, args...)
    plottable(sx::PersistenceInterval, args...)
    plottable(sx::Vector{<:AbstractCell}, args...)

Turn `sx` into a series that can be plotted. `args...`, should contain the data which tells
the plot where to place simplex vertices.
"""
function plottable(sx::AbstractCell, args...)
    return plottable([sx], args...)
end
function plottable(int::PersistenceInterval, args...)
    return plottable(simplex.(representative(int)), args...)
end
function plottable(rep::AbstractVector{<:AbstractChainElement}, args...)
    return plottable(simplex.(rep), args...)
end
function plottable(element::AbstractChainElement, args...)
    return plottable(simplex(element), args...)
end
function plottable(sxs::AbstractVector{<:AbstractCell{0}}, args...)
    indices = sort(only.(vertices.(sxs)))
    return index_data(indices, args...), [:seriestype => :scatter], 0
end
function plottable(sxs::AbstractVector{<:AbstractCell{1}}, args...)
    indices = mapreduce(vcat, vertices.(sxs)) do (u, v)
        [u, v, 0]
    end
    return index_data(indices, args...), [:seriestype => :path], 1
end
function plottable(sxs::AbstractVector{<:AbstractCell{D}}, args...) where {D}
    indices = mapreduce(vcat, vertices.(sxs)) do vs
        idxs = Int[]
        for (u, v, w) in subsets(vs, Val(3))
            append!(idxs, (u, v, w, u))
        end
        push!(idxs, 0)
        idxs
    end
    return index_data(indices, args...), [:seriestype => :path], D
end

function apply_threshold(sx::AbstractCell, thresh, thresh_strict)
    return birth(sx) ≤ thresh && birth(sx) < thresh_strict ? sx : nothing
end
function apply_threshold(sxs::AbstractVector{<:AbstractCell}, thresh, thresh_strict)
    return filter(sx -> birth(sx) ≤ thresh && birth(sx) < thresh_strict, sxs)
end
function apply_threshold(els::AbstractVector{<:AbstractChainElement}, thresh, thresh_strict)
    return apply_threshold(simplex.(els), thresh, thresh_strict)
end
function apply_threshold(int::PersistenceInterval, thresh, thresh_strict)
    if hasproperty(int, :representative)
        return apply_threshold(representative(int), thresh, thresh_strict)
    else
        error("interval has no representative. Run `ripserer` with `representatives=true`")
    end
end

@recipe function f(
    sx::Union{
        AbstractCell,
        AbstractChainElement,
        AbstractVector{<:AbstractCell},
        AbstractVector{<:AbstractChainElement},
        PersistenceInterval,
    },
    args...;
    threshold=Inf,
    threshold_strict=Inf,
)
    sx = apply_threshold(sx, threshold, threshold_strict)
    isnothing(sx) && return ()
    series, attrs, D = plottable(sx, args...)
    for (key, value) in attrs
        plotattributes[key] = get(plotattributes, key, value)
    end
    # splat colors and similar attributes over simplices
    for attr in keys(plotattributes)
        value = plotattributes[attr]
        if value isa Vector
            n = D < 2 ? D + 2 : 4 * binomial(D + 1, 3) + 1
            plotattributes[attr] = mapreduce(x -> fill(x, n), vcat, value)
        end
    end
    return series
end


text:
module FinEtoolsAcousticsTutorials

end # module


text:
#################### Flat Distribution ####################

struct Flat <: ContinuousUnivariateDistribution end

minimum(d::Flat) = -Inf
maximum(d::Flat) = Inf
insupport(d::Flat, x::Real) = true

logpdf(d::Flat, x::Real) = 0.0

#################### Null Distribution ####################

struct NullUnivariateDistribution <: UnivariateDistribution{ValueSupport} end

#################### SymUniform Distribution ####################

struct SymUniform <: ContinuousUnivariateDistribution
    SymUniform() = Uniform(-1.0, 1.0)
    SymUniform(μ::Real, σ::Real) = Uniform(μ - σ, μ + σ)
end

#################### Type Aliases ####################

const SymDistributionType = Union{
    Type{Biweight},
    Type{Cosine},
    Type{Epanechnikov},
    Type{Normal},
    Type{SymTriangularDist},
    Type{Triweight},
    Type{SymUniform},
}

const KernelDensityType = SymDistributionType

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
@testset "send_ack.jl" begin
    @testset begin
        isn = rand(UInt32(0):typemax(UInt32))
        sender = TCPSender(fixed_isn=WrappingInt32(isn))
        fill_window!(sender)
        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 1))
        write_bytes!(sender, "a")
        expect_seg(sender, no_flag=true, data="a")
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 1))
        expect_no_seg(sender)
    end

    @testset begin
        isn = rand(UInt32(0):typemax(UInt32))
        sender = TCPSender(fixed_isn=WrappingInt32(isn))
        fill_window!(sender)
        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 1))
        write_bytes!(sender, "a")
        expect_seg(sender, no_flag=true, data="a")
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 2))

        write_bytes!(sender, "b")
        expect_seg(sender, no_flag=true, data="b")
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 1))
        expect_no_seg(sender)
    end
    @testset begin
        isn = rand(UInt32(0):typemax(UInt32))
        sender = TCPSender(fixed_isn=WrappingInt32(isn))
        fill_window!(sender)
        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 1))
        write_bytes!(sender, "a")
        expect_seg(sender, no_flag=true, data="a")
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 2))

        write_bytes!(sender, "b")
        expect_seg(sender, no_flag=true, data="b")
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn + 1))
        expect_no_seg(sender)
    end

    @testset begin
        isn = rand(UInt32(0):typemax(UInt32))
        sender = TCPSender(fixed_isn=WrappingInt32(isn))
        fill_window!(sender)
        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        expect_state(sender, JLSponge.SYN_SENT)
        ack_received_test(sender, WrappingInt32(isn + 2), 1000)
        expect_state(sender, JLSponge.SYN_SENT)
    end

    @testset begin
        isn = rand(UInt32(0):typemax(UInt32))
        sender = TCPSender(fixed_isn=WrappingInt32(isn))
        fill_window!(sender)
        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        expect_state(sender, JLSponge.SYN_SENT)
        ack_received_test(sender, WrappingInt32(isn + 2), 1000)
        expect_state(sender, JLSponge.SYN_SENT)
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
using Sundials, Test

callback_f = function (du, u, p, t)
    du[1] = u[2]
    du[2] = -9.81
end

condtion = function (u, t, integrator) # Event when event_f(u,t,k) == 0
    u[1]
end

affect! = nothing
affect_neg! = function (integrator)
    integrator.u[2] = -integrator.u[2]
end

callback = ContinuousCallback(condtion, affect!, affect_neg!)

u0 = [50.0, 0.0]
tspan = (0.0, 15.0)
prob = ODEProblem(callback_f, u0, tspan)

sol = solve(prob, CVODE_Adams(), callback = callback)
@test sol(4.0)[1] > 0
sol = solve(prob, CVODE_BDF(), callback = callback)
@test sol(4.0)[1] > 0

condition = function (out, u, t, integrator)
    out[1] = u[1]
end

affect! = nothing
affect_neg! = function (integrator, idx)
    if idx == 1
        integrator.u[2] = -integrator.u[2]
    end
end

callback = VectorContinuousCallback(condition, affect!, affect_neg!, 1)
sol = solve(prob, CVODE_Adams(), callback = callback)
@test sol(4.0)[1] > 0
sol = solve(prob, CVODE_BDF(), callback = callback)
@test sol(4.0)[1] > 0

u0 = [1.0, 0.0]
function fun2(du, u, p, t)
    du[2] = -u[1]
    du[1] = u[2]
end
tspan = (0.0, 10.0)
prob = ODEProblem(fun2, u0, tspan)

function condition2(u, t, integrator)
    get_du(integrator)[1] > 0
end
affect2!(integrator) = terminate!(integrator)
cb = DiscreteCallback(condition2, affect2!)
sol = solve(prob, CVODE_BDF(), callback = cb)
@test sol.t[end] < 3.5

condition3(u, t, integrator) = u[2]
affect3!(integrator) = terminate!(integrator)
cb = ContinuousCallback(condition3, affect3!)
sol = solve(prob, CVODE_Adams(), callback = cb, abstol = 1e-12, reltol = 1e-12)
@test sol.t[end] ≈ pi


text:
const UIntPtr = Sys.WORD_SIZE == 64 ? UInt64 : UInt32

cookie = UInt32(0)

"""
    RemotePtr(ref::Ref{T})

Serializable pointer to a Julia object.

```
julia> c = RemotePtr(Ref(7))
julia> c[]
7
```

The Julia `serialize` function [stores pointers as zero](https://github.com/JuliaLang/julia/blob/08a233ddac2c0bcb07ca8a66792e8305a6124d3c/stdlib/Serialization/src/Serialization.jl#L727-L728).

RemotePtr converts a pointer to UInt64 so that it survives serialization.
A serialized RemotePtr can be stored in a file, or sent to a remote process.
A random cookie ensures that a pointer can only be dereferenced by
the process that created it.
"""
struct RemotePtr{T}
    p::UInt64
    cookie::UInt32
    function RemotePtr(o::Ref{T}) where T
        global cookie
        if cookie == 0
            cookie = rand(UInt32)
        end
        new{T}(convert(UInt64, pointer_from_objref(o)), cookie)
    end
end

function Base.getindex(rp::RemotePtr{T}) where T
    rp.cookie == cookie ||
        throw(ArgumentError("RemotePtr not valid on this node."))
    p = convert(Ptr{Ref{T}}, convert(UIntPtr, rp.p))
    ref = unsafe_pointer_to_objref(p)
    ref[]
end


text:
# this reproduces parsec's behaviour, by disallowing matched input to be
# used again.  to do this we need to:
# 1 - provide a source that allows input to be discarded
# 2 - discard input on success
# 3 - disable discarding input when inside Try()
# 4 - throw an exception when discarded input is accessed
# 5 - treat that exception as failure

# the source wraps an IO instance.  this is how julia manages files
# (which is presumably where this is needed most, since strings are
# already available in memory).  but strings can also be wrapped.

type TrySource{S}<:LineAt
    io::IO
    frozen::Int    # non-zero is frozen; count allows nested Try()
    zero::Int      # offset to lines (lines[x] contains line x+zero)
    right::Int     # rightmost expired column
    lines::Vector{S}
    TrySource(io::IO, line::S) = new(io, 0, 0, 0, S[line])
end

function TrySource(io::IO)
    line = readline(io)
    TrySource{typeof(line)}(io, line)
end

TrySource{S<:AbstractString}(s::S) = TrySource(IOBuffer(s))

function expire(s::TrySource, i::LineIter)
    if s.frozen == 0
        n = i.line - s.zero
        if n > 0
            s.lines = s.lines[n:end]
            s.zero += (n-1)
            if n > 1 || i.column > s.right
                s.right = i.column
            end
        end
    end
end

function line_at(f::TrySource, s::LineIter; check::Bool=true)
    if check
        if s.line <= f.zero || (s.line == f.zero+1 && s.column < f.right)
            throw(LineException())
        end
    end
    n = s.line - f.zero
    while length(f.lines) < n
        push!(f.lines, readline(f.io))
    end
    f.lines[n]
end

function next(f::TrySource, s::LineIter)
    # there's a subtlelty here.  the line is always correct for
    # reading more data (the check on done() comes *after* next).
    # this is so that getindex can access the line correctly if needed
    # (if we didn't have the line correct, getindex would take a slice
    # from the end of the previous line).
    line = line_at(f, s)
    c, col = next(line, s.column)
    if done(line, col)
        c, LineIter(s.line+1, 1)
    else
        c, LineIter(s.line, col)
    end
end

function done(s::TrySource, i::LineIter)
    line = line_at(s, i; check=false)
    done(line, i.column) && eof(s.io)
end

# the Try() matcher that enables backtracking

@auto_hash_equals type Try<:Delegate
    name::Symbol
    matcher::Matcher
    Try(matcher) = new(:Try, matcher)
end

@auto_hash_equals immutable TryState<:DelegateState
    state::State
end

execute(k::Config, m::Try, s::Clean, i) = error("use Try only with TrySource")

execute{S<:TrySource}(k::Config{S}, m::Try, s::Clean, i) = execute(k, m, TryState(CLEAN), i)

function execute{S<:TrySource}(k::Config{S}, m::Try, s::TryState, i)
    k.source.frozen += 1
    Execute(m, s, m.matcher, s.state, i)
end

function success{S<:TrySource}(k::Config{S}, m::Try, s::TryState, t, i, r::Value)
    k.source.frozen -= 1
    Success(TryState(t), i, r)
end

function failure{S<:TrySource}(k::Config{S}, m::Try, s::TryState)
    k.source.frozen -= 1
    FAILURE
end

function dispatch{S<:TrySource}(k::NoCache{S}, s::Success)
    (parent, parent_state) = pop!(k.stack)
    expire(k.source, s.iter)
    try
        success(k, parent, parent_state, s.child_state, s.iter, s.result)
    catch x
        isa(x, FailureException) ? FAILURE : rethrow()
    end
end

function dispatch{S<:TrySource}(k::Cache{S}, s::Success)
    parent, parent_state, key = pop!(k.stack)
    expire(k.source, s.iter)
    try
        k.cache[key] = s
    catch x
        isa(x, CacheException) ? nothing : rethrow()
    end
    try
        success(k, parent, parent_state, s.child_state, s.iter, s.result)
    catch x
        isa(x, FailureException) ? FAILURE : rethrow()
    end
end

# and a simple interface

parse_try(source, matcher; kargs...) = parse_one(TrySource(source), matcher; kargs...)
parse_try_dbg(source, matcher; kargs...) = parse_one_dbg(TrySource(source), matcher; kargs...)
parse_try_cache(source, matcher; kargs...) = parse_one_cache(TrySource(source), matcher; kargs...)
parse_try_cache_dbg(source, matcher; kargs...) = parse_one_cache_dbg(TrySource(source), matcher; kargs...)


text:
# Forwarding to Metric methods ## Does not seem to work for some reason?
# for Func in [:MetricPartials, :AutoMetricPartials, :ChristoffelSymbol, :ChristoffelPartials, :Riemann, :Ricci, :RicciScalar, :Weyl]
#     @eval ($Func(DM::AbstractDataModel, θ::AbstractVector{<:Number}; BigCalc::Bool=false, kwargs...) = $Func(FisherMetric(DM; kwargs...), θ; BicCalc=BigCalc))
# end

MetricPartials(DM::AbstractDataModel, θ::AbstractVector{<:Number}, Auto::Val=Val(:ForwardDiff); ADmode::Val=Auto, BigCalc::Bool=false, kwargs...) = MetricPartials(FisherMetric(DM; kwargs...), θ; ADmode=ADmode, BigCalc=BigCalc)
AutoMetricPartials(DM::AbstractDataModel, θ::AbstractVector{<:Number}, Auto::Val=Val(:ForwardDiff); ADmode::Val=Auto, BigCalc::Bool=false, kwargs...) = AutoMetricPartials(FisherMetric(DM; kwargs...), θ; ADmode=ADmode, BigCalc=BigCalc)
ChristoffelSymbol(DM::AbstractDataModel, θ::AbstractVector{<:Number}, Auto::Val=Val(:ForwardDiff); ADmode::Val=Auto, BigCalc::Bool=false, kwargs...) = ChristoffelSymbol(FisherMetric(DM; kwargs...), θ; ADmode=ADmode, BigCalc=BigCalc)
ChristoffelPartials(DM::AbstractDataModel, θ::AbstractVector{<:Number}, Auto::Val=Val(:ForwardDiff); ADmode::Val=Auto, BigCalc::Bool=false, kwargs...) = ChristoffelPartials(FisherMetric(DM; kwargs...), θ; ADmode=ADmode, BigCalc=BigCalc)
Riemann(DM::AbstractDataModel, θ::AbstractVector{<:Number}, Auto::Val=Val(:ForwardDiff); ADmode::Val=Auto, BigCalc::Bool=false, kwargs...) = Riemann(FisherMetric(DM; kwargs...), θ; ADmode=ADmode, BigCalc=BigCalc)
Ricci(DM::AbstractDataModel, θ::AbstractVector{<:Number}, Auto::Val=Val(:ForwardDiff); ADmode::Val=Auto, BigCalc::Bool=false, kwargs...) = Ricci(FisherMetric(DM; kwargs...), θ; ADmode=ADmode, BigCalc=BigCalc)
RicciScalar(DM::AbstractDataModel, θ::AbstractVector{<:Number}, Auto::Val=Val(:ForwardDiff); ADmode::Val=Auto, BigCalc::Bool=false, kwargs...) = RicciScalar(FisherMetric(DM; kwargs...), θ; ADmode=ADmode, BigCalc=BigCalc)
Weyl(DM::AbstractDataModel, θ::AbstractVector{<:Number}, Auto::Val=Val(:ForwardDiff); ADmode::Val=Auto, BigCalc::Bool=false, kwargs...) = Weyl(FisherMetric(DM; kwargs...), θ; ADmode=ADmode, BigCalc=BigCalc)

"""
Computes
```math
PDV[i,j,m] = (\\partial g)_{ij,m}
```
where the last index is the direction of the derivative.
"""
function MetricPartials(Metric::Function, θ::AbstractVector{<:Number}; ADmode::Val=Val(:ForwardDiff), kwargs...)
    PDV = Array{suff(θ), 3}(undef, length(θ), length(θ), length(θ))
    MetricPartials!(PDV, Metric, θ, ADmode; kwargs...);   PDV
end
function MetricPartials!(PDV::AbstractArray{<:Number,3}, Metric::Function, θ::AbstractVector{<:Number}, Auto::Val{false}; BigCalc::Bool=false)
    BigCalc && (θ = BigFloat.(θ))
    h = GetH(θ)
    for i in Base.OneTo(length(θ))
        PDV[:,:,i] = (1/(2*h))*(Metric(θ + h*BasisVector(i,length(θ))) - Metric(θ - h*BasisVector(i,length(θ))))
    end
    nothing
end

# For Float64, use AD to compute partial derivatives of metric, else use finite difference with BigFloat.
function AutoMetricPartials(Metric::Function, θ::AbstractVector{<:Number}; ADmode::Val=Val(:ForwardDiff), kwargs...)
    PDV = Array{suff(θ), 3}(undef, length(θ), length(θ), length(θ))
    MetricPartials!(PDV, Metric, θ, ADmode; kwargs...); PDV
end
function MetricPartials!(PDV::AbstractArray{<:Number,3}, Metric::Function, θ::AbstractVector{<:Number}, ADmode::Val; BigCalc::Bool=false)
    BigCalc && (θ = BigFloat.(θ))
    # PDV[:] = vec(GetJac(ADmode,Metric)(θ))
    # GetJac!(ADmode)(PDV, Metric, θ) ## no passthrough
    GetMatrixJac!(ADmode,Metric)(PDV, θ)
    nothing
end

# function ChristoffelSymbol(Metric::Function, θ::AbstractVector{<:Number}; BigCalc::Bool=false)
#     Finv = inv(Metric(θ));    FPDV = MetricPartials(Metric, θ; BigCalc=BigCalc)
#     if (suff(θ) == Float64) && BigCalc     FPDV = convert(Array{Float64,3}, FPDV)    end
#     @tullio Christoffels[a,i,j] := ((1/2) * Finv)[a,m] * (FPDV[j,m,i] + FPDV[m,i,j] - FPDV[i,j,m])
# end

# Accuracy ≈ 3e-11
# BigCalc for using BigFloat Calculation in finite differencing step but outputting Float64 again.
"""
    ChristoffelSymbol(DM::DataModel, θ::AbstractVector; BigCalc::Bool=false)
    ChristoffelSymbol(Metric::Function, θ::AbstractVector; BigCalc::Bool=false)
Calculates the components of the ``(1,2)`` Christoffel symbol ``\\Gamma`` at a point ``\\theta`` (i.e. the Christoffel symbol "of the second kind") through finite differencing of the `Metric`. Accurate to ≈ 3e-11.
`BigCalc=true` increases accuracy through `BigFloat` calculation.
"""
function ChristoffelSymbol(Metric::Function, θ::AbstractVector{<:Number}; ADmode::Val=Val(:ForwardDiff), kwargs...)
    # Need to initialize with zeros() instead of Array{}() due to tensor macro
    Γ = zeros(suff(θ), length(θ), length(θ), length(θ))
    ChristoffelSymbol!(Γ, Metric, θ, ADmode; kwargs...);  Γ
end
function ChristoffelSymbol!(Γ::AbstractArray{<:Number,3}, Metric::Function, θ::AbstractVector{<:Number}, ADmode::Val; kwargs...)
    PDV = zeros(suff(θ), length(θ), length(θ), length(θ))
    ChristoffelSymbol!!(PDV, Γ, Metric, θ, ADmode; kwargs...)
end
function ChristoffelSymbol!!(PDV::AbstractArray{<:Number,3}, Γ::AbstractArray{<:Number,3}, Metric::Function, θ::AbstractVector{<:Number}, ADmode::Val; kwargs...)
    MetricPartials!(PDV, Metric, θ, ADmode; kwargs...)
    # ((suff(θ) == Float64) && BigCalc) && (FPDV = convert(Array{Float64,3}, FPDV))
    PDV2Christoffel!(Γ, inv(Metric(θ)), PDV)
    nothing
end

function PDV2Christoffel!(Γ::AbstractArray{<:Number,3}, InvMetric::AbstractMatrix{<:Number}, PDV::AbstractArray{<:Number,3})
    @tullio Γ[a,i,j] = ((1/2) * InvMetric)[a,m] * (PDV[j,m,i] + PDV[m,i,j] - PDV[i,j,m])
    nothing
end
function PDV2Christoffel(InvMetric::AbstractMatrix{<:Number}, PDV::AbstractArray{<:Number,3})
    n = size(InvMetric)[1];    Γ = Array{suff(PDV), 3}(undef, n, n, n)
    PDV2Christoffel!(Γ, InvMetric, PDV);    Γ
end

function PDV2Christoffel2!(Γ::AbstractArray{<:Number,3}, InvMetric::AbstractMatrix{<:Number}, PDV::AbstractArray{<:Number,3})
    # Use symmetry in lower indices of Christoffel symbols due to Levi--Civita connection to eliminate one term in contraction
    @tullio Γ[a,i,j] = InvMetric[a,m] * (PDV[m,j,i] - (0.5*PDV)[i,j,m])
    nothing
end

ChristoffelTerm(Γ::AbstractArray{<:Number,3}, v::AbstractVector{<:Number}) = @tullio Res[a] := (-1*Γ)[a,b,c] * v[b] * v[c]

# function ChristoffelPartials(Metric::Function, θ::AbstractVector{<:Number}; BigCalc::Bool=false)
#     BigCalc && (θ = BigFloat.(θ))
#     DownUpDownDown = Array{suff(θ)}(undef,length(θ),length(θ),length(θ),length(θ))
#     h = GetH(θ)
#     for i in 1:length(θ)
#         DownUpDownDown[i,:,:,:] .= (ChristoffelSymbol(Metric,θ + h*BasisVector(i,length(θ))) .- ChristoffelSymbol(Metric,θ - h*BasisVector(i,length(θ))))
#     end;        (1/(2*h))*DownUpDownDown
# end

function ChristoffelPartials(Metric::Function, θ::AbstractVector{<:Number}; ADmode::Val=Val(:ForwardDiff), kwargs...)
    ∂Γ = Array{suff(θ), 4}(undef, length(θ), length(θ), length(θ), length(θ))
    ChristoffelPartials!(∂Γ, Metric, θ, ADmode; kwargs...);   ∂Γ
end

"""
Returns partial derivatives of Christoffel Symbols Γ with index structure up-down-down-down, i.e. (∂Γ)ᵇₑⱼₐ = Γᵇₑⱼ,ₐ where the last index denotes the direction of the derivative.
"""
function ChristoffelPartials!(∂Γ::AbstractArray{<:Number,4}, Metric::Function, θ::AbstractVector{<:Number}, ADmode::Val{false}; BigCalc::Bool=false)
    BigCalc && (θ = BigFloat.(θ))
    h = GetH(θ)
    FPDV = Array{suff(θ),3}(undef, length(θ), length(θ), length(θ))
    Γ₁ = Array{suff(θ),3}(undef, length(θ), length(θ), length(θ));  Γ₂ = Array{suff(θ),3}(undef, length(θ), length(θ), length(θ))
    for i in Base.OneTo(length(θ))
        # ∂Γ[i,:,:,:] .= (1/(2*h))*(ChristoffelSymbol(Metric,θ + h*BasisVector(i,length(θ))) .- ChristoffelSymbol(Metric,θ - h*BasisVector(i,length(θ))))
        ChristoffelSymbol!!(FPDV, Γ₂, Metric, θ + h*BasisVector(i,length(θ)), ADmode; BigCalc=BigCalc)
        ChristoffelSymbol!!(FPDV, Γ₁, Metric, θ - h*BasisVector(i,length(θ)), ADmode; BigCalc=BigCalc)
        ∂Γ[:,:,:,i] = (1/(2*h))*(Γ₂ - Γ₁)
    end;    nothing
end
function ChristoffelPartials!(∂Γ::AbstractArray{<:Number,4}, Metric::Function, θ::AbstractVector{<:Number}, ADmode::Val; BigCalc::Bool=false)
    BigCalc && (θ = BigFloat.(θ))
    GetMatrixJac!(ADmode, x->ChristoffelSymbol(Metric, x; ADmode=ADmode))(∂Γ, θ)
    nothing
end

"""
    Riemann(DM::DataModel, θ::AbstractVector; BigCalc::Bool=false)
    Riemann(Metric::Function, θ::AbstractVector; BigCalc::Bool=false)
Calculates the components of the ``(1,3)`` Riemann tensor by finite differencing of the `Metric`. `BigCalc=true` increases accuracy through BigFloat calculation.
"""
function Riemann(Metric::Function, θ::AbstractVector{<:Number}; ADmode::Val=Val(:ForwardDiff), BigCalc::Bool=false)
    ∂Γ = ChristoffelPartials(Metric, θ; ADmode=ADmode, BigCalc=BigCalc)
    if (suff(θ) == Float64) && BigCalc
        ∂Γ = convert(Array{Float64,4}, ∂Γ)
    end
    Γ = ChristoffelSymbol(Metric, θ; ADmode=ADmode, BigCalc=BigCalc)
    # @tullio Riem[m,i,k,p] := DownUpDownDown[k,m,i,p] - DownUpDownDown[p,m,i,k] + Gamma[a,i,p]*Gamma[m,a,k] - Gamma[a,i,k]*Gamma[m,a,p]
    Riem = Array{suff(θ), 4}(undef, length(θ), length(θ), length(θ), length(θ))
    RiemannLastInd!(Riem, ∂Γ, Γ);  Riem
end

# """
# Assumes index structure (∂Γ)ₐᵇₑⱼ where the FIRST index denotes direction of derivative.
# ```math
# \\tensor{\\mathrm{Riem}}{^i _j _k _l} = \\pdv{x^k}\\, \\tensor{\\Gamma}{^i _j _l} - \\pdv{x^l}\\, \\tensor{\\Gamma}{^i _j _k} + \\tensor{\\Gamma}{^i _a _k} \\, \\tensor{\\Gamma}{^a _j _l} - \\tensor{\\Gamma}{^i _a _l} \\, \\tensor{\\Gamma}{^a _j _k}
# ```
# """
# function RiemannFirstInd!(Riem::AbstractArray{<:Number,4}, ∂Γ::AbstractArray{<:Number,4}, Γ::AbstractArray{<:Number,3})
#     # @tensor Riem[i,j,k,l] = ∂Γ[k,i,j,l] - ∂Γ[l,i,j,k] + Γ[i,a,k]*Γ[a,j,l] - Γ[i,a,l]*Γ[a,j,k]
#     # With Tullio, this needs to be split over two lines!!! (Why?!)
#     @tullio Riem[i,j,k,l] = ∂Γ[k,i,j,l] - ∂Γ[l,i,j,k]
#     @tullio Riem[i,j,k,l] += Γ[i,a,k]*Γ[a,j,l] - Γ[i,a,l]*Γ[a,j,k]
#     nothing
# end
"""
Assumes index structure (∂Γ)ᵇₑⱼₐ = Γᵇₑⱼ,ₐ where the LAST index denotes direction of derivative.
```math
\\tensor{\\mathrm{Riem}}{^i _j _k _l} = \\pdv{x^k}\\, \\tensor{\\Gamma}{^i _j _l} - \\pdv{x^l}\\, \\tensor{\\Gamma}{^i _j _k} + \\tensor{\\Gamma}{^i _a _k} \\, \\tensor{\\Gamma}{^a _j _l} - \\tensor{\\Gamma}{^i _a _l} \\, \\tensor{\\Gamma}{^a _j _k}
```
"""
function RiemannLastInd!(Riem::AbstractArray{<:Number,4}, ∂Γ::AbstractArray{<:Number,4}, Γ::AbstractArray{<:Number,3})
    # @tensor Riem[i,j,k,l] = ∂Γ[i,j,l,k] - ∂Γ[i,j,k,l] + Γ[i,a,k]*Γ[a,j,l] - Γ[i,a,l]*Γ[a,j,k]
    # With Tullio, this needs to be split over two lines!!! (Why?!)
    @tullio Riem[i,j,k,l] = ∂Γ[i,j,l,k] - ∂Γ[i,j,k,l]
    @tullio Riem[i,j,k,l] += Γ[i,a,k]*Γ[a,j,l] - Γ[i,a,l]*Γ[a,j,k]
    nothing
end

"""
    Ricci(DM::DataModel, θ::AbstractVector; BigCalc::Bool=false)
    Ricci(Metric::Function, θ::AbstractVector; BigCalc::Bool=false)
Calculates the components of the ``(0,2)`` Ricci tensor by finite differencing of the `Metric`. `BigCalc=true` increases accuracy through `BigFloat` calculation.
"""
function Ricci(Metric::Function, θ::AbstractVector{<:Number}; kwargs...)
    Riem = Riemann(Metric, θ; kwargs...)
    # For some reason, it is necessary to prefill here.
    RIC = Array{suff(θ),2}(undef, length(θ), length(θ))
    @tullio RIC[a,b] = Riem[c,a,c,b]
end

"""
    RicciScalar(DM::DataModel, θ::AbstractVector; BigCalc::Bool=false) -> Real
    RicciScalar(Metric::Function, θ::AbstractVector; BigCalc::Bool=false) -> Real
Calculates the Ricci scalar by finite differencing of the `Metric`. `BigCalc=true` increases accuracy through `BigFloat` calculation.
"""
RicciScalar(Metric::Function, θ::AbstractVector{<:Number}; kwargs...) = tr(transpose(Ricci(Metric, θ; kwargs...)) * inv(Metric(θ)))

"""
(0,4) Weyl curvature tensor. NEEDS TESTING.
"""
function Weyl(Metric::Function, θ::AbstractVector{<:Number}; kwargs...)
    length(θ) < 4 && return zeros(length(θ),length(θ),length(θ),length(θ))
    Riem = Riemann(Metric, θ; kwargs...)
    g = BigCalc ? Metric(BigFloat.(θ)) : Metric(θ)
    @tullio Ric[a,b] := Riem[m,a,m,b]
    @tullio PartA[i,k,l,m] := Ric[i,m]*g[k,l] - Ric[i,l] * g[k,m] + Ric[k,l] * g[i,m] - Ric[k,m] * g[i,l]
    @tullio PartB[i,k,l,m] := g[i,l] * g[k,m] - g[i,m] * g[k,l]
    @tullio Rlow[a,b,c,d] := g[a,e] * Riem[e,b,c,d]
    Rlow .+ (length(θ) - 2)^(-1) .* PartA .+ ((length(θ) - 1)^(-1) * (length(θ) - 2)^(-1) * tr(inv(g) * Ric)) .* PartB
end