text:
export σ, sigmoid, hardσ, hardsigmoid, hardtanh, relu, leakyrelu, relu6, rrelu, elu, gelu, swish, selu, celu, softplus, softsign, logσ,
       logsigmoid, logcosh, mish, tanhshrink, softshrink, thresholdrelu, trelu, lisht

## Activation functions
#
# Some of activation functions have its wrapper function for GPU in CuArrays.jl.
# https://github.com/JuliaGPU/CuArrays.jl/issues/614

"""
    σ(x) = 1 / (1 + exp(-x))

Classic [sigmoid](https://en.wikipedia.org/wiki/Sigmoid_function) activation
function.
"""
function σ(x)
    t = exp(-abs(x))
    ifelse(x ≥ 0, inv(one(t) + t), t / (one(t) + t))
end
const sigmoid = σ

"""
    hardσ(x, a=0.2) = max(0, min(1.0, a * x + 0.5))

Segment-wise linear approximation of sigmoid.
See [BinaryConnect: Training Deep Neural Networks withbinary weights during propagations](https://arxiv.org/abs/1511.00363).
"""
hardσ(x, a=0.2) = oftype(x/1, max(zero(x/1), min(one(x/1), oftype(x/1,a) * x + oftype(x/1,0.5))))
const hardsigmoid = hardσ

"""
    logσ(x)

Return `log(σ(x))` which is computed in a numerically stable way.

    julia> logσ(0)
    -0.6931471805599453
    julia> logσ.([-100, -10, 100])
    3-element Array{Float64,1}:
     -100.0
      -10.000045398899218
       -3.720075976020836e-44
"""
logσ(x) = -softplus(-x)
const logsigmoid = logσ

"""
    hardtanh(x) = max(-1, min(1, x))

Segment-wise linear approximation of tanh. Cheaper  and  more  computational  efficient version of tanh.
See [Large Scale Machine Learning](https://ronan.collobert.com/pub/matos/2004_phdthesis_lip6.pdf).
"""
hardtanh(x) = max(-one(x), min( one(x), x))

"""
    relu(x) = max(0, x)

[Rectified Linear Unit](https://en.wikipedia.org/wiki/Rectifier_(neural_networks))
activation function.
"""
relu(x) = max(zero(x), x)

"""
    leakyrelu(x, a=0.01) = max(a*x, x)

Leaky [Rectified Linear Unit](https://en.wikipedia.org/wiki/Rectifier_(neural_networks))
activation function.
You can also specify the coefficient explicitly, e.g. `leakyrelu(x, 0.01)`.
"""
leakyrelu(x, a = oftype(x / 1, 0.01)) = max(a * x, x / one(x))

"""
    relu6(x) = min(max(0, x), 6)

[Rectified Linear Unit](https://en.wikipedia.org/wiki/Rectifier_(neural_networks))
activation function capped at 6.
See [Convolutional Deep Belief Networks on CIFAR-10](https://www.cs.toronto.edu/~kriz/conv-cifar10-aug2010.pdf)
"""
relu6(x) = min(relu(x), oftype(x, 6))

"""
    rrelu(x, l=1/8, u=1/3) = max(a*x, x)

    a = randomly sampled from uniform distribution U(l, u)

Randomized Leaky [Rectified Linear Unit](https://arxiv.org/abs/1505.00853)
activation function.
You can also specify the bound explicitly, e.g. `rrelu(x, 0.0, 1.0)`.
"""
function rrelu(x, l = 1 / 8.0, u = 1 / 3.0)
    a = oftype(x / 1, (u - l) * rand() + l)
    return leakyrelu(x, a)
end

"""
    elu(x, α=1) =
      x > 0 ? x : α * (exp(x) - 1)

Exponential Linear Unit activation function.
See [Fast and Accurate Deep Network Learning by Exponential Linear Units](https://arxiv.org/abs/1511.07289).
You can also specify the coefficient explicitly, e.g. `elu(x, 1)`.
"""
elu(x, α = one(x)) = ifelse(x ≥ 0, x / one(x), α * (exp(x) - one(x)))

"""
    gelu(x) = 0.5x * (1 + tanh(√(2/π) * (x + 0.044715x^3)))

[Gaussian Error Linear Unit](https://arxiv.org/abs/1606.08415)
activation function.
"""
function gelu(x)
    p = oftype(x / 1, π)
    λ = oftype(x / 1, √(2 / p))
    α = oftype(x / 1, 0.044715)
    h = oftype(x / 1, 0.5)
    h * x * (one(x) + tanh(λ * (x + α * x^3)))
end

"""
    swish(x) = x * σ(x)

Self-gated activation function.
See [Swish: a Self-Gated Activation Function](https://arxiv.org/abs/1710.05941).
"""
swish(x) = x * σ(x)

"""
    lisht(x) = x * tanh(x)

Non-Parametric Linearly Scaled Hyperbolic Tangent Activation Function.
See [LiSHT](https://arxiv.org/abs/1901.05894)
"""
lisht(x) = x * tanh(x)

"""
    selu(x) = λ * (x ≥ 0 ? x : α * (exp(x) - 1))

    λ ≈ 1.0507
    α ≈ 1.6733

Scaled exponential linear units.
See [Self-Normalizing Neural Networks](https://arxiv.org/abs/1706.02515).
"""
function selu(x)
  λ = oftype(x / 1, 1.0507009873554804934193349852946)
  α = oftype(x / 1, 1.6732632423543772848170429916717)
  λ * ifelse(x > 0, x / one(x), α * (exp(x) - one(x)))
end

"""
    celu(x, α=1) =
        (x ≥ 0 ? x : α * (exp(x/α) - 1))

Continuously Differentiable Exponential Linear Units
See [Continuously Differentiable Exponential Linear Units](https://arxiv.org/abs/1704.07483).
"""
celu(x, α = one(x)) = ifelse(x ≥ 0, x / one(x), α * (exp(x/α) - one(x)))

"""
    trelu(x, theta = 1.0) = x > theta ? x : 0

Threshold Gated Rectified Linear.
See [ThresholdRelu](https://arxiv.org/abs/1402.3337)
"""
trelu(x, theta = one(x)) = ifelse(x> theta, x, zero(x))
const thresholdrelu = trelu

"""
    softsign(x) = x / (1 + |x|)

See [Quadratic Polynomials Learn Better Image Features](http://www.iro.umontreal.ca/~lisa/publications2/index.php/attachments/single/205).
"""
softsign(x) = x / (one(x) + abs(x))

"""
    softplus(x) = log(exp(x) + 1)

See [Deep Sparse Rectifier Neural Networks](http://proceedings.mlr.press/v15/glorot11a/glorot11a.pdf).
"""
softplus(x) = ifelse(x > 0, x + log1p(exp(-x)), log1p(exp(x)))

"""
    logcosh(x)

Return `log(cosh(x))` which is computed in a numerically stable way.
"""
logcosh(x) = x + softplus(-2x) - log(oftype(x, 2))

"""
    mish(x) = x * tanh(softplus(x))

Self Regularized Non-Monotonic Neural Activation Function.
See [Mish: A Self Regularized Non-Monotonic Neural Activation Function](https://arxiv.org/abs/1908.08681).
"""
mish(x) = x * tanh(softplus(x))

"""
    tanhshrink(x) = x - tanh(x)

See [Tanhshrink Activation Function](https://www.gabormelli.com/RKB/Tanhshrink_Activation_Function).
"""
tanhshrink(x) = x - tanh(x)

"""
    softshrink(x, λ=0.5) =
        (x ≥ λ ? x - λ : (-λ ≥ x ? x + λ : 0))

See [Softshrink Activation Function](https://www.gabormelli.com/RKB/Softshrink_Activation_Function).
"""
softshrink(x, λ = oftype(x/1, 0.5)) = min(max(zero(x), x - λ), x + λ)

# Provide an informative error message if activation functions are called with an array
for f in (:σ, :hardσ, :logσ, :hardtanh, :relu, :leakyrelu, :relu6, :rrelu, :elu, :gelu, :swish, :lisht, :selu, :celu, :trelu, :softsign, :softplus, :logcosh, :mish, :tanhshrink, :softshrink)
    @eval $(f)(x::AbstractArray, args...) =
      error("Use broadcasting (`", $(string(f)), ".(x)`) to apply activation functions to arrays.")
end


text:
@testset "Transforms" begin
    println("VITUIX TRANSFORM TESTS")

    test_set = [
        (src="T/polar", target="T/polar_r20", freqrange=(1kHz, 20kHz), T=Dict(:r => 20.0), max_mag=3.2, mae_mag=0.4, max_phi=180.0, mae_phi=6.5),
        (src="T/polar", target="T/polar_r10_t15", freqrange=(1kHz, 20kHz), T=Dict(:r => 10.0, :t => 15.0), max_mag=4.7, mae_mag=0.3, max_phi=180.0, mae_phi=8.5),
        (src="T/polar_noXO", target="T/polar_noXO_r20", freqrange=(1kHz, 20kHz), T=Dict(:r => 20.0), max_mag=3.2, mae_mag=0.4, max_phi=180.0, mae_phi=6.3),
        (src="T/polar_noXO", target="T/polar_noXO_r-80_t-80", freqrange=(1kHz, 20kHz), T=Dict(:r => -80.0, :t => -80.0), max_mag=8.5, mae_mag=0.4, max_phi=180.0, mae_phi=6.0),

        (src="M/polar", target="M/polar_r10", freqrange=(0kHz, 20kHz), T=Dict(:r => 10.0), max_mag=4.3, mae_mag=0.25, max_phi=180.0, mae_phi=3.2),
        (src="M/polar", target="M/polar_r170", freqrange=(0kHz, 20kHz), T=Dict(:r => 170.0), max_mag=5.0, mae_mag=0.25, max_phi=180.0, mae_phi=4.0),
        (src="M/polar", target="M/polar_y-118", freqrange=(0kHz, 20kHz), T=Dict(:y => -0.118), max_mag=13.0, mae_mag=0.25, max_phi=180.0, mae_phi=0.7),
        (src="M/polar", target="M/polar_y-118_r170", freqrange=(0kHz, 20kHz), T=Dict(:y => -0.118, :r => 170.0), max_mag=13.0, mae_mag=0.25, max_phi=180.0, mae_phi=2.2),
        (src="M/polar", target="M/polar_y-118_r180", freqrange=(0kHz, 20kHz), T=Dict(:y => -0.118, :r => 180.0), max_mag=13.0, mae_mag=0.25, max_phi=180.0, mae_phi=2.2),
        (src="M/polar", target="M/polar_x95_y-118_z150_r12_t-15", freqrange=(0kHz, 20kHz), T=Dict(:x => 0.095, :y => -0.118, :z => 0.150, :r => 12.0, :t => -15.0), max_mag=6.0, mae_mag=0.3, max_phi=180.0, mae_phi=4.0),
        (src="M/polar_noXO", target="M/polar_noXO_r170", freqrange=(0kHz, 20kHz), T=Dict(:r => 170.0), max_mag=5.0, mae_mag=0.25, max_phi=180.0, mae_phi=4.0),
        (src="M/polar_noXO", target="M/polar_noXO_r-170", freqrange=(0kHz, 20kHz), T=Dict(:r => -170.0), max_mag=5.0, mae_mag=0.25, max_phi=180.0, mae_phi=3.5),
        (src="M/polar_noXO", target="M/polar_noXO_y-118_r170", freqrange=(0kHz, 20kHz), T=Dict(:y => -0.118, :r => 170.0), max_mag=13.0, mae_mag=0.25, max_phi=180.0, mae_phi=2.2),

        (src="W/polar", target="W/polar_y-348", freqrange=(0kHz, 2kHz), T=Dict(:y => -0.348), max_mag=3.5, mae_mag=0.15, max_phi=180.0, mae_phi=1.2),
        (src="W/polar", target="W/polar_y-348_r14", freqrange=(0kHz, 2kHz), T=Dict(:y => -0.348, :r => 14.0), max_mag=5.6, mae_mag=0.15, max_phi=180.0, mae_phi=1.2),
        (src="W/polar", target="W/polar_x-150_y-348_z-35_r-90_t3", freqrange=(0kHz, 2kHz), T=Dict(:x => -0.15, :y => -0.348, :z => -0.035, :r => -90.0, :t => 3.0), max_mag=4.5, mae_mag=0.15, max_phi=180.0, mae_phi=1.5),
        (src="W/polar_noXO", target="W/polar_noXO_y-348_r14", freqrange=(0kHz, 2kHz), T=Dict(:y => -0.348, :r => 14.0), max_mag=5.6, mae_mag=0.15, max_phi=180.0, mae_phi=1.2),
    ]

    for e in test_set
        println("  \"$(e.src)\" -> \"$(e.target)\"")
        src, target = joinpath(test_dir, e.src), joinpath(test_dir, e.target)
        print("\tA:")
        @test test_driver_transform(src, target, e.freqrange, compare_magdb, e.max_mag, e.mae_mag, false; e.T...)
        print("\tΦ:")
        @test test_driver_transform(src, target, e.freqrange, compare_phasedeg, e.max_phi, e.mae_phi, false; e.T...)
    end
end

ext: jl


text:
# Run like a user might run it
@testset "End-to-end" begin

    @testset "Particles with same shape" begin

        # 2D acoustics
        circle1 = Sphere((0.0,0.0),1.0)
        circle2 = Sphere((0.0,5.0),2.0)
        a_p = Acoustic(2; ρ = 1.0, c = 1.0)
        a = Acoustic(2; ρ = 0.3, c = 0.2)
        particles = [Particle(a_p,circle1), Particle(a_p,circle2)]

        source = plane_source(a; position = [0.0,0.0], direction = [1.0,0.0])
        sim = FrequencySimulation(particles,source)

        # show(sim);

        result = run(sim, SVector(1.0,2.0), 0.1)
        result = run(particles, source, SVector(1.0,2.0), 0.1)
        result = 3.2*result + result*4.0im + 0.3+4.0im # changes result.field
        3.2 + result

        @test field(result)[1] == result.field[1][1] # returns
        @test field(result,1,1) == result.field[1][1] # returns

        # run simulation over a region rather than a specif set of x points
        ω = 0.1
        region = bounding_box([p.shape for p in particles])
        result = run(particles, source, region, [ω]; res=10)

        # 3D acoustics
        s1 = Sphere((1.0,-1.0,2.0),1.0)
        s2 = Sphere((0.0,3.0,0.0),2.0)
        a_p = Acoustic(3; ρ = 1.0, c = 1.0)
        a = Acoustic(3; ρ = 0.3, c = 0.2)

        particles = [Particle(a_p,s1), Particle(a_p,s2)]
        source = plane_source(a; position = [0.0,0.0,0.0], direction = [0.0,0.0,1.0])

        x = SVector(-3.0,2.0,3.0)
        ω = 0.1
        sim = FrequencySimulation(particles,source)
        result = run(particles, source, x, ω)
        result = 3.2*result + result*4.0im + 0.3+4.0im # changes result.field
        3.2 + result

        @test field(result)[1] == result.field[1][1] # returns
        @test field(result,1,1) == result.field[1][1] # returns

    end

    @testset "Particles with different shape" begin
        circle = Sphere((0.0,0.0),1.0)
        rect = Box((2.0,2.0),(3.0,2.0))
        a_p = Acoustic(1.0,1.0,2)
        a = Acoustic(0.3,0.2,2)
        particles = [Particle(a_p,circle), Particle(a_p,rect)]
        source = plane_source(a,[0.0,0.0],[1.0,0.0])
        sim = FrequencySimulation(particles,source)
        @test true
    end

    include("time_one_particle.jl")

end


text:
using AcceleratedDubins

## Parameters of the vehicle
params = [5., 15., 2.6, 4] ## v_min v_max a_max -a_min
r_min = 1. # minimum turning radius
r_max = 5. # maximum turning radius

## Generate random configurations
start = [10*rand(), -10+5*rand(), 2*pi*rand()]
stop = [10+20*rand(), -10+5*rand(), 2*pi*rand()]

## try radii combinations and get resulting path time
    path_fastest, _ = fastest_path(start, stop, radii_samples_exp(r_min, r_max, 3), params)
    time_fastest = path_time(path_fastest, params)
## try to optimize path using Optim
    path_optimized, _ = optimized_path(start, stop, [r_min, r_max], params)
    time_optimized = path_time(path_optimized, params)
## compute shortest path (original Dubins with r_min)
    path_shortest, _ = retrieve_path(start, stop, [r_min, r_min], params)
    time_shortest = path_time(path_shortest, params)

@show time_fastest, time_optimized, time_shortest

## sample the path for plotting
times, speeds = speed_profile(path_fastest, params)


text:
#=
    The File contains tests for find_basis Algorithms
=#

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

# Tests the given algorithm for correctness by comparing
# it with the output of standard find_basis_1 on random samples
# The latter is assumed to be correct
#
# Yet some tests are hardcoded just in case
function test_find_basis(algorithm)

    # nilpotent space
    set0 = [
        from_dense([0 1; 0 0;], Nemo.QQ),
    ]
    ans0 = linear_span!([
        from_dense([0 1; 0 0;], Nemo.QQ)
    ])

    # almost nilpotent space
    set1 = [
        from_dense([0 1; 0 0;], Nemo.QQ),
        from_dense([0 0; 0 1;], Nemo.QQ),
    ]
    ans1 = linear_span!([
        from_dense([0 1; 0 0;], Nemo.QQ),
        from_dense([0 0; 0 1;], Nemo.QQ),
    ])

    # some space
    set2 = [
        from_dense([1 2 3; 0 1 0; 0 0 0;], Nemo.QQ),
        from_dense([1 0 0; 0 0 0; 0 2 1;], Nemo.QQ)
    ]
    ans2 = linear_span!([
        from_dense([1 0 0; 0 0 0; 0 0 1;], Nemo.QQ),
        from_dense([0 1 0; 0 0 0; 0 0 0;], Nemo.QQ),
        from_dense([0 0 1; 0 0 0; 0 0 -1//3;], Nemo.QQ),
        from_dense([0 0 0; 0 1 0; 0 0 0;], Nemo.QQ),
        from_dense([0 0 0; 0 0 0; 0 1 0;], Nemo.QQ)
    ])

    # complete space
    set3 = [
        from_dense([1 0 0; 2 1 3; 1 0 0;], Nemo.QQ),
        from_dense([0 1 0; 1 0 0; 0 0 1;], Nemo.QQ)
    ]
    ans3 = linear_span!([
        from_dense([1 0 0; 0 0 0; 0 0 0;], Nemo.QQ),
        from_dense([0 1 0; 0 0 0; 0 0 0;], Nemo.QQ),
        from_dense([0 0 1; 0 0 0; 0 0 0;], Nemo.QQ),
        from_dense([0 0 0; 1 0 0; 0 0 0;], Nemo.QQ),
        from_dense([0 0 0; 0 1 0; 0 0 0;], Nemo.QQ),
        from_dense([0 0 0; 0 0 1; 0 0 0;], Nemo.QQ),
        from_dense([0 0 0; 0 0 0; 1 0 0;], Nemo.QQ),
        from_dense([0 0 0; 0 0 0; 0 1 0;], Nemo.QQ),
        from_dense([0 0 0; 0 0 0; 0 0 1;], Nemo.QQ),
    ])

    for (set, ans) in zip([set0, set1, set2, set3], [ans0, ans1, ans2, ans3])
        V = algorithm(deepcopy(set))
        @test check_inclusion!(deepcopy(ans), deepcopy(V)) && check_inclusion!(V, deepcopy(ans))
    end

    # --------
    # a series of random tests

    ZZ  = GF(2^31 - 1)

    dims      = [3, 3, 4, 4, 5, 5, 5, 5, 8, 8, 8, 8, 10, 10, 10]
    sizes     = [1, 3, 10]
    densities = [0.05, 0.1, 0.2, 0.5]

    for n in dims
        for sz in sizes
            for d in densities
                As = [
                    random_sparsik((n, n), ZZ, density=d)
                    for _ in 1:sz
                ]
                ans = find_basis_1(deepcopy(As))
                V   = algorithm(deepcopy(As))
                @test check_inclusion!(ans, deepcopy(V)) && check_inclusion!(V, deepcopy(ans))
            end
        end
    end

end

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

@testset "basis, just in case" begin
    # This must always pass
    test_find_basis(find_basis_1)

end

@testset "basis, Wied. algorithm" begin

    # the algorithm is probabilistic, we should probably apply_vectors once again
    # at the end to obtain the answer, just as in find_basis_1_beta

    # test_find_basis(find_basis_2)

end

@testset "basis, Beta algorithm" begin

    test_find_basis(find_basis_1_β)

end

@testset "basis, Gamma-Beta algorithm" begin

    # test_find_basis(find_basis_1_γβ)

end

@info "OK"


text:
# Test variational information

using Base.Test
using Clustering

a1 = [1, 1, 1, 2, 2, 2, 3, 3, 3, 3]
a2 = [1, 1, 1, 1, 2, 2, 2, 2, 2, 2]

@test isapprox(varinfo(3, a1, 3, a1), 0.0, atol=1.0e-12)
@test isapprox(varinfo(2, a2, 2, a2), 0.0, atol=1.0e-12)

v = varinfo(3, a1, 2, a2)
v_ = varinfo(2, a2, 3, a1)
@test 0.0 < v < log(3)
@test isapprox(v, v_)

a1 = [1, 2, 3, 4, 5]
a2 = [1, 1, 1, 1, 1]

@test isapprox(varinfo(5, a1, 1, a2), log(5))
@test isapprox(varinfo(1, a2, 5, a1), log(5))

a1 = [1, 1, 1, 2, 2, 2]
a2 = [2, 2, 2, 1, 1, 1]
@test isapprox(varinfo(2, a1, 2, a2), 0.0, atol=1.0e-12)
@test isapprox(varinfo(2, a1, 3, a2), 0.0, atol=1.0e-12)
@test isapprox(varinfo(4, a1, 3, a2), 0.0, atol=1.0e-12)


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# ╔═╡ 45376f68-2552-11ec-32d7-a5de0d6b6882
begin

	# Setup the Julia environment -
	using PlutoUI
	using MAT

	# setup paths -
	const _PATH_TO_NOTEBOOK = pwd()
	const _PATH_TO_DATA = joinpath(_PATH_TO_NOTEBOOK,"data")
	const _PATH_TO_FIGS = joinpath(_PATH_TO_NOTEBOOK,"figs")

	# return -
	nothing
end

# ╔═╡ 16a831a8-14b2-4d1e-86e7-a28f21d0dbd4
md"""
## Introduction to Metabolic Engineering

Metabolic engineering is the practice of optimizing genetic and regulatory processes within cells to increase the cell's production of a desired molecule or protein of interest.

Metabolic engineers manipulate the biochemical networks used by cells to convert raw materials into molecules necessary for the cell's survival. Metabolic engineering specifically seeks to:

1. Mathematically model biochemical networks, calculate the yield (product divided substrate) of useful products, and identify parts of the network that constrain the production of the products of interest.
1. Use genetic engineering techniques to modify the biochemical network in order to relieve constraints limiting production. The modified network can then be modeled to calculate the new product yield, and to identify new constraints (back to 1).

Resources for biochemical network information:
* [Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000 Jan 1;28(1):27-30. doi: 10.1093/nar/28.1.27. PMID: 10592173; PMCID: PMC102409.](https://www.genome.jp/kegg/)

* [Karp, Peter D et al. “The BioCyc collection of microbial genomes and metabolic pathways.” Briefings in bioinformatics vol. 20,4 (2019): 1085-1093. doi:10.1093/bib/bbx085](https://pubmed.ncbi.nlm.nih.gov/29447345/)

* [Gama-Castro, Socorro et al. “RegulonDB version 9.0: high-level integration of gene regulation, coexpression, motif clustering and beyond.” Nucleic acids research vol. 44,D1 (2016): D133-43. doi:10.1093/nar/gkv1156](https://pubmed.ncbi.nlm.nih.gov/26527724/)

"""

# ╔═╡ 716d10ed-38bd-4654-9a54-33dfba020591
md"""
__Fig 1.__ The overall metabolic map from the KEGG database. Each dot (_node_) is a metabolite, each line (_edge_) is a metabolic reaction.
"""

# ╔═╡ ad2d1e2b-4222-43d0-819a-4a7d14746327
PlutoUI.LocalResource("./figs/KEGG-map01100.png")

# ╔═╡ cef733a0-5541-4338-bd28-0916b697fb19
md"""
## What is Flux Balance Analysis (FBA)?
Flux balance analysis (FBA) is a mathematical modeling and analysis approach which computes the flow (or _flux_) of carbon and energy throughout a metabolic network. FBA, a member of the constraint based family of mathematical modeling tools, is a widely used approach to compute metabolic flux. However, there are alternatives to FBA, such as metabolic flux analysis (MFA), but these alternatives vary more in the solution approach than the structure of the estimation problem.

Let's look at the following reference to better understand the different components of a flux balance anaysis problem:

* [Orth, J., Thiele, I. & Palsson, B. What is flux balance analysis?. Nat Biotechnol 28, 245–248 (2010). https://doi.org/10.1038/nbt.1614](https://www.ncbi.nlm.nih.gov/labs/pmc/articles/PMC3108565/)

Computational tools for working with flux balance analysis models:

* [Heirendt, Laurent et al. “Creation and analysis of biochemical constraint-based models using the COBRA Toolbox v.3.0.” Nature protocols vol. 14,3 (2019): 639-702. doi:10.1038/s41596-018-0098-2](https://pubmed.ncbi.nlm.nih.gov/30787451/)

#### Flux balance analysis problem structure
The FBA problem is typically encoded as a linear programming (LP) problem of the form:

$$\max_{v}\sum_{i=1}^{\mathcal{R}}c_{i}v_{i}$$

subject to the constraints:

$$\begin{eqnarray}
\mathbf{S}\mathbf{v} &=& \mathbf{0}\\
\mathcal{L}\leq&\mathbf{v}&\leq\mathcal{U}\\
\dots &=& \dots
\end{eqnarray}$$

where $\mathbf{S}$ denotes the stoichiometric matrix, $c_{i}$ denote the objective coefficients,
$\mathbf{v}$ denotes the metabolic flux (the unknown that we are trying to estimate), and
$\mathcal{L}$ ($\mathcal{U}$) denote the permissible lower (upper) bounds on the _unkown_ metabolic flux. The first set of constraints enforces conservation of mass, while the second imparts thermodynamic and kinetic information into the calculation. Finally, there are potentially other types of constraints (both linear and nonlinear) that can be used in this type of problem (we will not cover these here, but these additional constraints may be important in specific applications).

### What is a stoichiometric matrix?
The basis for flux balance calculations is the stoichiometric matrix.
The stochiometric matrix is a digitial representation of the biochemistry that can occur insides cells (such as that shown in Fig. 1). The stochiometric matrix, denoted by $\mathbf{S}$, is a $\mathcal{M}\times\mathcal{R}$ array of stoichiometric coefficients, denoted by the symbol $\sigma_{ij}$ ($i$ is the row index, $j$ is the col index). Each of the $\mathcal{M}$ rows of $\mathbf{S}$ describes a particular metabolite (_node_ in the metabolic network) while each of the $\mathcal{R}$ columns corresponds to a metabolic reaction (_edge_ in the metabolic network). Thus:

* A stoichiometric coefficient $\sigma_{ij}$ > 0 implies that metabolite $i$ is __produced__ by reaction $j$
* A stoichiometric coefficient $\sigma_{ij}$  = 0 implies that metabolite $i$ is __not connected__ to reaction $j$
* A stoichiometric coefficient $\sigma_{ij}$ < 0 implies that metabolite $i$ is __consumed__ by reaction $j$
"""

# ╔═╡ 679f15d4-e94d-43eb-8342-701111c46201
md"""

__Example__: Core stoichiometric matrix from _Escherichia coli_: [Orth, Jeffrey D et al. “Reconstruction and Use of Microbial Metabolic Networks: the Core Escherichia coli Metabolic Model as an Educational Guide.” EcoSal Plus vol. 4,1 (2010): 10.1128/ecosalplus.10.2.1. doi:10.1128/ecosalplus.10.2.1](https://pubmed.ncbi.nlm.nih.gov/26443778/)
"""

# ╔═╡ 0f683f10-71d7-477c-9b7c-49545e60fa8a
# load the stoichiometric matrix, and mess around with it -
begin

	# what is the model name -
	model_file_name = "modelReg.mat"
	model_name = "modelReg"

	# where is model file?
	_PATH_TO_MODEL_FILE = joinpath(_PATH_TO_DATA,model_file_name)

	# load the mat file -> get the cobra_dictionary
	file = matopen(_PATH_TO_MODEL_FILE)
    cobra_dictionary = read(file, model_name)
    close(file)
end

# ╔═╡ 095a2f1a-6ae3-451a-bc3f-7dad9a7834ea
keys(cobra_dictionary)

# ╔═╡ 5d7c052a-73d7-487a-939e-f813561daff9
# get the stoichiometric matrix
stm_sparse = cobra_dictionary["S"]

# ╔═╡ b3e23596-0b67-4ed1-8caa-6353307ee315
stm_full = Matrix(stm_sparse)

# ╔═╡ 5f19c57c-e4ac-4c8c-9bc0-f6a015c8fb3f
mets = cobra_dictionary["metFormulas"]

# ╔═╡ 1062535d-3453-4a00-89f5-66060c150dbf
rxns = cobra_dictionary["rxns"]

# ╔═╡ 5578b9b1-8fbd-4c22-a89d-cd39ea63fe87
md"""
Additional resources for the stoichiometric matrix

* [Lecture 3: The genome reconstruction process. Systems Biology: Constraint-based Reconstruction and Analysis, Cambridge University Press, 2015](https://www.youtube.com/watch?v=jcRVSJj3XP8&feature=youtu.be)

* [Lecture 9: Properties of the Stoichiometric Matrix. Systems Biology: Constraint-based Reconstruction and Analysis, Cambridge University Press, 2015](https://www.youtube.com/watch?v=0nsjKx6a7X4&feature=youtu.be)

"""

# ╔═╡ 6629d6a5-af02-40c4-a7a3-a7146bf6638b
md"""
### What are flux bounds constraints?
Flux bounds constraints control the permissible ranges that a metabolic reaction rate (flux) can have. Flux bounds constraints can incorporate thermodynamic information (whether a reaction is reversible) as well as kinetic information.

Resources for thermodynamics in the contraint based world:

* [Peres, Sabine, and Vincent Fromion. “Thermodynamic Approaches in Flux Analysis.” Methods in molecular biology (Clifton, N.J.) vol. 2088 (2020): 359-367. doi:10.1007/978-1-0716-0159-4_17](https://pubmed.ncbi.nlm.nih.gov/31893383/)

* [Flamholz, Avi et al. “eQuilibrator--the biochemical thermodynamics calculator.” Nucleic acids research vol. 40,Database issue (2012): D770-5. doi:10.1093/nar/gkr874](https://pubmed.ncbi.nlm.nih.gov/22064852/)

Let's do some sample calculations using [eQuilibrator](https://equilibrator.weizmann.ac.il)

Thermodynamic information tells you the directionality of the reaction (is the reaction reversible), but to get the permissible magnitudes of the reaction we need to compute the kinetics. We (and others) formulate flux bounds as the product of succesive corrections t the maximum possible rate e.g., for an irrereversible reaction:

$$0\leq{v_{i}}\leq{V_{max,j}^{\circ}}\left(\frac{\epsilon}{\epsilon^{\circ}}\right)\theta_{j}\left(\dots\right){f_{j}\left(\dots\right)}$$

where $V_{max,j}^{\circ}$ denotes the maximum reaction velocity computed for enzyme $j$ at some characteristic enzyme concetration (and full activity), the ratio $\epsilon/\epsilon^{\circ}$ is a correction for enzyme concentration, $\theta_{j}\left(\dots\right)\in\left[0,1\right]$ is an enzyme activity function (or measurement) and $f_{j}\left(\dots\right)$ is a function describing the substrate dependence of the reaction rate $j$. Both $\theta_{j}\left(\dots\right)$ and $f_{j}\left(\dots\right)$ could have associated parameters, e.g., saturation or binding constants, etc.

"""

# ╔═╡ baedaf07-7236-432c-a67b-bd7eaa99bb7c
md"""
### What is the objective of an _E. coli_ cell?

__Short answer__: maximize the growth rate

__Longer answer__: This question was raging in the metabolic engineering community twenty years ago, but is mostly forgotten today. First, there are long standing arguments from ecology (or even from the perspective of microeconomics) as to why _E.coli_ would be goal oriented. However, Palsson and coworkers largely put this question to rest (at least for _E.coli_) with the publication:

* [Ibarra, Rafael U et al. “Escherichia coli K-12 undergoes adaptive evolution to achieve in silico predicted optimal growth.” Nature vol. 420,6912 (2002): 186-9. doi:10.1038/nature01149](https://pubmed.ncbi.nlm.nih.gov/12432395/)

__What about cell-free networks__?
No cells, but most of the biochemistry. Is there an objective? This is an open question. However, initial work in this area suggests the maximization of translation, at least for single gene systems:

* [Vilkhovoy, Michael et al. “Sequence Specific Modeling of E. coli Cell-Free Protein Synthesis.” ACS synthetic biology vol. 7,8 (2018): 1844-1857. doi:10.1021/acssynbio.7b00465](https://pubmed.ncbi.nlm.nih.gov/29944340/)

"""

# ╔═╡ fd376fb1-1617-48a8-877b-31aa8f85476a
md"""
## Flux balance analysis case studies and examples

* [Wayman, Joseph A et al. “Improving designer glycan production in Escherichia coli through model-guided metabolic engineering.” Metabolic engineering communications vol. 9 e00088. 29 Mar. 2019, doi:10.1016/j.mec.2019.e00088](https://pubmed.ncbi.nlm.nih.gov/31008057/)

"""

# ╔═╡ 774dd760-cebf-44ff-b753-e54353c1317e
TableOfContents(title="CHEME 5430: Constraint Based Methods and Models", depth=4)

# ╔═╡ debd1859-15ce-49ec-b45a-4707695f5c8a
html"""
<style>
main {
    max-width: 1200px;
    margin: auto;
	width: 90%;
    font-family: "Roboto, monospace";
}

a {
    color: blue;
    text-decoration: none;
}

.H1 {
    padding: 0px 30px;
}
</style>
"""

# ╔═╡ 8904e3f1-88ad-424e-8ac2-e4c8dfd01c4a
html"""
<script>
var section = 0;
var subsection = 0;
var headers = document.querySelectorAll('h2, h3');
for (var i=0; i < headers.length; i++) {
    var header = headers[i];
    var text = header.innerText;
    var original = header.getAttribute("text-original");
    if (original === null) {
        // Save original header text
        header.setAttribute("text-original", text);
    } else {
        // Replace with original text before adding section number
        text = header.getAttribute("text-original");
    }
    var numbering = "";
    switch (header.tagName) {
        case 'H2':
            section += 1;
            numbering = section + ".";
            subsection = 0;
            break;
        case 'H3':
            subsection += 1;
            numbering = section + "." + subsection;
            break;
    }
    header.innerText = numbering + " " + text;
};
</script>
"""

# ╔═╡ d1bd21e8-86c0-4607-953b-d33a7b964ecc
html"""
<style>
main {
    max-width: 1200px;
  	width: 45%;
	margin: auto;
	font-family: "Roboto, monospace";
}

a {
	color: blue;
	text-decoration: none;
}

.H1 {
    padding: 0px 30px;
}
</style>
<script>
var section = 0;
var subsection = 0;
var headers = document.querySelectorAll('h2, h3');
for (var i=0; i < headers.length; i++) {
    var header = headers[i];
    var text = header.innerText;
    var original = header.getAttribute("text-original");
    if (original === null) {
        // Save original header text
        header.setAttribute("text-original", text);
    } else {
        // Replace with original text before adding section number
        text = header.getAttribute("text-original");
    }
    var numbering = "";
    switch (header.tagName) {
        case 'H2':
            section += 1;
            numbering = section + ".";
            subsection = 0;
            break;
        case 'H3':
            subsection += 1;
            numbering = section + "." + subsection;
            break;
    }
    header.innerText = numbering + " " + text;
};
</script>
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
MAT = "~0.10.1"
PlutoUI = "~0.7.14"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Blosc]]
deps = ["Blosc_jll"]
git-tree-sha1 = "84cf7d0f8fd46ca6f1b3e0305b4b4a37afe50fd6"
uuid = "a74b3585-a348-5f62-a45c-50e91977d574"
version = "0.7.0"

[[Blosc_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Lz4_jll", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "e747dac84f39c62aff6956651ec359686490134e"
uuid = "0b7ba130-8d10-5ba8-a3d6-c5182647fed9"
version = "1.21.0+0"

[[BufferedStreams]]
deps = ["Compat", "Test"]
git-tree-sha1 = "5d55b9486590fdda5905c275bb21ce1f0754020f"
uuid = "e1450e63-4bb3-523b-b2a4-4ffa8c0fd77d"
version = "1.0.0"

[[CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "31d0151f5716b655421d9d75b7fa74cc4e744df2"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.39.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[HDF5]]
deps = ["Blosc", "Compat", "HDF5_jll", "Libdl", "Mmap", "Random", "Requires"]
git-tree-sha1 = "83173193dc242ce4b037f0263a7cc45afb5a0b85"
uuid = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
version = "0.15.6"

[[HDF5_jll]]
deps = ["Artifacts", "JLLWrappers", "LibCURL_jll", "Libdl", "OpenSSL_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "fd83fa0bde42e01952757f01149dd968c06c4dba"
uuid = "0234f1f7-429e-5d53-9886-15a909be8d59"
version = "1.12.0+1"

[[HypertextLiteral]]
git-tree-sha1 = "72053798e1be56026b81d4e2682dbe58922e5ec9"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.0"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Lz4_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "5d494bc6e85c4c9b626ee0cab05daa4085486ab1"
uuid = "5ced341a-0733-55b8-9ab6-a4889d929147"
version = "1.9.3+0"

[[MAT]]
deps = ["BufferedStreams", "CodecZlib", "HDF5", "SparseArrays"]
git-tree-sha1 = "5c62992f3d46b8dce69bdd234279bb5a369db7d5"
uuid = "23992714-dd62-5051-b70f-ba57cb901cac"
version = "0.10.1"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OpenSSL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "15003dcb7d8db3c6c857fda14891a539a8f2705a"
uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
version = "1.1.10+0"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "9d8c00ef7a8d110787ff6f170579846f776133a9"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.4"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlutoUI]]
deps = ["Base64", "Dates", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "d1fb76655a95bf6ea4348d7197b22e889a4375f4"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.14"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "cc4bf3fdde8b7e3e9fa0351bdeedba1cf3b7f6e6"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.5.0+0"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─16a831a8-14b2-4d1e-86e7-a28f21d0dbd4
# ╟─716d10ed-38bd-4654-9a54-33dfba020591
# ╠═ad2d1e2b-4222-43d0-819a-4a7d14746327
# ╟─cef733a0-5541-4338-bd28-0916b697fb19
# ╟─679f15d4-e94d-43eb-8342-701111c46201
# ╠═0f683f10-71d7-477c-9b7c-49545e60fa8a
# ╠═095a2f1a-6ae3-451a-bc3f-7dad9a7834ea
# ╠═5d7c052a-73d7-487a-939e-f813561daff9
# ╠═b3e23596-0b67-4ed1-8caa-6353307ee315
# ╠═5f19c57c-e4ac-4c8c-9bc0-f6a015c8fb3f
# ╠═1062535d-3453-4a00-89f5-66060c150dbf
# ╟─5578b9b1-8fbd-4c22-a89d-cd39ea63fe87
# ╟─6629d6a5-af02-40c4-a7a3-a7146bf6638b
# ╟─baedaf07-7236-432c-a67b-bd7eaa99bb7c
# ╟─fd376fb1-1617-48a8-877b-31aa8f85476a
# ╠═45376f68-2552-11ec-32d7-a5de0d6b6882
# ╠═774dd760-cebf-44ff-b753-e54353c1317e
# ╠═debd1859-15ce-49ec-b45a-4707695f5c8a
# ╠═8904e3f1-88ad-424e-8ac2-e4c8dfd01c4a
# ╟─d1bd21e8-86c0-4607-953b-d33a7b964ecc
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
#

# Fallback
"""
    cache=alloc_mem(A,method)

Pre-allocates memory associated with matrix exponential function `method` and matrix `A`. To be used in combination with [`exponential!`](@ref).
"""
function alloc_mem(A,method)
    return nothing;
end

@deprecate _exp! exponential!
@deprecate exp_generic exponential!
exponential!(A)=exponential!(A,ExpMethodHigham2005(A));

## The diagonalization based
"""
    ExpMethodDiagonalization(enforce_real=true)

Matrix exponential method corresponding to the diagonalization with `eigen` possibly by removing imaginary part introduced by the numerical approximation.

"""
struct ExpMethodDiagonalization
    enforce_real::Bool
end
ExpMethodDiagonalization()=ExpMethodDiagonalization(true);

"""
    E=exponential!(A,[method [cache]])

Computes the matrix exponential with method specified in `method`. The contents of `A` is modified allowing for less allocations. The `method` parameter specifies the implementation and implementation parameters, e.g. [`ExpMethodNative`](@ref), [`ExpMethodDiagonalization`](@ref), [`ExpMethodGeneric`](@ref), [`ExpMethodHigham2005`](@ref). Memory
needed can be preallocated and provided in parameter `cache` such that the memory can recycled when calling `exponential!` several times. The preallocation is done with the command [`alloc_mem`](@ref): `cache=alloc_mem(A,method)`.

Example
```julia-repl
julia> A=randn(50,50);
julia> Acopy=B*2;
julia> method=ExpMethodHigham2005();
julia> cache=alloc_mem(A,method); # Main allocation done here
julia> E1=exponential!(A,method,cache) # Very little allocation here
julia> E2=exponential!(B,method,cache) # Very little allocation here
```

"""
function exponential!(A,method::ExpMethodDiagonalization,cache=nothing)
    F=eigen!(A)
    E=F.vectors*Diagonal(exp.(F.values))/F.vectors
    if (method.enforce_real && isreal(A))
        E=real.(E);
    end
    copyto!(A,E)
    return A
end

"""
    ExpMethodNative()

Matrix exponential method corresponding to calling `Base.exp`.

"""
struct ExpMethodNative
end
function exponential!(A,method::ExpMethodNative,cache=nothing)
    return exp(A)
end

max_stars_repo_name: JuliaDiffEq/ExponentialUtilities

max_issues_repo_name: JuliaDiffEq/ExponentialUtilities

max_forks_repo_name: JuliaDiffEq/ExponentialUtilities


text:
# name and module reflection
module_name(m::Module) = ccall(:jl_module_name, Any, (Any,), m)::Symbol
module_parent(m::Module) = ccall(:jl_module_parent, Any, (Any,), m)::Module
current_module() = ccall(:jl_get_current_module, Any, ())::Module

fullname(m::Module) = m===Main ? () : tuple(fullname(module_parent(m))...,
                                            module_name(m))

names(m::Module, all::Bool, imported::Bool) = ccall(:jl_module_names, Array{Symbol,1}, (Any,Int32,Int32), m, all, imported)
names(m::Module, all::Bool) = names(m, all, false)
names(m::Module) = names(m, false, false)
names(t::DataType) = collect(t.names)

function names(v)
    t = typeof(v)
    if isa(t,DataType)
        return names(t)
    else
        error("cannot call names() on a non-composite type")
    end
end

isconst(s::Symbol) =
    ccall(:jl_is_const, Int32, (Ptr{Void}, Any), C_NULL, s) != 0

isconst(m::Module, s::Symbol) =
    ccall(:jl_is_const, Int32, (Any, Any), m, s) != 0

# return an integer such that object_id(x)==object_id(y) if is(x,y)
object_id(x::ANY) = ccall(:jl_object_id, Uint, (Any,), x)

# type predicates
const isimmutable = x->(isa(x,Tuple) || !typeof(x).mutable)
isstructtype(t::DataType) = t.names!=() || (t.size==0 && !t.abstract)
isstructtype(x) = false
isbits(t::DataType) = !t.mutable & t.pointerfree & isleaftype(t)
isbits(t::Type) = false
isbits(x) = isbits(typeof(x))
isleaftype(t::ANY) = ccall(:jl_is_leaf_type, Int32, (Any,), t) != 0

typeintersect(a::ANY,b::ANY) = ccall(:jl_type_intersection, Any, (Any,Any), a, b)
typeseq(a::ANY,b::ANY) = a<:b && b<:a

function fieldoffsets(x::DataType)
    offsets = Array(Int, length(x.names))
    ccall(:jl_field_offsets, Void, (Any, Ptr{Int}), x, offsets)
    offsets
end

# subtypes
function _subtypes(m::Module, x::DataType, sts=Set(), visited=Set())
    push!(visited, m)
    for s in names(m,true)
        if isdefined(m,s)
            t = eval(m,s)
            if isa(t, DataType) && super(t).name == x.name
                push!(sts, t)
            elseif isa(t, Module) && !in(t, visited)
                _subtypes(t, x, sts, visited)
            end
        end
    end
    sts
end
subtypes(m::Module, x::DataType) = sort(collect(_subtypes(m, x)), by=string)
subtypes(x::DataType) = subtypes(Main, x)

subtypetree(x::DataType, level=-1) = (level == 0 ? (x, {}) : (x, {subtypetree(y, level-1) for y in subtypes(x)}))

# function reflection
isgeneric(f::ANY) = (isa(f,Function)||isa(f,DataType)) && isa(f.env,MethodTable)

function_name(f::Function) = isgeneric(f) ? f.env.name : (:anonymous)

code_lowered(f::Function,t::(Type...)) = map(m->uncompressed_ast(m.func.code), methods(f,t))
methods(f::ANY,t::ANY) = Any[m[3] for m in _methods(f,t,-1)]
_methods(f::ANY,t::ANY,lim) = _methods(f,{(t::Tuple)...},length(t::Tuple),lim,{})
function _methods(f::ANY,t::Array,i,lim::Integer,matching::Array{Any,1})
    if i == 0
        new = ccall(:jl_matching_methods, Any, (Any,Any,Int32), f, tuple(t...), lim)
        if new === false
            return false
        end
        append!(matching, new::Array{Any,1})
    else
        ti = t[i]
        if isa(ti, UnionType)
            for ty in (ti::UnionType).types
                t[i] = ty
                if _methods(f,t,i-1,lim,matching) === false
                    t[i] = ty
                    return false
                end
            end
            t[i] = ti
        else
            return _methods(f,t,i-1,lim,matching)
        end
    end
    matching
end

function methods(f::Function)
    if !isgeneric(f)
        error("not a generic function")
    end
    f.env
end

methods(t::DataType) = (_methods(t,Tuple,0);  # force constructor creation
                        t.env)

function length(mt::MethodTable)
    n = 0
    d = mt.defs
    while !is(d,())
        n += 1
        d = d.next
    end
    n
end

start(mt::MethodTable) = mt.defs
next(mt::MethodTable, m::Method) = (m,m.next)
done(mt::MethodTable, m::Method) = false
done(mt::MethodTable, i::()) = true

uncompressed_ast(l::LambdaStaticData) =
    isa(l.ast,Expr) ? l.ast : ccall(:jl_uncompress_ast, Any, (Any,Any), l, l.ast)

function _dump_function(f, t::ANY, native, wrapper)
    str = ccall(:jl_dump_function, Any, (Any,Any,Bool,Bool), f, t, native, wrapper)::ByteString
    if str == ""
        error("no method found for the specified argument types")
    end
    str
end

code_llvm  (f::Callable, types::(Type...)) = print(_dump_function(f, types, false, false))
code_native(f::Callable, types::(Type...)) = print(_dump_function(f, types, true, false))

function functionlocs(f::Callable, types=(Type...))
    locs = Any[]
    for m in methods(f, types)
        lsd = m.func.code::LambdaStaticData
        ln = lsd.line
        if ln > 0
            push!(locs, (find_source_file(string(lsd.file)), ln))
        end
    end
    if length(locs) == 0
       error("could not find function definition")
    end
    locs
end

functionloc(f::Callable, types=(Any...)) = functionlocs(f, types)[1]

function function_module(f::Function, types=(Any...))
    m = methods(f, types)
    if isempty(m)
        error("no matching methods")
    end
    m[1].func.code.module
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
module TestOneDimensionalRangeIterators

using Test
using MLJBase
using Random
import Distributions
using Statistics
using StableRNGs

rng = StableRNG(66600099)

const Dist = Distributions

mutable struct DummyModel <: Deterministic
    K::Int
    metric::Float64
    kernel::Char
end

dummy_model = DummyModel(4, 9.5, 'k')

mutable struct SuperModel <: Deterministic
    lambda::Float64
    model1::DummyModel
    model2::DummyModel
end

dummy1 = DummyModel(1, 9.5, 'k')
dummy2 = DummyModel(2, 9.5, 'k')
super_model = SuperModel(0.5, dummy1, dummy2)

z1 = range(dummy_model, :K, lower=1, upper=10)
z2 = range(dummy_model, :K, lower=10, origin=10^6, unit=10^5)
z3 = range(dummy_model, :K, upper=-10, origin=-10^6, unit=10^5)
z4 = range(super_model, :lambda, lower=1, upper=10)
z5 = range(dummy_model, :K, origin=10, unit=20)
p1 = range(dummy_model, :K, lower=1, upper=10, scale=:log10)
p2 = range(dummy_model, :kernel, values=['c', 'd'])
p3 = range(super_model, :lambda, lower=0.1, upper=1, scale=:log2)
p4 = range(dummy_model, :K, lower=1, upper=3, scale=x->2x)

[p4, p4]

@testset "scale transformations" begin
    @test transform(MLJBase.Scale, scale(:log), ℯ) == 1
    @test inverse_transform(MLJBase.Scale, scale(:log), 1) == float(ℯ)
end

@testset "inspecting scales of ranges" begin
    @test scale(p1) == :log10
    @test scale(p2) == :none
    @test scale(p3) == :log2
    @test scale(p4) == :custom
    @test scale(sin) === sin
end

@testset "iterators" begin

    @test iterator(p1, 5)  == [1, 2, 3, 6, 10]
    @test iterator(p2) == collect(p2.values)
    u = 2^(log2(0.1)/2)
    @test iterator(p3, 3) ≈ [0.1, u, 1]
    @test iterator(p4, 3) == [2, 4, 6]

    # semi-unbounded ranges:
    v = Int.(round.(exp.([(1-t)*log(10) + t*log(10+2e5)
                     for t in 0:(1/3):1]))) |> unique
    @test iterator(z2, 4) == v
    @test iterator(z3, 4) == reverse(-v)

    # doubly-unbounded ranges:
    @test iterator(z5, 4) ==
        iterator(range(Int, :foo, lower=-10, upper=30), 4)

    # truncated nominal range iteration:
    rr = range(Char, :foo, values=collect("abcdefg"))
    @test iterator(rr, 3) == ['a', 'b', 'c']

    # random:
    rng = StableRNG(66); @test iterator(rng, p1, 5) == [10, 2, 3, 6, 1]
    rng = StableRNG(22); @test iterator(rng, p1, 5) != [10, 2, 3, 6, 1]
    rng = StableRNG(33); @test iterator(rng, rr) == ['b', 'c', 'a', 'g',
                                    'f', 'd', 'e']
    rng = StableRNG(44); @test iterator(rng, rr) != ['b', 'c', 'a', 'g',
                                    'f', 'd', 'e']
    rng = StableRNG(88); @test iterator(rng, rr, 3) == ['b', 'c', 'a']
    rng = StableRNG(99); @test iterator(rng, rr, 3) != ['a', 'c', 'b']

    # with callable as scale:
    r = range(Int, :dummy, lower=1, upper=2, scale=x->10^x)
    expecting = map(x->round(Int,10^x), range(1, stop= 2, length=10))
    @test iterator(r, 10) == expecting
end

@testset "fitting distributions to NumericRange objects" begin

    # characterizations

    l = rand(rng)
    u = max(l, rand(rng)) + 1
    r = range(Int, :dummy, lower=l, upper=u)
    for D in [:Arcsine, :Uniform, :Biweight, :Cosine, :Epanechnikov,
              :SymTriangularDist, :Triweight]
        eval(quote
             d = Dist.fit(Dist.$D, $r)
             @test minimum(d) ≈ $l
             @test maximum(d) ≈ $u
             end
             )
    end

    o = randn(rng)
    s = rand(rng)
    r = range(Int, :dummy, lower=-Inf, upper=Inf, origin=o, unit=s)
    for D in [:Cauchy, :Gumbel, :Normal, :Laplace]
        eval(quote
             d = Dist.fit(Dist.$D, $r)
             @test Dist.location(d) ≈ $o
             @test Dist.scale(d) ≈ $s
             end
             )
    end

    o = rand(rng)
    s = o/(1 + rand(rng))
    r = range(Int, :dummy, lower=-Inf, upper=Inf, origin=o, unit=s)
    for D in [:Normal, :Gamma, :InverseGaussian, :LogNormal, :Logistic]
        eval(quote
             d = Dist.fit(Dist.$D, $r)
             @test mean(d) ≈ $o
             @test std(d) ≈ $s
             end
             )
    end

    r = range(Float64, :dummy, lower=-Inf, upper=Inf, unit=s, origin=o,)
    d = Dist.fit(Dist.Poisson, r)
    @test mean(d) ≈ s

    # truncation

    r = range(Int, :dummy, lower=l, upper=u)
    d = Dist.fit(Dist.Normal, r)
    @test minimum(d) == l
    @test maximum(d) == u

    # unsupported distributions

    @test_throws ArgumentError Dist.fit(Dist.Beta, r)

end

@testset "NumericSampler - distribution instance specified"  begin

    @testset  "integers" begin
        r = range(Int, :dummy, lower=11, upper=13)
        d = Dist.Uniform(1, 20)

        s = MLJBase.sampler(r, d)

        rng = StableRNG(0)
        dict = Dist.countmap(rand(rng, s, 1000))
        eleven, twelve, thirteen = map(x -> dict[x], 11:13)
        @test eleven == 252 && twelve == 514 && thirteen == 234

        rng = StableRNG(999);
        dict = Dist.countmap(rand(rng, s, 1000))
        eleven, twelve, thirteen = map(x -> dict[x], 11:13)
        @test eleven == 236 && twelve == 494 && thirteen == 270
    end

    @testset "right-unbounded floats" begin
        r = range(Float64, :dummy, lower=0.2, upper = Inf,
                  origin=5, unit=1) # origin and unit not relevant here
        s = MLJBase.sampler(r, Dist.Normal())

        Random.seed!(1);
        v = rand(s, 1000)
        @test all(x >= 0.2 for x in v)
        @test abs(minimum(v)/0.2 - 1) <= 0.01

        rng = Random.MersenneTwister(1);
        @test rand(rng, s, 1000) == v

        q = quantile(v, 0.0:0.1:1.0)
        Random.seed!(1);
        v2 = filter(x -> x>=0.2, rand(Dist.Normal(), 3000))[1:1000]
        q2 = quantile(v2, 0.0:0.1:1.0)
        @test all(x -> x≈1.0, q ./ q2)
    end

    @testset "sampler using callable scale" begin

        r = range(Int, :dummy, lower=1, upper=2, scale=x->10^x)
        s = sampler(r, Dist.Uniform)
        Random.seed!(123)
        v = rand(s, 100)
        @test issubset(v, 10:100)
        rng = MersenneTwister(123)
        @test rand(rng, s, 100) == v

        r = range(Float64, :dummy, lower=1, upper=2, scale=x->10^x)
        s = sampler(r, Dist.Uniform)
        Random.seed!(1)
        v = rand(s, 1000)
        @test abs(minimum(v) - 10) < 0.02
        @test abs(maximum(v) - 100) < 0.02
        rng = MersenneTwister(1)
        @test rand(rng, s, 1000) == v

    end

end

@testset "NumericSampler - distribution type specified"  begin

    r = range(Int, :k, lower=2, upper=6, origin=4.5, unit=1.2)
    s = MLJBase.sampler(r, Dist.Normal)
    v1 = rand(MersenneTwister(1), s, 50)
    d = Dist.truncated(Dist.Normal(r.origin, r.unit), r.lower, r.upper)
    v2 = map(x -> round(Int, x), rand(MersenneTwister(1), d, 50))
    @test v1 == v2

end

@testset "NominalSampler" begin

    r = range(Char, :(model.dummy), values=collect("cab"))

    @testset "probability vector specified" begin
        s = MLJBase.sampler(r, [0.1, 0.2, 0.7])
        rng = StableRNG(600)
        dict = Dist.countmap(rand(rng,s, 1000))
        c, a, b = map(x -> dict[x], collect("cab"))
        @test a == 201 && b == 714 && c == 85

        rng = StableRNG(89);
        dict = Dist.countmap(rand(rng, s, 1000))
        c, a, b = map(x -> dict[x], collect("cab"))
        @test a == 173 && b == 733 && c == 94
    end

    @testset "probability vector unspecified (uniform)" begin
        s = MLJBase.sampler(r)
        rng = StableRNG(55)
        dict = Dist.countmap(rand(rng,s, 1000))
        c, a, b = map(x -> dict[x], collect("cab"))
        @test a == 361 && b == 335 && c == 304

        rng = StableRNG(550)
        dict = Dist.countmap(rand(rng, s, 1000))
        c, a, b = map(x -> dict[x], collect("cab"))
        @test a == 332 && b == 356 && c == 312
    end

end

struct MySampler end
Base.rand(rng::AbstractRNG, ::MySampler) = rand(rng)

@testset "scale(s) for s a sampler" begin
    @test scale(MySampler()) == :none
    r = range(Char, :(model.dummy), values=collect("cab"))
    @test scale(MLJBase.sampler(r, [0.1, 0.2, 0.7])) == :none
    r1 = range(Int, :dummy, lower=1, upper=2, scale=x->10^x)
    @test scale(MLJBase.sampler(r1, Dist.Uniform)) == :custom
    r2 = range(Int, :k, lower=2, upper=6, origin=4.5, unit=1.2, scale=:log2)
    @test scale(MLJBase.sampler(r2, Dist.Normal)) == :log2
end

end

true


text:
function GaussianRandom(sys::System)
    return rand(sys.normal)
end

function RandomNumber()
    return rand()*2-1.0
end


text:
export KSP, SetOperators, GetOperators, SetFromOptions, KSPSolve,  SetUp, KSPSolveTranspose, KSP_NULL

export GetConvergedReason, PetscView, SetType, GetType, SetTolerances, GetTolerances, SetInitialGuessNonzero, GetInitialGuessNonzero, GetResidualNorm

"""
  KSP object
"""
type KSP
  pobj::Ptr{Void}

  """
    Constructor

    **Inputs**

     * comm: MPI_Comm

    **Outputs**

     * ksp object
  """
  function KSP(comm::MPI_Comm)  # constructor for KSP owned by the user
      ptr = Array{Ptr{Void}}(1)
      ierr = ccall((:KSPCreate,petsc),PetscErrorCode,(comm_type, Ptr{Void}),comm, ptr)
      @assert(ierr == 0)
      obj = new(ptr[1])
#      finalizer(obj, PetscDestroy)
      return obj
  end

  """
    KSP constructor from pointer to existing KSP object
  """
  function KSP(ptr::Ptr{Void})  # constructor for KSP *not* owned by the user
    return new(ptr)
  end
end

"""
  Null pointer KSP object
"""
global const KSP_NULL = KSP(C_NULL)

"""
  PetscDestroy for KSP object.  Safe to call multiple times.
"""
function PetscDestroy(ksp::KSP)

  if ksp.pobj != C_NULL
    err = ccall((:KSPDestroy,petsc),PetscErrorCode,(Ptr{Ptr{Void}},),&ksp.pobj)
  end

  ksp.pobj = C_NULL
end

"""
  PetscView

  **Inputs**

   * ksp: KSP object
   * viewer: PetscViewer, defaults to Petsc stdout
"""
function PetscView(ksp::KSP, viewer::PetscViewer=C_NULL)
    ccall((:KSPView,petsc),PetscErrorCode,(Ptr{Void}, Int64), ksp.pobj, viewer)
end

"""
  KSPSetOperators
"""
function SetOperators(ksp::KSP,Amat::PetscMat,Pmat::PetscMat)
   err = ccall((:KSPSetOperators,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void},Ptr{Void}), ksp.pobj, Amat.pobj, Pmat.pobj)
end

"""
  KSPGetOperators
"""
function GetOperators(ksp::KSP)
   Amat = Ref{Ptr{Void}}()
   Pmat = Ref{Ptr{Void}}()
   err = ccall((:KSPGetOperators,petsc),PetscErrorCode,(Ptr{Void},Ptr{Ptr{Void}},Ptr{Ptr{Void}}), ksp.pobj, Amat, Pmat)

   return PetscMat(Amat[]), PetscMat(Pmat[])
end

"""
  SetFromOptions
"""
function SetFromOptions(ksp::KSP)
    ccall((:KSPSetFromOptions,petsc),PetscErrorCode,(Ptr{Void},),ksp.pobj)
end

"""
  KSPSolve
"""
function KSPSolve(ksp::KSP, b::PetscVec, x::PetscVec)
    err = ccall((:KSPSolve,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void},Ptr{Void}), ksp.pobj, b.pobj, x.pobj)
end

"""
  KSPSolveTranspose
"""
function KSPSolveTranspose(ksp::KSP, b::PetscVec, x::PetscVec)
    err = ccall((:KSPSolveTranspose,petsc),PetscErrorCode,(Ptr{Void},Ptr{Void},Ptr{Void}), ksp.pobj, b.pobj, x.pobj)
end

"""
  KSPSetUp
"""
function SetUp(ksp::KSP)
    err = ccall((:KSPSetUp,petsc),PetscErrorCode,(Ptr{Void},), ksp.pobj)
end

"""
  KSPGetConvergedReason

  **Inputs**

   * ksp: KSP object

  **Outputs**

   * string containing reason
"""
function GetConvergedReason(ksp::KSP)
    reason = Array{KSPConvergedReason}(1)
    ccall((:KSPGetConvergedReason,petsc),PetscErrorCode,(Ptr{Void},Ptr{KSPConvergedReason}), ksp.pobj, reason)

    return reason[1]
end

"""
  KSPSetType
"""
function SetType(ksp::KSP, ksptype::KSPType)
    ccall((:KSPSetType,petsc),PetscErrorCode,(Ptr{Void},Cstring), ksp.pobj, ksptype)
end

"""
  KSPGetType

  **Inputs**

   * KSP

  **Outputs**

   * string containing KSP type
"""
function GetType(ksp::KSP)
    ksptype = Array{Ptr{UInt8}}(1)
    ccall((:KSPGetType,petsc),PetscErrorCode,(Ptr{Void}, Ptr{Ptr{UInt8}}), ksp.pobj, ksptype)

    return unsafe_string(ksptype[1])
end

"""
  KSPSetTolerances
"""
function SetTolerances(ksp::KSP,rtol::Number, abstol::Number, dtol::Number, maxits::Integer)

    ccall((:KSPSetTolerances,petsc),PetscErrorCode,(Ptr{Void}, PetscReal, PetscReal, PetscReal, PetscInt), ksp.pobj, rtol, abstol, dtol, maxits)
end

"""
  KSPGetTolerances

  **Inputs**

   * KSP

  **Outputs**

   * rtol
   * abstol
   * dtol
   * maxits
"""
function GetTolerances(ksp::KSP)
    rtol = Array{PetscReal}(1)
    abstol = Array{PetscReal}(1)
    dtol = Array{PetscReal}(1)
    maxits = Array{PetscInt}(1)

    ccall((:KSPGetTolerances,petsc),PetscErrorCode,(KSP,Ptr{PetscReal},Ptr{PetscReal},Ptr{PetscReal},Ptr{PetscInt}), ksp, rtol, abstol, dtol, maxits)

    return rtol[1], abstol[1], dtol[1], maxits[1]
end

"""
  KSPSetInitialGuessNonzero
"""
function SetInitialGuessNonzero(ksp::KSP, flg::PetscBool)
    ccall((:KSPSetInitialGuessNonzero,petsc),PetscErrorCode,(Ptr{Void},PetscBool), ksp.pobj, flg)
end

"""
  KSPGetInitialGuessNonzero

  **Inputs**

   * KSP

  **Outputs**

   * PetscBool
"""
function GetInitialGuessNonzero(ksp::KSP)
    flg_arr = Array{PetscBool}(1)
    ccall((:KSPGetInitialGuessNonzero,petsc),PetscErrorCode,(Ptr{Void}, Ptr{PetscBool}), ksp.pobj, flg_arr)
    return flg_arr[1]
end

"""
  KSPGetResidualNorm

  **Inputs**

   * KSP

  **Outputs**

   * PetscReal
"""
function GetResidualNorm(ksp::KSP)
   rnorm = Array{PetscReal}(1)
   ccall((:KSPGetResidualNorm,petsc),PetscErrorCode,(Ptr{Void}, Ptr{PetscReal}),ksp.pobj, rnorm)

    return rnorm[1]
end

### new function
# not tested

"""
  KSPSetReusePreconditioner
"""
function SetReusePreconditioner(arg1::KSP,arg2::PetscBool)
    ccall((:KSPSetReusePreconditioner,petsc),PetscErrorCode,(KSP,PetscBool),arg1,arg2)
end


text:
"""
Case 27:
This case study a three bus system with 1 machine (One d- One q-: 4th order model), a VSM of 19 states and an infinite source.
The test changes botht he voltage magnitude and phase angle of the source bus.
"""

##################################################
############### LOAD DATA ########################
##################################################

# Use the sme test data as Test 09
include(joinpath(TEST_FILES_DIR, "data_tests/test09.jl"))

##################################################
############### SOLVE PROBLEM ####################
##################################################

####### Changing magnitude of votlage at source bus #########

# time span
tspan = (0.0, 20.0);
case_source = collect(PSY.get_components(PSY.Source, threebus_sys))[1]

# Define Fault using Callbacks
V_source_change = SourceBusVoltageChange(1.0, case_source, :V_ref, 1.02)

@testset "Test 27 Source Bus Voltage Magnitude Perturbation ResidualModel" begin
    path = (joinpath(pwd(), "test-27"))
    !isdir(path) && mkdir(path)
    try
        # Define Simulation Problem
        sim = Simulation(
            ResidualModel,
            threebus_sys, # system
            path,
            tspan,
            V_source_change,
        )

        # Test Initial Condition
        diff_val = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test09_x0_init
            diff_val[1] += LinearAlgebra.norm(res[k] - v)
        end
        @test (diff_val[1] < 1e-3)

        # Solve problem
        execute!(sim, IDA(), dtmax = 0.02)
        results = read_results(sim)

        # Obtain data for angles
        series = get_state_series(results, ("generator-103-1", :θ_oc))
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

@testset "Test 27 Source Bus Voltage Magnitude Perturbation MassMatrixModel" begin
    path = (joinpath(pwd(), "test-27"))
    !isdir(path) && mkdir(path)
    try
        # Define Simulation Problem
        sim = Simulation(
            MassMatrixModel,
            threebus_sys, # system
            path,
            tspan,
            V_source_change,
        )

        # Test Initial Condition
        diff_val = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test09_x0_init
            diff_val[1] += LinearAlgebra.norm(res[k] - v)
        end
        @test (diff_val[1] < 1e-3)

        # Solve problem
        execute!(sim, Rodas4(), dtmax = 0.02)
        results = read_results(sim)

        # Obtain data for angles
        series = get_state_series(results, ("generator-103-1", :θ_oc))
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

####### Changing angle of voltage at source bus #########

#time span
tspan = (0.0, 20.0);
case_source = collect(PSY.get_components(PSY.Source, threebus_sys))[1]

#Define Fault using Callbacks
V_source_change = SourceBusVoltageChange(1.0, case_source, :θ_ref, 0.1)

@testset "Test 27 Source Bus Voltage Angle Perturbation ResidualModel" begin
    path = (joinpath(pwd(), "test-27"))
    !isdir(path) && mkdir(path)
    try
        # Define Simulation Problem
        sim = Simulation(
            ResidualModel,
            threebus_sys, # system
            path,
            tspan,
            V_source_change,
        )

        # Test Initial Condition
        diff_val = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test09_x0_init
            diff_val[1] += LinearAlgebra.norm(res[k] - v)
        end
        @test (diff_val[1] < 1e-3)

        # Solve problem
        execute!(sim, IDA(), dtmax = 0.02)
        results = read_results(sim)

        # Obtain data for angles
        series = get_state_series(results, ("generator-103-1", :θ_oc))
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

@testset "Test 27 Source Bus Voltage Angle Perturbation MassMatrixModel" begin
    path = (joinpath(pwd(), "test-27"))
    !isdir(path) && mkdir(path)
    try
        sim = Simulation(
            MassMatrixModel,
            threebus_sys, # system
            path,
            tspan,
            V_source_change,
        )

        # Test Initial Condition
        diff_val = [0.0]
        res = get_init_values_for_comparison(sim)
        for (k, v) in test09_x0_init
            diff_val[1] += LinearAlgebra.norm(res[k] - v)
        end
        @test (diff_val[1] < 1e-3)

        # Solve problem
        execute!(sim, Rodas4(), dtmax = 0.02)
        results = read_results(sim)

        # Obtain data for angles
        series = get_state_series(results, ("generator-103-1", :θ_oc))
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end


text:
# domain.jl

# Trait definition dispatcher
abstract type Dispatcher end
struct T0 <: Dispatcher end
struct T1 <: Dispatcher end
struct T2 <: Dispatcher end
struct T3 <: Dispatcher end

Dispatcher(::Type{<: String}) = T0()
Dispatcher(::Type{<: Array{AppliSales.Order, 1}}) = T1()
Dispatcher(::Type{<: Array{AppliGeneralLedger.JournalEntry,1}}) = T2()
Dispatcher(::Type{<: Array{AppliAR.BankStatement,1}}) = T3()

# dispatch(x::T) where {T} = dispatch(Dispatcher(T), x)
# end Trait defenition


text:
using Test, simgym
import simgym: leapfrog!
include("lineagetest.jl")

@test 1==1

testpos = [0.,0.]
testvel = [.0001,.0001]

cell = Cell[Cell(copy(testpos),copy(testvel),40,30)]
s = Simulation(cell,len_jones)

dt = .1
nsteps = 1

# this is what leapfrog SHOULD be doing
@. testpos += dt*testvel

leapfrog!(s,dt)

#@test testpos == s.cells[1].p

max_stars_repo_name: andi-bergeson/simgym

max_issues_repo_name: andi-bergeson/simgym

max_forks_repo_name: andi-bergeson/simgym


text:
# this allows a test of timing for an experiment without setting up any
# multimedia resources, so it can be run just about anywhere.
function find_timing(fn;keys...)
  empty!(Weber.null_record)
  exp = Experiment(;null_window=true,hide_output=true,keys...)
  setup(() -> fn(),exp)

  run(exp,await_input=false)

  nostarts = filter(x -> !endswith(string(x[:code]),"_start") &&
                    x[:code] != "terminated" &&
                    x[:code] != "closed",
                    Weber.null_record)
  map(x -> x[:code],nostarts),map(x -> x[:time],nostarts),nostarts
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct AdlsSecurityException <: SwaggerModel
    exception::Any # spec type: Union{ Nothing, String } # spec name: exception
    javaClassName::Any # spec type: Union{ Nothing, String } # spec name: javaClassName
    message::Any # spec type: Union{ Nothing, String } # spec name: message

    function AdlsSecurityException(;exception=nothing, javaClassName=nothing, message=nothing)
        o = new()
        validate_property(AdlsSecurityException, Symbol("exception"), exception)
        setfield!(o, Symbol("exception"), exception)
        validate_property(AdlsSecurityException, Symbol("javaClassName"), javaClassName)
        setfield!(o, Symbol("javaClassName"), javaClassName)
        validate_property(AdlsSecurityException, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        o
    end
end # type AdlsSecurityException

const _property_map_AdlsSecurityException = Dict{Symbol,Symbol}(Symbol("exception")=>Symbol("exception"), Symbol("javaClassName")=>Symbol("javaClassName"), Symbol("message")=>Symbol("message"))
const _property_types_AdlsSecurityException = Dict{Symbol,String}(Symbol("exception")=>"String", Symbol("javaClassName")=>"String", Symbol("message")=>"String")
Base.propertynames(::Type{ AdlsSecurityException }) = collect(keys(_property_map_AdlsSecurityException))
Swagger.property_type(::Type{ AdlsSecurityException }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_AdlsSecurityException[name]))}
Swagger.field_name(::Type{ AdlsSecurityException }, property_name::Symbol) =  _property_map_AdlsSecurityException[property_name]

function check_required(o::AdlsSecurityException)
    (getproperty(o, Symbol("exception")) === nothing) && (return false)
    true
end

function validate_property(::Type{ AdlsSecurityException }, name::Symbol, val)
end


text:
using Base.Test
using PGM

@testset "Graphs" begin

# @testset "Tree.Basics" begin
#
#     g1 = Tree(7, 2,
#               [3, 3, 4, 4, 4, 4, 6],  # parents
#               [1, 1, 2, 0, -1, -1, -2])  # parent link-types
#
#     @test nvertices(g1) == 7
#     @test nedgetypes(g1) == 2
#     @test nedges(g1) == 6
#
#     @test [parent(g1, v) for v in 1:7] == [3, 3, 4, 4, 4, 4, 6]
#     @test [parent_linktype(g1, v) for v in 1:7] == [1, 1, 2, 0, -1, -1, -2]
#     @test [children(g1, v) for v in 1:7] ==
#         [Int[], Int[], [1,2], [3,5,6], Int[], [7], Int[]]
#     @test [children_linktypes(g1, v) for v in 1:7] ==
#         [Int[], Int[], [-1,-1], [-2,1,1], Int[], [2], Int[]]
#     @test [degree(g1, v) for v in 1:7] == [1, 1, 3, 3, 1, 2, 1]
#
#     @test root(g1) == 4
#     @test [Int(isroot(g1, v)) for v in 1:7] == [0, 0, 0, 1, 0, 0, 0]
#     @test [Int(isleaf(g1, v)) for v in 1:7] == [1, 1, 0, 0, 1, 0, 1]
# end

@testset "UGraph.Basics" begin

    # shared_kind is false
    nv = 5
    g1 = UGraph(5, [(1, 2), (1, 3), (2, 3), (2, 4), (3, 4), (4, 5)])
    es1 = [Edge(1, 2, 1), Edge(1, 3, 2), Edge(2, 3, 3),
           Edge(2, 4, 4), Edge(3, 4, 5), Edge(4, 5, 6)]
    @test nvertices(g1) == 5
    @test nedgekinds(g1) == 6
    @test nedges(g1) == 6
    @test vertices(g1) == 1:5
    @test edges(g1) == es1
    @test [edges_of(g1, k) for k = 1:6] == [[e] for e in es1]
    @test [indegree_at(g1, v) for v=1:nv] == [0, 1, 2, 2, 1]
    @test [outdegree_at(g1, v) for v=1:nv] == [2, 2, 1, 1, 0]
    @test [degree_at(g1, v) for v=1:nv] == [2, 3, 3, 3, 1]
    @test [inedges_at(g1, v) for v=1:nv] ==
        [Edge[], [Edge(1, 2, 1)], [Edge(1, 3, 2), Edge(2, 3, 3)],
         [Edge(2, 4, 4), Edge(3, 4, 5)], [Edge(4, 5, 6)]]
    @test [outedges_at(g1, v) for v=1:nv] ==
        [[Edge(1, 2, 1), Edge(1, 3, 2)], [Edge(2, 3, 3), Edge(2, 4, 4)],
         [Edge(3, 4, 5)], [Edge(4, 5, 6)], Edge[]]

    # shared_kind is true
    nv = 5
    g2 = UGraph(5, [(1, 2), (1, 3), (2, 3), (2, 4), (3, 4), (4, 5)];
                shared_kind=true)
    es2 = [Edge(1, 2, 1), Edge(1, 3, 1), Edge(2, 3, 1),
           Edge(2, 4, 1), Edge(3, 4, 1), Edge(4, 5, 1)]
    @test nvertices(g2) == 5
    @test nedgekinds(g2) == 1
    @test nedges(g2) == 6
    @test vertices(g2) == 1:5
    @test edges(g2) == es2
    @test edges_of(g2, 1) == edges(g2)
    @test [indegree_at(g2, v) for v=1:nv] == [0, 1, 2, 2, 1]
    @test [outdegree_at(g2, v) for v=1:nv] == [2, 2, 1, 1, 0]
    @test [degree_at(g2, v) for v=1:nv] == [2, 3, 3, 3, 1]
    @test [inedges_at(g2, v) for v=1:nv] ==
     [Edge[], [Edge(1, 2, 1)], [Edge(1, 3, 1), Edge(2, 3, 1)],
      [Edge(2, 4, 1), Edge(3, 4, 1)], [Edge(4, 5, 1)]]
    @test [outedges_at(g2, v) for v=1:nv] ==
     [[Edge(1, 2, 1), Edge(1, 3, 1)], [Edge(2, 3, 1), Edge(2, 4, 1)],
      [Edge(3, 4, 1)], [Edge(4, 5, 1)], Edge[]]

    # grouped explicitly
    g3 = UGraph(5, [[(1, 2), (1, 3), (2, 4), (3, 4)],  # edges of 1st type
                    [(2, 3), (4, 5)]])                 # edges of 2nd type
    es3 = [Edge(1, 2, 1), Edge(1, 3, 1), Edge(2, 4, 1), Edge(3, 4, 1),
           Edge(2, 3, 2), Edge(4, 5, 2)]
    @test nvertices(g3) == 5
    @test nedgekinds(g3) == 2
    @test nedges(g3) == 6
    @test vertices(g3) == 1:5
    @test edges(g3) == es3
    @test edges_of(g3, 1) == es3[1:4]
    @test edges_of(g3, 2) == es3[5:6]
    @test [indegree_at(g3, v) for v=1:nv] == [0, 1, 2, 2, 1]
    @test [outdegree_at(g3, v) for v=1:nv] == [2, 2, 1, 1, 0]
    @test [degree_at(g3, v) for v=1:nv] == [2, 3, 3, 3, 1]
    @test [inedges_at(g3, v) for v=1:nv] ==
        [Edge[], [Edge(1, 2, 1)], [Edge(1, 3, 1), Edge(2, 3, 2)],
         [Edge(2, 4, 1), Edge(3, 4, 1)], [Edge(4, 5, 2)]]
    @test [outedges_at(g3, v) for v=1:nv] ==
        [[Edge(1, 2, 1), Edge(1, 3, 1)], [Edge(2, 4, 1), Edge(2, 3, 2)],
         [Edge(3, 4, 1)], [Edge(4, 5, 2)], Edge[]]
end

# @testset "UGraph.ToTree" begin
#
#     g1 = ugraph_with_tedges(7,
#         [[(1, 3), (2, 3), (4, 5), (4, 6)], [(3, 4), (6, 7)]])
#     @assert nvertices(g1) == 7
#     @assert nedgetypes(g1) == 2
#     @assert nedges(g1) == 6
#     @assert nedges(g1, 1) == 4
#     @assert nedges(g1, 2) == 2
#
#     t1a = to_tree(g1, 1)
#     @test nvertices(t1a) == 7
#     @test nedgetypes(t1a) == 2
#     @test nedges(t1a) == 6
#     @test [parent(t1a, v) for v in 1:7] == [1, 3, 1, 3, 4, 4, 6]
#     @test [parent_linktype(t1a, v) for v in 1:7] == [0, 1, -1, -2, -1, -1, -2]
#     @test [children(t1a, v) for v in 1:7] ==
#         [[3], Int[], [2,4], [5,6], Int[], [7], Int[]]
#     @test [children_linktypes(t1a, v) for v in 1:7] ==
#         [[1], Int[], [-1,2], [1,1], Int[], [2], Int[]]
#
#     t1b = to_tree(g1, 4)
#     @test nvertices(t1b) == 7
#     @test nedgetypes(t1b) == 2
#     @test nedges(t1b) == 6
#     @test [parent(t1b, v) for v in 1:7] == [3, 3, 4, 4, 4, 4, 6]
#     @test [parent_linktype(t1b, v) for v in 1:7] == [1, 1, 2, 0, -1, -1, -2]
#     @test [children(t1b, v) for v in 1:7] ==
#         [Int[], Int[], [1,2], [3,5,6], Int[], [7], Int[]]
#     @test [children_linktypes(t1b, v) for v in 1:7] ==
#         [Int[], Int[], [-1,-1], [-2,1,1], Int[], [2], Int[]]
#
#     t1c = to_tree(g1, 7)
#     @test nvertices(t1c) == 7
#     @test nedgetypes(t1c) == 2
#     @test nedges(t1c) == 6
#     @test [parent(t1c, v) for v in 1:7] == [3, 3, 4, 6, 4, 7, 7]
#     @test [parent_linktype(t1c, v) for v in 1:7] == [1, 1, 2, 1, -1, 2, 0]
#     @test [children(t1c, v) for v in 1:7] ==
#         [Int[], Int[], [1,2], [3,5], Int[], [4], [6]]
#     @test [children_linktypes(t1c, v) for v in 1:7] ==
#         [Int[], Int[], [-1,-1], [-2,1], Int[], [-1], [-2]]
# end

end # Graphs


text:
# functions to delete hybridization
# originally in functions.jl
# Claudia MArch 2015

# --------------------------------- delete hybridization -------------------------------

# function to identify inCycle (with priority queue)
# based on updateInCycle as it is the exact same code
# only without changing the incycle attribute
# only returning array of edges/nodes affected by the hybrid
# used when attempting to delete
# input: hybrid node around which we want to identify inCycle
# needs module "Base.Collections"
# returns tuple: nocycle, array of edges changed, array of nodes changed
# check: is this traversal much faster than a simple loop over
#        all edges/nodes and check if incycle==hybrid.number?
function identifyInCycle(net::Network,node::Node)
    node.hybrid || error("node $(node.number) is not hybrid, cannot identifyInCycle")
    start = node;
    hybedge = getHybridEdge(node);
    last = getOtherNode(hybedge,node);
    dist = 0;
    queue = PriorityQueue();
    path = Node[];
    net.edges_changed = Edge[];
    net.nodes_changed = Node[];
    push!(net.edges_changed,hybedge);
    push!(net.nodes_changed,node);
    found = false;
    net.visited = [false for i = 1:size(net.node,1)];
    enqueue!(queue,node,dist);
    while(!found)
        if(isempty(queue))
            return true, net.edges_changed, net.nodes_changed
        else
            curr = dequeue!(queue);
            if(isEqual(curr,last))
                found = true;
                push!(path,curr);
            else
                if(!net.visited[getIndex(curr,net)])
                    net.visited[getIndex(curr,net)] = true;
                    if(isEqual(curr,start))
                        for e in curr.edge
                            if(!e.hybrid || e.isMajor)
                                other = getOtherNode(e,curr);
                                other.prev = curr;
                                dist = dist+1;
                                enqueue!(queue,other,dist);
                            end
                        end
                    else
                        for e in curr.edge
                            if(!e.hybrid || e.isMajor)
                                other = getOtherNode(e,curr);
                                if(!other.leaf && !net.visited[getIndex(other,net)])
                                    other.prev = curr;
                                    dist = dist+1;
                                    enqueue!(queue,other,dist);
                                end
                            end
                        end
                    end
                end
            end
        end
    end # end while
    curr = pop!(path);
    while(!isEqual(curr, start))
        if(curr.inCycle == start.number)
            push!(net.nodes_changed, curr);
            edge = getConnectingEdge(curr,curr.prev);
            if(edge.inCycle == start.number)
                push!(net.edges_changed, edge);
            end
            curr = curr.prev;
        end
    end
    return false, net.edges_changed, net.nodes_changed
end

# aux function to traverse the network
# similar to traverseContainRoot but only
# identifying the edges that would be changed by a given
# hybridization
# warning: it does not go accross hybrid node/edge,
#          nor tree node with minor hybrid edge
function traverseIdentifyRoot(node::Node, edge::Edge, edges_changed::Array{Edge,1})
    if(!node.leaf && !node.hybrid)
        for e in node.edge
            if(!isEqual(edge,e) && e.isMajor && !e.hybrid)
                other = getOtherNode(e,node);
                push!(edges_changed, e);
                if(!other.hybrid)
                    #if(!other.hasHybEdge)
                        traverseIdentifyRoot(other,e, edges_changed);
                    #else
                    #    if(hybridEdges(other)[1].isMajor)
                    #        traverseIdentifyRoot(other,e, edges_changed);
                    #    end
                    #end
                end
            end
        end
    end
end

# function to identify containRoot
# depending on a hybrid node on the network
# input: hybrid node (can come from searchHybridNode)
# return array of edges affected by the hybrid node
function identifyContainRoot(net::HybridNetwork, node::Node)
    node.hybrid || error("node $(node.number) is not hybrid, cannot identify containRoot")
    net.edges_changed = Edge[];
    for e in node.edge
        if(!e.hybrid)
            other = getOtherNode(e,node);
            push!(net.edges_changed,e);
            traverseIdentifyRoot(other,e, net.edges_changed);
        end
    end
    return net.edges_changed
end

# function to undo the effect of a hybridization
# and then delete it
# input: network, hybrid node, random flag
#        random = true, deletes one hybrid egde at random
#        (minor with prob 1-gamma, major with prob gamma)
#        random = false, deletes the minor edge always
# warning: it uses the gamma of the hybrid edges even if
#          it is not identifiable like in bad diamond I (assumes undone by now)
# blacklist = true: add the edge as a bad choice to put a hybridization (not fully tested)
function deleteHybridizationUpdate!(net::HybridNetwork, hybrid::Node, random::Bool, blacklist::Bool)
    global DEBUG
    hybrid.hybrid || error("node $(hybrid.number) is not hybrid, so we cannot delete hybridization event around it")
    DEBUG && println("MOVE: delete hybridization on hybrid node $(hybrid.number)")
    nocycle, edgesInCycle, nodesInCycle = identifyInCycle(net,hybrid);
    !nocycle || error("the hybrid node $(hybrid.number) does not create a cycle")
    edgesRoot = identifyContainRoot(net,hybrid);
    edges = hybridEdges(hybrid);
    undoGammaz!(hybrid,net);
    othermaj = getOtherNode(edges[1],hybrid)
    edgesmaj = hybridEdges(othermaj)
    DEBUG && println("edgesmaj[3] $(edgesmaj[3].number) is the one to check if containRoot=false already: $(edgesmaj[3].containRoot)")
    if(edgesmaj[3].containRoot) #if containRoot=true, then we need to undo
        push!(edgesRoot, edges[1]) ## add hybrid edges to edgesRoot to undo containRoot
        push!(edgesRoot, edges[2])
        undoContainRoot!(edgesRoot);
    end
    if(DEBUG)
        edges[1].gamma >= 0.5 || println("strange major hybrid edge $(edges[1].number) with gamma $(edges[1].gamma) less than 0.5")
        edges[1].gamma != 1.0 || println("strange major hybrid edge $(edges[1].number) with gamma $(edges[1].gamma) equal to 1.0")
    end
    limit = edges[1].gamma
    if(random)
        minor = rand() < limit ? false : true
    else
        minor = true;
    end
    deleteHybrid!(hybrid,net,minor, blacklist)
    undoInCycle!(edgesInCycle, nodesInCycle); #moved after deleteHybrid to mantain who is incycle when deleteEdge and look for partition
    undoPartition!(net,hybrid, edgesInCycle)
end

deleteHybridizationUpdate!(net::HybridNetwork, hybrid::Node) = deleteHybridizationUpdate!(net, hybrid, true, false)

# function to delete a hybridization event
# input: hybrid node and network
#        minor: true (deletes minor edge), false (deletes major)
# warning: it is meant after undoing the effect of the
#          hybridization in deleteHybridizationUpdate!
#          by itself, it leaves things as is
# branch lengths of -1.0 are interpreted as missing.
function deleteHybrid!(node::Node,net::HybridNetwork,minor::Bool, blacklist::Bool)
    node.hybrid || error("node $(node.number) has to be hybrid for deleteHybrid")
    if(minor)
        hybedge1,hybedge2,treeedge1 = hybridEdges(node);
        other1 = getOtherNode(hybedge1,node);
        other2 = getOtherNode(hybedge2,node);
        other3 =  getOtherNode(treeedge1,node);
        if(hybedge1.number > treeedge1.number)
            setLength!(treeedge1, addBL(treeedge1.length, hybedge1.length));
            removeNode!(node,treeedge1);
            setNode!(treeedge1,other1);
            setEdge!(other1,treeedge1);
            removeEdge!(other1, hybedge1);
            deleteEdge!(net,hybedge1);
            #treeedge1.containRoot = (!treeedge1.containRoot || !hybedge1.containRoot) ? false : true #causes problems if hybrid.CR=false
            if(blacklist)
                println("put in blacklist edge $(treeedge1.number)")
                push!(net.blacklist, treeedge1.number)
            end
        else
            makeEdgeTree!(hybedge1,node)
            other1.hasHybEdge = false;
            setLength!(hybedge1, addBL(hybedge1.length, treeedge1.length));
            removeNode!(node,hybedge1);
            setNode!(hybedge1,other3);
            setEdge!(other3,hybedge1);
            removeEdge!(other3,treeedge1);
            deleteEdge!(net,treeedge1);
            hybedge1.containRoot = (!treeedge1.containRoot || !hybedge1.containRoot) ? false : true
            if(blacklist)
                println("put in blacklist edge $(hybedge1.number)")
                push!(net.blacklist, hybedge1.number)
            end
        end
        hybindex = getIndex(true,[e.hybrid for e in other2.edge]);
        if(hybindex == 1)
            treeedge1 = other2.edge[2];
            treeedge2 = other2.edge[3];
        elseif(hybindex == 2)
            treeedge1 = other2.edge[1];
            treeedge2 = other2.edge[3];
        elseif(hybindex == 3)
            treeedge1 = other2.edge[1];
            treeedge2 = other2.edge[2];
        else
            error("strange node has more than three edges")
        end
        treenode1 = getOtherNode(treeedge1,other2);
        treenode2 = getOtherNode(treeedge2,other2);
        if(abs(treeedge1.number) > abs(treeedge2.number))
            setLength!(treeedge2, addBL(treeedge2.length, treeedge1.length));
            removeNode!(other2,treeedge2);
            setNode!(treeedge2,treenode1);
            setEdge!(treenode1,treeedge2);
            removeEdge!(treenode1,treeedge1);
            deleteEdge!(net,treeedge1);
            treeedge2.containRoot = (!treeedge1.containRoot || !treeedge2.containRoot) ? false : true
            if(blacklist)
                println("put in blacklist edge $(treeedge2.number)")
                push!(net.blacklist, treeedge2.number)
            end
        else
            setLength!(treeedge1, addBL(treeedge2.length, treeedge1.length));
            removeNode!(other2,treeedge1);
            setNode!(treeedge1,treenode2);
            setEdge!(treenode2,treeedge1);
            removeEdge!(treenode2,treeedge2);
            deleteEdge!(net,treeedge2);
            treeedge1.containRoot = (!treeedge1.containRoot || !treeedge2.containRoot) ? false : true
            if(blacklist)
                println("put in blacklist edge $(treeedge1.number)")
                push!(net.blacklist, treeedge1.number)
            end
        end
        #removeHybrid!(net,node);
        deleteNode!(net,node);
        deleteNode!(net,other2);
        deleteEdge!(net,hybedge2);
    else
        hybedge1,hybedge2,treeedge1 = hybridEdges(node);
        other1 = getOtherNode(hybedge1,node);
        other2 = getOtherNode(hybedge2,node);
        setLength!(treeedge1, addBL(treeedge1.length, hybedge2.length))
        removeEdge!(other2,hybedge2)
        removeNode!(node,treeedge1)
        setEdge!(other2,treeedge1)
        setNode!(treeedge1,other2)
        #removeHybrid!(net,node)
        deleteNode!(net,node)
        deleteEdge!(net,hybedge1)
        deleteEdge!(net,hybedge2)
        removeEdge!(other1,hybedge1)
        size(other1.edge,1) == 2 || error("strange node $(other1.number) had 4 edges")
        if(abs(other1.edge[1].number) < abs(other1.edge[2].number))
            edge = other1.edge[1]
            otheredge = other1.edge[2]
        else
            edge = other1.edge[2]
            otheredge = other1.edge[1]
        end
        setLength!(other1.edge[1], addBL(other1.edge[1].length, other1.edge[2].length))
        other3 =  getOtherNode(otheredge,other1);
        removeNode!(other1,edge)
        removeEdge!(other3,otheredge)
        setEdge!(other3,edge)
        setNode!(edge,other3)
        deleteNode!(net,other1)
        deleteEdge!(net,otheredge)
    end
end

deleteHybrid!(node::Node,net::HybridNetwork,minor::Bool) = deleteHybrid!(node,net,minor, false)

"""
`deleteHybridEdge!(net::HybridNetwork,edge::Edge)`

Deletes a hybrid edge from a network. The network does not have to be of level 1,
and may contain some polytomies. Updates branch lengths, allowing for missing values.
Returns the network.

At each of the 2 junctions, the child edge is retained (i.e. the tree edge is retained,
below the hybrid node).

Warnings:

- does **not** update containRoot (could be implemented later)
- does **not** update attributes needed for snaq! (like containRoot, inCycle, edge.z, edge.y etc.)
- if the parent of edge is the root, the root will be moved to keep the network unrooted
  with a root of degree two.
"""
function deleteHybridEdge!(net::HybridNetwork,edge::Edge)
    edge.hybrid || error("edge $(edge.number) has to be hybrid for deleteHybridEdge!")
    n1 = (edge.isChild1 ? edge.node[1] : edge.node[2])  # child  of edge, to be deleted
    n1.hybrid || error("child node $(n1.number) of hybrid edge $(edge.number) should be a hybrid.")
    n2 = (edge.isChild1 ? edge.node[2] : edge.node[1])  # parent of edge, to be deleted too.
    # next: keep hybrid node n1 if it has 4+ edges (2 parents and 2+ children).
    #       2 or 1 edges should never occur.
    if length(n1.edge) < 3
        error("node $(n1.number) has $length(n1.edge) edges instead of 3+");
    elseif length(n1.edge) == 3
        pe = nothing # will be other parent (hybrid) edge of n1
        ce = nothing # will be child edge of n1, to be merged with pe
        for e in n1.edge
            if (e.hybrid && e!=edge) pe = e; end
            if !(e.hybrid)           ce = e; end
        end
        pn = getOtherNode(pe,n1); # parent node of n1, other than n2
        pn ≢ n2 || error("k=2 cycle: 2 hybrid edges from node $(n2.number) to node $(n1.number)")
        atRoot = (net.node[net.root] ≡ n1) # n1 should not be root, but if so, pn will be new root
        # next: replace ce by pe+ce, remove n1 and pe from network.
        ce.length = addBL(ce.length, pe.length)
        removeNode!(n1,ce) # ce now has 1 single node cn
        setNode!(ce,pn)    # ce now has 2 nodes in this order: cn, pn
        ce.isChild1 = true
        setEdge!(pn,ce)
        removeEdge!(pn,pe)
        # if (pe.number<ce.number) ce.number = pe.number; end # bad to match edges between networks
        deleteEdge!(net,pe,part=false) # decreases net.numEdges   by 1
        deleteNode!(net,n1) # decreases net.numHybrids by 1, numNodes too.
        # warning: containRoot could be updated in ce and down the tree.
        if (atRoot)
            try
                net.root = getIndex(pn,net)
            catch e
                if isa(e, ErrorException) error("node $(pn.number) not in net!"); end
            end
        end
    else # n1 has 4+ edges (polytomy): keep n2 but detach it from 'edge'
        removeEdge!(n1,edge) # does not update n1.hybrid at this time
        warn("polytomy at node $(n1.number). Assuming only 2 hybrid parents, though.")
        n1.hybrid = false
    end

    # next: keep n2 if it has 4+ edges. 2 or 1 edges should never occur.
    #       If root, would have no parent: treat network as unrooted and change the root.
    if length(n2.edge) < 2
        error("node $(n2.number) (parent of hybrid edge $(edge.number) to be deleted) has 1 edge only!")
    elseif length(n2.edge) == 2 # if n2=root or degree-2 node
        # remove n2 and both of its edges (including 'edge')
        pe = (edge ≡ n2.edge[1] ? n2.edge[2] : n2.edge[1]) #  <--edge-- n2 ---pe--- pn
        pn = (  n2 ≡ pe.node[1] ? pe.node[2] : pe.node[1])
        if net.node[net.root] ≡ n2 # if n2 was root, new root = pn
            net.root = getIndex(pn, net)
        end
        # remove n2 and pe
        removeEdge!(pn,pe)
        deleteEdge!(net,pe,part=false)
        deleteNode!(net,n2)
    elseif length(n2.edge) == 3
        oei = Int[] # n2's edges' indices, other than 'edge'.
        for i=1:length(n2.edge)
            if (n2.edge[i] != edge) push!(oei, i); end
        end
        length(oei)==2 || error("node $(n2.number) has 3 edges, but $(length(oei)) different from edge $(edge.number)")
        ce = n2.edge[oei[1]] # ce will be kept
        pe = n2.edge[oei[2]] # pe will be folded into new ce = pe+ce
        switch = false
        if getOtherNode(pe, n2).leaf
            !getOtherNode(ce, n2).leaf ||
                error("root $(n2.number) connected to 1 hybrid and 2 leaves: edges $(pe.number) and $(ce.number).")
            switch = true
        elseif (n2 ≡ ce.node[(ce.isChild1 ? 1 : 2)]) && (n2 ≡ pe.node[(pe.isChild1 ? 2 : 1)])
            switch = true
        end
        if switch # to ensure correct direction isChild1 for new edges if the original was up-to-date
            ce = n2.edge[oei[2]] # but no error if original isChild1 was outdated
            pe = n2.edge[oei[1]] # and to ensure that pn can be the new root if needed.
        end
        atRoot = (net.node[net.root] ≡ n2) # if n2=root, new root will be 'pn' = other node of pe
        # next: replace ce by pe+ce, remove n2 and pe from network.
        pn = getOtherNode(pe,n2) # parent node of n2 if n2 not root. Otherwise, pn will be new root.
        ce.length = addBL(ce.length, pe.length)
        removeNode!(n2,ce) # ce now has 1 single node cn
        setNode!(ce,pn)    # ce now has 2 nodes in this order: cn, pn
        ce.isChild1 = true
        setEdge!(pn,ce)
        removeEdge!(pn,pe)
        # if (pe.number<ce.number) ce.number = pe.number; end # bad to match edges between networks
        deleteEdge!(net,pe,part=false)
        deleteNode!(net,n2)
        if (atRoot)
            try
                net.root = getIndex(pn,net)
            catch e
                if isa(e, ErrorException) error("node $(pn.number) not in net!"); end
            end
        end
    else # n2 has 4+ edges (polytomy): keep n2 but detach it from 'edge'
        removeEdge!(n2,edge)
    end
    # finally: remove hybrid 'edge' from network
    deleteEdge!(net,edge,part=false)
    return net
end

# function to update net.partition after deleting a hybrid node
# needs a list of the edges in cycle
function undoPartition!(net::HybridNetwork, hybrid::Node, edgesInCycle::Vector{Edge})
    global DEBUG
    hybrid.hybrid || error("node $(hybrid.number) is not hybrid, and we need hybrid node inside deleteHybUpdate for undoPartition")
    if(net.numHybrids == 0)
        net.partition = Partition[]
    else
        cycles = Int[]
        edges = Edge[]
        N = length(net.partition)
        i = 1
        while(i <= N)
            DEBUG && println("hybrid number is $(hybrid.number) and partition is $([e.number for e in net.partition[i].edges]), with cycle $(net.partition[i].cycle)")
            if(in(hybrid.number,net.partition[i].cycle))
                DEBUG && println("hybrid number matches with partition.cycle")
                p = splice!(net.partition,i)
                DEBUG && println("after splice, p partition has edges $([e.number for e in p.edges]) and cycle $(p.cycle)")
                ind = getIndex(hybrid.number,p.cycle)
                deleteat!(p.cycle,ind) #get rid of that hybrid number
                cycles = vcat(cycles,p.cycle)
                edges = vcat(edges,p.edges)
                DEBUG && println("edges is $([e.number for e in edges]) and cycles is $(cycles)")
                N = length(net.partition)
            else
                i += 1
            end
        end
        for e in edgesInCycle
            DEBUG && println("edge in cycle is $(e.number)")
            if(isEdgeNumIn(e,net.edge)) #only include edge if still in net
                DEBUG && println("edge is in net still")
                push!(edges,e)
            end
        end
        newPartition = Partition(unique(cycles),edges)
        DEBUG && println("new partition with cycle $(newPartition.cycle), edges $([e.number for e in newPartition.edges])")
        push!(net.partition,newPartition)
    end
end


text:
"""
$(TYPEDEF)

Introduce direct observations on all dimensions of a Pose2 variable:

Example:
--------
```julia
PriorPose2( MvNormal([10; 10; pi/6.0], Matrix(Diagonal([0.1;0.1;0.05].^2))) )
```
"""
# mutable struct PriorPose2{T} <: IncrementalInference.AbstractPrior  where {T <: IncrementalInference.SamplableBelief}
#   Z::T
#   # empty constructor
#   PriorPose2{T}() where T = new{T}()
#   # regular constructor
#   PriorPose2{T}(x::T) where {T <: IncrementalInference.SamplableBelief}  = new{T}(x)
# end
# # convenience and default object helper
# PriorPose2(x::T) where {T <: IncrementalInference.SamplableBelief} = PriorPose2{T}(x)
struct PriorPose2{T <: SamplableBelief, P} <: IIF.AbstractPrior
  Z::T
  p::P
end

PriorPose2(z::IIF.SamplableBelief) = PriorPose2(z, getPointIdentity(Pose2))
PriorPose2{T,P}(z::IIF.SamplableBelief) where {T,P} = PriorPose2{T,P}(z, getPointIdentity(Pose2))

DFG.getManifold(::PriorPose2) = SpecialEuclidean(2)

function getSample(cf::CalcFactor{<:PriorPose2}, N::Int=1)
  Z = cf.factor.Z
  p = cf.factor.p
  M = getManifold(cf.factor)

  Xc = rand(Z)

  # X = get_vector.(Ref(M), Ref(p), Xc, Ref(DefaultOrthogonalBasis()))
  X = hat(M, p, Xc)
  points = exp(M, p, X)

  return points
end

function (cf::CalcFactor{<:PriorPose2})(m, p)
  #
  M = getManifold(cf.factor)
  return log(M, p, m)
  # iXihat = SE2(meas) \ SE2(wXi)
  # return se2vee(iXihat)
end

#TODO Serialization of reference point p
## Serialization support

"""
$(TYPEDEF)
"""
mutable struct PackedPriorPose2  <: IncrementalInference.PackedInferenceType
    str::String
    PackedPriorPose2() = new()
    PackedPriorPose2(x::String) = new(x)
end
function convert(::Type{PackedPriorPose2}, d::PriorPose2)
  return PackedPriorPose2(convert(PackedSamplableBelief, d.Z))
end
function convert(::Type{PriorPose2}, d::PackedPriorPose2)
  distr = convert(SamplableBelief, d.str)
  return PriorPose2(distr)
end

## NOTE likely deprecated comparitors, see DFG compareFields, compareAll instead
function compare(a::PriorPose2,b::PriorPose2; tol::Float64=1e-10)
  compareDensity(a.Z, b.Z)
end


text:
function do_all(filename, dir; Ny=32, Nz=32, yLength=2, zLength=2,
                               ν=1e-2, timestepper = :RungeKutta3, advection = WENO5(),
                               z_position=0.6, width=0.3, amplitude=-1e-1,
                               Δt=0.002, stop_time=40.0, iteration_interval=1000,
                               schedule_interval=0.5, output_fields = :all)

    grid = make_grid(Ny, Nz, yLength, zLength)
    v_bc, w_bc = no_slip_bc(grid)
    model = make_model(grid, ν)
    set_buoyancy!(model, z_position=z_position,
                         width=width,
                         amplitude=amplitude)
    sim = make_simulation(model, dir, filename, Δt=Δt,
                                                stop_time=stop_time,
                                                iteration_interval=iteration_interval,
                                                schedule_interval=schedule_interval,
                                                output_fields=output_fields)
    run!(sim)
    return nothing
end