text:
using Cairo
using ComplexVisual
@ComplexVisual.import_huge
using ComplexVisualGtk
@ComplexVisualGtk.import_huge

const fontface = cv_fontface("sans-serif")

function get_param_n()
    value_min, value_max = 1, 19
    slider_pos = CV_TranslateByOffset(Int)
    slider_pos.value = value_min
    set_slider_value_func = z -> begin
        slider_pos.value = round(UInt, max(min(z, value_max), value_min))
        return CV_Response(; redraw_flag=true)
    end
    label_style = cv_black → fontface → cv_fontsize(20)
    rulers=(CV_TickLabelAppearance(; label_style) ↦
                ("%.0f" ⇒ value_min:2:value_max),
            CV_Ruler("" ⇒  value_min+1:2:value_max),
            CV_TickLabelAppearance(; tick_length=0, gap=15, label_style) ↦
                (("n:" ⇒ 0.5),))

    return CV_ParamWithSlider(;
        slider_pos, value_min=value_min-0.1, value_max=value_max+0.1,
        rulers, set_slider_value_func)
end

param_n = get_param_n()

trafo1 = z -> 1/(1-z)
trafo2 = z -> abs(z) < 1 ?
    sum(k -> z^k, param_n.slider_pos.value:-1:1; init=1.0) :
    sum(k -> z^k, 1:param_n.slider_pos.value; init=1.0)

label_style = cv_black → fontface → cv_fontsize(20)
rulers=(CV_TickLabelAppearance(; label_style) ↦ ("%.0f" ⇒ [-1.0, 0.0, 1.0]),)
codomain1  = CV_Math2DCanvas(-1.5 + 1.5im, 1.5 - 1.5im, 150)
codomain2  = CV_Math2DCanvas(-1.5 + 1.5im, 1.5 - 1.5im, 150)

circle = (cv_black → cv_linewidth(3)) ↦ CV_LinePainter(
    cv_arc_lines(0.0, 2π, (1.0,)))

scene = cv_scene_comp_codomains_std((param_n, ), trafo1, trafo2,
    codomain1, codomain2; codomain1_re_rulers=rulers,
    painter1=CV_PortraitPainter(trafo1) → circle,
    painter2=CV_PortraitPainter(trafo2) → circle)
cv_get_redraw_func(scene)()

handler = cvg_visualize(scene)
cvg_wait_for_destroy(handler.window)

cvg_close(handler); cv_destroy(scene)

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
#
# Correctness Tests
#

fatalerrors = length(ARGS) > 0 && ARGS[1] == "-f"
anyerrors = false

using Base.Test
using Ook

my_tests = [
    "lexandparse.jl",
    "interpret.jl",
    "programs.jl"
]

println("Running tests:")

for my_test in my_tests
    try
        include(my_test)
        println("\t\033[1m\033[32mPASSED\033[0m: $(my_test)")
    catch
        anyerrors = true
        println("\t\033[1m\033[31mFAILED\033[0m: $(my_test)")
        if fatalerrors
            rethrow()
        end
    end
end

if anyerrors
    throw("Tests failed")
end

ext: jl


text:
#nb # %% slideshow [markdown] {"slideshow": {"slide_type": "slide"}}
#md # # [Lecture 2 - Types and Functions](@id lecture2)
#!md # # [Lecture 2 - Types and Functions]
#
# **Essential Skills For Computational Biology**
#
# *2020-06-08*
#
# Kevin Bonham, PhD

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "slide"}}
# ## A program is data and operations on data

#nb %% A slide [code] {"slideshow": {"slide_type": "fragment"}}
chromosomes = 23 # this is data

function diploid(n)
    return n * 2 # this is an operation
end
#-
#nb %% A slide [code] {"slideshow": {"slide_type": "fragment"}}
diploid(chromosomes) # this is a program

##nb # %% A slide [markdown] {"slideshow": {"slide_type": "subslide"}}
# ## How to run a julia program
#
# - From the REPL
# - From a script (Live demo)
# - From a file inside VS Code (Live demo)

#nb # %% A slide [markdown] {"slideshow": {"slide_type": "fragment"}}
# [Placeholder for video link](#)


text:
#####
##### Navier-Stokes and tracer advection equations
#####

""" Store previous value of the source term and calculate current source term. """
function calculate_interior_source_terms!(G, arch, grid, coriolis, buoyancy, surface_waves, closure, U, C, pHY′, K, F, parameters, time)
    # Manually choose thread-block layout here as it's ~20% faster.
    # See: https://github.com/climate-machine/Oceananigans.jl/pull/308
    Nx, Ny, Nz = grid.Nx, grid.Ny, grid.Nz

    if Nx == 1
        Tx, Ty = 1, min(256, Ny)
        Bx, By, Bz = Tx, floor(Int, Ny/Ty), Nz
    elseif Ny == 1
        Tx, Ty = min(256, Nx), 1
        Bx, By, Bz = floor(Int, Nx/Tx), Ty, Nz
    else
        Tx, Ty = 16, 16
        Bx, By, Bz = floor(Int, Nx/Tx), floor(Int, Ny/Ty), Nz
    end

    @launch(device(arch), threads=(Tx, Ty), blocks=(Bx, By, Bz),
            calculate_Gu!(G.u, grid, coriolis, surface_waves, closure, U, C, K, F, pHY′, parameters, time))

    @launch(device(arch), threads=(Tx, Ty), blocks=(Bx, By, Bz),
            calculate_Gv!(G.v, grid, coriolis, surface_waves, closure, U, C, K, F, pHY′, parameters, time))

    @launch(device(arch), threads=(Tx, Ty), blocks=(Bx, By, Bz),
            calculate_Gw!(G.w, grid, coriolis, surface_waves, closure, U, C, K, F, parameters, time))

    for tracer_index in 1:length(C)
        @inbounds Gc = G[tracer_index+3]
        @inbounds Fc = F[tracer_index+3]
        @inbounds  c = C[tracer_index]

        @launch(device(arch), threads=(Tx, Ty), blocks=(Bx, By, Bz),
                calculate_Gc!(Gc, grid, c, Val(tracer_index), closure, buoyancy, U, C, K, Fc, parameters, time))
    end

    return nothing
end

""" Calculate the right-hand-side of the u-momentum equation. """
function calculate_Gu!(Gu, grid, coriolis, surface_waves, closure, U, C, K, F, pHY′, parameters, time)
    @loop_xyz i j k grid begin
        @inbounds Gu[i, j, k] = ( - div_ũu(i, j, k, grid, U)
                                  - x_f_cross_U(i, j, k, grid, coriolis, U)
                                  - ∂xᶠᵃᵃ(i, j, k, grid, pHY′)
                                  + ∂ⱼ_2ν_Σ₁ⱼ(i, j, k, grid, closure, U, K)
                                  + x_curl_Uˢ_cross_U(i, j, k, grid, surface_waves, U, time)
                                  + ∂t_uˢ(i, j, k, grid, surface_waves, time)
                                  + F.u(i, j, k, grid, time, U, C, parameters))
    end
    return nothing
end

""" Calculate the right-hand-side of the v-momentum equation. """
function calculate_Gv!(Gv, grid, coriolis, surface_waves, closure, U, C, K, F, pHY′, parameters, time)
    @loop_xyz i j k grid begin
        @inbounds Gv[i, j, k] = ( - div_ũv(i, j, k, grid, U)
                                  - y_f_cross_U(i, j, k, grid, coriolis, U)
                                  - ∂yᵃᶠᵃ(i, j, k, grid, pHY′)
                                  + ∂ⱼ_2ν_Σ₂ⱼ(i, j, k, grid, closure, U, K)
                                  + y_curl_Uˢ_cross_U(i, j, k, grid, surface_waves, U, time)
                                  + ∂t_vˢ(i, j, k, grid, surface_waves, time)
                                  + F.v(i, j, k, grid, time, U, C, parameters))
    end
    return nothing
end

""" Calculate the right-hand-side of the w-momentum equation. """
function calculate_Gw!(Gw, grid, coriolis, surface_waves, closure, U, C, K, F, parameters, time)
    @loop_xyz i j k grid begin
        @inbounds Gw[i, j, k] = ( - div_ũw(i, j, k, grid, U)
                                  - z_f_cross_U(i, j, k, grid, coriolis, U)
                                  + ∂ⱼ_2ν_Σ₃ⱼ(i, j, k, grid, closure, U, K)
                                  + z_curl_Uˢ_cross_U(i, j, k, grid, surface_waves, U, time)
                                  + ∂t_wˢ(i, j, k, grid, surface_waves, time)
                                  + F.w(i, j, k, grid, time, U, C, parameters))
    end
    return nothing
end

""" Calculate the right-hand-side of the tracer advection-diffusion equation. """
function calculate_Gc!(Gc, grid, c, tracer_index, closure, buoyancy, U, C, K, Fc, parameters, time)
    @loop_xyz i j k grid begin
        @inbounds Gc[i, j, k] = (- div_uc(i, j, k, grid, U, c)
                                 + ∇_κ_∇c(i, j, k, grid, closure, c, tracer_index, K, C, buoyancy)
                                 + Fc(i, j, k, grid, time, U, C, parameters))
    end
    return nothing
end

""" Apply boundary conditions by adding flux divergences to the right-hand-side. """
function calculate_boundary_source_terms!(Gⁿ, arch, U, C, args...)

    # Velocity fields
    for i in 1:3
        apply_z_bcs!(Gⁿ[i], U[i], arch, args...)
        apply_y_bcs!(Gⁿ[i], U[i], arch, args...)
    end

    # Tracer fields
    for i in 4:length(Gⁿ)
        apply_z_bcs!(Gⁿ[i], C[i-3], arch, args...)
        apply_y_bcs!(Gⁿ[i], C[i-3], arch, args...)
    end

    return nothing
end

"""
Update the horizontal velocities u and v via

    `u^{n+1} = u^n + (Gu^{n+½} - δₓp_{NH} / Δx) Δt`

Note that the vertical velocity is not explicitly time stepped.
"""
function update_velocities!(U, grid, Δt, G, pNHS)
    @loop_xyz i j k grid begin
        @inbounds U.u[i, j, k] += (G.u[i, j, k] - ∂xᶠᵃᵃ(i, j, k, grid, pNHS)) * Δt
        @inbounds U.v[i, j, k] += (G.v[i, j, k] - ∂yᵃᶠᵃ(i, j, k, grid, pNHS)) * Δt
    end
    return nothing
end

"""
Update the horizontal velocities u and v via

    `u^{n+1} = u^n + Gu^{n+½} / Δt`

Note that the vertical velocity is not explicitly time stepped.
"""
function update_velocities!(U, grid, Δt, G, ::Nothing)
    @loop_xyz i j k grid begin
        @inbounds U.u[i, j, k] += G.u[i, j, k] * Δt
        @inbounds U.v[i, j, k] += G.v[i, j, k] * Δt
    end
    return nothing
end

"""
Update tracers via

    `c^{n+1} = c^n + Gc^{n+½} Δt`
"""
function update_tracer!(c, grid, Δt, Gc)
    @loop_xyz i j k grid begin
        @inbounds c[i, j, k] += Gc[i, j, k] * Δt
    end
    return nothing
end

"Update the solution variables (velocities and tracers)."
function update_solution!(U, C, arch, grid, Δt, G, pNHS)
    @launch device(arch) config=launch_config(grid, :xyz) update_velocities!(U, grid, Δt, G, pNHS)

    for i in 1:length(C)
        @inbounds c = C[i]
        @inbounds Gc = G[i+3]
        @launch device(arch) config=launch_config(grid, :xyz) update_tracer!(c, grid, Δt, Gc)
    end

    return nothing
end

#####
##### Vertical integrals
#####

"""
Update the hydrostatic pressure perturbation pHY′. This is done by integrating
the `buoyancy_perturbation` downwards:

    `pHY′ = ∫ buoyancy_perturbation dz` from `z=0` down to `z=-Lz`
"""
function update_hydrostatic_pressure!(pHY′, grid, buoyancy, C)
    @loop_xy i j grid begin
        @inbounds pHY′[i, j, grid.Nz] = - ℑzᵃᵃᶠ(i, j, grid.Nz+1, grid, buoyancy_perturbation, buoyancy, C) * ΔzF(i, j, grid.Nz+1, grid)
        @unroll for k in grid.Nz-1 : -1 : 1
            @inbounds pHY′[i, j, k] =
                pHY′[i, j, k+1] - ℑzᵃᵃᶠ(i, j, k+1, grid, buoyancy_perturbation, buoyancy, C) * ΔzF(i, j, k+1, grid)
        end
    end
    return nothing
end

"""
Compute the vertical velocity w by integrating the continuity equation from the bottom upwards

    `w^{n+1} = -∫ [∂/∂x (u^{n+1}) + ∂/∂y (v^{n+1})] dz`
"""
function compute_w_from_continuity!(model)
    @launch(device(model.architecture), config=launch_config(model.grid, :xy),
            _compute_w_from_continuity!(datatuple(model.velocities), model.grid))
    return nothing
end

function _compute_w_from_continuity!(U, grid)
    @loop_xy i j grid begin
        # U.w[i, j, 1] = 0 is enforced via halo regions.
        @unroll for k in 2:grid.Nz
            @inbounds U.w[i, j, k] = U.w[i, j, k-1] - ΔzC(i, j, k, grid) * hdivᶜᶜᵃ(i, j, k-1, grid, U.u, U.v)
        end
    end
    return nothing
end


text:
using Test, StatsModels
using DataFrames

using Unfold
include("test_utilities.jl")

data,evts = loadtestdata("test_case_3a") #
##

f_spl  = @formula 0~1+conditionA+spl(continuousA,4) # 1
f  = @formula 0~1+conditionA+continuousA # 1
data_r = reshape(data,(1,:))
data_e,times = Unfold.epoch(data=data_r,tbl=evts,τ=(-1.,1.),sfreq=10) # cut the data into epochs

m_mul,m_mul_results = fit(UnfoldLinearModel,f,evts,data_e,times)
m_mul_spl,m_mul_results_spl = fit(UnfoldLinearModel,f_spl,evts,data_e,times)

# asking for 4 splines should generate 4 splines
@test length(unique(m_mul_results_spl.term)) == 6 # XXX check back with Unfold whether this is the same! could be n-1 splines in Unfold. We should keep that comparable I guess

basisfunction = firbasis(τ=(-1,1),sfreq=10,name="A")
m_tul,m_tul_results = fit(UnfoldLinearModel,f,evts,data_r,basisfunction)
m_tul_spl,m_tul_results_spl = fit(UnfoldLinearModel,f_spl,evts,data_r,basisfunction)

evts_grid = gridexpand()

# results from timeexpanded and non should be equal
#yhat_tul  = predict(m_tul_spl,evts_grid)
#yhat_mul  = predict(m_mul_spl,evts_grid)
if 1 == 0
using AlgebraOfGraphics
yhat_mul.conditionA = categorical(yhat_mul.conditionA)
yhat_mul.continuousA = categorical(yhat_mul.continuousA)
m = mapping(:times,:yhat,color=:continuousA,linestyle=:conditionA)
df = yhat_mul
AlgebraOfGraphics.data(df) * visual(Lines) * m  |> draw
end

# test much higher number of splines
f_spl_many  = @formula 0~1+spl(continuousA,131) # 1
m_mul_spl_many,m_mul_results_spl_many = fit(UnfoldLinearModel,f_spl_many,evts,data_e,times)
@test length(unique(m_mul_results_spl_many.term)) == 132


text:
################################################################################
#
#   IMPLEMENTATIONS OF DIFFERENT PLAQUETTE RELATED FUNCTIONS
#
#   STRUCTURE OF THE FILE
#
#   1) FINDING PLAQUETTES IN LATTICE OBJECTS
#   - of a site
#   - of complete lattice
#
#   2) PRINTING PLAQUETTE STATISTICS
#   - of a site
#   - of complete lattice
#
#
#   BUG: Plaquettes dont implement periodic boundaries so far
#
################################################################################

#-----------------------------------------------------------------------------------------------------------------------------
#
#   Find Plaquettes of some site in the lattice with desired length len
#
#-----------------------------------------------------------------------------------------------------------------------------

"""
    getPlaquettesOfSite(lattice::Lattice, site::Int64, len::Int64)
    getPlaquettesOfSite(lattice::Lattice, site::Int64)

Function to find the plaquettes that loop through a given site with index `site`
of a passed `Lattice` object and have length `len`.

If no length is specified, but only a site, it is searched for the plaquettes with the
shortest length.

NOTE: Plaquettes cannot wrap around periodic boundaries correctly so far (in the sense that
they cannot visit sites with the same label in different periodic copies of the lattice)

The plaquettes are returned as a list of Arrays, where each element in the array is a site-index.

# Examples

```julia-repl
julia> plaquettes = getPlaquettesOfSite(lattice, 1, 6)
Array{Array{Int64, 1}}[
    [1, 5, 7, 22, 10, 17, 1],
    ...
]
```
"""
function getPlaquettesOfSite(lattice::Lattice, site::Int64, len::Int64)
    # get the connectivity of the lattice
    cl = getConnectionList(lattice)
    # history includes all sites that are already visited, including the current site
    function givePlaquettePaths(history::Array{Int64,1}, steps::Int64, patharray::Array{Array{Int64,1}, 1})
        # check if the maximum depth of the search has been reached
        if length(history) == steps+1
            # check if a closed loop exists for the maximum depth
            if history[end] == history[1]
                push!(patharray, history)
            end
            return
        elseif history[end] in history[1:end-1]
            # building a closed loop of not proper length
            return
        end
        # current site
        currentSite = history[end]
        # get all options
        #options = []
        #for c in cl[currentSite]
        #    push!(options, Int(c[2]))
        #end
        options = Int64[Int(c[2]) for c in cl[currentSite]]
        # distinguish between origin and further apart sites
        if length(history) == 1
            # go for all options
            for o in options
                # get the new history
                history_new = copy(history)
                push!(history_new, o)
                # get all subpaths
                givePlaquettePaths(history_new, steps, patharray)
            end
        else
            # go for all options that dont go back
            for o in options
                # check if going back
                if o == history[end-1]
                    continue
                end
                # get the new history
                history_new = copy(history)
                push!(history_new, o)
                # get all subpaths
                givePlaquettePaths(history_new, steps, patharray)
            end
        end
        # return
        return
    end
    # call the function to generate array of plaquettes
    paths = Array{Int64,1}[]
    givePlaquettePaths([site], len, paths)
    # check all paths if they are plaquettes
    plaquettes = Array{Int64,1}[]
    for p in paths
        if p[1] == p[end]
            # check if the inverse path is already in the plaquettes list
            if !(p[end:-1:1] in plaquettes)
                push!(plaquettes, p)
            end
        end
    end
    # return the array
    return plaquettes
end
function getPlaquettesOfSite(lattice::Lattice, site::Int64)
    # start with length = 3
    len = 3
    # search for plaquettes
    plaquettes = getPlaquettesOfSite(lattice, site, len)
    # iterate while there are no plaquettes found
    while length(plaquettes) == 0
        # increment the length
        len = len + 1
        # research for plaquettes
        plaquettes = getPlaquettesOfSite(lattice, site, len)
    end
    # return the plaquettes
    return plaquettes
end
export getPlaquettesOfSite

"""
    getPlaquettesOfLattice(lattice::Lattice, len::Int64)
    getPlaquettesOfLattice(lattice::Lattice)

Function to find all plaquettes which have length `len` of a given `Lattice` object.

If not length is specified, the plaquettes with shortest available length will be returned.

NOTE: Plaquettes cannot wrap around periodic boundaries correctly so far (in the sense that
they cannot visit sites with the same label in different periodic copies of the lattice)

The plaquettes are returned as a list of Arrays, where each element in the array is a site-index.

# Examples

```julia-repl
julia> plaquettes = getPlaquettesOfLattice(lattice, 6)
Array{Array{Int64, 1}}[
    [1, 5, 7, 22, 10, 17, 1],
    ...
]
```
"""
function getPlaquettesOfLattice(lattice::Lattice, len::Int64)
    # get all plaquettes of all sites
    plaquettes = Array{Int64,1}[]
    plaquettes_sorted = Array{Int64,1}[]
    for i in 1:size(lattice.positions,1)
        # obtain plaquettes
        plaquettes_site = getPlaquettesOfSite(lattice, i, len)
        # push plaquettes into array
        for p in plaquettes_site
            # check if present
            p_sorted = copy(p[1:end-1])
            sort!(p_sorted)
            # search if present
            found = false
            for ps in plaquettes_sorted
                # check if lists equal
                if findfirst([ps[i] == p_sorted[i] for i in 1:len], false) == 0
                    found = true
                    break
                end
            end
            if !(found)
                push!(plaquettes, p)
                push!(plaquettes_sorted, p_sorted)
            end
        end
    end
    # return the list
    return plaquettes
end
function getPlaquettesOfLattice(lattice::Lattice)
    # start with length = 3
    len = 3
    # search for plaquettes
    plaquettes = getPlaquettesOfLattice(lattice, len)
    # iterate while there are no plaquettes found
    while length(plaquettes) == 0
        # increment the length
        len = len + 1
        # research for plaquettes
        plaquettes = getPlaquettesOfLattice(lattice, len)
    end
    # return the plaquettes
    return plaquettes
end
export getPlaquettesOfLattice

#-----------------------------------------------------------------------------------------------------------------------------
#
#   Obtain and print information on plaquettes
#
#-----------------------------------------------------------------------------------------------------------------------------

"""
    printPlaquetteStatisticsOfSite(lattice::Lattice, site::Int64 [ ; detailed::Bool=false, max_length::Int64=12])

Function to find and print information about all plaquettes going through a given site of a passed `Lattice` object.
If `detailed=true`, the individual plaquettes are printed as well.
The option `max_length` specifies what maximal length of plaquettes should be detected.

NOTE: Plaquettes cannot wrap around periodic boundaries correctly so far (in the sense that
they cannot visit sites with the same label in different periodic copies of the lattice)

# Examples

```julia-repl
julia> plaquettes = printPlaquetteStatisticsOfSite(lattice, 1)
...

julia> plaquettes = printPlaquetteStatisticsOfSite(lattice, 1, detailed=true, max_length=10)
...
```
"""
function printPlaquetteStatisticsOfSite(lattice::Lattice, site::Int64; detailed::Bool=false, max_length::Int64=12)
    # print header information
    println("printing plaquette information for site $(site) of lattice:")
    # number of total plaquettes
    plaquettes_total = 0
    # iterate over all lengths
    for l in 3:max_length
        # get plaquettes
        plaquettes = getPlaquettesOfSite(lattice, site, l)
        # print statistics
        if detailed
            # definately print the number of plaquettes
            println("- length $(l) bonds: $(size(plaquettes,1)) plaquettes")
            # print more information
            for p in plaquettes
                println("    - $(p)")
            end
        else
            # maybe print the number of plaquettes
            if size(plaquettes,1) > 0
                println("- length $(l) bonds: $(size(plaquettes,1)) plaquettes")
            end
        end
        # add up how many plaquettes were found
        plaquettes_total += size(plaquettes,1)
    end
    # print total number
    println("In total: $(plaquettes_total) plaquettes found")
end
export printPlaquetteStatisticsOfSite

"""
    printPlaquetteStatisticsOfLattice(lattice::Lattice [ ; detailed::Bool=false, max_length::Int64=12])

Function to find and print information about all plaquettes of a passed `Lattice` object.
If `detailed=true`, the individual plaquettes are printed as well.
The option `max_length` specifies what maximal length of plaquettes should be detected.

NOTE: Plaquettes cannot wrap around periodic boundaries correctly so far (in the sense that
they cannot visit sites with the same label in different periodic copies of the lattice)

# Examples

```julia-repl
julia> plaquettes = printPlaquetteStatisticsOfLattice(lattice)
...

julia> plaquettes = printPlaquetteStatisticsOfLattice(lattice, detailed=true, max_length=10)
...
```
"""
function printPlaquetteStatisticsOfLattice(lattice::Lattice; detailed::Bool=false, max_length::Int64=12)
    # print header information
    println("printing plaquette information for lattice:")
    # number of total plaquettes
    plaquettes_total = 0
    # iterate over all lengths
    for l in 3:max_length
        # get plaquettes
        plaquettes = getPlaquettesOfLattice(lattice, l)
        # print statistics
        if detailed
            # definately print the number of plaquettes
            println("- length $(l) bonds: $(size(plaquettes,1)) plaquettes")
            # print more information
            for p in plaquettes
                println("    - $(p)")
            end
        else
            # maybe print the number of plaquettes
            if size(plaquettes,1) > 0
                println("- length $(l) bonds: $(size(plaquettes,1)) plaquettes")
            end
        end
        # add up how many plaquettes were found
        plaquettes_total += size(plaquettes,1)
    end
    # print total number
    println("In total: $(plaquettes_total) plaquettes found")
end
export printPlaquetteStatisticsOfLattice


text:
using GDAL
using Test

@testset "GDAL" begin

    # drivers
    # before being able to use any drivers, they must be registered first
    GDAL.gdalallregister()

    version = GDAL.gdalversioninfo("--version")
    n_gdal_driver = GDAL.gdalgetdrivercount()
    n_ogr_driver = GDAL.ogrgetdrivercount()
    @info """$version
    $n_gdal_driver GDAL drivers found
    $n_ogr_driver OGR drivers found
    """

    @test n_gdal_driver > 0
    @test n_ogr_driver > 0

    srs = GDAL.osrnewspatialreference(C_NULL)
    GDAL.osrimportfromepsg(srs, 4326) # fails if GDAL_DATA is not set correctly

    xmlnode_pointer = GDAL.cplparsexmlstring("hi")
    @test GDAL.cplgetxmlvalue(xmlnode_pointer, "b", "") == "hi"
    # load into Julia struct, mutate, and put back as Ref into GDAL
    xmlnode = unsafe_load(xmlnode_pointer)
    @test GDAL.cplserializexmltree(Ref(xmlnode)) == "\n  hi\n\n"
    GDAL.cpldestroyxmlnode(xmlnode_pointer)

    # ref https://github.com/JuliaGeo/GDAL.jl/pull/41#discussion_r143345433
    gfld = GDAL.ogr_gfld_create("name-a", GDAL.wkbPoint)
    @test gfld isa GDAL.OGRGeomFieldDefnH
    @test GDAL.ogr_gfld_getnameref(gfld) == "name-a"
    @test GDAL.ogr_gfld_gettype(gfld) == GDAL.wkbPoint
    # same as above but for the lower level C API
    gfld = GDAL.ogr_gfld_create("name-b", GDAL.wkbPolygon)
    @test gfld isa Ptr{GDAL.OGRGeomFieldDefnHS}
    @test GDAL.ogr_gfld_getnameref(gfld) == "name-b"
    @test GDAL.ogr_gfld_gettype(gfld) == GDAL.wkbPolygon

    cd(dirname(@__FILE__)) do
        rm("tmp", recursive = true, force = true)
        mkpath("tmp") # ensure it exists
        include("tutorial_raster.jl")
        include("tutorial_vector.jl")
        include("tutorial_vrt.jl")
        include("gdal_utils.jl")
        include("gdal_jll_utils.jl")
        include("drivers.jl")
        include("error.jl")
    end

    GDAL.gdaldestroydrivermanager()

end


text:
using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the linear advection equation

advectionvelocity = (0.2, -0.3)
equations = LinearScalarAdvectionEquation2D(advectionvelocity)

initial_condition = initial_condition_gauss

boundary_condition = BoundaryConditionDirichlet(initial_condition)
boundary_conditions = Dict(
  :all => boundary_condition
)

solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

# Deformed rectangle that looks like a waving flag,
# lower and upper faces are sinus curves, left and right are vertical lines.
f1(s) = SVector(-5.0, 5 * s - 5.0)
f2(s) = SVector( 5.0, 5 * s + 5.0)
f3(s) = SVector(5 * s, -5.0 + 5 * sin(0.5 * pi * s))
f4(s) = SVector(5 * s,  5.0 + 5 * sin(0.5 * pi * s))
faces = (f1, f2, f3, f4)

# This creates a mapping that transforms [-1, 1]^2 to the domain with the faces defined above.
# It generally doesn't work for meshes loaded from mesh files because these can be meshes
# of arbitrary domains, but the mesh below is specifically built on the domain [-1, 1]^2.
Trixi.validate_faces(faces)
mapping_flag = Trixi.transfinite_mapping(faces)

# Unstructured mesh with 24 cells of the square domain [-1, 1]^n
mesh_file = joinpath(@__DIR__, "square_unstructured_2.inp")
isfile(mesh_file) || download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp",
                              mesh_file)

mesh = P4estMesh{2}(mesh_file, polydeg=3,
                    mapping=mapping_flag,
                    initial_refinement_level=1)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 10.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
                                     extra_analysis_integrals=(entropy,))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_restart = SaveRestartCallback(interval=100,
                                   save_final_restart=true)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first),
                                      base_level=1,
                                      med_level=2, med_threshold=0.1,
                                      max_level=3, max_threshold=0.6)
amr_callback = AMRCallback(semi, amr_controller,
                           interval=5,
                           adapt_initial_condition=true,
                           adapt_initial_condition_only_refine=true)

stepsize_callback = StepsizeCallback(cfl=0.7)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_restart, save_solution,
                        amr_callback, stepsize_callback);

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);

summary_callback() # print the timer summary


text:
abstract type abstractDeflationFactor end

struct CustomDist{T}
	dist::T
end
@inline (cd::CustomDist)(u, v) = cd.dist(u, v)

"""
$(TYPEDEF)

This operator allows to handle the following situation. Assume you want to solve `F(x)=0` with a Newton algorithm but you want to avoid the process to return some already known solutions ``roots_i``. The deflation operator penalizes these roots ; the algorithm works very well despite its simplicity. You can create a `DeflationOperator` to define a scalar function `M(u)` used to find, with Newton iterations, the zeros of the following function
``F(u) \\cdot Π_i(\\|u - root_i\\|^{-2p} + \\alpha) := F(u) \\cdot M(u)`` where ``\\|u\\|^2 = dot(u, u)``. The fields of the struct `DeflationOperator` are as follows:

$(TYPEDFIELDS)

Given `defOp::DeflationOperator`, one can access its roots as `defOp[n]` as a shortcut for `defOp.roots[n]`. Note that you can also use `defOp[end]`.

Also, one can add (resp.remove) a new root by using `push!(defOp, newroot)` (resp. `pop!(defOp)`). Finally `length(defOp)` is a shortcut for `length(defOp.roots)`

## Constructors

- `DeflationOperator(p::Real, α::Real, roots::Vector{vectype}; autodiff = false)`
- `DeflationOperator(p::Real, dt, α::Real, roots::Vector{vectype}; autodiff = false)`
- `DeflationOperator(p::Real, α::Real, roots::Vector{vectype}, v::vectype; autodiff = false)`

The option `autodiff` triggers the use of automatic differentiation for the computation of the gradient of the scalar function `M`. This works only on `AbstractVector` for now.

## Custom distance

You are asked to pass a scalar product like `dot` to build a `DeflationOperator`. However, in some cases, you may want to pass a custom distance `dist(u,v)`. You can do this using

	`DeflationOperator(p, CustomDist(dist), α, roots)`

Note that passing `CustomDist(dist, true)` will trigger the use of automatic differentiation for the gradient of `M`.
"""
struct DeflationOperator{Tp <: Real, Tdot, T <: Real, vectype} <: abstractDeflationFactor
	"power `p`. You can use an `Int` for example"
	power::Tp

	"function, this function has to be bilinear and symmetric for the linear solver to work well"
	dot::Tdot

	"shift"
	α::T

	"roots"
	roots::Vector{vectype}

	# internal, to reduce allocations during computation
	tmp::vectype

	# internal, to reduce allocations during computation
	autodiff::Bool
end

# constructors
DeflationOperator(p::Real, α::Real, roots::Vector{vectype}; autodiff = false) where vectype = DeflationOperator(p, dot, α, roots, _copy(roots[1]), autodiff)
DeflationOperator(p::Real, dt, α::Real, roots::Vector{vectype}; autodiff = false) where vectype = DeflationOperator(p, dt, α, roots, _copy(roots[1]), autodiff)
DeflationOperator(p::Real, α::Real, roots::Vector{vectype}, v::vectype; autodiff = false) where vectype = DeflationOperator(p, dot, α, roots, v, autodiff)

# methods to deal with DeflationOperator
Base.eltype(df::DeflationOperator{Tp, Tdot, T, vectype}) where {Tp, Tdot, T, vectype} = T
Base.push!(df::DeflationOperator{Tp, Tdot, T, vectype}, v::vectype) where {Tp, Tdot, T, vectype} = push!(df.roots, v)
Base.pop!(df::DeflationOperator) = pop!(df.roots)
Base.getindex(df::DeflationOperator, inds...) = getindex(df.roots, inds...)
Base.length(df::DeflationOperator) = length(df.roots)
Base.deleteat!(df::DeflationOperator, id) = deleteat!(df.roots, id)
Base.empty!(df::DeflationOperator) = empty!(df.roots)
Base.firstindex(df::DeflationOperator) = 1
Base.lastindex(df::DeflationOperator) = length(df)
Base.copy(df::DeflationOperator) = DeflationOperator(df.power, df.dot, df.α, deepcopy(df.roots), copy(df.tmp), df.autodiff)

function Base.show(io::IO, df::DeflationOperator)
	println(io, "┌─ Deflation operator with ", length(df.roots)," root(s)")
	println(io, "├─ power    = ", df.power)
	println(io, "├─ α        = ", df.α)
	println(io, "├─ dist     = ", df.dot)
	println(io, "└─ autodiff = ", df.autodiff)
end

# Compute M(u)
# optimised version which does not allocate much
function (df::DeflationOperator{Tp, Tdot, T, vectype})(::Val{:inplace}, u, tmp) where {Tp, Tdot, T, vectype}
	length(df.roots) == 0 && return T(1)
	M(u) = T(1) / df.dot(u, u)^df.power + df.α
	# compute u - df.roots[1]
	copyto!(tmp, u); axpy!(T(-1), df.roots[1], tmp)
	out = M(tmp)
	for ii in 2:length(df.roots)
		copyto!(tmp, u); axpy!(T(-1), df.roots[ii], tmp)
		out *= M(tmp)
	end
	return out
end

# Compute M(u), efficient and do not allocate
(df::DeflationOperator{Tp, Tdot, T, vectype})(u::vectype) where {Tp, Tdot, T, vectype} = df(Val(:inplace), u, df.tmp)
(df::DeflationOperator{Tp, Tdot, T, vectype})(u) where {Tp, Tdot, T, vectype} = df(Val(:inplace), u, similar(u))

# version when a custom distance is passed
function (df::DeflationOperator{Tp, Tdot, T, vectype})(::Val{:inplace}, u, tmp) where {Tp, Tdot <: CustomDist, T, vectype}
	length(df.roots) == 0 && return T(1)
	M(u, v) = T(1) / df.dot(u, v)^df.power + df.α
	out = M(u, df.roots[1])
	for ii in 2:length(df.roots)
		out *= M(u, df.roots[ii])
	end
	return out
end

# Compute dM(u)⋅du. We use a tmp for storing
function (df::DeflationOperator{Tp, Tdot, T, vectype})(::Val{:dMwithTmp}, tmp, u, du) where {Tp, Tdot, T, vectype}
	length(df) == 0 && return T(0)
	if df.autodiff
		return ForwardDiff.derivative(t -> df(u .+ t .* du), 0)
	else
		δ = T(1e-8)
		copyto!(tmp, u); axpy!(δ, du, tmp)
		return (df(tmp) - df(u)) / δ
	end
end
(df::DeflationOperator)(u, du) = df(Val(:dMwithTmp), similar(u), u, du)

"""
	pb = DeflatedProblem(F, J, M::DeflationOperator)

This creates a deflated problem ``M(u) \\cdot F(u) = 0`` where `M` is a `DeflationOperator` which encodes the penalization term. `J` is the jacobian of `F`. Can be used to call `newton` and `continuation`.
"""
struct DeflatedProblem{Tp, Tdot, T, vectype, TF, TJ}
	F::TF
	J::TJ
	M::DeflationOperator{Tp, Tdot, T, vectype}
end
@inline length(prob::DeflatedProblem) = length(prob.M)

"""
Return the deflated function M(u) * F(u) where M(u) ∈ R
"""
function (df::DeflatedProblem{Tp, Tdot, T, vectype, TF, TJ})(u, par) where {Tp, Tdot, T, vectype, TF, TJ}
	out = df.F(u, par)
	rmul!(out, df.M(u))
	return out
end

"""
Return the jacobian of the deflated function M(u) * F(u) where M(u) ∈ R
"""
function (df::DeflatedProblem{Tp, Tdot, T, vectype, TF, TJ})(u::vectype, par, du) where {Tp, Tdot, T, vectype, TF, TJ}
	# dF(u)⋅du * M(u) + F(u) dM(u)⋅du
	# out = dF(u)⋅du * M(u)
	out = apply(df.J(u, par), du)
	M = df.M(u)
	rmul!(out, M)
	# we add the remaining part
	if length(df) > 0
		F = df.F(u, par)
		# F(u) dM(u)⋅du, out .+= df.M(u, du) .* F
		axpy!(df.M(u, du), F, out)
	end
	return out
end
###################################################################################################
# Implement the Jacobian operator of the Deflated problem

# this is used to define a custom linear solver
struct DeflatedLinearSolver <: AbstractLinearSolver end

"""
Implement the linear solver for the deflated problem.
"""
function (dfl::DeflatedLinearSolver)(J, rhs)
	# the expression of the Functional is now
	# F(u) * Π_i(dot(u - root_i, u - root_i)^{-power} + shift) := F(u) * M(u)
	# the expression of the differential is
	# dF(u)⋅du * M(u) + F(u) dM(u)⋅du

	# the point at which to compute the Jacobian
	u = J[1]
	p = J[2]

	# deflated Problem composite type
	defPb = J[3]
	linsolve = J[4]

	Fu = defPb.F(u, p)
	Mu = defPb.M(u)
	Ju = defPb.J(u, p)

	if length(defPb.M) == 0
		h1, _, it1 = linsolve(Ju, rhs)
		return h1, true, (it1, 0)
	end

	# linear solve for the deflated problem. We note that Mu ∈ R
	# hence dM(u)⋅du is a scalar. We now solve the following linear problem
	# M(u) * dF(u)⋅h + F(u) dM(u)⋅h = rhs
	h1, h2, _, (it1, it2) = linsolve(Ju, rhs, Fu)

	# We look for the expression of dM(u)⋅h
	# the solution is then h = Mu * h1 - z h2 where z has to be determined
	# z1 = dM(h)⋅h1
	tmp = similar(u)
	z1 = defPb.M(Val(:dMwithTmp), tmp, u, h1)

	# z2 = dM(h)⋅h2
	z2 = defPb.M(Val(:dMwithTmp), tmp, u, h2)

	z = z1 / (Mu + z2)

	# we extract the type of defPb
	_T = eltype(defPb.M)

	# return (h1 - z * h2) / Mu, true, (it1, it2)
	copyto!(tmp, h1)
	axpy!(-z, h2, tmp)
	rmul!(tmp, _T(1) / Mu)
	return tmp, true, (it1, it2)
end
####################################################################################################
"""
$(SIGNATURES)

This is the deflated version of the Krylov-Newton Solver for `F(x, p0) = 0` with Jacobian `J(x, p0)`. We refer to [`newton`](@ref) for more information. It penalises the roots saved in `defOp.roots`. The other arguments are as for `newton`. See [`DeflationOperator`](@ref) for more information.

# Arguments
Compared to [`newton`](@ref), the only different arguments are
- `defOp::DeflationOperator` deflation operator
- `linsolver` linear solver used to invert the Jacobian of the deflated functional.
    - We have a custom solver `DeflatedLinearSolver()` with requires solving two linear systems `J⋅x = rhs`.
    - For other linear solvers `<: AbstractLinearSolver`, a matrix free method is used for the deflated functional.
    - if passed `Val(:autodiff)`, then `ForwardDiff.jl` is used to compute the jacobian of the deflated problem

# Output:
- solution:
- history of residuals
- flag of convergence
- number of iterations

# Simplified call
When `J` is not passed. It then computed with finite differences. The call is as follows:

	newton(F, x0, p0, options, defOp; kwargs...)
"""
function newton(F, J, x0::vectype, p0, options::NewtonPar{T, L, E}, defOp::DeflationOperator{Tp, Tdot, T, vectype}, linsolver::DeflatedLinearSolver = DeflatedLinearSolver(); kwargs...) where {T, L, E, Tp, Tdot, vectype}
	# we create the new functional
	deflatedPb = DeflatedProblem(F, J, defOp)

	# and its jacobian
	Jacdf = (u, p) -> (u, p, deflatedPb, options.linsolver)

	# change the linear solver
	opt_def = @set options.linsolver = linsolver

	return newton(deflatedPb, Jacdf, x0, p0, opt_def; kwargs...)
end

function newton(F, J, x0::vectype, p0, options::NewtonPar{T, L, E}, defOp::DeflationOperator{Tp, Tdot, T, vectype}, linsolver::AbstractLinearSolver; kwargs...) where {Tp, T, Tdot, vectype, L, E}
	# we create the new functional
	deflatedPb = DeflatedProblem(F, J, defOp)

	# change the linear solver
	opt_def = @set options.linsolver = linsolver

	return newton(deflatedPb, (x, p) -> (dx -> deflatedPb(x, p, dx)), x0, p0, opt_def; kwargs...)
end

function newton(F, J, x0::vectype, p0, options::NewtonPar{T, L, E}, defOp::DeflationOperator{Tp, Tdot, T, vectype}, ::Val{:autodiff}; kwargs...) where {Tp, T, Tdot, vectype, L, E}
	# we create the new functional
	deflatedPb = DeflatedProblem(F, J, defOp)

	# we create the jacobian
	Jac = (x, p) -> ForwardDiff.jacobian(z -> deflatedPb(z, p), x)

	return newton(deflatedPb, Jac, x0, p0, options; kwargs...)
end

# simplified call when no Jacobian is given
function newton(F, x0::vectype, p0, options::NewtonPar{T, L, E}, defOp::DeflationOperator{Tp, Tdot, T, vectype}, linsolver = DeflatedLinearSolver(); kwargs...) where {Tp, T, Tdot, vectype, L, E}
	J = (u, p) -> finiteDifferences(z -> F(z,p), u)
	return newton(F, J, x0, p0, options, defOp, linsolver; kwargs...)
end

"""
$(TYPEDEF)

This specific Newton-Krylov method first tries to converge to a solution `sol0` close the guess `x0`. It then attempts to converge from the guess `x1` while avoiding the previous coonverged solution close to `sol0`. This is very handy for branch switching. The method is based on a deflated Newton-Krylov solver.
"""
function newton(F, J, x0::vectype, x1::vectype, p0, options::NewtonPar{T, L, E}, defOp::DeflationOperator = DeflationOperator(2, dot, 1.0, Vector{vectype}(), _copy(x0), true), linsolver = DeflatedLinearSolver(); kwargs...) where {T, vectype, L, E}
	res0 = newton(F, J, x0, p0, options; kwargs...)
	@assert res0[3] "Newton did not converge to the trivial solution x0."
	push!(defOp, res0[1])
	res1 = newton(F, J, x1, p0, (@set options.maxIter = 10options.maxIter), defOp; kwargs...)
	@assert res1[3] "Deflated Newton did not converge to the non-trivial solution ( i.e. on the bifurcated branch)."
	return res1, res0
end


text:
# modified by N Laohakunakorn
# University of Edinburgh 2020

using Statistics
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2017 Varnerlab
# School of Chemical Engineering Purdue University
# W. Lafayette IN 46907 USA

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #
function Control(t,x,rate_vector,data_dictionary)
# ---------------------------------------------------------------------- #
# Control.jl was generated using the Kwatee code generation system.
# Username: nicholas
# Type: CFPS-JULIA
# Version: 1.0
# Generation timestamp: 03-24-2017 13:01:36
#
# Arguments:
# t  - current time
# x  - state vector
# rate_vector - vector of reaction rates
# data_dictionary  - Data dictionary instance (holds model parameters)
# ---------------------------------------------------------------------- #

# Set a default value for the allosteric control variables -
number_of_reactions = length(rate_vector);
control_vector = ones(number_of_reactions);
control_parameter_array = data_dictionary["CONTROL_PARAMETER_ARRAY"];

# Alias the species vector -
M_g6p_c = x[1];
M_f6p_c = x[2];
M_fdp_c = x[3];
M_dhap_c = x[4];
M_g3p_c = x[5];
M_13dpg_c = x[6];
M_3pg_c = x[7];
M_2pg_c = x[8];
M_oaa_c = x[9];
M_coa_c = x[10];
M_accoa_c = x[11];
M_6pgl_c = x[12];
M_6pgc_c = x[13];
M_ru5p_D_c = x[14];
M_xu5p_D_c = x[15];
M_r5p_c = x[16];
M_s7p_c = x[17];
M_e4p_c = x[18];
M_2ddg6p_c = x[19];
M_cit_c = x[20];
M_icit_c = x[21];
M_akg_c = x[22];
M_succoa_c = x[23];
M_q8_c = x[24];
M_fum_c = x[25];
M_q8h2_c = x[26];
M_mql8_c = x[27];
M_mqn8_c = x[28];
M_h_e = x[29];
M_ppi_c = x[30];
M_glx_c = x[31];
M_actp_c = x[32];
M_etoh_c = x[33];
M_for_c = x[34];
M_nh3_c = x[35];
M_arg_L_c = x[36];
M_h2s_c = x[37];
M_thf_c = x[38];
M_mlthf_c = x[39];
M_aicar_c = x[40];
M_5mthf_c = x[41];
M_chor_c = x[42];
M_h2o2_c = x[43];
M_mglx_c = x[44];
M_prop_c = x[45];
M_indole_c = x[46];
M_cadav_c = x[47];
M_gaba_c = x[48];
M_glycoA_c = x[49];
M_78mdp_c = x[50];
M_4adochor_c = x[51];
M_4abz_c = x[52];
M_78dhf_c = x[53];
M_dhf_c = x[54];
M_methf_c = x[55];
M_10fthf_c = x[56];
M_prpp_c = x[57];
M_hco3_c = x[58];
M_clasp_c = x[59];
M_or_c = x[60];
M_omp_c = x[61];
M_5pbdra = x[62];
M_gar_c = x[63];
M_fgar_c = x[64];
M_fgam_c = x[65];
M_air_c = x[66];
M_cair_c = x[67];
M_saicar_c = x[68];
M_faicar_c = x[69];
M_imp_c = x[70];
M_xmp_c = x[71];
GENE_CAT = x[72];
RNAP = x[73];
OPEN_GENE_CAT = x[74];
mRNA_CAT = x[75];
RIBOSOME = x[76];
RIBOSOME_START_CAT = x[77];
M_ala_L_c_tRNA = x[78];
M_arg_L_c_tRNA = x[79];
M_asn_L_c_tRNA = x[80];
M_asp_L_c_tRNA = x[81];
M_cys_L_c_tRNA = x[82];
M_glu_L_c_tRNA = x[83];
M_gln_L_c_tRNA = x[84];
M_gly_L_c_tRNA = x[85];
M_his_L_c_tRNA = x[86];
M_ile_L_c_tRNA = x[87];
M_leu_L_c_tRNA = x[88];
M_lys_L_c_tRNA = x[89];
M_met_L_c_tRNA = x[90];
M_phe_L_c_tRNA = x[91];
M_pro_L_c_tRNA = x[92];
M_ser_L_c_tRNA = x[93];
M_thr_L_c_tRNA = x[94];
M_trp_L_c_tRNA = x[95];
M_tyr_L_c_tRNA = x[96];
M_val_L_c_tRNA = x[97];
PROTEIN_CAT = x[98];
tRNA = x[99];
M_glc_D_c = x[100];
M_pep_c = x[101];
M_pyr_c = x[102];
M_ac_c = x[103];
M_lac_D_c = x[104];
M_mal_L_c = x[105];
M_atp_c = x[106];
M_adp_c = x[107];
M_amp_c = x[108];
M_gtp_c = x[109];
M_gdp_c = x[110];
M_gmp_c = x[111];
M_utp_c = x[112];
M_udp_c = x[113];
M_ump_c = x[114];
M_ctp_c = x[115];
M_cdp_c = x[116];
M_cmp_c = x[117];
M_succ_c = x[118];
M_asp_L_c = x[119];
M_gly_L_c = x[120];
M_ile_L_c = x[121];
M_asn_L_c = x[122];
M_cys_L_c = x[123];
M_lys_L_c = x[124];
M_his_L_c = x[125];
M_ala_L_c = x[126];
M_phe_L_c = x[127];
M_pro_L_c = x[128];
M_ser_L_c = x[129];
M_thr_L_c = x[130];
M_trp_L_c = x[131];
M_tyr_L_c = x[132];
M_val_L_c = x[133];
M_met_L_c = x[134];
M_leu_L_c = x[135];
M_glu_L_c = x[136];
M_gln_L_c = x[137];
M_o2_c = x[138];
M_co2_c = x[139];
M_pi_c = x[140];
M_nh4_c = x[141];
M_so4_c = x[142];
M_h_c = x[143];
M_h2o_c = x[144];
M_nad_c = x[145];
M_nadh_c = x[146];
M_nadp_c = x[147];
M_nadph_c = x[148];
E_R_glk_atp = x[149];
E_R_pgi = x[150];
E_R_pgi_R = x[151];
E_R_pfk = x[152];
E_R_fdp = x[153];
E_R_fbaA = x[154];
E_R_fbaA_R = x[155];
E_R_tpiA = x[156];
E_R_tpiA_R = x[157];
E_R_gapA = x[158];
E_R_gapA_R = x[159];
E_R_pgk = x[160];
E_R_pgk_R = x[161];
E_R_gpm = x[162];
E_R_gpm_R = x[163];
E_R_eno = x[164];
E_R_eno_R = x[165];
E_R_pyk = x[166];
E_R_pck = x[167];
E_R_ppc = x[168];
E_R_pdh = x[169];
E_R_pps = x[170];
E_R_zwf = x[171];
E_R_zwf_R = x[172];
E_R_pgl = x[173];
E_R_gnd = x[174];
E_R_rpe = x[175];
E_R_rpe_R = x[176];
E_R_rpi = x[177];
E_R_rpi_R = x[178];
E_R_talAB = x[179];
E_R_talAB_R = x[180];
E_R_tkt1 = x[181];
E_R_tkt1_R = x[182];
E_R_tkt2 = x[183];
E_R_tkt2_R = x[184];
E_R_edd = x[185];
E_R_eda = x[186];
E_R_gltA = x[187];
E_R_acn = x[188];
E_R_acn_R = x[189];
E_R_icd = x[190];
E_R_icd_R = x[191];
E_R_sucAB = x[192];
E_R_sucCD = x[193];
E_R_sdh = x[194];
E_R_frd = x[195];
E_R_fum = x[196];
E_R_fum_R = x[197];
E_R_mdh = x[198];
E_R_mdh_R = x[199];
E_R_cyd = x[200];
E_R_cyo = x[201];
E_R_app = x[202];
E_R_atp = x[203];
E_R_nuo = x[204];
E_R_pnt1 = x[205];
E_R_pnt2 = x[206];
E_R_ndh1 = x[207];
E_R_ndh2 = x[208];
E_R_ppa = x[209];
E_R_aceA = x[210];
E_R_aceB = x[211];
E_R_maeA = x[212];
E_R_maeB = x[213];
E_R_pta = x[214];
E_R_pta_R = x[215];
E_R_ackA = x[216];
E_R_ackA_R = x[217];
E_R_acs = x[218];
E_R_adhE = x[219];
E_R_adhE_R = x[220];
E_R_ldh = x[221];
E_R_ldh_R = x[222];
E_R_pflAB = x[223];
E_R_alaAC = x[224];
E_R_alaAC_R = x[225];
E_R_arg = x[226];
E_R_aspC = x[227];
E_R_asnB = x[228];
E_R_asnA = x[229];
E_R_cysEMK = x[230];
E_R_gltBD = x[231];
E_R_gdhA = x[232];
E_R_gdhA_R = x[233];
E_R_glnA = x[234];
E_R_glyA = x[235];
E_R_his = x[236];
E_R_ile = x[237];
E_R_leu = x[238];
E_R_lys = x[239];
E_R_met = x[240];
E_R_phe = x[241];
E_R_pro = x[242];
E_R_serABC = x[243];
E_R_thr = x[244];
E_R_trp = x[245];
E_R_tyr = x[246];
E_R_val = x[247];
E_R_arg_deg = x[248];
E_R_asp_deg = x[249];
E_R_asn_deg = x[250];
E_R_gly_deg = x[251];
E_R_mglx_deg = x[252];
E_R_ser_deg = x[253];
E_R_pro_deg = x[254];
E_R_thr_deg1 = x[255];
E_R_thr_deg2 = x[256];
E_R_thr_deg3 = x[257];
E_R_trp_deg = x[258];
E_R_cys_deg = x[259];
E_R_lys_deg = x[260];
E_R_gln_deg = x[261];
E_R_glu_deg = x[262];
E_R_gaba_deg1 = x[263];
E_R_gaba_deg2 = x[264];
E_R_chor = x[265];
E_R_fol_e = x[266];
E_R_fol_1 = x[267];
E_R_fol_2a = x[268];
E_R_fol_2b = x[269];
E_R_fol_3 = x[270];
E_R_fol_4 = x[271];
E_R_gly_fol = x[272];
E_R_gly_fol_R = x[273];
E_R_mthfd = x[274];
E_R_mthfd_R = x[275];
E_R_mthfc = x[276];
E_R_mthfc_R = x[277];
E_R_mthfr2a = x[278];
E_R_mthfr2b = x[279];
E_R_prpp_syn = x[280];
E_R_or_syn_1 = x[281];
E_R_or_syn_2 = x[282];
E_R_omp_syn = x[283];
E_R_ump_syn = x[284];
E_R_ctp_1 = x[285];
E_R_ctp_2 = x[286];
E_R_A_syn_1 = x[287];
E_R_A_syn_2 = x[288];
E_R_A_syn_3 = x[289];
E_R_A_syn_4 = x[290];
E_R_A_syn_5 = x[291];
E_R_A_syn_6 = x[292];
E_R_A_syn_7 = x[293];
E_R_A_syn_8 = x[294];
E_R_A_syn_9 = x[295];
E_R_A_syn_10 = x[296];
E_R_A_syn_12 = x[297];
E_R_xmp_syn = x[298];
E_R_gmp_syn = x[299];
E_R_atp_amp = x[300];
E_R_utp_ump = x[301];
E_R_ctp_cmp = x[302];
E_R_gtp_gmp = x[303];
E_R_atp_adp = x[304];
E_R_utp_adp = x[305];
E_R_ctp_adp = x[306];
E_R_gtp_adp = x[307];
E_R_udp_utp = x[308];
E_R_cdp_ctp = x[309];
E_R_gdp_gtp = x[310];
E_R_atp_ump = x[311];
E_R_atp_cmp = x[312];
E_R_atp_gmp = x[313];
E_R_adk_atp = x[314];
E_Import_o2 = x[315];
E_Import_co2 = x[316];
E_Import_pi = x[317];
E_Import_nh4 = x[318];
E_Import_so4 = x[319];
E_Import_h2o = x[320];
E_Export_o2 = x[321];
E_Export_co2 = x[322];
E_Export_pi = x[323];
E_Export_nh4 = x[324];
E_Export_so4 = x[325];
E_Export_h2o = x[326];
E_Proton_gradient = x[327];
E_transcriptional_initiation_CAT = x[328];
E_transcription_CAT = x[329];
E_mRNA_degradation_CAT = x[330];
E_translation_initiation_CAT = x[331];
E_translation_CAT = x[332];
E_tRNA_charging_M_ala_L_c_CAT = x[333];
E_tRNA_charging_M_arg_L_c_CAT = x[334];
E_tRNA_charging_M_asn_L_c_CAT = x[335];
E_tRNA_charging_M_asp_L_c_CAT = x[336];
E_tRNA_charging_M_cys_L_c_CAT = x[337];
E_tRNA_charging_M_glu_L_c_CAT = x[338];
E_tRNA_charging_M_gln_L_c_CAT = x[339];
E_tRNA_charging_M_gly_L_c_CAT = x[340];
E_tRNA_charging_M_his_L_c_CAT = x[341];
E_tRNA_charging_M_ile_L_c_CAT = x[342];
E_tRNA_charging_M_leu_L_c_CAT = x[343];
E_tRNA_charging_M_lys_L_c_CAT = x[344];
E_tRNA_charging_M_met_L_c_CAT = x[345];
E_tRNA_charging_M_phe_L_c_CAT = x[346];
E_tRNA_charging_M_pro_L_c_CAT = x[347];
E_tRNA_charging_M_ser_L_c_CAT = x[348];
E_tRNA_charging_M_thr_L_c_CAT = x[349];
E_tRNA_charging_M_trp_L_c_CAT = x[350];
E_tRNA_charging_M_tyr_L_c_CAT = x[351];
E_tRNA_charging_M_val_L_c_CAT = x[352];

transfer_function_vector = Float64[];
# type: inhibition actor: M_pep_c target: R_pfk
push!(transfer_function_vector,1.0 - control_parameter_array[1,1]*M_pep_c^control_parameter_array[1,2]/(1+control_parameter_array[1,1]*M_pep_c^control_parameter_array[1,2]));
control_vector[4] = mean(transfer_function_vector);
transfer_function_vector = 0;

transfer_function_vector = Float64[];
# type: activation actor: M_pep_c target: R_fdp
push!(transfer_function_vector,control_parameter_array[2,1]*M_pep_c^control_parameter_array[2,2]/(1+control_parameter_array[2,1]*M_pep_c^control_parameter_array[2,2]));
control_vector[5] = mean(transfer_function_vector);
transfer_function_vector = 0;

transfer_function_vector = Float64[];
# type: activation actor: M_fdp_c target: R_pyk
push!(transfer_function_vector,control_parameter_array[3,1]*M_fdp_c^control_parameter_array[3,2]/(1+control_parameter_array[3,1]*M_fdp_c^control_parameter_array[3,2]));
control_vector[18] = mean(transfer_function_vector);
transfer_function_vector = 0;

transfer_function_vector = Float64[];
# type: inhibition actor: M_pep_c target: R_pck
push!(transfer_function_vector,1.0 - control_parameter_array[4,1]*M_pep_c^control_parameter_array[4,2]/(1+control_parameter_array[4,1]*M_pep_c^control_parameter_array[4,2]));
control_vector[19] = mean(transfer_function_vector);
transfer_function_vector = 0;

transfer_function_vector = Float64[];
# type: activation actor: M_fdp_c target: R_ppc
push!(transfer_function_vector,control_parameter_array[5,1]*M_fdp_c^control_parameter_array[5,2]/(1+control_parameter_array[5,1]*M_fdp_c^control_parameter_array[5,2]));
control_vector[20] = mean(transfer_function_vector);
transfer_function_vector = 0;

transfer_function_vector = Float64[];
# type: inhibition actor: M_pyr_c target: R_pdh
push!(transfer_function_vector,1.0 - control_parameter_array[6,1]*M_pyr_c^control_parameter_array[6,2]/(1+control_parameter_array[6,1]*M_pyr_c^control_parameter_array[6,2]));
control_vector[21] = mean(transfer_function_vector);
transfer_function_vector = 0;

transfer_function_vector = Float64[];
# type: inhibition actor: M_pep_c target: R_pps
push!(transfer_function_vector,1.0 - control_parameter_array[7,1]*M_pep_c^control_parameter_array[7,2]/(1+control_parameter_array[7,1]*M_pep_c^control_parameter_array[7,2]));
control_vector[22] = mean(transfer_function_vector);
transfer_function_vector = 0;

transfer_function_vector = Float64[];
# type: inhibition actor: M_akg_c target: R_gltA
push!(transfer_function_vector,1.0 - control_parameter_array[8,1]*M_akg_c^control_parameter_array[8,2]/(1+control_parameter_array[8,1]*M_akg_c^control_parameter_array[8,2]));
control_vector[39] = mean(transfer_function_vector);
transfer_function_vector = 0;

transfer_function_vector = Float64[];
# type: inhibition actor: M_pep_c target: R_icd
push!(transfer_function_vector,1.0 - control_parameter_array[9,1]*M_pep_c^control_parameter_array[9,2]/(1+control_parameter_array[9,1]*M_pep_c^control_parameter_array[9,2]));
control_vector[42] = mean(transfer_function_vector);
transfer_function_vector = 0;

transfer_function_vector = Float64[];
# type: inhibition actor: M_akg_c target: R_aceA
push!(transfer_function_vector,1.0 - control_parameter_array[10,1]*M_akg_c^control_parameter_array[10,2]/(1+control_parameter_array[10,1]*M_akg_c^control_parameter_array[10,2]));
# type: inhibition actor: M_pep_c target: R_aceA
push!(transfer_function_vector,1.0 - control_parameter_array[11,1]*M_pep_c^control_parameter_array[11,2]/(1+control_parameter_array[11,1]*M_pep_c^control_parameter_array[11,2]));
# type: inhibition actor: M_3pg_c target: R_aceA
push!(transfer_function_vector,1.0 - control_parameter_array[12,1]*M_3pg_c^control_parameter_array[12,2]/(1+control_parameter_array[12,1]*M_3pg_c^control_parameter_array[12,2]));
control_vector[62] = mean(transfer_function_vector);
transfer_function_vector = 0;

transfer_function_vector = Float64[];
# type: inhibition actor: M_akg_c target: R_aceB
push!(transfer_function_vector,1.0 - control_parameter_array[13,1]*M_akg_c^control_parameter_array[13,2]/(1+control_parameter_array[13,1]*M_akg_c^control_parameter_array[13,2]));
# type: inhibition actor: M_pep_c target: R_aceB
push!(transfer_function_vector,1.0 - control_parameter_array[14,1]*M_pep_c^control_parameter_array[14,2]/(1+control_parameter_array[14,1]*M_pep_c^control_parameter_array[14,2]));
# type: inhibition actor: M_3pg_c target: R_aceB
push!(transfer_function_vector,1.0 - control_parameter_array[15,1]*M_3pg_c^control_parameter_array[15,2]/(1+control_parameter_array[15,1]*M_3pg_c^control_parameter_array[15,2]));
control_vector[63] = mean(transfer_function_vector);
transfer_function_vector = 0;

transfer_function_vector = Float64[];
# type: inhibition actor: M_accoa_c target: R_maeB
push!(transfer_function_vector,1.0 - control_parameter_array[16,1]*M_accoa_c^control_parameter_array[16,2]/(1+control_parameter_array[16,1]*M_accoa_c^control_parameter_array[16,2]));
control_vector[65] = mean(transfer_function_vector);
transfer_function_vector = 0;

transfer_function_vector = Float64[];
# type: activation actor: M_pyr_c target: R_ldh_R
push!(transfer_function_vector,control_parameter_array[17,1]*M_pyr_c^control_parameter_array[17,2]/(1+control_parameter_array[17,1]*M_pyr_c^control_parameter_array[17,2]));
control_vector[74] = mean(transfer_function_vector);
transfer_function_vector = 0;

# Modify the rate_vector with the control variables -
rate_vector = rate_vector.*control_vector;

# Return the modified rate vector -
return rate_vector;
end

max_stars_repo_name: Laohakunakorn-Group/nllab-CFPSmetabolic

max_issues_repo_name: Laohakunakorn-Group/nllab-CFPSmetabolic

max_forks_repo_name: Laohakunakorn-Group/nllab-CFPSmetabolic


text:
abstract type Basis{T} <: LazyQuasiMatrix{T} end
abstract type Weight{T} <: LazyQuasiVector{T} end

const WeightedBasis{T, A<:AbstractQuasiVector, B<:Basis} = BroadcastQuasiMatrix{T,typeof(*),<:Tuple{A,B}}

struct WeightLayout <: AbstractQuasiLazyLayout end
abstract type AbstractBasisLayout <: AbstractQuasiLazyLayout end
struct BasisLayout <: AbstractBasisLayout end
struct SubBasisLayout <: AbstractBasisLayout end
struct MappedBasisLayout <: AbstractBasisLayout end
struct WeightedBasisLayout <: AbstractBasisLayout end
struct SubWeightedBasisLayout <: AbstractBasisLayout end
struct MappedWeightedBasisLayout <: AbstractBasisLayout end

SubBasisLayouts = Union{SubBasisLayout,SubWeightedBasisLayout}
WeightedBasisLayouts = Union{WeightedBasisLayout,SubWeightedBasisLayout,MappedWeightedBasisLayout}
MappedBasisLayouts = Union{MappedBasisLayout,MappedWeightedBasisLayout}

abstract type AbstractAdjointBasisLayout <: AbstractQuasiLazyLayout end
struct AdjointBasisLayout <: AbstractAdjointBasisLayout end
struct AdjointSubBasisLayout <: AbstractAdjointBasisLayout end
struct AdjointMappedBasisLayout <: AbstractAdjointBasisLayout end

MemoryLayout(::Type{<:Basis}) = BasisLayout()
MemoryLayout(::Type{<:Weight}) = WeightLayout()

adjointlayout(::Type, ::AbstractBasisLayout) = AdjointBasisLayout()
adjointlayout(::Type, ::SubBasisLayout) = AdjointSubBasisLayout()
adjointlayout(::Type, ::MappedBasisLayouts) = AdjointMappedBasisLayout()
broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::AbstractBasisLayout) = WeightedBasisLayout()
broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::SubBasisLayout) = WeightedBasisLayout()
broadcastlayout(::Type{typeof(*)}, ::WeightLayout, ::MappedBasisLayouts) = MappedWeightedBasisLayout()

# A sub of a weight is still a weight
sublayout(::WeightLayout, _) = WeightLayout()

## Weighted basis interface
unweightedbasis(P::BroadcastQuasiMatrix{<:Any,typeof(*),<:Tuple{AbstractQuasiVector,AbstractQuasiMatrix}}) = last(P.args)
unweightedbasis(V::SubQuasiArray) = view(unweightedbasis(parent(V)), parentindices(V)...)

# Default is lazy
ApplyStyle(::typeof(pinv), ::Type{<:Basis}) = LazyQuasiArrayApplyStyle()
pinv(J::Basis) = apply(pinv,J)

function ==(A::Basis, B::Basis)
    axes(A) == axes(B) && throw(ArgumentError("Override == to compare bases of type $(typeof(A)) and $(typeof(B))"))
    false
end

@inline copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(+)}}) = +(broadcast(\,Ref(L.A),arguments(L.B))...)
@inline copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(+)},<:Any,<:AbstractQuasiVector}) =
    transform_ldiv(L.A, L.B)

@inline function copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(-)}})
    a,b = arguments(L.B)
    (L.A\a)-(L.A\b)
end

@inline copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(-)},<:Any,<:AbstractQuasiVector}) =
    transform_ldiv(L.A, L.B)

@inline function copy(P::Ldiv{<:AbstractBasisLayout,<:AbstractBasisLayout})
    A, B = P.A, P.B
    A == B || throw(ArgumentError("Override copy for $(typeof(A)) \\ $(typeof(B))"))
    SquareEye{eltype(P)}((axes(A,2),))
end
@inline function copy(P::Ldiv{<:SubBasisLayouts,<:SubBasisLayouts})
    A, B = P.A, P.B
    parent(A) == parent(B) ||
        throw(ArgumentError("Override copy for $(typeof(A)) \\ $(typeof(B))"))
    Eye{eltype(P)}((axes(A,2),axes(B,2)))
end

@inline function copy(P::Ldiv{<:MappedBasisLayouts,<:MappedBasisLayouts})
    A, B = P.A, P.B
    demap(A)\demap(B)
end

@inline copy(L::Ldiv{<:AbstractBasisLayout,<:SubBasisLayouts}) = apply(\, L.A, ApplyQuasiArray(L.B))
@inline function copy(L::Ldiv{<:SubBasisLayouts,<:AbstractBasisLayout})
    P = parent(L.A)
    kr, jr = parentindices(L.A)
    layout_getindex(apply(\, P, L.B), jr, :) # avoid sparse arrays
end

for Bas1 in (:Basis, :WeightedBasis), Bas2 in (:Basis, :WeightedBasis)
    @eval ==(A::SubQuasiArray{<:Any,2,<:$Bas1}, B::SubQuasiArray{<:Any,2,<:$Bas2}) =
        parentindices(A) == parentindices(B) && parent(A) == parent(B)
end

# multiplication operators, reexpand in basis A
@inline function _broadcast_mul_ldiv(::Tuple{Any,AbstractBasisLayout}, A, B)
    a,b = arguments(B)
    @assert a isa AbstractQuasiVector # Only works for vec .* mat
    ab = (A * (A \ a)) .* b # broadcasted should be overloaded
    MemoryLayout(ab) isa BroadcastLayout && error("Overload broadcasted(_, ::typeof(*), ::$(typeof(ab.args[1])), ::$(typeof(b)))")
    A \ ab
end

_broadcast_mul_ldiv(_, A, B) = copy(Ldiv{typeof(MemoryLayout(A)),UnknownLayout}(A,B))

copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(*)}}) = _broadcast_mul_ldiv(map(MemoryLayout,arguments(L.B)), L.A, L.B)
copy(L::Ldiv{<:AbstractBasisLayout,BroadcastLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = _broadcast_mul_ldiv(map(MemoryLayout,arguments(L.B)), L.A, L.B)

# expansion
_grid(_, P) = error("Overload Grid")
_grid(::MappedBasisLayout, P) = igetindex.(Ref(parentindices(P)[1]), grid(demap(P)))
_grid(::SubBasisLayout, P) = grid(parent(P))
_grid(::WeightedBasisLayouts, P) = grid(unweightedbasis(P))
grid(P) = _grid(MemoryLayout(typeof(P)), P)

struct TransformFactorization{T,Grid,Plan,IPlan} <: Factorization{T}
    grid::Grid
    plan::Plan
    iplan::IPlan
end

TransformFactorization(grid, plan) =
    TransformFactorization{promote_type(eltype(grid),eltype(plan)),typeof(grid),typeof(plan),Nothing}(grid, plan, nothing)

TransformFactorization(grid, ::Nothing, iplan) =
    TransformFactorization{promote_type(eltype(grid),eltype(iplan)),typeof(grid),Nothing,typeof(iplan)}(grid, nothing, iplan)

grid(T::TransformFactorization) = T.grid

\(a::TransformFactorization{<:Any,<:Any,Nothing}, b::AbstractQuasiVector{T}) where T = a.iplan \  convert(Array{T}, b[a.grid])
\(a::TransformFactorization, b::AbstractQuasiVector) = a.plan * convert(Array, b[a.grid])

\(a::TransformFactorization{<:Any,<:Any,Nothing}, b::AbstractVector) = a.iplan \  b
\(a::TransformFactorization, b::AbstractVector) = a.plan * b

function _factorize(::AbstractBasisLayout, L)
    p = grid(L)
    TransformFactorization(p, nothing, factorize(L[p,:]))
end

struct ProjectionFactorization{T, FAC<:Factorization{T}, INDS} <: Factorization{T}
    F::FAC
    inds::INDS
end

\(a::ProjectionFactorization, b::AbstractQuasiVector) = (a.F \ b)[a.inds]
\(a::ProjectionFactorization, b::AbstractVector) = (a.F \ b)[a.inds]

_factorize(::SubBasisLayout, L) = ProjectionFactorization(factorize(parent(L)), parentindices(L)[2])
# function _factorize(::MappedBasisLayout, L)
#     kr, jr = parentindices(L)
#     P = parent(L)
#     ProjectionFactorization(factorize(view(P,:,jr)), parentindices(L)[2])
# end

transform_ldiv(A, B, _) = factorize(A) \ B
transform_ldiv(A, B) = transform_ldiv(A, B, size(A))

copy(L::Ldiv{<:AbstractBasisLayout,<:Any,<:Any,<:AbstractQuasiVector}) =
    transform_ldiv(L.A, L.B)

copy(L::Ldiv{<:AbstractBasisLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) =
    transform_ldiv(L.A, L.B)

##
# Algebra
##

# struct ExpansionLayout <: MemoryLayout end
# applylayout(::Type{typeof(*)}, ::BasisLayout, _) = ExpansionLayout()

const Expansion{T,Space<:AbstractQuasiMatrix,Coeffs<:AbstractVector} = ApplyQuasiVector{T,typeof(*),<:Tuple{Space,Coeffs}}

basis(v::Expansion) = v.args[1]

for op in (:*, :\)
    @eval function broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), x::Number, f::Expansion)
        S,c = arguments(f)
        S * broadcast($op, x, c)
    end
end
for op in (:*, :/)
    @eval function broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), f::Expansion, x::Number)
        S,c = arguments(f)
        S * broadcast($op, c, x)
    end
end

function _broadcastbasis(::typeof(+), _, _, a, b)
    try
        a ≠ b && error("Overload broadcastbasis(::typeof(+), ::$(typeof(a)), ::$(typeof(b)))")
    catch
        error("Overload broadcastbasis(::typeof(+), ::$(typeof(a)), ::$(typeof(b)))")
    end
    a
end

_broadcastbasis(::typeof(+), ::MappedBasisLayouts, ::MappedBasisLayouts, a, b) = broadcastbasis(+, demap(a), demap(b))[basismap(a), :]

broadcastbasis(::typeof(+), a, b) = _broadcastbasis(+, MemoryLayout(a), MemoryLayout(b), a, b)

broadcastbasis(::typeof(-), a, b) = broadcastbasis(+, a, b)

for op in (:+, :-)
    @eval function broadcasted(::LazyQuasiArrayStyle{1}, ::typeof($op), f::Expansion, g::Expansion)
        S,c = arguments(f)
        T,d = arguments(g)
        ST = broadcastbasis($op, S, T)
        ST * $op((ST \ S) * c , (ST \ T) * d)
    end
end

function broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(*), a::Expansion, f::Expansion)
    axes(a,1) == axes(f,1) || throw(DimensionMismatch())
    P,c = arguments(f)
    (a .* P) * c
end

@eval function ==(f::Expansion, g::Expansion)
    S,c = arguments(f)
    T,d = arguments(g)
    ST = broadcastbasis(+, S, T)
    (ST \ S) * c == (ST \ T) * d
end

function show(io::IO, f::Expansion)
    P,c = arguments(f)
    show(io, P)
    print(io, " * ")
    show(io, c)
end
show(io::IO, ::MIME"text/plain", f::Expansion) = show(io, f)

## materialize views

# materialize(S::SubQuasiArray{<:Any,2,<:ApplyQuasiArray{<:Any,2,typeof(*),<:Tuple{<:Basis,<:Any}}}) =
#     *(arguments(S)...)

# mass matrix
# y = p(x), dy = p'(x) * dx
# \int_a^b f(y) g(y) dy = \int_{-1}^1 f(p(x))*g(p(x)) * p'(x) dx

_sub_getindex(A, kr, jr) = A[kr, jr]
_sub_getindex(A, ::Slice, ::Slice) = A

@simplify function *(Ac::QuasiAdjoint{<:Any,<:SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}},
             B::SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
    A = Ac'
    PA,PB = parent(A),parent(B)
    kr,jr = parentindices(B)
    _sub_getindex((PA'PB)/kr.A,parentindices(A)[2],jr)
end

# Differentiation of sub-arrays
@simplify function *(A::Derivative, B::SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:Inclusion,<:Any}})
    P = parent(B)
    (Derivative(axes(P,1))*P)[parentindices(B)...]
end

@simplify function *(A::Derivative, B::SubQuasiArray{<:Any,2,<:AbstractQuasiMatrix,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
    P = parent(B)
    kr,jr = parentindices(B)
    (Derivative(axes(P,1))*P*kr.A)[kr,jr]
end

# we represent as a Mul with a banded matrix
sublayout(::AbstractBasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractUnitRange}}) = SubBasisLayout()
sublayout(::AbstractBasisLayout, ::Type{<:Tuple{<:AbstractAffineQuasiVector,<:AbstractUnitRange}}) = MappedBasisLayout()
sublayout(::WeightedBasisLayouts, ::Type{<:Tuple{<:AbstractAffineQuasiVector,<:AbstractUnitRange}}) = MappedWeightedBasisLayout()
sublayout(::WeightedBasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractUnitRange}}) = SubWeightedBasisLayout()
sublayout(::MappedWeightedBasisLayout, ::Type{<:Tuple{<:Inclusion,<:AbstractUnitRange}}) = MappedWeightedBasisLayout()

@inline sub_materialize(::AbstractBasisLayout, V::AbstractQuasiArray) = V
@inline sub_materialize(::AbstractBasisLayout, V::AbstractArray) = V

demap(x) = x
demap(x::BroadcastQuasiArray) = BroadcastQuasiArray(x.f, map(demap, arguments(x))...)
demap(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}}) = parent(V)
demap(V::SubQuasiArray{<:Any,1,<:Any,<:Tuple{<:AbstractAffineQuasiVector}}) = parent(V)
function demap(V::SubQuasiArray{<:Any,2})
    kr, jr = parentindices(V)
    demap(parent(V)[kr,:])[:,jr]
end

basismap(x::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}) = parentindices(x)[1]
basismap(x::SubQuasiArray{<:Any,1,<:Any,<:Tuple{<:AbstractAffineQuasiVector}}) = parentindices(x)[1]
basismap(x::BroadcastQuasiArray) = basismap(x.args[1])

##
# SubLayout behaves like ApplyLayout{typeof(*)}

combine_mul_styles(::SubBasisLayouts) = combine_mul_styles(ApplyLayout{typeof(*)}())
_mul_arguments(::SubBasisLayouts, A) = _mul_arguments(ApplyLayout{typeof(*)}(), A)
arguments(::SubBasisLayouts, A) = arguments(ApplyLayout{typeof(*)}(), A)
call(::SubBasisLayouts, ::SubQuasiArray) = *

combine_mul_styles(::AdjointSubBasisLayout) = combine_mul_styles(ApplyLayout{typeof(*)}())
_mul_arguments(::AdjointSubBasisLayout, A) = _mul_arguments(ApplyLayout{typeof(*)}(), A)
arguments(::AdjointSubBasisLayout, A) = arguments(ApplyLayout{typeof(*)}(), A)
call(::AdjointSubBasisLayout, ::SubQuasiArray) = *

copy(M::Mul{AdjointSubBasisLayout,<:SubBasisLayouts}) = apply(*, arguments(M.A)..., arguments(M.B)...)

function arguments(::ApplyLayout{typeof(*)}, V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:Inclusion,<:AbstractUnitRange}})
    A = parent(V)
    _,jr = parentindices(V)
    first(jr) ≥ 1 || throw(BoundsError())
    P = _BandedMatrix(Ones{Int}(1,length(jr)), axes(A,2), first(jr)-1,1-first(jr))
    A,P
end

####
# sum
####

_sum(V::AbstractQuasiArray, dims) = __sum(MemoryLayout(typeof(V)), V, dims)
_sum(V::AbstractQuasiArray, ::Colon) = __sum(MemoryLayout(typeof(V)), V, :)
sum(V::AbstractQuasiArray; dims=:) = _sum(V, dims)

__sum(L, Vm, _) = error("Override for $L")
function __sum(::SubBasisLayout, Vm, dims)
    @assert dims == 1
    sum(parent(Vm); dims=dims)[:,parentindices(Vm)[2]]
end
function __sum(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, ::Colon)
    a = arguments(LAY, V)
    first(apply(*, sum(a[1]; dims=1), tail(a)...))
end

function __sum(::MappedBasisLayouts, V::AbstractQuasiArray, dims)
    kr = basismap(V)
    @assert kr isa AbstractAffineQuasiVector
    sum(demap(V); dims=dims)/kr.A
end

include("splines.jl")


text:
function copper_plate(canonical_model::CanonicalModel, expression::Symbol, bus_count::Int64)

    time_steps = model_time_steps(canonical_model)
    devices_netinjection = _remove_undef!(canonical_model.expressions[expression])

    canonical_model.constraints[:CopperPlateBalance] = JuMPConstraintArray(undef, time_steps)

    for t in time_steps
        canonical_model.constraints[:CopperPlateBalance][t] = JuMP.@constraint(canonical_model.JuMPmodel, sum(canonical_model.expressions[expression].data[1:bus_count, t]) == 0)
    end

    return

end


text:
# Free convection in Oceananigans.jl

# Grid setup

using Oceananigans

grid = RegularRectilinearGrid(size = (64, 64, 64),
                              x = (0, 256), y = (0, 256), z = (-128, 0),
                              topology = (Periodic, Periodic, Bounded))

# Boundary conditions and model setup

Q₀ = 1e-8 # Surface buoyancy flux (m² s⁻³)
N² = 1e-5 # Initial and bottom buoyancy gradient (s⁻²)

buoyancy_bcs = TracerBoundaryConditions(grid,
                                        top = BoundaryCondition(Flux, Q₀),
                                        bottom = BoundaryCondition(Gradient, N²))

model = IncompressibleModel(architecture = CPU(),
                            grid = grid,
                            advection = UpwindBiasedFifthOrder(),
                            tracers = :b,
                            buoyancy = BuoyancyTracer(),
                            coriolis = FPlane(f=1e-4),
                            closure = IsotropicDiffusivity(ν=1e-4, κ=1e-4),
                            boundary_conditions = (b=buoyancy_bcs,))

# Set initial condition
noise(z) = 1e-3 * (Q₀ * grid.Δz)^(1/3) * rand()

initial_buoyancy(x, y, z) = N² * z + noise(z)

set!(model, b=initial_buoyancy)

# Simulation setup

using Oceananigans.Units: minute, minutes, hours
using Oceananigans.Utils: prettytime

simulation_time(s) = @info "Time: $(prettytime(s.model.clock.time))"

simulation = Simulation(model,
                        Δt = 1minute,
                        stop_time = 48hours,
                        iteration_interval = 100,
                        progress = simulation_time)

# Output / diagnostics setup

using Oceananigans.AbstractOperations

u, v, w = model.velocities

k = (u^2 + v^2 + w^2) / 2
ξ = ∂y(w) - ∂z(v)

outputs = (kinetic_energy = ComputedField(k),
           x_vorticity = ComputedField(ξ))

simulation.output_writers[:fields] =
    NetCDFOutputWriter(model, outputs,
                       schedule = TimeInterval(30minutes),
                       mode = "c",
                       filepath = "oceananigans_convection.nc")

# Time-step!

run!(simulation)


text:
mystring = "Hello, World!"
println(mystring)


text:
struct DFA
    states
    aps
    transitions
    accepting_state
    sink_state
    initial_state
end

"""
Transition function for any DFA.
"""
function δ(transitions, q, label)
    for relation in transitions
        if relation[1] == q && (relation[2] == label || relation[2] == "true")
            return relation[3]
        end
    end

    # TODO: If no transitions, stay in the current state - this is a workaround to use Spot.jl
    return q
end

"""
Calculate Transition reward
"""
function distance_from_accept_state(dfa, dfa_state)

    @assert dfa_state ∈ dfa.states

    if dfa_state == dfa.sink_state
        return Inf
    end

    if dfa_state == dfa.accepting_state
        return 0
    end

    states_to_try = [dfa.accepting_state]
    dist = 1

    while !isempty(states_to_try)
        for state in states_to_try
            for relation in dfa.transitions
                if relation[3] == state
                    if relation[1] == dfa_state
                        # @info relation[1], state
                        return dist
                    else
                        push!(states_to_try, relation[1])
                        # @info states_to_try
                    end
                end
            end
            dist += 1
            setdiff!(states_to_try, state)
        end
    end

    return dist
end

function create_dot_graph(dfa::DFA, filename::String)
    open(filename, "w") do f
        println(f, "digraph G {")
        println(f, "  rankdir=LR\n  node [shape=\"circle\"]\n  fontname=\"Lato\"\n  node [fontname=\"Lato\"]\n  edge [fontname=\"Lato\"]")
        println(f, "  size=\"8.2,8.2\" node[style=filled,fillcolor=\"#FDEDD3\"] edge[arrowhead=vee, arrowsize=.7]")

        # Initial State
        @printf(f, "  I [label=\"\", style=invis, width=0]\n  I -> %d\n", dfa.initial_state)

        for state in dfa.states
            if state == dfa.accepting_state
                @printf(f, "  %d [shape=\"doublecircle\", label=<z%d>]\n", state, state)
                @printf(f, "  %d -> %d [label=<true>]\n", state, state)
            else
                @printf(f, "  %d [label=<z%d>]\n", state, state)
                for transition in dfa.transitions
                    if transition[1] == state
                        # Assume 1 label for now
                        @printf(f, "  %d -> %d [label=<%s>]\n", state, transition[3], transition[2])
                    end
                end

            end

        end
        println(f, "}")
    end
end

"""
Given the current dfa state, return labels that will induce a positive transition
"""
function get_next_dfa_labels(dfa, dfa_state)
    labels = []
    for relation in dfa.transitions
        if relation[1] == dfa_state &&  relation[1] != relation[3] && relation[3] != dfa.sink_state
            push!(labels, relation[2])
        end
    end
    return labels
end

"""
Get all labels that transition to sink state
"""
function get_dfa_sink_labels(dfa)
    labels = []
    for relation in dfa.transitions
        if relation[1] != relation[3] && relation[3] == dfa.sink_state
            push!(labels, relation[2])
        end
    end
    return labels
end

"""
Check if current state can transition to accept state
"""
function possible_accept_transition(dfa, dfa_state, symbol)
    for relation in dfa.transitions
        if dfa_state == relation[1] && symbol == relation[2] && relation[3] == dfa.accepting_state
            return true
        end
    end
    return false
end

"""
Construct DFA object using Spot
"""
function construct_DFA_from_LTL(ltl_form)
    # ltl_form = ltl"$ltl_string"
    translator = LTLTranslator(deterministic=true)
    spot_automaton = Spot.translate(translator, ltl_form)

    transitions = []
    for (edge, label) in zip(get_edges(spot_automaton), get_labels(spot_automaton))
        label = string(label) == "1" ? "true" :  create_full_label(atomic_propositions(spot_automaton), label)
        new_transition = (edge[1], label, edge[2])
        push!(transitions, new_transition)
    end

    # TODO: add transition to the/a violating state
    # ! It doesn't exist here currently
    accept_states = get_rabin_acceptance(spot_automaton)
    accept_state = collect(accept_states[1][2])[1]

    dfa = DFA(collect(1:num_states(spot_automaton)), string.(atomic_propositions(spot_automaton)), transitions, accept_state, nothing, get_init_state_number(spot_automaton))

    return dfa
end

function create_full_label(props, label)
    # all props
    label_props = atomic_prop_collect(label)
    if props == label_props
        return string(label)
    end

    # rebuild
    label = string(label)
    new_label = ""
    for prop in string.(props)
        if occursin("!$prop", label)
            new_label = "$new_label & !$prop"
        elseif occursin("$prop", label)
            new_label = "$new_label & $prop"
        else
            new_label = "$new_label & !$prop"
        end
    end
    new_label = new_label[4:end]
    return new_label
end


text:
module Libcint

using libcint_jll

#####
##### Functions names
#####

const dim = [
    "CINTgto_cart",
    "CINTgto_spheric",
    "CINTgto_spinor",
    ]

const f1e = [
    "cint1e_ovlp_cart",
    "cint1e_nuc_cart",
    "cint1e_kin_cart",
    "cint1e_ia01p_cart",
    "cint1e_irixp_cart",
    "cint1e_ircxp_cart",
    "cint1e_iking_cart",
    "cint1e_iovlpg_cart",
    "cint1e_inucg_cart",
    "cint1e_ipovlp_cart",
    "cint1e_ipkin_cart",
    "cint1e_ipnuc_cart",
    "cint1e_iprinv_cart",
    "cint1e_rinv_cart",
    "cint2e_cart",
    "cint2e_ig1_cart",
    "cint2e_ip1_cart",
    "cint1e_ovlp_sph",
    "cint1e_nuc_sph",
    "cint1e_kin_sph",
    "cint1e_ia01p_sph",
    "cint1e_irixp_sph",
    "cint1e_ircxp_sph",
    "cint1e_iking_sph",
    "cint1e_iovlpg_sph",
    "cint1e_inucg_sph",
    "cint1e_ipovlp_sph",
    "cint1e_ipkin_sph",
    "cint1e_ipnuc_sph",
    "cint1e_iprinv_sph",
    "cint1e_rinv_sph",
    "cint2e_sph",
    "cint2e_ig1_sph",
    "cint2e_ip1_sph",
    "cint1e_ovlp",
    "cint1e_nuc",
    "cint1e_nucg",
    "cint1e_srsr",
    "cint1e_sr",
    "cint1e_srsp",
    "cint1e_spsp",
    "cint1e_sp",
    "cint1e_spspsp",
    "cint1e_spnuc",
    "cint1e_spnucsp",
    "cint1e_srnucsr",
    "cint1e_sa10sa01",
    "cint1e_ovlpg",
    "cint1e_sa10sp",
    "cint1e_sa10nucsp",
    "cint1e_sa01sp",
    "cint1e_spgsp",
    "cint1e_spgnucsp",
    "cint1e_spgsa01",
    "cint1e_ipovlp",
    "cint1e_ipkin",
    "cint1e_ipnuc",
    "cint1e_iprinv",
    "cint1e_ipspnucsp",
    "cint1e_ipsprinvsp",
    "cint2e",
    "cint2e_spsp1",
    "cint2e_spsp1spsp2",
    "cint2e_srsr1",
    "cint2e_srsr1srsr2",
    "cint2e_sa10sp1",
    "cint2e_sa10sp1spsp2",
    "cint2e_g1",
    "cint2e_spgsp1",
    "cint2e_g1spsp2",
    "cint2e_spgsp1spsp2",
    "cint2e_ip1",
    "cint2e_ipspsp1",
    "cint2e_ip1spsp2",
    "cint2e_ipspsp1spsp2",
    "cint2e_ssp1ssp2",
    ]

#####
##### Function wrapper macros
#####

macro make_dim(name::String)
    fn = Symbol(name)
    _fn = Symbol('_'*name)
    quote
        function $(esc(fn))(
            i::Cint,
            bas::AbstractVector{Cint}
            )
            n = ccall(
                ($_fn, libcint),
                Cint,
                (Cint, Ref{Cint}),
                i, bas,
                )
        end
    end
end

macro make_fn(name::String)
    fn = Symbol(name)
    _fn = Symbol('_'*name)
    quote
        function $(esc(fn))(
            buf::AbstractArray{Cdouble},
            shls::AbstractVector{Cint},
            atm::AbstractVector{Cint},
            natm::Cint,
            bas::AbstractVector{Cint},
            nbas::Cint,
            env::AbstractVector{Cdouble},
            )
            not0 = ccall(
                ($_fn, libcint),
                Cint,
                (Ref{Cdouble}, Ref{Cint}, Ref{Cint}, Cint, Ref{Cint}, Cint, Ref{Cdouble}, Ptr{Cvoid}),
                buf, shls, atm, natm, bas, nbas, env, C_NULL,
                )
        end
    end
end

#####
##### Define functions
#####

for s in dim
    _s = Meta.parse('_'*s);
    @eval const $_s = $s
    @eval @make_dim($s)
end

for s in f1e
    _s = Meta.parse('_'*s);
    @eval const $_s = $s
    @eval @make_fn($s)
end

end # module


text:
using LoggingFacilities
using Logging

logger = compose(
    current_logger(),
    logger -> FixedMessageWidthTransformerLogger(logger, 40),
    OneLineTransformerLogger,
)

with_logger(logger) do
    for i in 1:10
        x = round(rand(), digits = 10)
        y = round(rand(), digits = 10)
        @info "Iteration #$i" x y
    end
end

#=
[ Info: Iteration #1                             x=0.5674366266778048 y=0.9024632811266229
[ Info: Iteration #2                             x=0.7070619211896252 y=0.44200367657577777
[ Info: Iteration #3                             x=0.034904964067742794 y=0.09434400347104588
[ Info: Iteration #4                             x=0.38486209773503965 y=0.023043539047058914
[ Info: Iteration #5                             x=0.8135371107442477 y=0.9173422607084372
[ Info: Iteration #6                             x=0.9445325770356989 y=0.3864405581804362
[ Info: Iteration #7                             x=0.850683727155505 y=0.6299730919535942
[ Info: Iteration #8                             x=0.09683784334946521 y=0.9662550162409316
[ Info: Iteration #9                             x=0.7733711241925625 y=0.7778339378610857
[ Info: Iteration #10                            x=0.13435678232322879 y=0.546560581934818
=#


text:
# Contains common functions for all grid types.

"""
Types allowed as input to `vtk_point_data()` and `vtk_cell_data()`.

Either (abstract) arrays or tuples of arrays are allowed.
In the second case, the length of the tuple determines the number of components
of the input data (e.g. if N = 3, it corresponds to a 3D vector field).
"""
const InputDataType =
    Union{AbstractArray,
          NTuple{N, T} where {N, T <: AbstractArray}}

"""
Determine number of components of input data.
"""
function num_components(data::AbstractArray, num_points_or_cells::Int)
    Nc = div(length(data), num_points_or_cells)
    if Nc * num_points_or_cells != length(data)
        throw(ArgumentError("Incorrect dimensions of input array."))
    end
    Nc
end

num_components(data::NTuple, ::Int) = length(data)

"Union of data types allowed by VTK (see file-formats.pdf, page 15)."
const VTKDataType = Union{Int8, UInt8, Int16, UInt16, Int32, UInt32,
                          Int64, UInt64, Float32, Float64}

"""
Return the VTK string representation of a numerical data type.
"""
function datatype_str(::Type{T}) where T <: VTKDataType
    # Note: the VTK type names are exactly the same as the Julia type names
    # (e.g.  Float64 -> "Float64"), so that we can simply use the `string`
    # function.
    string(T)
end
datatype_str(::Type{T}) where T =
    throw(ArgumentError("Data type not supported by VTK: $T"))
datatype_str(::AbstractArray{T}) where T = datatype_str(T)
datatype_str(::NTuple{N, T} where N) where T <: AbstractArray =
    datatype_str(eltype(T))

"""
Total size of data in bytes.
"""
sizeof_data(x::Array) = sizeof(x)
sizeof_data(x::AbstractArray) = length(x) * sizeof(eltype(x))
sizeof_data(x::NTuple) = sum(sizeof_data, x)

"Write array of numerical data to stream."
function write_array(io, data::AbstractArray)
    write(io, data)
    nothing
end

function write_array(io, data::NTuple{N, T} where T <: AbstractArray) where N
    # We assume that all arrays in the tuple have the same size.
    L = length(data[1])
    for i = 1:L, x in data
        write(io, x[i])
    end
    nothing
end

"""
Add numerical data to VTK XML file.

Data is written under the `xParent` XML node.

`Nc` corresponds to the number of components of the data.
"""
function data_to_xml(vtk::DatasetFile, xParent::XMLElement, data::InputDataType,
                     varname::AbstractString, Nc::Integer=1)
    @assert name(xParent) in ("Points", "PointData", "Coordinates", "Cells",
                              "CellData")
    func :: Function = vtk.appended ? data_to_xml_appended : data_to_xml_inline
    func(vtk, xParent, data, varname, Nc) :: XMLElement
end

"""
Add appended raw binary data to VTK XML file.

Data is written to the `vtk.buf` buffer.

When compression is enabled:

  * the data array is written in compressed form (obviously);

  * the header, written before the actual numerical data, is an array of UInt32
    values:
        `[num_blocks, blocksize, last_blocksize, compressed_blocksizes]`
    All the sizes are in bytes. The header itself is not compressed, only the
    data is.
    For more details, see:
        http://public.kitware.com/pipermail/paraview/2005-April/001391.html
        http://mathema.tician.de/what-they-dont-tell-you-about-vtk-xml-binary-formats
    (This is not really documented in the VTK specification...)

Otherwise, if compression is disabled, the header is just a single UInt32 value
containing the size of the data array in bytes.

"""
function data_to_xml_appended(vtk::DatasetFile, xParent::XMLElement,
                              data::InputDataType, varname::AbstractString,
                              Nc::Integer)
    @assert vtk.appended

    buf = vtk.buf    # append buffer
    compress = vtk.compression_level > 0

    # DataArray node
    xDA = new_child(xParent, "DataArray")
    set_attribute(xDA, "type", datatype_str(data))
    set_attribute(xDA, "Name", varname)
    set_attribute(xDA, "format", "appended")
    set_attribute(xDA, "offset", position(buf))
    set_attribute(xDA, "NumberOfComponents", string(Nc))

    # Size of data array (in bytes).
    nb = sizeof_data(data)

    if compress
        initpos = position(buf)

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

        # Write compressed data.
        zWriter = ZlibCompressorStream(buf, level=vtk.compression_level)
        write_array(zWriter, data)
        write(zWriter, TranscodingStreams.TOKEN_END)
        flush(zWriter)

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

    xDA::XMLElement
end

"Add inline, base64-encoded data to VTK XML file."
function data_to_xml_inline(vtk::DatasetFile, xParent::XMLElement,
                            data::InputDataType, varname::AbstractString,
                            Nc::Integer)
    @assert !vtk.appended
    compress = vtk.compression_level > 0

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

    # Number of bytes of data.
    nb = sizeof_data(data)

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

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

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

    if compress
        close(zWriter)
    end
    close(buf)

    xDA::XMLElement
end

"""
Add either point or cell data to VTK file.

Here `Nc` is the number of components of the data (Nc >= 1).
"""
function vtk_point_or_cell_data(vtk::DatasetFile, data::InputDataType,
                                name::AbstractString, nodetype::AbstractString,
                                Nc::Integer)
    # Find Piece node.
    xroot = root(vtk.xdoc)
    xGrid = find_element(xroot, vtk.grid_type)
    xPiece = find_element(xGrid, "Piece")

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

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

    vtk
end

function vtk_point_data(vtk::DatasetFile, data::InputDataType, name::AbstractString)
    Nc = num_components(data, vtk.Npts)
    vtk_point_or_cell_data(vtk, data, name, "PointData", Nc)
end

function vtk_cell_data(vtk::DatasetFile, data::InputDataType, name::AbstractString)
    Nc = num_components(data, vtk.Ncls)
    vtk_point_or_cell_data(vtk, data, name, "CellData", Nc)
end

function vtk_save(vtk::DatasetFile)
    if isopen(vtk)
        if vtk.appended
            save_with_appended_data(vtk)
        else
            save_file(vtk.xdoc, vtk.path)
        end
        @assert isopen(vtk)
        close(vtk)
    end
    return [vtk.path] :: Vector{String}
end

"""
Write VTK XML file containing appended binary data to disk.

In this case, the XML file is written manually instead of using the `save_file`
function of `LightXML`, which doesn't allow to write raw binary data.
"""
function save_with_appended_data(vtk::DatasetFile)
    @assert vtk.appended
    @assert isopen(vtk.buf)

    # Convert XML document to a string, and split the last two lines.
    lines = rsplit(string(vtk.xdoc), '\n', limit=3, keepempty=true)

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

    open(vtk.path, "w") do io
        # Write everything but the last two lines.
        write(io, lines[1])
        write(io, "\n")

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

        close(vtk.buf)
    end

    nothing
end

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

"""
Return the "extent" attribute required for structured (including rectilinear)
grids.
"""
extent_attribute(Ni, Nj, Nk, ::Nothing=nothing) =
    "0 $(Ni - 1) 0 $(Nj - 1) 0 $(Nk - 1)"

function extent_attribute(Ni, Nj, Nk, extent::Array{T}) where T <: Integer
    length(extent) == 6 || throw(ArgumentError("Extent must have length 6."))
    (extent[2] - extent[1] + 1 == Ni) &&
    (extent[4] - extent[3] + 1 == Nj) &&
    (extent[6] - extent[5] + 1 == Nk) ||
    throw(ArgumentError("Extent is not consistent with dataset dimensions."))
    join(extent, " ")
end

"""
Number of cells in structured grids.

In 3D, all cells are hexahedrons (i.e. VTK_HEXAHEDRON), and the number of
cells is (Ni-1)*(Nj-1)*(Nk-1). In 2D, they are quadrilaterals (VTK_QUAD), and in
1D they are line segments (VTK_LINE).
"""
function num_cells_structured(Ni, Nj, Nk)
    Ncls = one(Ni)
    for N in (Ni, Nj, Nk)
        Ncls *= max(1, N - 1)
    end
    Ncls
end

# Accept passing coordinates as a tuple.
vtk_grid(filename, xyz::NTuple{3, T} where T; kwargs...) =
    vtk_grid(filename, xyz...; kwargs...)
vtk_grid(filename, xy::NTuple{2, T} where T; kwargs...) =
    vtk_grid(filename, xy...; kwargs...)


text:
using StanVariational
using Test

bernoulli_model = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

bernoulli_data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

stanmodel = VariationalModel("bernoulli", bernoulli_model)

(sample_file, log_file) = stan_variational(stanmodel; data=bernoulli_data)

if sample_file !== Nothing

  @testset "Bernoulli variational example" begin
    # Read sample summary (in ChainDataFrame format)
    sdf = read_summary(stanmodel)

    @test sdf[:theta, :mean][1] ≈ 0.32 atol=0.1
  end

end


text:
using DashBootstrapComponents, DashHtmlComponents

pagination = html_div(dbc_pagination(max_value = 10))

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
const IGC_TIME_FMT = "HHMMSS"
const IGC_DATE_FMT = "ddmmyy"

struct IGCTimeFormat
end

struct IGCDateFormat
end