text:
"""
A model for determining the memory of agents traversing a network given a list
of observed paths. For example, it can answer the question, "how many webpages
back do we have to look to predict where someone surfing the web will go next?"

Based on "When is a Network a Network? Multi-Order Graphical Model Selection in
Pathways and Temporal Networks" by Ingo Scholtes. The original paper may be
found at https://arxiv.org/pdf/1702.05499.pdf.

Also of note is Scholtes' argument that traditional network analytic methods are
only justified if the optimal memory for modeling paths is 1, i.e. the Markov
property holds.
"""

using MatrixNetworks # Used for realizing the network implied by observed paths
using SparseArrays # Used for going from MatrixNetwork to a matrix
using StatsFuns # Used for the χ² CDF
using Dates # Used for runtime calculations

"""
`HO_Markov`
=========

Represents a higher-order Markov model on a set of observed paths. (pp. 2)

Fields
------
    - `k::Int64`: The order of the model, which denotes the length of the paths
    - `P::Dict{Vector{Int64}, Float64}`: The log probabilities that particular paths
                                         will be generated. More specifically,
                                         given the first `k-1` steps of the path,
                                         the log probability that step `k` occurs

While traditional Markov models are memoryless, in that where they go next is
based solely on their current location, "higher-order" Markov models determine
the probability of each possible next location based on the `k` previous locations.
"""
struct HO_Markov
    k::Int64
    P::Dict{Vector{Int64}, Float64}
end

"""
`MultiOrderModel`
=================

Represents a multi-order model consisting of many layers of `HO_Markov` models.
(pp. 4)

Fields
------
    - `K::Int64`: The order of the highest `HO_Markov` model
    - `layers::Vector{HO_Markov}`: Array of `HO_Markov` models, going from orders
                                   0 to K
    - `P::Dict{Vector{Int64}, Float64}`: The cumulative log probabilities for the
                                         realization of each path based on the
                                         layers. For example, P(a,b,c) = P(a)P(b|a)P(c|a,b)
"""
struct MultiOrderModel
    K::Int64
    layers::Vector{HO_Markov}
    P::Dict{Vector{Int64}, Float64}
end

"""
`optimal_HO_Markov`
=================

Given a multiset of observed paths `S`, constructs the `HO_Markov` that is most
likely to generate `S`. (pp. 3)

Arguments
---------
    - `S::Vector{Vector{Int64}}`: Observed paths
    - `k::Int64`: Order of the desired model

Example
-------
~~~~
optimal_HO_Markov([[1,2], [1,3], [1,2,3], [2,3], [3,2], [3,1,2]], 1)
~~~~
"""
function optimal_HO_Markov(S::Vector{Vector{Int64}}, k::Int64)
    P = Dict{Vector{Int64}, Float64}() # Log probabilities
    p_counts = Dict{Vector{Int64}, Int64}() # Number of times a path of length `k` appears in S
    w_counts = Dict{Vector{Int64}, Int64}() # Number of times a path of length `k-1` appears in S
    S_k = []
    for q in S
        for i = 1:size(q,1)-k # Get all k-subpaths of q
            p_k = q[i:i+k] # k-subpath
            w_k = p_k[1:k] # p_k excluding the last node

            push!(S_k, p_k)

            if !haskey(p_counts, p_k)
                p_counts[p_k] = 0
            end
            p_counts[p_k] += 1

            if !haskey(w_counts, w_k)
                w_counts[w_k] = 0
            end
            w_counts[w_k] += 1
        end
    end

    for p in Set(S_k)
        P[p] = log(p_counts[p] / w_counts[p[1:k]]) # Pr = proportion of paths identical to p[1:k] that end in p[k]
                                                   # E.g., P(a,b,c) = |{x = (a,b,c) for x in S}| / |{x = (a,b,d) for x in S for all observed nodes d}|
    end
    return HO_Markov(k, P)
end

"""
`optimal_multiorder`
====================

Generates the optimal multi-order model of maximum order `K` for a set of
observed paths. (pp. 4)

Arguments
---------
    - `S::Vector{Vector{Int64}}`: The observed paths
    - `K::Int64`: The maximum order of the model
"""
function optimal_multiorder(S::Vector{Vector{Int64}}, K::Int64)
    t1 = now()
    layers = [optimal_HO_Markov(S, k-1) for k = 1:K+1] # Generate each layer
    println("Layer models generated for K=$K. Took $(now() - t1).")
    P = Dict{Vector{Int64}, Float64}()

    for p in S
        P[p] = 0.0
        for k = 1:min(K,size(p,1)) # Multiply by probability in each layer. E.g., P(a,b,c) = P(a) * P(b|a) * P(c|a,b)
            P[p] += layers[k].P[p[1:k]]
        end
        for i = K+1:size(p,1) # Multiply by M_K probability to produce path. For paths longer than K
            P[p] += layers[K+1].P[p[i-K:i]]
        end
    end

    return MultiOrderModel(K, layers, P)
end

"""
`increment_model`
=================

Generate a multi-order model of one order higher than the one given.

Arguments
---------
    - `model::MultiOrderModel`: The model to increment

This method takes advantage of the nestedness of multi-order models: that a model
of order K+1 contains all the layers of a model of order K.
"""
function increment_model(model::MultiOrderModel)
    S = collect(keys(model.P))
    K = model.K + 1
    layers = [model.layers; [optimal_HO_Markov(S, K)]]
    P = Dict{Vector{Int64}, Float64}()

    for p in S
        P[p] = 0.0
        for k = 1:min(K,size(p,1)) # Multiply by probability in each layer. E.g., P(a,b,c) = P(a)P(b|a)P(c|a,b)
            P[p] += layers[k].P[p[1:k]]
        end
        for i = K+1:size(p,1) # Multiply by M_K probability to produce path. For paths longer than K
            P[p] += layers[K+1].P[p[i-K:i]]
        end
    end

    return MultiOrderModel(model.K+1, layers, P)
end

"""
`multiorder_likelihood`
=======================

Compute the likelihood that a multi-order model `M_k` produces the observed paths `S`
by taking the product of all path probabilities. (pp. 4)

Precondition
------------
    - S ⊆ M_k.P.keys
"""
function multiorder_likelihood(M_k::MultiOrderModel, S::Vector{Vector{Int64}})
    @assert issubset(Set(S), keys(M_k.P))
    return sum([M_k.P[p] for p in S])
end

"""
`K_opt`
=======

Compute the memory length `K` by which the observed paths are optimally modeled.
(pp. 4,5)

Arguments
---------
    - `S::Vector{Vector{Int64}}`: The observed paths
    - `A::MatrixNetwork`: The graph underlying `S`
    - `p::Float64 (=0.01)`: The significance threshold that must be crossed to
                            prefer a simpler model to a more complex one

Returns
-------
    - The optimal value for `K`
    - The corresponding multi-order model of order `K`
"""
function K_opt(S::Vector{Vector{Int64}}, A::MatrixNetwork; p::Float64=0.01)
    K_opt = 0
    t1 = now()
    M_k = optimal_multiorder(S, K_opt)

    println("Initial K=0 model generated. Took $(now() - t1).")

    ### Compute the degrees of freedom in a binary adjacency matrix ###
    function matrix_df(A::Array{Bool,2}, k::Int64)
        Ak = A^k
        return sum(Ak) - sum([1 for i = 1:size(A,1) if any(Ak[i,:] .!= 0)])
    end

    ### Compute the degrees of freedom of a multi-order model ###
    function df(M::MultiOrderModel)
        return (A.n - 1) + sum([matrix_df(Array(sparse(A)), k) for k = 1:M.K])
    end

    M_kl = multiorder_likelihood(M_k, S)
    dfMk = df(M_k)

    ### Calculate the p-value of M_k over M_K###
    function p_value(M_K::MultiOrderModel) # pp. 5
        M_Kl = multiorder_likelihood(M_K, S)
        dfMK = df(M_K)

        return chisqcdf(dfMK-dfMk, 2(M_Kl-M_kl)), M_Kl, dfMK
    end

    while true # Keep increasing maximum order until the significance threshold
        t2 = now()
        M_K = increment_model(M_k)
        println("Higher K=$(K_opt+1) model generated. Took $(now() - t1).")
        t3 = now()
        cp, M_kl, dfMk = p_value(M_K)
        println("Comparative p-value between K=$K_opt and K=$(K_opt+1): $cp. Took $(now() - t3).\n")
        if cp < p
            break
        else
            M_k = deepcopy(M_K)
            K_opt += 1
        end
    end
    println("Done. K_opt=$K_opt. Total analysis took $(now() - t1).")
    return K_opt, M_k
end

"""
`read_paths`
============

Read a collection of paths from a file.

Each line must contain a single path with an optional frequency at the end to
indicate the number of times to count the path.

Arguments
---------
    - `filepath::String`: Where to read paths from
    - `separator::String (=",")`: The string that separates nodes on each line
    - `frequency::Bool (=false)`: Whether the file has frequencies at the end of
                                  each line

Returns
-------
    - An array of paths, where each path is an array of integers, in which each
      node from the file has been assigned a unique number
    - The `MatrixNetwork` implied by the paths
    - The dictionary relating node strings to node numbers

Example
-------
~~~~
read_paths("data.txt"; separator=" ", frequency=true)
~~~~
"""
function read_paths(filepath::String; separator::String=",", frequency::Bool=false)
    println("Begin reading $filepath.")
    nodes = Dict() # Maps strings in the file to positive integers
    paths = Array{Array{Int64,1},1}()
    edges = Set{Tuple{Int64,Int64}}() # Edges observed in paths
    t1 = now()
    open(filepath) do file
        lines = [ln for ln in eachline(file)]
        for ln in lines
            separated = split(ln, separator)
            p = []
            for i = 1:(frequency ? length(separated)-1 : length(separated)) # Parse strings in line
                n = separated[i]
                if !haskey(nodes,n)
                    nodes[n] = length(nodes) + 1
                end
                if i < length(separated)-1 && !haskey(nodes,separated[i+1]) # Get following string
                    nodes[separated[i+1]] = length(nodes) + 1
                end
                push!(p, nodes[n])
                if i < length(separated)-1
                    push!(edges, (nodes[n], nodes[separated[i+1]])) # Add edge between current and next
                end
            end
            if frequency
                for i = 1:Int(parse(Float64, separated[length(separated)]))
                    push!(paths, p)
                end
            else
                push!(paths, p)
            end
        end
    end
    println("$(length(paths)) paths read, implying a network of $(length(nodes)) nodes and $(length(edges)) edges. Took $(now() - t1).")
    return paths, MatrixNetwork(collect(edges), length(nodes)), nodes
end


text:
#
# This file is part of the Guards.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

using Guards, Test, Distributed

length(procs()) == 1 && addprocs(1)

@everywhere using Actors, Guards

@everywhere did(x) = x isa Guard ? (@grd x) : x

act1 = Actors.spawn(Bhv(did), pid=2)

gd = guard([1,2,3])

@test call(act1, gd) == [1,2,3]
@grd push!(gd, 4)
@test call(act1, gd) == [1,2,3,4]

ext: jl


text:
using StaticCompiler
using Test
using LLVM
using Libdl

cd(@__DIR__)
@testset "ccalls" begin
    include("ccalls.jl")
end

@testset "globals" begin
    include("globals.jl")
end

@testset "others" begin
    include("others.jl")
end

@testset "standalone" begin
    include("standalone-exe.jl")
end


text:
module PuzzleGame

using Images, ImageView, TestImages
using Gtk, Reactive, GtkReactive
using Random

export main

function rules()
    println("please swap tiles to complete the image")
end

function init_game(img::Array{T,2}; nof_row_blocks::Int=4, nof_col_blocks::Int=4 ) where T  <: Colorant
    global block_size
    block_size = Int.(size(img) ./ (nof_row_blocks, nof_col_blocks))
    row_slices = [ i*block_size[1] + 1 : (i+1) * block_size[1] for i in 0:nof_row_blocks-1]
    col_slices = [ i*block_size[2] + 1 : (i+1) * block_size[2] for i in 0:nof_col_blocks-1]

    blocks = [(x,y) for y in col_slices for x in row_slices]

    shuf = zeros(eltype(img), size(img))
    shuf_blocks = shuffle(blocks)

    for b in 1:length(blocks)
        shuf[shuf_blocks[b][1], shuf_blocks[b][2]] = img[blocks[b][1], blocks[b][2]]
    end
    shuf
end

"""
get the margins of the block from (x,y) pixel location
"""
function get_block_coords(coord::Tuple{AbstractFloat, AbstractFloat}, block_size::Tuple{Int,Int})
    offset = rem.(coord, block_size)
    left = coord[1] - offset[1] + 1
    top = coord[2] - offset[2] + 1

    right = left + block_size[1] - 1
    bottom = top + block_size[2] - 1

    margins = Int.(round.([left, top, right, bottom]))
    return margins
end

"""
box = Int.([left, top, right, bottom])
"""
function swap_blocks(boxA, boxB, img)
    tmp = img[boxA[2]:boxA[4], boxA[1]:boxA[3]]
    img[boxA[2]:boxA[4], boxA[1]:boxA[3]] = img[boxB[2]:boxB[4], boxB[1]:boxB[3]]
    img[boxB[2]:boxB[4], boxB[1]:boxB[3]] = tmp
    return img
end

function main()
    img = testimage("fabio")
    new = init_game(img)
    guidict = imshow(new)
    rules()

    @guarded signal_connect(guidict["gui"]["window"], :destroy) do widget
        isAcitve = false
        println("you closded the window")
    end

    coord_sig = @guarded map(guidict["gui"]["canvas"].mouse.buttonpress, init=(0.0,0.0)) do btn
        global ann_handler, firstClick
        row = Float64(btn.position.x)
        col = Float64(btn.position.y)
        println("you clicked @ [x=$(row), y=$(col)] with $(btn.button)")
        margins = get_block_coords((row, col), block_size)

        if firstClick
            println("highlight the selected block")
            ann_handler = annotate!(guidict, AnnotationBox(margins[1], margins[2], margins[3], margins[4],
                                    linewidth=2.5, color=colorant"blue"))
            firstClick = false

        else
            println("swap the blocks")
            prev_margins = Int.([ann_handler.ann.data.left, ann_handler.ann.data.top, ann_handler.ann.data.right, ann_handler.ann.data.bottom])
            swap_blocks(margins, prev_margins, new)
            delete!(guidict, ann_handler)
            firstClick = true
        end

        return row,col
    end

    return guidict
end

#Top level code
firstClick = true
isActive = true
if !isinteractive()
    main()
    while isAcitve
        sleep(1)
    end
    println("code done..")
end

end # module PuzzleGame


text:
# ## Axis
# ### axis

#md # ![FILE_NAME.png](images/FILE_NAME.png)

using Measures
gr()

plot(sin, 0, 2π,
    xticks=0:0.5:2π,
    xrotation=60,
    xtickfont=font(30),
    bottom_margin=15mm)


text:
macro FZ_FR()#(N, F_yf, F_yr, V, Ax, R)
	code=quote
		# define the lateral tire forces
		F_yf = @F_YF()
		F_yr = @F_YR()
		FZ_fr = zeros(length(Ax),1)
		for ii in eachindex(Ax)
			FZ_fr[ii] = 0.5*(FzF0 - KZX*(Ax[ii] - V[ii]*R[ii])) + KZYF*((F_yf[ii] + F_yr[ii])/m)
		end
		FZ_fr
	end
	return esc(code)
end


text:
"""
Measure the roles of the species in a network (given as a `Modular` object),
according to two parameters. *z* is the z-score of the degree within a modules
(which species connect *within* their modules). *c* is the degree across
networks (which species connect different networks).

The result is given as a `data.frame` object.
"""
function network_roles(M::Modular)
   roles = zeros(Float64, (sum(size(M.A)), 2))
   # Step 1 -- measure of z
   for m in 1:size(M.S)[2]
      blues = M.S[collect(1:size(M.A)[2]),m]
      reds = M.S[[size(M.A)[2]+[1:size(M.A)[1]]],m]
      mod_only = M.A[vec(reds.==1),vec(blues.==1)]
      degree_seq = append!(vec(sum(mod_only, 1)), vec(sum(mod_only, 2)))
      if std(degree_seq) == 0
         # If all species have equal degree within a module
         # the z-score is 0
         z = (degree_seq .- mean(degree_seq))
      else
         z = (degree_seq .- mean(degree_seq))./std(degree_seq)
      end
      idx = [1:size(M.S)[1]][append!(vec(blues.==1),vec(reds.==1))]
      roles[idx,1] = z
   end
   # Step 2 -- measure of c
   nr, nc = size(M.A)
   for col = 1:nc
      k = sum(M.A, 1)[col]
      km = zeros(Int64, size(M.S)[2])
      for row = 1:nr
         if M.A[row,col] == 1
            km[[1:(size(M.S)[2])][vec(M.S[(nc+row),:] .== 1)]] += 1
         end
      end
      c = 1-sum((km./k).^2)
      roles[col,2] = c
   end
   for row = 1:nr
      k = sum(M.A, 2)[row]
      km = zeros(Int64, size(M.S)[2])
      for col = 1:nc
         if M.A[row,col] == 1
            km[[1:(size(M.S)[2])][vec(M.S[col,:] .== 1)]] += 1
         end
      end
      c = 1-sum((km./k).^2)
      roles[(nc+row),2] = c
   end
   roles_df = DataFrame(id = append!(collect(1:nc), collect(1:nr)),
                        margin = append!(rep("col", nc), rep("row", nr)),
                        z = roles[:,1],
                        c = roles[:,2]
                        )
   return roles_df
end


text:
#----------------------------------------------------------------------------------------------#
#                                 Plain Abstract Type Factory                                  #
#----------------------------------------------------------------------------------------------#

"""
`function tyArchy(t::Union{DataType,UnionAll})`\n
Returns a string suitable for documenting the hierarchy of an abstract type.
"""
function tyArchy(t::Union{DataType,UnionAll})
    h = Any[t]; while h[end] != Any; append!(h, [supertype(h[end])]); end
    H = Tuple(string(nameof(i)) for i in h)
    join(H, " <: ")
end

"""
`function mkNonPAbs(TY::Symbol, TP::Symbol, what::AbstractString, xp::Bool=true)`\n
Declares exactly one new, non-parametric, abstract type `TY <: TP`. Argument `what` is inserted
in the new type documentation, and `xp` controls whether or not the new abstract type is
exported (default `true`).
"""
function mkNonPAbs(TY::Symbol, TP::Symbol, what::AbstractString, xp::Bool=true)
    if !(eval(TP) isa DataType)
        error("Type parent must be a DataType. Got $(string(TP)).")
    end
    hiStr = tyArchy(eval(TP))
    dcStr = """
`abstract type $(TY) <: $(TP) end`\n
Abstract supertype for $(what).\n
## Hierarchy\n
`$(TY) <: $(hiStr)`
    """
    @eval begin
        # Abstract type definition
        abstract type $TY <: $TP end
        # Type documentation
        @doc $dcStr $TY
        # Type exporting
        if $(xp); export $TY; end
    end
end

#----------------------------------------------------------------------------------------------#
#                        EngTherm Plain Abstract Supertype Definitions                         #
#----------------------------------------------------------------------------------------------#

# EngTherm root abstract type
mkNonPAbs(:AbstractTherm   , :Any          , "thermodynamic entities"                         )

# BASE branch
mkNonPAbs(  :BASES         , :AbstractTherm, "quantity bases"                                 )
mkNonPAbs(    :ThermBase   , :BASES        , "thermodynamic bases"                            )
mkNonPAbs(      :IntBase   , :ThermBase    , "intensive bases"                                )
mkNonPAbs(        :MA      , :IntBase      , "the MAss base"                                  )
mkNonPAbs(        :MO      , :IntBase      , "the MOlar base"                                 )
mkNonPAbs(      :ExtBase   , :ThermBase    , "non-intensive bases"                            )
mkNonPAbs(        :SY      , :ExtBase      , "the SYstem (extensive) base"                    )
mkNonPAbs(        :DT      , :ExtBase      , "the Time Derivative (rate) base"                )
mkNonPAbs(    :ExactBase   , :BASES        , "type-exactness bases"                           )
mkNonPAbs(      :EX        , :ExactBase    , "the EXact base"                                 )
mkNonPAbs(      :MM        , :ExactBase    , "the MeasureMent base"                           )

# AUX branch
mkNonPAbs(  :AUX           , :AbstractTherm, "ancillary EngTherm types"                       )
mkNonPAbs(    :AuxFunc     , :AUX          , "ancillary functions"                            )

#----------------------------------------------------------------------------------------------#
#                Concrete type unions for allowed abstract supertype parameters                #
#----------------------------------------------------------------------------------------------#

"""
`const PREC = Union{Float16,Float32,Float64,BigFloat}`\n
Concrete PREC type union for parametric abstract types.
"""
const PREC = Union{Float16,Float32,Float64,BigFloat}

"""
`const EXAC = Union{EX,MM}`\n
Concrete EXAC type union for parametric abstract types.
"""
const EXAC = Union{EX,MM}

"""
`const BASE = Union{MA,MO,SY,DT}`\n
Concrete BASE type union for parametric abstract types.
"""
const BASE = Union{MA,MO,SY,DT}

#----------------------------------------------------------------------------------------------#
#                   {PREC[,EXAC[,BASE]]} Parametric Abstract Type Factories                    #
#----------------------------------------------------------------------------------------------#

"""
`function mk1ParAbs(TY::Symbol, TP::Symbol, what::AbstractString, pp::Integer=1,
xp::Bool=true)`\n
Declares a new, 1-parameter abstract type. Parent type parameter count is a function of `pp`, so
that declarations are as follows:\n
- `TY{𝗽} <: TP{𝗽}` for `pp >= 1` (default);
- `TY{𝗽<:PREC} <: TP` for `pp <= 0`.\n
Argument `what` is inserted in the new type documentation, and `xp` controls whether or not the
new abstract type is exported (default `true`).
"""
function mk1ParAbs(TY::Symbol, TP::Symbol, what::AbstractString,
                   pp::Integer=1, xp::Bool=true)
    #if !(eval(TP) isa DataType)
    #    error("Type parent must be a DataType. Got $(string(TP)).")
    #end
    hiStr = tyArchy(eval(TP))
    ppStr = pp>=1 ? "{𝗽}" : ""
    dcStr = """
`abstract type $(TY){𝗽<:PREC} <: $(TP)$(ppStr) end`\n
Abstract supertype for $(what).\n
## Hierarchy\n
`$(TY) <: $(hiStr)`
    """
    if      pp>=1   @eval (abstract type $TY{𝗽} <: $TP{𝗽} end)
    elseif  pp<=0   @eval (abstract type $TY{𝗽<:PREC} <: $TP end)
    end
    @eval begin
        # Type documentation
        @doc $dcStr $TY
        # Type exporting
        if $(xp); export $TY; end
    end
end

"""
`function mk2ParAbs(TY::Symbol, TP::Symbol, what::AbstractString, pp::Integer=2,
xp::Bool=true)`\n
Declares a new, 2-parameter abstract type. Parent type parameter count is a function of `pp`, so
that declarations are as follows:\n
- `TY{𝗽,𝘅} <: TP{𝗽,𝘅}` for `pp >= 2` (default);
- `TY{𝗽,𝘅<:EXAC} <: TP{𝗽}` for `pp = 1`;
- `TY{𝗽<:PREC,𝘅<:EXAC} <: TP` for `pp <= 0`.\n
Argument `what` is inserted in the new type documentation, and `xp` controls whether or not the
new abstract type is exported (default `true`).
"""
function mk2ParAbs(TY::Symbol, TP::Symbol, what::AbstractString,
                   pp::Integer=2, xp::Bool=true)
    #if !(eval(TP) isa DataType)
    #    error("Type parent must be a DataType. Got $(string(TP)).")
    #end
    hiStr = tyArchy(eval(TP))
    ppStr = pp>=2 ? "{𝗽,𝘅}" : pp==1 ? "{𝗽}" : ""
    dcStr = """
`abstract type $(TY){𝗽<:PREC,𝘅<:EXAC} <: $(TP)$(ppStr) end`\n
Abstract supertype for $(what).\n
## Hierarchy\n
`$(TY) <: $(hiStr)`
    """
    if      pp>=2   @eval (abstract type $TY{𝗽,𝘅} <: $TP{𝗽,𝘅} end)
    elseif  pp==1   @eval (abstract type $TY{𝗽,𝘅<:EXAC} <: $TP{𝗽} end)
    elseif  pp<=0   @eval (abstract type $TY{𝗽<:PREC,𝘅<:EXAC} <: $TP end)
    end
    @eval begin
        # Type documentation
        @doc $dcStr $TY
        # Type exporting
        if $(xp); export $TY; end
    end
end

"""
`function mk3ParAbs(TY::Symbol, TP::Symbol, what::AbstractString, pp::Integer=3,
xp::Bool=true)`\n
Declares a new, 3-parameter abstract type. Parent type parameter count is a function of `pp`, so
that declarations are as follows:\n
- `TY{𝗽,𝘅,𝗯} <: TP{𝗽,𝘅,𝗯}` for `pp >= 3` (default);
- `TY{𝗽,𝘅,𝗯<:BASE} <: TP{𝗽,𝘅}` for `pp == 2`;
- `TY{𝗽,𝘅<:EXAC,𝗯<:BASE} <: TP{𝗽}` for `pp = 1`;
- `TY{𝗽<:PREC,𝘅<:EXAC,𝗯<:BASE} <: TP` for `pp <= 0`.\n
Argument `what` is inserted in the new type documentation, and `xp` controls whether or not the
new abstract type is exported (default `true`).
"""
function mk3ParAbs(TY::Symbol, TP::Symbol, what::AbstractString,
                   pp::Integer=3, xp::Bool=true)
    #if !(eval(TP) isa DataType)
    #    error("Type parent must be a DataType. Got $(string(TP)).")
    #end
    hiStr = tyArchy(eval(TP))
    ppStr = pp>=3 ? "{𝗽,𝘅,𝗯}" : pp==2 ? "{𝗽,𝘅}" : pp==1 ? "{𝗽}" : ""
    dcStr = """
`abstract type $(TY){𝗽<:PREC,𝘅<:EXAC,𝗯<:BASE} <: $(TP)$(ppStr) end`\n
Abstract supertype for $(what).\n
## Hierarchy\n
`$(TY) <: $(hiStr)`
    """
    if pp>=3
        @eval (abstract type $TY{𝗽,𝘅,𝗯} <: $TP{𝗽,𝘅,𝗯} end)
    elseif pp==2
        @eval (abstract type $TY{𝗽,𝘅,𝗯<:BASE} <: $TP{𝗽,𝘅} end)
    elseif pp==1
        @eval (abstract type $TY{𝗽,𝘅<:EXAC,𝗯<:BASE} <: $TP{𝗽} end)
    elseif pp<=0
        @eval (abstract type $TY{𝗽<:PREC,𝘅<:EXAC,𝗯<:BASE} <: $TP end)
    end
    @eval begin
        # Type documentation
        @doc $dcStr $TY
        # Type exporting
        if $(xp); export $TY; end
    end
end

#----------------------------------------------------------------------------------------------#
#                      EngTherm Parametric Abstract Supertype Definitions                      #
#----------------------------------------------------------------------------------------------#

# AMOUNT branch — Pars are (i) precision, and (ii) exactness
mk2ParAbs(  :AMOUNTS       , :AbstractTherm, "thermodynamic amounts"                       , 0)
mk2ParAbs(    :WholeAmt    , :AMOUNTS      , "whole, unbased amounts of fixed units"       , 2)
mk2ParAbs(      :WProperty , :WholeAmt     , "whole, unbased properties"                   , 2)
mk2ParAbs(      :WInteract , :WholeAmt     , "whole, unbased interactions"                 , 2)
mk2ParAbs(      :WUnranked , :WholeAmt     , "whole, unbased unranked amounts"             , 2)
mk3ParAbs(    :BasedAmt    , :AMOUNTS      , "based amount groups of fixed units"          , 2)
mk3ParAbs(      :BProperty , :BasedAmt     , "based property groups"                       , 3)
mk3ParAbs(      :BInteract , :BasedAmt     , "based interaction groups"                    , 3)
mk3ParAbs(      :BUnranked , :BasedAmt     , "based unranked amount groups"                , 3)
mk2ParAbs(    :GenerAmt    , :AMOUNTS      , "generic, arbitrary unit amounts"             , 2)

Property{𝗽,𝘅} = Union{WProperty{𝗽,𝘅},BProperty{𝗽,𝘅,𝗯} where 𝗯} where {𝗽,𝘅}
Interact{𝗽,𝘅} = Union{WInteract{𝗽,𝘅},BInteract{𝗽,𝘅,𝗯} where 𝗯} where {𝗽,𝘅}
Unranked{𝗽,𝘅} = Union{WUnranked{𝗽,𝘅},BUnranked{𝗽,𝘅,𝗯} where 𝗯} where {𝗽,𝘅}

export Property, Interact, Unranked

# STATE branch — Pars are (i) precision, and (ii) exactness
mk2ParAbs(  :STATES        , :AbstractTherm, "thermodynamic states"                        , 0)
mk2ParAbs(    :PropPair    , :STATES       , "propery pairs"                               , 2)
mk2ParAbs(    :PropTrio    , :STATES       , "propery trios"                               , 2)
mk2ParAbs(    :PropQuad    , :STATES       , "propery quads"                               , 2)

# MODEL branch — Pars are (i) precision, and (ii) exactness
mk2ParAbs(  :MODELS        , :AbstractTherm, "thermodynamic models"                        , 0)
mk2ParAbs(    :Heat        , :MODELS       , "specific heat models"                        , 2)
mk2ParAbs(      :ConstHeat , :Heat         , "constant specific heat models"               , 2)
mk2ParAbs(      :UnvarHeat , :Heat         , "univariate specific heat models"             , 2)
mk2ParAbs(      :BivarHeat , :Heat         , "bivariate specific heat models"              , 2)
mk2ParAbs(    :Medium      , :MODELS       , "substance/medium models"                     , 2)
mk2ParAbs(      :Substance , :Medium       , "substance model by Equation of State"        , 2)
mk2ParAbs(    :System      , :MODELS       , "system models"                               , 2)
mk2ParAbs(      :Closed    , :System       , "closed systems"                              , 2)
mk2ParAbs(      :Open      , :System       , "open systems"                                , 2)


text:
using CPUTime
using Plots
using DifferentialEquations

f = @ode_def begin
    dx = v
    dv = -4(-D*v + ω^2*l*m*sin(θ))/(-4M + 3m*cos(θ)^2 - 4m) +
    3cos(θ)*(-d + ω + g*l*m*sin(θ))/(l*(-4M + 3m*cos(θ)^2 -4m)) +
    -4/(-4M + 3m*cos(θ)^2 -4m)
    dθ = ω
    dω = 3cos(θ)*(D*v + ω^2*l*m*sin(θ))/(l*(-4M + 3m*cos(θ)^2 -4*m)) +
    -3(M + m)*(-d + ω + g*l*m*sin(θ))/(l^2*(-4M*m + 3m^2*cos(θ)^2 - 4m^2)) +
    3cos(θ)/(l*(-4M + 3m*cos(θ)^2 - 4m))
end M m L l D d g

function do_3_1()
    # パラメータ
    M = 5.0  # 車両質量[kg]
    m = 1.0  # 振子質量[kg]
    L = 1.5  # 振子長さ[m]
    l = L/2
    D = 0.01  # 車両と地面の動摩擦係数
    d = 0.01  # 振子と車両の動摩擦係数
    g = 9.80665  # 重力加速度[m/s²]
    p = [M, m, L, l, D, d, g]

    x0 = [0.0, 0.0, π/6, 0.0]  # 初期値を設定
    tspan = (0.0, 3.0)  # 解く範囲を設定

    prob = ODEProblem(f ,x0 ,tspan, p)
    sol = solve(prob)  # ソルバで解く

    anim = @animate for i in 1:length(sol.t)
        now = round(sol.t[i], digits=1)
        s = "t = " * string(now) * " [s]"
        x, v, θ, ω = sol.u[i]
        plot(
            [x, x+L*cos(-θ+π/2)],
            [0, L*sin(-θ+π/2)],
            #size=(400,400),
            linewidth=5,
            xlims = (-3L, 3L),
            ylims = (-1.1L, 1.1L),
            aspect_ratio = 1,
            legend=false,
            annotate=(0, -1, s),
        )
    end
    gif(anim, "example_3_1_jl.gif", fps = 30)

    return sol, anim

end

@time sol, anim = do_3_1()


text:
using CUDA: CUDA, @cuda
using DataFrames: DataFrame
using LoopVectorization: @avx

export model_var_prob, satisfies_flows_down, satisfies_flows, count_downflow, downflow_all

"Compute the probability of each variable for a random satisfying assignment of the logical circuit"
function model_var_prob(root::LogicCircuit)
    nvars = num_variables(root)
    v, f, _ = satisfies_flows(root, DataFrame(fill(0.5, 1, nvars), :auto))
    f[3:2+nvars]./v[end]
end

#####################
# Bit circuit values and flows (up and downward pass)
#####################

"Compute the value and flow of each node"
function satisfies_flows(circuit::LogicCircuit, data,
                         reuse_values=nothing, reuse_flows=nothing; on_node=noop, on_edge=noop,
                         weights=nothing)
    bc = same_device(BitCircuit(circuit, data), data)
    v, f = satisfies_flows(bc, data, reuse_values, reuse_flows; on_node, on_edge, weights)
    v,f, bc.node2id
end

function satisfies_flows(circuit::BitCircuit, data,
                         reuse_values=nothing, reuse_flows=nothing; on_node=noop, on_edge=noop,
                         weights=nothing)
    @assert isgpu(data) == isgpu(circuit) "BitCircuit and data need to be on the same device"
    values = satisfies_all(circuit, data, reuse_values)
    flows = satisfies_flows_down(circuit, values, reuse_flows; on_node, on_edge, weights)
    return values, flows
end

#####################
# Bit circuit flows downward pass
#####################

"When values of nodes have already been computed, do a downward pass computing the flows at each node"
function satisfies_flows_down(circuit::BitCircuit, values, reuse=nothing; on_node=noop, on_edge=noop, weights=nothing)
    flows = similar!(reuse, typeof(values), size(values)...)
    if !(weights isa Nothing) && values isa CuMatrix
        weights = convert(CuVector{Float32}, weights)
    end
    satisfies_flows_down_layers(circuit, flows, values, on_node, on_edge; weights = weights)
    return flows
end

# downward pass helpers on CPU

"Evaluate the layers of a bit circuit on the CPU (SIMD & multi-threaded)"
function satisfies_flows_down_layers(circuit::BitCircuit, flows::Matrix, values::Matrix, on_node, on_edge; weights = nothing)
    els = circuit.elements
    for layer in Iterators.reverse(circuit.layers)
        Threads.@threads for dec_id in layer
            par_start = @inbounds circuit.nodes[3,dec_id]
            if iszero(par_start)
                if dec_id == num_nodes(circuit)
                    # populate root flow from values
                    @inbounds @views @. flows[:, dec_id] = values[:, dec_id]
                end
                # no parents, ignore (can happen for false/true node and root)
            else
                par_end = @inbounds circuit.nodes[4,dec_id]
                for j = par_start:par_end
                    par = @inbounds circuit.parents[j]
                    grandpa = @inbounds els[1,par]
                    sib_id = sibling(els, par, dec_id)
                    # TODO: check if things would speed up if we also allowed the first parent to write flows in lower-down layers
                    single_child = has_single_child(circuit.nodes, grandpa)
                    if single_child
                        if j == par_start
                            @inbounds @views @. flows[:, dec_id] = flows[:, grandpa]
                        else
                            accum_flow(flows, dec_id, grandpa)
                        end
                    else
                        if j == par_start
                            assign_flow(flows, values, dec_id, grandpa, sib_id)
                        else
                            accum_flow(flows, values, dec_id, grandpa, sib_id)
                        end
                    end
                    # report edge flow only once:
                    sib_id > dec_id && on_edge(flows, values, dec_id, sib_id, par, grandpa, single_child, weights)
                end
            end
            on_node(flows, values, dec_id, weights)
        end
    end
end

assign_flow(f::Matrix{<:Unsigned}, v, d, g, s) =
    @inbounds @views @. f[:, d] = v[:, s] & v[:, d] & f[:, g]

function assign_flow(f::Matrix{<:AbstractFloat}, v, d, g, s)
    @simd for j in 1:size(f,1) # adding @avx here gives incorrect results
        edge_flow = v[j, s] * v[j, d] / v[j, g] * f[j, g]
        edge_flow = ifelse(isfinite(edge_flow), edge_flow, zero(eltype(f)))
        f[j, d] = edge_flow
    end
end

accum_flow(f::Matrix{<:Unsigned}, d, g) =
    @inbounds @views @. f[:, d] |= f[:, g]

accum_flow(f::Matrix{<:AbstractFloat}, d, g) =
    @inbounds @views @. f[:, d] += f[:, g]

accum_flow(f::Matrix{<:Unsigned}, v, d, g, s) =
    @inbounds @views @. f[:, d] |= v[:, s] & v[:, d] & f[:, g]

function accum_flow(f::Matrix{<:AbstractFloat}, v, d, g, s)
    @simd for j in 1:size(f,1) # adding @avx here gives incorrect results
        edge_flow = v[j, s] * v[j, d] / v[j, g] * f[j, g]
        edge_flow = ifelse(isfinite(edge_flow), edge_flow, zero(eltype(f)))
        f[j, d] += edge_flow
    end
end

# # downward pass helpers on GPU

"Pass flows down the layers of a bit circuit on the GPU"
function satisfies_flows_down_layers(circuit::BitCircuit, flows::CuMatrix, values::CuMatrix, on_node, on_edge;
            weights = nothing)
    CUDA.@sync for layer in Iterators.reverse(circuit.layers)
        num_examples = size(values, 1)
        num_decision_sets = length(layer)
        kernel = @cuda name="satisfies_flows_down_layers_cuda" launch=false satisfies_flows_down_layers_cuda(layer, circuit.nodes, circuit.elements, circuit.parents, flows, values, on_node, on_edge, weights)
        config = launch_configuration(kernel.fun)
        threads, blocks =  balance_threads_2d(num_examples, num_decision_sets, config.threads)
        kernel(layer, circuit.nodes, circuit.elements, circuit.parents, flows, values, on_node, on_edge, weights; threads, blocks)
    end
end

"CUDA kernel for passing flows down circuit"
function satisfies_flows_down_layers_cuda(layer, nodes, elements, parents, flows, values, on_node, on_edge, weights::Nothing)
    index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    index_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stride_x = blockDim().x * gridDim().x
    stride_y = blockDim().y * gridDim().y
    for k = index_x:stride_x:size(values,1)
        for i = index_y:stride_y:length(layer)
            dec_id = @inbounds layer[i]
            if dec_id == size(nodes,2)
                # populate root flow from values
                flow = values[k, dec_id]
            else
                par_start = @inbounds nodes[3,dec_id]
                flow = zero(eltype(flows))
                if !iszero(par_start)
                    par_end = @inbounds nodes[4,dec_id]
                    for j = par_start:par_end
                        par = @inbounds parents[j]
                        grandpa = @inbounds elements[1,par]
                        v_gp = @inbounds values[k, grandpa]
                        prime = elements[2,par]
                        sub = elements[3,par]
                        if !iszero(v_gp) # edge flow only gets reported when non-zero
                            f_gp = @inbounds flows[k, grandpa]
                            single_child = has_single_child(nodes, grandpa)
                            if single_child
                                edge_flow = f_gp
                            else
                                v_prime = @inbounds values[k, prime]
                                v_sub = @inbounds values[k, sub]
                                edge_flow = compute_edge_flow( v_prime, v_sub, v_gp, f_gp)
                            end
                            flow = sum_flows(flow, edge_flow)
                            # report edge flow only once:
                            dec_id == prime && on_edge(flows, values, prime, sub, par, grandpa, k, edge_flow, single_child, nothing)
                        end
                    end
                end
            end
            @inbounds flows[k, dec_id] = flow
            on_node(flows, values, dec_id, k, flow, nothing)
        end
    end
    return nothing
end
function satisfies_flows_down_layers_cuda(layer, nodes, elements, parents, flows, values, on_node, on_edge, weights::CUDA.CuDeviceArray{<:AbstractFloat,1,1})
    index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    index_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stride_x = blockDim().x * gridDim().x
    stride_y = blockDim().y * gridDim().y
    for k = index_x:stride_x:size(values,1)
        start_idx = 64 * k - 63
        end_idx = CUDA.min(64 * k, CUDA.length(weights))
        for i = index_y:stride_y:length(layer)
            dec_id = @inbounds layer[i]
            if dec_id == size(nodes,2)
                # populate root flow from values
                flow = values[k, dec_id]
            else
                par_start = @inbounds nodes[3,dec_id]
                flow = zero(eltype(flows))
                if !iszero(par_start)
                    par_end = @inbounds nodes[4,dec_id]
                    for j = par_start:par_end
                        par = @inbounds parents[j]
                        grandpa = @inbounds elements[1,par]
                        v_gp = @inbounds values[k, grandpa]
                        prime = elements[2,par]
                        sub = elements[3,par]
                        if !iszero(v_gp) # edge flow only gets reported when non-zero
                            f_gp = @inbounds flows[k, grandpa]
                            single_child = has_single_child(nodes, grandpa)
                            if single_child
                                edge_flow = f_gp
                            else
                                v_prime = @inbounds values[k, prime]
                                v_sub = @inbounds values[k, sub]
                                edge_flow = compute_edge_flow( v_prime, v_sub, v_gp, f_gp)
                            end
                            flow = sum_flows(flow, edge_flow)
                            # report edge flow only once:
                            if dec_id == prime
                                for bit_idx::UInt32 = 0 : (end_idx - start_idx)
                                    weight::Float32 = @inbounds weights[start_idx + bit_idx]
                                    on_edge(flows, values, prime, sub, par, grandpa, bit_idx, edge_flow, single_child, weight)
                                end
                            end
                        end
                    end
                end
            end
            @inbounds flows[k, dec_id] = flow
            for bit_idx::UInt32 = 0 : (end_idx - start_idx)
                weight::Float32 = @inbounds weights[start_idx + bit_idx]
                on_node(flows, values, dec_id, bit_idx, flow, weight)
            end
        end
    end
    return nothing
end
function satisfies_flows_down_layers_cuda(layer, nodes, elements, parents, flows, values::CUDA.CuDeviceArray{<:AbstractFloat,2,1}, on_node, on_edge, weights::CUDA.CuDeviceArray{<:AbstractFloat,1,1})
    index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    index_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stride_x = blockDim().x * gridDim().x
    stride_y = blockDim().y * gridDim().y
    for k = index_x:stride_x:size(values,1)
        for i = index_y:stride_y:length(layer)
            dec_id = @inbounds layer[i]
            if dec_id == size(nodes,2)
                # populate root flow from values
                flow = values[k, dec_id]
            else
                par_start = @inbounds nodes[3,dec_id]
                flow = zero(eltype(flows))
                if !iszero(par_start)
                    par_end = @inbounds nodes[4,dec_id]
                    for j = par_start:par_end
                        par = @inbounds parents[j]
                        grandpa = @inbounds elements[1,par]
                        v_gp = @inbounds values[k, grandpa]
                        prime = elements[2,par]
                        sub = elements[3,par]
                        if !iszero(v_gp) # edge flow only gets reported when non-zero
                            f_gp = @inbounds flows[k, grandpa]
                            single_child = has_single_child(nodes, grandpa)
                            if single_child
                                edge_flow = f_gp
                            else
                                v_prime = @inbounds values[k, prime]
                                v_sub = @inbounds values[k, sub]
                                edge_flow = compute_edge_flow( v_prime, v_sub, v_gp, f_gp)
                            end
                            flow = sum_flows(flow, edge_flow)
                            # report edge flow only once:
                            dec_id == prime && on_edge(flows, values, prime, sub, par, grandpa, k, edge_flow, single_child, weights[k])
                        end
                    end
                end
            end
            @inbounds flows[k, dec_id] = flow
            on_node(flows, values, dec_id, k, flow, weights[k])
        end
    end
    return nothing
end

@inline sum_flows(a,b::AbstractFloat) = a + b
@inline sum_flows(a,b::Unsigned) = a | b

@inline compute_edge_flow(p_up::AbstractFloat, s_up, n_up, n_down) = p_up * s_up / n_up * n_down
@inline compute_edge_flow(p_up::Unsigned, s_up, n_up, n_down) = p_up & s_up & n_down

# API to get flows

#TODO: these functions need documentation!

function count_downflow(values::Matrix{UInt64}, flows::Matrix{UInt64}, N, n::LogicCircuit, node2id)
    dec_id = node2id[n].node_id
    sum(1:size(flows,1)) do i
        count_ones(flows[i, dec_id])
    end
end

function downflow_all(values::Matrix{UInt64}, flows::Matrix{UInt64}, N, n::LogicCircuit, node2id)
    dec_id = node2id[n].node_id
    indices = map(1:size(flows,1)) do i
        digits(Bool, flows[i, dec_id], base=2, pad=64)
    end
    BitArray(vcat(indices...)[1:N])
end

function count_downflow(values::Matrix{UInt64}, flows::Matrix{UInt64}, N, n::LogicCircuit, c::LogicCircuit, node2id)
    grandpa = node2id[n].node_id
    if isleafgate(c)
        par = node2id[c].node_id
        return sum(1:size(flows,1)) do i
            count_ones(values[i, par] & flows[i, grandpa])
        end
    else
        ids = [node2id[x].node_id for x in children(c)]
        return sum(1:size(flows,1)) do i
            indices = flows[i, grandpa]
            for id in ids
                indices &= values[i, id]
            end
            count_ones(indices)
        end
    end
end

function downflow_all(values::Matrix{UInt64}, flows::Matrix{UInt64}, N, n::LogicCircuit, c::LogicCircuit, node2id)
    grandpa = node2id[n].node_id
    if isleafgate(c)
        par = node2id[c].node_id
        edge = map(1:size(flows,1)) do i
            digits(Bool, values[i, par] & flows[i, grandpa], base=2, pad=64)
        end
    else
        ids = [node2id[x].node_id for x in children(c)]
        edge = map(1:size(flows,1)) do i
            indices = flows[i, grandpa]
            for id in ids
                indices = indices & values[i, id]
            end
            digits(Bool, indices, base=2, pad=64)
        end
    end
    BitArray(vcat(edge...)[1:N])
end

function count_downflow(values::Matrix{<:AbstractFloat}, flows::Matrix{<:AbstractFloat}, N, n::LogicCircuit, node2id)
    return sum(downflow_all(values, flows, N, n, node2id))
end

function downflow_all(values::Matrix{<:AbstractFloat}, flows::Matrix{<:AbstractFloat}, N, n::LogicCircuit, node2id)
    dec_id = node2id[n].node_id
    map(1:size(flows, 1)) do i
        flows[i, dec_id]
    end
end

function count_downflow(values::Matrix{<:AbstractFloat}, flows::Matrix{<:AbstractFloat}, N, n::LogicCircuit, c::LogicCircuit, node2id)
    sum(downflow_all(values, flows, N, n, c, node2id))
end

function downflow_all(values::Matrix{<:AbstractFloat}, flows::Matrix{<:AbstractFloat}, N, n::LogicCircuit, c::LogicCircuit, node2id)
    grandpa = node2id[n].node_id
    if isleafgate(c)
        par = node2id[c].node_id
        edge_flows = map(1:size(flows,1)) do i
            values[i, par] * flows[i, grandpa] / values[i, grandpa]
        end
        return edge_flows
    else
        ids = [node2id[x].node_id for x in children(c)]
        edge_flows = map(1:size(flows,1)) do i
            n_down = flows[i, grandpa]
            n_up = values[i, grandpa]
            edge_flow = 1.0
            for id in ids
                edge_flow = edge_flow * values[i, id]
            end
            edge_flow * n_down / n_up
        end
        return edge_flows
    end
end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
# -------------------------------------------
function forward!(net::Network, X::Signal)
    for i in 1:length(net.layers)
        reset!(net.layers[i], X)
        X = forward!(net.layers[i], X)
    end
    X
end

function forward!(l::Layer, X::Signal)
    copy!(l.X, X)
    @into! l.U = l.W * X
    u = @in1! l.U .+ l.b

    forward_activation!(l)

    l.Z
end

function forward_activation!(l::Layer)
    forward!(l.activation, l.U, l.Z)
    forward!(l.post_filters, l)
end

function forward_activation!(l::SimpleRecurrentLayer)
    forward!(l.activation, sub(l.U,:,1), sub(l.Z,:,1))
    forward!(l.post_filters, l; time_step=1)

    @inbounds for t in 2:size(l.U, 2)
        l.Uh[:,t] = l.Wh * l.Z[:,t-1]
        @inbounds for i in 1:size(l.Uh,1)
            l.U[i,t] += l.Uh[i,t]
        end
        forward!(l.activation, sub(l.U,:,t), sub(l.Z,:,t))
        forward!(l.post_filters, l; time_step=t)
    end

    if maximum(l.Z) > 100000.0f0
        @show maximum(l.Z), maximum(l.U), maximum(l.Uh)
        error("Large recurrent forward")
    end
end

# -----------------------------------------------------------
function backprop!(net::Network, X::Signal, Y::Signal)
    forward!(net, X)
    for i in length(net.layers):-1:1
        backprop!(net.layers[i], Y)
        Y = net.layers[i]
    end
end

function backprop!(l::OutputLayer, Y::Signal)
    copy!(l.ΔE, l.loss(:diff, l, Y, l.Z))
    backprop!(l.post_filters, l)
    backprop!(l)
end

function backprop!(l::HiddenLayer, l2::Layer)
    Δout = l2.W' * l2.ΔE
    backward!(l.activation, l.U, l.Zd)
    copy!(l.ΔE, (l.Zd .* Δout))
    backprop!(l.post_filters, l)
    backprop!(l)
end

function backprop!(l::SimpleRecurrentLayer, l2::Layer)
    Δout = l2.W' * l2.ΔE
    backward!(l.activation, l.U, l.Zd)

    # Recurrent
    l.ΔE[:,end:end] = l.Zd[:,end:end] .* Δout[:,end:end]
    backprop!(l.post_filters, l; time_step=size(l.ΔE,2))

    Δfuture = zeros(Float32, size(l.Wh,1))

    @inbounds for t in (size(l.ΔE,2)-1):-1:1
        @into! Δfuture = l.Wh' * l.ΔE[:,t+1]
        @inbounds for i in 1:size(l.Zd,1)
            l.ΔE[i,t] = l.Zd[i,t] .* (Δout[i,t] + Δfuture[i])
        end
        backprop!(l.post_filters, l; time_step=t)
    end

    @into! l.ΔWh = l.ΔE * l.Z'
    _ = @in1! l.ΔWh ./ l.batch_size
    backprop!(l)
end

function backprop!(l::Layer)
    copy!(l.ΔW, ((l.ΔE * l.X')               ./ l.batch_size))
    copy!(l.Δb, ((l.ΔE * ones(l.batch_size)) ./ l.batch_size))
end


text:
# Affinity Propagation clustering
# https://juliastats.org/Clustering.jl/stable/affprop.html#Affinity-Propagation-1

function AffProp(Sim)
    # If it fails to cluster, then they all stay in the same cluster
    Ind = ones(size(Sim,1))
    try
        AFFPROP = affinityprop(100*Sim) # Affinity propagation clustering
        Ind = AFFPROP.assignments
    catch
    end

    return Ind
end

max_stars_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_issues_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_forks_repo_name: gonzalez-alvaro/MultidimensionalGAPM


text:
module ModuleTypes
using StaticArrays
using GeometryTypes

export HomogeneousPoint, ProjectiveEntity, FundamentalMatrix, ProjectionMatrix
export HomographyMatrix, HomographyMatrices
export TotalViews, TwoViews
export ProjectionMatrices, FactorisationAlgorithm, Chojnacki, LatentVariables
export HomogeneousCoordinates, EssentialMatrix
export CameraModel, Pinhole, CanonicalLens
export EstimationAlgorithm, DirectLinearTransform, Taubin, FundamentalNumericalScheme
export BundleAdjustment
export CostFunction, ApproximateMaximumLikelihood, AML
export CoordinateSystemTransformation, CanonicalToHartley, HartleyToCanonical
export CovarianceMatrices
export Point2D, Point2DH, Point3D, Point3DH
export HessianApproximation, CanonicalApproximation, CovarianceEstimationScheme
export NoiseModel, GaussianNoise
export RasterSystem, CartesianSystem, OpticalSystem
include("types.jl")
end


text:
using ITensors.NDTensors,
      Test

@testset "BlockSparseTensor basic functionality" begin

  # Indices
  indsA = ([2,3],[4,5])

  # Locations of non-zero blocks
  locs = [(1,2),(2,1)]

  A = BlockSparseTensor(locs,indsA...)
  randn!(A)

  @test blockdims(A,(1,2)) == (2,5)
  @test blockdims(A,(2,1)) == (3,4)
  @test nnzblocks(A) == 2
  @test nnz(A) == 2*5+3*4
  @test inds(A) == ([2,3],[4,5])
  @test isblocknz(A,(2,1))
  @test isblocknz(A,(1,2))
  @test !isblocknz(A,(1,1))
  @test !isblocknz(A,(2,2))
  @test findblock(A,(2,1))==1
  @test findblock(A,(1,2))==2
  @test isnothing(findblock(A,(1,1)))
  @test isnothing(findblock(A,(2,2)))

  # Test different ways of getting nnz
  @test nnz(blockoffsets(A),inds(A)) == nnz(A)

  A[1,5] = 15
  A[2,5] = 25

  @test A[1,1] == 0
  @test A[1,5] == 15
  @test A[2,5] == 25

  D = dense(A)

  @test D == A

  for I in eachindex(A)
    @test D[I] == A[I]
  end

  A12 = blockview(A,(1,2))

  @test dims(A12) == (2,5)

  for I in eachindex(A12)
    @test A12[I] == A[I+CartesianIndex(0,4)]
  end

  B = BlockSparseTensor(undef,locs,indsA)
  randn!(B)

  C = A+B

  for I in eachindex(C)
    @test C[I] == A[I]+B[I]
  end

  Ap = permutedims(A,(2,1))

  @test blockdims(Ap,(1,2)) == (4,3)
  @test blockdims(Ap,(2,1)) == (5,2)
  @test nnz(A) == nnz(Ap)
  @test nnzblocks(A) == nnzblocks(Ap)

  for I in eachindex(C)
    @test A[I] == Ap[NDTensors.permute(I,(2,1))]
  end

  A = BlockSparseTensor(ComplexF64,locs,indsA)
  randn!(A)
  @test conj(data(store(A))) == data(store(conj(A)))
  @test typeof(conj(A)) <: BlockSparseTensor

  @testset "BlockSparseTensor setindex! add block" begin
    T = BlockSparseTensor([2,3],[4,5])

    for I in eachindex(C)
      @test T[I] == 0.0
    end
    @test nnz(T) == 0
    @test nnzblocks(T) == 0
    @test !isblocknz(T,(1,1))
    @test !isblocknz(T,(2,1))
    @test !isblocknz(T,(1,2))
    @test !isblocknz(T,(2,2))

    T[1,1] = 1.0

    @test T[1,1] == 1.0
    @test nnz(T) == 8
    @test nnzblocks(T) == 1
    @test isblocknz(T,(1,1))
    @test !isblocknz(T,(2,1))
    @test !isblocknz(T,(1,2))
    @test !isblocknz(T,(2,2))

    T[4,8] = 2.0

    @test T[4,8] == 2.0
    @test nnz(T) == 8+15
    @test nnzblocks(T) == 2
    @test isblocknz(T,(1,1))
    @test !isblocknz(T,(2,1))
    @test !isblocknz(T,(1,2))
    @test isblocknz(T,(2,2))

    T[1,6] = 3.0

    @test T[1,6] == 3.0
    @test nnz(T) == 8+15+10
    @test nnzblocks(T) == 3
    @test isblocknz(T,(1,1))
    @test !isblocknz(T,(2,1))
    @test isblocknz(T,(1,2))
    @test isblocknz(T,(2,2))

    T[4,2] = 4.0

    @test T[4,2] == 4.0
    @test nnz(T) == 8+15+10+12
    @test nnzblocks(T) == 4
    @test isblocknz(T,(1,1))
    @test isblocknz(T,(2,1))
    @test isblocknz(T,(1,2))
    @test isblocknz(T,(2,2))
  end

  @testset "Add with different blocks" begin
    # Indices
    inds = ([2,3],[4,5])

    # Locations of non-zero blocks
    locsA = [(1,1),(1,2),(2,2)]
    A = BlockSparseTensor(locsA,inds...)
    randn!(A)

    locsB = [(1,2),(2,1)]
    B = BlockSparseTensor(locsB,inds...)
    randn!(B)

    R = A+B

    @test nnz(R) == dim(R)
    for I in eachindex(R)
      @test R[I] == A[I] + B[I]
    end
  end

  @testset "permutedims!! with different blocks" begin
    # Indices
    indsA = ([2,3],[4,5])

    # Locations of non-zero blocks
    locsA = [(1,2),(2,1)]
    A = BlockSparseTensor(locsA,indsA...)
    randn!(A)

    perm = (2,1)

    locsB = [(2,1)]
    indsB = NDTensors.permute(indsA, perm)
    B = BlockSparseTensor(locsB,indsB...)
    randn!(B)

    R = permutedims!!(B,A,perm)

    @test nnz(R) == nnz(A)
    for I in eachindex(A)
      @test R[NDTensors.permute(I,perm)] == A[I]
    end
  end

  @testset "Contract" begin
    indsA = ([2,3],[4,5])
    locsA = [(1,1),(2,2),(2,1),(1,2)]
    A = BlockSparseTensor(locsA,indsA...)
    randn!(A)

    indsB = ([4,5],[3,2])
    locsB = [(1,2),(2,1),(1,1)]
    B = BlockSparseTensor(locsB,indsB...)
    randn!(B)

    R = contract(A,(1,-1),B,(-1,2))

    DA = dense(A)
    DB = dense(B)
    DR = contract(DA,(1,-1),DB,(-1,2))

    for I in eachindex(R)
      @test R[I] ≈ DR[I]
    end
  end

  @testset "reshape" begin
    indsA = ([2,3],[4,5])
    locsA = [(2,1),(1,2)]
    A = BlockSparseTensor(locsA,indsA...)
    randn!(A)

    indsB = ([8,12,10,15],)
    B = reshape(A,indsB)

    @test nnzblocks(A)==nnzblocks(B)
    @test nnz(A)==nnz(B)
    for i in 1:nnzblocks(B)
      blockA = blockview(A,i)
      blockB = blockview(B,i)
      @test reshape(blockA,size(blockB))==blockB
    end
  end

  @testset "permute_combine" begin
    indsA = ([2,3],[4,5],[6,7,8])
    locsA = [(2,1,1),(1,2,1),(2,2,3)]
    A = BlockSparseTensor(locsA,indsA...)
    randn!(A)

    B = NDTensors.permute_combine(A,3,(2,1))

    @test nnzblocks(A)==nnzblocks(B)
    @test nnz(A)==nnz(B)

    Ap = permutedims(A,(3,2,1))

    for i in 1:nnzblocks(A)
      blockAp = blockview(Ap,i)
      blockB = blockview(B,i)
      @test reshape(blockAp,size(blockB))==blockB
    end
  end

  @testset "svd" begin

    @testset "svd example 1" begin
      A = BlockSparseTensor([(2,1),(1,2)],[2,2],[2,2])
      randn!(A)
      U,S,V = svd(A)
      @test isapprox(norm(array(U)*array(S)*array(V)'-array(A)),0; atol=1e-14)
    end

    @testset "svd example 2" begin
      A = BlockSparseTensor([(1,2),(2,3)],[2,2],[3,2,3])
      randn!(A)
      U,S,V = svd(A)
      @test isapprox(norm(array(U)*array(S)*array(V)'-array(A)),0.0; atol=1e-14)
    end

    @testset "svd example 3" begin
      A = BlockSparseTensor([(2,1),(3,2)],[3,2,3],[2,2])
      randn!(A)
      U,S,V = svd(A)
      @test isapprox(norm(array(U)*array(S)*array(V)'-array(A)),0.0; atol=1e-14)
    end

    @testset "svd example 4" begin
      A = BlockSparseTensor([(2,1),(3,2)],[2,3,4],[5,6])
      randn!(A)
      U,S,V = svd(A)
      @test isapprox(norm(array(U)*array(S)*array(V)'-array(A)),0.0; atol=1e-13)
    end

    @testset "svd example 5" begin
      A = BlockSparseTensor([(1,2),(2,3)],[5,6],[2,3,4])
      randn!(A)
      U,S,V = svd(A)
      @test isapprox(norm(array(U)*array(S)*array(V)'-array(A)),0.0; atol=1e-14)
    end
  end

end


text:
@testset "Core Test" begin
  agent = Agent()

  # Layer
  l1 = addnode!(agent, GenericLayer(8))
  @test length(agent.nodes) == 1
  l2 = addnode!(agent, GenericLayer(8))
  addedge!(agent, l2, (
    (:input, l1),
  ))
  @test length(agent.edges) == 1
  run!(agent)

  lt = addnode!(agent, GenericLayer(8))
  addedge!(agent, lt, (
    (:input, l1),
  ))
  deledge!(agent, lt, l1, :input)
  @test length(agent.edges) == 1
  delnode!(agent, lt)
  @test length(agent.nodes) == 2

  # Deledges!
  lt = addnode!(agent, GenericLayer(8))
  addedge!(agent, lt, (
    (:input, l1),
  ))
  deledges!(agent, lt)
  @test length(root(agent.nodes[lt]).synapses.conns) == 0

  # PatchClamp
  p1 = addnode!(agent, PatchClamp(Bool, 8))
  p2 = addnode!(agent, PatchClamp(8))
  run!(agent)

  # Layer and PatchClamp edges
  addedge!(agent, l1, (
    (:input, p1),
  ))
  addedge!(agent, p2, (
    (:input, l2),
  ))
  run!(agent)

  # RewardModulator
  p3 = addnode!(agent, PatchClamp(Bool, 8))
  addedge!(agent, l2, (
    (:input, p3, (modulator=RewardModulator(),)),
  ))
  run!(agent)

  # Merge
  agent2 = Agent()
  l3 = addnode!(agent2, GenericLayer(16))
  merge!(agent, agent2)
  addedge!(agent, l3, (
    (:input, l2),
  ))
  run!(agent)

  # Clear
  reinit!(agent)

  # Deepcopy
  agent3 = deepcopy(agent)
  run!(agent3)

  # Size
  for (idx, node) in agent.nodes
    @test size(root(node)) != nothing
  end
end


text:
module test_02_OneScalarSignalWithUnit

using ModiaPlot

include("test_base/test_02_OneScalarSignalWithUnit.jl")

end


text:
@enum ProjectionEnum Perspective Orthographic

struct Camera3D <: AbstractCamera
    rotationspeed::Node{Float32}
    translationspeed::Node{Float32}
    eyeposition::Node{Vec3f0}
    lookat::Node{Vec3f0}
    upvector::Node{Vec3f0}
    fov::Node{Float32}
    near::Node{Float32}
    far::Node{Float32}
    projectiontype::Node{ProjectionEnum}
    pan_button::Node{ButtonTypes}
    rotate_button::Node{ButtonTypes}
    move_key::Node{ButtonTypes}
end

"""
    cam3d_cad!(scene; kw_args...)

Creates a 3D camera for `scene` which rotates around
the _viewer_'s "up" axis - similarly to how it's done
in CAD software cameras.
"""
function cam3d_cad!(scene; kw_args...)
    cam_attributes = merged_get!(:cam3d, scene, Attributes(kw_args)) do
        Attributes(
            rotationspeed = 0.01,
            translationspeed = 1.0,
            eyeposition = Vec3f0(3),
            lookat = Vec3f0(0),
            upvector = Vec3f0(0, 0, 1),
            fov = 45f0,
            near = 0.01f0,
            far = 100f0,
            projectiontype = Perspective,
            pan_button = Mouse.right,
            rotate_button = Mouse.left,
            move_key = nothing
        )
    end
    cam = from_dict(Camera3D, cam_attributes)
    # remove previously connected camera
    disconnect!(scene.camera)
    add_translation!(scene, cam, cam.pan_button, cam.move_key, false)
    add_rotation!(scene, cam, cam.rotate_button, cam.move_key, false)
    cameracontrols!(scene, cam)
    on(camera(scene), scene.px_area) do area
        # update cam when screen ratio changes
        update_cam!(scene, cam)
    end
    cam
end

"""
    cam3d_turntable!(scene; kw_args...)

Creates a 3D camera for `scene`, which rotates around
the plot's axis.
"""
function cam3d_turntable!(scene; kw_args...)
    cam_attributes = merged_get!(:cam3d, scene, Attributes(kw_args)) do
        Attributes(
            rotationspeed = 0.01,
            translationspeed = 1.0,
            eyeposition = Vec3f0(3),
            lookat = Vec3f0(0),
            upvector = Vec3f0(0, 0, 1),
            fov = 45f0,
            near = 0.01f0,
            far = 100f0,
            projectiontype = Perspective,
            pan_button = Mouse.right,
            rotate_button = Mouse.left,
            move_key = nothing
        )
    end
    cam = from_dict(Camera3D, cam_attributes)
    # remove previously connected camera
    disconnect!(scene.camera)
    add_translation!(scene, cam, cam.pan_button, cam.move_key, true)
    add_rotation!(scene, cam, cam.rotate_button, cam.move_key, true)
    cameracontrols!(scene, cam)
    on(camera(scene), scene.px_area) do area
        # update cam when screen ratio changes
        update_cam!(scene, cam)
    end
    cam
end

"""
    cam3d!(scene; kwargs...)

An alias to [`cam3d_turntable!`](@ref).
Creates a 3D camera for `scene`, which rotates around
the plot's axis.
"""
const cam3d! = cam3d_turntable!

function projection_switch(
        wh::Rect2D,
        fov::T, near::T, far::T,
        projectiontype::ProjectionEnum, zoom::T
    ) where T <: Real
    aspect = T((/)(widths(wh)...))
    h = T(tan(fov / 360.0 * pi) * near)
    w = T(h * aspect)
    projectiontype == Perspective && return frustum(-w, w, -h, h, near, far)
    h, w = h * zoom, w * zoom
    orthographicprojection(-w, w, -h, h, near, far)
end

function rotate_cam(
        theta::Vec{3, T},
        cam_right::Vec{3, T}, cam_up::Vec{3, T}, cam_dir::Vec{3, T}
    ) where T
    rotation = Quaternion{T}(0, 0, 0, 1)
    if !all(isfinite.(theta))
        # We can only rotate for finite values
        # Makie#338
        return rotation
    end
    # first the rotation around up axis, since the other rotation should be relative to that rotation
    if theta[1] != 0
        rotation *= qrotation(cam_up, theta[1])
    end
    # then right rotation
    if theta[2] != 0
        rotation *= qrotation(cam_right, theta[2])
    end
    # last rotation around camera axis
    if theta[3] != 0
        rotation *= qrotation(cam_dir, theta[3])
    end
    rotation
end

# TODO switch button and key because this is the wrong order
function add_translation!(scene, cam, key, button, zoom_shift_lookat::Bool)
    last_mousepos = RefValue(Vec2f0(0, 0))
    dragging = RefValue(false)
    on(camera(scene), scene.events.mousebutton) do event
        if event.button == key[] && ispressed(scene, button[])
            if event.action == Mouse.press && is_mouseinside(scene)
                last_mousepos[] = mouseposition_px(scene)
                dragging[] = true
                return true
            elseif event.action == Mouse.release && dragging[]
                mousepos = mouseposition_px(scene)
                dragging[] = false
                diff = (last_mousepos[] - mousepos) * cam.translationspeed[]
                last_mousepos[] = mousepos
                translate_cam!(scene, cam, Vec3f0(0f0, diff[1], diff[2]))
                return true
            end
        end
        return false
    end

    on(camera(scene), scene.events.mouseposition) do mp
        if dragging[] && ispressed(scene, button[])
            mousepos = screen_relative(scene, mp)
            diff = (last_mousepos[] .- mousepos) * cam.translationspeed[]
            last_mousepos[] = mousepos
            translate_cam!(scene, cam, Vec3f0(0f0, diff[1], diff[2]))
            return true
        end
        return false
    end

    on(camera(scene), scene.events.scroll) do scroll
        if ispressed(scene, button[]) && is_mouseinside(scene)
            cam_res = Vec2f0(widths(scene.px_area[]))
            mouse_pos_normalized = mouseposition_px(scene) ./ cam_res
            mouse_pos_normalized = 2*mouse_pos_normalized .- 1f0
            zoom_step = scroll[2]
            zoom!(scene, mouse_pos_normalized, zoom_step, zoom_shift_lookat)
            return true
        end
        return false
    end
end

function add_rotation!(scene, cam, button, key, fixed_axis::Bool)
    last_mousepos = RefValue(Vec2f0(0, 0))
    dragging = RefValue(false)
    e = events(scene)

    on(camera(scene), e.mousebutton) do event
        if event.button == button[] && ispressed(scene, key[])
            if event.action == Mouse.press && is_mouseinside(scene)
                last_mousepos[] = mouseposition_px(scene)
                dragging[] = true
                return true
            elseif event.action == Mouse.release && dragging[]
                mousepos = mouseposition_px(scene)
                dragging[] = false
                rot_scaling = cam.rotationspeed[] * (e.window_dpi[] * 0.005)
                mp = (last_mousepos[] - mousepos) * rot_scaling
                last_mousepos[] = mousepos
                rotate_cam!(scene, cam, Vec3f0(mp[1], -mp[2], 0f0), fixed_axis)
                return true
            end
        end
        return false
    end

    on(camera(scene), e.mouseposition) do mp
        if dragging[] && ispressed(scene, key[])
            mousepos = screen_relative(scene, mp)
            rot_scaling = cam.rotationspeed[] * (e.window_dpi[] * 0.005)
            mp = (last_mousepos[] .- mousepos) * rot_scaling
            last_mousepos[] = mousepos
            rotate_cam!(scene, cam, Vec3f0(mp[1], -mp[2], 0f0), fixed_axis)
            return true
        end
        return false
    end
end

"""
    translate_cam!(scene::Scene. translation::VecTypes)

Translate the camera to the given coordinates.
"""
translate_cam!(scene::Scene, translation::VecTypes) = translate_cam!(scene, cameracontrols(scene), translation)
function translate_cam!(scene::Scene, cam::Camera3D, _translation::VecTypes)
    translation = Vec3f0(_translation)
    translation == Vec3f0(0) && return
    @extractvalue cam (projectiontype, lookat, eyeposition, upvector)

    dir = eyeposition - lookat
    dir_len = norm(dir)
    cam_res = Vec2f0(widths(scene.px_area[]))
    z, x, y = translation
    z *= 0.1f0 * dir_len

    x, y = (Vec2f0(x, y) ./ cam_res) .* dir_len

    dir_norm = normalize(dir)
    right = normalize(cross(dir_norm, upvector))
    z_trans = dir_norm * z
    side_trans = right * (-x) + normalize(upvector) * y
    newpos = eyeposition + side_trans + z_trans
    cam.eyeposition[] = newpos
    cam.lookat[] = lookat + side_trans
    update_cam!(scene, cam)
    return
end

"""
    zoom!(scene, point, zoom_step, shift_lookat::Bool)

Zooms the camera of `scene` in towards `point` by a factor of `zoom_step`. A positive
`zoom_step` zooms in while a negative `zoom_step` zooms out.
"""
function zoom!(scene, point, zoom_step, shift_lookat::Bool)
    cam = cameracontrols(scene)
    @extractvalue cam (projectiontype, lookat, eyeposition, upvector, projectiontype)

    # split zoom into two components:
    # the offset perpendicular to `eyeposition - lookat`, based on mouse offset ~ ray_dir
    # the offset parallel to `eyeposition - lookat` ~ dir
    ray_eye = inv(scene.camera.projection[]) * Vec4f0(point[1],point[2],0,0)
    ray_eye = Vec4f0(ray_eye[1:2]...,0,0)
    ray_dir = Vec3f0((inv(scene.camera.view[]) * ray_eye))

    dir = eyeposition - lookat

    if shift_lookat
        # This results in the point under the cursor remaining stationary
        if projectiontype == Perspective
            ray_dir *= norm(dir)
        end
        cam.eyeposition[] = eyeposition + (ray_dir - dir) * (1f0 - 0.9f0 ^ zoom_step)
        cam.lookat[] = lookat + (1f0 - 0.9f0 ^ zoom_step) * ray_dir
    else
        # Rotations need more extreme eyeposition shifts
        step = zoom_step
        while abs(step) > 1f0
            cam.eyeposition[] = cam.eyeposition[] + sign(zoom_step) * (ray_dir - dir * 0.1f0)
            dir = cam.eyeposition[] - lookat
            step -= sign(step)
        end
    end

    update_cam!(scene, cam)
end

"""
    rotate_cam!(scene::Scene, theta_v::Number...)
    rotate_cam!(scene::Scene, theta_v::VecTypes)

Rotate the camera of the Scene by the given rotation. Passing `theta_v = (α, β, γ)` will rotate
the camera according to the Euler angles (α, β, γ).
"""
rotate_cam!(scene::Scene, theta_v::Number...) = rotate_cam!(scene, cameracontrols(scene), theta_v)
rotate_cam!(scene::Scene, theta_v::VecTypes) = rotate_cam!(scene, cameracontrols(scene), theta_v)
function rotate_cam!(scene::Scene, cam::Camera3D, _theta_v::VecTypes, fixed_axis::Bool = true)
    theta_v = Vec3f0(_theta_v)
    theta_v == Vec3f0(0) && return #nothing to do!
    @extractvalue cam (eyeposition, lookat, upvector)

    dir = normalize(eyeposition - lookat)
    right_v = normalize(cross(upvector, dir))
    upvector = normalize(cross(dir, right_v))
    axis = fixed_axis ? Vec3f0(0, 0, sign(upvector[3])) : upvector
    rotation = rotate_cam(theta_v, right_v, axis, dir)
    r_eyepos = lookat + rotation * (eyeposition - lookat)
    r_up = normalize(rotation * upvector)
    cam.eyeposition[] = r_eyepos
    cam.upvector[] = r_up
    update_cam!(scene, cam)
    return
end

function update_cam!(scene::Scene, cam::Camera3D)
    @extractvalue cam (fov, near, projectiontype, lookat, eyeposition, upvector)

    zoom = norm(lookat - eyeposition)
    # TODO this means you can't set FarClip... SAD!
    # TODO use boundingbox(scene) for optimal far/near
    far = max(zoom * 5f0, 30f0)
    proj = projection_switch(scene.px_area[], fov, near, far, projectiontype, zoom)
    view = Makie.lookat(eyeposition, lookat, upvector)

    scene.camera.projection[] = proj
    scene.camera.view[] = view
    scene.camera.projectionview[] = proj * view
    scene.camera.eyeposition[] = cam.eyeposition[]
end

function update_cam!(scene::Scene, camera::Camera3D, area3d::Rect)
    @extractvalue camera (fov, near, lookat, eyeposition, upvector)
    bb = FRect3D(area3d)
    width = widths(bb)
    half_width = width/2f0
    lower_corner = minimum(bb)
    middle = maximum(bb) - half_width
    old_dir = normalize(eyeposition .- lookat)
    camera.lookat[] = middle
    neweyepos = middle .+ (1.2*norm(width) .* old_dir)
    camera.eyeposition[] = neweyepos
    camera.upvector[] = Vec3f0(0,0,1)
    camera.near[] = 0.1f0 * norm(widths(bb))
    camera.far[] = 3f0 * norm(widths(bb))
    update_cam!(scene, camera)
    return
end

"""
    update_cam!(scene::Scene, eyeposition, lookat, up = Vec3f0(0, 0, 1))

Updates the camera's controls to point to the specified location.
"""
update_cam!(scene::Scene, eyeposition, lookat, up = Vec3f0(0, 0, 1)) = update_cam!(scene, cameracontrols(scene), eyeposition, lookat, up)

function update_cam!(scene::Scene, camera::Camera3D, eyeposition, lookat, up = Vec3f0(0, 0, 1))
    camera.lookat[] = Vec3f0(lookat)
    camera.eyeposition[] = Vec3f0(eyeposition)
    camera.upvector[] = Vec3f0(up)
    update_cam!(scene, camera)
    return
end


text:
# An N-dimensional light-dark problem, originally used to test MCVI
# A very simple POMDP with continuous state and observation spaces.
# maintained by @zsunberg

import Base: ==, +, *, -

"""
    LightDarkNDState

## Fields
- `loc`: position
- `status`: 0 = normal, negative = terminal
"""
struct LightDarkNDState{N}
    status::Int64
    loc::SVector{N, Float64}
end

LightDarkNDState(status, loc) = LightDarkNDState{length(loc)}(status, loc)

*(n::Number, s::LightDarkNDState) = LightDarkNDState(s.status, n*s.loc)

# The light source is at (5,5,...,5)
default_sigma(x) = fill(norm(x .- 5.)/sqrt(2 * length(x)) + 1e-2, length(x))

# Each dimension is independently observed
independent_sigma(x) = abs.(x .- 5.)./sqrt(2) .+ 1e-2

"""
    LightDarkND

An N-dimensional light dark problem. The goal is to be near [0,..,0]∈R^N. Observations are noisy measurements of the position.

Model
-----

   -3-2-1 0 1 2 3
...| | | | | | | | ...
          G   S

Here G is the goal. S is the starting location
"""
@with_kw struct LightDarkND{N,F<:Function} <: POMDPs.POMDP{LightDarkNDState{N},Int,Vector{Float64}}
    discount_factor::Float64 = 0.95
    correct_r::Float64 = 10.0
    incorrect_r::Float64 = -10.0
    step_size::Matrix{Float64} = Matrix(1.0I, N, N)
    movement_cost::Float64 = 0.0
    sigma::F = default_sigma
end

discount(p::LightDarkND) = p.discount_factor

isterminal(::LightDarkND, act::Int64) = act == 0

isterminal(::LightDarkND, s::LightDarkNDState) = s.status < 0

actions(::LightDarkND{N}) where N = -N:N

struct LDNormalNDStateDist{N}
    mean::SVector{N,Float64}
    std::SVector{N,Float64}
end

sampletype(::Type{LDNormalNDStateDist{N}}) where N = LightDarkNDState{N}
rand(rng::AbstractRNG, d::LDNormalNDStateDist{N}) where N = LightDarkNDState(0, d.mean + randn(rng, N).*d.std)
initialstate(pomdp::LightDarkND{N}) where N = LDNormalNDStateDist{N}(fill(2.0, N), fill(3.0, N))
initialobs(m::LightDarkND, s) = observation(m, s)

observation(p::LightDarkND, sp::LightDarkNDState) = MvNormal(sp.loc, p.sigma(sp.loc))

function transition(p::LightDarkND, s::LightDarkNDState{N}, a::Int) where N
    if a == 0
        return Deterministic(LightDarkNDState{N}(-1, s.loc))
    else
        sgn = sign(a)
        return Deterministic(LightDarkNDState{N}(s.status, s.loc .+ sgn .* p.step_size[sgn * a, :]))
    end
end

function reward(p::LightDarkND, s::LightDarkNDState, a::Int)
    if s.status < 0
        return 0.0
    elseif a == 0
        if norm(s.loc) < 1.0
            return p.correct_r
        else
            return p.incorrect_r
        end
    else
        return -p.movement_cost
    end
end

convert_s(::Type{A}, s::LightDarkNDState{N}, p::LightDarkND) where {N, A<:AbstractArray} = eltype(A)[s.status, s.loc...]
convert_s(::Type{LightDarkNDState{N}}, s::AbstractArray, p::LightDarkND) where N = LightDarkNDState{N}(Int64(s[1]), s[2:end])


text:
_t₁(x, df::Float64) = tdistpdf(df, x)
_t₁(x::Vector, df::Float64) = map(u->tdistpdf(df, u), x)
_logt₁(x, df::Float64) = tdistlogpdf(df, x)
_logt₁(x::Vector, df::Float64) = map(u->tdistlogpdf(df, u), x)
_T₁(x, df::Float64) = tdistcdf(df, x)
_T₁(x::Vector, df::Float64) = map(u->tdistcdf(df, u), x)
_logT₁(x, df::Float64) = tdistlogcdf(df, x)
_logT₁(x::Vector, df::Float64) = map(u->tdistlogcdf(df, u), x)

function _log_g(Q, ν::Float64, k::Int)
    # Adapted from mvtdist_consts in Distributions package
    hdf = 0.5 * ν
    hdim = 0.5 * k
    shdfhdim = hdf + hdim
    v = lgamma(shdfhdim) - lgamma(hdf) - hdim*log(ν) - hdim*log(pi)
    return v - (shdfhdim * log(1 + Q/ν))
end

function scale_and_cor(Ω::Matrix{Float64})
    Ωz = similar(Ω)
    d = size(Ω,1)
    ω = sqrt(diag(Ω))
    for i in 1:d, j in 1:d
        Ωz[i,j] = Ω[i,j]/(ω[i]*ω[j])
    end
    diagm(ω), Ωz
end

function scale_and_cor(Ω::PDMat)
    Ωz = similar(Ω.mat)
    d = size(Ω,1)
    ω = sqrt(diag(Ω))
    for i in 1:d, j in 1:d
        Ωz[i,j] = Ω.mat[i,j]/(ω[i]*ω[j])
    end
    Diagonal(ω), PDMat(Ωz)
end

δ(α::Float64) = α/sqrt(1.0 + α^2)

function δ(Ω::Matrix{Float64}, α::Vector{Float64})
    ω, Ωz = scale_and_cor(Ω)
    Ωzα = Ωz*α
    αΩzα = dot(α, Ωzα)
    Ωzα/sqrt(1+αΩzα)
end

function δ(Ω::PDMat, α::Vector{Float64})
    ω, Ωz = scale_and_cor(Ω)
    (Ωz * α)/sqrt(1.0 + quad(Ωz, α))
end


text:
using CellularPotts
using Test
using Random

# @testset "constructor methods" begin
#     # We want to create a model in a variety of ways
# end

# Default 2D
M = ModelParameters()

M = ModelParameters(
    graphDimension = (120,120),
    cellCounts = [200],
    cellVolumes = [50],
    penalties = [AdhesionPenalty([0 50; 50 100]), VolumePenalty(fill(50,200),[5])]
)

CPM = CellPotts(M)

MHStep!(CPM)

using BenchmarkToolsd

@benchmark MHStep!(CPM)

#3D
M = ModelParameters(
    graphDimension = (50,50,50),
    cellCounts = [50],
    cellVolumes = [1000],
    penalties = [AdhesionPenalty([0 50; 50 100]), VolumePenalty(fill(200,50),[5])])

CellGUI(CPM)


text:
module Arr

using LinearAlgebra
using Formatting: format

arange(start; kwargs...) = collect(range(start; kwargs...))

unpack_arr(x::AbstractArray) = [x[:,i] for i in range(1,stop=size(x,2))]
hstack(x) = reduce(hcat, x)
vstack(x) = reduce(vcat, x)
eye(d) = Matrix(I, d, d)

arr2str(X::Vector{T}; digits=2) where T <: Number = "[" * join(map(x->format("{:." * string(digits) * "f}", x), X), ",") * "]"
arr2str(X::Number; digits=2) = "(scalar) " * format("{:." * string(digits) * "f}", X)

@inline zero_arrays!(x::Array{T, 1}) where T <: AbstractArray = for y in x; zero_arrays!(y); end
@inline zero_arrays!(x::Array{T, 1}) where T <: Real = (x .= 0.);
@inline zero_arrays!(x::Array{T, 2}) where T <: Real = (x .= 0.);
end


text:
function update_board(;key = key[])
    token = String( trim_padding_PKCS5(decrypt( Decryptor("AES256", key), Vector{UInt8}( configuration["token"] ) )) )
    header = ["Content-Type" => "application/json", "Accept" => "application/json", "Authorization" => token]
    board_id = Glo.BoardID( configuration["board_id"] )
    cards = Glo.boards_cards( board_id, header = header )
    columns = get(Glo.boards( board_id, header = header ),"columns",Dict{String,Any}[])
    todos = find_todos()
    # TODO: sort todos by project?
    for todo in todos
        project = last( splitpath(todo.project) )
        card_ind = findfirst(c->c["name"] == todo.text, cards)
        if card_ind === nothing
            column_ind = findfirst(c->c["name"] == project, columns)
            if column_ind === nothing
                column = Glo.boards_columns(Dict("name"=>project), board_id, header = header )
                push!( columns, column )
            else
                column = columns[column_ind]
            end
            Glo.boards_cards(Dict(
                "name"=>todo.text,
                "column_id"=>column["id"],
                "description"=>("text" => repr(todo))
                ), board_id, header = header)
        else
            card = cards[card_ind]
            deleteat!(cards, card_ind)
            if get(card,"description",Dict("text"=>""))["text"] != repr(todo)
                card["description"] = "text" => repr(todo)
            end
            Glo.boards_cards(card, board_id, Glo.CardID(card["id"]), header = header)
        end
    end
    for card in cards
        Glo.boards_cards(!,board_id,Glo.CardID(card["id"]), header = header)
    end
    # TODO: detect untouched columns
    # for column in columns
    #     Glo.boards_columns(!,board_id,column["id"], header = header)
    # end
    return nothing
end # function


text:
get_params_kwargs(::GPE; kwargs...) = delete!(Dict(kwargs), :lik)
get_params_kwargs(::GPA; kwargs...) = delete!(Dict(kwargs), :noise)

"""
    optimize!(gp::GPBase, args...; kwargs...)

Optimise the hyperparameters of Gaussian process `gp` based on type II maximum likelihood estimation. This function performs gradient based optimisation using the Optim pacakge to which the user is referred to for further details.

# Keyword arguments:
    * `domean::Bool`: Mean function hyperparameters should be optmized
    * `kern::Bool`: Kernel function hyperparameters should be optmized
    * `noise::Bool`: Observation noise hyperparameter should be optimized (GPE only)
    * `lik::Bool`: Likelihood hyperparameters should be optimized (GPA only)
    * `meanbounds`: [lowerbounds, upperbounds] for the mean hyperparameters
    * `kernbounds`: [lowerbounds, upperbounds] for the kernel hyperparameters
    * `noisebounds`: [lowerbound, upperbound] for the noise hyperparameter
    * `args/kwargs`: Arguments and keyword arguments for the optimize function from the Optim package https://julianlsolvers.github.io/Optim.jl/stable/#user/config/
"""
function optimize!(gp::GPBase, method = LBFGS(), args...; domean::Bool = true, kern::Bool = true,
                   noise::Bool = true, lik::Bool = true,
                   meanbounds = nothing, kernbounds = nothing,
                   noisebounds = nothing, likbounds = nothing, kwargs...)
    params_kwargs = get_params_kwargs(gp; domean=domean, kern=kern, noise=noise, lik=lik)
    # println(params_kwargs)
    func = get_optim_target(gp; params_kwargs...)
    init = get_params(gp; params_kwargs...)  # Initial hyperparameter values
    if meanbounds == kernbounds == noisebounds == likbounds == nothing
        results = optimize(func, init, method, args...; kwargs...)     # Run optimizer
    else
        lb, ub = bounds(gp, noisebounds, meanbounds, kernbounds, likbounds;
                        domean = domean, kern = kern, noise = noise, lik = lik)
        results = optimize(func.f, func.df, lb, ub, init, Fminbox(method), args...)
    end
    set_params!(gp, Optim.minimizer(results); params_kwargs...)
    update_target!(gp)
    return results
end

function get_optim_target(gp::GPBase; params_kwargs...)
    function ltarget(hyp::AbstractVector)
        prev = get_params(gp; params_kwargs...)
        try
            set_params!(gp, hyp; params_kwargs...)
            update_target!(gp)
            return -gp.target
        catch err
            # reset parameters to remove any NaNs
            set_params!(gp, prev; params_kwargs...)

            if !all(isfinite.(hyp))
                println(err)
                return Inf
            elseif isa(err, ArgumentError)
                println(err)
                return Inf
            elseif isa(err, LinearAlgebra.PosDefException)
                println(err)
                return Inf
            else
                throw(err)
            end
        end
    end

    function ltarget_and_dltarget!(grad::AbstractVector, hyp::AbstractVector)
        prev = get_params(gp; params_kwargs...)
        try
            set_params!(gp, hyp; params_kwargs...)
            update_target_and_dtarget!(gp; params_kwargs...)
            grad[:] = -gp.dtarget
            return -gp.target
        catch err
            # reset parameters to remove any NaNs
            set_params!(gp, prev; params_kwargs...)
            if !all(isfinite.(hyp))
                println(err)
                return Inf
            elseif isa(err, ArgumentError)
                println(err)
                return Inf
            elseif isa(err, LinearAlgebra.PosDefException)
                println(err)
                return Inf
            else
                throw(err)
            end
        end
    end

    function dltarget!(grad::AbstractVector, hyp::AbstractVector)
        ltarget_and_dltarget!(grad::AbstractVector, hyp::AbstractVector)
    end

    xinit = get_params(gp; params_kwargs...)
    func = OnceDifferentiable(ltarget, dltarget!, ltarget_and_dltarget!, xinit)
    return func
end


text:
# Contains common functions for all grid types.

ZlibCompressStream(buf::IO) =
    ZlibDeflateOutputStream(buf; gzip=false, level=COMPRESSION_LEVEL)

function data_to_xml{T<:Real}(
        vtk::DatasetFile, xParent::XMLElement, data::AbstractArray{T},
        varname::AbstractString, Nc::Integer=1)
    #==========================================================================
    This variant of data_to_xml should be used when writing appended data.
      * buf is the buffer where the appended data is written.
      * xParent is the XML node under which the DataArray node will be created.
        It is either a "Points" or a "PointData" node.

    When vtk.compressed == true:
      * the data array is written in compressed form (obviously);
      * the header, written before the actual numerical data, is an array of
        UInt32 values:
            [num_blocks, blocksize, last_blocksize, compressed_blocksizes]
        All the sizes are in bytes. The header itself is not compressed, only
        the data is.
        See also:
            http://public.kitware.com/pipermail/paraview/2005-April/001391.html
            http://mathema.tician.de/what-they-dont-tell-you-about-vtk-xml-binary-formats

    Otherwise, the header is just a single UInt32 value containing the size of
    the data array in bytes.
    ==========================================================================#

    if !vtk.appended
        # Redirect to the inline version of this function.
        return data_to_xml_inline(vtk, xParent, data, varname, Nc)::XMLElement
    end

    const buf = vtk.buf    # append buffer
    const compress = vtk.compressed

    @assert name(xParent) in ("Points", "PointData", "Coordinates",
                              "Cells", "CellData")

    sType = T === Float32 ? "Float32" :
            T === Float64 ? "Float64" :
            T === Int32   ? "Int32"   :
            T === Int64   ? "Int64"   :
            T === UInt8   ? "UInt8"   :
            error("Real subtype not supported: $T")

    # DataArray node
    xDA = new_child(xParent, "DataArray")
    set_attribute(xDA, "type",   sType)
    set_attribute(xDA, "Name",   varname)
    set_attribute(xDA, "format", "appended")
    set_attribute(xDA, "offset", string(buf.position - 1))
    set_attribute(xDA, "NumberOfComponents", string(Nc))

    # Size of data array (in bytes).
    const nb = length(data) * sizeof(eltype(data))

    if compress
        const initpos = buf.position

        # Write temporary array that will be replaced later by the real header.
        header = zeros(UInt32, 4)
        write(buf, header)

        # Write compressed data.
        zWriter = ZlibCompressStream(buf)
        for val in data
            write(zWriter, val)
        end
        flush(zWriter)

        # Go back to `initpos` and write real header.
        endpos = buf.position
        compbytes = buf.position - initpos - sizeof(header)
        header[:] = [1, nb, nb, compbytes]
        buf.position = initpos
        write(buf, header)
        buf.position = endpos
    else
        write(buf, UInt32(nb))       # header (uncompressed version)
        for val in data
            write(buf, val)
        end
    end

    return xDA::XMLElement
end

function data_to_xml_inline{T<:Real}(
        vtk::DatasetFile, xParent::XMLElement, data::AbstractArray{T},
        varname::AbstractString, Nc::Integer=1)
    #==========================================================================
    This variant of data_to_xml should be used when writing data "inline" into
    the XML file (not appended at the end).

    See the other variant of this function for more info.
    ==========================================================================#

    @assert !vtk.appended
    @assert name(xParent) in ("Points", "PointData", "Coordinates",
                              "Cells", "CellData")

    const compress = vtk.compressed

    sType = T === Float32 ? "Float32" :
            T === Float64 ? "Float64" :
            T === Int32   ? "Int32"   :
            T === Int64   ? "Int64"   :
            T === UInt8   ? "UInt8"   :
            error("Real subtype not supported: $T")

    # DataArray node
    xDA = new_child(xParent, "DataArray")
    set_attribute(xDA, "type",   sType)
    set_attribute(xDA, "Name",   varname)
    set_attribute(xDA, "format", "binary")   # here, binary means base64-encoded
    set_attribute(xDA, "NumberOfComponents", "$Nc")

    # Number of bytes of data.
    const nb = length(data) * sizeof(eltype(data))

    # Write data to a buffer, which is then base64-encoded and added to the
    # XML document.
    buf = BufferedOutputStream()

    # NOTE: in the compressed case, the header and the data need to be
    # base64-encoded separately!!
    # That's why we don't use a single buffer that contains both, like in the
    # other data_to_xml function.
    if compress
        # Write compressed data.
        zWriter = ZlibCompressStream(buf)
        for val in data
            write(zWriter, val)
        end
        flush(zWriter)
    else
        for val in data
            write(buf, val)
        end
    end

    # Write buffer with data to XML document.
    add_text(xDA, "\n")
    if compress
        header = UInt32[1, nb, nb, buf.position - 1]
        add_text(xDA, base64encode(header))
    else
        add_text(xDA, base64encode(UInt32(nb)))     # header (uncompressed version)
    end
    add_text(xDA, base64encode(takebuf_string(buf)))
    add_text(xDA, "\n")

    close(buf)

    return xDA::XMLElement
end

function vtk_point_or_cell_data{T<:Real}(
        vtk::DatasetFile, data::AbstractArray{T}, name::AbstractString,
        nodetype::AbstractString, Nc::Integer)
    # Nc: number of components (Nc >= 1)

    # Find Piece node.
    xroot = root(vtk.xdoc)
    xGrid = find_element(xroot, vtk.gridType_str)
    xPiece = find_element(xGrid, "Piece")

    # Find or create "nodetype" (PointData or CellData) node.
    xtmp = find_element(xPiece, nodetype)
    xPD = (xtmp === nothing) ? new_child(xPiece, nodetype) : xtmp

    # DataArray node
    xDA = data_to_xml(vtk, xPD, data, name, Nc)

    return
end

function vtk_point_data(vtk::DatasetFile, data::AbstractArray, name::AbstractString)
    # Number of components.
    Nc = div(length(data), vtk.Npts)
    @assert Nc*vtk.Npts == length(data)
    return vtk_point_or_cell_data(vtk, data, name, "PointData", Nc)
end

function vtk_cell_data(vtk::DatasetFile, data::AbstractArray, name::AbstractString)
    # Number of components.
    Nc = div(length(data), vtk.Ncls)
    @assert Nc*vtk.Ncls == length(data)
    return vtk_point_or_cell_data(vtk, data, name, "CellData", Nc)
end

function vtk_save(vtk::DatasetFile)
    outfiles = [vtk.path]::Vector{UTF8String}
    if !isopen(vtk)
        return outfiles
    end

    if !vtk.appended
        save_file(vtk.xdoc, vtk.path)
        close(vtk)
        return outfiles
    end

    ## if vtk.appended:
    # Write XML file manually, including appended data.

    # NOTE: this is not very clean, but I can't write raw binary data with the
    # LightXML package.
    # Using raw data is way more efficient than base64 encoding in terms of
    # filesize AND time (especially time).

    # Convert XML document to a string, and split it by lines.
    lines = split(string(vtk.xdoc), '\n')

    # Verify that the last two lines are what they're supposed to be.
    @assert lines[end-1] == "</VTKFile>"
    @assert lines[end] == ""

    io = open(vtk.path, "w")

    # Write everything but the last two lines.
    for line in lines[1:end-2]
        write(io, line)
        write(io, "\n")
    end

    # Write raw data (contents of buffer vtk.buf).
    # An underscore "_" is needed before writing appended data.
    write(io, "  <AppendedData encoding=\"raw\">")
    write(io, "\n_")
    write(io, takebuf_string(vtk.buf))
    write(io, "\n  </AppendedData>")
    write(io, "\n</VTKFile>")

    close(vtk.buf)
    close(vtk)
    close(io)

    return outfiles
end

function vtk_xml_write_header(vtk::DatasetFile)
    xroot = create_root(vtk.xdoc, "VTKFile")
    set_attribute(xroot, "type", vtk.gridType_str)
    set_attribute(xroot, "version", "1.0")
    if IS_LITTLE_ENDIAN
        set_attribute(xroot, "byte_order", "LittleEndian")
    else
        set_attribute(xroot, "byte_order", "BigEndian")
    end
    if vtk.compressed
        set_attribute(xroot, "compressor", "vtkZLibDataCompressor")
        set_attribute(xroot, "header_type", "UInt32")
    end
    return xroot::XMLElement
end

# Returns the "extent" attribute required for structured (including rectilinear)
# grids.
function extent_attribute(Ni, Nj, Nk, extent::Void=nothing)
    return @sprintf("%d %d %d %d %d %d", 0, Ni-1, 0, Nj-1, 0, Nk-1)
end

function extent_attribute{T<:Integer}(Ni, Nj, Nk, extent::Array{T})
    length(extent) == 6 || error("extent must have length 6.")
    (extent[2] - extent[1] + 1 == Ni) &&
    (extent[4] - extent[3] + 1 == Nj) &&
    (extent[6] - extent[5] + 1 == Nk) ||
    error("extent is not consistent with dataset dimensions")
    ext = string(extent[1])
    for n = 2:6
        ext *= " " * string(extent[n])
    end
    return ext
end

# Number of cells in structured grids (includes structured and rectilinear
# VTK files).
# In 3D, all cells are hexahedrons (i.e. VTK_HEXAHEDRON), and the number of
# cells is (Ni-1)*(Nj-1)*(Nk-1).
# In 2D, they are quadrilaterals (VTK_QUAD), and in 1D they are line segments
# (VTK_LINE).
function num_cells_structured(Ni, Nj, Nk)
    Ncls = 1
    for N in (Ni, Nj, Nk)
        Ncls *= (N == 1) ? 1 : (N - 1)
    end
    return Ncls
end


text:
@testset "Schwefel CMA-ES" begin

    # Schwefel's ellipsoid test function:
    # a moderately conditioned ellipsoid with a dominating isolated eigenvalue
    function schwefel(x::AbstractVector{T}) where {T<:AbstractFloat}
        out = zero(T)
        for i = 1:length(x)
            out += sum(abs2, x[1:i])
        end
        return out
    end

    # Parameters
    N = 30

    # Testing: CMA-ES
    result, fitness, cnt = optimize(schwefel, CMAES(N=N; μ = 3, λ = 12),iterations = 1000)
    println("(3/3,12)-CMA-ES (Schwefel) => F: $(fitness), C: $(cnt)")

    @test ≈(result, zeros(N), atol=1e-5)
    @test ≈(fitness, 0.0, atol=1e-5)
end