text:
using LinearAlgebra, Plots, JLD2
using Flux, Flux.Zygote, Optim, DiffEqFlux
using KitBase, ProgressMeter
import KitML

maxwell_boltzmann(x) = x * log(x) - x

maxwell_boltzmann_prime(x) = log(x)

function kinetic_entropy(α, m, weights)
    B = KitBase.maxwell_boltzmann_dual_prime.(α' * m)[:]
    return sum(maxwell_boltzmann.(B) .* weights)
end

function is_absorb(x::T, y::T) where {T<:Real}
    cds = Array{Bool}(undef, 11) # conditions

    cds[1] = -2.5<x<-1.5 && 1.5<y<2.5
    cds[2] = -2.5<x<-1.5 && -0.5<y<0.5
    cds[3] = -2.5<x<-1.5 && -2.5<y<-1.5
    cds[4] = -1.5<x<-0.5 && 0.5<y<1.5
    cds[5] = -1.5<x<-0.5 && -1.5<y<-0.5
    cds[6] = -0.5<x<0.5 && -2.5<y<-1.5
    cds[7] = 0.5<x<1.5 && 0.5<y<1.5
    cds[8] = 0.5<x<1.5 && -1.5<y<-0.5
    cds[9] = 1.5<x<2.5 && 1.5<y<2.5
    cds[10] = 1.5<x<2.5 && -0.5<y<0.5
    cds[11] = 1.5<x<2.5 && -2.5<y<-1.5

    if any(cds) == true
        return true
    else
        return false
    end
end

cd(@__DIR__)
include("math.jl")
model = KitML.load_model("best_model.h5"; mode = :tf)

begin
    # space
    x0 = -3.5
    x1 = 3.5
    y0 = -3.5
    y1 = 3.5
    nx = 80#100
    ny = 80#100
    dx = (x1 - x0) / nx
    dy = (y1 - y0) / ny

    pspace = KitBase.PSpace2D(x0, x1, nx, y0, y1, ny)

    # time
    tEnd = 1.0
    cfl = 0.5
    dt = cfl / 2 * (dx * dy) / (dx + dy)

    # quadrature
    quadratureorder = 5
    points, triangulation = KitBase.octa_quadrature(quadratureorder)
    weights = KitBase.quadrature_weights(points, triangulation)
    nq = size(points, 1)

    # moments
    L = 1
    #ne = (L + 1)^2
    ne = GetBasisSize(L,3)
    phi = zeros(Float32, ne, nx, ny)
    phi[1, :, :] .= 1e-6
    α = zeros(Float32, ne, nx, ny)
    #m = KitBase.eval_spherharmonic(points, L)
    m = ComputeSphericalBasisAnalytical( points)
end

begin
    σs = zeros(Float64, nx, ny)
    σa = zeros(Float64, nx, ny)
    for i = 1:nx, j = 1:ny
        if is_absorb(pspace.x[i, j], pspace.y[i, j])
            σs[i, j] = 0.0
            σa[i, j] = 10.0
        else
            σs[i, j] = 1.0
            σa[i, j] = 0.0
        end
    end
    σt = σs + σa
    σq = zeros(Float64, nx, ny)
    for i = 1:nx, j = 1:ny
        if -0.5<pspace.x[i, j]<0.5 && -0.5<pspace.y[i, j]<0.5
            σq[i, j] = 1.0 / (4.0 * π)
        else
            σq[i, j] = 0.0
        end
    end
end

#contourf(pspace.x[1:nx, 1], pspace.y[1, 1:ny], σs')

global t = 0.0
flux1 = zeros(Float32, ne, nx + 1, ny)
flux2 = zeros(Float32, ne, nx, ny + 1)

αT = zeros(Float32, nx*ny, ne)
phiT = zeros(Float32, nx*ny, ne)

global X = zeros(Float32, 1, ne)
global Y = zeros(Float32, 1, 1)

res = KitBase.optimize_closure(X[1, :], m, weights, phi[:, nx÷2, ny÷2], KitBase.maxwell_boltzmann_dual)
X[1, :] .= res.minimizer
Y[1, 1] = kinetic_entropy(X[1, :], m, weights)

#=
# compare mathematical & NN optimizer
α0 = similar(α)
@inbounds Threads.@threads for j = 1:ny
    for i = 1:nx
        res = KitBase.optimize_closure(α[:, i, j], m, weights, phi[:, i, j], KitBase.maxwell_boltzmann_dual)
        α0[:, i, j] .= res.minimizer
    end
end

@inbounds for i = 1:ne
    phiT[:, i] .= phi[i, :, :][:]
end
αT = KitML.neural_closure(model, phiT)
α1 = similar(α)
@inbounds for i = 1:ne
    α1[i, :, :] .= reshape(αT[:, i], nx, ny)
end

plot(pspace.x[1:nx, 1], α0[1, :, 40], label="Newton", xlabel="x", ylabel="alpha_0")
plot!(pspace.x[1:nx, 1], α1[1, :, 40], label="NN")
=#

phi_temp = zero(phi)
phi_old = zero(phi)
anim = @animate for iter = 1:1000#20
    println("Iteration $(iter)")

    phi_old .= phi

    # regularization
    @inbounds for i = 1:ne
        phiT[:, i] .= phi[i, :, :][:]
    end
    αT .= KitML.neural_closure(model, phiT)
    @inbounds for i = 1:ne
        α[i, :, :] .= reshape(αT[:, i], nx, ny)
    end
    @inbounds Threads.@threads for j = 1:ny
        for i = 1:nx
            phi_temp[:, i, j] .= KitBase.realizable_reconstruct(α[:, i, j], m, weights, KitBase.maxwell_boltzmann_dual_prime)
        end
    end

    @inbounds for j = 1:ny
        for i = 1:nx
            @show norm(phi_temp[:, i, j] - phi_old[:, i, j], 2)

            if norm(phi_temp[:, i, j] - phi_old[:, i, j], 2) > 1e-1
                res = KitBase.optimize_closure(α[:, i, j], m, weights, phi[:, i, j], KitBase.maxwell_boltzmann_dual)
                α[:, i, j] .= res.minimizer
                phi[:, i, j] .= KitBase.realizable_reconstruct(res.minimizer, m, weights, KitBase.maxwell_boltzmann_dual_prime)

                if phi[1, i, j] > 0.01
                    X = vcat(X, permutedims(phi[:, i, j]))
                    Y = vcat(Y, kinetic_entropy(α[:, i, j], m, weights))
                end
            else
                phi[:, i, j] .= phi_temp[:, i, j]
            end
        end
    end

    # flux
    fη1 = zeros(nq)
    @inbounds for j = 1:ny
        for i = 2:nx
            KitBase.flux_kfvs!(fη1, KitBase.maxwell_boltzmann_dual.(α[:, i-1, j]' * m)[:], KitBase.maxwell_boltzmann_dual.(α[:, i, j]' * m)[:], points[:, 1], dt)

            for k in axes(flux1, 1)
                flux1[k, i, j] = sum(m[k, :] .* weights .* fη1)
            end
        end
    end

    fη2 = zeros(nq)
    @inbounds for i = 1:nx
        for j = 2:ny
            KitBase.flux_kfvs!(fη2, KitBase.maxwell_boltzmann_dual.(α[:, i, j-1]' * m)[:], KitBase.maxwell_boltzmann_dual.(α[:, i, j]' * m)[:], points[:, 2], dt)

            for k in axes(flux2, 1)
                flux2[k, i, j] = sum(m[k, :] .* (weights .* fη2))
            end
        end
    end

    # update
    @inbounds for j = 2:ny-1
        for i = 2:nx-1
            for q = 1:1
                phi[q, i, j] =
                    phi[q, i, j] +
                    (flux1[q, i, j] - flux1[q, i+1, j]) / dx +
                    (flux2[q, i, j] - flux2[q, i, j+1]) / dy +
                    (σs[i, j] * phi[q, i, j] - σt[i, j] * phi[q, i, j]) * dt +
                    σq[i, j] * dt * 100.0
            end

            for q = 2:ne
                phi[q, i, j] =
                    phi[q, i, j] +
                    (flux1[q, i, j] - flux1[q, i+1, j]) / dx +
                    (flux2[q, i, j] - flux2[q, i, j+1]) / dy +
                    (-σt[i, j] * phi[q, i, j]) * dt
            end
        end
    end

    global t += dt
    contourf(pspace.x[1:nx, 1], pspace.y[1, 1:ny], phi[1, :, :])

    if iter%9 == 0
        model.fit(X, Y, epochs=2)

        X = zeros(Float32, 1, ne)
        Y = zeros(Float32, 1, 1)
        res = KitBase.optimize_closure(X[1, :], m, weights, phi[:, nx÷2, ny÷2], KitBase.maxwell_boltzmann_dual)
        X[1, :] .= res.minimizer
        Y[1, 1] = kinetic_entropy(X[1, :], m, weights)
    end
end

cd(@__DIR__)
gif(anim, "lattice_neural.gif")

#contourf(pspace.x[1:nx, 1], pspace.y[1, 1:ny], phi[1, :, :])


text:
using Test, BlockArrays

@testset "blocks" begin
    @testset "blocks(::BlockVector)" begin
        vector_blocks = [[1, 2], [3, 4, 5], Int[]]
        @test blocks(mortar(vector_blocks)) === vector_blocks
        @test collect(blocks(mortar(vector_blocks))) == vector_blocks
    end

    @testset "blocks(::BlockMatrix)" begin
        matrix_blocks = permutedims(reshape([
            1ones(1, 3), 2ones(1, 2),
            3ones(2, 3), 4ones(2, 2),
        ], (2, 2)))
        @test blocks(mortar(matrix_blocks)) === matrix_blocks
    end

    @testset "blocks(::PseudoBlockVector)" begin
        v0 = rand(3)
        vb = PseudoBlockArray(v0, [1, 2])
        @test size(blocks(vb)) == (2,)
        blocks(vb)[1] = [123]
        @test v0[1] == 123
        @test parent(blocks(vb)[1]) === v0

        # toplevel = true:
        str = sprint(show, "text/plain", blocks(vb))
        @test occursin("blocks of PseudoBlockArray of", str)

        # toplevel = false:
        str = sprint(show, "text/plain", view(blocks(vb), 1:1))
        @test occursin("::BlocksView{…,::PseudoBlockArray{…,", str)
    end

    @testset "blocks(::PseudoBlockMatrix)" begin
        m0 = rand(2, 4)
        mb = PseudoBlockArray(m0, [1, 1], [2, 1, 1])
        @test size(blocks(mb)) == (2, 3)
        blocks(mb)[1, 1] = [123 456]
        @test m0[1, 1] == 123
        @test m0[1, 2] == 456
        @test parent(blocks(mb)[1, 1]) === m0

        # linear indexing
        @test blocks(mb)[1] == m0[1:1, 1:2]
        blocks(mb)[1] = [111 222]
        @test mb[Block(1, 1)] == [111 222]

        # toplevel = true:
        str = sprint(show, "text/plain", blocks(mb))
        @test occursin("blocks of PseudoBlockArray of", str)

        # toplevel = false:
        str = sprint(show, "text/plain", view(blocks(mb), 1:1, 1:1))
        @test occursin("::BlocksView{…,::PseudoBlockArray{…,", str)
    end

    @testset "blocks(::Vector)" begin
        v = rand(3)
        @test size(blocks(v)) == (1,)
        blocks(v)[1][1] = 123
        @test v[1] == 123
        @test parent(blocks(v)[1]) === v
    end

    @testset "blocks(::Matrix)" begin
        m = rand(2, 4)
        @test size(blocks(m)) == (1, 1)
        blocks(m)[1, 1][1, 1] = 123
        @test m[1, 1] == 123
        @test parent(blocks(m)[1, 1]) === m
    end

    @testset "blocks(::Adjoint|Transpose)" begin
        m = BlockArray([rand(ComplexF64, 2, 2) for _ in 1:3, _ in 1:5], [1, 2], [2, 3])
        @testset for i in 1:2, j in 1:2
            @test blocks(m')[i, j] == m'[Block(i), Block(j)]
            @test blocks(transpose(m))[i, j] == transpose(m)[Block(i), Block(j)]
        end
    end

    @testset "blocks(::SubArray)" begin
        vector_blocks = [[1, 2], [3, 4, 5], Int[]]
        b = view(mortar(vector_blocks), Block(1):Block(2))
        v = blocks(b)
        @test v == vector_blocks[1:2]
        v[1][1] = 111
        @test b[1] == 111
        @test parent(v) === parent(b).blocks  # special path works
    end

    @testset "blocks(::SubArray)" begin
        matrix_blocks = permutedims(reshape([
            1ones(1, 3), 2ones(1, 2),
            3ones(2, 3), 4ones(2, 2),
        ], (2, 2)))
        b = view(mortar(matrix_blocks), Block(1):Block(2), Block(2):Block(2))
        m = blocks(b)
        @test m == matrix_blocks[1:2, 2:2]
        m[1, 1][1, 1] = 111
        @test b[1, 1] == 111
        @test parent(m) === parent(b).blocks  # special path works
    end
end

ext: jl


text:
# ForwardDiff integration

for f in libdevice
  if haskey(ForwardDiff.DiffRules.DEFINED_DIFFRULES, (:Base,f,1))
    f == :tanh && continue
    diffrule = ForwardDiff.DiffRules.DEFINED_DIFFRULES[(:Base,f,1)]
    ForwardDiff.DiffRules.DEFINED_DIFFRULES[(:CUDAnative,f,1)] =
      (args...) -> replace_device(diffrule(args...))
    eval(ForwardDiff.unary_dual_definition(:CUDAnative, f))
  end
end

ForwardDiff.DiffRules.DEFINED_DIFFRULES[(:CUDAnative, :tanh, 1)] = x ->
  replace_device(:(1-tanh(x)^2))
eval(ForwardDiff.unary_dual_definition(:CUDAnative, :tanh))


text:
# __precompile__(false)

module Routing

include("Routing_include.jl")

export  Label, ESPPRC_Instance, VRPTW_Instance,
        Node, Fleet, Request, Solomon,
        solveESPPRC, solveESPPRC_vrp, solveVRP,
        read_solomon_data, calculate_solomon_cost,
        calculate_path_time, calculate_path_cost,
        load_solomon, generate_solomon_vrptw_instance,
        show_details, solomon_to_espprc, print_pct_neg_arcs,
        MAX_INT

end # module


text:
using .HyperDualNumbers, FFTW
import FFTW.fft!;

function fft!(x::AbstractArray{T}) where {T <: Hyper{cpx}} where {cpx <: Complex{num}} where {num <: Number}
    Xcomplex, Xder1, Xder2, Xder12 = HyperDualNumbers.value.(x), HyperDualNumbers.epsilon1.(x), HyperDualNumbers.epsilon2.(x), HyperDualNumbers.epsilon12.(x)
    planned_fft = plan_fft!(similar(Xcomplex))
    Xcomplex .= planned_fft * Xcomplex
    Xder1 .= planned_fft * Xder1
    Xder2 .= planned_fft * Xder2
    Xder12 .= planned_fft * Xder12
    @. x = hyper(Xcomplex, Xder1, Xder2, Xder12)
    nothing
end

function real_mod(x::Hyper)
    return hyper(real(x.value), real(x.epsilon1), real(x.epsilon2), real(x.epsilon12))
end


text:
"""
    binary_tournament(population)
Apply binary tournament to obtain a solution from from population.
"""
function binary_tournament(population)
    a, b = rand(1:length(population), 2)

    return is_better(population[a], population[b]) ? population[a] : population[b]
end

function binary_tournament(population, fitness::Vector{T})  where T <: Real
    a, b = rand(1:length(population), 2)

    return fitness[a] < fitness[b] ? population[a] : population[b]
end


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

function set_terminal(term::String)
    global gaston_config

    if validate_terminal(term)
        gaston_config.terminal = term
    else
        error(string("Terminal type ", term, " not supported."))
    end
    return term
end

# For terminals that support file output, set the file name
function set_filename(s::String)
    global gaston_config

    gaston_config.outputfile = s
end

# Set default values for CurveConf and AxesConf. Return set value upon
# success, error out otherwise.
function set_default_legend(s::String)
    global gaston_config

    gaston_config.legend = s
end

function set_default_plotstyle(s::String)
    global gaston_config

    if validate_2d_plotstyle(s) || validate_3d_plotstyle(s) ||
        validate_image_plotstyle(s)
        gaston_config.plotstyle = s
    else
        error(string("Plotstyle ", s, " not supported."))
    end
end

function set_default_color(s::String)
    global gaston_config

    # We can't validate colors yet
    gaston_config.color = s
end

function set_default_marker(s::String)
    global gaston_config

    if validate_marker(s)
        gaston_config.marker = s
    else
        error(string("Marker ", s, " not supported."))
    end
end

function set_default_linewidth(i::Real)
    global gaston_config

    gaston_config.linewidth = i
end

function set_default_pointsize(i::Real)
    global gaston_config

    gaston_config.pointsize = i
end

function set_default_title(s::String)
    global gaston_config

    gaston_config.title = s
end

function set_default_xlabel(s::String)
    global gaston_config

    gaston_config.xlabel = s
end

function set_default_ylabel(s::String)
    global gaston_config

    gaston_config.ylabel = s
end

function set_default_zlabel(s::String)
    global gaston_config

    gaston_config.zlabel = s
end

function set_default_box(s::String)
    global gaston_config

    # We can't validate box yet
    gaston_config.box = s
end

function set_default_axis(s::String)
    global gaston_config

    if validate_axis(s)
        gaston_config.marker = s
    else
        error(string("Axis ", s, " not supported."))
    end
end

# functions to set print parameters
function set_print_color(s::String)
    global gaston_config

    gaston_config.print_color = b
end

function set_print_fontface(s::String)
    global gaston_config

    gaston_config.print_fontface = s
end

function set_print_fontsize(s::Real)
    global gaston_config

    gaston_config.print_fontsize = s
end

function set_print_fontscale(r::Real)
    global gaston_config

    gaston_config.print_fontscale = r
end

function set_print_linewidth(r::Real)
    global gaston_config

    gaston_config.print_linewidth = r
end

function set_print_size(s::String)
    global gaston_config

    gaston_config.print_size = s
end

function set_default_grid(s::String)
    global gaston_config

    if s=="on" || s=="off"
        gaston_config.grid=s
    else
        error(string("grid ", s, " not supported."))
    end
end


text:
# define common PETSc constants
# this excludes configurable constants (e.g. PetscScalar) which are set in lib.jl

const PetscErrorCode = Cint

struct PetscError <: Exception
    code::PetscErrorCode
end

macro chk(expr)
    :((errcode = $(esc(expr))) == 0 || throw(PetscError(errcode)))
end

const PETSC_DEFAULT = -2
const PETSC_DECIDE = PETSC_DETERMINE = -1

@enum PetscBool PETSC_FALSE PETSC_TRUE

@enum PetscDataType begin
    PETSC_DATATYPE_UNKNOWN = 0
    PETSC_DOUBLE = 1
    PETSC_COMPLEX = 2
    PETSC_LONG = 3
    PETSC_SHORT = 4
    PETSC_FLOAT = 5
    PETSC_CHAR = 6
    PETSC_BIT_LOGICAL = 7
    PETSC_ENUM = 8
    PETSC_BOOL = 9
    PETSC_FLOAT128 = 10
    PETSC_OBJECT = 11
    PETSC_FUNCTION = 12
    PETSC_STRING = 13
    PETSC___FP16 = 14
    PETSC_STRUCT = 15
    PETSC_INT = 16
    PETSC_INT64 = 17
end

const Petsc64bitInt = Int64
const PetscLogDouble = Cdouble

@enum InsertMode NOT_SET_VALUES INSERT_VALUES ADD_VALUES MAX_VALUES MIN_VALUES INSERT_ALL_VALUES ADD_ALL_VALUES INSERT_BC_VALUES ADD_BC_VALUES

@enum NormType begin
    NORM_1 = 0
    NORM_2 = 1
    NORM_FROBENIUS = 2
    NORM_INFINITY = 3
    NORM_1_AND_2 = 4
end
const NORM_MAX = NORM_INFINITY

@enum MatAssemblyType MAT_FLUSH_ASSEMBLY = 1 MAT_FINAL_ASSEMBLY = 0

@enum MatFactorType begin
    MAT_FACTOR_NONE = 0
    MAT_FACTOR_LU = 1
    MAT_FACTOR_CHOLESKY = 2
    MAT_FACTOR_ILU = 3
    MAT_FACTOR_ICC = 4
    MAT_FACTOR_ILUDT = 5
end

@enum ScatterMode begin
    SCATTER_FORWARD = 0
    SCATTER_REVERSE = 1
    SCATTER_FORWARD_LOCAL = 2
    SCATTER_REVERSE_LOCAL = 3
end
const SCATTER_LOCAL = SCATTER_FORWARD_LOCAL

@enum PetscViewerFormat begin
    PETSC_VIEWER_DEFAULT
    PETSC_VIEWER_ASCII_MATLAB
    PETSC_VIEWER_ASCII_MATHEMATICA
    PETSC_VIEWER_ASCII_IMPL
    PETSC_VIEWER_ASCII_INFO
    PETSC_VIEWER_ASCII_INFO_DETAIL
    PETSC_VIEWER_ASCII_COMMON
    PETSC_VIEWER_ASCII_SYMMODU
    PETSC_VIEWER_ASCII_INDEX
    PETSC_VIEWER_ASCII_DENSE
    PETSC_VIEWER_ASCII_MATRIXMARKET
    PETSC_VIEWER_ASCII_VTK
    PETSC_VIEWER_ASCII_VTK_CELL
    PETSC_VIEWER_ASCII_VTK_COORDS
    PETSC_VIEWER_ASCII_PCICE
    PETSC_VIEWER_ASCII_PYTHON
    PETSC_VIEWER_ASCII_FACTOR_INFO
    PETSC_VIEWER_ASCII_LATEX
    PETSC_VIEWER_ASCII_XML
    PETSC_VIEWER_ASCII_GLVIS
    PETSC_VIEWER_ASCII_CSV
    PETSC_VIEWER_DRAW_BASIC
    PETSC_VIEWER_DRAW_LG
    PETSC_VIEWER_DRAW_LG_XRANGE
    PETSC_VIEWER_DRAW_CONTOUR
    PETSC_VIEWER_DRAW_PORTS
    PETSC_VIEWER_VTK_VTS
    PETSC_VIEWER_VTK_VTR
    PETSC_VIEWER_VTK_VTU
    PETSC_VIEWER_BINARY_MATLAB
    PETSC_VIEWER_NATIVE
    PETSC_VIEWER_HDF5_PETSC
    PETSC_VIEWER_HDF5_VIZ
    PETSC_VIEWER_HDF5_XDMF
    PETSC_VIEWER_HDF5_MAT
    PETSC_VIEWER_NOFORMAT
    PETSC_VIEWER_LOAD_BALANCE
end

@enum MatOperation begin
    MATOP_SET_VALUES = 0
    MATOP_GET_ROW = 1
    MATOP_RESTORE_ROW = 2
    MATOP_MULT = 3
    MATOP_MULT_ADD = 4
    MATOP_MULT_TRANSPOSE = 5
    MATOP_MULT_TRANSPOSE_ADD = 6
    MATOP_SOLVE = 7
    MATOP_SOLVE_ADD = 8
    MATOP_SOLVE_TRANSPOSE = 9
    MATOP_SOLVE_TRANSPOSE_ADD = 10
end

@enum DMBoundaryType begin
    DM_BOUNDARY_NONE = 0
    DM_BOUNDARY_GHOSTED = 1
    DM_BOUNDARY_MIRROR = 2
    DM_BOUNDARY_PERIODIC = 3
    DM_BOUNDARY_TWIST = 4
end

@enum DMDAStencilType begin
    DMDA_STENCIL_STAR = 0
    DMDA_STENCIL_BOX = 1
end


text:
#=
Copyright 2022 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

"""
    solveConjunction(b :: Blood, v1 :: Vessel, v2 :: Vessel)

Solve the non-linear system at the conjunction node between two vessels.
"""
function solveConjunction(b :: Blood, v1 :: Vessel, v2 :: Vessel)
    @fastmath @inbounds U0 = @SArray [v1.u[end],
         v2.u[1],
         sqrt(sqrt(v1.A[end])),
         sqrt(sqrt(v2.A[1]))]

    @inbounds k1 = v1.s_15_gamma[end]
    @inbounds k2 = v2.s_15_gamma[1]
    @inbounds k3 = b.rho
    k  = @SArray [k1, k2, k3]

    J = calculateJacobianConjunction(v1, v2, U0, k)
    U = newtonRaphson([v1, v2], J, U0, k, calculateWstarConjunction, calculateFConjunction)

    updateConjunction(U, v1, v2)
end

"""
    calculateJacobianConjunction(v1 :: Vessel, v2 :: Vessel, U, k)

Return the Jacobian for conjunction equations.
"""
function calculateJacobianConjunction(v1 :: Vessel, v2 :: Vessel, U, k)
    @fastmath @inbounds U33 = U[3]*U[3]*U[3]
    @fastmath @inbounds U43 = U[4]*U[4]*U[4]

    @fastmath @inbounds J13 =  4.0*k[1]
    @fastmath @inbounds J24 = -4.0*k[2]

    @fastmath @inbounds J31 =  U33*U[3]
    @fastmath @inbounds J32 = -U43*U[4]
    @fastmath @inbounds J33 =  4.0*U[1]*U33
    @fastmath @inbounds J34 = -4.0*U[2]*U43

    @fastmath @inbounds J41 =  k[3]*U[1]
    @fastmath @inbounds J42 = -k[3]*U[2]
    @fastmath @inbounds J43 =  2.0*v1.beta[end]*U[3]*v1.s_inv_A0[end]
    @fastmath @inbounds J44 = -2.0*v2.beta[1]*U[4]*v2.s_inv_A0[1]

    return @SArray [1.0 0.0 J13 0.0;
                    0.0 1.0 0.0 J24;
                    J31 J32 J33 J34;
                    J41 J42 J43 J44]
end

"""
    calculateWstarConjunction(U, k)

Return the Riemann invariants at the conjunction node.
"""
function calculateWstarConjunction(U, k)
    @fastmath @inbounds W1 = U[1] + 4.0*k[1]*U[3]
    @fastmath @inbounds W2 = U[2] - 4.0*k[2]*U[4]

    return @SArray [W1, W2]
end

"""
    calculateFConjunction(vessels :: Array{Vessel,1}, U, k, W)
"""
function calculateFConjunction(vessels :: Array{Vessel,1}, U, k, W)
    v1 = vessels[1]
    v2 = vessels[2]

    @fastmath @inbounds U32 = U[3]*U[3]
    @fastmath @inbounds U42 = U[4]*U[4]

    @fastmath @inbounds f1 = U[1] + 4.0*k[1]*U[3] - W[1]

    @fastmath @inbounds f2 = U[2] - 4.0*k[2]*U[4] - W[2]

    @fastmath @inbounds f3 = U[1]*U32*U32 - U[2]*U42*U42

    f4 = 0.5*k[3]*U[1]*U[1] + v1.beta[end]*(U32*v1.s_inv_A0[end] - 1.0) -
        (0.5*k[3]*U[2]*U[2] + v2.beta[1]*(U42*v2.s_inv_A0[1] - 1.0))

    return @SArray [f1, f2, f3, f4]
end

"""
    updateConjunction(U, v1 :: Vessel, v2 :: Vessel)

Update the values at the conjunction node for the two vessels.
"""
function updateConjunction(U, v1 :: Vessel, v2 :: Vessel)
    @inbounds v1.u[end] = U[1]
    @inbounds v2.u[1] = U[2]

    @fastmath @inbounds v1.A[end] = U[3]*U[3]*U[3]*U[3]
    @fastmath @inbounds v1.Q[end] = v1.u[end]*v1.A[end]

    @fastmath @inbounds v2.A[1] = U[4]*U[4]*U[4]*U[4]
    @fastmath @inbounds v2.Q[1] = v2.u[1]*v2.A[1]

    @inbounds v1.P[end] = pressure(v1.A[end], v1.A0[end], v1.beta[end], v1.Pext)
    @inbounds v2.P[1] = pressure(v2.A[1], v2.A0[1], v2.beta[1], v2.Pext)

    @fastmath @inbounds v1.c[end] = waveSpeed(v1.A[end], v1.gamma[end])
    @fastmath @inbounds v2.c[1] = waveSpeed(v2.A[1], v2.gamma[1])
end


text:
# TensorMap & Tensor:
# general tensor implementation with arbitrary symmetries
#==========================================================#
"""
    struct TensorMap{S<:IndexSpace, N₁, N₂, ...} <: AbstractTensorMap{S, N₁, N₂}

Specific subtype of [`AbstractTensorMap`](@ref) for representing tensor maps (morphisms in
a tensor category) whose data is stored in blocks of some subtype of `DenseMatrix`.
"""
struct TensorMap{S<:IndexSpace, N₁, N₂, G<:Sector, A<:Union{<:DenseMatrix,SectorDict{G,<:DenseMatrix}}, F₁, F₂} <: AbstractTensorMap{S, N₁, N₂}
    data::A
    codom::ProductSpace{S,N₁}
    dom::ProductSpace{S,N₂}
    rowr::SectorDict{G,FusionTreeDict{F₁,UnitRange{Int}}}
    colr::SectorDict{G,FusionTreeDict{F₂,UnitRange{Int}}}
    function TensorMap{S, N₁, N₂, G, A, F₁, F₂}(data::A,
                codom::ProductSpace{S,N₁}, dom::ProductSpace{S,N₂},
                rowr::SectorDict{G,FusionTreeDict{F₁,UnitRange{Int}}},
                colr::SectorDict{G,FusionTreeDict{F₂,UnitRange{Int}}}) where
                    {S<:IndexSpace, N₁, N₂, G<:Sector, A<:SectorDict{G,<:DenseMatrix},
                     F₁<:FusionTree{G,N₁}, F₂<:FusionTree{G,N₂}}
        eltype(valtype(data)) ⊆ field(S) ||
            @warn("eltype(data) = $(eltype(data)) ⊈ $(field(S)))", maxlog=1)
        new{S, N₁, N₂, G, A, F₁, F₂}(data, codom, dom, rowr, colr)
    end
    function TensorMap{S, N₁, N₂, Trivial, A, Nothing, Nothing}(data::A,
                codom::ProductSpace{S,N₁}, dom::ProductSpace{S,N₂}) where
                    {S<:IndexSpace, N₁, N₂, A<:DenseMatrix}
        eltype(data) ⊆ field(S) ||
            @warn("eltype(data) = $(eltype(data)) ⊈ $(field(S)))", maxlog=1)
        new{S, N₁, N₂, Trivial, A, Nothing, Nothing}(data, codom, dom)
    end
end

const Tensor{S<:IndexSpace, N, G<:Sector, A, F₁, F₂} = TensorMap{S, N, 0, G, A, F₁, F₂}
const TrivialTensorMap{S<:IndexSpace, N₁, N₂, A<:DenseMatrix} = TensorMap{S, N₁, N₂, Trivial, A, Nothing, Nothing}

# Basic methods for characterising a tensor:
#--------------------------------------------
codomain(t::TensorMap) = t.codom
domain(t::TensorMap) = t.dom

blocksectors(t::TrivialTensorMap) = (Trivial(),)
blocksectors(t::TensorMap) = keys(t.data)

Base.@pure storagetype(::Type{<:TensorMap{<:IndexSpace,N₁,N₂,Trivial,A}}) where
    {N₁,N₂,A<:DenseMatrix} = A
Base.@pure storagetype(::Type{<:TensorMap{<:IndexSpace,N₁,N₂,G,<:SectorDict{G,A}}}) where
    {N₁,N₂,G<:Sector,A<:DenseMatrix} = A

dim(t::TensorMap) = mapreduce(x->length(x[2]), +, blocks(t); init = 0)

# General TensorMap constructors
#--------------------------------
# with data
function TensorMap(data::DenseArray, codom::ProductSpace{S,N₁}, dom::ProductSpace{S,N₂}) where {S<:IndexSpace, N₁, N₂}
    if sectortype(S) === Trivial # For now, we can only accept array data for Trivial sectortype
        (d1, d2) = (dim(codom), dim(dom))
        if !(length(data) == d1*d2 || size(data) == (d1, d2) ||
            size(data) == (dims(codom)..., dims(dom)...))
            throw(DimensionMismatch())
        end
        data2 = reshape(data, (d1, d2))
        A = typeof(data2)
        return TensorMap{S, N₁, N₂, Trivial, A, Nothing, Nothing}(data2, codom, dom)
    else
        # TODO: allow to start from full data (a single DenseArray) and create the dictionary, in the first place for Abelian sectors, or for e.g. SU₂ using Wigner 3j symbols
        throw(SectorMismatch())
    end
end

function TensorMap(data::AbstractDict{<:Sector,<:DenseMatrix}, codom::ProductSpace{S,N₁}, dom::ProductSpace{S,N₂}) where {S<:IndexSpace, N₁, N₂}
    G = sectortype(S)
    G == keytype(data) || throw(SectorMismatch())
    G == Trivial && return TensorMap(data[Trivial()], codom, dom)
    F₁ = fusiontreetype(G, StaticLength(N₁))
    F₂ = fusiontreetype(G, StaticLength(N₂))
    rowr = SectorDict{G, FusionTreeDict{F₁, UnitRange{Int}}}()
    colr = SectorDict{G, FusionTreeDict{F₂, UnitRange{Int}}}()
    for c in blocksectors(codom ← dom)
        rowrc = FusionTreeDict{F₁, UnitRange{Int}}()
        colrc = FusionTreeDict{F₂, UnitRange{Int}}()
        offset1 = 0
        for s1 in sectors(codom)
            for f1 in fusiontrees(s1, c, map(isdual, codom.spaces))
                r = (offset1 + 1):(offset1 + dim(codom, s1))
                push!(rowrc, f1 => r)
                offset1 = last(r)
            end
        end
        offset2 = 0
        for s2 in sectors(dom)
            for f2 in fusiontrees(s2, c, map(isdual, dom.spaces))
                r = (offset2 + 1):(offset2 + dim(dom, s2))
                push!(colrc, f2 => r)
                offset2 = last(r)
            end
        end
        (haskey(data, c) && size(data[c]) == (offset1, offset2)) ||
            throw(DimensionMismatch())
        push!(rowr, c=>rowrc)
        push!(colr, c=>colrc)
    end
    if !isreal(G) && eltype(valtype(data)) <: Real
        b = valtype(data)(undef, (0,0))
        V = typeof(complex(b))
        K = keytype(data)
        data2 = SectorDict{K,V}((c=>complex(b)) for (c,b) in data)
        A = typeof(data2)
        return TensorMap{S, N₁, N₂, G, A, F₁, F₂}(data2, codom, dom, rowr, colr)
    else
        V = valtype(data)
        K = keytype(data)
        data2 = SectorDict{K,V}((c=>d) for (c,d) in data)
        A = typeof(data2)
        return TensorMap{S, N₁, N₂, G, A, F₁, F₂}(data2, codom, dom, rowr, colr)
    end
end

# without data: generic constructor from callable:
function TensorMap(f, codom::ProductSpace{S,N₁}, dom::ProductSpace{S,N₂}) where {S<:IndexSpace, N₁, N₂}
    G = sectortype(S)
    if G == Trivial
        d1 = dim(codom)
        d2 = dim(dom)
        data = f((d1,d2))
        A = typeof(data)
        return TensorMap{S, N₁, N₂, Trivial, A, Nothing, Nothing}(data, codom, dom)
    else
        F₁ = fusiontreetype(G, StaticLength(N₁))
        F₂ = fusiontreetype(G, StaticLength(N₂))
        # TODO: the current approach is not very efficient and somewhat wasteful
        sampledata = f((1,1))
        if !isreal(G) && eltype(sampledata) <: Real
            A = typeof(complex(sampledata))
        else
            A = typeof(sampledata)
        end
        data = SectorDict{G,A}()
        rowr = SectorDict{G, FusionTreeDict{F₁, UnitRange{Int}}}()
        colr = SectorDict{G, FusionTreeDict{F₂, UnitRange{Int}}}()
        for c in blocksectors(codom ← dom)
            rowrc = FusionTreeDict{F₁, UnitRange{Int}}()
            colrc = FusionTreeDict{F₂, UnitRange{Int}}()
            offset1 = 0
            for s1 in sectors(codom)
                for f1 in fusiontrees(s1, c, map(isdual, codom.spaces))
                    r = (offset1 + 1):(offset1 + dim(codom, s1))
                    push!(rowrc, f1 => r)
                    offset1 = last(r)
                end
            end
            dim1 = offset1
            offset2 = 0
            for s2 in sectors(dom)
                for f2 in fusiontrees(s2, c, map(isdual, dom.spaces))
                    r = (offset2 + 1):(offset2 + dim(dom, s2))
                    push!(colrc, f2 => r)
                    offset2 = last(r)
                end
            end
            dim2 = offset2
            push!(data, c=>f((dim1, dim2)))
            push!(rowr, c=>rowrc)
            push!(colr, c=>colrc)
        end
        return TensorMap{S, N₁, N₂, G, SectorDict{G,A}, F₁, F₂}(data, codom, dom, rowr, colr)
    end
end

TensorMap(f,
            ::Type{T},
            codom::ProductSpace{S},
            dom::ProductSpace{S}) where {S<:IndexSpace, T<:Number} =
    TensorMap(d->f(T, d), codom, dom)

TensorMap(::Type{T},
            codom::ProductSpace{S},
            dom::ProductSpace{S}) where {S<:IndexSpace, T<:Number} =
    TensorMap(d->Array{T}(undef, d), codom, dom)

TensorMap(::UndefInitializer,
            ::Type{T},
            codom::ProductSpace{S},
            dom::ProductSpace{S}) where {S<:IndexSpace, T<:Number} =
    TensorMap(d->Array{T}(undef, d), codom, dom)

TensorMap(::UndefInitializer,
            codom::ProductSpace{S},
            dom::ProductSpace{S}) where {S<:IndexSpace} =
    TensorMap(undef, Float64, codom, dom)

function TensorMap(I::UniformScaling,
                    ::Type{T},
                    codom::ProductSpace{S},
                    dom::ProductSpace{S}) where {S<:IndexSpace, T<:Number}
    Base.depwarn("`TensorMap(I, T, codomain, domain)`  using `LinearAlgebra.I` is deprecated, use `id([A,], space)` for creating the identity tensor on a space, and `isomorphism([A,], codomain, domain)` to construct a fixed invertible map between two isomorphic spaces. When `spacetype(domain)<:EuclideanSpace`, one can also use `unitary([A,], codomain, domain)` which is then equivalent to `isomorphism`. `A<:AbstractMatrix` is the type of storage, and is by default chosen equal to `Matrix{Float64}`.", ((Base.Core).Typeof(TensorMap)).name.mt.name)
    isomorphism(Matrix{T}, codom, dom)
end
TensorMap(I::UniformScaling,
            codom::ProductSpace{S},
            dom::ProductSpace{S}) where {S<:IndexSpace} =
    TensorMap(I, Float64, codom, dom)

TensorMap(::Type{T},
            codom::TensorSpace{S},
            dom::TensorSpace{S}) where {T<:Number, S<:IndexSpace} =
    TensorMap(T, convert(ProductSpace, codom), convert(ProductSpace, dom))

TensorMap(dataorf, codom::TensorSpace{S}, dom::TensorSpace{S}) where {S<:IndexSpace} =
    TensorMap(dataorf, convert(ProductSpace, codom), convert(ProductSpace, dom))

TensorMap(dataorf, ::Type{T},
            codom::TensorSpace{S},
            dom::TensorSpace{S}) where {T<:Number, S<:IndexSpace} =
    TensorMap(dataorf, T, convert(ProductSpace, codom), convert(ProductSpace, dom))

TensorMap(codom::TensorSpace{S}, dom::TensorSpace{S}) where {S<:IndexSpace} =
    TensorMap(Float64, convert(ProductSpace, codom), convert(ProductSpace, dom))

TensorMap(dataorf, T::Type{<:Number}, P::TensorMapSpace{S}) where {S<:IndexSpace} =
    TensorMap(dataorf, T, codomain(P), domain(P))

TensorMap(dataorf, P::TensorMapSpace{S}) where {S<:IndexSpace} =
    TensorMap(dataorf, codomain(P), domain(P))

TensorMap(T::Type{<:Number}, P::TensorMapSpace{S}) where {S<:IndexSpace} =
    TensorMap(T, codomain(P), domain(P))

TensorMap(P::TensorMapSpace{S}) where {S<:IndexSpace} = TensorMap(codomain(P), domain(P))

Tensor(dataorf, T::Type{<:Number}, P::TensorSpace{S}) where {S<:IndexSpace} =
    TensorMap(dataorf, T, P, one(P))

Tensor(dataorf, P::TensorSpace{S}) where {S<:IndexSpace} = TensorMap(dataorf, P, one(P))

Tensor(T::Type{<:Number}, P::TensorSpace{S}) where {S<:IndexSpace} = TensorMap(T, P, one(P))

Tensor(P::TensorSpace{S}) where {S<:IndexSpace} = TensorMap(P, one(P))

# Efficient copy constructors
#-----------------------------
function Base.copy(t::TrivialTensorMap{S, N₁, N₂, A}) where {S, N₁, N₂, A}
    return TrivialTensorMap{S, N₁, N₂, A}(copy(t.data), t.codom, t.dom)
end
function Base.copy(t::TensorMap{S, N₁, N₂, G, A, F₁, F₂}) where {S, N₁, N₂, G, A, F₁, F₂}
    return TensorMap{S, N₁, N₂, G, A, F₁, F₂}(deepcopy(t.data), t.codom, t.dom, t.rowr, t.colr)
end

# Similar
#---------
Base.similar(t::AbstractTensorMap, T::Type, codomain::VectorSpace, domain::VectorSpace) =
    similar(t, T, codomain←domain)
Base.similar(t::AbstractTensorMap, codomain::VectorSpace, domain::VectorSpace) =
    similar(t, codomain←domain)

Base.similar(t::AbstractTensorMap{S}, ::Type{T},
                P::TensorMapSpace{S} = (domain(t) → codomain(t))) where {T,S} =
    TensorMap(d->similarstoragetype(t, T)(undef, d), P)
Base.similar(t::AbstractTensorMap{S}, ::Type{T}, P::TensorSpace{S}) where {T,S} =
    Tensor(d->similarstoragetype(t, T)(undef, d), P)
Base.similar(t::AbstractTensorMap{S},
                P::TensorMapSpace{S} = (domain(t) → codomain(t))) where {S} =
    TensorMap(d->storagetype(t)(undef, d), P)
Base.similar(t::AbstractTensorMap{S}, P::TensorSpace{S}) where {S} =
    Tensor(d->storagetype(t)(undef, d), P)

function Base.complex(t::AbstractTensorMap)
    if eltype(t) <: Complex
        return t
    elseif t.data isa AbstractArray
        return TensorMap(complex(t.data), codomain(t), domain(t))
    else
        data = SectorDict(c=>complex(d) for (c,d) in t.data)
        return TensorMap(data, codomain(t), domain(t))
    end
end

# Conversion between TensorMap and Dict, for read and write purpose
#------------------------------------------------------------------
function Base.convert(::Type{Dict}, t::AbstractTensorMap)
    d = Dict{Symbol,Any}()
    d[:codomain] = repr(codomain(t))
    d[:domain] = repr(domain(t))
    data = Dict{String,Any}()
    for (c,b) in blocks(t)
        data[repr(c)] = Array(b)
    end
    d[:data] = data
    return d
end
function Base.convert(::Type{TensorMap}, d::Dict{Symbol,Any})
    codomain = eval(Meta.parse(d[:codomain]))
    domain = eval(Meta.parse(d[:domain]))
    data = SectorDict(eval(Meta.parse(c))=>b for (c,b) in d[:data])
    return TensorMap(data, codomain, domain)
end

# Getting and setting the data
#------------------------------
hasblock(t::TrivialTensorMap, ::Trivial) = true
hasblock(t::TensorMap, s::Sector) = haskey(t.data, s)

block(t::TrivialTensorMap, ::Trivial) = t.data
function block(t::TensorMap, s::Sector)
    sectortype(t) == typeof(s) || throw(SectorMismatch())
    A = valtype(t.data)
    if haskey(t.data, s)
        return t.data[s]
    else # at least one of the two matrix dimensions will be zero
        return storagetype(t)(undef, (blockdim(codomain(t),s), blockdim(domain(t), s)))
    end
end

blocks(t::TensorMap{<:IndexSpace,N₁,N₂,Trivial}) where {N₁,N₂} =
    SingletonDict(Trivial()=>t.data)
blocks(t::TensorMap) = t.data

fusiontrees(t::TrivialTensorMap) = ((nothing, nothing),)
fusiontrees(t::TensorMap) = TensorKeyIterator(t.rowr, t.colr)

@inline function Base.getindex(t::TensorMap{<:IndexSpace,N₁,N₂,G},
                                sectors::Tuple{Vararg{G}}) where {N₁,N₂,G<:Sector}

    FusionStyle(G) isa Abelian ||
        throw(SectorMismatch("Indexing with sectors only possible if abelian"))
    s1 = TupleTools.getindices(sectors, codomainind(t))
    s2 = TupleTools.getindices(sectors, domainind(t))
    c1 = length(s1) == 0 ? one(G) : (length(s1) == 1 ? s1[1] : first(⊗(s1...)))
    @boundscheck begin
        c2 = length(s2) == 0 ? one(G) : (length(s2) == 1 ? s2[1] : first(⊗(s1...)))
        c2 == c1 || throw(SectorMismatch())
        hassector(codomain(t), s1) && hassector(domain(t), s2)
    end
    f1 = FusionTree(s1,c1)
    f2 = FusionTree(s2,c1)
    @inbounds begin
        return t[f1,f2]
    end
end
@propagate_inbounds Base.getindex(t::TensorMap, sectors::Tuple) =
    t[map(sectortype(t), sectors)]

@inline function Base.getindex(t::TensorMap{<:IndexSpace,N₁,N₂,G},
        f1::FusionTree{G,N₁}, f2::FusionTree{G,N₂}) where {N₁,N₂,G<:Sector}

    c = f1.coupled
    @boundscheck begin
        c == f2.coupled || throw(SectorMismatch())
        haskey(t.rowr[c], f1) || throw(SectorMismatch())
        haskey(t.colr[c], f2) || throw(SectorMismatch())
    end
    @inbounds begin
        d = (dims(codomain(t), f1.uncoupled)..., dims(domain(t), f2.uncoupled)...)
        return sreshape(StridedView(t.data[c])[t.rowr[c][f1], t.colr[c][f2]], d)
    end
end
@propagate_inbounds Base.setindex!(t::TensorMap{<:IndexSpace,N₁,N₂,G}, v, f1::FusionTree{G,N₁}, f2::FusionTree{G,N₂}) where {N₁,N₂,G<:Sector} = copyto!(getindex(t, f1, f2), v)

# For a tensor with trivial symmetry, allow no argument indexing
@inline Base.getindex(t::TrivialTensorMap) =
    sreshape(StridedView(t.data), (dims(codomain(t))..., dims(domain(t))...))
@inline Base.setindex!(t::TrivialTensorMap, v) = copyto!(getindex(t), v)

# For a tensor with trivial symmetry, fusiontrees returns (nothing,nothing)
@inline Base.getindex(t::TrivialTensorMap, ::Tuple{Nothing,Nothing}) = getindex(t)
@inline Base.setindex!(t::TrivialTensorMap, v, ::Tuple{Nothing,Nothing}) = setindex!(t, v)

# For a tensor with trivial symmetry, allow direct indexing
@inline function Base.getindex(t::TrivialTensorMap, I::Vararg{Int})
    data = t[]
    @boundscheck checkbounds(data, I...)
    @inbounds v = data[I...]
    return v
end
@inline function Base.setindex!(t::TrivialTensorMap, v, I::Vararg{Int})
    data = t[]
    @boundscheck checkbounds(data, I...)
    @inbounds data[I...] = v
    return v
end

# Show
#------
function Base.summary(t::TensorMap)
    print("TensorMap(", codomain(t), " ← ", domain(t), ")")
end
function Base.show(io::IO, t::TensorMap{S}) where {S<:IndexSpace}
    if get(io, :compact, false)
        print(io, "TensorMap(", codomain(t), " ← ", domain(t), ")")
        return
    end
    println(io, "TensorMap(", codomain(t), " ← ", domain(t), "):")
    if sectortype(S) == Trivial
        Base.print_array(io, t[])
        println(io)
    elseif FusionStyle(sectortype(S)) isa Abelian
        for (f1,f2) in fusiontrees(t)
            println(io, "* Data for sector ", f1.uncoupled, " ← ", f2.uncoupled, ":")
            Base.print_array(io, t[f1,f2])
            println(io)
        end
    else
        for (f1,f2) in fusiontrees(t)
            println(io, "* Data for fusiontree ", f1, " ← ", f2, ":")
            Base.print_array(io, t[f1,f2])
            println(io)
        end
    end
end

# Conversion and promotion:
#---------------------------
# TODO


text:
push!(LOAD_PATH,"../../src/")
import UNSflow

alphadef = UNSflow.ConstDef(6. *pi/180)

hdef = UNSflow.ConstDef(0.)

udef = UNSflow.SinDef(0., 0.5, 0.25, 0.)

full_kinem = UNSflow.KinemDef(alphadef, hdef, udef)

pvt = 0.25

geometry = "FlatPlate"

lespcrit = [0.2;]

surf = UNSflow.TwoDSurf(geometry, pvt, full_kinem, lespcrit)

curfield = UNSflow.TwoDFlowField(UNSflow.ConstDef(1.), UNSflow.ConstDef(0.))

dtstar = 0.05

t_tot = 60.

nsteps = Int(round(t_tot/dtstar))+1

startflag = 0

writeflag = 1

writeInterval = t_tot/5.

#delvort = delSpalart(500, 12, 1e-5)
delvort = UNSflow.delNone()

mat, surf, curfield = UNSflow.ldvmLin(surf, curfield, nsteps, dtstar,startflag, writeflag, writeInterval, delvort)

UNSflow.makeVortPlots2D()

UNSflow.makeForcePlots2D()

UNSflow.cleanWrite()


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2016 Varnerlab
# Robert Frederick Smith School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850
#
# 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: calculate_transcription_rates
# Description: Calculate the transcriptional rate array at time t
# Generated on: 2016-11-02T07:00:01.333
#
# Input arguments:
# t::Float64 => Current time value (scalar)
# x::Array{Float64,1} => State array (number_of_species x 1)
# data_dictionary::Dict{AbstractString,Any} => Dictionary holding model parameters
#
# Output arguments:
# transcription_rate_array::Array{Float64,1} => Transcriptional rate array (number_of_genes x 1) at time t
# ----------------------------------------------------------------------------------- #
function calculate_transcription_rates(t::Float64,x::Array{Float64,1},data_dictionary::Dict{AbstractString,Any})

	# Alias the species -
	gene_AP1 = x[1]
	gene_AhR = x[2]
	gene_CD11b = x[3]
	gene_CD14 = x[4]
	gene_CD38 = x[5]
	gene_CEBPa = x[6]
	gene_E2F = x[7]
	gene_EGR1 = x[8]
	gene_GFI1 = x[9]
	gene_IRF1 = x[10]
	gene_OCT1 = x[11]
	gene_OCT4 = x[12]
	gene_P21 = x[13]
	gene_P47Phox = x[14]
	gene_PPARg = x[15]
	gene_PU1 = x[16]
	gene_Trigger = x[17]
	gene_cRAF = x[18]

	# Initialize the transcription rate -
	transcription_rate_array = zeros(18)
	KSAT = data_dictionary["saturation_constant_transcription"]
	kcat_transcription = data_dictionary["kcat_transcription"]
	rnapII_concentration = data_dictionary["rnapII_concentration"]
	average_transcript_length = data_dictionary["average_transcript_length"]
	gene_coding_length_array = data_dictionary["gene_coding_length_array"]

	# Populate the transcription rate array -
	# Gene: gene_AP1
	gene_length = gene_coding_length_array[1]
	scale_factor = (average_transcript_length/gene_length)
	transcription_rate_array[1] = scale_factor*kcat_transcription*(rnapII_concentration)*((gene_AP1)/(KSAT+gene_AP1))

	# Gene: gene_AhR
	gene_length = gene_coding_length_array[2]
	scale_factor = (average_transcript_length/gene_length)
	transcription_rate_array[2] = scale_factor*kcat_transcription*(rnapII_concentration)*((gene_AhR)/(KSAT+gene_AhR))

	# Gene: gene_CD11b
	gene_length = gene_coding_length_array[3]
	scale_factor = (average_transcript_length/gene_length)
	transcription_rate_array[3] = scale_factor*kcat_transcription*(rnapII_concentration)*((gene_CD11b)/(KSAT+gene_CD11b))

	# Gene: gene_CD14
	gene_length = gene_coding_length_array[4]
	scale_factor = (average_transcript_length/gene_length)
	transcription_rate_array[4] = scale_factor*kcat_transcription*(rnapII_concentration)*((gene_CD14)/(KSAT+gene_CD14))

	# Gene: gene_CD38
	gene_length = gene_coding_length_array[5]
	scale_factor = (average_transcript_length/gene_length)
	transcription_rate_array[5] = scale_factor*kcat_transcription*(rnapII_concentration)*((gene_CD38)/(KSAT+gene_CD38))

	# Gene: gene_CEBPa
	gene_length = gene_coding_length_array[6]
	scale_factor = (average_transcript_length/gene_length)
	transcription_rate_array[6] = scale_factor*kcat_transcription*(rnapII_concentration)*((gene_CEBPa)/(KSAT+gene_CEBPa))

	# Gene: gene_E2F
	gene_length = gene_coding_length_array[7]
	scale_factor = (average_transcript_length/gene_length)
	transcription_rate_array[7] = scale_factor*kcat_transcription*(rnapII_concentration)*((gene_E2F)/(KSAT+gene_E2F))

	# Gene: gene_EGR1
	gene_length = gene_coding_length_array[8]
	scale_factor = (average_transcript_length/gene_length)
	transcription_rate_array[8] = scale_factor*kcat_transcription*(rnapII_concentration)*((gene_EGR1)/(KSAT+gene_EGR1))

	# Gene: gene_GFI1
	gene_length = gene_coding_length_array[9]
	scale_factor = (average_transcript_length/gene_length)
	transcription_rate_array[9] = scale_factor*kcat_transcription*(rnapII_concentration)*((gene_GFI1)/(KSAT+gene_GFI1))

	# Gene: gene_IRF1
	gene_length = gene_coding_length_array[10]
	scale_factor = (average_transcript_length/gene_length)
	transcription_rate_array[10] = scale_factor*kcat_transcription*(rnapII_concentration)*((gene_IRF1)/(KSAT+gene_IRF1))

	# Gene: gene_OCT1
	gene_length = gene_coding_length_array[11]
	scale_factor = (average_transcript_length/gene_length)
	transcription_rate_array[11] = scale_factor*kcat_transcription*(rnapII_concentration)*((gene_OCT1)/(KSAT+gene_OCT1))

	# Gene: gene_OCT4
	gene_length = gene_coding_length_array[12]
	scale_factor = (average_transcript_length/gene_length)
	transcription_rate_array[12] = scale_factor*kcat_transcription*(rnapII_concentration)*((gene_OCT4)/(KSAT+gene_OCT4))

	# Gene: gene_P21
	gene_length = gene_coding_length_array[13]
	scale_factor = (average_transcript_length/gene_length)
	transcription_rate_array[13] = scale_factor*kcat_transcription*(rnapII_concentration)*((gene_P21)/(KSAT+gene_P21))

	# Gene: gene_P47Phox
	gene_length = gene_coding_length_array[14]
	scale_factor = (average_transcript_length/gene_length)
	transcription_rate_array[14] = scale_factor*kcat_transcription*(rnapII_concentration)*((gene_P47Phox)/(KSAT+gene_P47Phox))

	# Gene: gene_PPARg
	gene_length = gene_coding_length_array[15]
	scale_factor = (average_transcript_length/gene_length)
	transcription_rate_array[15] = scale_factor*kcat_transcription*(rnapII_concentration)*((gene_PPARg)/(KSAT+gene_PPARg))

	# Gene: gene_PU1
	gene_length = gene_coding_length_array[16]
	scale_factor = (average_transcript_length/gene_length)
	transcription_rate_array[16] = scale_factor*kcat_transcription*(rnapII_concentration)*((gene_PU1)/(KSAT+gene_PU1))

	# Gene: gene_Trigger
	gene_length = gene_coding_length_array[17]
	scale_factor = (average_transcript_length/gene_length)
	transcription_rate_array[17] = scale_factor*kcat_transcription*(rnapII_concentration)*((gene_Trigger)/(KSAT+gene_Trigger))

	# Gene: gene_cRAF
	gene_length = gene_coding_length_array[18]
	scale_factor = (average_transcript_length/gene_length)
	transcription_rate_array[18] = scale_factor*kcat_transcription*(rnapII_concentration)*((gene_cRAF)/(KSAT+gene_cRAF))

	# return transcription_rate_array -
	return transcription_rate_array
end

# ----------------------------------------------------------------------------------- #
# Function: calculate_background_transcription_rates
# Description: Calculate the leak transcriptional rate array at time t
# Generated on: 2016-11-02T07:00:01.334
#
# Input arguments:
# t::Float64 => Current time value (scalar)
# x::Array{Float64,1} => State array (number_of_species x 1)
# data_dictionary::Dict{AbstractString,Any} => Dictionary holding model parameters
#
# Output arguments:
# background_transcription_rate_array::Array{Float64,1} => Background transcriptional rate array (number_of_genes x 1) at time t
# ----------------------------------------------------------------------------------- #
function calculate_background_transcription_rates(t::Float64,x::Array{Float64,1},transcription_rate_array::Array{Float64,1},data_dictionary::Dict{AbstractString,Any})
	return zeros(length(x))
end

# ----------------------------------------------------------------------------------- #
# Function: calculate_translation_rates
# Description: Calculate the translation rate array at time t
# Generated on: 2016-11-02T07:00:01.334
#
# Input arguments:
# t::Float64 => Current time value (scalar)
# x::Array{Float64,1} => State array (number_of_species x 1)
# data_dictionary::Dict{AbstractString,Any} => Dictionary holding model parameters
#
# Output arguments:
# translation_rate_array::Array{Float64,1} => Translation rate array (number_of_genes x 1) at time t
# ----------------------------------------------------------------------------------- #
function calculate_translation_rates(t::Float64,x::Array{Float64,1},data_dictionary::Dict{AbstractString,Any})

	# Alias the species -
	mRNA_gene_AP1 = x[19]
	mRNA_gene_AhR = x[20]
	mRNA_gene_CD11b = x[21]
	mRNA_gene_CD14 = x[22]
	mRNA_gene_CD38 = x[23]
	mRNA_gene_CEBPa = x[24]
	mRNA_gene_E2F = x[25]
	mRNA_gene_EGR1 = x[26]
	mRNA_gene_GFI1 = x[27]
	mRNA_gene_IRF1 = x[28]
	mRNA_gene_OCT1 = x[29]
	mRNA_gene_OCT4 = x[30]
	mRNA_gene_P21 = x[31]
	mRNA_gene_P47Phox = x[32]
	mRNA_gene_PPARg = x[33]
	mRNA_gene_PU1 = x[34]
	mRNA_gene_Trigger = x[35]
	mRNA_gene_cRAF = x[36]

	# Initialize the translation rate -
	translation_rate_array = zeros(18)
	KSAT = data_dictionary["saturation_constant_translation"]
	kcat_translation = data_dictionary["kcat_translation"]
	ribosome_concentration = data_dictionary["ribosome_concentration"]
	average_protein_length = data_dictionary["average_protein_length"]
	protein_coding_length_array = data_dictionary["protein_coding_length_array"]

	# Populate the translation rate array -
	# Transcript: mRNA_gene_AP1
	protein_length = protein_coding_length_array[1]
	scale_factor = (average_protein_length/protein_length)
	translation_rate_array[1] = scale_factor*kcat_translation*(ribosome_concentration)*((mRNA_gene_AP1)/(KSAT+mRNA_gene_AP1))

	# Transcript: mRNA_gene_AhR
	protein_length = protein_coding_length_array[2]
	scale_factor = (average_protein_length/protein_length)
	translation_rate_array[2] = scale_factor*kcat_translation*(ribosome_concentration)*((mRNA_gene_AhR)/(KSAT+mRNA_gene_AhR))

	# Transcript: mRNA_gene_CD11b
	protein_length = protein_coding_length_array[3]
	scale_factor = (average_protein_length/protein_length)
	translation_rate_array[3] = scale_factor*kcat_translation*(ribosome_concentration)*((mRNA_gene_CD11b)/(KSAT+mRNA_gene_CD11b))

	# Transcript: mRNA_gene_CD14
	protein_length = protein_coding_length_array[4]
	scale_factor = (average_protein_length/protein_length)
	translation_rate_array[4] = scale_factor*kcat_translation*(ribosome_concentration)*((mRNA_gene_CD14)/(KSAT+mRNA_gene_CD14))

	# Transcript: mRNA_gene_CD38
	protein_length = protein_coding_length_array[5]
	scale_factor = (average_protein_length/protein_length)
	translation_rate_array[5] = scale_factor*kcat_translation*(ribosome_concentration)*((mRNA_gene_CD38)/(KSAT+mRNA_gene_CD38))

	# Transcript: mRNA_gene_CEBPa
	protein_length = protein_coding_length_array[6]
	scale_factor = (average_protein_length/protein_length)
	translation_rate_array[6] = scale_factor*kcat_translation*(ribosome_concentration)*((mRNA_gene_CEBPa)/(KSAT+mRNA_gene_CEBPa))

	# Transcript: mRNA_gene_E2F
	protein_length = protein_coding_length_array[7]
	scale_factor = (average_protein_length/protein_length)
	translation_rate_array[7] = scale_factor*kcat_translation*(ribosome_concentration)*((mRNA_gene_E2F)/(KSAT+mRNA_gene_E2F))

	# Transcript: mRNA_gene_EGR1
	protein_length = protein_coding_length_array[8]
	scale_factor = (average_protein_length/protein_length)
	translation_rate_array[8] = scale_factor*kcat_translation*(ribosome_concentration)*((mRNA_gene_EGR1)/(KSAT+mRNA_gene_EGR1))

	# Transcript: mRNA_gene_GFI1
	protein_length = protein_coding_length_array[9]
	scale_factor = (average_protein_length/protein_length)
	translation_rate_array[9] = scale_factor*kcat_translation*(ribosome_concentration)*((mRNA_gene_GFI1)/(KSAT+mRNA_gene_GFI1))

	# Transcript: mRNA_gene_IRF1
	protein_length = protein_coding_length_array[10]
	scale_factor = (average_protein_length/protein_length)
	translation_rate_array[10] = scale_factor*kcat_translation*(ribosome_concentration)*((mRNA_gene_IRF1)/(KSAT+mRNA_gene_IRF1))

	# Transcript: mRNA_gene_OCT1
	protein_length = protein_coding_length_array[11]
	scale_factor = (average_protein_length/protein_length)
	translation_rate_array[11] = scale_factor*kcat_translation*(ribosome_concentration)*((mRNA_gene_OCT1)/(KSAT+mRNA_gene_OCT1))

	# Transcript: mRNA_gene_OCT4
	protein_length = protein_coding_length_array[12]
	scale_factor = (average_protein_length/protein_length)
	translation_rate_array[12] = scale_factor*kcat_translation*(ribosome_concentration)*((mRNA_gene_OCT4)/(KSAT+mRNA_gene_OCT4))

	# Transcript: mRNA_gene_P21
	protein_length = protein_coding_length_array[13]
	scale_factor = (average_protein_length/protein_length)
	translation_rate_array[13] = scale_factor*kcat_translation*(ribosome_concentration)*((mRNA_gene_P21)/(KSAT+mRNA_gene_P21))

	# Transcript: mRNA_gene_P47Phox
	protein_length = protein_coding_length_array[14]
	scale_factor = (average_protein_length/protein_length)
	translation_rate_array[14] = scale_factor*kcat_translation*(ribosome_concentration)*((mRNA_gene_P47Phox)/(KSAT+mRNA_gene_P47Phox))

	# Transcript: mRNA_gene_PPARg
	protein_length = protein_coding_length_array[15]
	scale_factor = (average_protein_length/protein_length)
	translation_rate_array[15] = scale_factor*kcat_translation*(ribosome_concentration)*((mRNA_gene_PPARg)/(KSAT+mRNA_gene_PPARg))

	# Transcript: mRNA_gene_PU1
	protein_length = protein_coding_length_array[16]
	scale_factor = (average_protein_length/protein_length)
	translation_rate_array[16] = scale_factor*kcat_translation*(ribosome_concentration)*((mRNA_gene_PU1)/(KSAT+mRNA_gene_PU1))

	# Transcript: mRNA_gene_Trigger
	protein_length = protein_coding_length_array[17]
	scale_factor = (average_protein_length/protein_length)
	translation_rate_array[17] = scale_factor*kcat_translation*(ribosome_concentration)*((mRNA_gene_Trigger)/(KSAT+mRNA_gene_Trigger))

	# Transcript: mRNA_gene_cRAF
	protein_length = protein_coding_length_array[18]
	scale_factor = (average_protein_length/protein_length)
	translation_rate_array[18] = scale_factor*kcat_translation*(ribosome_concentration)*((mRNA_gene_cRAF)/(KSAT+mRNA_gene_cRAF))

	# return translation array -
	return translation_rate_array
end

# ----------------------------------------------------------------------------------- #
# Function: calculate_mRNA_degradation_rates
# Description: Calculate the mRNA degradation rate array at time t
# Generated on: 2016-11-02T07:00:01.337
#
# Input arguments:
# t::Float64 => Current time value (scalar)
# x::Array{Float64,1} => State array (number_of_species x 1)
# data_dictionary::Dict{AbstractString,Any} => Dictionary holding model parameters
#
# Output arguments:
# mRNA_degradation_rate_array::Array{Float64,1} => mRNA degradation rate array (number_of_genes x 1) at time t
# ----------------------------------------------------------------------------------- #
function calculate_mRNA_degradation_rates(t::Float64,x::Array{Float64,1},data_dictionary::Dict{AbstractString,Any})

	# Alias the species -
	mRNA_gene_AP1 = x[19]
	mRNA_gene_AhR = x[20]
	mRNA_gene_CD11b = x[21]
	mRNA_gene_CD14 = x[22]
	mRNA_gene_CD38 = x[23]
	mRNA_gene_CEBPa = x[24]
	mRNA_gene_E2F = x[25]
	mRNA_gene_EGR1 = x[26]
	mRNA_gene_GFI1 = x[27]
	mRNA_gene_IRF1 = x[28]
	mRNA_gene_OCT1 = x[29]
	mRNA_gene_OCT4 = x[30]
	mRNA_gene_P21 = x[31]
	mRNA_gene_P47Phox = x[32]
	mRNA_gene_PPARg = x[33]
	mRNA_gene_PU1 = x[34]
	mRNA_gene_Trigger = x[35]
	mRNA_gene_cRAF = x[36]

	# Initialize the degrdation array -
	degradation_rate_array = zeros(18)
	mRNA_degrdation_constant = data_dictionary["degradation_constant_mRNA"]

	# Calculate the degradation_rate_array -
	degradation_rate_array[1] = (mRNA_degrdation_constant)*mRNA_gene_AP1
	degradation_rate_array[2] = (mRNA_degrdation_constant)*mRNA_gene_AhR
	degradation_rate_array[3] = (mRNA_degrdation_constant)*mRNA_gene_CD11b
	degradation_rate_array[4] = (mRNA_degrdation_constant)*mRNA_gene_CD14
	degradation_rate_array[5] = (mRNA_degrdation_constant)*mRNA_gene_CD38
	degradation_rate_array[6] = (mRNA_degrdation_constant)*mRNA_gene_CEBPa
	degradation_rate_array[7] = (mRNA_degrdation_constant)*mRNA_gene_E2F
	degradation_rate_array[8] = (mRNA_degrdation_constant)*mRNA_gene_EGR1
	degradation_rate_array[9] = (mRNA_degrdation_constant)*mRNA_gene_GFI1
	degradation_rate_array[10] = (mRNA_degrdation_constant)*mRNA_gene_IRF1
	degradation_rate_array[11] = (mRNA_degrdation_constant)*mRNA_gene_OCT1
	degradation_rate_array[12] = (mRNA_degrdation_constant)*mRNA_gene_OCT4
	degradation_rate_array[13] = (mRNA_degrdation_constant)*mRNA_gene_P21
	degradation_rate_array[14] = (mRNA_degrdation_constant)*mRNA_gene_P47Phox
	degradation_rate_array[15] = (mRNA_degrdation_constant)*mRNA_gene_PPARg
	degradation_rate_array[16] = (mRNA_degrdation_constant)*mRNA_gene_PU1
	degradation_rate_array[17] = (mRNA_degrdation_constant)*mRNA_gene_Trigger
	degradation_rate_array[18] = (mRNA_degrdation_constant)*mRNA_gene_cRAF

	# return the degrdation rate array -
	return degradation_rate_array
end

# ----------------------------------------------------------------------------------- #
# Function: calculate_protein_degradation_rates
# Description: Calculate the protein degradation rate array at time t
# Generated on: 2016-11-02T07:00:01.34
#
# Input arguments:
# t::Float64 => Current time value (scalar)
# x::Array{Float64,1} => State array (number_of_species x 1)
# data_dictionary::Dict{AbstractString,Any} => Dictionary holding model parameters
#
# Output arguments:
# protein_degradation_rate_array::Array{Float64,1} => protein degradation rate array (number_of_proteins x 1) at time t
# ----------------------------------------------------------------------------------- #
function calculate_protein_degradation_rates(t::Float64,x::Array{Float64,1},data_dictionary::Dict{AbstractString,Any})

	# Alias the species -
	protein_gene_AP1 = x[37]
	protein_gene_AhR = x[38]
	protein_gene_CD11b = x[39]
	protein_gene_CD14 = x[40]
	protein_gene_CD38 = x[41]
	protein_gene_CEBPa = x[42]
	protein_gene_E2F = x[43]
	protein_gene_EGR1 = x[44]
	protein_gene_GFI1 = x[45]
	protein_gene_IRF1 = x[46]
	protein_gene_OCT1 = x[47]
	protein_gene_OCT4 = x[48]
	protein_gene_P21 = x[49]
	protein_gene_P47Phox = x[50]
	protein_gene_PPARg = x[51]
	protein_gene_PU1 = x[52]
	protein_gene_Trigger = x[53]
	protein_gene_cRAF = x[54]

	# Initialize the degrdation array -
	degradation_rate_array = zeros(18)
	protein_degrdation_constant = data_dictionary["degradation_constant_protein"]

	# Calculate the degradation_rate_array -
	degradation_rate_array[1] = (protein_degrdation_constant)*protein_gene_AP1
	degradation_rate_array[2] = (protein_degrdation_constant)*protein_gene_AhR
	degradation_rate_array[3] = (protein_degrdation_constant)*protein_gene_CD11b
	degradation_rate_array[4] = (protein_degrdation_constant)*protein_gene_CD14
	degradation_rate_array[5] = (protein_degrdation_constant)*protein_gene_CD38
	degradation_rate_array[6] = (protein_degrdation_constant)*protein_gene_CEBPa
	degradation_rate_array[7] = (protein_degrdation_constant)*protein_gene_E2F
	degradation_rate_array[8] = (protein_degrdation_constant)*protein_gene_EGR1
	degradation_rate_array[9] = (protein_degrdation_constant)*protein_gene_GFI1
	degradation_rate_array[10] = (protein_degrdation_constant)*protein_gene_IRF1
	degradation_rate_array[11] = (protein_degrdation_constant)*protein_gene_OCT1
	degradation_rate_array[12] = (protein_degrdation_constant)*protein_gene_OCT4
	degradation_rate_array[13] = (protein_degrdation_constant)*protein_gene_P21
	degradation_rate_array[14] = (protein_degrdation_constant)*protein_gene_P47Phox
	degradation_rate_array[15] = (protein_degrdation_constant)*protein_gene_PPARg
	degradation_rate_array[16] = (protein_degrdation_constant)*protein_gene_PU1
	degradation_rate_array[17] = (protein_degrdation_constant)*protein_gene_Trigger
	degradation_rate_array[18] = (protein_degrdation_constant)*protein_gene_cRAF

	# return the degrdation rate array -
	return degradation_rate_array
end


text:
using JuMP, Cbc

pModel = Model(Cbc.Optimizer)

setupTime = [45 45 45 45; 60 60 60 60; 100 100 100 100]
price = [50 50 50 50; 70 70 70 70; 120 120 120 120]
cost = [18 18 18 18; 15 15 15 15; 25 25 25 25]
rate = [5 5 5 5; 4 4 4 4 ;2 2 2 2]
demand = [400 600 200 800;
          240 440 100 660;
          80 120 40 100]

@variable(pModel, time[1:3,1:4]<= 435)
@variable(pModel, setUp[1:3,1:4], binary=true)
@variable(pModel, production[1:3,1:4])
@variable(pModel, inventory[1:3,1:4])

@constraint(pModel, time.<=480*setUp)
@constraint(pModel, sum(((setupTime.*setUp)+time),dims = 1) .<= 480)
@constraint(pModel, production .== (time.*rate))

d = 1
    @constraint(pModel, inventory[:,d] .== production[:,d]-demand[:,d])
for d = 2:4
    @constraint(pModel, inventory[:,d] .== inventory[:,d-1]+production[:,d]-demand[:,d])
end

@constraint(pModel, inventory .>= 0)
@constraint(pModel, inventory[:,4] .<= 0)

@objective(pModel, Max, sum(price.*production)-sum(cost.*inventory))

result = JuMP.optimize!(pModel)

max_stars_repo_name: valenmgama/optimization-tec

max_issues_repo_name: valenmgama/optimization-tec

max_forks_repo_name: valenmgama/optimization-tec


text:
module TestLogConfig

using ..Fixtures
using FactCheck
using Log4jl
import Log4jl: LoggerConfig, name, level, isadditive, isenabled, parent, parent!, level!

facts("Logger Configuration") do
    context("specifies `root` logger configuration by default") do
        lc = LoggerConfig()
        @fact level(lc) --> Log4jl.LOG4JL_DEFAULT_STATUS_LEVEL
        @fact name(lc) --> isempty
        @fact lc --> not(isadditive)
        @fact parent(lc) --> isnull
    end

    context("created with name and level parameters") do
        lc = LoggerConfig(Fixtures.TESTNAME)
        @fact level(lc) --> Log4jl.LOG4JL_DEFAULT_STATUS_LEVEL
        @fact name(lc) --> Fixtures.TESTNAME
        @fact lc --> not(isadditive)
    end

    context("could have a different log level and a parent") do
        lc = LoggerConfig(Fixtures.TESTNAME, level=Fixtures.TESTLEVELN)
        @fact level(lc) --> Fixtures.TESTLEVEL
        level!(lc, Log4jl.Level.OFF)
        @fact level(lc) --> Log4jl.Level.OFF
        plc = LoggerConfig()
        parent!(lc, plc)
        @fact parent(lc) --> not(isnull)
        @fact get(parent(lc)) --> plc
    end
end

end


text:
using BCTRNN
using DiffEqSensitivity
using OrdinaryDiffEq
import DiffEqFlux: FastChain, FastDense
import Flux: ClipValue, ADAM

# Not in Project.toml
using Plots
gr()

include("sine_wave_dataloader.jl")

function train_sine_fc(epochs, solver=Tsit5(); model_size=2,
  sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)),
  T=Float32, mtkize=false, lr=0.02, kwargs...)

  train_dl = generate_2d_data(T)

  f_in = 2
  f_out = 1

  im = BCTRNN.InputAllToAll()
  sm = BCTRNN.SynsFullyConnected()
  om = BCTRNN.OutputAll()
  wiring = BCTRNN.WiringConfig(f_in, model_size, im,sm,om)

  model = FastChain(BCTRNN.Mapper(f_in),
                    BCTRNN.LTCGJ(wiring, solver, sensealg; T, mtkize, kwargs...),
                    FastDense(wiring.n_out, f_out))

  hs = []
  for (k,v) in wiring.matrices
    push!(hs, heatmap(v, title=k))
  end
  display(plot(hs..., layout=length(hs)))

  cb = BCTRNN.MyCallback(T; ecb=mycb, nepochs=epochs, nsamples=length(train_dl))
  #opt = GalacticOptim.Flux.Optimiser(ClipValue(0.5), ADAM(0.02))
  opt = BCTRNN.ClampBoundOptim(BCTRNN.get_bounds(model,T)..., ClipValue(T(1.0)), ADAM(T(lr)))
  BCTRNN.optimize(model, BCTRNN.loss_seq, cb, opt, train_dl, epochs, T), model
end

@time train_sine_fc(100, Tsit5(); model_size=7)

function train_sine(epochs, solver=Tsit5();
  n_sensory=2, n_inter=5, n_command=5, n_motor=1,
  sensory_inter=2, inter_command=3, command_command=2, command_motor=2,
  mtkize=false, lr=0.02,
  sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)), T=Float32, kwargs...)

  train_dl = generate_2d_data(T)

  f_in = 2
  f_out = 1

  im = BCTRNN.InputAllToFirstN(n_sensory)
  sm = BCTRNN.SynsNCP(; n_sensory, n_inter, n_command, n_motor,
    sensory_inter, inter_command, command_command, command_motor
  )
  om = BCTRNN.OutputAll()
  wiring = BCTRNN.WiringConfig(f_in, im,sm,om)

  model = FastChain(BCTRNN.Mapper(f_in),
                    BCTRNN.LTCGJ(wiring, solver, sensealg; T, mtkize, kwargs...),
                    FastDense(wiring.n_out, f_out))

  hs = []
  for (k,v) in wiring.matrices
    push!(hs, heatmap(v, title=k))
  end
  display(plot(hs..., layout=length(hs)))

  cb = BCTRNN.MyCallback(T; ecb=mycb, nepochs=epochs, nsamples=length(train_dl))
  #opt = GalacticOptim.Flux.Optimiser(ClipValue(0.5), ADAM(0.02))
  opt = BCTRNN.ClampBoundOptim(BCTRNN.get_bounds(model,T)..., ClipValue(T(1.0)), ADAM(T(lr)))
  BCTRNN.optimize(model, BCTRNN.loss_seq, cb, opt, train_dl, epochs, T), model
end

@time train_sine(100, Tsit5())


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2019 Varnerlab
# Robert Frederick School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850

# 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.
# ----------------------------------------------------------------------------------- #

# includes -
include("Types.jl")
include("Kinetics.jl")
include("Control.jl")
include("Inputs.jl")
include("Data.jl")
include("SolveBalances.jl")
include("Balances.jl")
include("Utility.jl")
include("Error.jl")

# system packages - these are required to be installed -
# check - if they are installed, all is good with the world.
# if not then install them
using LinearAlgebra # pre-installed w/Julia
using Statistics    # pre-installed w/Julia
using Pkg           # pre-installed w/Julia
installed_package_set = keys(Pkg.installed())

# Do we have DifferentialEquations?
if (in("DifferentialEquations",installed_package_set) == false)
    Pkg.add("DifferentialEquations")
end

# Do we have DelimitedFiles?
if (in("DelimitedFiles",installed_package_set) == false)
    Pkg.add("DelimitedFiles")
end

# Do we have JSON?
if (in("JSON",installed_package_set) == false)
    Pkg.add("JSON")
end

# Do we have ProgressMeter?
if (in("ProgressMeter",installed_package_set) == false)
    Pkg.add("ProgressMeter")
end

if (in("Optim",installed_package_set) == false)
    Pkg.add("Optim")
end

if (in("DataFrames",installed_package_set) == false)
    Pkg.add("DataFrames")
end

if (in("CSV",installed_package_set) == false)
    Pkg.add("CSV")
end

if (in("Interpolations",installed_package_set) == false)
    Pkg.add("Interpolations")
end

#if (in("PyPlot",installed_package_set) == false)
#    Pkg.add("PyPlot")
#end

# Ok to use now ...
using DifferentialEquations
using DelimitedFiles
using JSON
using ProgressMeter
using Optim
using DataFrames
using CSV
using Interpolations
#using PyPlot

# List any custom includes here ...
# ...


text:
function pointbefore(firstpoints::P, scale=ReflectionScale) where P
    initialx = reflectback(first.(firstpoints[1:2])..., scale=scale)
    initialys = thiele4(firstpoints..., initialx)
    if isa(eltype(firstpoints), Tuple)
        initialpoint = (initialx, initialys[1:length(firstpoints[1])-1]...,)
    else
        initialpoint = [initialx, initialys[1:length(firstpoints[1])-1]...,]
    end
    return convert(eltype(firstpoints), initialpoint)
end

function pointafter(lastpoints::P, scale=ReflectionScale) where P
    finalx = reflectforward(first.(lastpoints[end-1:end])..., scale=scale)
    finalys = thiele4(lastpoints..., finalx)
    finalpoint = (finalx, finalys...,)
    if isa(eltype(lastpoints), Tuple)
        finalpoint = (finalx, finalys[1:length(lastpoints[1])-1]...,)
    else
        finalpoint = [finalx, finalys[1:length(lastpoints[1])-1]...,]
    end
    return convert(eltype(lastpoints), finalpoint)
end

allfinite(result) = all(Iterators.flatten(map(x->isfinite.(x),result)))


text:
#===  ===#
function writeinfo(
    model::Union{elmod2d,acmod2d},
    sou::Array{MTsource};
    PartextPath="SimuPara.txt",
    ParbinPath="SimuPara.bin")

    sounumber = length(sou)
#=== Simu&OutputInfo.txt ===#
    if PartextPath != nothing
        if PartextPath[end-2:end] != "txt"
            error("Check file extensions. It has to be in 'txt' format. ")
        else
            fid = open(PartextPath,"a+")
            wMT(fid,model)
            wST(fid,sounumber,sou,model.medium.pkf)
            wOPT(fid)
            close(fid)
        end
    end
#=== SimuPara.bin ===#
    # disabled, in garbage folder
    # wMSB(model,sou,sounumber,ParbinPath)
end
#======#

#===  ===#
function writeinfo(
    model::Union{elmod2d,acmod2d},
    sou::Array{MTsource},
    rec::receiver,
    OptCpnt::Array{String};
    PartextPath="SimuPara.txt",
    ParbinPath="SimuPara.bin")

    sounumber = length(sou)
#=== Simu&OutputInfo.txt ===#
    if PartextPath != nothing
        if PartextPath[end-2:end] != "txt"
            error("Check file extensions. It has to be in 'txt' format. ")
        else
            fid = open(PartextPath,"a+")
            wMT(fid,model)
            wST(fid,sounumber,sou,model.medium.pkf)
            wOPT(fid,model,rec,OptCpnt)
            close(fid)
        end
    end
# #=== SimuPara.bin ===#
#     wMSB(model,sou,sounumber,ParbinPath)
end
#======#

#===  ===#
function writeinfo(
    model::Union{elmod2d,acmod2d},
    sou::Array{MTsource},
    Slices::Union{StepRange{Int64,Int64},Array{Int64,1}},
    OptCpnt::Array{String};
    PartextPath="SimuPara.txt",
    ParbinPath="SimuPara.bin")

    sounumber = length(sou)
#=== Simu&OutputInfo.txt ===#
    if PartextPath != nothing
        if PartextPath[end-2:end] != "txt"
            error("Check file extensions. It has to be in 'txt' format. ")
        else
            fid = open(PartextPath,"a+")
            wMT(fid,model)
            wST(fid,sounumber,sou,model.medium.pkf)
            wOPT(fid,model,Slices,OptCpnt)
            close(fid)
        end
    end
# #=== SimuPara.bin ===#
#     wMSB(model,sou,sounumber,ParbinPath)
end
#======#

#===  ===#
function writeinfo(
    model::Union{elmod2d,acmod2d},
    sou::Array{MTsource},
    Slices::Union{StepRange{Int64,Int64},Array{Int64,1}},
    rec::receiver,
    OptCpnt::Array{String};
    PartextPath="SimuPara.txt",
    ParbinPath="SimuPara.bin")

    sounumber = length(sou)
#=== Simu&OutputInfo.txt ===#
    if PartextPath != nothing
        if PartextPath[end-2:end] != "txt"
            error("Check file extensions. It has to be in 'txt' format. ")
        else
            fid = open(PartextPath,"a+")
            wMT(fid,model)
            wST(fid,sounumber,sou,model.medium.pkf)
            wOPT(fid,model,Slices,rec,OptCpnt)
            close(fid)
        end
    end
# #=== SimuPara.bin ===#
#     wMSB(model,sou,sounumber,ParbinPath)
end
#======#

#===  ===#
function writeinfo(
    model::Union{elmod2d},
    sou::Array{MTsource},
    Slices::Union{StepRange{Int64,Int64},Array{Int64,1}};
    PartextPath="SimuPara.txt",
    ParbinPath="SimuPara.bin")

    sounumber = length(sou)
#=== Simu&OutputInfo.txt ===#
    if PartextPath != nothing
        if PartextPath[end-2:end] != "txt"
            error("Check file extensions. It has to be in 'txt' format. ")
        else
            fid = open(PartextPath,"a+")
            wMT(fid,model)
            wST(fid,sounumber,sou,model.medium.pkf)
            wOPT(fid,model,Slices)
            close(fid)
        end
    end
# #=== SimuPara.bin ===#
#     wMSB(model,sou,sounumber,ParbinPath)
end
#======#

#===  ===#
function writeinfo(
    model::Union{acmod2d},
    sou::Array{MTsource},
    Slices::Union{StepRange{Int64,Int64},Array{Int64,1}};
    PartextPath="SimuPara.txt",
    ParbinPath="SimuPara.bin")

    sounumber = length(sou)
#=== Simu&OutputInfo.txt ===#
    if PartextPath != nothing
        if PartextPath[end-2:end] != "txt"
            error("Check file extensions. It has to be in 'txt' format. ")
        else
            fid = open(PartextPath,"a+")
            wMT(fid,model)
            wST(fid,sounumber,sou,model.medium.pkf)
            wOPT(fid,model,Slices)
            close(fid)
        end
    end
# #=== SimuPara.bin ===#
#     wMSB(model,sou,sounumber,ParbinPath)
end
#======#


text:
"""
	MMCK(μ, λ, c, k)

Tạo mô hình M/M/c/k
"""
struct MMCK{T} <:AbstractMMCK
	λ::T
	μ::T
	c::Union{T, Integer}
	k::Union{T, Integer}
	r::Union{T, Real}
	ρ::Union{T, Real}
	function MMCK(λ, μ, c, k)
		T = Union{typeof(μ), typeof(λ)}
		new{T}(λ, μ, c,	k, λ/μ,	λ/μ/c)
	end
end

function pn(m::MMCK, n::Integer)
	r = m.r
	k = m.k
	c = m.c
	if n == 0
		1 / (
		     sum(r^i / factorial(i) for i = 0:c-1) +
		     sum(r^i / (c^(i-c) * factorial(c)) for i = c:k)
		     )
	elseif n < c
		pn(m, 0) * r^n / factorial(n)
	elseif n <= k
		pn(m, 0) * r^n / (c^(n-c) * factorial(c))
	else 0
	end
end

function L(m::MMCK)
	sum(i * pn(m, i) for i=0:m.k)
end

function Lq(m::MMCK)
	sum((i - m.c) * pn(m, i) for i=m.c:m.k)
end

function W(m::MMCK)
	L(m) / λe(m)
end

function Wq(m::MMCK)
	Lq(m) / λe(m)
end


text:
# Imports
using DelimitedFiles
using Combinatorics
using Permutations

# Run a simple search through words data
# Getting the word dataset
# words_path = "./english-words/words_alpha.txt"
words_path = "./oxford_words.txt"
words_alpha = readdlm(words_path);

# Declare sets of variables
letters = (('b', 'e', 'p'), ('u', 'i', 'n'), ('c', 't', 'a'), ('l', 'm', 'o'));
# letters = (('h', 'i', 'g'), ('a', 's', 'b'), ('m', 'x', 'c'), ('o', 'e', 't'));
# Construct the rules for the words search
# Letters that should be in the words
all_letters = collect(Iterators.flatten(letters));

# Step-1: Get the words containing only the basis letters
min_word_size = 10

# Declare empty word list
contains_letters_list = String[]
for word in words_alpha
    split_word = collect(word)
    word_length = size(split_word)
    # Add rules to filter only required words
    if all(in(all_letters).(split_word)) && (word_length[1] >= min_word_size)
        # Add it to the new word list
        push!(contains_letters_list, word)
    end
end

# Letter pairings that should not be in the words
bad_combos = String[]
for lett in letters
    for l in lett
        push!(bad_combos, string(l,l))
    end
    for perm in collect(permutations(lett,2))
        push!(bad_combos, join(perm))
    end
end

# Step-2: Get the candidate words

candidate_words = String[]
for word in contains_letters_list
    word_split = collect(word)
    word_length = (size(word_split))[1]
    combo_check = true
    for ii = 1:word_length-1
        combo_check = combo_check && !in(bad_combos).(join(word_split[ii:ii+1]))
    end
    if combo_check == true
        push!(candidate_words, word)
    end
end

# Now get candidate solutions
candidate_solutions = []
for first_word in candidate_words
    first_word_split = collect(first_word)
    last_letter = (first_word_split)[end]
    for second_word in candidate_words
        second_word_split = collect(second_word)
        first_letter = (second_word_split)[1]
        if last_letter == first_letter
            two_word_letters = union(first_word_split, second_word_split)
            if all(in(two_word_letters).(all_letters))
                push!(candidate_solutions, (first_word, second_word))
            end
        end
    end
end

println(candidate_solutions)


text:
# This file includes macro tools to define new components types

#=
    Returns true if `ex` is of the form

```julia
julia> ex = :(
       struct Foo <: Bar
       x::Int
       end
       )
```
Note that `Foo` is a subtype of `Bar`.
=#
issubtype(ex::Expr) = (nex = ex.args[2]; nex isa Expr && nex.head == :(<:))

#=
    Returns true if `ex` is of the form

```julia
julia> ex = :(
       struct Foo{T} <: Bar
       x::T
       end
       )
:(struct Foo{T} <: Bar
      #= REPL[17]:3 =#
      x::T
  end)
```
Note that `Foo` is a parametric type.
=#
isparametric(ex::Expr) = (nex = ex.args[2].args[1]; nex isa Expr && nex.head == :curly)

#=
   Returns true if `ex` has default values.

```julia
ex = :(x::T = 1.)
```
=#
hasdefault(ex::Expr) = ex.head == :(=)

#=
    Appends the `nex` to `ex`. A mininmum working example:

```julia
julia> ex = :(
       struct Foo <: Bar
       x::Int
       end
       )
:(struct Foo <: Bar
      #= REPL[15]:3 =#
      x::Int
  end)

julia> Causal.Components.ComponentsBase.appendex!(ex, nex)
:(struct Foo{T} <: Bar
      #= REPL[15]:3 =#
      x::Int
      y::T = 1.0
  end)
```
=#
function appendex!(ex::Expr, nex::Expr)
    ex.args[3].head == :block || error("Invalid expression $ex")

    # Append expression
    push!(ex.args[3].args, nex)
    # field::T
    if nex.head == :(::)
        typesymbol = nex.args[2]
    # field::T = default or field = T
    elseif nex.head == :(=)
        lhs = nex.args[1]
        typesymbol  =
            # field::T
            if lhs isa Expr && lhs.head == :(::)
                typesymbol = lhs.args[2]
            # field
            elseif lhs isa Symbol
                 typesymbol = lhs
            else
                error("Invalid usage.")
            end
    end
    appendtype!(ex, typesymbol)
    ex
end

function appendtype!(ex::Expr, typesymbol::Symbol)
    if isparametric(ex)
        push!(ex.args[2].args[1].args, typesymbol)
    else
        ex.args[2].args[1] = :($(ex.args[2].args[1]){$typesymbol})
    end
    ex
end

# function _append_common_fields!(ex, newbody, newparamtypes)
#     # Append body
#     body = ex.args[3]
#     append!(body.args, newbody.args)

#     # Append struct type parameters
#     name = ex.args[2]
#     if name isa Expr && name.head === :(<:)
#         name = name.args[1]
#     end

#     if name isa Expr && name.head === :curly
#         append!(name.args, newparamtypes)
#     elseif name isa Symbol
#         ex.args[2] = Expr(:curly, name, newparamtypes...)  # parametrize ex
#     end
# end

# function deftype(ex)
#     # Check ex head
#     ex isa Expr && ex.head == :struct || error("Invalid source defition")

#     # Get struct name
#     name = ex.args[2]
#     if name isa Expr && name.head === :(<:)
#         name = name.args[1]
#     end

#     # Process struct body
#     body = ex.args[3]
#     kwargs = Expr(:parameters)
#     callargs = Symbol[]

#     _def!(body, kwargs, callargs)

#     # struct has no fields
#     isempty(kwargs.args) && return quote
#         Base.@__doc__($(esc(ex)))
#     end

#     if name isa Symbol
#         return quote
#             Base.@__doc__($(esc(ex)))
#             $(esc(name))($kwargs) = $(esc(name))($(callargs...))
#         end
#     elseif name isa Expr && name.head === :curly
#         _name = name.args[1]
#         _param_types = name.args[2:end]
#         __param_types = [_type_ isa Symbol ? _type_  : _type_.args[1] for _type_ in _param_types]
#         return quote
#             Base.@__doc__($(esc(ex)))
#             $(esc(_name))($kwargs) = $(esc(_name))($(callargs...))
#             $(esc(_name)){$(esc.(__param_types)...)}($kwargs) where {$(esc.(_param_types)...)} =
#                 $(esc(_name)){$(esc.(__param_types)...)}($(callargs...))
#         end
#     end
# end

# function _def!(body, kwargs, callargs)
#     for i in 1 : length(body.args)
#         bodyex = body.args[i]
#         if bodyex isa Symbol # var
#             push!(kwargs.args, bodyex)
#             push!(callargs, bodyex)
#         elseif bodyex isa Expr
#             if bodyex.head === :(=)
#                 rhs = bodyex.args[2]
#                 lhs = bodyex.args[1]
#                 if lhs isa Expr && lhs.head === :(::) # var::T = 1
#                     var = lhs.args[1]
#                 elseif lhs isa Symbol # var = 1
#                     var = lhs
#                 elseif lhs isa Expr && lhs.head == :call # inner constructors
#                     continue
#                 end
#                 push!(kwargs.args, Expr(:kw, var, esc(rhs)))
#                 push!(callargs, var)
#                 body.args[i] = lhs
#             elseif bodyex.head === :(::)  # var::T
#                 var = bodyex.args[1]
#                 push!(kwargs.args, var)
#                 push!(callargs, var)
#             end
#         end
#     end
# end


text:
@testset "Client" begin

  @test sizeof(Jib.Client.HEADTYPE) == 4

  @test Jib.Client.MAX_LEN < typemax(Jib.Client.HEADTYPE)

  # isascii
  @test Jib.Client.isascii([0x80, 0x79, 0x79], 1)

  @test !Jib.Client.isascii([0x80, 0x80, 0x79], 1)

  # buffer
  buf = Jib.Client.buffer(true)
  @test ismarked(buf)
  reset(buf)
  @test String(take!(buf)) == Jib.Client.API_SIGN * "\0\0\0\0"

  # write_one
  buf = Jib.Client.buffer(true)
  write(buf, "ABC")
  bo = IOBuffer()
  Jib.Client.write_one(bo, buf)

  @test buf.size == 0
  @test String(take!(bo)) == Jib.Client.API_SIGN * "\0\0\0\x03ABC"

  # Round trip
  buf = Jib.Client.buffer(false)
  write(buf, "ABC")

  bo = IOBuffer()
  Jib.Client.write_one(bo, buf)

  seekstart(bo) # Rewind
  @test String(Jib.Client.read_one(bo)) == "ABC"
  @test eof(bo)
  @test bo.size == 7

end


text:
function QuantumOptics.expect(Obs::Union{GraphOperator, SparseOperator},
                     net::NeuralNetwork, sampler::Sampler=FullSumSampler())

    oprob = ObservablesProblem(NeuralQuantum.input_type(net), Obs)
    ov          = state(oprob, net)
    mosampl     = multithread(sampler);
    oitercache  = SamplingCache(oprob)
    moitercache = [SamplingCache(oprob) for i=1:Threads.nthreads()]
    moc         = init_sampler!(mosampl, net, ov)

    sample_net = sampler isa FullSumSampler ? sample_network_wholespace! : sample_network!
    #println(sample_net)
    init_sampler!(mosampl, net, ov, moc)
    #Threads.@threads
    for i=1:Threads.nthreads()
        s = sampler_list(mosampl)[i]
        n = moc.nets[i]
        v = moc.σs[i]
        _ic = moitercache[i]
        zero!(_ic)
        while true
            sample_net(_ic, oprob, n, v)
            !samplenext!(v, s, n, moc.caches[i]) && break
        end
    end
    zero!(oitercache)
    for ic=moitercache
        add!(oitercache, ic)
    end
    ObsEval = evaluation_post_sampling(oitercache)
    if sampler isa FullSumSampler
        return ObsEval.ObsAve[1], ObsEval.ObsVals[1]
    else
        return ObsEval.ObsAve[1], ObsEval.ObsVals[1]
    end
end


text:
doc"""
    Beta(α,β)

The *Beta distribution* has probability density function

$f(x; \alpha, \beta) = \frac{1}{B(\alpha, \beta)}
 x^{\alpha - 1} (1 - x)^{\beta - 1}, \quad x \in [0, 1]$

The Beta distribution is related to the [`Gamma`](:func:`Gamma`) distribution via the
property that if $X \sim \operatorname{Gamma}(\alpha)$ and $Y \sim \operatorname{Gamma}
(\beta)$ independently, then $X / (X + Y) \sim \operatorname{Beta}(\alpha, \beta)$.

```julia
Beta()        # equivalent to Beta(1.0, 1.0)
Beta(a)       # equivalent to Beta(a, a)
Beta(a, b)    # Beta distribution with shape parameters a and b

params(d)     # Get the parameters, i.e. (a, b)
```

External links

* [Beta distribution on Wikipedia](http://en.wikipedia.org/wiki/Beta_distribution)

"""
immutable Beta <: ContinuousUnivariateDistribution
    α::Float64
    β::Float64

    function Beta(α::Real, β::Real)
        @check_args(Beta, α > zero(α) && β > zero(β))
        new(α, β)
    end
    Beta(α::Real) = Beta(α, α)
    Beta() = new(1.0, 1.0)
end

@distr_support Beta 0.0 1.0

#### Parameters

params(d::Beta) = (d.α, d.β)

#### Statistics

mean(d::Beta) = ((α, β) = params(d); α / (α + β))

function mode(d::Beta)
    (α, β) = params(d)
    (α > 1.0 && β > 1.0) || error("mode is defined only when α > 1 and β > 1.")
    return (α - 1.0) / (α + β - 2.0)
end

modes(d::Beta) = [mode(d)]

function var(d::Beta)
    (α, β) = params(d)
    s = α + β
    return (α * β) / (abs2(s) * (s + 1.0))
end

meanlogx(d::Beta) = ((α, β) = params(d); digamma(α) - digamma(α + β))

varlogx(d::Beta) = ((α, β) = params(d); trigamma(α) - trigamma(α + β))
stdlogx(d::Beta) = sqrt(varlogx(d))

function skewness(d::Beta)
    (α, β) = params(d)
    if α == β
        return 0.0
    else
        s = α + β
        (2.0 * (β - α) * sqrt(s + 1.0)) / ((s + 2.0) * sqrt(α * β))
    end
end

function kurtosis(d::Beta)
    α, β = params(d)
    s = α + β
    p = α * β
    6.0 * (abs2(α - β) * (s + 1.0) - p * (s + 2.0)) / (p * (s + 2.0) * (s + 3.0))
end

function entropy(d::Beta)
    α, β = params(d)
    s = α + β
    lbeta(α, β) - (α - 1.0) * digamma(α) - (β - 1.0) * digamma(β) +
        (s - 2.0) * digamma(s)
end

#### Evaluation

@_delegate_statsfuns Beta beta α β

gradlogpdf(d::Beta, x::Float64) =
    ((α, β) = params(d); 0.0 <= x <= 1.0 ? (α - 1.0) / x - (β - 1.0) / (1 - x) : 0.0)

#### Sampling

rand(d::Beta) = StatsFuns.Rmath.betarand(d.α, d.β)

#### Fit model

# TODO: add MLE method (should be similar to Dirichlet)

# This is a moment-matching method (not MLE)
#
function fit{T<:Real}(::Type{Beta}, x::AbstractArray{T})
    x_bar = mean(x)
    v_bar = varm(x, x_bar)
    α = x_bar * (((x_bar * (1.0 - x_bar)) / v_bar) - 1.0)
    β = (1.0 - x_bar) * (((x_bar * (1.0 - x_bar)) / v_bar) - 1.0)
    Beta(α, β)
end


text:
@test begin
    b = Box(Point(-1, -1), Point(1, 1))
    [Box(Point(-1, -1), Point(0, 0)),
     Box(Point(-1, 0), Point(0, 1)),
     Box(Point(0, -1), Point(1, 0)),
     Box(Point(0, 0), Point(1, 1))] == split(b)
end

@test begin
    b = Box(Point(-1, -1), Point(1, 1))
    [Segment(Point(-1, -1), Point(-1, 1)),
     Segment(Point(-1, 1), Point(1, 1)),
     Segment(Point(1, 1), Point(1, -1)),
     Segment(Point(1, -1), Point(-1, -1))] == surface(b)
end

@test begin
    b = Box(Point(0, 0, 0), Point(1, 1, 1))
    [Box(Point(0, 0, 0), Point(0.5, 0.5, 0.5)),
     Box(Point(0, 0.5, 0), Point(0.5, 1, 0.5)),
     Box(Point(0.5, 0, 0), Point(1, 0.5, 0.5)),
     Box(Point(0.5, 0.5, 0), Point(1, 1, 0.5)),
     Box(Point(0, 0, 0.5), Point(0.5, 0.5, 1)),
     Box(Point(0, 0.5, 0.5), Point(0.5, 1, 1)),
     Box(Point(0.5, 0, 0.5), Point(1, 0.5, 1)),
     Box(Point(0.5, 0.5, 0.5), Point(1, 1, 1))] == split(b)
end

@test begin
    x = [-1, 2]
    y = [-2, 3]
    z = [-3, 4]

    b = Box(Point(-1, -2, -3), Point(2, 3, 4))
    faces0 = [[f(0, 0), f(1, 0), f(0, 1), f(1, 1)] for f in surface(b)]

    pts = [Point(xi, yi, zi) for zi in z for yi in y for xi in x]

    faces = [[1, 2, 3, 4], [5, 6, 7, 8], [1, 3, 7, 5], [2, 4, 6, 8],
             [1, 2, 5, 6], [3, 4, 7, 8]]

    faces = [pts[f] for f in faces]

    function match_face(points0, points)
        for p0 in points0
            found = false
            for p in points
                found = norm(p0 - p) < 1E-13
                found && break
            end
            !found && return false
        end
        true
    end

    for f0 in faces0
        found = false
        for f in faces
            found = match_face(f0, f)
            found && break
        end
        @assert found
    end
    true
end

max_stars_repo_name: MiroK/fractional-dim


text:
using CImGui
using CImGui.CSyntax
using CImGui.CSyntax.CStatic
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using CImGui.GLFWBackend.GLFW
using CImGui.OpenGLBackend.ModernGL
using Printf

using GR
using LaTeXStrings

const glsl_version = 150
GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac

function draw(phi)
    formulas = (
        L"- \frac{{\hbar ^2}}{{2m}}\frac{{\partial ^2 \psi (x,t)}}{{\partial x^2 }} + U(x)\psi (x,t) = i\hbar \frac{{\partial \psi (x,t)}}{{\partial t}}",
        L"\zeta \left({s}\right) := \sum_{n=1}^\infty \frac{1}{n^s} \quad \sigma = \Re(s) > 1",
        L"\zeta \left({s}\right) := \frac{1}{\Gamma(s)} \int_{0}^\infty \frac{x^{s-1}}{e^x-1} dx" )

    selntran(0)
    settextfontprec(232, 3)
    settextalign(2, 3)
    settextcolorind(0) # white
    chh = 0.036

    clearws()
    setcharheight(chh)
    setcharup(sin(phi), cos(phi))
    y = 0.6
    for s in formulas
        mathtex(0.5, y, s)
        tbx, tby = inqmathtex(0.5, y, s)
        fillarea(tbx, tby)
        y -= 0.2
    end
    updatews()
end

# setup GLFW error callback
error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
GLFW.SetErrorCallback(error_callback)

# create window
window = GLFW.CreateWindow(960, 720, "Demo")
@assert window != C_NULL
GLFW.MakeContextCurrent(window)
GLFW.SwapInterval(1)  # enable vsync

# setup Dear ImGui context
ctx = CImGui.CreateContext()

# create texture for image drawing
img_width, img_height = 500, 500
image_id = ImGui_ImplOpenGL3_CreateImageTexture(img_width, img_height)

ENV["GKS_WSTYPE"] = "100"
image = rand(GLuint, img_width, img_height) # allocate the memory
mem = Printf.@sprintf("%p",pointer(image))
conid = Printf.@sprintf("!%dx%d@%s.mem",  img_width, img_height, mem[3:end])

# setup Platform/Renderer bindings
ImGui_ImplGlfw_InitForOpenGL(window, true)
ImGui_ImplOpenGL3_Init(glsl_version)

try
    clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]
    while !GLFW.WindowShouldClose(window)
        GLFW.PollEvents()
        # start the Dear ImGui frame
        ImGui_ImplOpenGL3_NewFrame()
        ImGui_ImplGlfw_NewFrame()
        CImGui.NewFrame()

        # show image example
        CImGui.Begin("GR Demo")

        @cstatic phi = Cfloat(0.0) begin
            @c CImGui.SliderFloat("Angle", &phi, 0, 360, "%.4f")
            beginprint(conid)
            draw(phi * π/180)
            endprint()
        end

        ImGui_ImplOpenGL3_UpdateImageTexture(image_id, image, img_width, img_height)
        CImGui.Image(Ptr{Cvoid}(image_id), (img_width, img_height))
        CImGui.End()

        # rendering
        CImGui.Render()
        GLFW.MakeContextCurrent(window)
        display_w, display_h = GLFW.GetFramebufferSize(window)
        glViewport(0, 0, display_w, display_h)
        glClearColor(clear_color...)
        glClear(GL_COLOR_BUFFER_BIT)
        ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())

        GLFW.MakeContextCurrent(window)
        GLFW.SwapBuffers(window)
    end
catch e
    @error "Error in renderloop!" exception=e
    Base.show_backtrace(stderr, catch_backtrace())
finally
    ImGui_ImplOpenGL3_Shutdown()
    ImGui_ImplGlfw_Shutdown()
    CImGui.DestroyContext(ctx)
    GLFW.HideWindow(window)
    GLFW.DestroyWindow(window)
end


text:
using Base.Test
using IndexedTables

let a = NDSparse([12,21,32], [52,41,34], [11,53,150]), b = NDSparse([12,23,32], [52,43,34], [56,13,10])
    p = collect(IndexedTables.product(a, b))
    @test p == [(11,56), (11,13), (11,10), (53,56), (53,13), (53,10), (150,56), (150,13), (150,10)]

    p = collect(IndexedTables.product(a, b, a))
    @test p == [(11,56,11),(11,56,53),(11,56,150),(11,13,11),(11,13,53),(11,13,150),(11,10,11),(11,10,53),(11,10,150),
                (53,56,11),(53,56,53),(53,56,150),(53,13,11),(53,13,53),(53,13,150),(53,10,11),(53,10,53),(53,10,150),
                (150,56,11),(150,56,53),(150,56,150),(150,13,11),(150,13,53),(150,13,150),(150,10,11),(150,10,53),(150,10,150)]
end

let a = [1:10;]
    @test IndexedTables._sizehint!(1:10, 20) == 1:10
    @test IndexedTables._sizehint!(a, 20) === a
end

@test Columns([1,2], [3,4]) == Columns([1,2], [3.0,4.0])
@test Columns([1,2], [3,4]) != Columns([1,2], [3.0,4.1])
@test Columns([1,2], [3,4]) != Columns(a=[1,2], b=[3,4])

function roundtrips(x)
    b = IOBuffer()
    serialize(b, x)
    seekstart(b)
    return deserialize(b) == x
end

@test roundtrips(Columns(rand(5), rand(5)))
@test roundtrips(Columns(c1 = rand(5), c2 = rand(5)))
@test roundtrips(convert(NDSparse, rand(3,3)))
@test roundtrips(NDSparse(Columns(y=rand(3), x=rand(3)), rand(3)))

let x = rand(3), y = rand(3), v = rand(3), w = rand(3)
    @test vcat(Columns(x,y), Columns(v,w)) == Columns(vcat(x,v), vcat(y,w))
    @test vcat(Columns(x=x,y=y), Columns(x=v,y=w)) == Columns(x=vcat(x,v), y=vcat(y,w))
end

import IndexedTables: _promote_op

let
    f = x -> 1/x
    @test _promote_op(f, Int) == Float64

    f = x -> (x,1/x)
    @test _promote_op(f, Int) == Tuple{Int, Float64}

    foo(i) = [1, 1//2, "x"][i]
    f = x -> (x,foo(x))
    @test _promote_op(f, Int) == Tuple{Int, Any}

    bar(i) = Union{Int, String}[1, "x"][i]
    g = x -> (x,bar(x))
    @test _promote_op(g, Int) == Tuple{Int, Any}

    h = x -> rand(Bool) ? @NT(x=1) : @NT(y=2)
    @test _promote_op(h, Int) == NamedTuple

    #101
    s = (1,)
    @test _promote_op(i -> Tuple(getfield(i, j) for j in s), Tuple{Int}) == Any

    # 97
    x = ndsparse(@NT(t=[0.01, 0.05]), @NT(x=[1,2], y=[3,4]))
    @test map(p->@NT(r = sum(p)), x).data == Columns([4,6], names=[:r])
end


text:
using Distributed
using JLD2

include("../BioClim/src/required.jl")
include("expandlayers_vulpes_vulpes.jl")

taxa_occ[3]

##Presence maps
# create new empty matrix for predator presence, dimensions of predator newlayers[1]
pred_presence = zeros(Int64, size(newlayers[1]))

# fill pred_presence matrix with binary values (1 if present, 0 if not)
#length of columns
for i in 1:size(newlayers[1].grid,2)
    #length of rows
    for j in 1:size(newlayers[1].grid,1)
        if isnan((newlayers[1].grid)[j,i])
            pred_presence[j,i] = 0
        else
            pred_presence[j,i] = 1
        end
    end
end

# view binary presence-absence result
pred_presence

# create new empty matrix for prey1 presence, Armadillidium opacum
prey1_pres = zeros(Int64, 2031, 4223)

# fill prey1_pres matrix with binary values (1 if present, 0 if not)
for i in 1:size(newlayers[16].grid,2)
    for j in 1:size(newlayers[16].grid,1)
        if isnan((newlayers[16].grid)[j,i])
            prey1_pres[j,i] = 0
        else
            prey1_pres[j,i] = 1
        end
    end
end
prey1_pres

#add predator and prey1 matrices
sum1 = pred_presence .+ prey1_pres

# #add only where predator is present
# for k in 1:length(sum1)
#     if pred_presence[k] == 1
#         sum1 = pred_presence[k] .+ prey1_pres[k]
#     end
# end

#filter only occurence data (ie. values of 2)
for l in 1:length(sum1)
    if sum1[l] == 2
        sum1[l] = 2
    else
        sum1[l] = 0
    end
end

#replace zeros to NaN for clearer image when mapping
sum1 = replace(sum1, 0 => NaN)

# Convert to SDMLayer
sum1_lay = SDMLayer(sum1, newlayers[1].left, newlayers[1].right, newlayers[1].bottom, newlayers[1].top)

plot_sum1 = plotSDM(sum1_lay)
plot_sum1 = plot!(title="Vulpes vulpes and Cynomys ludovicianus")
savefig("../SDM/figures/vulpes_vulpes/sdm_map-vulpes_vulpes_p15.pdf")

plot_sum1

####################
#just predator
worldmap = clip(worldshape(50), pred)

sdm_plot = plot([0.0], lab = " ", msw = 0.0, ms = 0.0, size = (900, 450), frame=:box,
                title = "Vulpes vulpes")
xaxis!(sdm_plot, (pred.left, pred.right), "Longitude")
yaxis!(sdm_plot, (pred.bottom, pred.top), "Latitude")

for p in worldmap
    sh = Shape([pp.x for pp in p.points], [pp.y for pp in p.points])
    plot!(sdm_plot, sh, c = :lightgrey, lab = " ")
end

heatmap!(
    sdm_plot,
    longitudes(pred), latitudes(pred), pred.grid,
    aspectratio = 92.60/60.75, c = :BuPu,
    clim = (0.0, maximum(filter(!isnan, pred.grid)))
    )

for p in worldmap
    xy = map(x -> (x.x, x.y), p.points)
    plot!(sdm_plot, xy, c=:grey, lab = " ", leg=false)
end
sdm_plot

savefig("../SDM/figures/vulpes_vulpes/dm_map-Vulpes_vulpes.pdf")

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM


text:
let
    model0 = iJR.load_model("max_model")

    bins = 50
    Fd_rxns_map = iJR.load_rxns_map()
    offsetf = 1.1

    exglcidx = ChU.rxnindex(model0, iJR.EX_GLC_IDER)
    biomidx = ChU.rxnindex(model0, iJR.BIOMASS_IDER)
    exglcL, exglcU = ChU.bounds(model0, exglcidx)

    maxD = maximum(Fd.val(:D)) * offsetf
    max_cgD_X = -maximum(Fd.ciD_X(:GLC)) * offsetf

    Ds = range(0.01, maxD; length = bins)
    cgD_Xs = range(max_cgD_X, exglcU; length = bins)

    box_vols = zeros(bins, bins)

    Fd_ids = ["GLC", "PYR", "SUCC", "LAC", "FORM", "AC"]
    model_ids = [Fd_rxns_map[Fd_id] for Fd_id in Fd_ids]
    model_idxs = [ChU.rxnindex(model0, model_id) for model_id in model_ids]

    # feeding task
    Ch = Channel(nthreads()) do Ch_
        @showprogress for (Di, D) in enumerate(Ds)
            for (cgD_Xi, cgD_X) in enumerate(cgD_Xs)
                put!(Ch_, (Di, D, cgD_Xi, cgD_X))
            end
        end
    end

    # compute volume map
    @threads for _ in 1:nthreads()
        thid = threadid()
        thid == 1 && continue
        model = deepcopy(model0)
        for (Di, D, cgD_Xi, cgD_X) in Ch

            # Reduce Pol
            ChU.lb!(model, exglcidx, cgD_X)
            ChU.bounds!(model, biomidx, D, D)
            try
                L, U = ChLP.fva(model, model_idxs)
                vol = prod(abs.(L .- U))
                box_vols[Di, cgD_Xi] = max(0.0, log10(vol + 1e-50))
            catch
                box_vols[Di, cgD_Xi] = NaN
            end
        end
    end

    # vol map
    p = heatmap(Ds, cgD_Xs, box_vols';
        title = "Polytope Box volume", label = "",
        xlabel = "D", ylabel = "cgD_X"
    )

    # exp vals
    scatter!(p, Fd.val(:D), -Fd.ciD_X(:GLC);
        label = "exp", color = :white,
        m = 8
    )

    mysavefig(p, "pol_box_volume"; bins)
end