text:
module TestMap
    using Base.Test
    using NullableArrays

    # create m random arrays each with N dimensions and length dims[i] along
    # dimension i for i=1:N
    N = rand(2:5)
    dims = Int[ rand(3:8) for i in 1:N ]
    m = rand(3:5)
    As = [ rand(dims...) for i in 1:m ]

    Ms = [ rand(Bool, dims...) for i in 1:m ]
    Xs = Array{NullableArray{Float64, N}, 1}()
    for i in 1:m
        push!(Xs, NullableArray(As[i]))
    end
    Ys = Array{NullableArray{Float64, N}, 1}()
    for i in 1:m
        push!(Ys, NullableArray(As[i], Ms[i]))
    end

    C = Array(Float64, dims...)
    Z = NullableArray(Float64, dims...)

    R = map(|, Ms...)

    f(x...) = sum(x)
    g(x::Float64...) = prod(x)

    dests = (C, Z)

    # 1 arg
    for (args, masks) in (
        ((As, Xs), fill((), m)), ((As, Ys), [ (Ms[i],) for i in 1:m ])
    )
        for i in 1:m
            # map!
            map!(f, args[1][i]) # map!(f, As[i])
            map!(f, args[2][i]) # map!(f, Xs[i])
            @test isequal(args[2][i], NullableArray(args[1][i], masks[i]...))
            # lifted map!
            map!(g, args[1][i])
            map!(g, args[2][i]; lift=true)
            @test isequal(args[2][i], NullableArray(args[1][i], masks[i]...))
            # map
            A = map(f, args[1][i])
            X = map(f, args[2][i])
            @test isequal(X, NullableArray(A, masks[i]...))
            # lifted map
            A = map(f, args[1][i])
            X = map(g, args[2][i]; lift=true)
            @test isequal(X, NullableArray(A, masks[i]...))
        end
    end
    # 2 arg
    i, j = rand(1:m), rand(1:m)
    S = map(|, Ms[i], Ms[j])
    for (args, mask) in (
        ((As, Xs), ()), ((As, Ys), (S,))
    )
        # map!
        map!(f, dests[1], args[1][i], args[1][j])
        map!(f, dests[2], args[2][i], args[2][j])
        @test isequal(dests[2], NullableArray(dests[1], mask...))
        # lifted map!
        map!(g, dests[1], args[1][i], args[1][j])
        map!(g, dests[2], args[2][i], args[2][j]; lift=true)
        @test isequal(dests[2], NullableArray(dests[1], mask...))
        # map
        map(f, args[1][i], args[1][j])
        map(f, args[2][i], args[2][j])
        @test isequal(dests[2], NullableArray(dests[1], mask...))
        # lifted map
        map(g, args[1][i], args[1][j])
        map(g, args[2][i], args[2][j]; lift=true)
        @test isequal(dests[2], NullableArray(dests[1], mask...))
    end
    # n arg
    for (args, mask) in (
        ((As, Xs), ()), ((As, Ys), (R,))
    )
        # map!
        map!(f, dests[1], args[1]...)
        map!(f, dests[2], args[2]...)
        @test isequal(dests[2], NullableArray(dests[1], mask...))
        # lifted map!
        map!(g, dests[1], args[1]...)
        map!(g, dests[2], args[2]...; lift=true)
        @test isequal(dests[2], NullableArray(dests[1], mask...))
        # map
        map(f, args[1]...)
        map(f, args[2]...)
        @test isequal(dests[2], NullableArray(dests[1], mask...))
        # lifted map
        map(g, args[1]...)
        map(g, args[2]...; lift=true)
        @test isequal(dests[2], NullableArray(dests[1], mask...))
    end

    # test map over empty NullableArrays
    X = NullableArray(Int[])
    h1(x) = 5.0*x
    h2(x::Nullable) = x
    h2(x::Nullable...) = prod(x)

    Z1 = map(h1, X)
    Z2 = map(h2, X)
    if VERSION >= v"0.5.0-dev+3294"
        @test isequal(Z1, NullableArray(Union{}[])) # no type inference
    end
    if VERSION >= v"0.5.0-dev+5297"
        @test isequal(Z2, NullableArray(Float64[])) # use type inference
    end

    Z1 = map(h1, X; lift=true) # use type inference
    @test isequal(Z1, NullableArray(Int64[]))

    # if a function has no method for inner_eltype of empty NullableArray,
    # result should be empty NullableArray{Union{}}()
    h3(x::Float64...) = prod(x)
    Z3 = map(h3, X)
    @test isequal(Z3, NullableArray(Union{}[]))
    Z3 = map(h3, X, X)
    @test isequal(Z3, NullableArray(Union{}[]))
    Z3 = map(h3, X, X, X)
    @test isequal(Z3, NullableArray(Union{}[]))
    Z3 = map(h3, X; lift=true)
    @test isequal(Z3, NullableArray(Union{}[]))
    Z3 = map(h3, X, X; lift=true)
    @test isequal(Z3, NullableArray(Union{}[]))
    Z3 = map(h3, X, X, X; lift=true)
    @test isequal(Z3, NullableArray(Union{}[]))

    # test map over all null NullableArray
    n = rand(10:100)
    Ys = [ NullableArray(rand(Int, n), fill(true, n)) for i in 1:rand(3:5) ]

    Z2 = map(h2, Ys[1]) # behavior should be determined solely by h2, i.e. no type inference
    @test isequal(Z2, NullableArray(Int, n))
    Z2 = map(h2, Ys[1], Ys[2])
    @test isequal(Z2, NullableArray(Int, n))
    Z2 = map(h2, Ys...)
    @test isequal(Z2, NullableArray(Int, n))

    Z1 = map(h1, Ys[1]; lift=true) # fall back on type inference
    @test isequal(Z1, NullableArray(Float64, n))
    Z1 = map(h1, Ys[1], Ys[2]; lift=true)
    @test isequal(Z1, NullableArray(Float64, n))
    Z1 = map(h1, Ys...; lift=true)
    @test isequal(Z1, NullableArray(Float64, n))

    # if a function has no method for inner_eltype of all null NullableArray,
    # result of lifted map should be NullableArray{Union{}}()
    h4(x::Int...) = prod(x)
    Z4 = map(h4, Ys[1]; lift=true)
    @test isequal(Z4, NullableArray(Union{}, n))
    Z4 = map(h4, Ys[1], Ys[2]; lift=true)
    @test isequal(Z4, NullableArray(Union{}, n))
    Z4 = map(h4, Ys...; lift=true)
    @test isequal(Z4, NullableArray(Union{}, n))
end # module TestMap


text:
using BinaryBuilder

cuda_version = v"9.0.176"   # NOTE: could be less specific

sources_linux = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/7.6.5.32/Production/9.0_20191031/cudnn-9.0-linux-x64-v7.6.5.32.tgz",
                  "bd0a4c0090d5b02feec3f195738968690cc2470b9bc6026e6fe8ff245cd261c8")
]
sources_windows = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/7.6.5.32/Production/9.0_20191031/cudnn-9.0-windows10-x64-v7.6.5.32.zip",
                  "c7401514a6d7d24e8541f88c12e4328f165b5c5afd010ee462d356cac2158268")
]

include("../common.jl")

ext: jl


text:
using LinearAlgebra

mutable struct Nystrom{T}
	Ω::AbstractMatrix{T}
	S::AbstractMatrix{T}
	function Nystrom{T}(n, R) where {T}
		new(randn(T, n, R), zeros(n, R))
	end
end

function update!(sketch, v, η)
	sketch.S .*= (1 - η)
	sketch.S .+= η * v * (v' * sketch.Ω)
end

function reconstruct(sketch)
	n = size(sketch.S, 1)

	σ = sqrt(n) * eps() * norm(sketch.S)
	Sσ = sketch.S + σ * sketch.Ω
	L = cholesky(Hermitian(sketch.Ω' * Sσ); check=false)
	U, Σ, = svd(Sσ / L.U)
	Λ = max.(0, Σ .^ 2 .- σ)

	U, Diagonal(Λ)
end


text:
# @doc """A DReal context manages all other DReal objects, global configuration options, etc.
#     Objects created in one context cannot be used in another one.""" ->
type Context
  ctx::Ptr{Void}
  vars::Set{String}
end

Context(ctx::Ptr{Void}) = Context(ctx,Set{String}())
Context(l::Logic) = Context(opensmt_mk_context(@compat UInt32(l.val)))

add_vars!(ctx::Context, v::Set{String}) = union!(ctx.vars,v)
add_vars!(ctx::Context, v::String) = push!(ctx.vars,v)

create_global_ctx!(l::Logic = qf_nra) =
  (global default_global_context; default_global_context = Context(l))
global_context() = (global default_global_context; default_global_context)

function set_global_ctx!(ctx::Context)
  global default_global_context
  default_global_context = ctx
end

@doc """push creates a new scope by saving the current stack size.
  A `pop` following push will undo all assertions in between""" ->
push_ctx!(ctx::Context) = opensmt_push(ctx.ctx)
push_ctx!() = push_ctx!(global_context())

@doc "Pop removes any assertion or declaration performed between it and the matching push." ->
pop_ctx!(ctx::Context) = opensmt_pop(ctx.ctx)
pop_ctx!() = pop_ctx!(global_context())

reset_ctx!(ctx::Context) = (opensmt_reset(ctx.ctx);ctx.vars = Set{String}())
reset_ctx!() = reset_ctx!(global_context())

delete_ctx!(ctx::Context) = opensmt_del_context(ctx.ctx)
delete_ctx!() = delete_ctx!(global_context())

function reset_global_ctx!(l::Logic = qf_nra)
  delete_ctx!(global_context())
  create_global_ctx!(l)
end


text:
using Checkpointing

function main(steps, checkpoints; verbose=0)
    store = function f() end
    revolve = Revolve(steps, checkpoints, store, store; verbose=verbose)
    guessed_checkpoints = guess(revolve)
    @show guessed_checkpoints
    println("Revolve suggests : $guessed_checkpoints checkpoints for a factor of $(factor(revolve, steps,guessed_checkpoints))")
    println("actions: ")
    while true
        next_action = next_action!(revolve)
        if next_action.actionflag == Checkpointing.done
            break
        end
    end
    return revolve
end


text:
"""
    HMCSampling(;ϵ::T=1e-5,nBurnin::Int=100,thinning::Int=10)

Draw samples from the true posterior via Hamiltonian Monte Carlo.

## Keywords arguments
- `ϵ::T` : convergence criteria
- `nBurnin::Int` : Number of samples discarded before starting to save samples
- `thinning::Int` : Frequency of sampling
"""
mutable struct HMCSampling{T<:Real,N,Tx,Ty} <: SamplingInference{T}
    nBurnin::Int # Number of burnin samples
    thinning::Integer # Frequency at which samples are saved
    ϵ::T #Convergence criteria
    n_iter::Int #Number of samples computed
    HyperParametersUpdated::Bool #To know if the inverse kernel matrix must updated
    opt::NTuple{N,SOptimizer}
    xview::Tx
    yview::Ty
    sample_store::Array{T,3}
    function HMCSampling{T}(nBurnin::Int, thinning::Int, ϵ::Real) where {T}
        nBurnin >= 0 || error("nBurnin should be positive")
        thinning >= 0 || error("thinning should be positive")
        return new{T,1,Vector{T},Vector{T}}(nBurnin, thinning, ϵ)
    end
    function HMCSampling{T}(
        nBurnin::Int,
        thinning::Int,
        ϵ::Real,
        nFeatures::Vector{<:Int},
        nSamples::Int,
        nLatent::Int,
        xview::Tx,
        yview::Ty,
    ) where {T,Tx,Ty}
        opts = ntuple(_ -> SOptimizer{T}(nothing), nLatent)
        return new{T,nLatent,Tx,Ty}(
            nBurnin, thinning, ϵ, 0, nSamples, true, opts, xview, yview
        )
    end
end

ρ(::HMCSampling{T}) where {T} = one(T)
nMinibatch(i::HMCSampling) = i.nSamples

function HMCSampling(; ϵ::T=1e-5, nBurnin::Int=100, thinning::Int=10) where {T<:Real}
    return HMCSampling{Float64}(nBurnin, thinning, ϵ)
end

function Base.show(io::IO, inference::HMCSampling{T}) where {T<:Real}
    return print(io, "Hamilton Monte Carlo Sampler")
end

function tuple_inference(
    i::HMCSampling{T},
    nLatent::Int,
    nFeatures::Vector{<:Int},
    nSamples::Int,
    nMinibatch::Int, #unused
    xview,
    yview,
) where {T}
    return HMCSampling{T}(
        i.nBurnin, i.thinning, i.ϵ, nFeatures, nSamples, nLatent, xview, yview
    )
end

function grad_log_joint_pdf(model::MCGP{T,L,HMCSampling{T}}, x) where {T,L}
    return grad_log_likelihood(likelihood(model), yview(model), x) +
           sum(grad_log_gp_prior.(model.f, x))
end

function init_sampler!(
    inference::HMCSampling{T},
    nLatent::Integer,
    nFeatures::Integer,
    nSamples::Integer,
    cat_samples::Bool,
) where {T<:Real}
    if n_iter(inference) == 0 || !cat_samples
        inference.sample_store = zeros(T, nSamples, nFeatures, nLatent)
    else
        inference.sample_store = cat(
            inference.sample_store, zeros(T, nSamples, nFeatures, nLatent); dims=1
        )
    end
    return inference
end

function sample_parameters(
    model::MCGP{T,L,<:HMCSampling{T}}, nSamples::Int, callback, cat_samples::Bool
) where {T,L}
    f_init = vcat(means(model)...)
    logpdf(x) = log_joint_model(model, x)
    gradlogpdf(x) = (log_joint_model(model, x), grad_log_joint_model(model, x))
    metric = DiagEuclideanMetric(length(f_init))
    h = Hamiltonian(metric, logpdf, gradlogpdf)
    int = Leapfrog(find_good_eps(h, f_init))
    prop = NUTS{MultinomialTS,GeneralisedNoUTurn}(int)
    adaptor = StanHMCAdaptor(
        n_adapts, Preconditioner(metric), NesterovDualAveraging(0.8, int.ϵ)
    )
    return samples, stats = sample(
        h, prop, f_init, nSamples, adaptor, n_adapts; progress=true
    )
end


text:
#=

Author: Shunsuke Hori

=#

using PyPlot

"""
If we expand the implicit function for h_j(n_k) then we find that
it is a quadratic. We know that h_j(n_k) > 0 so we can get its
value by using the quadratic form
"""

function h_j(j::Integer, nk::Real, s1::Real, s2::Real,
             θ::Real, δ::Real, ρ::Real)
    # Find out who's h we are evaluating
    if j == 1
        sj = s1
        sk = s2
    else
        sj = s2
        sk = s1
    end

    # Coefficients on the quadratic a x^2 + b x + c = 0
    a = 1.0
    b = ((ρ + 1 / ρ) * nk - sj - sk)
    c = (nk * nk - (sj * nk) / ρ - sk * ρ * nk)

    # Positive solution of quadratic form
    root = (-b + sqrt(b * b - 4 * a * c)) / (2 * a)

    return root
end

"""
Determine whether (n1, n2) is in the set DLL
"""
DLL(n1::Real, n2::Real,
    s1_ρ::Real, s2_ρ::Real,
    s1::Real, s2::Real,
    θ::Real, δ::Real, ρ::Real) =
        (n1 <= s1_ρ) && (n2 <= s2_ρ)

"""
Determine whether (n1, n2) is in the set DHH
"""
DHH(n1::Real, n2::Real,
    s1_ρ::Real, s2_ρ::Real,
    s1::Real, s2::Real,
    θ::Real, δ::Real, ρ::Real) =
        (n1 >= h_j(1, n2, s1, s2, θ, δ, ρ)) && (n2 >= h_j(2, n1, s1, s2, θ, δ, ρ))

"""
Determine whether (n1, n2) is in the set DHL
"""
DHL(n1::Real, n2::Real,
    s1_ρ::Real, s2_ρ::Real,
    s1::Real, s2::Real,
    θ::Real, δ::Real, ρ::Real) =
        (n1 >= s1_ρ) && (n2 <= h_j(2, n1, s1, s2, θ, δ, ρ))

"""
Determine whether (n1, n2) is in the set DLH
"""
DLH(n1::Real, n2::Real,
    s1_ρ::Real, s2_ρ::Real,
    s1::Real, s2::Real,
    θ::Real, δ::Real, ρ::Real) =
        (n1 <= h_j(1, n2, s1, s2, θ, δ, ρ)) && (n2 >= s2_ρ)

"""
Takes a current value for (n_{1, t}, n_{2, t}) and returns the
values (n_{1, t+1}, n_{2, t+1}) according to the law of motion.
"""
function one_step(n1::Real, n2::Real,
                  s1_ρ::Real, s2_ρ::Real,
                  s1::Real, s2::Real,
                  θ::Real, δ::Real, ρ::Real)
    # Depending on where we are, evaluate the right branch
    if DLL(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ)
        n1_tp1 = δ * (θ * s1_ρ + (1 - θ) * n1)
        n2_tp1 = δ * (θ * s2_ρ + (1 - θ) * n2)
    elseif DHH(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ)
        n1_tp1 = δ * n1
        n2_tp1 = δ * n2
    elseif DHL(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ)
        n1_tp1 = δ * n1
        n2_tp1 = δ * (θ * h_j(2, n1, s1, s2, θ, δ, ρ) + (1 - θ) * n2)
    elseif DLH(n1, n2, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ)
        n1_tp1 = δ * (θ * h_j(1, n2, s1, s2, θ, δ, ρ) + (1 - θ) * n1)
        n2_tp1 = δ * n2
    end

    return n1_tp1, n2_tp1
end

"""
Given an initial condition, continues to yield new values of `n1` and `n2`
"""
new_n1n2(n1_0::Real, n2_0::Real,
         s1_ρ::Real, s2_ρ::Real,
         s1::Real, s2::Real,
         θ::Real, δ::Real, ρ::Real) =
            one_step(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ)

"""
Takes initial values and iterates forward to see whether
the histories eventually end up in sync.

If countries are symmetric then as soon as the two countries have the
same measure of firms then they will by synchronized -- However, if
they are not symmetric then it is possible they have the same measure
of firms but are not yet synchronized. To address this, we check whether
firms stay synchronized for `npers` periods with Euclidean norm

##### Parameters
----------
- `n1_0` : `Real`,
Initial normalized measure of firms in country one
- `n2_0` : `Real`,
Initial normalized measure of firms in country two
- `maxiter` : `Integer`,
Maximum number of periods to simulate
- `npers` : `Integer`,
Number of periods we would like the countries to have the same measure for

##### Returns
-------
- `synchronized` : `Bool`,
Did they two economies end up synchronized
- `pers_2_sync` : `Integer`,
The number of periods required until they synchronized
"""
function pers_till_sync(n1_0::Real, n2_0::Real,
                        s1_ρ::Real, s2_ρ::Real,
                        s1::Real, s2::Real,
                        θ::Real, δ::Real, ρ::Real,
                        maxiter::Integer, npers::Integer)

    # Initialize the status of synchronization
    synchronized = false
    pers_2_sync = maxiter
    iters = 0

    nsync = 0

    while (~synchronized) && (iters < maxiter)
        # Increment the number of iterations and get next values
        iters += 1

        n1_t, n2_t = new_n1n2(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ)

        # Check whether same in this period
        if abs(n1_t - n2_t) < 1e-8
            nsync += 1
        # If not, then reset the nsync counter
        else
            nsync = 0
        end

        # If we have been in sync for npers then stop and countries
        # became synchronized nsync periods ago
        if nsync > npers
            synchronized = true
            pers_2_sync = iters - nsync
        end
        n1_0, n2_0 = n1_t, n2_t
    end
    return synchronized, pers_2_sync
end

function create_attraction_basis{TR <: Real}(s1_ρ::TR, s2_ρ::TR,
                                 s1::TR, s2::TR, θ::TR, δ::TR, ρ::TR,
                                 maxiter::Integer, npers::Integer, npts::Integer)
    # Create unit range with npts
    synchronized, pers_2_sync = false, 0
    unit_range = linspace(0.0, 1.0, npts)

    # Allocate space to store time to sync
    time_2_sync = Matrix{TR}(npts, npts)
    # Iterate over initial conditions
    for (i, n1_0) in enumerate(unit_range)
        for (j, n2_0) in enumerate(unit_range)
            synchronized, pers_2_sync = pers_till_sync(n1_0, n2_0, s1_ρ, s2_ρ,
                                                        s1, s2, θ, δ, ρ,
                                                        maxiter, npers)
            time_2_sync[i, j] = pers_2_sync
        end
    end

    return time_2_sync
end

# == Now we define a type for the model == #

"""
The paper "Globalization and Synchronization of Innovation Cycles" presents
a two country model with endogenous innovation cycles. Combines elements
from Deneckere Judd (1985) and Helpman Krugman (1985) to allow for a
model with trade that has firms who can introduce new varieties into
the economy.

We focus on being able to determine whether two countries eventually
synchronize their innovation cycles. To do this, we only need a few
of the many parameters. In particular, we need the parameters listed
below

##### Parameters
----------
- `s1` : `Real`,
Amount of total labor in country 1 relative to total worldwide labor
- `θ` : `Real`,
A measure of how mcuh more of the competitive variety is used in
production of final goods
- `δ` : `Real`,
Percentage of firms that are not exogenously destroyed every period
- `ρ` : `Real`,
Measure of how expensive it is to trade between countries
"""
struct MSGSync{TR <: Real}
    s1::TR
    s2::TR
    s1_ρ::TR
    s2_ρ::TR
    θ::TR
    δ::TR
    ρ::TR
end

function MSGSync(s1::Real=0.5, θ::Real=2.5,
                 δ::Real=0.7, ρ::Real=0.2)
    # Store other cutoffs and parameters we use
    s2 = 1 - s1
    s1_ρ = min((s1 - ρ * s2) / (1 - ρ), 1)
    s2_ρ = 1 - s1_ρ

    model=MSGSync(s1, s2, s1_ρ, s2_ρ, θ, δ, ρ)

    return model
end

unpack_params(model::MSGSync) =
    model.s1, model.s2, model.θ, model.δ, model.ρ, model.s1_ρ, model.s2_ρ

"""
Simulates the values of `n1` and `n2` for `T` periods

##### Parameters
----------
- `n1_0` : `Real`, Initial normalized measure of firms in country one
- `n2_0` : `Real`, Initial normalized measure of firms in country two
- `T` : `Integer`, Number of periods to simulate

##### Returns
-------
- `n1` : `Vector{TR}(ndim=1) where TR <: Real`,
A history of normalized measures of firms in country one
- `n2` : `Vector{TR}(ndim=1) where TR <: Real`,
A history of normalized measures of firms in country two
"""
function simulate_n{TR <: Real}(model::MSGSync, n1_0::TR, n2_0::TR, T::Integer)
    # Unpack parameters
    s1, s2, θ, δ, ρ, s1_ρ, s2_ρ = unpack_params(model)

    # Allocate space
    n1 = Vector{TR}(T)
    n2 = Vector{TR}(T)

    # Simulate for T periods
    for t in 1:T
        # Get next values
        n1[t], n2[t] = n1_0, n2_0
        n1_0, n2_0 = new_n1n2(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2, θ, δ, ρ)
    end

    return n1, n2
end

"""
Takes initial values and iterates forward to see whether
the histories eventually end up in sync.

If countries are symmetric then as soon as the two countries have the
same measure of firms then they will by synchronized -- However, if
they are not symmetric then it is possible they have the same measure
of firms but are not yet synchronized. To address this, we check whether
firms stay synchronized for `npers` periods with Euclidean norm

##### Parameters
----------
- `n1_0` : `Real`,
Initial normalized measure of firms in country one
- `n2_0` : `Real`,
Initial normalized measure of firms in country two
- `maxiter` : `Integer`,
Maximum number of periods to simulate
- `npers` : `Integer`,
Number of periods we would like the countries to have the same measure for

##### Returns
-------
- `synchronized` : `Bool`,
Did they two economies end up synchronized
- `pers_2_sync` : `Integer`,
The number of periods required until they synchronized
"""
function pers_till_sync(model::MSGSync, n1_0::Real, n2_0::Real,
                        maxiter::Integer=500, npers::Integer=3)
    # Unpack parameters
    s1, s2, θ, δ, ρ, s1_ρ, s2_ρ = unpack_params(model)

    return pers_till_sync(n1_0, n2_0, s1_ρ, s2_ρ, s1, s2,
                          θ, δ, ρ, maxiter, npers)
end

"""
Creates an attraction basis for values of n on [0, 1] X [0, 1] with npts in each dimension
"""
function create_attraction_basis(model::MSGSync;
                                 maxiter::Integer=250,
                                 npers::Integer=3,
                                 npts::Integer=50)
    # Unpack parameters
    s1, s2, θ, δ, ρ, s1_ρ, s2_ρ = unpack_params(model)

    ab = create_attraction_basis(s1_ρ, s2_ρ, s1, s2, θ, δ,
                                  ρ, maxiter, npers, npts)
    return ab
end


text:
@testset "load lcmc w/ no errors" begin
    lcmc_cf = charfreq(SimplifiedLCMC())
    @test length(lcmc_cf) != 0
end

@testset "load junda w/ no errors" begin
    junda_cf = charfreq(SimplifiedJunDa())
    @test length(junda_cf) != 0
end


text:
# Julia interface when not defined in source files

Base.eltype(x::DistMultiVec{T}) where {T} = T

function Base.size(A::ElementalMatrix, i::Integer)
    if i < 1
        error("dimension out of range")
    elseif i == 1
        return Int(height(A))
    elseif i == 2
        return Int(width(A))
    else
        return 1
    end
end

Base.size(A::ElementalMatrix) = (size(A, 1), size(A, 2))

Base.copyto!(A::T, B::T) where {T<:ElementalMatrix} = _copy!(B, A)
# copy(A::ElementalMatrix) = copy!(similar(A), A)
Base.length(A::ElementalMatrix) = prod(size(A))

## Current mutating Julia multiplication API
LinearAlgebra.mul!(C::T, A::T, B::T) where {T<:ElementalMatrix} = gemm(NORMAL, NORMAL, one(eltype(T)), A, B, zero(eltype(T)), C)
LinearAlgebra.mul!(C::T, adjA::Adjoint{<:Any,T}, B::T) where {T<:ElementalMatrix} = gemm(ADJOINT, NORMAL, one(eltype(T)), parent(adjA), B, zero(eltype(T)), C)
LinearAlgebra.mul!(C::T, trA::Transpose{<:Any,T}, B::T) where {T<:ElementalMatrix} = gemm(TRANSPOSE, NORMAL, one(eltype(T)), parent(trA), B, zero(eltype(T)), C)
LinearAlgebra.mul!(C::T, A::T, adjB::Adjoint{<:Any,T}) where {T<:ElementalMatrix} = gemm(NORMAL, ADJOINT, one(eltype(T)), A, parent(adjB), zero(eltype(T)), C)
LinearAlgebra.mul!(C::T, A::T, trB::Transpose{<:Any,T}) where {T<:ElementalMatrix} = gemm(NORMAL, TRANSPOSE, one(eltype(T)), A, parent(trB), zero(eltype(T)), C)

## BLAS like multiplication API (i.e. with α and β)
LinearAlgebra.mul!(C::S, A::S, B::S, α::Number, β::Number) where {S<:ElementalMatrix{T}} where {T} =
    gemm(NORMAL, NORMAL, T(α), A, B, T(β), C)
LinearAlgebra.mul!(C::S, adjA::Adjoint{<:Any,S}, B::S, α::Number, β::Number) where {S<:ElementalMatrix{T}} where {T} =
    gemm(ADJOINT, NORMAL, T(α), parent(adjA), B, T(β), C)
LinearAlgebra.mul!(C::S, trA::Transpose{<:Any,S}, B::S, α::Number, β::Number) where {S<:ElementalMatrix{T}} where {T} =
    gemm(TRANSPOSE, NORMAL, T(α), parent(trA), B, T(β), C)
LinearAlgebra.mul!(C::S, A::S, adjB::Adjoint{<:Any,S}, α::Number, β::Number) where {S<:ElementalMatrix{T}} where {T} =
    gemm(NORMAL, ADJOINT, T(α), A, parent(adjB), T(β), C)
LinearAlgebra.mul!(C::S, A::S, trB::Transpose{<:Any,S}, α::Number, β::Number) where {S<:ElementalMatrix{T}} where {T} =
    gemm(NORMAL, TRANSPOSE, T(α), A, parent(trB), T(β), C)

## Linear solve API
LinearAlgebra.ldiv!(A::LowerTriangular{T,S}, B::S)  where {T,S<:ElementalMatrix} =
    trsm(LEFT, LOWER, NORMAL   , NON_UNIT, one(T), A.data, B)
LinearAlgebra.ldiv!(adjA::Adjoint{T,LowerTriangular{T,S}}, B::S) where {T,S<:ElementalMatrix} =
    trsm(LEFT, LOWER, ADJOINT  , NON_UNIT, one(T), parent(adjA).data, B)
LinearAlgebra.ldiv!(trA::Transpose{T,LowerTriangular{T,S}}, B::S) where {T,S<:ElementalMatrix} =
    trsm(LEFT, LOWER, TRANSPOSE, NON_UNIT, one(T), parent(trA).data, B)
LinearAlgebra.ldiv!(A::UpperTriangular{T,S}, B::S)  where {T,S<:ElementalMatrix} =
    trsm(LEFT, UPPER, NORMAL   , NON_UNIT, one(T), A.data, B)
LinearAlgebra.ldiv!(adjA::Adjoint{T,UpperTriangular{T,S}}, B::S) where {T,S<:ElementalMatrix} =
    trsm(LEFT, UPPER, ADJOINT  , NON_UNIT, one(T), parent(adjA).data, B)
LinearAlgebra.ldiv!(trA::Transpose{T,UpperTriangular{T,S}}, B::S) where {T,S<:ElementalMatrix} =
    trsm(LEFT, UPPER, TRANSPOSE, NON_UNIT, one(T), parent(trA).data, B)

LinearAlgebra.rdiv!(A::S, B::LowerTriangular{T,S})  where {T,S<:ElementalMatrix} =
    trsm(RIGHT, LOWER, NORMAL   , NON_UNIT, one(T), B.data, A)
LinearAlgebra.rdiv!(A::S, adjB::Adjoint{T,LowerTriangular{T,S}}) where {T,S<:ElementalMatrix} =
    trsm(RIGHT, LOWER, ADJOINT  , NON_UNIT, one(T), parent(adjB).data, A)
LinearAlgebra.rdiv!(A::S, trB::Transpose{T,LowerTriangular{T,S}}) where {T,S<:ElementalMatrix} =
    trsm(RIGHT, LOWER, TRANSPOSE, NON_UNIT, one(T), B.data, A)
LinearAlgebra.rdiv!(A::S, B::UpperTriangular{T,S})  where {T,S<:ElementalMatrix} =
    trsm(RIGHT, UPPER, NORMAL   , NON_UNIT, one(T), B.data, A)
LinearAlgebra.rdiv!(A::S, adjB::Adjoint{T,UpperTriangular{T,S}}) where {T,S<:ElementalMatrix} =
    trsm(RIGHT, UPPER, ADJOINT  , NON_UNIT, one(T), parent(adjB).data, A)
LinearAlgebra.rdiv!(A::S, trB::Transpose{T,UpperTriangular{T,S}}) where {T,S<:ElementalMatrix} =
    trsm(RIGHT, UPPER, TRANSPOSE, NON_UNIT, one(T), parent(trB).data, A)

# Spectral
LinearAlgebra.svd(A::ElementalMatrix) = svd!(copy(A))
LinearAlgebra.svd(A::ElementalMatrix, ctrl::SVDCtrl) = svd!(copy(A), ctrl)
LinearAlgebra.svdvals(A::ElementalMatrix, ctrl::SVDCtrl) = svdvals!(copy(A), ctrl)

# conversions to and from julia arrays

# function copy!(dest::Matrix{T}, src::Base.VecOrMat{T}) where {T}
#     m, n = size(src, 1), size(src, 2)
#     resize!(dest, m, n)
#     Base.unsafe_copy!(pointer(dest), pointer(src), m*n)
#     return dest
# end
# function copy!(dest::Base.Matrix{T}, src::Matrix{T}) where {T}
#     m, n = size(dest)
#     if m != size(src, 1) || n != size(src, 2)
#         throw(DimensionMisMatch("source and destination must have same shape"))
#     end
#     Base.unsafe_copy!(pointer(dest), pointer(src), m*n)
#     return dest
# end

function Base.copyto!(dest::DistMatrix{T}, src::Base.AbstractVecOrMat) where {T}
    m, n = size(src, 1), size(src, 2)
    zeros!(dest, m, n)
    if MPI.commRank(comm(dest)) == 0
        for j = 1:n
            for i = 1:m
                queueUpdate(dest, i, j, src[i,j])
            end
        end
    end
    processQueues(dest)
    return dest
end

Base.copyto!(dest::DistMatrix, src::ElementalMatrix) = _copy!(src, dest)

function Base.copyto!(dest::Base.VecOrMat, src::DistMatrix{T}) where {T}
    m, n = size(src, 1), size(src, 2)
    # if MPI.commRank(comm(src)) == 0
        for j = 1:n
            for i = 1:m
                queuePull(src, i, j)
            end
        end
    # end
    dest_mat = ndims(dest) == 1 ? reshape(dest, :, 1) : dest
    processPullQueue(src, dest_mat)
    return dest
end

function Base.convert(::Type{Matrix{T}}, A::Base.VecOrMat{T}) where {T}
    m, n = size(A, 1), size(A, 2)
    B = Matrix(T)
    resize!(B, m, n)
    Base.unsafe_copyto!(pointer(B), pointer(A), m*n)
    return B
end
function Base.convert(::Type{Base.Matrix{T}}, A::Matrix{T}) where {T}
    m, n = size(A)
    B = Base.Matrix{T}(undef, m, n)
    Base.unsafe_copyto!(pointer(B), pointer(A), m*n)
    return B
end

function Base.convert(::Type{DistMatrix{T}}, A::Base.VecOrMat{T}) where {T}
    m, n = size(A, 1), size(A, 2)
    B = DistMatrix(T)
    zeros!(B, m, n)
    if MPI.commRank(comm(B)) == 0
        for j = 1:n
            for i = 1:m
                queueUpdate(B, i, j, A[i,j])
            end
        end
    end
    processQueues(B)
    return B
end

function Base.convert(::Type{DistMultiVec{T}}, A::Base.VecOrMat{T}) where {T}
    m, n = size(A, 1), size(A, 2)
    B = DistMultiVec(T)
    zeros!(B, m, n)
    if MPI.commRank(comm(B)) == 0
        for j = 1:n
            for i = 1:m
                queueUpdate(B, i, j, A[i,j])
            end
        end
    end
    processQueues(B)
    return B
end

function Base.convert(::Type{DistMatrix{T}}, A::DistMultiVec{T}) where {T}
    m, n = size(A)
    B = DistMatrix(T)
    zeros!(B, m, n)
    for j = 1:n
        for i = 1:localHeight(A)
            xij = getLocal(A, i, j)
            queueUpdate(B, globalRow(A, i), j, xij)
        end
    end
    processQueues(B)
    return B
end

Base.convert(::Type{Array}, xd::DistMatrix{T}) where {T} =
    copyto!(Base.zeros(T, size(xd)), xd)

Base.Array(xd::DistMatrix) = convert(Array, xd)

LinearAlgebra.norm(x::ElementalMatrix) = nrm2(x)
# function LinearAlgebra.norm(x::ElementalMatrix)
#     if size(x, 2) == 1
#         return nrm2(x)
#     else
#         return twoNorm(x)
#     end
# end

LinearAlgebra.cholesky!(A::Hermitian{<:Union{Real,Complex},<:ElementalMatrix}) = LinearAlgebra.Cholesky(_cholesky!(A.uplo == 'U' ? UPPER : LOWER, A.data), A.uplo, 0)
LinearAlgebra.cholesky(A::Hermitian{<:Union{Real,Complex},<:ElementalMatrix}) = cholesky!(copy(A))

LinearAlgebra.lu(A::ElementalMatrix) = _lu!(copy(A))

# Mixed multiplication with Julia Arrays
(*)(A::DistMatrix{T}, B::StridedVecOrMat{T}) where {T} = A*convert(DistMatrix{T}, B)
(*)(A::DistMultiVec{T}, B::StridedVecOrMat{T}) where {T} = convert(DistMatrix{T}, A)*convert(DistMatrix{T}, B)
(*)(A::DistSparseMatrix{T}, B::StridedVecOrMat{T}) where {T} = A*convert(DistMultiVec{T}, B)
(*)(A::Adjoint{T,DistMatrix{T}}, B::StridedVecOrMat{T}) where {T} = A*convert(DistMatrix{T}, B)
(*)(A::Adjoint{T,DistMultiVec{T}}, B::Base.VecOrMat{T}) where {T} = convert(DistMatrix{T}, parent(A))'*convert(DistMatrix{T}, B)
(*)(A::Adjoint{T,DistSparseMatrix{T}}, B::Base.VecOrMat{T}) where {T} = A*convert(DistMultiVec{T}, B)


text:
###Equations
eq=equation(:x)
@test eq==Equation(:x,0)

eq1=:a≖:b*:c;eq2=:c≖:d*:e;
eq3=eq1&eq2
@test eq3==(:a≖:b*:d*:e)

energy=@equ E=m*c^2
c=@equ c=299792458
m=@equ m=3*n
n=@equ n=9
@test (energy&c&m&n).rhs==2426638982589407628

tri=@equ c^2=a^2+b^2
@test randeval(sqrt(tri).rhs)==randeval(sqrt(:a^2+:b^2))

rel=@equ Oneable(a)*x*z=y
@test (:q*:t)&rel==:y

#evaluate
eq=equation(1/:x-sqrt(:x))
@test evaluate(eq,Dict(:x=>1))==(0,0)
m=matches(eq)
for tm in m
	l,r=evaluate(tm,Dict(:x=>1))
	@test l==r
end

#matches
ex=3*:x^2-5*:x+1.5
m=matches(ex,quadratic)
@test evaluate(ex,Dict(:x=>m[1].rhs))<1e-9 && evaluate(ex,Dict(:x=>m[2].rhs))<1e-9
ex=:x*3*:x-5*:x*3+1.5
m=matches(ex,quadratic)
@test evaluate(ex,Dict(:x=>m[1].rhs))<1e-9 && evaluate(ex,Dict(:x=>m[2].rhs))<1e-9

#solve
x=@equ(x+3=0)&:x
@test x.lhs==:x
@test x.rhs==-3
eq=equation(:x+:z+:t)
sol=solve(eq)
for s in sol
	res=evaluate(eq,Dict(s.lhs=>s.rhs))
	@test res[1]==res[2]
end
eq=simplify(@equ m*V=m*v1+M*v2);
seq=eq&:v1
@test randeval(seq.rhs)==randeval(((eq-:M*:v2)/:m).lhs)

eq=relations["Der"][3]
@test eval(parse(ps(eq)))==eq


text:
module Neurons
using Requires

include("base.jl")
include("to_TPN.jl")
include("to_TA.jl")

function __init__()
    # Optional dependency for makie-based plotting
    @require GLMakie="e9467ef8-e4e7-5192-8a1a-b1aee30e663a" begin
        include("plot_network.jl")
    end
end

end


text:
"""
Function heure_indiv

"""
function heure_indiv(statut::Int32)

    heure_depart = Journee.Semaine.Deb_Type ;
    heure_arrivee = Journee.Semaine.Fin_Type ;

    return(heure_depart,heure_arrivee)

end

"""
Function choix_resto

"""
function choix_resto(Prob::Array{Float64,1},id_aliment::Array{Int32,1})

    if(sum(Prob) > 0.5)
        id_lieu_loc = id_aliment[rand(Categorical(Prob))] ;
    else
        id_lieu_loc = id_aliment[floor(Int,rand(Uniform(1,length(id_aliment))))] ;
    end

    return(id_lieu_loc) ;

end


text:
# ========== Single Step Stationary Wavelet Transform ==========
# ----- 1D SWT -----
"""
    sdwt_step(v, d, h, g)

Perform one level of the stationary discrete wavelet transform (SDWT) on the vector `v`,
which is the `d`-th level scaling coefficients (Note the 0th level scaling coefficients is
the raw signal). The vectors `h` and `g` are the detail and scaling filters.

# Arguments
- `v::AbstractArray{T} where T<:Number`: Array of coefficients from a node at level `d`.
- `d::Integer`: Depth level of `v`.
- `h::Vector{S} where S<:Number`: High pass filter.
- `g::Vector{S} where S<:Number`: Low pass filter.

# Returns
- `w₁::AbstractVector{T} where T<:Number` or `w₁::AbstractMatrix{T} where T<:Number`: Output
  from low pass filter (1D case); or output from low + low pass filter (2D case).
- `w₂::AbstractVector{T} where T<:Number` or `w₂::AbstractMatrix{T} where T<:Number`: Output
  from high pass filter (1D case); or output from low + high pass filter (2D case).
- `w₃::AbstractVector{T} where T<:Number` or `w₃::AbstractMatrix{T} where T<:Number`: Output
  from high + low pass filter (2D case).
- `w₄::AbstractVector{T} where T<:Number` or `w₄::AbstractMatrix{T} where T<:Number`: Output
  from high + high pass filter (2D case).

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8)
wt = wavelet(WT.haar)
g, h = WT.makereverseqmfpair(wt, true)

# One step of SDWT
SWT.sdwt_step(v, 0, h, g)
```

**See also:** [`sdwt_step!`](@ref)
"""
function sdwt_step(v::AbstractVector{T}, d::Integer, h::Array{S,1}, g::Array{S,1}) where
                  {T<:Number, S<:Number}
    n = length(v)
    w₁ = zeros(T, n)
    w₂ = zeros(T, n)

    sdwt_step!(w₁, w₂, v, d, h, g)
    return w₁, w₂
end

"""
    sdwt_step!(w₁, w₂, v, d, h, g)

Same as `sdwt_step` but without array allocation.

# Arguments
- `w₁::AbstractVector{T} where T<:Number` or `w₁::AbstractMatrix{T} where T<:Number`: Vector
  allocation for output from low pass filter (1D case); or matrix allocation for output from
  low + low pass filter (2D case).
- `w₂::AbstractVector{T} where T<:Number` or `w₂::AbstractMatrix{T} where T<:Number`: Vector
  allocation for output from high pass filter (1D case); or matrix allocation for output
  from low + high pass filter (2D case).
- `w₃::AbstractVector{T} where T<:Number` or `w₃::AbstractMatrix{T} where T<:Number`: Matrix
  allocation for output from high + low pass filter (2D case).
- `w₄::AbstractVector{T} where T<:Number` or `w₄::AbstractMatrix{T} where T<:Number`: Matrix
  allocation for output from high + high pass filter (2D case).
- `v::AbstractArray{T} where T<:Number`: Array of coefficients from a node at level `d`.
- `d::Integer`: Depth level of `v`.
- `h::Vector{S} where S<:Number`: High pass filter.
- `g::Vector{S} where S<:Number`: Low pass filter.

# Returns
- `w₁::AbstractVector{T} where T<:Number` or `w₁::AbstractMatrix{T} where T<:Number`: Output
  from low pass filter (1D case); or output from low + low pass filter (2D case).
- `w₂::AbstractVector{T} where T<:Number` or `w₂::AbstractMatrix{T} where T<:Number`: Output
  from high pass filter (1D case); or output from low + high pass filter (2D case).
- `w₃::AbstractVector{T} where T<:Number` or `w₃::AbstractMatrix{T} where T<:Number`: Output
  from high + low pass filter (2D case).
- `w₄::AbstractVector{T} where T<:Number` or `w₄::AbstractMatrix{T} where T<:Number`: Output
  from high + high pass filter (2D case).

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8)
wt = wavelet(WT.haar)
g, h = WT.makereverseqmfpair(wt, true)
w₁ = zeros(8)
w₂ = zeros(8)

# One step of SDWT
SWT.sdwt_step!(w₁, w₂, v, 0, h, g)
```

**See also:** [`sdwt_step`](@ref)
"""
function sdwt_step!(w₁::AbstractVector{T},
                    w₂::AbstractVector{T},
                    v::AbstractVector{T},
                    d::Integer,
                    h::Array{S,1},
                    g::Array{S,1}) where {T<:Number, S<:Number}
    # Sanity check
    @assert length(w₁) == length(w₂) == length(v)
    @assert length(h) == length(g)

    # Setup
    n = length(v)                   # Signal length
    filtlen = length(h)             # Filter length

    # One step of stationary transform
    for i in 1:n
        k₁ = mod1(i-(1<<d),n)       # Start index for low pass filtering
        k₂ = i                      # Start index for high pass filtering
        @inbounds w₁[i] = g[end] * v[k₁]
        @inbounds w₂[i] = h[1] * v[k₂]
        for j in 2:filtlen
            k₁ = k₁+(1<<d) |> k₁ -> k₁>n ? mod1(k₁,n) : k₁
            k₂ = k₂-(1<<d) |> k₂ -> k₂≤0 ? mod1(k₂,n) : k₂
            @inbounds w₁[i] += g[end-j+1] * v[k₁]
            @inbounds w₂[i] += h[j] * v[k₂]
        end
    end
    return w₁, w₂
end

"""
    isdwt_step(w₁, w₂, d, h, g)
    isdwt_step(w₁, w₂, d, sv, sw, h, g)

Perform one level of the inverse stationary discrete wavelet transform (ISDWT) on the
vectors `w₁` and `w₂`, which are the `d+1`-th level scaling coefficients (Note the 0th level
scaling coefficients is the raw signal). The vectors `h` and `g` are the detail and scaling
filters.

!!! note
    One can decide to choose the average-based inverse transform or the shift-based inverse
    transform. For shift based, one needs to specify the shifts of the parent and children
    nodes; for the average based, the output is the average of all possible shift-based
    inverse transform outputs.

# Arguments
- `w₁::AbstractVector{T} where T<:Number` or `w₁::AbstractMatrix{T} where T<:Number`:
  Coefficients of left child node (1D case); or coefficients from top left child node (2D
  case).
- `w₂::AbstractVector{T} where T<:Number` or `w₂::AbstractMatrix{T} where T<:Number`:
  Coefficients of right child node (1D case); or coefficients from top right child node (2D
  case).
- `w₃::AbstractVector{T} where T<:Number` or `w₃::AbstractMatrix{T} where T<:Number`:
  Coefficients from bottom left child node (2D case).
- `w₄::AbstractVector{T} where T<:Number` or `w₄::AbstractMatrix{T} where T<:Number`:
  Coefficients from bottom right child node (2D case).
- `d::Integer`: Depth level of parent node of `w₁` and `w₂`.
- `sv::Integer`: Shift of parent node `v`.
- `sw::Integer`: Shift of children nodes `w₁` and `w₂`.
- `h::Vector{S} where S<:Number`: High pass filter.
- `g::Vector{S} where S<:Number`: Low pass filter.

# Returns
- `v::Vector{T}`: Reconstructed coefficients.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8)
wt = wavelet(WT.haar)
g, h = WT.makereverseqmfpair(wt, true)

# One step of SDWT
w₁, w₂ = SWT.sdwt_step(v, 0, h, g)

# One step of ISDWT
v̂ = SWT.isdwt_step(w₁, w₂, 0, h, g)         # Average based
ṽ = SWT.isdwt_step(w₁, w₂, 0, 0, 1, h, g)   # Shift based
```

**See also:** [`isdwt_step!`](@ref)
"""
function isdwt_step(w₁::AbstractVector{T},
                    w₂::AbstractVector{T},
                    d::Integer,
                    h::Array{S,1},
                    g::Array{S,1}) where {T<:Number, S<:Number}
    v = similar(w₁)
    isdwt_step!(v, w₁, w₂, d, h, g)
    return v
end

function isdwt_step(w₁::AbstractVector{T},
                    w₂::AbstractVector{T},
                    d::Integer,
                    sv::Integer,
                    sw::Integer,
                    h::Array{S,1},
                    g::Array{S,1}) where {T<:Number, S<:Number}
    v = similar(w₁)
    isdwt_step!(v, w₁, w₂, d, sv, sw, h, g)
    return v
end

"""
    isdwt_step!(v, w₁, w₂, d, h, g)
    isdwt_step!(v, w₁, w₂, d, sv, sw, h, g[; add2out])

Same as `isdwt_step` but without array allocation.

# Arguments
- `v::AbstractArray{T} where T<:Number`: Array allocation for reconstructed coefficients.
- `w₁::AbstractVector{T} where T<:Number` or `w₁::AbstractMatrix{T} where T<:Number`:
  Coefficients of left child node (1D case); or coefficients from top left child node (2D
  case).
- `w₂::AbstractVector{T} where T<:Number` or `w₂::AbstractMatrix{T} where T<:Number`:
  Coefficients of right child node (1D case); or coefficients from top right child node (2D
  case).
- `w₃::AbstractVector{T} where T<:Number` or `w₃::AbstractMatrix{T} where T<:Number`:
  Coefficients from bottom left child node (2D case).
- `w₄::AbstractVector{T} where T<:Number` or `w₄::AbstractMatrix{T} where T<:Number`:
  Coefficients from bottom right child node (2D case).
- `d::Integer`: Depth level of parent node of `w₁` and `w₂`.
- `sv::Integer`: Shift of parent node `v`.
- `sw::Integer`: Shift of children nodes `w₁` and `w₂`.
- `h::Vector{S} where S<:Number`: High pass filter.
- `g::Vector{S} where S<:Number`: Low pass filter.

# Keyword Arguments
- `add2out::Bool`: (Default: `false`) Whether to add computed result directly to output `v`
  or rewrite computed result to `v`.

# Returns
- `v::Array{T}`: Reconstructed coefficients.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8)
v̂ = similar(v)
ṽ = similar(v)
wt = wavelet(WT.haar)
g, h = WT.makereverseqmfpair(wt, true)

# One step of SDWT
w₁, w₂ = SWT.sdwt_step(v, 0, h, g)

# One step of ISDWT
SWT.isdwt_step!(v̂, w₁, w₂, 0, h, g)          # Average based
SWT.isdwt_step!(ṽ, w₁, w₂, 0, 0, 1, h, g)    # Shift based
```

**See also:** [`isdwt_step`](@ref)
"""
function isdwt_step!(v::AbstractVector{T},
                     w₁::AbstractVector{T},
                     w₂::AbstractVector{T},
                     d::Integer,
                     h::Array{S,1},
                     g::Array{S,1}) where {T<:Number, S<:Number}
    nd = 1 << d             # Number of blocks for parent node

    # isdwt_step for each shift
    for sv in 0:(nd-1)
        sw₁ = sv            # Shift of w₁, w₂ with no additional shift
        sw₂ = sv + 1<<d     # Shift of w₁, w₂ with addtional shift
        isdwt_step!(v, w₁, w₂, d, sv, sw₁, h, g)                 # Without shift
        isdwt_step!(v, w₁, w₂, d, sv, sw₂, h, g, add2out=true)   # With shift
    end

    # Get average output of shift vs no shift
    for i in eachindex(v)
        @inbounds v[i] /= 2
    end
end

function isdwt_step!(v::AbstractVector{T},
                     w₁::AbstractVector{T},
                     w₂::AbstractVector{T},
                     d::Integer,
                     sv::Integer,
                     sw::Integer,
                     h::Array{S,1},
                     g::Array{S,1};
                     add2out::Bool = false) where {T<:Number, S<:Number}
    # Sanity check
    @assert 0 ≤ sv < 1<<d
    @assert sv ≤ sw < 1<<(d+1)

    # Setup
    n = length(v)               # Signal length
    filtlen = length(h)         # filter length
    ip = sv + 1                 # Parent start index
    sp = 1 << d                 # Parent step size
    ic = sw + 1                 # Child start index
    sc = 1 << (d+1)             # Child step size

    # One step of inverse SDWT
    for (t, m) in enumerate(ip:sp:n)
        i₀ = mod1(t,2)                                  # Pivot point for filter
        i₁ = filtlen-i₀+1                               # Index for low pass filter g
        i₂ = mod1(t+1,2)                                # Index for high pass filter h
        j = sw==sv ? mod1(m-1<<d,n) : mod1(m+sp-1<<d,n) # Position of v, shift needed if sw≠sv
        k₁ = ((t-1)>>1) * sc + ic                       # Index for approx coefs w₁
        k₂ = ((t-1)>>1) * sc + ic                       # Index for detail coefs w₂
        @inbounds v[j] = add2out ? v[j]+g[i₁]*w₁[k₁]+h[i₂]*w₂[k₂] : g[i₁]*w₁[k₁]+h[i₂]*w₂[k₂]
        for i in (i₀+2):2:filtlen
            i₁ = filtlen-i+1
            i₂ = i + isodd(i) - iseven(i)
            k₁ = k₁-(1<<(d+1)) |> k₁ -> k₁≤0 ? mod1(k₁,n) : k₁
            k₂ = k₂+(1<<(d+1)) |> k₂ -> k₂>n ? mod1(k₂,n) : k₂
            @inbounds v[j] += g[i₁]*w₁[k₁]+h[i₂]*w₂[k₂]
        end
    end
    return v
end

# ----- 2D SWT -----
# Forward transform with allocation
function sdwt_step(v::AbstractMatrix{T}, d::Integer, h::Vector{S}, g::Vector{S}) where
                  {T<:Number, S<:Number}
    n, m = size(v)
    w₁ = Matrix{T}(undef, (n,m))
    w₂ = Matrix{T}(undef, (n,m))
    w₃ = Matrix{T}(undef, (n,m))
    w₄ = Matrix{T}(undef, (n,m))
    temp = Array{T,3}(undef, (n,m,2))
    sdwt_step!(w₁, w₂, w₃, w₄, v, d, h, g, temp)
    return w₁, w₂, w₃, w₄
end
# Forward transform without allocation
function sdwt_step!(w₁::AbstractMatrix{T}, w₂::AbstractMatrix{T},
                    w₃::AbstractMatrix{T}, w₄::AbstractMatrix{T},
                    v::AbstractMatrix{T},
                    d::Integer,
                    h::Vector{S}, g::Vector{S},
                    temp::AbstractArray{T,3}) where {T<:Number, S<:Number}
    # Sanity check
    @assert size(v) == size(w₁) == size(w₂) == size(w₃) == size(w₄)
    @assert ndims(temp) == 3
    @assert size(temp,1) == size(v,1)
    @assert size(temp,2) == size(v,2)

    # Setup
    n, m = size(w₁)

    # --- Transform ---
    # Compute sdwt for all columns
    for j in 1:m
        @inbounds temp₁ⱼ = @view temp[:,j,1]
        @inbounds temp₂ⱼ = @view temp[:,j,2]
        @inbounds vⱼ = @view v[:,j]
        @inbounds sdwt_step!(temp₁ⱼ, temp₂ⱼ, vⱼ, d, h, g)
    end
    # Compute sdwt for all rows
    for i in 1:n
        @inbounds temp₁ᵢ = @view temp[i,:,1]
        @inbounds w₁ᵢ = @view w₁[i,:]
        @inbounds w₂ᵢ = @view w₂[i,:]
        @inbounds sdwt_step!(w₁ᵢ, w₂ᵢ, temp₁ᵢ, d, h, g)

        @inbounds temp₂ᵢ = @view temp[i,:,2]
        @inbounds w₃ᵢ = @view w₃[i,:]
        @inbounds w₄ᵢ = @view w₄[i,:]
        @inbounds sdwt_step!(w₃ᵢ, w₄ᵢ, temp₂ᵢ, d, h, g)
    end
    return w₁, w₂, w₃, w₄
end

# Inverse transform (Average based) with allocation
function isdwt_step(w₁::AbstractMatrix{T}, w₂::AbstractMatrix{T},
                    w₃::AbstractMatrix{T}, w₄::AbstractMatrix{T},
                    d::Integer,
                    h::Array{S,1}, g::Array{S,1}) where {T<:Number, S<:Number}
    n, m = size(w₁)
    v = Matrix{T}(undef, (n,m))
    temp = Array{T,3}(undef, (n,m,2))
    isdwt_step!(v, w₁, w₂, w₃, w₄, d, h, g, temp)
    return v
end
# Inverse transform (Shift based) with allocation
function isdwt_step(w₁::AbstractMatrix{T}, w₂::AbstractMatrix{T},
                    w₃::AbstractMatrix{T}, w₄::AbstractMatrix{T},
                    d::Integer, sv::Integer, sw::Integer,
                    h::Vector{S}, g::Vector{S}) where {T<:Number, S<:Number}
    n, m = size(w₁)
    v = Matrix{T}(undef, (n,m))
    temp = Array{T,3}(undef, (n,m,2))
    isdwt_step!(v, w₁, w₂, w₃, w₄, d, sv, sw, h, g, temp)
    return v
end
# Inverse transform (Average based) without allocation
function isdwt_step!(v::AbstractMatrix{T},
                     w₁::AbstractMatrix{T}, w₂::AbstractMatrix{T},
                     w₃::AbstractMatrix{T}, w₄::AbstractMatrix{T},
                     d::Integer,
                     h::Vector{S}, g::Vector{S},
                     temp::AbstractArray{T,3}) where {T<:Number, S<:Number}
    # Sanity check
    @assert size(v) == size(w₁) == size(w₂) == size(w₃) == size(w₄)
    @assert ndims(temp) == 3
    @assert size(temp,1) == size(v,1)
    @assert size(temp,2) == size(v,2)

    # Setup
    n, m = size(w₁)

    # --- Transform ---
    # Compute isdwt for all rows
    for i in 1:n
        @inbounds temp₁ᵢ = @view temp[i,:,1]
        @inbounds w₁ᵢ = @view w₁[i,:]
        @inbounds w₂ᵢ = @view w₂[i,:]
        @inbounds isdwt_step!(temp₁ᵢ, w₁ᵢ, w₂ᵢ, d, h, g)

        @inbounds temp₂ᵢ = @view temp[i,:,2]
        @inbounds w₃ᵢ = @view w₃[i,:]
        @inbounds w₄ᵢ = @view w₄[i,:]
        @inbounds isdwt_step!(temp₂ᵢ, w₃ᵢ, w₄ᵢ, d, h, g)
    end
    # Compute isdwt for all columns
    for j in 1:m
        @inbounds temp₁ⱼ = @view temp[:,j,1]
        @inbounds temp₂ⱼ = @view temp[:,j,2]
        @inbounds vⱼ = @view v[:,j]
        @inbounds isdwt_step!(vⱼ, temp₁ⱼ, temp₂ⱼ, d, h, g)
    end
    return v
end
# Inverse transform (Shift based) without array allocation
function isdwt_step!(v::AbstractMatrix{T},
                     w₁::AbstractMatrix{T}, w₂::AbstractMatrix{T},
                     w₃::AbstractMatrix{T}, w₄::AbstractMatrix{T},
                     d::Integer, sv::Integer, sw::Integer,
                     h::Vector{S}, g::Vector{S},
                     temp::AbstractArray{T,3}) where {T<:Number, S<:Number}
    # Sanity check
    @assert size(v) == size(w₁) == size(w₂) == size(w₃) == size(w₄)
    @assert ndims(temp) == 3
    @assert size(temp,1) == size(v,1)
    @assert size(temp,2) == size(v,2)

    # Setup
    n, m = size(w₁)

    # --- Transform ---
    # Compute isdwt for all rows
    for i in 1:n
        @inbounds temp₁ᵢ = @view temp[i,:,1]
        @inbounds w₁ᵢ = @view w₁[i,:]
        @inbounds w₂ᵢ = @view w₂[i,:]
        @inbounds isdwt_step!(temp₁ᵢ, w₁ᵢ, w₂ᵢ, d, sv, sw, h, g)

        @inbounds temp₂ᵢ = @view temp[i,:,2]
        @inbounds w₃ᵢ = @view w₃[i,:]
        @inbounds w₄ᵢ = @view w₄[i,:]
        @inbounds isdwt_step!(temp₂ᵢ, w₃ᵢ, w₄ᵢ, d, sv, sw, h, g)
    end
    # Compute isdwt for all columns
    for j in 1:m
        @inbounds temp₁ⱼ = @view temp[:,j,1]
        @inbounds temp₂ⱼ = @view temp[:,j,2]
        @inbounds vⱼ = @view v[:,j]
        @inbounds isdwt_step!(vⱼ, temp₁ⱼ, temp₂ⱼ, d, sv, sw, h, g)
    end
    return v
end


text:
using RowTables
using DataFrames
using BenchmarkTools

## DataFrames will not run on Jupyter at the moment

ds = [
    Dict(:a=>2,:b=>1,:c=>7)
    Dict(:a=>1,:b=>4,:c=>8)
    Dict(:a=>9,:b=>6,:c=>4)
    Dict(:a=>10,:b=>8,:c=>7)
    Dict(:a=>10,:b=>4,:c=>3)];

rt = RowTable(ds,[:a,:b,:c])

df = DataFrame(rt)

df[3,:]

rt[3,:]

rt[3:3,:]

@btime rt[3,:]

@btime df[3,:]

@btime rt[3:3,:]

@btime rt[3:4,:]

@btime df[3:4,:]

@btime df[3,:]

df[:,1]

rt[:,1]

@btime rt[:,1];

@btime df[:,1];

@btime DataFrame(rt)

@btime RowTable(df)


text:
using Sundials
pyplot(size=(700,700))

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  ## Define the system residual function.
  function resrob(tres, y, yp, r)
      r[1]  = -0.04*y[1] + 1.0e4*y[2]*y[3]
      r[2]  = -r[1] - 3.0e7*y[2]*y[2] - yp[2]
      r[1] -=  yp[1]
      r[3]  =  y[1] + y[2] + y[3] - 1.0
  end

  #t = [0.0, 4.0 * logspace(-1., 5., 7)]
  t = [0.0, 0.4, 4.0, 40.0, 400.0, 4000.0, 40000.0, 400000.0]
  yout, ypout = Sundials.idasol(resrob, [1.0, 0, 0], [-0.04, 0.04, 0.0], t)
  yout |> display
  println()
  ypout |> display
  println()

  p1 = plot(t, yout[:,1], color=:darkblue,
    xlabel="x", ylabel="y[1]",
    title="Example Sundials01A (using Sundials.jl)"
  )
  plot!(t, yout[:,2], color=:red)
  plot!(t, yout[:,3], color=:green)

  p2 = plot(t, ypout[:,1], color=:darkred,
    xlabel="x", ylabel="y[2]"
  )
  plot!(t, ypout[:,2], color=:red)
  plot!(t, ypout[:,3], color=:green)

  plot(p1, p2, layout=(1, 2))
  savefig("Sundials01B.png")
  gui()

end


text:
# Sample network and compute observables for a set of operators stored as
# transposed sparse matrices on |ψ⟩ Neural Quantum States (KetNet).
function sample_network!(
  res::MCMCObsEvaluationCache,
  problem::ObservablesProblem,
  net::KetNet, σ, wholespace=false)
  # Load tmp caches to avoid allocating
  σp = res.σ

  # The denominator of this state
  lnψ = net(σ)
  prob = wholespace ? exp(2*real(lnψ)) : 1.0
  res.Zave += prob

  # Iterate through all elements in row i_σ of the matrix computing
  # Since Julia has no CSR matrices, I iterate on columns of the
  # transpose matrix
  i_σ = index(σ)
  for (obs_id, Oᵗ) = enumerate(problem.ObservablesTransposed)
    O_loc = 0.0+0.0im
    for row_id = Oᵗ.colptr[i_σ]:(Oᵗ.colptr[i_σ+1]-1)
      # Find nonzero elements s by doing <i_sp|Oᵗ|i_σ>
      i_σp = Oᵗ.rowval[row_id]
      set_index!(σp, i_σp)

      # Compute the log(ψ(σ)/ψ(σ')), by only computing differences.
      log_ratio = net(σp) - lnψ
      O_loc += conj(Oᵗ.nzval[row_id]) * exp(log_ratio) # conj because transpose
    end
    res.ObsAve[obs_id] += prob * O_loc
    push!(res.ObsVals[obs_id], prob * O_loc)
  end
  return res
end

# Sample network and compute observables for a set of operators stored as a
# KLocalOperator on |ψ⟩ Neural Quantum States (KetNet).
function sample_network!(
  res::MCMCObsEvaluationCache,
  problem::ObservablesProblem{B,SM},
  net::KetNet, σ, wholespace=false) where {B, SM<:AbsLinearOperator}
  # Load tmp caches to avoid allocating
  σp = res.σ

  # The denominator of this state
  lnψ = net(σ)
  prob = wholespace ? exp(2*real(lnψ)) : 1.0
  res.Zave += prob

  i_σ = index(σ)
  for (obs_id, O) = enumerate(problem.ObservablesTransposed)
    O_loc = 0.0+0.0im
    #diffs_O = row_valdiff(O, σ)
    for op=operators(O)
      r = local_index(σ, sites(op))
      for (mel, changes)=op.op_conns[r]
        set_index!(σp, i_σ)
        for (site, val)=changes
          setat!(σp, site, val)
        end

        # Compute the log(ψ(σ)/ψ(σ')), by only computing differences.
        log_ratio = net(σp) - lnψ
        O_loc += mel * exp(log_ratio)
      end
    end

    res.ObsAve[obs_id] += prob * O_loc
    push!(res.ObsVals[obs_id], prob * O_loc)
  end
  return res
end


text:
using Parameters, StatsBase, Random
import Distributions: logpdf, rand, loglikelihood

struct Siegler{T1,T2,T3,T4} <: ContinuousUnivariateDistribution
    δ::T1
    τ::T2
    s::T3
    parms::T4
end

Siegler(;δ, τ, s, parms) = Siegler(δ, τ, s, parms)

loglikelihood(d::Siegler, data::Array{<:NamedTuple,1}) = logpdf(d, data)

function logpdf(d::Siegler, data::Array{<:NamedTuple,1})
    LL = computeLL(d.parms, data; δ=d.δ, τ=d.τ, s=d.s)
    return LL
end

function simulate(stimuli, parms; args...)
    # populate chunks
    chunks = populate_memory()
    # set the base level constant based on sum slot
    set_baselevels!(chunks)
    # create a declarative memory object
    memory = Declarative(;memory=chunks)
    # create an ACTR model object
    actr = ACTR(;declarative=memory, parms..., args...)
    shuffle!(stimuli)
    N = length(stimuli)
    data = Array{NamedTuple,1}(undef, N)
    # loop through trials
    for (i,s) in enumerate(stimuli)
        # retrieve a chunk based on retrieval request s
        chunk = retrieve(actr; s...)
        # compute reaction time
        rt = compute_RT(actr, chunk) + parms.ter
        if isempty(chunk)
            # retrieval failure
            data[i] = (s...,resp = -100,rt = rt)
        else
            # chunk retrieved.
            data[i] = (s...,resp = chunk[1].slots.sum,rt = rt)
        end
    end
    return data
end

function sim_fun(actr, chunk; request...)
    slots = chunk.slots
    p = 0.0; δ = actr.parms.δ
    for (c,v) in request
        p += .1 * δ * abs(slots[c] - v)
    end
    return p
end

function populate_memory(act=0.0)
    chunks = [Chunk(;num1=num1,num2=num2,
        sum=num1 + num2,act=act) for num1 in 0:5
        for num2 in 0:5]
    pop!(chunks)
    return chunks
end

function set_baselevels!(chunks)
    for chunk in chunks
        if chunk.slots.sum < 5
            chunk.bl = 0.65
        end
    end
    return nothing
end

function computeLL(parms, data; δ, τ, s)
    type = typeof(δ)
    # populate chunks
    chunks = populate_memory(zero(δ))
    # set base level constant based on sum slot
    set_baselevels!(chunks)
    # create declarative memory object
    memory = Declarative(;memory=chunks)
    # create ACTR model object
    actr = ACTR(;declarative=memory, parms..., δ, τ, s)
    # remove random values from activation
    actr.parms.noise = false
    N = length(chunks) + 1
    @unpack s,ter,τ = actr.parms
    LL = 0.0; idx = 0; μ = Array{type,1}(undef, N)
    σ = s * pi / sqrt(3)
    ϕ = ter
    for (num1,num2,resp,rt) in data
        # compute activation
        compute_activation!(actr; num1, num2)
        # extract mean activation values
        map!(x -> x.act, μ, chunks)
        # last mean activation is for retrieval failure
        μ[end] = τ
        # log normal race distribution
        dist = LNR(;μ=-μ, σ, ϕ)
        # no retrieval error
        if resp != -100
            # get all chunk indices such that sum = response
            indices = find_indices(actr; sum=resp)
            log_probs = zeros(type, length(indices))
            # loop over each component of mixture and compute log likelihood
            for (c,idx) in enumerate(indices)
                log_probs[c] = logpdf(dist, idx, rt)
            end
            # compute marginal likelihood
            LL += logsumexp(log_probs)
        else
            # retrieval failure
            LL += logpdf(dist, N, rt)
        end
    end
    return LL
end

function rt_histogram(df, stimuli; kwargs...)
    vals = NamedTuple[]
    for (num1,num2) in stimuli
        idx = @. ((df[:,:num1] == num1 ) & (df[:,:num2] == num2)) | ((df[:,:num1] == num2) & (df[:,:num2] == num1))
        subdf = df[idx,:]
        str = string(num1, "+", num2)
        temp = filter(x -> x[:resp] != -100, subdf)
        g = groupby(temp, :resp)
        rt_resp = combine(g, :rt => mean)
        push!(vals, (title = str, data = rt_resp))
    end
    p = bar(layout=(2,3), leg=false, xlims=(0,10), xlabel="Response", ylabel="Mean RT",
        size=(600,600), xaxis=font(6), yaxis=font(6), titlefont=font(7), grid=false; kwargs...)
    for (i,v) in enumerate(vals)
        @df v.data bar!(p, :resp, :rt_mean, subplot=i, title=v.title, bar_width=1, color=:grey,
            grid=false, ylims=(0,4.5), xlims=(-.5,9))
    end
    return p
end

function response_histogram(df, stimuli; kwargs...)
    vals = NamedTuple[]
    for (num1,num2) in stimuli
        idx = @. ((df[:,:num1] == num1 ) & (df[:,:num2] == num2)) | ((df[:,:num1] == num2) & (df[:,:num2] == num1))
        subdf = df[idx,:]
        str = string(num1, "+", num2)
        v = filter(x -> x != -100, subdf[:,:resp])
        push!(vals, (title = str, data = v))
    end
    p = histogram(layout=(2,3), leg=false, xlims=(0,10), xlabel="Response",ylabel="Proportion",
        size=(600,600), xaxis=font(6), yaxis=font(6), titlefont=font(7), grid=false; kwargs...)
    for (i,v) in enumerate(vals)
        histogram!(p, v.data, subplot=i, title=v.title, bar_width=1, color=:grey, grid=false,
        normalize=:probability, ylims=(0,1))
    end
    return p
end


text:
module SeisSolXDMF

export XDMFFile, timesteps_of, grid_of, data_of

using XMLDict
using HDF5

    struct XDMFFile
        xml         :: XMLDict.XMLDictElement
        base_path   :: AbstractString
        timesteps   :: Array{XMLDict.XMLDictElement,1}
    end

    function expect(expression, error_msg)
        if !expression
            error(error_msg)
        end
    end

    function XDMFFile(file_path::AbstractString)
        xml_file = open(file_path) do file; read(file, String) end
        xml_file = parse_xml(xml_file)

        try
            # The root node "Xdmf" is assumed
            timesteps = xml_file["Domain"]["Grid"]["Grid"]
            num_timesteps = length(timesteps)

            expect(num_timesteps >  0, "No timesteps found in XDMF file.")

            xdmf_base_path = splitdir(file_path)[1]
            return XDMFFile(xml_file, xdmf_base_path, timesteps)
        catch
            error("Input file has no valid SeisSol XDMF structure.")
        end
    end

    function timesteps_of(xdmf::XDMFFile) :: Array{Float64, 1}
        return Array(map(data_item -> parse(Float64, data_item["Time"][:Value]), xdmf.timesteps))
    end

    function grid_of(xdmf::XDMFFile) :: Tuple{AbstractArray{Integer, 2}, AbstractArray{AbstractFloat, 2}}
        # Read geometry (points) and topology (simplices (triangles/tetrahedra)) files.
        # Since the ids of the points referred to in the topology array start at 0 (and Julia counts from 1) we have to add 1.
        simplices = read_dataset(xdmf, xdmf.timesteps[1]["Topology"]["DataItem"]) .+ 1
        expect(ndims(simplices) == 2, "Topology should have 2 dimensions but has $(ndims(simplices)).")
        expect(size(simplices, 1) ∈ (3, 4), "Topology entries must have 3 (triangles) or 4 (tetrahedra) points but have $(size(simplices, 1)).")

        points = read_dataset(xdmf, xdmf.timesteps[1]["Geometry"]["DataItem"])
        expect(ndims(points) == 2, "Geometry should have 2 dimensions but has $(ndims(simplices)).")
        expect(size(points, 1) == 3, "Geometry entries must have 3 coordinates but have $(size(simplices, 1)).")

        return (simplices, points)
    end

    function data_of(xdmf::XDMFFile, timestep::Integer, var_name::AbstractString)
        expect(timestep >= 1, "Timestep $(timestep) is smaller than 1.")
        expect(timestep <= length(xdmf.timesteps), "Timestep $(timestep) exceeds number of timesteps in XDMF file ($(length(xdmf.timesteps))).")

        attr_or_list = xdmf.timesteps[timestep]["Attribute"]
        var_entry = if isa(attr_or_list, AbstractArray)
            var_index = findfirst(attr -> attr[:Name] == var_name, attr_or_list)
            expect(!isnothing(var_index), """Variable $var_name not found in timestep $timestep.
                Found variables: $(join(map(attr -> attr[:Name], attr_or_list), ", ")).""")

            attr_or_list[var_index]["DataItem"]
        else
            attr = attr_or_list
            expect(attr[:Name] == var_name, """Variable $var_name not found in timestep $timestep. Found variable: $(attr[:Name]).""")

            attr["DataItem"]
        end
        return reshape(read_hyperslab(xdmf, var_entry), :)
    end

    function read_dataset(xdmf::XDMFFile, data_item::XMLDict.XMLDictElement; indices=nothing)
        filename = data_item[""]

        file_range = data_item[:Dimensions]
        file_range = split(file_range, ' ')
        file_range = parse.(Int, file_range)
        file_period = file_range[2]

        number_type = get_number_type(data_item)

        if data_item[:Format] == "HDF"
            path_parts = split(filename, ':', limit=2)
            expect(!isempty(path_parts[1]) && !isempty(path_parts[2]), "HDF5 group path is invalid.")

            filename = joinpath(xdmf.base_path, path_parts[1])
            hdf_dataset_path = String(path_parts[2])

            dataset =
                if isnothing(indices)
                    h5read(filename, hdf_dataset_path, (1:file_range[2], 1:file_range[1]))
                else
                    h5read(filename, hdf_dataset_path, (
                        indices[1,2]:indices[2,2]:(indices[1,2] + indices[3,2] - 1),
                        indices[1,1]:indices[2,1]:(indices[1,1] + indices[3,1] - 1)))
                end
        else
            filename = joinpath(xdmf.base_path, filename)

            if !isnothing(indices)
                start_index = (indices[1,1]-1) * file_period + (indices[1,2]-1) # Starting at 0
                stride = (indices[2,1]-1) * file_period + (indices[2,2]-1)
                expect(stride == 0, "Stride not yet supported for POSIX output.") # (no need (?))
                read_size = (indices[3,2], indices[3,1])
            else
                start_index = 0
                read_size = (file_range[2], file_range[1])
            end
            dataset = Array{number_type, 2}(undef, read_size)
            open(filename, "r") do file
                seek(file, start_index * sizeof(number_type))
                read!(file, dataset)
            end
        end

        return dataset
    end

    function read_hyperslab(xdmf::XDMFFile, data_item::XMLDict.XMLDictElement)
        #=
        Example DataItem:
        <DataItem ItemType="HyperSlab" Dimensions="4561037">
            <DataItem NumberType="UInt" Precision="4" Format="XML" Dimensions="3 2">0 0 1 1 1 4561037</DataItem>
            <DataItem NumberType="Float" Precision="8" Format="Binary" Dimensions="1 5242880">out_cell/mesh0/u.bin</DataItem>
        </DataItem>
        =#
        expect(data_item[:ItemType] == "HyperSlab", "Unexpected item type '$(data_item[:ItemType])', expected 'HyperSlab'.")
        range_item = data_item["DataItem"][1]
        binary_item = data_item["DataItem"][2]

        expect(range_item[:Format] == "XML", "Unexpected format '$(range_item[:Format])', expected 'XML'.")
        expect(range_item[:Dimensions] == "3 2", "Unexpected dimensions '$(range_item[:Dimensions])', expected '3 2'.")
        expect(binary_item[:Format] ∈ ["Binary", "HDF"], "Unexpected format '$(binary_item[:Format])', expected 'Binary' or 'HDF'.")

        rng = range_item[""]
        rng = split(rng, ' ')
        rng = parse.(Int, rng)
        rng = reshape(rng, (2, 3))'

        rng[1,:] .+= 1 # Julia start indices

        return read_dataset(xdmf, binary_item, indices=rng)
    end

    function get_number_type(data_item::XMLDict.XMLDictElement) :: Type
        bytes_of_precision = parse(UInt, data_item[:Precision])
        number_type_name = data_item[:NumberType]

        print_error(type_name, precision) =
            error("""Type "$type_name" of size $precision bytes is not recognized by the parser.""")

        if bytes_of_precision == 8
            if number_type_name == "Int"
                return Int64
            elseif number_type_name == "UInt"
                return UInt64
            elseif number_type_name == "Float"
                return Float64
            else
                print_error(number_type_name, bytes_of_precision)
            end
        elseif bytes_of_precision == 4
            if number_type_name == "Int"
                return Int32
            elseif number_type_name == "UInt"
                return UInt32
            elseif number_type_name == "Float"
                return Float32
            else
                print_error(number_type_name, bytes_of_precision)
            end
        else
            print_error(number_type_name, bytes_of_precision)
        end
    end

end # module

max_stars_repo_name: mtronicsDev/SeisSolXDMF

max_issues_repo_name: mtronicsDev/SeisSolXDMF

max_forks_repo_name: mtronicsDev/SeisSolXDMF


text:
using Random
using MLJModelInterface
using Test

mutable struct Foo <: MLJType
    rng::AbstractRNG
    x::Int
    y::Int
end

mutable struct Bar{names} <: MLJType
    rng::AbstractRNG
    v::Tuple{Int,Int}
    Bar{names}(rng, x, y) where names = new{names}(rng, (x, y))
end

Bar(rng, x, y) = Bar{(:x, :y)}(rng, x, y)

# overload `getproperty` so that components of `v` are accessed with
# the names given in `names` (which will be (:x, :y) when using
# the above constructor):
Base.propertynames(::Bar{names}) where names = (:rng, names...)

function Base.getproperty(b::Bar{names}, name::Symbol) where names
    name === :rng && return getfield(b, :rng)
    v = getfield(b, :v)
    name === names[1] && return v[1]
    return v[2]
end
function Base.setproperty!(b::Bar{names}, name::Symbol, value) where names
    name === :rng && return setfield!(b, :rng, value)
    v = getfield(b, :v)
    if name === names[1]
        setfield!(b, :v, (value, v[2]))
        return value
    end
    setfield!(b, :v, (v[1], value))
    return value
end

mutable struct Super <: MLJType
    sub::Foo
    z::Int
end

mutable struct Partial <: MLJType
    x::Int
    y::Vector{Int}
    Partial(x) = new(x)
end

mutable struct Sub
    x::Int
end

mutable struct Deep
    x::Int
    s::Union{Sub,Int}
end

mutable struct Super2 <: MLJType
    sub::Sub
    z::Int
end

MLJModelInterface.deep_properties(::Type{<:Super2}) = (:sub,)

@testset "_isdefined" begin
    b = Bar(MersenneTwister(), 2, 3)
    @test MLJModelInterface._isdefined(b, :x)
end

@testset "_equal_to_depth_one" begin
    d1 = Deep(1, 2)
    d2 = Deep(1, 2)
    @test MLJModelInterface._equal_to_depth_one(d1, d2)
    d2.x = 3
    @test !MLJModelInterface._equal_to_depth_one(d1, d2)

    d1 = Deep(1, Sub(2))
    d2 = Deep(1, Sub(2))
    @test !MLJModelInterface._equal_to_depth_one(d1, d2)
end

@testset "equality for MLJType" begin
    f1 = Foo(MersenneTwister(7), 1, 2)
    f2 = Foo(MersenneTwister(8), 1, 2)

    @test f1.rng != f2.rng
    @test f1 == f2
    f1.x = 10
    @test f1 != f2
    b = Bar(MersenneTwister(7), 1, 2)
    @test f2 != b

    @test is_same_except(f1, f2, :x)
    f1.y = 20
    @test f1 != f2
    @test is_same_except(f1, f2, :x, :y)

    f1 = Foo(MersenneTwister(7), 1, 2)
    f2 = Foo(MersenneTwister(8), 1, 2)
    s1 = Super(f1, 20)
    s2 = Super(f2, 20)
    @test s1 == s2
    s2.sub.x = 10
    @test f1 != f2

    @test !(f1 == Super(f1, 4))

    @test !isequal(
        Foo(MersenneTwister(1), 1, 2),
        Foo(MersenneTwister(1), 1, 2)
    )

    p1 = Partial(1)
    p2 = Partial(1)
    @test p1 == p2
    p2.y = [1,2]
    @test !(p1 == p2)

    # test of "deep" properties
    s1 = Super2(Sub(1), 2)
    s2 = Super2(Sub(1), 2)
    @test s1 == s2

    # replacing value of a "strap-on" property is detected:
    b = Bar(MersenneTwister(), 1, 2)
    b1 = deepcopy(b)
    b.x = 42
    @test b != b1
end

@testset "in(x, collection) for MLJType" begin
    f1 = Foo(MersenneTwister(7), 1, 2)
    f2 = Foo(MersenneTwister(7), 1, 2)
    f3 = Super(f1, 20)

    tv = (f1, f3)
    tk = (f2, f3)
    tw = (f3, f3)
    v = [tv...]
    k = [tk...]
    w = [tw...]
    @test f1 in tv
    @test !(f1 in tk)
    @test !(f1 in tw)
    @test f1 in v
    @test !(f1 in k)
    @test !(f1 in w)
    @test f1 in Set(v)
    @test !(f1 in Set(k))
    @test !(f1 in Set(w))

end

@testset "isrepresented" begin
    m = Foo(MersenneTwister(7), 1, 2)
    m2 = Foo(MersenneTwister(8), 1, 2)
    n = Foo(MersenneTwister(7), 3, 2)
    p = Foo(MersenneTwister(9), 4, 5)
    models = [m, n]

    @test isrepresented(m, nothing) == false
    @test isrepresented(m, models)
    @test isrepresented(m2, models)
    @test !isrepresented(p, models)
end

true


text:
function initIT_range(N::Int, components::Vector{Int}, f::Vector{Float64} = Float64[]; MAXITER::Int=100)
	#components = Int[]
	nelements = zeros(Int, TOTALCOMPONENTS)
	elements = Vector{Vector{Int}}(TOTALCOMPONENTS)
	obs = N

	if STEP in components
		#push!(components,STEP)
		nelements[STEP] = N-1
		elements[STEP] = collect( 1:(N-1) )
	else
		nelements[STEP] = 0
		elements[STEP] = Int[]
	end

	if SPIKE in components
		#push!(components,SPIKE)
		if ( STEP in components || SLOPE in components ) # last step and last spike are the same
			nelements[SPIKE] = N-1
			elements[SPIKE] = collect( 1:(N-1) )
		else
			nelements[SPIKE] = N
			elements[SPIKE] = collect( 1:N )
		end
	else
		nelements[SPIKE] = 0
		elements[SPIKE] = Int[]
	end

	if SLOPE in components
		#push!(components,SLOPE)
		nelements[SLOPE] = N-1
		elements[SLOPE] = collect( 1:(N-1) )
	else
		nelements[SLOPE] = 0
		elements[SLOPE] = Int[]
	end

	if SIN in components
		#push!(components,SIN)
		nelements[SIN] = length(f)
		elements[SIN] = collect( 1:length(f))
	else
		nelements[SIN] = 0
		elements[SIN] = Int[]
	end

	if COS in components
		#push!(components,COS)
		nelements[COS] = length(f)
		elements[COS] = collect( 1:length(f))
	else
		nelements[COS] =  0
		elements[COS] = Int[]
	end

	return iterator(obs, components, elements, nelements, ones(Bool,TOTALCOMPONENTS), sum(nelements),MAXITER)
end

function initSparse(IT::iterator)

	BCD = Array{SparseMatrixCSC{Float64,Int},1}[]

	beta_tilde = SparseMatrixCSC{Float64,Int}[]
	beta = SparseMatrixCSC{Float64,Int}[]

	activeSet = SparseMatrixCSC{Bool,Int}[]

	for i in 1:TOTALCOMPONENTS

		if i in IT.components

			push!(beta_tilde,spzeros(IT.nelements[i],1))
			push!(beta,spzeros(IT.nelements[i],1))
			push!(activeSet,spzeros(Bool,IT.nelements[i],1))

		else

			push!(beta_tilde,spzeros(0,0))
			push!(beta,spzeros(0,0))
			push!(activeSet,spzeros(Bool,0,0))

		end

	end

	return BCD, beta_tilde, beta, activeSet
end

function initDense(IT)
	BCD = Array{Vector{Float64},1}[]

	beta_tilde = Vector{Float64}[]
	beta = Vector{Float64}[]

	activeSet = Vector{Bool}[]

	for i in 1:TOTALCOMPONENTS

		if i in IT.components

			push!(beta_tilde, zeros(IT.nelements[i]))
			push!(beta, zeros(IT.nelements[i]))
			push!(activeSet,zeros(Bool,IT.nelements[i]))

		else

			push!(beta_tilde,zeros(0))
			push!(beta, zeros(0))
			push!(activeSet,zeros(Bool,0))

		end

	end

	return BCD, beta_tilde, beta, activeSet
end

function initXDY(IT,y,data)

	xdy0 = Vector{Float64}[]

	for i in 1:TOTALCOMPONENTS
		if in(i,IT.components)

			push!(xdy0, xdy[i](IT,y,data) )

		else
			push!(xdy0, Float64[] )
		end
		#temp)
	end

	return xdy0
end

#change for no preallocation and simply edit field
#requires changing normalizations.jl (returns)
function initData(IT::iterator, fs::Vector{Float64}=Float64[], fc::Vector{Float64}=Float64[])

	d = dataCD()
	for i in IT.components
		d.σ[i],d.μ[i] = getData[i](IT,fs)
	end
	d.fs = fs
	d.fc = fc

	return d
end

max_stars_repo_name: joaquimg/ComponentFilterCD

max_issues_repo_name: joaquimg/ComponentFilterCD

max_forks_repo_name: joaquimg/ComponentFilterCD


text:
#=
   From Allen Downey
   https://twitter.com/AllenDowney/status/1063460117029535746
   """
   Today's Bayesian problem of the week: Suppose we visit a wild animal preserve where we
   know that the only animals are lions and tigers and bears, but we don't know how
   many of each there are.

   During the tour, we see 3 lions, 2 tigers, and 1 bear. Assuming that every animal had an equal
   chance to appear in our sample, estimate the prevalence of each species.

   What is the probability that the next animal we see is a bear?
   """

   Also see: https://towardsdatascience.com/estimating-probabilities-with-bayesian-modeling-in-python-7144be007815
    From https://towardsdatascience.com/estimating-probabilities-with-bayesian-modeling-in-python-7144be007815
    """
    Species: lions    Prevalence: 44.44%.
    Species: tigers   Prevalence: 33.33%.
    Species: bears    Prevalence: 22.22%.
    """

   Distributions of variable animal7
   lion       =>    4526  (0.452600)
   tiger      =>    3251  (0.325100)
   bear       =>    2223  (0.222300)

   Cf ~/blog/lions_tigers_and_bears.blog
      ~/webppl/lions_tigers_and_bears.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function lions_tigers_and_bears()
    lion = 1
    tiger = 2
    bear = 3

    # Prior
    # probLion  ~ Beta(1,1)
    # probTiger ~ Beta(1,1)
    # probBear  ~ Beta(1,1)

    # The Dirichlet distribution ensures that the sum of probabilities is 1
    # i.e. we don't have to ensure this via some specific constraint.
    v = [3, 2, 1] # We saw three lions, two tigers, and one bear
    # v = [2,2,2]
    x ~ Dirichlet(v)

    # The probabilities to calculate ("aliased" for simplicity)
    probLion  ~ Dirac(x[1])
    probTiger ~ Dirac(x[2])
    probBear  ~ Dirac(x[3])

    n = 7
    o = tzeros(n)
    for i in 1:n
        o[i] ~ Categorical([probLion,probTiger,probBear])
    end

    # It shouldn't matter in what order we see the different animals.
    true ~ Dirac(o[1] == lion)
    true ~ Dirac(o[2] == lion)
    true ~ Dirac(o[3] == lion)
    true ~ Dirac(o[4] == tiger)
    true ~ Dirac(o[5] == tiger)
    true ~ Dirac(o[6] == bear)

    animal7 ~ Dirac(o[7])

    # probbear7 ~ Dirac(o[7] == bear)
    # problion7 ~ Dirac(o[7] == lion)
    # probtiger7 ~ Dirac(o[7] == tiger)

end

model = lions_tigers_and_bears()
num_chns = 4

# chns = sample(model, Prior(), 10_000)

# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(15), 10_000)
# chns = sample(model, SMC(1000), 10_000)
# chns = sample(model, IS(), 10_000)

#
display(chns)

show_var_dist_pct(chns,:animal7,["lion","tiger","bear"])

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
export SymSemiseparableChol

struct SymSemiseparableChol <: SymSemiseparableCholesky
    n::Int64
    p::Int64
    U::AbstractArray
    W::AbstractArray
end

# Constuctors
function SymSemiseparableChol(U::AbstractArray, W::AbstractArray)
	if size(U,1) == size(W,1) && size(U,2) == size(W,2)
		return SymSemiseparableChol(size(U,1),size(U,2),U,W)
	else
		error("Dimension mismatch between the generators U and W")
	end
end

# Mappings
mul!(y::AbstractArray, L::SymSemiseparableChol, 		         b::AbstractArray) =
	ss_tri_mul!(y, L.U,   L.W,   b)
mul!(y::AbstractArray, L::AdjointOperator{SymSemiseparableChol}, b::AbstractArray) =
	ssa_tri_mul!(y, L.A.U, L.A.W, b)
inv!(y::AbstractArray, L::SymSemiseparableChol, 		    	 b::AbstractArray) =
	ss_forward!(y, L.U,   L.W,   b);
inv!(y::AbstractArray, L::AdjointOperator{SymSemiseparableChol}, b::AbstractArray) =
	ssa_backward!(y, L.A.U, L.A.W, b);
newlogdet(L::SymSemiseparableChol) = ss_logdet(L.U, L.W)
newlogdet(L::AdjointOperator{SymSemiseparableChol}) = ss_logdet(L.A.U, L.A.W)

#### Inverse of a SymSemiseparableChol using ####
function inv(L::SymSemiseparableChol)
	return L'\(L\Diagonal(ones(L.n)))
end
function inv(L::SymSemiseparableChol, b::AbstractArray)
	return L'\(L\b)
end

########################################################################
#### Cholesky factoriaztion of:                                     ####
####    Extended generator representable {p}-semiseperable matrices ####
########################################################################

# Creating W, s.t. L = tril(UW')
function ss_create_w(U,V)
    n,p = size(U);
    wTw = zeros(p,p);
    W = zeros(n,p)
    for j = 1:n
        tmpu = U[j,:];
        tmp  = V[j,:] - wTw*tmpu;
        w = tmp/sqrt(abs(tmpu'*tmp))
        W[j,:] = w;
        wTw = wTw + w*w';
    end
    return W
end
#### Triangular product ####
function ss_tri_mul!(Y::AbstractArray,U::AbstractArray,
                     W::AbstractArray,X::AbstractArray)
     n, m = size(U);
     mx = size(X,2);
     Wbar = zeros(m,mx);
     for i = 1:n
         tmpW   = W[i,:]
         tmpU   = U[i,:]
         tmpX   = X[i:i,:]
         Wbar  += tmpW .* tmpX;
         Y[i,:] = Wbar'*tmpU;
     end
end
#### Adjoint triangular product ####
function ssa_tri_mul!(Y::AbstractArray,U::AbstractArray,
                      W::AbstractArray,X::AbstractArray)
     n, m = size(U);
	 mx = size(X,2);
     Ubar = zeros(m,mx);
     Ubar = Ubar + U'*X;
     for i = 1:n
         tmpW = W[i,:]
         tmpU = U[i,:]
         tmpX = X[i:i,:]
         Y[i,:] = Ubar'*tmpW;
         Ubar  -= tmpU .* tmpX
     end
end
#### Forward substitution (solve Lx = b) ####
function ss_forward!(X::AbstractArray, U::AbstractArray,
                     W::AbstractArray, B::AbstractArray)
    n, m = size(U)
    mx = size(B,2)
    Wbar = zeros(m, mx);
    for i = 1:n
        tmpU = U[i,:]
        tmpW = W[i,:]
        X[i,:] = (B[i:i,:] - tmpU'*Wbar)./(tmpU'*tmpW);
        Wbar += tmpW .* X[i:i,:];
    end
end
#### Backward substitution (solve L'x = b) ####
function ssa_backward!(X::AbstractArray, U::AbstractArray,
                       W::AbstractArray, B::AbstractArray)
    n, m = size(U)
    mx = size(B,2)
    Ubar = zeros(m,mx);
    for i = n:-1:1
        tmpU = U[i,:];
        tmpW = W[i,:];
        X[i,:] = (B[i:i,:] - tmpW'*Ubar)/(tmpU'*tmpW);
        Ubar += tmpU .* X[i:i,:];
    end
end
#### Logarithm of determinant ####
function ss_logdet(U::AbstractArray, W::AbstractArray)
    n = size(U, 1);
    a = 0;
    for i = 1:n
        a += log(dot(U[i,:],W[i,:]));
    end
    return a
end


text:
include("header.jl")

#=
# Fluid properties, flow quantities, and units
In this notebook, we will discuss *fluid properties* and *flow quantities*,
and the systems of units we use for them. First, some basic definitions:
- Fluid properties are properties of the material (the fluid)
- Flow quantities are characteristics of the flow of this fluid
We will also introduce some syntax we will use in this notebook and those that follow for dealing with units.
=#

# ### Set up the module
using MAE103

#=
We will generally focus on SI units, and these will be the default system
for displaying quantities. However, we also need to be familiar
with Imperial (sometimes called "English") units, which arise in many
situations in engineering. In the examples below, we will show that the tools
in these notebooks allow us to easily convert from one system to another.

For any quantity, we can see what the default units are by using `default_unit`.
For example, for pressure,
=#
default_unit(Pressure)
#-
#=
### Fluid properties
Let's start by discussing the basic properties of a fluid.

#### Density
The **density** provides a measure of the amount of fluid material per unit volume.
It is measured in units of mass/volume:
=#
default_unit(Density)

#=
For example, the density of water at a reference temperature of 15.6 degrees C is
=#
Density(Water)

#=
and for air at temperature of 15 degrees C and pressure of 1 atmosphere is
=#
Density(Air)

#=
Notice that water is around 800 times denser than air. This fact is very important
in fluid mechanics!
=#

#=
Remember that the density of gases, like air, depend on pressure and temperature (by
the ideal gas law), so the density of a gas may be sensitive to the local conditions.
This will happen in flows traveling relatively fast ([Q: compared to what?]). In such a case, density is a flow quantity,
not a fluid property, and we have to determine its value as part of the problem.

However, when a gas is traveling relatively slow [Q: compared to what?], we can
often treat a gas as having constant density, so it can be treated as a fluid property.

And liquid density is generally not sensitive. In fact, usually we can assume that the
density in a liquid is constant and uniform. We can treat
density in a liquid as a fluid property rather than a flow quantity.
=#
#-
#=
Also note that the density of seawater is larger (at average salt concentration)
than that of freshwater
=#
Density(Seawater)

#=
What if we want to set the density in other units? For example, in imperial
units, we would usually set it with lbm/ft^3. (In the notebooks, we use `lbm` for
pound (mass), and `lbf` is used for pound (force).) To set the units of a quantity,
rather than rely on the default units, we use the syntax to follow the number with
 `u"units"`. It will automatically convert it to the default units. For example, 2 lb/ft^3:
=#
Density(2u"lbm/ft^3")

#=
and to report the value of a quantity in other units, we use the `value`
function. For example, to get the density of water in lb/ft^3:
=#
value(Density(Water),u"lbm/ft^3")

#=
#### Viscosity
The internal friction in the fluid, called the *viscosity* (or, more specifically, the *dynamic viscosity*), given
by the Greek symbol $\mu$. [We can get this symbol by typing `\mu+TAB`.] The
viscosity controls the relationship between *shear stress*  $\tau$ (the frictional
force per unit area) and the *strain rate*, given by the gradient of velocity,
$\mathrm{d}u/\mathrm{d}y$. This latter quantity measures the difference in speeds in adjacent layers
of fluid, and larger differences suggest more shear stress. Viscosity is
the proportionality constant

$$\tau = \mu \dfrac{\mathrm{d}u}{\mathrm{d}y}$$

We will learn
much more about viscosity later, but for now, it is sufficient to
know that viscosity has units of kg/m/s:
=#
μw = Viscosity(Water)
#-
μa = Viscosity(Air)

#=
Note that water is much more viscous the air:
=#
μw/μa
#=
This is probably intuitive to you.
=#
#=
Note that both of these viscosities are much smaller than that of glycerin:
=#
Viscosity(Glycerin)
#-
#=
We will also occasionally make use of the ratio between viscosity and density.
This is called the *kinematic viscosity*. We use the symbol $\nu$ for this
[obtained by typing `\mu+TAB`.]
=#
νw = KinematicViscosity(Water)
#-
νa = KinematicViscosity(Air)

#=
#### Surface tension
The surface tension is a property associated with liquid interfaces, and particularly,
liquid interfaces with gases. The surface tension has units of force per unit length,
because if we imagine "cutting" a bit of the interface from the rest of it,
then this cut would form a perimeter of the snipped part of the interface, and
surface tension would act along this perimeter, representing how much the rest
of the liquid interface was pulling on it.

Surface tension also happens to have units of energy per unit area, so it
is sometimes referred to as "surface energy". Different liquids have different
surface energies, depending on the strength of their inter-molecular forces.
=#
SurfaceTension(Water)
#-
SurfaceTension(Glycerin)

#=
### Flow quantities
Now, let's discuss quantities that describe the fluid flow, or
at least, the *state* of the fluid. It is important to understand that
these quantities are, in general, *field quantities*: they vary from location
to location, and perhaps vary over time. So each of them should be
thought of as a *function* of the spatial coordinates, $(x,y,z)$ and time $t$.
Finding these functions is often our ultimate goal in solving a problem.
=#

#=
#### Pressure
Pressure represents the average force that the molecules exert per unit
area of surface. It is important to understand that pressure acts the same in every
direction. That is, it is an *isotropic* quantity. This means that, no matter
what the orientation of the surface, the pressure acts the same on it, and further,
it only acts *perpendicular* to the surface.

The default SI unit of pressure is the *Pascal* (Pa, equal to 1 N/m^2, or 1 kg/m/s^2).
=#
Pressure(20)

#=
But there are many other units for pressure in use. Some common ones. The
atmosphere (atm) represents the standard ambient pressure.
=#
p = Pressure(1u"atm")
#=
So 1 atm is about 101325 Pa, or 101.325 kPa. Let's see this written in other units
=#
value(p,u"psi")
#-
value(p,u"mmHg")
#=
These are all *absolute* pressures. For many of the flows we will
study, the absolute pressure will not be important. Consider,
for example, the flow of water through a pipe. Only the *difference* between
the pressure at the inlet and the pressure at the outlet matters for
driving this flow; the absolute pressure is irrelevant. So
we can also define a `PressureDifference`, measured in the
same units:
=#
PressureDifference(50)

#=
Pressure is a flow quantity. In general, we need to find its values
as part of the problem.
=#

#=
#### Velocity
The velocity of the fluid describes how fast it is moving in each direction,
and has units of length/time.
=#
u = Velocity(20)
#-
#=
In other units,
=#
value(u,u"cm/s")
#-
value(u,u"ft/s")

#=
We know that velocity is a vector, so it has, in general, three components $(u,v,w)$.
It is also a field quantity. This means that each of the three components
depends on $(x,y,z)$ and $t$. That's a lot of detail to determine, but once
we determine it, we know *everything* about a fluid flow.
=#