text:
const PYJLGCCACHE = Dict{PyPtr,Any}()

@kwdef struct PyJuliaValueObject
    ob_base::PyObject = PyObject()
    value::Ptr{Cvoid} = C_NULL
    weaklist::PyPtr = C_NULL
end

pyjlbase_new(t::PyPtr, ::PyPtr, ::PyPtr) = begin
    o = ccall(UnsafePtr{PyTypeObject}(t).alloc[!], PyPtr, (PyPtr, Py_ssize_t), t, 0)
    if !isnull(o)
        PyJuliaValue_SetValue(o, nothing)
        UnsafePtr{PyJuliaValueObject}(o).weaklist[] = C_NULL
    end
    o
end

pyjlbase_init(o::PyPtr, args::PyPtr, kwargs::PyPtr) = begin
    ism1(PyArg_CheckNumArgsEq("__init__", args, 1)) && return Cint(-1)
    ism1(PyArg_CheckNoKwargs("__init__", kwargs)) && return Cint(-1)
    ism1(PyArg_GetArg(Any, "__init__", args, 0)) && return Cint(-1)
    PyJuliaValue_SetValue(o, takeresult(Any))
    Cint(0)
end

pyjlbase_dealloc(o::PyPtr) = begin
    delete!(PYJLGCCACHE, o)
    isnull(UnsafePtr{PyJuliaValueObject}(o).weaklist[!]) || PyObject_ClearWeakRefs(o)
    ccall(UnsafePtr{PyTypeObject}(Py_Type(o)).free[!], Cvoid, (PyPtr,), o)
    nothing
end

const PyJuliaBaseValue_Type__ref = Ref(PyNULL)
PyJuliaBaseValue_Type() = begin
    ptr = PyJuliaBaseValue_Type__ref[]
    if isnull(ptr)
        c = []
        t = fill(
            PyType_Create(
                c,
                name = "juliaaa.ValueBase",
                basicsize = sizeof(PyJuliaValueObject),
                new = pyglobal(:PyType_GenericNew),
                init = @cfunction(pyjlbase_init, Cint, (PyPtr, PyPtr, PyPtr)),
                dealloc = @cfunction(pyjlbase_dealloc, Cvoid, (PyPtr,)),
                flags = Py_TPFLAGS_BASETYPE |
                        Py_TPFLAGS_HAVE_VERSION_TAG |
                        (CONFIG.isstackless ? Py_TPFLAGS_HAVE_STACKLESS_EXTENSION : 0x00),
                weaklistoffset = fieldoffset(PyJuliaValueObject, 3),
                getattro = pyglobal(:PyObject_GenericGetAttr),
                setattro = pyglobal(:PyObject_GenericSetAttr),
            ),
        )
        ptr = PyPtr(pointer(t))
        err = PyType_Ready(ptr)
        ism1(err) && return PyNULL
        PYJLGCCACHE[ptr] = push!(c, t)
        PyJuliaBaseValue_Type__ref[] = ptr
    end
    ptr
end

PyJuliaValue_Check(o) = Py_TypeCheck(o, PyJuliaBaseValue_Type())

PyJuliaValue_GetValue(__o) = begin
    _o = Base.cconvert(PyPtr, __o)
    GC.@preserve _o begin
        o = Base.unsafe_convert(PyPtr, _o)
        p = UnsafePtr{PyJuliaValueObject}(o)
        Base.unsafe_pointer_to_objref(p.value[!])
    end
end

PyJuliaValue_SetValue(__o, v) = begin
    _o = Base.cconvert(PyPtr, __o)
    GC.@preserve _o begin
        o = Base.unsafe_convert(PyPtr, _o)
        p = UnsafePtr{PyJuliaValueObject}(o)
        p.value[!], PYJLGCCACHE[o] = Python.pointer_from_obj(v)
    end
end

PyJuliaValue_New(t, v) = begin
    if isnull(t)
        if !PyErr_IsSet()
            PyErr_SetString(PyExc_Exception(), "Got NULL type with no error set")
        end
        return PyNULL
    end
    bt = PyJuliaBaseValue_Type()
    isnull(bt) && return PyNULL
    PyType_IsSubtype(t, bt) != 0 || (
        PyErr_SetString(PyExc_TypeError(), "Expecting a subtype of 'juliaaa.ValueBase'"); return PyNULL
    )
    o = _PyObject_New(t)
    isnull(o) && return PyNULL
    UnsafePtr{PyJuliaValueObject}(o).weaklist[] = C_NULL
    PyJuliaValue_SetValue(o, v)
    o
end

PyJuliaBaseValue_New(v) = begin
    t = PyJuliaBaseValue_Type()
    isnull(t) && return PyNULL
    PyJuliaValue_New(t, v)
end

PyJuliaValue_TryConvert_any(o, ::Type{S}) where {S} = begin
    x = PyJuliaValue_GetValue(o)
    putresult(tryconvert(S, x))
end


text:
import CurveFit
using Printf
using stoosc
rc("font", family="arial")

# Maximal considered input intensity
I_max = 0.3

## Deterministic Arnold tongue
filename = joinpath("arnold_tongues", "model01", "ode.csv")
arnold_tongue = load_data(filename)
entrainment_area_ode = estimate_entrainment_area(arnold_tongue, I_max)
println("ode entrainment area = $(entrainment_area_ode)")
fig2, ax2 = subplots(figsize=(2.5, 2.25))
plot_arnold_tongue(arnold_tongue, I_max, ax=ax2)
ax2.set_title("Deterministic")
fig2.tight_layout()
filename = joinpath("arnold_tongues", "model01", "arnold_deterministic.svg")
save_figure(fig2, filename)

## Arnold tongue for a stochastic model with many oscillators
filename = joinpath("arnold_tongues", "model01", "sde_n=1000_σ=0.005.csv")
arnold_tongue = load_data(filename)
entrainment_area = estimate_entrainment_area(arnold_tongue, I_max)
println("entrainment area = $(entrainment_area)")
fig2, ax2 = subplots(figsize=(2.5, 2.25))
plot_arnold_tongue(arnold_tongue, I_max, ax=ax2)
ax2.set_title("n = 1000, σ = 0.005")
fig2.tight_layout()
filename = joinpath("arnold_tongues", "model01", "arnold_many_oscillators.svg")
save_figure(fig2, filename)

## Arnold tongue for a stochastic model with few oscillators
filename = joinpath("arnold_tongues", "model01", "sde_n=10_σ=0.005.csv")
arnold_tongue = load_data(filename)
entrainment_area = estimate_entrainment_area(arnold_tongue, I_max)
println("entrainment area = $(entrainment_area)")
fig2, ax2 = subplots(figsize=(2.5, 2.25))
plot_arnold_tongue(arnold_tongue, I_max, ax=ax2)
ax2.set_title("n = 10, σ = 0.005")
fig2.tight_layout()
filename = joinpath("arnold_tongues", "model01", "arnold_few_oscillators.svg")
save_figure(fig2, filename)

## Increasing noise for constant population size
fig, ax = subplots(figsize=(7, 3))
ax.plot(log10.([0.001, 0.03]), [entrainment_area_ode, entrainment_area_ode], "--",
    color="black", label="Deterministic")

# Save the values for the maximal entrainment area
maximum_entrainment_area = Dict(
    "n" => Float64[],
    "σ" => Float64[],
    "entrainment_area" => Float64[]
)

x = [0.001, 0.002, 0.003, 0.005, 0.01, 0.02, 0.03]  # noise intensitiess
for n = [1000, 100, 10, 1]  # number of oscillators
    # for each population size

    y = Float64[]  # entrainment area

    for σ = x
        # for each noise value

        # Load the arnold tongue for given n and σ
        filename = joinpath("arnold_tongues", "model01", "sde_n=$(n)_σ=$(σ).csv")
        arnold_tongue = load_data(filename)

        # Save arnold tongues for the individual points in the figure
        fig2, ax2 = subplots(figsize=(2.5, 2.25))
        plot_arnold_tongue(arnold_tongue, I_max, ax=ax2)
        ax2.set_title("n = $(n); σ = $(σ)")
        fig2.tight_layout()
        filename = joinpath("arnold_tongues", "model01", "sde_n_$(n)_σ_$(σ).svg")
        save_figure(fig2, filename)

        # Calculate and save the entrainment area
        entrainment_area = estimate_entrainment_area(arnold_tongue, I_max)
        push!(y, entrainment_area)

    end

    # Plot entrainment area versus noise intensity
    x_log = log10.(x)
    f_fit = CurveFit.curve_fit(CurveFit.Polynomial, x_log, y, 2)
    h = ax.plot(x_log, y, "o", label="n = $n")

    # Fit a polynomial for a smooth curve interpolation and plot it
    x_new = range(x_log[1], x_log[end], length=100)
    y_new = f_fit.(x_new)
    ax.plot(x_new, y_new, "-", color=h[1].get_color())

    # Find point of maximal entrainment area and save it into a dictionary
    index = argmax(y_new)
    push!(maximum_entrainment_area["σ"], x_new[index])
    push!(maximum_entrainment_area["n"], n)
    push!(maximum_entrainment_area["entrainment_area"], y_new[index])
    println("Maximum for n = $n is EA = $(y_new[index]) at σ ≈ $(10^x_new[index])")

end

# Set labels for the figure
ax.set_xlabel("log\$_{10}\$(Noise intensity)")
ax.set_ylabel("Entrainment area")
ax.set_title("Entrainment area for constant number of oscillators (n)")
ax.legend(edgecolor="black", framealpha=1.0, loc="best", ncol=1)
fig.tight_layout()

# Save the main figure
filename = joinpath("arnold_tongues", "model01", "overview_n.svg")
save_figure(fig, filename)

## Plot maximal entrainment area and corresponding noise intenisty for each n
fig, ax_arr = subplots(1, 2, figsize=(7, 3))
n = log10.(maximum_entrainment_area["n"])
σ = maximum_entrainment_area["σ"]
entrainment_area = maximum_entrainment_area["entrainment_area"]

ax_arr[1].plot(n, σ, "o", color="black")
ax_arr[1].set_xlabel("log\$_{10}\$(Number of oscillators, n)")
ax_arr[1].set_ylabel("log\$_{10}\$(Noise intensity, σ)")
ax_arr[1].set_title("Noise intensity for maximal entrainment area")

ax_arr[2].plot(n, entrainment_area, "o", color="black")
ax_arr[2].set_xlabel("log\$_{10}\$(Number of oscillators, n)")
ax_arr[2].set_ylabel("Entrainment area")
ax_arr[2].set_title("Maximal entrainment area")

fig.tight_layout()
save_figure(fig, "arnold_tongues", "model01", "maximal_entrainment_n.svg")

## Increasing number of oscillators for constant system volume
fig, ax = subplots(figsize=(7, 3))
ax.plot(log10.([1, 1000]), [entrainment_area_ode, entrainment_area_ode], "--",
    color="black", label="Deterministic")

# Save the values for the maximal entrainment area
maximum_entrainment_area = Dict(
    "V" => Float64[],
    "n" => Float64[],
    "entrainment_area" => Float64[]
)

n_arr = [1, 4, 10, 40, 100, 400, 1_000]
for i = [10, 1, 0.4, 0.2, 0.1]
    Ω_arr = convert.(Int, [1_000_000, 250_000, 100_000, 25_000, 10_000, 2_500, 1_000] .* i)
    y = Float64[]
    x = n_arr
    for (n, Ω) = zip(n_arr, Ω_arr)

        filename = joinpath("arnold_tongues", "model01", "sde_n=$(n)_Ω=$(Ω).csv")
        arnold_tongue = load_data(filename)

        # Save arnold tongues for the individual points in the figure
        fig2, ax2 = subplots(figsize=(2.5, 2.25))
        plot_arnold_tongue(arnold_tongue, I_max, ax=ax2)
        ax2.set_title("n = $(n); Ω = $(Ω)")
        fig2.tight_layout()
        filename = joinpath("arnold_tongues", "model01", "sde_n_$(n)_Ω_$(Ω).svg")
        save_figure(fig2, filename)

        entrainment_area = estimate_entrainment_area(arnold_tongue, I_max)
        push!(y, entrainment_area)

    end
    x_log = log10.(x)
    f_fit = CurveFit.curve_fit(CurveFit.Polynomial, x_log, y, 3)
    x_new = range(x_log[1], x_log[end], length=100)
    y_new = f_fit.(x_new)

    index = argmax(y_new)
    push!(maximum_entrainment_area["V"], Ω_arr[1])
    push!(maximum_entrainment_area["n"], x_new[index])
    push!(maximum_entrainment_area["entrainment_area"], y_new[index])
    println("Maximum for i = $i is EA = $(y_new[index]) at n ≈ $(10^x_new[index])")

    plot_label = @sprintf("V = %.0e", Ω_arr[1])
    h = ax.plot(x_log, y, "o", label=plot_label)
    ax.plot(x_new, y_new, "-", color=h[1].get_color())

end

# Set labels for the figure
ax.set_xlabel("Number of oscillators")
ax.set_xticks(log10.(n_arr))
ax.set_xticklabels(n_arr)
ax.set_ylabel("Entrainment area")
ax.set_title("Entrainment area for constant total volume (V)")
ax.legend(edgecolor="black", framealpha=1.0, loc="lower left", ncol=3)
fig.tight_layout()

filename = joinpath("arnold_tongues", "overview_V.svg")
save_figure(fig, filename)

## Plot maximal entrainment area and corresponding noise intenisty for each n
fig, ax_arr = subplots(1, 2, figsize=(7, 3))
V = log10.(maximum_entrainment_area["V"])[2:end]
n = maximum_entrainment_area["n"][2:end]
entrainment_area = maximum_entrainment_area["entrainment_area"][2:end]

ax_arr[1].plot(V, n, "o", color="black")
ax_arr[1].set_xlabel("log\$_{10}\$(Total volume, V)")
ax_arr[1].set_ylabel("log\$_{10}\$(Number of oscillators, n)")
ax_arr[1].set_title("Total volume for maximal entrainment area")

ax_arr[2].plot(V, entrainment_area, "o", color="black")
ax_arr[2].set_xlabel("log\$_{10}\$(Total volume, V)")
ax_arr[2].set_ylabel("Entrainment area")
ax_arr[2].set_title("Maximal entrainment area")

fig.tight_layout()
save_figure(fig, "arnold_tongues", "model01", "maximal_entrainment_V.svg")

ext: jl

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
using WeightedArrays
using Test

R = Weighted(rand(3,10), rand(10))

@test size(R) == (3,10)
@test lastlength(R) == 10

@test sum(R.weights) ≈ 1 # normalised in construction

R.weights .+= 1
@test sum(weights(R)) ≈ 1 # normalised on reading

normalise!(R)
@test sum(R.weights) ≈ 1 # restored

@test tanh(R).opt.aname == "tanh-θ"
@test tanh(R).opt.lo == -1

S = sobol(2,9)

@test log(S).opt.clamp == false

@test size(unique(S + S)) == size(S) # hcat & unique

@test (2S).weights[1] ≈ 2/9
@test (2 .* S).array ≈ 2 .* (S.array)

X = xgrid(2,0:0.1:1)
X = near(X, S, 0.2)
@test size(X) == (2,78) # near

Y = wrandnp(Float32, 2,5)
@test eltype(Y.array) == eltype(weights(Y)) == Float32

Z = wrandn(big, 2,5)
@test eltype(Z.array) == eltype(weights(Z)) == BigFloat

pca = wPCA(R)
@test size(pca(sobol(3,10))) == (2,10) # pca can be applied

sPCA(R)
@test size(rPCA(sobol(3,20))) == (2,20) # was saved

using SliceMap

@test mapslices(x -> 2 .* x, S).array ≈ 2 .* S.array

@test mapcols(x -> 2 .* x, S).array ≈ 2 .* S.array
@test MapCols(x -> 2 .* x, S).array ≈ 2 .* S.array

mktempdir() do path
    w1 = Weighted([1 2; 3 4], [0.5, 1.0], WeightOpt(aname="a", wname="w"))

    f1csv = joinpath(path, "w1.csv")
    save(w1, f1csv)
    w1csv = load(f1csv)

    @test w1csv.array == w1.array
    @test w1csv.weights == weights(w1) # normalised
    @test_broken w1csv.opt.aname == "w1" # not saved

    f1json = joinpath(path, "w1.json")
    save(w1, f1json)
    w1json = load(f1json)

    @test w1json.array == w1.array
    @test w1json.weights == w1.weights
    @test w1json.opt.aname == "a"
end


text:
using Random

const rng = MersenneTwister()

abstract type Profiler end

Base.empty!(p::Profiler) = nothing
log_dispatch(p::Profiler, fn, type) = nothing

typefreqs(prof::Profiler) = values(prof.typefreqs)

function profileexpr(calledfn, argname)
    return quote
        Catwalk.log_dispatch(Catwalk.profiler(jitctx.callctxs.$calledfn), $(calledfn), typeof($argname))
    end
end

struct NoProfiler <: Profiler end

struct FullProfiler <: Profiler
    typefreqs::DataTypeFrequencies
    FullProfiler() = new(DataTypeFrequencies())
end

Base.empty!(p::FullProfiler) = empty!(p.typefreqs)

@inline function log_dispatch(p::FullProfiler, fn, type)
    increment!(p.typefreqs, type)
end

abstract type ProfileStrategy end
struct SparseProfile <: ProfileStrategy
    sparsity::Float16
    profiler::FullProfiler
    SparseProfile(sparsity = 0.01) = new(sparsity, FullProfiler())
end

function getprofiler(strategy::SparseProfile, round::Integer)
    if round <= 2 || rand(rng) < strategy.sparsity # First round is exploration only, second always profiled
        return strategy.profiler
    end
    return NoProfiler()
end


text:
function multi_binomial(a, b)
    # take two vectors a and b as arguments and calculate
    # a product of binomials computed from the corresponding vector components
    # C^a_b = binomial(a₁,b₁) * binomial(a₂,b₂) * ... * binomial(a_n, β_n)
    if length(a) != length(b)
        throw(ErrorException("vector lengths are different for some reason"))
    else
        C = prod([binomial(a[i], b[i]) for i in 1:length(a)])
        return C
    end
end

function derivative_matching(sys::MomentEquations, binary_vars::Array{Int,1}=Int[])

    closure_exp = OrderedDict()
    closure = OrderedDict() # additional dict to hold not expanded symbolic expressions

    N = sys.N

    if isempty(binary_vars)
        isbernoulli = false
    else
        sys = bernoulli_moment_eqs(sys, binary_vars)
        isbernoulli = true
    end

    # construct the raw moments up to mth order from the solved-for central moments
    if typeof(sys) == CentralMomentEquations
        closed_μ = central_to_raw_moments(N, sys.m_order)
        μ = central_to_raw_moments(N, sys.q_order)
    else
        μ = sys.μ
        closed_μ = copy(μ)
    end
    μ_symbolic = define_μ(sys.N, sys.q_order)

    # note that derivative matching is originally constructed by truncating at order of m_order+1
    # so if q_order > m_order + 1, we have to consider m_order+1, m_order+2, and so on in sequence
    # to build up the truncated raw moment expressions

    iter_k = vcat(sys.iter_1, sys.iter_m)

    for order in sys.m_order+1:sys.q_order

        # iterator through all moments of lower order
        iter_k = vcat(iter_k, filter(x -> sum(x) == order-1, sys.iter_q))

        length_k = length(iter_k)

        # iterator through all moments of the current truncation order
        iter_order = filter(x -> sum(x) == order, sys.iter_q)

        # initialise matrix needed for linear system defined below
        A = Matrix{Float64}(undef, length_k, length_k)

        for m in iter_order

            # each higher order moment μ_m is a product of lower order moments ∏ᵣ₌₁ᵏ(μ_mᵣ)^γᵣ
            # the coefficients γᵣ (contained in vector γ) are obtained by solving a system
            # of k equations (k - no. of all lower order moments):
            # C^(m)_(mᵤ) = ∑ᵣ₌₁ᵏ γᵣ * C^(mᵣ)_(mᵤ), s = 1, ..., k
            # where C^(a)_(b) is a product of binomials, so that
            # C^(a)_(b) = binomial(a₁, b₁)*binomial(a₂, b₂)*...*binomial(a_n, b_n)
            # hence for each higher order moment we set up a linear system A*γ=b and solve it

            b = [multi_binomial(m, iter_k[i]) for i in 1:length_k]
            for i in 1:length_k
                for j in 1:length_k
                    A[i, j] = multi_binomial(iter_k[j], iter_k[i])
                end
            end

            # The solution γ appears to always be integer-valued (also mentioned in the paper).
            # While there is no proof that it holds generall, we assume so in order
            # to use integer powers instead of float ones for simpler symbolic manipulations.
            # The code will break otherwise, throwing an InexactError attempting to convert Float to Int.
            # If this ever happens we know what to fix. TODO: add exception handling?
            γ = A\b

            closed_μ[m] = prod([closed_μ[iter_k[i]]^Int(γ[i]) for i in 1:length_k])
            closed_μ[m] = simplify(closed_μ[m])
            closure[μ_symbolic[m]] = prod([μ[iter_k[i]]^Int(γ[i]) for i in 1:length_k])
            closure[μ_symbolic[m]] = simplify(closure[μ_symbolic[m]])
        end

    end

    if typeof(sys) == CentralMomentEquations
        # construct the corresponding truncated expressions of higher order
        # central moments from the obtained raw moment expressions
        raw_to_central = raw_to_central_moments(N, sys.q_order, closed_μ, bernoulli=isbernoulli)
        central_to_raw = central_to_raw_moments(N, sys.q_order)
        closure_M = OrderedDict()
        for i in sys.iter_q
            closure_exp[sys.M[i]] = raw_to_central[i]
            expr = simplify(central_to_raw[i]-sys.M[i])
            closure_M[sys.M[i]] = simplify(closure[μ_symbolic[i]]-expr)
        end
        closure = closure_M
    else
        for i in sys.iter_q
            closure_exp[sys.μ[i]] = closed_μ[i]
        end
    end

    close_eqs(sys, closure_exp, closure, false)

end


text:
include("generator_utils.jl")

# repl_emoji_symbols_170
# repl_emoji_symbols_below_170
include("repl_emoji_symbols.jl")

emoji_data     = get_emoji_data()
dict_170       = get_diff_dict(repl_emoji_symbols_170,       emoji_data)
dict_below_170 = get_diff_dict(repl_emoji_symbols_below_170, emoji_data)

write_dict_to_file("generated_emoji_symbols_170.jl",       "emoji_symbols", dict_170)
write_dict_to_file("generated_emoji_symbols_below_170.jl", "emoji_symbols", dict_below_170)


text:
int(x) = (Int ∘ round)(x)
fint(x) = (Int ∘ floor)(x)
cint(x) = (Int ∘ ceil)(x)

"""
  or(x, y)

Return `y` if `x` is `Nothing`
"""
or(x, y) = isnothing(x) ? y : x

"""
  loop_last(v)

  Returns an iterable yielding tuples (is_last, value).
"""
function loop_last(v)
  is_last =  1:length(v) .== length(v)
    return zip(is_last, v)
end

function loop_firstlast(v)
  is_first =  1:length(v) .== 1
  is_last =  1:length(v) .== length(v)
  return zip(is_first, is_last, v)
end

"""
  get_lr_widths(width::Int)

To split something with `width` in 2, get the lengths
of the left/right widths.

When width is even that's easy, when it's odd we need to
be careful.
"""
function get_lr_widths(width::Int)::Tuple{Int, Int}
  iseven(width) && return (int(width/2), int(width/2))
  return (fint(width/2), cint(width/2))
end

"""
Get a clean string representation of an expression
"""
function expr2string(e::Expr)::String
  return replace_multi(
    string(e),
    '\n'=>"", ' '=>"", r"#=.*=#"=>"", "begin"=>"", "end"=>"")
end

"""
  get_file_format(nbytes; suffix="B")

Return a string with formatted file size.
"""
function get_file_format(nbytes; suffix="B")
  for unit in ["", "K", "M", "G", "T", "P", "E", "Z", "Y"]
      if nbytes < 1024.0
          _nd = round(nbytes; digits=2)
          return string(_nd) * " $unit$suffix"
      end
      nbytes = nbytes / 1024.0
  end
  return "cacca"
end


text:
using NumericalMethodsforEngineers, Test

a = [10. 1. -5.; -20. 3. 20.; 5. 3. 5.]
b = [1., 2., 6.]
n = size(A, 1)
x = 0.0

# Convert to upper triang;e
function to_upper(a::Matrix, b::Vector)
  at = deepcopy(a)
  bt = deepcopy(b)
  r, c = size(at)

  for k in 1:r-1
    if abs(at[k, k]) > 1.0e-6
    for i in k+1:r
      x = at[i, k] / at[k, k]
      at[i, k] = 0.0
      for j in k+1:r
      at[i, j] -= at[k, j] * x
      end
      bt[i] -= bt[k] * x
    end
    else
    md"Zero pivot found in row $k"
    end
  end
  (at, bt)
end

function direct_solve(a::Matrix, b::Vector)
  at = deepcopy(a)
  bt = deepcopy(b)
  r, c = size(at)

  for i in r:-1:1
    x = bt[i]
    if i < r
    for j in i+1:r
      x -= at[i, j] * bt[j]
    end
    end
    bt[i] = x / at[i, i]
  end
  bt
end

println("Modified lhs matrix and rhs")
at, bt = to_upper(a, b)
at |> display
println()
bt |> display
println()

println("Solution vector:")
c_nmfe = direct_solve(at, bt)
c_nmfe |> display
println()

c = a \ [1., 2., 6.]
c |> display

println()
@test ≈(c, c_nmfe, atol=1e-10)
@test round.(a * c; digits=14) == b


text:
struct SanchezLacombeParam <: EoSParam
    Mw::SingleParam{Float64}
    segment::SingleParam{Float64}
    epsilon::PairParam{Float64}
    vol::PairParam{Float64}
end

abstract type SanchezLacombeModel <: LatticeFluidModel end
include("mixing/mixing.jl")
struct SanchezLacombe{T <: SLMixingRule,I<:IdealModel} <:SanchezLacombeModel
    components::Array{String,1}
    icomponents::UnitRange{Int}
    mixing::T
    params::SanchezLacombeParam
    idealmodel::I
    references::Array{String,1}
end
@registermodel SanchezLacombe
const SL = SanchezLacombe

function SanchezLacombe(components;
    idealmodel=BasicIdeal,
    mixing = SLk0k1lMixingRule,
    userlocations=String[],
    ideal_userlocations=String[],
    mixing_userlocations = String[],
    verbose=false)
    params = getparams(components, ["LatticeFluid/SanchezLacombe","properties/molarmass.csv"]; userlocations=userlocations, verbose=verbose)

    segment = params["segment"]
    unmixed_epsilon = params["epsilon"]
    unmixed_vol = params["vol"]
    unmixed_epsilon.values #.*= k_B #to convert from temperature to eps
    unmixed_vol.values .*= 1e-6 #convert from cm3/mol to m3/mol
    Mw = params["Mw"]
    mixmodel = init_model(mixing,components,mixing_userlocations,verbose)
    ideal = init_model(idealmodel,components,ideal_userlocations,verbose)
    premixed_vol,premixed_epsilon = sl_mix(unmixed_vol,unmixed_epsilon,mixmodel)
    packagedparams = SanchezLacombeParam(Mw, segment, premixed_epsilon, premixed_vol)
    references = ["10.1016/S0378-3812(02)00176-0"]
    icomponents = 1:length(components)
    model = SanchezLacombe(components,icomponents,mixmodel,packagedparams,ideal,references)
    return model
end

include("mixing/SLk0k1lrule.jl")

function a_res(model::SanchezLacombe,V,T,z=SA[1.0])
    Σz = sum(z)
    r = model.params.segment.values
    mixing = model.mixing
    r̄ = dot(z,r)
    r̄ = r̄/Σz
    v_r,ε_r = mix_vε(model,V,T,z,mixing,r̄,Σz)
    v = V/Σz
    ρ̃ = r̄*v_r/v
    T̃ = R̄*T/ε_r
    _1 = one(V+T+first(z))
    return r̄*(-ρ̃ /T̃ + (_1/ρ̃  - _1)*log(1-ρ̃ )+_1)
end

function lb_volume(model::SanchezLacombe,z=SA[1.0])
    Σz = sum(z)
    r = model.params.segment.values
    v = model.params.vol.diagvalues
    r̄ = dot(z,r)
    #v_r,ε_r = mix_vε(model,0.0,0.0,z,model.mixing,r̄,Σz)
    return sum(r[i]*z[i]*v[i] for i in @comps)
end

function T_scale(model::SanchezLacombe,z=SA[1.0])
    Σz = sum(z)
    r = model.params.segment.values
    r̄ = dot(z,r)
    v_r,ε_r = mix_vε(model,0.0,0.0,z,model.mixing,r̄,Σz)
    return ε_r/R̄
end

function p_scale(model::SanchezLacombe,z=SA[1.0])
    Σz = sum(z)
    r = model.params.segment.values
    r̄ = dot(z,r)
    v_r,ε_r = mix_vε(model,0.0,0.0,z,model.mixing,r̄,Σz)
    Ts = ε_r/R̄
    vs = v_r*r̄
    return R̄*Ts/vs
end

function x0_volume_liquid(model::SanchezLacombe,T,z)
    v_lb = lb_volume(model,z)
    return v_lb*1.1
end
#SL does not work with the virial coefficient
function x0_volume_gas(model::SanchezLacombe,p,T,z)
    return sum(z)*R̄*T/p
end

function x0_sat_pure(model::SanchezLacombe,T,z=SA[1.0])
    Σz = sum(z)
    r = model.params.segment.values
    r̄ = dot(z,r)
    v_r,ε_r = mix_vε(model,0.0,T,z,model.mixing,r̄,Σz)
    Ts = ε_r/R̄
    vs = v_r*r̄
    Ps = R̄*Ts/vs
    Tr = T/Ts
    nan = zero(Tr)/zero(Tr)
    Tr > 1 && return [nan,nan]
    Tstar = Tr*PropaneRef_consts.T_c
    rhov = _propaneref_rhovsat(Tstar)
    vv = 1/rhov
    psat = pressure(model,vv,T)
    #vv = R̄*T/Ps
    vl = volume(model,psat,T,phase =:l)
    if isnan(vl)
        vv = nan
    end
    return [log10(vl),log10(vv)]
end

export SL,SanchezLacombe,SLk0k1lMixingRule


text:
# coding: utf-8

# In[1]:

#its also called levenshtein distance
#it can be done via recursion too
#the recursion version is much more elegant yet less efficient
# https://github.com/je-suis-tm/recursion/blob/master/edit%20distance%20recursion.jl

#edit distance is to minimize the steps transforming one string to another
#the way to solve this problem is very similar is to knapsack
#assume we have two strings text1 and text2
#we build a matrix with the size of (length(text1)+1)*(length(text2)+1)

#there are three different ways to transform a string
#insert,delete and replace
#we can use any of them or combined
#lets take a look at the best case first
#assume string text1 is string text2
#we dont need to do anything
#so 0 steps would be the answer
#for the worst case
#when string text1 has nothing in common with string text2
#we would have to replace the whole string text1
#the steps become the maximum step which is max(length(text1),length(text2))
#for general case,the number of steps would fall between the worst and the best case
#assume we are at i th letter of string text1 and j th letter string text2
#if we wanna get the optimal steps of transforming string text1 to string text2
#we have to make sure at each letter transformation
#text1[1:i] and text2[1:j] have reached their optimal status
#otherwise,we could always find another combination of insert,delete and replace
#to get a "real" optimal text1[1:i] and text2[1:j]
#it would make our string transformation not so optimal any more
#it is the same logic as the optimization of knapsack problem
#after we set our logic straight
#we would take a look at three different approaches
#lets take a look at insertion
#basically we need to insert j th letter from string text2 into string text1 at i th position
#the cumulated steps we have taken should be matrix[i][j-1]+1
#matrix[i][j-1] is the steps for text1[1:i] to text2[1:j]
#for delete,it is vice versa
#the cumulated steps we have taken should be matrix[i-1][j]+1
#for replacement,it is a lil bit tricky
#there are two scenarios
#if text1[i-1]==text2[j-1]
#it should be matrix[i-1][j-1]
#we dont need any replacement at all
#else,it should be matrix[i-1][j-1]+1
#we replace i th letter of string text1 with j th letter of string text2
#after we managed to understand three different approaches
#we want to take the minimum number of steps among these three approaches
#throughout the iteration of different positions of both strings
#in the end,we would get the optimal steps to transform one string to another,YAY

# In[2]:

function edit_distance(text1,text2)

    len1=length(text1)+1
    len2=length(text2)+1

    #this part is to create a matrix of (length(text1)+1)*(length(text2)+1)
    matrix=[[0 for _ in 1:len2] for _ in 1:len1]

    for i in 1:len1

        matrix[i][1]=i-1

    end

    for i in 1:len2

        matrix[1][i]=i-1

    end

    #we take iterations on both string text1 and text2
    #next,we check if text1[i-1]==text2[j-1]
    #if yes,no replacement needed
    #if no,replacement needed
    #we take a minimum function to see which combination would give the minimum steps
    #eventually we got what we are after
    for i in 2:len1

        for j in 2:len2

            if text1[i-1]==text2[j-1]

                matrix[i][j]=min(matrix[i-1][j]+1,
                matrix[i][j-1]+1,
                matrix[i-1][j-1])

            else

                matrix[i][j]=min(matrix[i-1][j]+1,
                matrix[i][j-1]+1,
                matrix[i-1][j-1]+1)

            end

        end

    end

    return matrix[len1][len2]

end

# In[3]:

println(edit_distance("baiseé","bas"))

max_stars_repo_path: edit distance dynamic programming.jl

max_issues_repo_path: edit distance dynamic programming.jl

max_forks_repo_path: edit distance dynamic programming.jl


text:
struct WaveReg <: AbstractRegister{1, ComplexF64}
    nbits::Int
    nt::Int
    wave::Matrix{Float64}
end

function WaveReg(nbits::Int, nt::Int)
    WaveReg(nbits, zeros(Float64, nt, nbits))
end

function instruct!(wr::WaveReg, loc::Int)
    t = t0(wr, loc)
    wavelet!(view(wr, t+1:t+wavelet_length(wr), loc))
end

wavelet_length(::Val{:X}) = 100
function wavelet!(v::AbstractVector, ::Val{:X})
end


text:
using ControlSystems;
using DifferentialEquations;
using PyPlot;
using Polynomials;

A=[0 1;1 0];
B=[0;1];
C=[0 0];
D=0;

sys=ss(A,B,C,D,0);

u(x,t)=[-2.11803 -2.11803 ]*x;
x0=[1;1];
tspan=0:0.01:20;
y, t, x, uout = lsim(sys,u,tspan,x0=x0);

plot(t,x[:,1],label="x1")
plot(t,x[:,2],label="x2")
xlabel("time (t)")
legend()
grid(true)

max_stars_repo_name: akshayhiregoudar/Control-System-Design

max_issues_repo_name: akshayhiregoudar/Control-System-Design

max_forks_repo_name: akshayhiregoudar/Control-System-Design


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

using Markdown
using InteractiveUtils

# ╔═╡ 81d44ebc-396d-11ec-2c9d-f16ce59cf2f3
md"""
# Customized application configuration with initializers

Initializers are plain Julia files which are loaded early in the application cycle (before routes, controller, or models). They are designed to expose configuration code which might be needed later on in the application life cycle.

Initializers should be placed within the `config/initializers/` folder and they will be automatically loaded (included) by Genie.

**If your configuration is environment dependent (ie a database connection which is different between dev and prod), it should be added to the corresponding `config/env/*.jl` file.**

## Best practices

* You can name the initializers as you wish (ideally a descriptive name, like maybe `redis.jl`).
* Don't use uppercase names unless you define a module (in order to respect Julia's naming practices).
* Keep your initializer files small and focused, so they serve only one purpose.
* You can add as many initializers as you need.
* Do not abuse them, they are not meant to host complex code.

## Load order

The initializers are loaded in the order they are read from the file system. If you have initializers which depend on other initializers, this is most likely a sign that you need to refactor using a model or a library file.

---
**HEADS UP**

Library files are Julia files which provide distinct functionality and can be placed in the `lib/` folder where they are also automatically loaded by Genie. If the `lib/` folder does not exist, you can create it yourself.

---

## Scope

All the definitions (variables, constants, functions, modules, etc) added to initializer files are loaded into your app's module. So if your app is called `MyGenieApp`, the definitions will be available under the `MyGenieApp` module.

---
**HEADS UP**

Given that your app's name is variable, you can also access your app's module through the `UserApp` constant. So all the definitions added to initializers can also be accessed through the `UserApp` module (`UserApp === MyGenieApp`).

---

## Example

This is Genie's default initializer for loading the SearchLight ORM. Please notice that the initializers does contain database configuration information, which is defined in a dedicated, environment dependent file. The initializer simply delegates the configuration loading and setup to the SearchLight package."""

# ╔═╡ c7dbff04-629f-4466-ab8a-7e4906cbbf6e
md"""
```julia
using SearchLight, SearchLight.QueryBuilder

SearchLight.Configuration.load()

if SearchLight.config.db_config_settings["adapter"] !== nothing
  SearchLight.Database.setup_adapter()
  SearchLight.Database.connect()
  SearchLight.load_resources()
end
```
"""

# ╔═╡ 22a66247-e1dc-43d7-8230-d7159bbf3d81
md"""
For more examples of initializers check [MVC App](https://github.com/essenciary/genie-watch-tonight/tree/main/config/initializers)
"""

# ╔═╡ Cell order:
# ╟─81d44ebc-396d-11ec-2c9d-f16ce59cf2f3
# ╟─c7dbff04-629f-4466-ab8a-7e4906cbbf6e
# ╟─22a66247-e1dc-43d7-8230-d7159bbf3d81


text:
using DrWatson
@quickactivate "INN_Velocity-Migration"

using JLD, Statistics, PyPlot, ImageFiltering, InvertibleNetworks, LinearAlgebra
include(srcdir("RTM_64x200.jl"))
using .RTM_64x200

X_orig = load(datadir("vel_4k_samples_64x200_lin_vel.jld"), "m_all")
X_orig = Float32.(X_orig)
n1, n2, nc, nsamples = size(X_orig)
ntrain = Int(nsamples*.9)
X_train_orig=X_orig[:, :, :, 1:ntrain]

AN = ActNorm(ntrain)
X_train_orig = AN.forward(X_train_orig) # zero mean and unit std of the training data X

# define data and figure path
datapath = plotsdir("chint/test3_6000_16_16")
figfolder = "rtm/test3_6000_16_16/posterior_samples2"
mkpath(plotsdir(figfolder))

X_fixed = load(joinpath(datapath, "posterior_samples2.jld"), "X_fixed")
X_post  = load(joinpath(datapath, "posterior_samples2.jld"), "X_post")
Y_fixed = load(joinpath(datapath, "posterior_samples2.jld"), "Y_fixed")
X_fixed = AN.inverse(X_fixed) # unnormarlize based on the same parameter of the training data X
X_post = AN.inverse(X_post)   # unnormarlize based on the same parameter of the training data X
X_post_mean = mean(X_post; dims=4)
Y_fixed = AN.inverse(Y_fixed) # unnormarlize based on the same parameter of the training data X

m = X_fixed[:, :, 1, 1] # true velocity in squared slowness
rtm0 = Y_fixed[:, :, 1, 1] # rtm images used in training INN
d = (25., 25.)
idx_wb = 9 # index of the water bottom

figure(); plot_rtm(adjoint(rtm0), d; new_fig=false, aspect = "auto"); title("RTM used in INN training")
savefig(plotsdir(figfolder, "rtm_INN.png"))

####################################################################################################
## Scenario 1: migration velocity m0 = all posterior samples
println(string("Scenario 1: migration velocity m0 = posterior samples"))

test_size = size(X_post)[4] #100
m0_all = zeros(Float32, n1, n2, test_size)
rtm_all = zeros(Float32, n1, n2, test_size)

for i = 1:3 #test_size

    m0 = X_post[:, :, 1, i]
    m0[1:idx_wb, :] .= m[1,1]
    m0[idx_wb+1:end, :] = imfilter(m0[idx_wb+1:end, :], Kernel.gaussian(3)) # smoothed velocity
    rtm = rtm_isic(m, m0)

    m0_all[:,:,i] = m0
    rtm_all[:,:,i] = adjoint(reshape(rtm, (n2, n1) ))

    println(string("RTM using the ", i, "th sample of X posterior finished"))
end

save(plotsdir(figfolder, "rtm_post_samples.jld"), "m0_all", m0_all, "rtm_all", rtm_all)

####################################################################################################
## Scenario 2: migration velocity m0 = posterior mean
println(string("Scenario 2: migration velocity m0 = posterior mean"))

m0 = X_post_mean[:, :, 1, 1]
m0[1:idx_wb, :] .= m[1,1]
m0[idx_wb+1:end, :] = imfilter(m0[idx_wb+1:end, :], Kernel.gaussian(3)) # smoothed velocity

rtm2 = rtm_isic(m, m0)
rtm2 = adjoint(reshape(rtm2, (n2, n1) ))
save(plotsdir(figfolder, "rtm_post_mean.jld"), "m", m, "m0", m0, "rtm", rtm2)

####################################################################################################
## Scenario 3: migration velocity m0 = migration velocity used in rtm
println(string("Scenario 3: migration velocity m0 = migration velocity used in rtm"))

# build a migration velocity m0 whose value increases linearly with depth
v0 = ones(n1, 1)
v0[idx_wb+1:end] = collect(range(3.5, 5.5, length = n1-idx_wb))
v0 = (1f0 ./ v0).^2 # convert to slowness squared
v0[1:idx_wb] .= m[1,1] # water layer velocity
m_lin = repeat(v0, 1, n2)

rtm3 = rtm_isic(m, m_lin)
rtm3 = adjoint(reshape(rtm3, (n2, n1) ))
save(plotsdir(figfolder, "rtm_lin_vel.jld"), "m", m, "m0", m_lin, "rtm", rtm3)

####################################################################################################
## Scenario 4: migration velocity m0 = smoothed true velocity
println(string("Scenario 4: migration velocity m0 = smoothed true velocity"))

m_smooth = m[1,1] .* ones(Float32, size(m)) # smoothed velocity
m_smooth[idx_wb+1:end, :] = imfilter(m[idx_wb+1:end, :], Kernel.gaussian(3)) # smoothed velocity

rtm4 = rtm_isic(m, m_smooth)
rtm4 = adjoint(reshape(rtm4, (n2, n1) ))
save(plotsdir(figfolder, "rtm_smooth_vel.jld"), "m", m, "m0", m_smooth, "rtm", rtm4)

####################################################################################################
## Scenario 5: migration velocity m0 = true velocity
println(string("Scenario 5: migration velocity m0 = true velocity"))

rtm5 = rtm_isic(m, m)
rtm5 = adjoint(reshape(rtm5, (n2, n1) ))
save(plotsdir(figfolder, "rtm_true_vel.jld"), "m", m, "m0", m, "rtm", rtm5)

####################################################################################################
# Load data for Scenario 1-5
m0_all, rtm_all = load(plotsdir(figfolder, "rtm_post_samples.jld"), "m0_all", "rtm_all")
m, m0, rtm2 = load(plotsdir(figfolder, "rtm_post_mean.jld"), "m", "m0", "rtm")
m, m_lin, rtm3 = load(plotsdir(figfolder, "rtm_lin_vel.jld"), "m", "m0", "rtm")
m, m_smooth, rtm4 = load(plotsdir(figfolder, "rtm_smooth_vel.jld"), "m", "m0", "rtm")
m, m_true, rtm5 = load(plotsdir(figfolder, "rtm_true_vel.jld"), "m", "m0", "rtm")

# Figures in slowness squared (s^2/km^2)
figure(); imshow(m, cmap="jet", aspect = "auto"); plt.colorbar(); title(L"True velocity ($s^2/km^2$)")
savefig(plotsdir(figfolder, "vel_m.png"))

figure(figsize=[20,8])
ax1 = subplot(2, 4, 1); imshow(m0_all[:,:,1], cmap="jet", aspect = "auto"); title("Smoothed posterior sample1")
ax2 = subplot(2, 4, 2); imshow(m0_all[:,:,2], cmap="jet", aspect = "auto"); title("Smoothed posterior sample2")
ax3 = subplot(2, 4, 3); imshow(m0_all[:,:,3], cmap="jet", aspect = "auto"); title("Smoothed posterior sample3")
ax4 = subplot(2, 4, 4); imshow(m0, cmap="jet", aspect = "auto"); title("Smoothed posterior mean")
ax5 = subplot(2, 4, 5); plot_rtm(adjoint(rtm_all[:,:,1]), d; new_fig=false, aspect = "auto"); title("RTM using posterior sample1")
ax6 = subplot(2, 4, 6); plot_rtm(adjoint(rtm_all[:,:,2]), d; new_fig=false, aspect = "auto"); title("RTM using posterior sample2")
ax7 = subplot(2, 4, 7); plot_rtm(adjoint(rtm_all[:,:,3]), d; new_fig=false, aspect = "auto"); title("RTM using posterior sample3")
ax8 = subplot(2, 4, 8); plot_rtm(adjoint(rtm2), d; new_fig=false, aspect = "auto"); title("RTM using posterior mean")
savefig(plotsdir(figfolder, "rtm_m_post.png"))

figure(figsize=[20,8])
ax1 = subplot(2, 4, 1); imshow(m0_all[:,:,1], cmap="jet", aspect = "auto"); title("Smoothed posterior sample")
ax2 = subplot(2, 4, 2); imshow(m0, cmap="jet", aspect = "auto"); title("Smoothed posterior mean")
ax3 = subplot(2, 4, 3); imshow(m_lin, cmap="jet", aspect = "auto"); title("Migration velocity in RTM")
# ax4 = subplot(2, 4, 4); imshow(m_smooth, cmap="jet", aspect = "auto"); title("Smoothed true velocity")
ax4 = subplot(2, 4, 4); imshow(m_true, cmap="jet", aspect = "auto"); title("Smoothed true velocity")
ax5 = subplot(2, 4, 5); plot_rtm(adjoint(rtm_all[:,:,1]), d; new_fig=false, aspect = "auto"); title("RTM using posterior sample")
ax6 = subplot(2, 4, 6); plot_rtm(adjoint(rtm2), d; new_fig=false, aspect = "auto"); title("RTM using posterior mean")
ax7 = subplot(2, 4, 7); plot_rtm(adjoint(rtm3), d; new_fig=false, aspect = "auto"); title("RTM using linearly varying velocity")
# ax8 = subplot(2, 4, 8); plot_rtm(adjoint(rtm4), d; new_fig=false, aspect = "auto"); title("RTM using smootehd true velocity")
ax8 = subplot(2, 4, 8); plot_rtm(adjoint(rtm5), d; new_fig=false, aspect = "auto"); title("RTM using true velocity")
savefig(plotsdir(figfolder, "rtm_m.png"))

# Figures in velocity (km/s)
v = sqrt.(1f0 ./ m)
v0_all = sqrt.(1f0 ./ m0_all[:,:,1:3])
v0 = sqrt.(1f0 ./ m0)
v_lin = sqrt.(1f0 ./ m_lin)
v_smooth = sqrt.(1f0 ./ m_smooth)
v_true = sqrt.(1f0 ./ m_true)

figure(); imshow(v, vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); plt.colorbar(); title("True velocity (km/s)")
savefig(plotsdir(figfolder, "vel.png"))

figure(figsize=[20,8])
ax1 = subplot(2, 4, 1); imshow(v0_all[:,:,1], vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); title("Smoothed posterior sample1")
ax2 = subplot(2, 4, 2); imshow(v0_all[:,:,2], vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); title("Smoothed posterior sample2")
ax3 = subplot(2, 4, 3); imshow(v0_all[:,:,3], vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); title("Smoothed posterior sample3")
ax4 = subplot(2, 4, 4); imshow(v0, vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); title("Smoothed posterior mean")
ax5 = subplot(2, 4, 5); plot_rtm(adjoint(rtm_all[:,:,1]), d; new_fig=false, aspect = "auto"); title("RTM using posterior sample1")
ax6 = subplot(2, 4, 6); plot_rtm(adjoint(rtm_all[:,:,2]), d; new_fig=false, aspect = "auto"); title("RTM using posterior sample2")
ax7 = subplot(2, 4, 7); plot_rtm(adjoint(rtm_all[:,:,3]), d; new_fig=false, aspect = "auto"); title("RTM using posterior sample3")
ax8 = subplot(2, 4, 8); plot_rtm(adjoint(rtm2), d; new_fig=false, aspect = "auto"); title("RTM using posterior mean")
savefig(plotsdir(figfolder, "rtm_v_post.png"))

figure(figsize=[20,8])
ax1 = subplot(2, 4, 1); imshow(v0_all[:,:,1], vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); title("Smoothed posterior sample")
ax2 = subplot(2, 4, 2); imshow(v0, vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); title("Smoothed posterior mean")
ax3 = subplot(2, 4, 3); imshow(v_lin, vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); title("Migration velocity in RTM")
# ax4 = subplot(2, 4, 4); imshow(v_smooth, vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); title("Smoothed true velocity")
ax4 = subplot(2, 4, 4); imshow(v_true, vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); title("True velocity")
ax5 = subplot(2, 4, 5); plot_rtm(adjoint(rtm_all[:,:,1]), d; new_fig=false, aspect = "auto"); title("RTM using posterior sample")
ax6 = subplot(2, 4, 6); plot_rtm(adjoint(rtm2), d; new_fig=false, aspect = "auto"); title("RTM using posterior mean")
ax7 = subplot(2, 4, 7); plot_rtm(adjoint(rtm3), d; new_fig=false, aspect = "auto"); title("RTM using used for INN training")
# ax8 = subplot(2, 4, 8); plot_rtm(adjoint(rtm4), d; new_fig=false, aspect = "auto"); title("RTM using smootehd true velocity")
ax8 = subplot(2, 4, 8); plot_rtm(adjoint(rtm5), d; new_fig=false, aspect = "auto"); title("RTM using true velocity")
savefig(plotsdir(figfolder, "rtm_v.png"))

figure(figsize=[20,4])
ax1 = subplot(1, 3, 1); imshow(v, cmap="jet", aspect = "auto"); plt.colorbar(); title("True velocity (km/s)")
ax2 = subplot(1, 3, 2); imshow(v0_all[:,:,1], cmap="jet", aspect = "auto"); plt.colorbar(); title("Migration velocity (km/s)")
ax3 = subplot(1, 3, 3); plot_rtm(adjoint(rtm_all[:,:,1]), d; new_fig=false, aspect = "auto"); title("RTM using posterior samples")
savefig(plotsdir(figfolder, "rtm_post_samples.png"))

figure(figsize=[20,4])
ax1 = subplot(1, 3, 1); imshow(v, cmap="jet", aspect = "auto"); plt.colorbar(); title("True velocity (km/s)")
ax2 = subplot(1, 3, 2); imshow(v0, cmap="jet", aspect = "auto"); plt.colorbar(); title("Migration velocity (km/s)")
ax3 = subplot(1, 3, 3); plot_rtm(adjoint(rtm2), d; new_fig=false, aspect = "auto"); title("RTM using posterior mean")
savefig(plotsdir(figfolder, "rtm_post_mean.png"))

figure(figsize=[20,4])
ax1 = subplot(1, 3, 1); imshow(v, cmap="jet", aspect = "auto"); plt.colorbar(); title("True velocity (km/s)")
ax2 = subplot(1, 3, 2); imshow(v_lin, cmap="jet", aspect = "auto"); plt.colorbar(); title("Migration velocity (km/s)")
ax3 = subplot(1, 3, 3); plot_rtm(adjoint(rtm3), d; new_fig=false, aspect = "auto"); title("RTM using linearly varying velocity")
savefig(plotsdir(figfolder, "rtm_lin_vel.png"))

figure(figsize=[20,4])
ax1 = subplot(1, 3, 1); imshow(v, cmap="jet", aspect = "auto"); plt.colorbar(); title("True velocity (km/s)")
ax2 = subplot(1, 3, 2); imshow(v_smooth, cmap="jet", aspect = "auto"); plt.colorbar(); title("Migration velocity (km/s)")
ax3 = subplot(1, 3, 3); plot_rtm(adjoint(rtm4), d; new_fig=false, aspect = "auto"); title("RTM using smoothed true velocity")
savefig(plotsdir(figfolder, "rtm_smooth_vel.png"))

figure(figsize=[20,4])
ax1 = subplot(1, 3, 1); imshow(v, cmap="jet", aspect = "auto"); plt.colorbar(); title("True velocity (km/s)")
ax2 = subplot(1, 3, 2); imshow(v_true, cmap="jet", aspect = "auto"); plt.colorbar(); title("Migration velocity (km/s)")
ax3 = subplot(1, 3, 3); plot_rtm(adjoint(rtm5), d; new_fig=false, aspect = "auto"); title("RTM using true velocity")
savefig(plotsdir(figfolder, "rtm_true_vel.png"))

figure()
plot( hcat(rtm_all[idx_wb+1:end,100,1], rtm2[idx_wb+1:end,100], rtm3[idx_wb+1:end,100], rtm5[idx_wb+1:end,100]) );
legend(["posterior samples", "posterior mean", "linearly varying velocity", "true velocity"])
title("Vertical profiles of RTM images")
savefig(plotsdir(figfolder, "rtm_vertical_profiles.png"))


text:
export binary_quadratic_form, can_solve, discriminant, isdiscriminant, QuadBin,
       isfundamental_discriminant, prime_form, prime_power_form, cycle,
       isindefinite, ispositive_definite, isnegative_definite, isreduced

###############################################################################
#
#   Constructor
#
###############################################################################

@doc Markdown.doc"""
     binary_quadratic_form(a, b, c)
Constructs the integral binary quadratic form $a$ * x^2 + $b$ * x * y + $c$ * y^2.
"""
binary_quadratic_form(a, b, c) = QuadBin(a, b, c)

###############################################################################
#
#   Coefficients
#
###############################################################################

function getindex(f::QuadBin, i::Int)
  if i == 1
    return f.a
  elseif i == 2
    return f.b
  elseif i == 3
    return f.c
  else
    throw(error("Index must be between 1 and 3"))
  end
end

###############################################################################
#
#   Values
#
###############################################################################

(f::QuadBin)(x,y) = f[1]*x^2 + f[2] * x * y + f[3] * y^2

@doc Markdown.doc"""
     can_solve(f::QuadBin, n::Int64)
For a binary quadratic form $f$ with negative discriminant and an integer $n$,
returns tuple (`true`, (`x`,`y`)) if $f$(`x`,`y`) = $n$ for integers `x`,`y`.
If no such integers exist return (`false`,(0,0))
"""
function can_solve(f::QuadBin, n::Int64)
    if discriminant(f) >= 0
        @error("$f has non-negative discriminant")
    end
    for y in 1:Int64(floor(sqrt(Int64(4*f[1]*n)//abs(Int64(discriminant(f))))))
        #now f(x,y) quadratic in one variable -> use quadratic formula
        aq = f[1]
        bq = f[2] * y
        cq = f[3] * y^2 - n
        d = bq^2 - 4*aq*cq
        if issquare(d)
            if divides(-bq + sqrt(d), 2*aq)[1]
                return(true, (divexact(-bq + sqrt(d), 2*aq), ZZ(y)))
            end
            if divides(-bq - sqrt(d), 2*aq)[1]
                return(true, (divexact(-bq - sqrt(d), 2*aq), ZZ(y)))
            end
        end
    end
    return (false,(ZZ(0),ZZ(0)))
end

###############################################################################
#
#   Discriminant
#
###############################################################################

discriminant(f::QuadBin) = f[2]^2 - 4 * f[1] * f[3]

@doc Markdown.doc"""
     isdiscriminant(D)
Returns `true` if $D$ is the discriminant of an integral binary
quadratic form, otherwise returns `false`.
"""
function isdiscriminant(D)
    if D == 0
        return false
    end
    m = mod(D,4)
    if m == 0 || m == 1
        return true
    end
    return false
end

@doc Markdown.doc"""
     isfundamental_discriminant(D)
Returns `true` if $D$ is a fundamental discriminant i.e. has conductor `1`,
otherwise returns `false`.
"""
function isfundamental_discriminant(D)
    m = mod(D,4)
    if m == 1 && issquarefree(D)
        return true
    end
    if m == 0
        h = divexact(D,4)
        c = mod(h,4)
        if (c == 2 || c == 3) && issquarefree(h)
            return true
        end
    end
    return false
end

@doc Markdown.doc"""
     conductor(D)
Returns the conductor of the discriminant $D$ i.e. the largest
positive integer `c` such that $\frac{D}{c^2}$ is a discriminant.
"""
function conductor(D)
  !isdiscriminant(D) && error("$D is no discriminant")
  d = divexact(D, fundamental_discriminant(D))
  return isqrt(d)
end

function fundamental_discriminant(D)
  fac = factor(D)
  sqf = one(FlintZZ)
  for (p, e) in fac
    if isodd(e)
      sqf = sqf * p
    end
  end
  # sqf = is the squarefree-part, so D = sqf * square and sqf square-free
  if mod(sqf, 4) == 1
    return sqf
  else
    return sqf * 4
  end
end

###############################################################################
#
#   Equality
#
###############################################################################

function ==(f1::QuadBin, f2::QuadBin)
  return f1[1] == f2[1] && f1[2] == f2[2] && f1[3] == f2[3]
end

###############################################################################
#
#   Arithmetic
#
###############################################################################

conjugate(f::QuadBin) = binary_quadratic_form(f[1], -f[2], f[3])

-(f::QuadBin) = binary_quadratic_form(-f[1], -f[2], -f[3])

Generic.content(f::QuadBin) = gcd([f[1], f[2], f[3]])

isindefinite(f::QuadBin) = discriminant(f) > 0 ? true : false

isnegative_definite(f::QuadBin) = (discriminant(f) < 0 && f[1] < 0)

ispositive_definite(f::QuadBin) = (discriminant(f) < 0 && f[1] > 0)

Base.iszero(f::QuadBin) = (f[1] == 0 && f[2] == 0 && f[3] == 0)

isprimitive(f::QuadBin) = (content(f) == 1)

isreducible(f::QuadBin) = issquare(discriminant(f))

function isreduced(f::QuadBin)
    if ispositive_definite(f)
        return (abs(f[2]) <= f[1] <= f[3]) && (f[2] >= 0 || (f[1] != f[2] && f[1] != f[3]))
    end

    if isnegative_definite(f)
        return isreduced(-conjugate(f))
    end

    D = Int64(discriminant(f))

    if isindefinite(f)
        h = (abs(sqrt(D) - float(Int64(2*abs(f[1])))) < f[2] < sqrt(D)) || f[1] == 0 && -f[2] < 2*f[3] <= f[2]
        return (h || f[3]==0 && -f[2] < 2*f[1] <= f[2])
    end
end

###############################################################################
#
#   Reduction
#
###############################################################################

@doc Markdown.doc"""
     reduce(f::QuadBin)
Returns for the positive (or negative) definite binary quadratic form $f$ an equivalent
reduced form.
"""
function Base.reduce(f::QuadBin)
    if isindefinite(f)
        @error("Not implemented for $f indefinite")
    end
    if isnegative_definite(f)
        return reduce(-conjugate(f))
    end
    if ispositive_definite(f)
        return _reduce(f)[1]
    end
end

function _rho(f::QuadBin, T::fmpz_mat)
    s = _round_erec(f[2]//(2*f[3]))
    res_f = binary_quadratic_form(f[3], -f[2] + 2*s*f[3], f[3]*s^2 - f[2]*s + f[1])
    return(res_f, matrix(ZZ, 2, 2, [T[2], T[1] + s*T[2], T[4], T[3] + s*T[4]]))
end

function _normalize(f::QuadBin)
    s = _round_erec((f[1]-f[2])//(2*f[1]), dir = "down")
    res_f = binary_quadratic_form(f[1], f[2] + 2*s*f[1], f[1]*s^2 + f[2]*s + f[3])
    return(res_f, matrix(ZZ, 2, 2, [1, s, 0, 1]))
end

function _reduce(f::QuadBin)
    (g, T) = _normalize(f)
    while !_isreduced(g)
        (g, T) = _rho(g, T)
    end
    return (g, T)
end

function _isreduced(f::QuadBin)
    bol1 = f[1] <= f[3]
    bol2 = true
    if f[1] == f[3]
        bol2 = 0 <= f[2]
    end
    #reduced is checked, could be omitted
    bol3 = -f[1] < f[2] && f[2] <= f[1]

    return bol1 && bol2 && bol3
end

###############################################################################
#
#   Other Utilities
#
###############################################################################

function _round_erec(x::fmpq; dir = "none")
    if dir == "up"
        return fmpz(Int64(ceil(Rational{Int64}(x.num, x.den))))
    end
    if dir == "down"
        return fmpz(Int64(floor(Rational{Int64}(x.num, x.den))))
    end
    if dir == "none"
        if abs(x-_round_erec(x, dir = "up")) <= abs(x-_round_erec(x, dir = "down"))
            return _round_erec(x, dir = "up")
        else
            return _round_erec(x, dir = "down")
        end
    end
end

###############################################################################
#
#   Composition
#
###############################################################################

#COR7.3.17: f1,f2  same discriminant
# function Hecke.compose(f1::QuadBin, f2::QuadBin)
#     d = discriminant(f1)
#     Bm = numerator((b(f1) + b(f2))//2)
#     m = gcd(gcd(a(f1), a(f2)), Bm)
#     A = numerator((a(f1)*a(f2))//m^2)
#     #better?
#     j=1
#     k=1
#     l=1
#     while true
#         j = rand(-50:50)
#         k = rand(-50:50)
#         l = rand(-50:50)
#         if j*a(f2) + k*a(f1) + l*Bm == m
#             break
#         end
#     end
#     s = numerator(j*a(f2)*b(f1) + k*a(f1)*b(f2) + numerator((l*(b(f1)*b(f2)+d))//2)//m)
#     B = mod(s,2*A)
#     C = numerator((B^2-d)//(4*A))
#     return(QuadBin(A, B, C))
# end

###COMPOSITION###(Dirichlet)
#using Dirichlet composition, compare warwick Def 2.10

@doc Markdown.doc"""
     compose(f1::QuadBin, f2::QuadBin)
Returns the composition of the binary quadratic forms $f1$ and $f2$.
The result is not reduced, uses Dirichlet Composition.
"""
function compose(f1::QuadBin, f2::QuadBin)
    #discriminants have to match
    D = discriminant(f1)
    D2 = discriminant(f2)
    if !(D==D2)
        @error("discriminants do not match")
    end
    (h, n_1, n_2) = gcdx(f1[1], f2[1])
    (e, n_3, n3) = gcdx(h, divexact(f1[2]+f2[2], 2))

    (n1, n2) = (n_3 * n_1, n_3 * n_2)
    B = n1 * divexact(f1[1] * f2[2], e) + n2 * divexact(f1[2] * f2[1], e) + n3 * divexact(f1[2] * f2[2] + D, 2*e)
    return QuadBin(divexact(f1[1]*f2[1], e^2), B, divexact(e^2*(B^2-D), 4*f1[1]*f2[1]))
end

# cmpr Buchmann Alg 6.1
@doc Markdown.doc"""
     cycle(f::QuadBin)
Returns the cycle of the irreducible indefinite reduced form $f$.
"""
function cycle(f::QuadBin)
    if !isindefinite(f)
        @error("$f not indefinite")
    elseif !isreduced(f)
        @error("$f not reduced")
    elseif isreducible(f)
        @error("$f reducible")
    end

    Cycle = Vector([f])
    g = _rhoTau(f)
    while !(g==f)
        push!(Cycle, g)
        g = _rhoTau(g)
    end
    return Cycle
end

function _rhoTau(f::QuadBin)
    #s = round(f[2]//(2*f[3]))
    num = Int64(f[2])+Int64(floor(sqrt(Int64(discriminant(f)))))
    den = Int64(2*abs(f[3]))
    s =  Int64(Int64(floor(num/den)))
    return QuadBin(abs(f[3]), -f[2] + 2*s*abs(f[3]), -(f[1] + f[2]*s + f[3]*s^2))
end

###############################################################################
#
#   Prime Forms
#
###############################################################################

function _sqrtmod4P(d::fmpz, p::fmpz)
    if jacobi(mod(d, p), p) == -1
        @error("$d is no square modulo $p")
    end
    if p == 2
        if iseven(d)
            return 2 * mod(divexact(d, 4), 2)
        else
            return 1
        end
    else
        r = sqrtmod(d, p)
        if mod(r,2) == mod(d,2)
            return r
        else
            return p-r
        end
    end
end

function _number_of_primeforms(d::fmpz, p::fmpz)
    return jacobi(mod(d, p), p) + 1
end

@doc Markdown.doc"""
     prime_form(d::fmpz, p::fmpz)
Returns an integral binary quadratic form of discriminant $d$ and leading coefficient
$p$ where $p$ is a prime number.
"""
function prime_form(d::fmpz, p::fmpz)
    if !isdiscriminant(d)
        @error("$d is no discriminant")
    end
    if _number_of_primeforms(d, p) == 0
        @error("prime form does not exist")
    end
    b = _sqrtmod4P(d, p)
    return QuadBin(p, b, divexact(b^2-d, 4*p))
end

function _sqrtmod4PE(d::fmpz, p::fmpz, e)
    if mod(d,p) == 0
        n = 1
        while (2*n + 2  <= e && mod(p^(2n+2), d) == 0)
            n = n + 1
        end
        if (p == 2 && !mod(2^(-2*n)*d, 4) in [0,1])
            n = n -2
        end
        d0 = divexact(d, p^(2*n))
        e0 = e - 2*n
        if e0 == 0
            return p^n*mod(d0, 2)
        else return p^n * _sqrtmod4PE(d0, p^e0, e)
        end
    else
        b = _sqrtmod4P(d, p)
        f = 1
        while f < e
            k = mod(divexact(d-b^2, b*4*p^f), p)
            b = mod(b + 2*k*p^f, 2*p^(f+1))
            f = f + 1
        end
        if b > p^e
            return b-2*p^e
        else
            return b
        end
    end
end

@doc Markdown.doc"""
     prime_power_form(d::fmpz, p::fmpz, e)
Returns an integral binary quadratic form of discriminant $d$ and leading coefficient
$p$^$e$ where $p$ is a prime number.
"""
function prime_power_form(d::fmpz, p::fmpz, e)
    if !isdiscriminant(d)
        @error("$d is no discriminant")
    end
    b = _sqrtmod4PE(d, p, e)
    c = divexact(b^2-d, 4*p^e)
    f = QuadBin(p^e, b, c)
    if discriminant(f) != d
        @error("prime power form does not exist")
    end
    return f
end


text:
# conversion routines between different sparse and dense storage formats

SparseArrays.sparse(::DenseCuArray, args...) = error("CUSPARSE supports multiple sparse formats, use specific constructors instead (e.g. CuSparseMatrixCSC)")

## CSR to CSC

# by flipping rows and columns, we can use that to get CSC to CSR too

for (fname,elty) in ((:cusparseScsr2csc, :Float32),
                     (:cusparseDcsr2csc, :Float64),
                     (:cusparseCcsr2csc, :ComplexF32),
                     (:cusparseZcsr2csc, :ComplexF64))
    @eval begin
        function CuSparseMatrixCSC{$elty}(csr::CuSparseMatrixCSR{$elty}; inda::SparseChar='O')
            m,n = csr.dims
            colPtr = CUDA.zeros(Cint, n+1)
            rowVal = CUDA.zeros(Cint, nnz(csr))
            nzVal = CUDA.zeros($elty, nnz(csr))
            if version() >= v"10.2"
                # TODO: algorithm configuratibility?
                @workspace size=@argout(
                        cusparseCsr2cscEx2_bufferSize(handle(), m, n, nnz(csr), nonzeros(csr),
                            csr.rowPtr, csr.colVal, nzVal, colPtr, rowVal,
                            $elty, CUSPARSE_ACTION_NUMERIC, inda,
                            CUSPARSE_CSR2CSC_ALG1, out(Ref{Csize_t}(1)))
                    )[] buffer->begin
                        cusparseCsr2cscEx2(handle(), m, n, nnz(csr), nonzeros(csr),
                            csr.rowPtr, csr.colVal, nzVal, colPtr, rowVal,
                            $elty, CUSPARSE_ACTION_NUMERIC, inda,
                            CUSPARSE_CSR2CSC_ALG1, buffer)
                    end
            else
                $fname(handle(), m, n, nnz(csr), nonzeros(csr),
                    csr.rowPtr, csr.colVal, nzVal, rowVal,
                    colPtr, CUSPARSE_ACTION_NUMERIC, inda)
            end
            CuSparseMatrixCSC(colPtr,rowVal,nzVal,csr.dims)
        end

        function CuSparseMatrixCSR{$elty}(csc::CuSparseMatrixCSC{$elty}; inda::SparseChar='O')
            m,n    = csc.dims
            rowPtr = CUDA.zeros(Cint,m+1)
            colVal = CUDA.zeros(Cint,nnz(csc))
            nzVal  = CUDA.zeros($elty,nnz(csc))
            if version() >= v"10.2"
                # TODO: algorithm configuratibility?
                @workspace size=@argout(
                        cusparseCsr2cscEx2_bufferSize(handle(), n, m, nnz(csc), nonzeros(csc),
                            csc.colPtr, rowvals(csc), nzVal, rowPtr, colVal,
                            $elty, CUSPARSE_ACTION_NUMERIC, inda,
                            CUSPARSE_CSR2CSC_ALG1, out(Ref{Csize_t}(1)))
                    )[] buffer->begin
                        cusparseCsr2cscEx2(handle(), n, m, nnz(csc), nonzeros(csc),
                            csc.colPtr, rowvals(csc), nzVal, rowPtr, colVal,
                            $elty, CUSPARSE_ACTION_NUMERIC, inda,
                            CUSPARSE_CSR2CSC_ALG1, buffer)
                    end
            else
                $fname(handle(), n, m, nnz(csc), nonzeros(csc),
                    csc.colPtr, rowvals(csc), nzVal, colVal,
                    rowPtr, CUSPARSE_ACTION_NUMERIC, inda)
            end
            CuSparseMatrixCSR(rowPtr,colVal,nzVal,csc.dims)
        end
    end
end

## CSR to BSR and vice-versa

for (fname,elty) in ((:cusparseScsr2bsr, :Float32),
                     (:cusparseDcsr2bsr, :Float64),
                     (:cusparseCcsr2bsr, :ComplexF32),
                     (:cusparseZcsr2bsr, :ComplexF64))
    @eval begin
        function CuSparseMatrixBSR{$elty}(csr::CuSparseMatrixCSR{$elty}, blockDim::Integer;
                                          dir::SparseChar='R', inda::SparseChar='O',
                                          indc::SparseChar='O')
            m,n = csr.dims
            nnz_ref = Ref{Cint}(1)
            mb = div((m + blockDim - 1),blockDim)
            nb = div((n + blockDim - 1),blockDim)
            bsrRowPtr = CUDA.zeros(Cint,mb + 1)
            cudesca = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, inda)
            cudescc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, indc)
            cusparseXcsr2bsrNnz(handle(), dir, m, n, cudesca, csr.rowPtr,
                                csr.colVal, blockDim, cudescc, bsrRowPtr, nnz_ref)
            bsrNzVal = CUDA.zeros($elty, nnz_ref[] * blockDim * blockDim )
            bsrColInd = CUDA.zeros(Cint, nnz_ref[])
            $fname(handle(), dir, m, n,
                   cudesca, nonzeros(csr), csr.rowPtr, csr.colVal,
                   blockDim, cudescc, bsrNzVal, bsrRowPtr,
                   bsrColInd)
            CuSparseMatrixBSR{$elty}(bsrRowPtr, bsrColInd, bsrNzVal, csr.dims, blockDim, dir, nnz_ref[])
        end
    end
end

for (fname,elty) in ((:cusparseSbsr2csr, :Float32),
                     (:cusparseDbsr2csr, :Float64),
                     (:cusparseCbsr2csr, :ComplexF32),
                     (:cusparseZbsr2csr, :ComplexF64))
    @eval begin
        function CuSparseMatrixCSR{$elty}(bsr::CuSparseMatrixBSR{$elty};
                                          inda::SparseChar='O', indc::SparseChar='O')
            m,n = bsr.dims
            mb = div(m,bsr.blockDim)
            nb = div(n,bsr.blockDim)
            nnzVal = nnz(bsr) * bsr.blockDim * bsr.blockDim
            cudesca = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, inda)
            cudescc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, indc)
            csrRowPtr = CUDA.zeros(Cint, m + 1)
            csrColInd = CUDA.zeros(Cint, nnzVal)
            csrNzVal  = CUDA.zeros($elty, nnzVal)
            $fname(handle(), bsr.dir, mb, nb,
                   cudesca, nonzeros(bsr), bsr.rowPtr, bsr.colVal,
                   bsr.blockDim, cudescc, csrNzVal, csrRowPtr,
                   csrColInd)
            CuSparseMatrixCSR(csrRowPtr, csrColInd, csrNzVal, bsr.dims)
        end
    end
end

## CSR to COO and vice-versa

function CuSparseMatrixCSR(coo::CuSparseMatrixCOO{Tv}, ind::SparseChar='O') where {Tv}
    m,n = coo.dims
    nnz = coo.nnz
    csrRowPtr = CUDA.zeros(Cint, nnz)
    cusparseXcoo2csr(handle(), coo.rowInd, nnz, m, csrRowPtr, ind)
    CuSparseMatrixCSR{Tv}(csrRowPtr, coo.colInd, coo.nzVal, coo.dims)
end

function CuSparseMatrixCOO(csr::CuSparseMatrixCSR{Tv}, ind::SparseChar='O') where {Tv}
    m,n = csr.dims
    nnz = csr.nnz
    cooRowInd = CUDA.zeros(Cint, nnz)
    cusparseXcsr2coo(handle(), csr.rowPtr, nnz, m, cooRowInd, ind)
    CuSparseMatrixCOO{Tv}(cooRowInd, csr.colVal, csr.nzVal, csr.dims, nnz)
end

## sparse to dense, and vice-versa

# TODO: use cusparseSparseToDense & cusparseDenseToSparse (3x faster, CUDA 11.1 Update 1)

for (cname,rname,elty) in ((:cusparseScsc2dense, :cusparseScsr2dense, :Float32),
                           (:cusparseDcsc2dense, :cusparseDcsr2dense, :Float64),
                           (:cusparseCcsc2dense, :cusparseCcsr2dense, :ComplexF32),
                           (:cusparseZcsc2dense, :cusparseZcsr2dense, :ComplexF64))
    @eval begin
        function CUDA.CuMatrix{$elty}(csr::CuSparseMatrixCSR{$elty}; ind::SparseChar='O')
            m,n = csr.dims
            denseA = CUDA.zeros($elty,m,n)
            lda = max(1,stride(denseA,2))
            cudesc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, ind)
            $rname(handle(), m, n, cudesc, nonzeros(csr),
                   csr.rowPtr, csr.colVal, denseA, lda)
            denseA
        end
        function CUDA.CuMatrix{$elty}(csc::CuSparseMatrixCSC{$elty}; ind::SparseChar='O')
            m,n = csc.dims
            denseA = CUDA.zeros($elty,m,n)
            lda = max(1,stride(denseA,2))
            cudesc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, ind)
            $cname(handle(), m, n, cudesc, nonzeros(csc),
                   rowvals(csc), csc.colPtr, denseA, lda)
            denseA
        end
    end
end

for (nname,cname,rname,elty) in ((:cusparseSnnz, :cusparseSdense2csc, :cusparseSdense2csr, :Float32),
                                 (:cusparseDnnz, :cusparseDdense2csc, :cusparseDdense2csr, :Float64),
                                 (:cusparseCnnz, :cusparseCdense2csc, :cusparseCdense2csr, :ComplexF32),
                                 (:cusparseZnnz, :cusparseZdense2csc, :cusparseZdense2csr, :ComplexF64))
    @eval begin
        function CuSparseMatrixCSR(A::CuMatrix{$elty}; ind::SparseChar='O')
            m,n = size(A)
            lda = max(1, stride(A,2))
            cudesc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL,
                                        CUSPARSE_FILL_MODE_LOWER,
                                        CUSPARSE_DIAG_TYPE_NON_UNIT, ind)
            nnzRowCol = CUDA.zeros(Cint, m)
            nnzTotal = Ref{Cint}(1)
            $nname(handle(),
                   'R', m, n, cudesc, A, lda, nnzRowCol,
                   nnzTotal)
            nzVal = CUDA.zeros($elty,nnzTotal[])

            rowPtr = CUDA.zeros(Cint,m+1)
            colInd = CUDA.zeros(Cint,nnzTotal[])
            $rname(handle(), m, n, cudesc, A,
                    lda, nnzRowCol, nzVal, rowPtr, colInd)
            return CuSparseMatrixCSR(rowPtr,colInd,nzVal,size(A))
        end

        function CuSparseMatrixCSC(A::CuMatrix{$elty}; ind::SparseChar='O')
            m,n = size(A)
            lda = max(1, stride(A,2))
            cudesc = CuMatrixDescriptor(CUSPARSE_MATRIX_TYPE_GENERAL,
                                        CUSPARSE_FILL_MODE_LOWER,
                                        CUSPARSE_DIAG_TYPE_NON_UNIT, ind)
            nnzRowCol = CUDA.zeros(Cint, n)
            nnzTotal = Ref{Cint}(1)
            $nname(handle(),
                   'C', m, n, cudesc, A, lda, nnzRowCol,
                   nnzTotal)
            nzVal = CUDA.zeros($elty,nnzTotal[])

            colPtr = CUDA.zeros(Cint,n+1)
            rowInd = CUDA.zeros(Cint,nnzTotal[])
            $cname(handle(), m, n, cudesc, A,
                    lda, nnzRowCol, nzVal, rowInd, colPtr)
            return CuSparseMatrixCSC(colPtr,rowInd,nzVal,size(A))
        end
    end
end

function CUDA.CuMatrix{T}(bsr::CuSparseMatrixBSR{T}; inda::SparseChar='O',
                          indc::SparseChar='O') where {T}
    CuMatrix{T}(CuSparseMatrixCSR{T}(bsr; inda, indc))
end

function CuSparseMatrixBSR(A::CuMatrix; ind::SparseChar='O')
    m,n = size(A)   # TODO: always let the user choose, or provide defaults for other methods too
    CuSparseMatrixBSR(CuSparseMatrixCSR(A; ind), gcd(m,n))
end


text:
using Luxor

Drawing(600, 600, "hello-world.png")
origin()

sethue(30 / 100, 43 / 100, 96 / 100)
rad = 200
circle(0, 0, rad, :fill)

sethue("red")

cent_ellipse = Point(0, 0)
ellipse(cent_ellipse, 2 * rad, rad / 10, :stroke)

pole = Point(0, -rad)
arrow(cent_ellipse, pole, linewidth = 2.0)

x = Point(rad, 0)
arrow(cent_ellipse, x, linewidth = 2.0)

θ = rand() * 2π
input_state = Point(rad * cos(θ), rad * sin(θ))
line(cent_ellipse, input_state, :stroke)
circle(input_state, rad / 20, :fill)

θ = rand() * 2π
final_state = Point(-rad * cos(θ), rad * sin(θ))
line(cent_ellipse, final_state, :stroke)
circle(final_state, rad / 20, :fill)

# random state between the two other ones

curve(input_state, x, final_state)
sethue("black")
setdash("shortdashed")
# arrow(input_state, final_state, [15, 15])
strokepath()

finish()
preview()


text:
module TestTUI
#=
Remaining problems:
- bug in TerminalUserInterfaces.terminal_size()
- printing of outputs is cut off. See `paragraph.jl` in TerminalUserInterfaces

ToDo:
- also track files in `test` directory
=#

using Glob
using Revise
using Parameters
using Pkg
using Suppressor
using TerminalUserInterfaces
const TUI = TerminalUserInterfaces

const DEVDIR = get(ENV, "JULIA_PKG_DEVDIR", joinpath(homedir(), ".julia/dev"))

include("./utils.jl")
include("./test.jl")
include("./textwidget.jl")
include("./tui.jl")

function testtui(
        pattern = "*";
        dir = DEVDIR)
    testtui(findpackages(pattern, dir))
end

#=
Tasks:
1. watch package directories and requeue changed packages
2. empty queue and update testresults
3. update UI
=#

function startwatching(packagepaths, results, channel)
    for (packagepath, result) in zip(packagepaths, results)
        paths = sourcepaths(packagepath)
        @async begin
            entr(paths, pause = .5) do
                try
                    if result.state != Queued
                        result.state = Queued
                        push!(channel, packagepath)
                    end
                catch e
                    @error e
                end
            end
        end
    end
end

function starttesting(resultdict, channel)
    for packagepath in channel
        try
            result = resultdict[packagepath]
            result.state = Running
            newresult = runtests(packagepath)
            if result.state != Queued
                result.output = newresult.output
                result.state = newresult.state
            end
        catch e
            @error e pkg=packagepath
        end
    end
end

function testtui(packagepaths::Vector{String})
    @assert all(istestablepackage.(packagepaths))
    packagepaths = sort(packagepaths)
    names = [splitpath(path)[end] for path in packagepaths]
    results = [TestResult(name) for name in names]
    resultdict = Dict(zip(packagepaths, results))

    channel = Channel{String}(100)
    startwatching(packagepaths, results, channel)
    Threads.@spawn starttesting(resultdict, channel)

    tuiloop(results)

    return results, channel
end

export testtui, runtests

end


text:
function PoissonDistribution(arg0::jdouble)
    return PoissonDistribution((jdouble,), arg0)
end

function PoissonDistribution(arg0::jdouble, arg1::jdouble)
    return PoissonDistribution((jdouble, jdouble), arg0, arg1)
end

function PoissonDistribution(arg0::jdouble, arg1::jdouble, arg2::jint)
    return PoissonDistribution((jdouble, jdouble, jint), arg0, arg1, arg2)
end

function PoissonDistribution(arg0::jdouble, arg1::jint)
    return PoissonDistribution((jdouble, jint), arg0, arg1)
end

function cumulative_probability(obj::PoissonDistribution, arg0::jint)
    return jcall(obj, "cumulativeProbability", jdouble, (jint,), arg0)
end

function get_mean(obj::PoissonDistribution)
    return jcall(obj, "getMean", jdouble, ())
end

function get_numerical_mean(obj::PoissonDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::PoissonDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

function get_support_lower_bound(obj::PoissonDistribution)
    return jcall(obj, "getSupportLowerBound", jint, ())
end

function get_support_upper_bound(obj::PoissonDistribution)
    return jcall(obj, "getSupportUpperBound", jint, ())
end

function is_support_connected(obj::PoissonDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end

function log_probability(obj::PoissonDistribution, arg0::jint)
    return jcall(obj, "logProbability", jdouble, (jint,), arg0)
end

function normal_approximate_probability(obj::PoissonDistribution, arg0::jint)
    return jcall(obj, "normalApproximateProbability", jdouble, (jint,), arg0)
end

function probability(obj::PoissonDistribution, arg0::jint)
    return jcall(obj, "probability", jdouble, (jint,), arg0)
end


text:
using Parameters
using SpecialFunctions
using LinearAlgebra
using Tracker
using ReverseDiff
using RecursiveArrayTools

@inline function calcgibbs(ph::U,T::W) where {U<:IdealPhase,W<:Real}
    return getGibbs.(getfield.(ph.species,:thermo),T)
end
export calcgibbs

@inline function calcenthalpyinternalgibbs(ph::U,T::W,P::Z,V::Q) where {U<:IdealPhase,W,Z,Q<:Real}
    Hs = getEnthalpy.(getfield.(ph.species,:thermo),T)
    Us = Hs .- P*V
    Gs = Hs .- T.*getEntropy.(getfield.(ph.species,:thermo),T)
    return Hs,Us,Gs
end

@inline function calcenthalpyinternalgibbs(ph::Union{IdealGas,IdealSurface},T::W,P::Z,V::Q) where {W,Z,Q<:Real}
    Hs = getEnthalpy.(getfield.(ph.species,:thermo),T)
    Us = Hs .- R*T
    Gs = Hs .- T.*getEntropy.(getfield.(ph.species,:thermo),T)
    return Hs,Us,Gs
end

export calcenthalpyinternalgibbs

function makespcsvector(phase,spcdict)
    y0 = zeros(length(phase.species))
    spnames = [x.name for x in phase.species]
    for (key,val) in spcdict
        if key == "T"
            continue
        elseif key == "P"
            continue
        elseif key == "V"
            continue
        else
            ind = findfirst(isequal(key),spnames)
            @assert typeof(ind)<: Integer  "$key not found in species list: $spnames"
            y0[ind] = val
        end
    end
    return y0
end

export makespcsvector

@inline function getkf(rxn::ElementaryReaction,ph,T,P,C,ns,V,phi)
    if isdefined(rxn.kinetics,:efficiencies) && length(rxn.kinetics.efficiencies) > 0
        @views @inbounds @fastmath C += sum([ns[i]*val for (i,val) in rxn.kinetics.efficiencies])/V
    end
    return rxn.kinetics(T=T,P=P,C=C,phi=phi)
end
export getkf

@inline function getkfs(ph::U,T::W1,P::W2,C::W3,ns::Q,V::W4,phi) where {U,W1,W2,W3,W4<:Real,Q<:AbstractArray}
    kfs = similar(ns,length(ph.reactions))
    i = 1
    oldind = 1
    ind = 0
    while i <= length(ph.veckineticsinds) #vectorized kinetics
        @inbounds ind = ph.veckineticsinds[i]
        @inbounds kfs[oldind:ind] = ph.veckinetics[i](;T=T,P=P,C=C)
        oldind = ind+1
        i += 1
    end
    @simd for i in ind+1:length(ph.reactions)
        @inbounds kfs[i] = getkf(ph.reactions[i],ph,T,P,C,ns,V,phi)
    end
    return kfs
end

export getkfs

"""
Calculates the diffusion limited rate coefficient
for 1 spc returns Inf
for 2 spc calculates using the Smolchowski equation
for >2 spc calculates using the Generalized Smolchowski equation
Equations from Flegg 2016
"""
@inline function getDiffusiveRate(spcs::Q,diffs::Array{W,1}) where {Q<:AbstractArray,W<:Real}
    if length(spcs) == 1
        return Inf
    elseif length(spcs) == 2
        @fastmath @inbounds kf = 4.0*Base.pi*(diffs[spcs[1].index]+diffs[spcs[2].index])*(spcs[1].radius+spcs[2].radius)*Na
    else
        @views @inbounds diffusivity = diffs[getfield.(spcs,:index)]
        N = length(spcs)
        @fastmath a = (3.0*length(spcs)-5.0)/2.0
        @fastmath Dinv = 1.0./diffusivity
        @fastmath Dbar = 1.0./reverse(cumsum(Dinv))
        @fastmath Dhat = diffusivity .+ Dbar
        @fastmath @inbounds deltaN = sum(Dinv)/sum([sum([1.0/(diffusivity[i]*diffusivity[m]) for m in 1:N-1 if i>m]) for i in 2:N])
        @views @fastmath @inbounds kf = prod(Dhat[2:end].^1.5)*4*Base.pi^(a+1)/gamma(a)*(sum(getfield.(spcs,:radius))/sqrt(deltaN))^(2*a)*Na^(N-1)
    end
    return kf
end
export getDiffusiveRate

@inline function getKc(rxn::ElementaryReaction,ph::U,T::Z,Gs::Q,phi::V=0.0) where {U<:AbstractPhase,V,Q,Z<:Real}
    Nreact = length(rxn.reactantinds)
    Nprod = length(rxn.productinds)
    dGrxn = 0.0
    if Nreact == 1
        @fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]]
    elseif Nreact == 2
        @fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]]+Gs[rxn.reactantinds[2]]
    elseif Nreact == 3
        @fastmath @inbounds dGrxn -= Gs[rxn.reactantinds[1]]+Gs[rxn.reactantinds[2]]+Gs[rxn.reactantinds[3]]
    end
    if Nprod == 1
        @fastmath @inbounds dGrxn += Gs[rxn.productinds[1]]
    elseif Nprod == 2
        @fastmath @inbounds dGrxn += Gs[rxn.productinds[1]]+Gs[rxn.productinds[2]]
    elseif Nprod == 3
        @fastmath @inbounds dGrxn += Gs[rxn.productinds[1]]+Gs[rxn.productinds[2]]+Gs[rxn.productinds[3]]
    end
    return @inbounds @fastmath exp(-(dGrxn+rxn.electronchange*phi)/(R*T))*(getC0(ph,T))^(Nprod-Nreact)
end

@inline function getKc(rxn::ElementaryReaction,phase1,phase2,Gs1,Gs2,T,phi=0.0) #for constant k interfaces
    dGrxn = 0.0
    dN1 = 0
    dN2 = 0
    for r in rxn.reactants
        isfirst = true
        ind = findfirst(isequal(r),phase1.species)
        if ind === nothing
            isfirst = false
            ind = findfirst(isequal(r),phase2.species)
            dGrxn -= Gs2[ind]
            dN2 -= 1
        else
            dGrxn -= Gs1[ind]
            dN1 -= 1
        end
    end
    for r in rxn.products
        isfirst = true
        ind = findfirst(isequal(r),phase1.species)
        if ind === nothing
            isfirst = false
            ind = findfirst(isequal(r),phase2.species)
            dGrxn += Gs2[ind]
            dN2 += 1
        else
            dGrxn += Gs1[ind]
            dN1 += 1
        end
    end
    return @inbounds @fastmath exp(-(dGrxn+rxn.electronchange*phi)/(R*T))*getC0(phase1,T)^dN1*getC0(phase2,T)^dN2
end
export getKc

@inline function getKcs(ph::U,T::Z,Gs::Q) where {U<:AbstractPhase,Q,Z<:Real}
    return @fastmath @inbounds exp.(ph.stoichmatrix*(Gs./(R*T)) .+ ph.Nrp.*log(getC0(ph,T)));
end

@inline function getKcs(ph::U,T::Z,Gs::Q,phi::V) where {U<:AbstractPhase,Q,Z<:Real,V<:Real}
    return @fastmath @inbounds exp.(ph.stoichmatrix*(Gs./(R*T)).+ph.electronchange.*(phi/(R*T)) .+ ph.Nrp.*log(getC0(ph,T)));
end

@inline function getKcs(ph,T,Gs1,Gs2)
    Gpart = ArrayPartition(Gs1,Gs2)
    return @fastmath @inbounds exp.(ph.stoichmatrix*(Gpart./(R*T)) .+ ph.Nrp1.*log(getC0(ph.domain1.phase,T)) .+ ph.Nrp2.*log(getC0(ph.domain2.phase,T)));
end

export getKcs

"""
Calculates the forward and reverse rate coefficients for a given reaction, phase and state
Maintains diffusion limitations if the phase has diffusionlimited=true
"""
@inline function getkfkrev(rxn::ElementaryReaction,ph::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W8;kf::W6=-1.0,f::W7=-1.0) where {U<:AbstractPhase,W8,W6,W7,W5,W4,W1,W2,W3<:Real,Q1,Q2,Q3<:AbstractArray}
    if signbit(kf)
        if signbit(f)
            kf = getkf(rxn,ph,T,P,C,ns,V,phi)
        else
            kf = getkf(rxn,ph,T,P,C,ns,V,phi)*f
        end
    end
    Kc = getKc(rxn,ph,T,Gs,phi)
    @fastmath krev = kf/Kc
    if ph.diffusionlimited
        if length(rxn.reactants) == 1
            if length(rxn.products) > 1
                krevdiff = getDiffusiveRate(rxn.products,diffs)
                @fastmath krev = krev*krevdiff/(krev+krevdiff)
                @fastmath kf = Kc*krev
            end
        elseif length(rxn.products) == 1
            kfdiff = getDiffusiveRate(rxn.reactants,diffs)
            @fastmath kf = kf*kfdiff/(kf+kfdiff)
            @fastmath krev = kf/Kc
        elseif length(rxn.products) == length(rxn.reactants)
            kfdiff = getDiffusiveRate(rxn.reactants,diffs)
            krevdiff = getDiffusiveRate(rxn.products,diffs)
            @fastmath kff = kf*kfdiff/(kf+kfdiff)
            @fastmath krevr = krev*krevdiff/(krev+krevdiff)
            @fastmath kfr = Kc*krevr
            if kff > kfr
                kf = kfr
                krev = krevr
            else
                kf = kff
                @fastmath krev = kf/Kc
            end
        end
    end
    if rxn.reversible
        return (kf,krev)
    else
        return (kf,0.0)
    end
end
export getkfkrev

@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W7;kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2,Q3<:AbstractArray}
    if !phase.diffusionlimited && kfs === nothing
        kfs = getkfs(phase,T,P,C,ns,V,phi)
        if phi == 0.0
            krev = @fastmath kfs./getKcs(phase,T,Gs).*phase.reversibility
        else
            krev = @fastmath kfs./getKcs(phase,T,Gs,phi).*phase.reversibility
        end
    elseif !phase.diffusionlimited && !(kfs === nothing)
        if phi == 0.0
            krev = @fastmath kfs./getKcs(phase,T,Gs).*phase.reversibility
        else
            krev = @fastmath kfs./getKcs(phase,T,Gs,phi).*phase.reversibility
        end
    elseif phase.diffusionlimited && !(kfs === nothing)
        len = length(phase.reactions)
        krev = zeros(typeof(N),len)
        @simd for i = 1:len
           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;kf=kfs[i])
        end
    else
        len = length(phase.reactions)
        kfs = zeros(typeof(N),len)
        krev = zeros(typeof(N),len)
        @simd for i = 1:len
           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi)
        end
    end
    return kfs,krev
end

@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Q2,diffs::Q3,V::W5,phi::W7;kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2<:Union{ReverseDiff.TrackedArray,Tracker.TrackedArray},Q3<:AbstractArray} #autodiff p
    if !phase.diffusionlimited && kfs === nothing
        kfs = getkfs(phase,T,P,C,ns,V,phi)
        if phi == 0.0
            krev = @fastmath kfs./getKcs(phase,T,Gs)
        else
            krev = @fastmath kfs./getKcs(phase,T,Gs,phi)
        end
    elseif !phase.diffusionlimited && !(kfs === nothing)
        if phi == 0.0
            krev = @fastmath kfs./getKcs(phase,T,Gs)
        else
            krev = @fastmath kfs./getKcs(phase,T,Gs,phi)
        end
    elseif phase.diffusionlimited && !(kfs === nothing)
        len = length(phase.reactions)
        krev = similar(kfs)
        kfsdiff = similar(kfs)
        @simd for i = 1:len
           @fastmath @inbounds kfsdiff[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;kf=kfs[i])
        end
        return kfsdiff, krev
    else
        len = length(phase.reactions)
        kfs = zeros(typeof(Gs[1]),len)
        krev = zeros(typeof(Gs[1]),len)
        @simd for i = 1:len
           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi)
        end
    end
    return kfs,krev
end

@inline function getkfkrevs(phase::U,T::W1,P::W2,C::W3,N::W4,ns::Q1,Gs::Array{Q2,1},diffs::Q3,V::W5,phi::W7;kfs::W6=nothing) where {U,W7,W6,W5<:Real,W1<:Real,W2<:Real,W3,W4,Q1<:AbstractArray,Q2<:ForwardDiff.Dual,Q3<:AbstractArray} #autodiff p
    if !phase.diffusionlimited && kfs === nothing
        kfs = getkfs(phase,T,P,C,ns,V,phi)
        if phi == 0.0
            krev = @fastmath kfs./getKcs(phase,T,Gs)
        else
            krev = @fastmath kfs./getKcs(phase,T,Gs,phi)
        end
    elseif !phase.diffusionlimited && !(kfs === nothing)
        if phi == 0.0
            krev = @fastmath kfs./getKcs(phase,T,Gs)
        else
            krev = @fastmath kfs./getKcs(phase,T,Gs,phi)
        end
    elseif phase.diffusionlimited && !(kfs === nothing)
        len = length(phase.reactions)
        krev = similar(kfs)
        @simd for i = 1:len
           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi;kf=kfs[i])
        end
    else
        len = length(phase.reactions)
        kfs = zeros(typeof(Gs[1]),len)
        krev = zeros(typeof(Gs[1]),len)
        @simd for i = 1:len
           @fastmath @inbounds kfs[i],krev[i] = getkfkrev(phase.reactions[i],phase,T,P,C,N,ns,Gs,diffs,V,phi)
        end
    end
    return kfs,krev
end

export getkfkrevs


text:
using Distributions, Plots; pyplot()

lambda, N = 1/3, 10^5
bulbs = [1,10,50]
xGrid = 0:0.1:10
C = [:blue :red :green]
dists = [Gamma(n,1/(n*lambda)) for n in bulbs]

function normalizedData(d::Gamma)
    sh = Int64(shape(d))
    data = [sum(-(1/(sh*lambda))*log.(rand(sh))) for _ in 1:N]
end

L = [ "Shape = "*string.(shape.(i))*", Scale = "*
    string.(round.(scale.(i),digits=2)) for i in dists ]

stephist( normalizedData.(dists), bins=50,
    normed=:true, c=C, xlims=(0,maximum(xGrid)),ylims=(0,1),
    xlabel="x", ylabel="Density", label="")
plot!(xGrid, [pdf.(i,xGrid) for i in dists], c=C, label=reshape(L, 1,:))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using MyPkgDemo
using Documenter

makedocs(;
    modules=[MyPkgDemo],
    authors="MegamindHenry",
    repo="https://github.com/MegamindHenry/MyPkgDemo.jl/blob/{commit}{path}#L{line}",
    sitename="MyPkgDemo.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://MegamindHenry.github.io/MyPkgDemo.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Add Nums" => "add_nums.md",
        "Cat Str" => "cat_str.md"
    ],
)

deploydocs(;
    repo="github.com/MegamindHenry/MyPkgDemo.jl",
)