text:
using FewShotAnomalyDetection
using UCI
using Flux
using MLDataPattern
using IPMeasures
using IPMeasures: crit_mmd2_var, crit_mxy_over_mltpl, crit_mxy_over_sum

hiddenDim = 32
latentDim = 3
numLayers = 3
nonlinearity = "swish"
layerType = "Dense"
β = 1
num_pseudoinputs = 1

batchSize = 100
numBatches = 1000

# dataset = "breast-cancer-wisconsin" # name of the UCI dataset. You can find the names in the e.g. in the LODA paper http://agents.fel.cvut.cz/stegodata/pdfs/Pev15-Loda.pdf - last page
# # dataset = "statlog-vehicle"
#
# data, normal_labels, anomaly_labels = UCI.get_umap_data(dataset)
# # This method iterates over all subdataset in a multiclass dataset or it only runs once for a single class dataset
# subdatasets = UCI.create_multiclass(data, normal_labels, anomaly_labels)
# (subdata, class_label) = subdatasets[1]
# println(dataset*" "*class_label)
# _X_tr, _y_tr, _X_tst, _y_tst = UCI.split_data(subdata, 0.8) #train/test split
#
# # dirty data normalization
# _X_tr .-= minimum(_X_tr)
# _X_tr ./= maximum(_X_tr)
# _X_tst .-= minimum(_X_tst)
# _X_tst ./= maximum(_X_tst)
#
# # no need to have it this way but we are used to it from legacy code
# train = (_X_tr, _y_tr)
# test = (_X_tst, _y_tst)

M = 2
N = 1000
s = Float32.([30 10; 10 1])
X = s*randn(Float32,M,N) .+ 100

svae = SVAEvampmeans(size(X, 1), hiddenDim, latentDim, numLayers, nonlinearity, layerType, num_pseudoinputs)

criterion = crit_mxy_over_mltpl
criterion = crit_mxy_over_sum
criterion = crit_mmd2_var

opt = Flux.Optimise.RMSProp(1e-4)
for i in 1:15
# Find the best width of the mmd kernel
z = FewShotAnomalyDetection.samplez(svae, zparams(svae, X)...).data
zp = FewShotAnomalyDetection.sampleVamp(svae, size(z, 2)).data
γs = -10:0.05:2
cs = [criterion(IPMeasures.IMQKernel(10.0 ^ γ), z, zp, IPMeasures.pairwisecos) for γ in γs]
γ = 10 ^ γs[argmax(cs)]
println("We chose kernal size $γ")

learnRepresentation(data) = wloss(svae, data, β, (x, y) -> FewShotAnomalyDetection.mmd_imq(x, y, γ))
cb = Flux.throttle(() -> println("SVAE: $(learnRepresentation(X))"), 10)
# there is a hack with RandomBatches because so far I can't manage to get them to work without the tuple - I have to find a different sampling iterator
Flux.train!(learnRepresentation, Flux.params(svae), RandomBatches((X,), size = batchSize, count = numBatches), opt, cb = cb)
end

s = [700, 700]
#Plotting

using Plots
plotlyjs()
scatter(X[1, :], X[2, :], size = s)
mus = Flux.Tracker.data(zparams(svae, X)[1])
xgz = Flux.Tracker.data(svae.g(mus))
scatter!(xgz[1, :], xgz[2, :], size = s)

scatter3d(mus[1, :], mus[2, :], mus[3, :], size = s)
zs = Flux.Tracker.data(FewShotAnomalyDetection.zfromx(svae, X))
scatter3d(zs[1, :], zs[2, :], zs[3, :], size = s, title = "sampled Z")

(μ, k) = zparams(svae, svae.pseudo_inputs)
lkhs = FewShotAnomalyDetection.pz(svae, X)
scatter3d(mus[1, :], mus[2, :], mus[3, :], zcolor = vec(lkhs), size = s)
scatter3d!(μ[1, :].data, μ[2, :].data, μ[3, :].data, color = "green", m = :xcross, markersize = 5, title = "Z means")

xs = minimum(X[1, :]):1:maximum(X[1, :])
ys = minimum(X[2, :]):1:maximum(X[2, :])

score = (x, y) -> FewShotAnomalyDetection.as_jacobian(svae, [x, y])
fillc = true
nlevels = 20
csvae = contour(xs, ys, score, fill = fillc, levels = nlevels, title = "Jacobian score")
scatter!(csvae, X[1, :], X[2, :], color = "green", alpha = 0.5)

score = (x, y) -> FewShotAnomalyDetection.pxexpectedz(svae, [x, y])[1]
csvae = contour(xs, ys, score, fill = fillc, levels = nlevels, title = "P(X) Vita")
scatter!(csvae, X[1, :], X[2, :], alpha = 0.5)

score = (x, y) -> log(FewShotAnomalyDetection.pz(svae, [x, y])[1])
csvae = contour(xs, ys, score, fill = fillc, levels = nlevels, title = "P(z)")
scatter!(csvae, X[1, :], X[2, :], alpha = 0.5)

score = (x, y) -> log(FewShotAnomalyDetection.pz(svae, [x, y])[1]) + FewShotAnomalyDetection.pxexpectedz(svae, [x, y])[1]
csvae = contour(xs, ys, score, fill = fillc, levels = nlevels, title = "P(X) + P(z)")
scatter!(csvae, X[1, :], X[2, :], alpha = 0.5)

# # one learning step
# for i in 1:1000
# @show l = learnRepresentation(X)
# Flux.Tracker.back!(l)
# for p in params(svae)
#    if any(isnan.(p.grad))
#       println("$p has gradient NaN")
#    end
#    Δ = Flux.Optimise.apply!(opt, p.data, p.grad)
#    p.data .-= Δ
#    p.grad .= 0
# end
# end

max_stars_repo_name: honzabim/FewShotAnomalyDetection

max_issues_repo_name: honzabim/FewShotAnomalyDetection

max_forks_repo_name: honzabim/FewShotAnomalyDetection


text:
typealias L1orL2HingeLoss Union{L1HingeLoss, L2HingeLoss}

# ==========================================================================
# Predicton type to dispatch on univariate/multivariate prediction

abstract PredictionType{INTERCEPT}

immutable UvPredicton{INTERCEPT} <: PredictionType{INTERCEPT}
  bias::Float64
  function UvPredicton(bias::Real)
    ((typeof(INTERCEPT) <: Bool) && (bias != 0.) == INTERCEPT) || throw(MethodError())
    new(Float64(bias))
  end
end
UvPredicton(bias::Real) = UvPredicton{bias!=0.}(bias)

immutable MvPredicton{INTERCEPT} <: PredictionType{INTERCEPT}
  bias::Float64
  function MvPredicton(bias::Real)
    ((typeof(INTERCEPT) <: Bool) && (bias != 0.) == INTERCEPT) || throw(MethodError())
    new(Float64(bias))
  end
end
MvPredicton(bias::Real) = MvPredicton{bias!=0.}(bias)

function prediction_type{T<:Real}(loss::MarginBasedLoss, y::AbstractVector{T}, bias::Real)
  for i in y
    if i < 1
      return UvPredicton(bias)
    elseif i > 1
      return MvPredicton(bias)
    end
  end
  return UvPredicton(bias)
end

function prediction_type{T<:Real}(loss::DistanceBasedLoss, y::AbstractVector{T}, bias::Float64)
  UvPredicton(bias)
end

function prediction_type{T<:Real}(loss::DistanceBasedLoss, y::AbstractMatrix{T}, bias::Float64)
  MvPredicton(bias)
end

# ==========================================================================
# Special solver for SVM problems

abstract SvmSolver <: AbstractSolver

# ==========================================================================
# Typed problem definition for returning the primal or dual solution

abstract SvmProblem
immutable PrimalProblem <: SvmProblem end
immutable DualProblem <: SvmProblem end

# ==========================================================================
# Internal exception that is thrown if an algorithm doesn't handle the
# multivariate prediction case itself. This will result in a generic handleing
# such as a one-vs-all in the case of classification

type MvNotNativelyHandled <: Exception end

# ==========================================================================
# These functions make sure that if no callback is defined, then the empty
# callback function will be optimized away through inlining

# Primal problem only
function _docallback(::Any, ::Void, iter, w, w⃗, f, G) false end
function _docallback(::PredictionType{false}, cb::Function, iter, w, w⃗, f, G) cb(iter, w⃗, f, G) == :Exit end
function _docallback(::PredictionType{true}, cb::Function, iter, w, w⃗, f, G) cb(iter, w, f, G) == :Exit end

# Dual problem only
function _docallback(::Any, ::Void, iter, α, f, G) false end
function _docallback(::PredictionType, cb::Function, iter, α, f, G) cb(iter, α, f, G) == :Exit end

# Primal or dual problem
function _docallback(::Any, ::Bool, ::Void, iter, w, w⃗, α, f, G) false end
function _docallback(::PredictionType{false}, dual::Bool, cb::Function, iter, w, w⃗, α, f, G) dual ? cb(iter, α, f, G) == :Exit : cb(iter, w⃗, f, G) == :Exit end
function _docallback(::PredictionType{true}, dual::Bool, cb::Function, iter, w, w⃗, α, f, G) dual ? cb(iter, α, f, G) == :Exit : cb(iter, w, f, G) == :Exit end

# ==========================================================================
# Custom print functions for the training output

function print_iter_head()
  @printf("%5s   %12s   %12s\n",
          "Iter", "f.value", "g.value")
  println("======================================================")
end

function print_iter(t::Int, J::Real, G::Real)
  @printf("%5d   %12.4e   %12.4e\n", t, J, G)
end

ext: jl


text:
include("src/DiscretisedFluidQueues.jl")
using .DiscretisedFluidQueues
using LinearAlgebra, SparseArrays, StableRNGs
using Test

T = [-2.5 2 0.5; 1 -2 1; 1 2 -3];

C = [0.0; 2.0; -3.0];
m = -1 .+ 2*Int.(DiscretisedFluidQueues._strictly_pos.(C));
S = DiscretisedFluidQueues.PhaseSet(C);

model = DiscretisedFluidQueues.BoundedFluidQueue(T,S);
am = DiscretisedFluidQueues.augment_model(model);
nodes = 0.0:4.0:12.0;
nodes = [0.0;3.0;4.0;12.0]
nbases = 3
dgmesh = DiscretisedFluidQueues.DGMesh(nodes,nbases)

# fv_order = 3
# fvmesh = DiscretisedFluidQueues.FVMesh(nodes,fv_order)

# order = 3;
# frapmesh = DiscretisedFluidQueues.FRAPMesh(nodes,order);

dq = DiscretisedFluidQueues.DiscretisedFluidQueue(am,dgmesh);
B = DiscretisedFluidQueues.build_lazy_generator(dq)
DiscretisedFluidQueues.fast_mul(B,Matrix{Float64}(I(40)))-B
DiscretisedFluidQueues.fast_mul(Matrix{Float64}(I(40)),B)-B

max_stars_repo_name: angus-lewis/StochasticFluidQueues

max_issues_repo_name: angus-lewis/StochasticFluidQueues

max_forks_repo_name: angus-lewis/StochasticFluidQueues


text:
module AbstractPaulis

using ..AbstractOps
import ..AbstractOps: _show_op_plain
import Random
import SparseArraysN
const SparseArrays = SparseArraysN
using StaticArrays
import LinearAlgebra
import IsApprox
import ..Utils: pow_of_minus_one
import ..AbstractOps: op_symbols

export AbstractPauli, pauli_vector

# """
#     AbstractPauli{T}
# Represents a single-qubit Pauli operator, one
# of `I`, `X`, `Y`, `Z`.
# """
abstract type AbstractPauli{T} <: AbstractOp end

const (Iop, Xop, Yop, Zop) = (0, 1, 2, 3)

@inline is_pauli_y(p::AbstractPauli) = op_index(p) === Yop
@inline is_pauli_I(p::AbstractPauli) = op_index(p) === Iop

####
#### Constructors
####

const _pauli_chars = ('I', 'X', 'Y', 'Z')
const _pauli_string_chars = string.(_pauli_chars)
const _pauli_symbol_chars = Symbol.(_pauli_chars)

op_symbols(::Type{<:AbstractPauli}, ::Type{<:AbstractChar}) = _pauli_chars
op_symbols(::Type{<:AbstractPauli}, ::Type{<:AbstractString}) = _pauli_string_chars
op_symbols(::Type{<:AbstractPauli}, ::Type{Symbol}) = _pauli_symbol_chars

# TODO: Test and remove this. The generic method is the same
# """
#     Vector{<:AbstractPauli}(ps::AbstractString)

# Construct a `Vector{<:AbstractPauli}` by parsing `ps` containing
# characters I, X, Y, Z.
#"""
#Vector{T}(ps::AbstractString) where {T <: AbstractPauli} = [T(s) for s in ps]

"""
    pauli_vector(::Type{PauliT}, op_index, n_qubits, indices=Vector{Int}(undef, n_qubits))

Return a `Vector` of `PauliT` representing a multi-qubit Pauli indexed by the base-4
representation of `op_index`. A temporary array `indices` may be passed to avoid
allocation.
"""
function pauli_vector(::Type{PauliT}, op_index::Integer, n_qubits::Integer,
                     indices=Vector{Int}(undef, n_qubits)) where PauliT
    digits!(indices, op_index; base=4)
    reverse!(indices)  # agree with sort order
    return PauliT.(indices)
end

AbstractOps.rand_ind_range(::Type{<:AbstractPauli}) = 0:3

####
#### IO
####

op_chars(::Type{AbstractPauli}) = _pauli_chars

## I think all of this was due to a bug elsewhere. Perhaps it can be reverted.
## Bizarre, some Int prints as 132, is converted to float -123.0
## and is not > 4, nor < 1
## So, we enumerate the values of valid indices.
function Base.show(io::IO, p::AbstractPauli)
    print(io, typeof(p), ": ")
    _show_op_plain(io, p)
end

function _show_op_plain(io::IO, p::AbstractPauli)
    i::Int = op_index(p) + 1
    if i != 1 && i != 2 && i != 3 && i != 4
        char = '0'
    else
        char = _pauli_chars[i]
    end
    print(io, char)
end

## Use this if `AbstractPauli <: AbstractMatrix`
## We need both definitions.
Base.display(io::IO, p::AbstractPauli) = show(io, p)
Base.display(p::AbstractPauli) = show(stdout, p)

## We may want to define `display` as well.
function Base.show(io::IO, ps::AbstractArray{<:AbstractPauli})
    for p in ps
        _show_op_plain(io, p)
    end
end

####
#### Conversion
####

## Use static matrices when constructing large matrices in loops.
## But, I don't seem to get any performance gain.
const _Imat = float.([1 0; 0 1])
const _Xmat = float.([0 1; 1 0])
const _Ymat = float.([0 -im; im 0])
const _Zmat = float.([1 0; 0 -1])

const _Isparse = SparseArrays.sparse(_Imat)
const _Xsparse = SparseArrays.sparse(_Xmat)
const _Ysparse = SparseArrays.sparse(_Ymat)
const _Zsparse = SparseArrays.sparse(_Zmat)

const _SImat = @SMatrix [1.0 0; 0 1]
const _SXmat = @SMatrix [0.0 1; 1 0]
const _SYmat = @SMatrix [0.0 -im; im 0]
const _SZmat = @SMatrix [1.0 0; 0 -1]

function Base.Matrix(p::AbstractPauli)
    matrices = (_Imat, _Xmat, _Ymat, _Zmat)
    return matrices[op_index(p) + 1]
end

function SparseArrays.sparse(p::AbstractPauli)
    matrices = (_Isparse, _Xsparse, _Ysparse, _Zsparse)
    return matrices[op_index(p) + 1]
end

####
#### Container interface
####

Base.size(::AbstractPauli) = (2, 2)

function smatrix(p::AbstractPauli)
    static = (_SImat, _SXmat, _SYmat, _SZmat)
    return static[op_index(p) + 1]
end

Base.getindex(p::AbstractPauli, i, j) = smatrix(p)[i, j]

####
#### Compare / predicates
####

Base.one(::Type{PauliT}) where PauliT <: AbstractPauli = PauliT(:I)

Base.one(::PauliT) where PauliT <: AbstractPauli = one(PauliT)

Base.:(==)(p1::AbstractPauli, p2::AbstractPauli) = op_index(p1) == op_index(p2)

"""
    isless(p1::AbstractPauli, p2::AbstractPauli)

Canonical (lexical) order of `AbstractPauli` is I < X < Y < Z.
"""
Base.isless(p1::AbstractPauli, p2::AbstractPauli) = op_index(p1) < op_index(p2)

IsApprox.ishermitian(::AbstractPauli) = true

IsApprox.isposdef(p::PauliT) where PauliT <: AbstractPauli = p === PauliT(:I)

IsApprox.isdiag(p::PauliT) where PauliT <: AbstractPauli = (p === PauliT(:I) || p === PauliT(:Z))

IsApprox.issymmetric(p::PauliT) where PauliT <: AbstractPauli = p != PauliT(:Y)

IsApprox.isunitary(::AbstractPauli) = true

####
#### Algebra / mathematical operations
####

Base.:*(p::AbstractPauli, m::AbstractMatrix) = smatrix(p) * m

Base.:*(m::AbstractMatrix, p::AbstractPauli) = m * smatrix(p)

Base.:^(p::AbstractPauli, n::Integer) = iseven(n) ? one(p) : p

Base.inv(p::AbstractPauli) = p

Base.adjoint(p::AbstractPauli) = p

IsApprox.commutes(p1::AbstractPauli, p2::AbstractPauli) = isone(p1) || isone(p2) || p1 === p2

IsApprox.anticommutes(p1::AbstractPauli, p2::AbstractPauli) = ! IsApprox.commutes(p1, p2)

"""
    phase(p1::AbstractPauli, p2::AbstractPauli)

Return a `NamedTuple` of two `Bool`s representing the phase of the product of `p1` and
`p2`. The elements are `:has_sign_flip` and `:has_imag_unit`. The phase can be reconstructed
as `im^has_imag_unit * (-1)^has_sign_flip`.
"""
function phase(p1::AbstractPauli, p2::AbstractPauli)
    if isone(p1) || isone(p2)
        has_imag_unit = false
        has_sign_flip = false
    else
        d = op_index(p1) - op_index(p2)
        if d == 0
            has_imag_unit = false
            has_sign_flip = false
        elseif d == 1 || d == -2  # Depends on stadard ordering of Paulis
            has_imag_unit = true
            has_sign_flip = true
        else
            has_imag_unit = true
            has_sign_flip = false
        end
    end
    return PauliPhase(has_sign_flip, has_imag_unit)
end

### Linear algebra

LinearAlgebra.tr(p::AbstractPauli) = is_pauli_I(p) ? 2 : 0

## This would otherwise be computed by LinearAlgebra calling getindex
LinearAlgebra.eigvals(p::AbstractPauli) = isone(p) ? [1.0, 1.0] : [-1.0, 1.0]

Base.length(::AbstractPauli) = 1
Base.iterate(p::AbstractPauli) = (p, nothing)
Base.iterate(::AbstractPauli, ::Any) = nothing
Base.isempty(::AbstractPauli) = false

Base.iszero(::AbstractPauli) = false

abstract type AbstractPhaseData end

struct PauliPhase
    sign_flip::Bool
    imag_unit::Bool
end

mutable struct PauliPhaseCounts <: AbstractPhaseData
    n_sign_flips::Int
    n_imag_units::Int
end

PauliPhaseCounts() = PauliPhaseCounts(0, 0)

AbstractOps.phase_data(::Type{<:AbstractPauli}) = PauliPhaseCounts()

function AbstractOps.accumulate_phase!(p_counts::PauliPhaseCounts, op1::T, op2::T) where {T <: AbstractPauli}
    pauli_phase = phase(op1, op2)
    return AbstractOps.accumulate_phase!(p_counts, pauli_phase)
end

function AbstractOps.accumulate_phase!(p_counts::PauliPhaseCounts, p_phase::PauliPhase)
    p_counts.n_sign_flips += p_phase.sign_flip
    p_counts.n_imag_units += p_phase.imag_unit
    return p_counts
end

function AbstractOps.accumulate_phase(old_phase_data, op1::T, op2::T) where {T <: AbstractPauli}
    phase_info = phase(op1, op2)
    new_phase_data = (old_phase_data[1] + phase_info.sign_flip, old_phase_data[2] + phase_info.imag_unit)
    return new_phase_data
end

AbstractOps.compute_phase(phase_counts::PauliPhaseCounts, _1, _2) =
    AbstractOps.compute_phase(AbstractPauli, (phase_counts.n_sign_flips, phase_counts.n_imag_units))

function AbstractOps.compute_phase(::Type{<:AbstractPauli}, phase_data)
    (n_sign_flips, n_imag_units) = phase_data
    _sign = pow_of_minus_one(n_sign_flips)
    n = n_imag_units % 4
    if n == 0
        cim = complex(1)
    elseif n == 1
        cim = complex(im)
    elseif n == 2
        cim = complex(-1)
    else
        cim = complex(-im)
    end
    return _sign * cim
end

####
#### Interface for subtypes
####

## Subtypes of AbstractPauli must implement these.

## The function `*` must return only the bare Pauli. That is, the phase +-im must be discarded.

# TODO: Is this MethodError worth anything? Won't a missing method throw a MethodError in any case?
"""
    p1::AbstractPauli * p2::AbstractPauli

Multiply single-qubit operators ignoring phase.
"""
function Base.:*(p1::AbstractPauli, p2::AbstractPauli)
    throw(MethodError(*, (p1, p2)))
end

"""
    abstract type AbstractPauli

Represents single-qubit Pauli matrices and their algebra.

    AbstractPauli(ind::Integer)::AbstractPauli

This function must be called as a concrete subtype of `AbstractPauli`.

Create an instance of (a subtype of) `AbstractPauli` given
the Pauli index in `[0, 3]`.

    AbstractPauli(s::Union{Symbol, AbstractString, AbstractChar})::AbstractPauli

Create an instance of (a subtype of) `AbstractPauli` from the `Symbol`s `:I, :X, :Y, :Z`,
or the `AbstractString`s `"I", "X", "Y", "Z"` or the `AbstractChar`s `'I', 'X', 'Y', 'Z'`.
In general creating `AbstractPauli`s with `String`s is slower than with the other types.
"""
function AbstractPauli end

end # module AbstractPaulis


text:
function MA(k,S)
    return k*S
end

function MAcom2(k,S1,S2)
    return k*S1*S2
end

function MM(k,Km,S)
    return k*S/(Km+S)
end

function Hill2(A,n)
    return A^n
end

function libPulseDelay(t,X,term)
    # signal parameter
    sinput = 1.0 - (1.0 - 0.01)*X
    sbase = 0.01
    slate = 0.05 - (0.05 - 0.01)*X
    tpulse = 0.25
    traise = 0.5
    tdecay = 0.5
    tdelay = 0.0

    if t <= tdelay || term == 1.0
        return sbase
    elseif tdelay < t <= (traise + tdelay)
        return (t - tdelay)*(sinput - sbase)/traise + sbase
    elseif (traise + tdelay) < t <= (tpulse + traise + tdelay)
        return sinput
    elseif (tpulse + traise + tdelay) < t
        return (sinput - slate)*exp(-(t - tpulse - traise - tdelay)/tdecay) + slate
    else
        return 0.0
    end
end

function diffeq!(du,u,h,p,t)

    #IKK activity
    IKKa = u[V.IKKp] + u[V.IKKpC] + u[V.IKKppC] + u[V.IKKpp]
    if IKKa == 0.0
        IKKa = 1.0
    end

    ##### time delay
    NFKBnDelay = h(p,t-p[C.delayrnae])[V.NFKBn]

    ###################----CBM module----#######################
    signal = libPulseDelay(t,p[C.X],p[C.term])

    v = Dict{Int64,Float64}()

    v[1] = MA(p[C.kCp0], u[V.C])                                                 # (*basal C phosphorylation*)
    v[2] = signal * MM(p[C.kCpS], p[C.kmCpS], u[V.C])                            # (*signal-dependent C phosphorylation*)
    v[3] = MM(p[C.kCpu], p[C.kmCpu], u[V.Cp])                                    # (* Cp dephosphorylation*)
    v[4] = MA(p[C.kCpB0], u[V.CB])                                               # (*basal CB phosphorylation*)
    v[5] = signal * MM(p[C.kCpBS], p[C.kmCpBS], u[V.CB])                         # (*signal-dependent CB phosphorylation*)
    v[6] = MM(p[C.kCpBu], p[C.kmCpBu], u[V.CpB])                                 # (* CpB dephosphorylation*)
    v[7] = MA(p[C.kCpM0], u[V.CM])                                               # (*basal CM phosphorylation*)
    v[8] = signal * MM(p[C.kCpMS], p[C.kmCpMS], u[V.CM])                         # (*signal-dependent CM phosphorylation*)
    v[9] = MM(p[C.kCpMu], p[C.kmCpMu], u[V.CpM])                                 # (* CpM dephosphorylation*)
    v[10] = MA(p[C.kCpBM0], u[V.CBM])                                            # (*basal CBM phosphorylation*)
    v[11] = signal * MM(p[C.kCpBMS], p[C.kmCpBMS], u[V.CBM])                     # (*signal-dependent CBM phosphorylation*)
    v[12] = MM(p[C.kCpBMu], p[C.kmCpBMu], u[V.CpBM])                             # (* CpBM dephosphorylation*)
    v[13] = MAcom2(p[C.kBaM], u[V.B], u[V.M])                                    # (*association between B and M*)
    v[14] = MA(p[C.kBdM], u[V.BM])                                               # (*disassociation between B and M*)
    v[15] = MAcom2(p[C.kCaB], u[V.C], u[V.B])                                    # (*association between C and B*)
    v[16] = MA(p[C.kCdB], u[V.CB])                                               # (*disassociation between C and B*)
    v[17] = MAcom2(p[C.kCpaB], u[V.Cp], u[V.B])                                  # (*association between Cp and B*)
    v[18] = MA(p[C.kCpdB], u[V.CpB])                                             # (*disassociation between Cp and B*)
    v[19] = MAcom2(p[C.kCaM], u[V.C], u[V.M])                                    # (*association between C and M*)
    v[20] = MA(p[C.kCdM], u[V.CM])                                               # (*disassociation between C and M*)
    v[21] = MAcom2(p[C.kCpaM], u[V.Cp], u[V.M])                                  # (*association between Cp and M*)
    v[22] = MA(p[C.kCpdM], u[V.CpM])                                             # (*disassociation between Cp and M*)
    v[23] = MAcom2(p[C.kCBaM], u[V.CB], u[V.M])                                  # (*association between CB and M*)
    v[24] = MA(p[C.kCBdM], u[V.CBM])                                             # (*disassociation between CB and M*)
    v[25] = MAcom2(p[C.kCpBaM], u[V.CpB], u[V.M])                                # (*association between CpB and M*)
    v[26] = MA(p[C.kCpBdM], u[V.CpBM])                                           # (*disassociation between CpB and M*)
    v[27] = MAcom2(p[C.kCMaB], u[V.CM], u[V.B])                                  # (*association between CM and B*)
    v[28] = MA(p[C.kCMdB], u[V.CBM])                                             # (*disassociation between CM and B*)
    v[29] = MAcom2(p[C.kCpMaB], u[V.CpM], u[V.B])                                # (*association between CpM and B*)
    v[30] = MA(p[C.kCpMdB], u[V.CpBM])                                           # (*disassociation between CpM and B*)
    v[31] = MAcom2(p[C.kCaBM], u[V.C], u[V.BM])                                  # (*association between C and BM*)
    v[32] = MA(p[C.kCdBM], u[V.CBM])                                             # (*disassociation between C and BM*)
    v[33] = MAcom2(p[C.kCpaBM], u[V.Cp], u[V.BM])                                # (*association between Cp and BM*)
    v[34] = MA(p[C.kCpdBM], u[V.CpBM])                                           # (*disassociation between Cp and BM*)
    ###################----CBM module----#######################

    ###################----TAK1 module----######################
    v[35] = MA(p[C.kTp0], u[V.TAK1])                                             # (*basal TAK1 phosphorylation*)
    v[36] = u[V.Cp] * MM(p[C.kCTpS], p[C.kmCTpS], u[V.TAK1])                     # (*Cp-dependent TAK1 phosphorylation*)
    v[37] = u[V.CpB] * MM(p[C.kCBTpS], p[C.kmCBTpS], u[V.TAK1])                  # (*CpB-dependent TAK1 phosphorylation*)
    v[38] = u[V.CpM] * MM(p[C.kCMTpS], p[C.kmCMTpS], u[V.TAK1])                  # (*CpM-dependent TAK1 phosphorylation*)
    v[39] = u[V.CpBM] * MM(p[C.kCBMTpS], p[C.kmCBMTpS], u[V.TAK1])               # (*CpBM-dependent TAK1 phosphorylation*)
    v[40] = u[V.IKKp] * MM(p[C.kTpIKK1], p[C.kmTpIKK1], u[V.TAK1])               # (*IKKp-dependent TAK1 phosphorylation*)
    v[41] = MM(p[C.kTpu], p[C.kmTpu], u[V.TAK1p])                                # (*TAK1 dephosphorylation*)
    v[42] = MA(p[C.kTpC0], u[V.TAK1C])                                           # (*basal TAK1C phosphorylation*)
    v[43] = u[V.Cp] * MM(p[C.kCTpCS], p[C.kmCTpCS], u[V.TAK1C])                  # (*Cp-dependent TAK1C phosphorylation*)
    v[44] = u[V.CpB] * MM(p[C.kCBTpCS], p[C.kmCBTpCS], u[V.TAK1C])               # (*CpB-dependent TAK1C phosphorylation*)
    v[45] = u[V.CpM] * MM(p[C.kCMTpCS], p[C.kmCMTpCS], u[V.TAK1C])               # (*CpM-dependent TAK1C phosphorylation*)
    v[46] = u[V.CpBM] * MM(p[C.kCBMTpCS], p[C.kmCBMTpCS], u[V.TAK1C])            # (*CpBM-dependent TAK1C phosphorylation*)
    v[47] = u[V.IKKp] * MM(p[C.kTpCIKK1], p[C.kmTpCIKK1], u[V.TAK1C])            # (*IKKp-dependent TAK1C phosphorylation*)
    v[48] = u[V.IKKpC] * MM(p[C.kTpCIKK2], p[C.kmTpCIKK2], u[V.TAK1C])           # (*IKKpC-dependent TAK1C phosphorylation*)
    v[49] = u[V.IKKppC] * MM(p[C.kTpCIKK3], p[C.kmTpCIKK3], u[V.TAK1C])          # (*IKKppC-dependent TAK1C phosphorylation*)
    v[50] = MM(p[C.kTpCu], p[C.kmTpCu], u[V.TAK1pC])                             # (*basal TAK1C dephosphorylation*)
    v[51] = MAcom2(p[C.kTaB], (u[V.CB] + u[V.CpB]), u[V.TAK1])                   # (*association between TAK1 and CB*)
    v[52] = MAcom2(p[C.kTaM], (u[V.CM] + u[V.CpM]), u[V.TAK1])                   # (*association between TAK1 and CM*)
    v[53] = MAcom2(p[C.kTaBM], (u[V.CBM] + u[V.CpBM]), u[V.TAK1])                # (*association between TAK1 and CBM*)
    v[54] = MA(p[C.kTCd], u[V.TAK1C])                                            # (*dissociation between TAK1 and CBM*)
    v[55] = MAcom2(p[C.kTpaB], (u[V.CB] + u[V.CpB]), u[V.TAK1p])                 # (*association between TAK1p and CBM*)
    v[56] = MAcom2(p[C.kTpaM], (u[V.CM] + u[V.CpM]), u[V.TAK1p])                 # (*association between TAK1p and CBM*)
    v[57] = MAcom2(p[C.kTpaBM], (u[V.CBM] + u[V.CpBM]), u[V.TAK1p])              # (*association between TAK1p and CBM*)
    v[58] = MA(p[C.kTpCd], u[V.TAK1pC])                                          # (*dissociation between TAK1 and CBM*)
    ###################----TAK1 module----######################

    ###################----IKK module----#######################
    v[59] = MA(p[C.kIp0], u[V.IKK])                                              # (*basal IKK phosphorylation*)
    v[60] = u[V.TAK1p] * MM(p[C.kIpTAKp], p[C.kmIpTAKp], u[V.IKK])               # (*TAK1p-dependent IKK phosphorylation*)
    v[61] = u[V.TAK1pC] * MM(p[C.kIpTAKpC], p[C.kmIpTAKpC], u[V.IKK])            # (*TAKpC-dependent IKK phosphorylation*)
    v[62] = MA(p[C.kICp0], u[V.IKKC])                                            # (*basal IKKC phosphorylation*)
    v[63] = u[V.TAK1p] * MM(p[C.kICpTAKp], p[C.kmICpTAKp], u[V.IKKC])            # (*TAK1p-dependent IKKC phosphorylation*)
    v[64] = u[V.TAK1pC] * MM(p[C.kICpTAKpC], p[C.kmICpTAKpC], u[V.IKKC])         # (*TAK1pC-dependent IKKC phosphorylation*)
    v[65] = MM(p[C.kIpu], p[C.kmIpu], u[V.IKKp])                                 # (*IKKp dephosphorylation*)
    v[66] = MM(p[C.kIpCu], p[C.kmIpCu], u[V.IKKpC])                              # (*IKKpC dephosphorylation*)
    v[67] = MAcom2(p[C.kIaB], (u[V.CB] + u[V.CpB]), u[V.IKK])                    # (*association between IKK and CB*)
    v[68] = MAcom2(p[C.kIaM], (u[V.CM] + u[V.CpM]), u[V.IKK])                    # (*association between IKK and CM*)
    v[69] = MAcom2(p[C.kIaBM], (u[V.CBM] + u[V.CpBM]), u[V.IKK])                 # (*association between IKK and CBM*)
    v[70] = MA(p[C.kICd], u[V.IKKC])                                             # (*dissociation between IKK and CBM*)
    v[71] = MAcom2(p[C.kIpaB], (u[V.CB] + u[V.CpB]), u[V.IKKp])                  # (*association between IKKp and CB*)
    v[72] = MAcom2(p[C.kIpaM], (u[V.CM] + u[V.CpM]), u[V.IKKp])                  # (*association between IKKp and CM*)
    v[73] = MAcom2(p[C.kIpaBM], (u[V.CBM] + u[V.CpBM]), u[V.IKKp])               # (*association between IKKp and CBM*)
    v[74] = MA(p[C.kIpCd], u[V.IKKpC])                                           # (*dissociation between IKKp and CBM*)
    v[75] = MM(p[C.kIpCfaIKKpC], p[C.kmIpCfaIKKpC], u[V.IKKpC])                  # (*basal IKKpC phosphorylation*)
    v[76] = u[V.IKKppC] * MM(p[C.kIpCfaIKKppC], p[C.kmIpCfaIKKppC], u[V.IKKpC])  # (*IKKppC-dependent IKKpC phosphorylation*)
    v[77] = MM(p[C.kIppCu], p[C.kmIppCu], u[V.IKKppC])                           # (*basal IKKppC dephosphorylation*)
    v[78] = MA(p[C.kIppCd], u[V.IKKppC])                                         # (*dissociation between IKKpp and CBM (IKKppC->IKKpp*)
    v[79] = MM(p[C.kIpphf], p[C.kmIpphf], u[V.IKKpp])                            # (*basal IKKpp inactivation (IKKpp->IKKi)*)
    v[80] = MM(p[C.kIppChf], p[C.kmIppChf], u[V.IKKppC])                         # (*basal IKKppC inactivation (IKKppC->IKKi)*)
    v[81] = MM(p[C.kIir], p[C.kmIir], u[V.IKKi])                                 # (* recycling (IKKi->IKK)*)
    v[82] = u[V.A20c] * MM(p[C.kIpA20], p[C.kmIpA20], u[V.IKKp])                 # (* A20-dependent IKKp phosphorylation inhibition *)
    v[83] = u[V.A20c] * MM(p[C.kIpCA20], p[C.kmIpCA20], u[V.IKKpC])              # (* A20-dependent IKKpC phosphorylation inhibition *)
    ###################----IKK module----#######################

    ########----IKKNFkBIkB module----################
    v[84] = MAcom2(p[C.kassanfkbikk], IKKa, u[V.NFKBIKBac])                      # (*association between IKKa and NFkBIkBac*)
    v[85] = MAcom2(p[C.kassbnfkbikk], IKKa, u[V.NFKBIKBbc])                      # (*association between IKKa and NFkBIkBbc*)
    v[86] = MAcom2(p[C.kassenfkbikk], IKKa, u[V.NFKBIKBec])                      # (*association between IKKa and NFkBIkBec*)
    v[87] = MA(p[C.kdisanfkbikk], u[V.IKKNFKBIKBac])                             # (*dissociation between IKKa and NFkBIkBac*)
    v[88] = MA(p[C.kdisbnfkbikk], u[V.IKKNFKBIKBbc])                             # (*dissociation between IKKa and NFkBIkBbc*)
    v[89] = MA(p[C.kdisenfkbikk], u[V.IKKNFKBIKBec])                             # (*dissociation between IKKa and NFkBIkBec*)
    v[90] = MA(p[C.kdegboundaIKK], u[V.IKKNFKBIKBac])                            # (*dissociation between IKKa and NFkBc and degradation of IkBac*)
    v[91] = MA(p[C.kdegboundbIKK], u[V.IKKNFKBIKBbc])                            # (*dissociation between IKKa and NFkBc and degradation of IkBbc*)
    v[92] = MA(p[C.kdegboundeIKK], u[V.IKKNFKBIKBec])                            # (*dissociation between IKKa and NFkBc and degradation of IkBec*)
    v[93] = MAcom2(p[C.kassaikk], u[V.IKBac], IKKa)                              # (*association between IKKa and IkBac*)
    v[94] = MAcom2(p[C.kassbikk], u[V.IKBbc], IKKa)                              # (*association between IKKa and IkBbc*)
    v[95] = MAcom2(p[C.kasseikk], u[V.IKBec], IKKa)                              # (*association between IKKa and IkBec*)
    v[96] = MA(p[C.kdisaikk], u[V.IKKIKBac])                                     # (*dissociation between IKKa and IkBac*)
    v[97] = MA(p[C.kdisbikk], u[V.IKKIKBbc])                                     # (*dissociation between IKKa and IkBbc*)
    v[98] = MA(p[C.kdiseikk], u[V.IKKIKBec])                                     # (*dissociation between IKKa and IkBec*)
    v[99] = MA(p[C.kdegfreeaIKK], u[V.IKKIKBac])                                 # (*dissociation of IKKa and degradation of IkBac*)
    v[100] = MA(p[C.kdegfreebIKK], u[V.IKKIKBbc])                                # (*dissociation of IKKa and degradation of IkBbc*)
    v[101] = MA(p[C.kdegfreeeIKK], u[V.IKKIKBec])                                # (*dissociation of IKKa and degradation of IkBec*)
    ########----IKKNFkBIkB module----################

    ###################----NFkB module----######################
    v[102] = MAcom2(p[C.kassa], u[V.IKBac], u[V.NFKBc])                          # (*association between NFkBc and IkBac*)
    v[103] = MA(p[C.kdisa], u[V.NFKBIKBac])                                      # (*dissociation between NFkBc and IkBac*)
    v[104] = MA(p[C.kdegbounda], u[V.NFKBIKBac])                                 # (*dissociation of NFkBc and degradation of IkBac*)
    v[105] = MAcom2(p[C.kassaikknfkb], u[V.IKKIKBac], u[V.NFKBc])                # (*association between NFkBc and IKKaIkBac*)
    v[106] = MA(p[C.kdisaikknfkb], u[V.IKKNFKBIKBac])                            # (*dissociation between NFkBc and IKKaIkBac*)
    v[107] = MAcom2(p[C.kassa], u[V.IKBan], u[V.NFKBn])                          # (*association between NFkBn and IkBan*)
    v[108] = MA(p[C.kdisa], u[V.NFKBIKBan])                                      # (*dissociation between NFkBn and IkBan*)
    v[109] = MA(p[C.kdegbounda], u[V.NFKBIKBan])                                 # (*dissociation of NFkBn and degradation of IkBan*)

    v[110] = MAcom2(p[C.kassb], u[V.IKBbc], u[V.NFKBc])                          # (*association between NFkBc and IkBbc*)
    v[111] = MA(p[C.kdisb], u[V.NFKBIKBbc])                                      # (*dissociation between NFkBc and IkBbc*)
    v[112] = MA(p[C.kdegboundb], u[V.NFKBIKBbc])                                 # (*dissociation of NFkBc and degradation of IkBbc*)
    v[113] = MAcom2(p[C.kassbikknfkb], u[V.IKKIKBbc], u[V.NFKBc])                # (*association between NFkBc and IKKaIkBbc*)
    v[114] = MA(p[C.kdisbikknfkb], u[V.IKKNFKBIKBbc])                            # (*dissociation between NFkBc and IKKaIkBbc*)
    v[115] = MAcom2(p[C.kassb], u[V.IKBbn], u[V.NFKBn])                          # (*association between NFkBn and IkBbn*)
    v[116] = MA(p[C.kdisb], u[V.NFKBIKBbn])                                      # (*dissociation between NFkBn and IkBbn*)
    v[117] = MA(p[C.kdegboundb], u[V.NFKBIKBbn])                                 # (*dissociation of NFkBn and degradation of IkBbn*)

    v[118] = MAcom2(p[C.kasse], u[V.IKBec], u[V.NFKBc])                          # (*association between NFkBc and IkBec*)
    v[119] = MA(p[C.kdise], u[V.NFKBIKBec])                                      # (*dissociation between NFkBc and IkBec*)
    v[120] = MA(p[C.kdegbounde], u[V.NFKBIKBec])                                 # (*dissociation of NFkBc and degradation of IkBec*)
    v[121] = MAcom2(p[C.kasseikknfkb], u[V.IKKIKBec], u[V.NFKBc])                # (*association between NFkBc and IKKaIkBec*)
    v[122] = MA(p[C.kdiseikknfkb], u[V.IKKNFKBIKBec])                            # (*dissociation between NFkBc and IKKaIkBec*)
    v[123] = MAcom2(p[C.kasse], u[V.IKBen], u[V.NFKBn])                          # (*association between NFkBn and IkBen*)
    v[124] = MA(p[C.kdise], u[V.NFKBIKBen])                                      # (*dissociation between NFkBn and IkBen*)
    v[125] = MA(p[C.kdegbounde], u[V.NFKBIKBen])                                 # (*dissociation of NFkBn and degradation of IkBen*)

    v[126] = MA(p[C.kshutboundikbain], u[V.NFKBIKBac])                           # (*transport NFkBIkBac into nucleos*)
    v[127] = MA(p[C.kshutboundikbaout], u[V.NFKBIKBan])                          # (*transport NFkBIkBan into cytoplasm*)
    v[128] = MA(p[C.kshutboundikbbin], u[V.NFKBIKBbc])                           # (*transport NFkBIkBbc into nucleos*)
    v[129] = MA(p[C.kshutboundikbbout], u[V.NFKBIKBbn])                          # (*transport NFkBIkBbn into cytoplasm*)
    v[130] = MA(p[C.kshutboundikbein], u[V.NFKBIKBec])                           # (*transport NFkBIkBec into nucleos*)
    v[131] = MA(p[C.kshutboundikbeout], u[V.NFKBIKBen])                          # (*transport NFkBIkBen into cytoplasm*)
    v[132] = MA(p[C.kshutfreeikbain], u[V.IKBac])                                # (*transport IkBac into nucleos*)
    v[133] = MA(p[C.kshutfreeikbaout], u[V.IKBan])                               # (*transport IkBan into cytoplasm*)
    v[134] = MA(p[C.kshutfreeikbbin], u[V.IKBbc])                                # (*transport IkBbc into nucleos*)
    v[135] = MA(p[C.kshutfreeikbbout], u[V.IKBbn])                               # (*transport IkBbn into cytoplasm*)
    v[136] = MA(p[C.kshutfreeikbein], u[V.IKBec])                                # (*transport IkBec into nucleos*)
    v[137] = MA(p[C.kshutfreeikbeout], u[V.IKBen])                               # (*transport IkBen into cytoplasm*)
    v[138] = MA(p[C.kshutnfkbin], u[V.NFKBc])                                    # (*transport NFkBc into nucleos*)
    v[139] = MA(p[C.kshutnfkbout], u[V.NFKBn])                                   # (*transport NFkBn into cytoplasm*)

    v[140] = p[C.k0mrnaikba]                                                     # (*basal mRNA(IkBa) transcription*)
    v[141] = p[C.kprodmrnaikba] * Hill2(u[V.NFKBn], p[C.khillprodmrnaikba])      # (*NFkB-induced IkBa transcription*)
    v[142] = MA(p[C.kdegmrnaikba], u[V.mRNAac])                                  # (*mRNA(IkBa) degradation*)
    v[143] = MA(p[C.kpikba], u[V.mRNAac])                                        # (*mRNA(IkBa) translation*)
    v[144] = MA(p[C.kdegfreea], u[V.IKBac])                                      # (*IkBac degradation*)
    v[145] = MA(p[C.kdegfreea], u[V.IKBan])                                      # (*IkBan degradation*)

    v[146] = p[C.k0mrnaikbb]                                                     # (*basal mRNA(IkBb) transcription*)
    v[147] = MA(p[C.kdegmrnaikbb], u[V.mRNAbc])                                  # (*mRNA(IkBb) degradation*)
    v[148] = MA(p[C.kpikbb], u[V.mRNAbc])                                        # (*mRNA(IkBb) translation*)
    v[149] = MA(p[C.kdegfreeb], u[V.IKBbc])                                      # (*IkBbc degradation*)
    v[150] = MA(p[C.kdegfreeb], u[V.IKBbn])                                      # (*IkBbn degradation*)

    v[151] = p[C.k0mrnaikbe]                                                     # (*basal mRNA(IkBe) transcription*)
    v[152] = p[C.kprodmrnaikbe] * Hill2(NFKBnDelay, p[C.khillprodmrnaikbe])      # (*NFkB-induced IkBe transcription*)
    v[153] = MA(p[C.kdegmrnaikbe], u[V.mRNAec])                                  # (*mRNA(IkBe) degradation*)
    v[154] = MA(p[C.kpikbe], u[V.mRNAec])                                        # (*mRNA(IkBe) translation*)
    v[155] = MA(p[C.kdegfreee], u[V.IKBec])                                      # (*IkBec degradation*)
    v[156] = MA(p[C.kdegfreee], u[V.IKBen])                                      # (*IkBen degradation*)

    v[157] = p[C.k0mrnaa20]                                                      # (*basal mRNA(A20) transcription*)
    v[158] = p[C.kprodmrnaa20] * Hill2(u[V.NFKBn], p[C.khillprodmrnaa20])        # (*NFkB-induced A20 transcription*)
    v[159] = MA(p[C.kdegmrnaa20], u[V.mRNAa20c])                                 # (*mRNA(A20) degradation*)
    v[160] = MA(p[C.kpa20], u[V.mRNAa20c])                                       # (*mRNA(A20) translation*)
    v[161] = MA(p[C.kdega20], u[V.A20c])                                         # (*A20 degradation*)

    IKKaNFKBIKB = - v[93] - v[94] - v[95] + v[99] + v[96] + v[100] + v[97] + v[101] + v[98] + v[90] + v[87] +
                    v[91] + v[88] + v[92] + v[89] - v[84] - v[85] - v[86]

    ###################----CBM module----######################
    du[V.B] = (- v[15] + v[16]) + (- v[17] + v[18]) + (- v[27] + v[28]) + (- v[29] + v[30]) + (- v[13] + v[14])
    du[V.M] = (- v[19] + v[20]) + (- v[21] + v[22]) + (- v[23] + v[24]) + (- v[25] + v[26]) + (- v[13] + v[14])
    du[V.BM] = ( v[13] - v[14]) + (- v[31] + v[32]) + (- v[33] + v[34])
    du[V.C] = (- v[15] + v[16]) + (- v[19] + v[20]) + (- v[31] + v[32]) + ( v[3] -(v[1] + v[2]) )
    du[V.CB] = (v[15] - v[16]) + (- v[23] + v[24]) + ( v[6] -(v[4] + v[5]) )
    du[V.CM] = (v[19] - v[20]) + (- v[27] + v[28]) + ( v[9] -(v[7] + v[8]) )
    du[V.Cp] = (- v[17] + v[18]) + (- v[21] + v[22]) + (- v[33] + v[34]) + (- v[3] +(v[1] + v[2]) )
    du[V.CpB] = (v[17] - v[18]) + (- v[25] + v[26]) + (- v[6] +(v[4] + v[5]) )
    du[V.CpM] = (v[21] - v[22]) + (- v[29] + v[30]) + (- v[9] +(v[7] + v[8]) )
    du[V.CBM] = (v[27] - v[28]) + (v[23] - v[24]) + ( v[31] - v[32]) + ( v[12] -(v[10] + v[11]) )
    du[V.CpBM] = (v[29] - v[30]) + (v[25] - v[26]) + ( v[33] - v[34]) + ( - v[12] +(v[10] + v[11]) )
    ###################----CBM module----######################

    ###################----TAK1 module---######################
    du[V.TAK1] =  -(v[35] + v[36] + v[37] + v[38] + v[39] + v[40])  + v[41] - (v[51] + v[52] + v[53]) + v[54]
    du[V.TAK1p] =  (v[35] + v[36] + v[37] + v[38] + v[39] + v[40]) - v[41] - (v[55] + v[56] + v[57]) + v[58]
    du[V.TAK1C] = -(v[42] + v[43] + v[44] + v[45] + v[46] + v[47] + v[48] + v[49] ) + v[50] + (v[51] + v[52] + v[53]) - v[54]
    du[V.TAK1pC] =  (v[42] + v[43] + v[44] + v[45] + v[46] + v[47] + v[48] + v[49] ) - v[50] + (v[55] + v[56] + v[57]) - v[58]
    ###################----TAK1 module---######################

    ###################----IKK module----######################
    du[V.IKK] =  -(v[59] + v[60] + v[61]) + v[65] + v[81] - (v[67] + v[68] + v[69]) + v[70] + v[82]
    du[V.IKKC] = -(v[62] + v[63] + v[64]) + v[66] + (v[67] + v[68] + v[69]) - v[70] + v[83]
    du[V.IKKp] = (v[59] + v[60] + v[61]) - v[65] - (v[71] + v[72] + v[73]) + v[74] - v[82] + u[V.IKKp]*IKKaNFKBIKB/IKKa
    du[V.IKKpC] = (v[71] + v[72] + v[73]) - v[74] - ( v[75] + v[76]) + v[77]+ (v[62] + v[63] + v[64]) - v[66] - v[83] + u[V.IKKpC]*IKKaNFKBIKB/IKKa
    du[V.IKKppC] =  ( v[75] + v[76]) - v[77] - v[78] - v[80] + u[V.IKKppC]*IKKaNFKBIKB/IKKa
    du[V.IKKpp] = v[78] - v[79] + u[V.IKKpp]*IKKaNFKBIKB/IKKa
    du[V.IKKi] = v[79] + v[80] - v[81]
    ###################----IKK module----######################

    ###################----NFkB module----######################
    du[V.NFKBc] = v[90] + v[91] + v[92] - v[102] + v[103]+ v[104] - v[105] + v[106] - v[110] + v[111] +
                v[112] - v[113] + v[114] - v[118] + v[119] + v[120] - v[121] + v[122] - v[138] + v[139]
    du[V.NFKBn] = - v[107] + v[108] + v[109] - v[115] + v[116] + v[117] - v[123] + v[124] + v[125] + v[138] - v[139]
    du[V.IKBac] = - v[93] + v[96] - v[102] + v[103] - v[132] + v[133] + v[143] - v[144]
    du[V.IKBbc] = - v[94] + v[97] - v[110] + v[111] - v[134] + v[135] + v[148] - v[149]
    du[V.IKBec] = - v[95] + v[98] - v[118] + v[119] - v[136] + v[137] + v[154] - v[155]
    du[V.IKBan] = - v[107] + v[108] + v[132] - v[133] - v[145]
    du[V.IKBbn] = - v[115] + v[116] + v[134] - v[135] - v[150]
    du[V.IKBen] = - v[123] + v[124] + v[136] - v[137] - v[156]

    du[V.NFKBIKBac] = - v[84] + v[87] + v[102] - v[103] - v[104] - v[126] + v[127]
    du[V.NFKBIKBbc] = - v[85] + v[88] + v[110] - v[111] - v[112] - v[128] + v[129]
    du[V.NFKBIKBec] = - v[86] + v[89] + v[118] - v[119] - v[120] - v[130] + v[131]
    du[V.NFKBIKBan] = v[107] - v[108] - v[109] + v[126] - v[127]
    du[V.NFKBIKBbn] = v[115] - v[116] - v[117] + v[128] - v[129]
    du[V.NFKBIKBen] = v[123] - v[124] - v[125] + v[130] - v[131]
    du[V.IKKIKBac] = v[93] - v[96] - v[99] - v[105] + v[106]
    du[V.IKKIKBbc] = v[94] - v[97] - v[100] - v[113] + v[114]
    du[V.IKKIKBec] = v[95] - v[98] - v[101] - v[121] + v[122]
    du[V.IKKNFKBIKBac] =  v[84] - v[87] - v[90] + v[105] - v[106]
    du[V.IKKNFKBIKBbc] =  v[85] - v[88] - v[91] + v[113] - v[114]
    du[V.IKKNFKBIKBec] =  v[86] - v[89] - v[92] + v[121] - v[122]

    du[V.mRNAac] = v[140] + v[141] - v[142]
    du[V.mRNAbc] = v[146] - v[147]
    du[V.mRNAec] = v[151] + v[152] - v[153]
    du[V.mRNAa20c] = v[157] + v[158] - v[159]
    du[V.A20c] = v[160] - v[161]
    ###################----NFkB module----######################
end

function param_values()

    p = zeros(C.NUM)

    p[C.kCp0] = 1.000000e-04
    p[C.kCpB0] = 1.000000e-04
    p[C.kCpM0] = 1.000000e-04
    p[C.kCpBM0] = 1.000000e-04
    p[C.kCpS] = 7.000000e+00
    p[C.kCpBS] = 2.000000e+01
    p[C.kCpMS] = 1.700000e+01
    p[C.kCpBMS] = 1.000000e+00
    p[C.kmCpS] = 1.000000e+00
    p[C.kmCpBS] = 1.000000e+00
    p[C.kmCpMS] = 1.000000e+00
    p[C.kmCpBMS] = 5.000000e+00
    p[C.kCpu] = 7.500000e-01
    p[C.kCpBu] = 1.950000e+00
    p[C.kCpMu] = 3.750000e+00
    p[C.kCpBMu] = 2.500000e-01
    p[C.kmCpu] = 1.000000e+00
    p[C.kmCpBu] = 1.000000e+00
    p[C.kmCpMu] = 1.000000e+00
    p[C.kmCpBMu] = 2.500000e+00

    p[C.kBaM] = 5.000000e-01
    p[C.kCaB] = 5.000000e-01
    p[C.kCaM] = 5.000000e-01
    p[C.kCBaM] = 1.500000e+00
    p[C.kCMaB] = 5.000000e-01
    p[C.kCaBM] = 5.000000e-01
    p[C.kCpaBM] = 1.885000e+00
    p[C.kCpaB] = 7.885000e+00
    p[C.kCpaM] = 4.585000e+00
    p[C.kCpBaM] = 2.500000e+00
    p[C.kCpMaB] = 1.500000e+00
    p[C.kBdM] = 2.500000e+00
    p[C.kCdB] = 1.000000e+00
    p[C.kCdM] = 2.500000e+00
    p[C.kCBdM] = 1.500000e+00
    p[C.kCMdB] = 2.500000e+00
    p[C.kCpdB] = 3.000000e-01
    p[C.kCpdM] = 3.000000e-01
    p[C.kCpBdM] = 1.500000e-01
    p[C.kCpMdB] = 4.500000e-01
    p[C.kCdBM] = 2.000000e-01
    p[C.kCpdBM] = 2.500000e-01

    p[C.kTp0] = 2.700000e-01
    p[C.kCTpS] = 3.200000e-01
    p[C.kCBTpS] = 1.320000e+00
    p[C.kCMTpS] = 1.320000e+00
    p[C.kCBMTpS] = 1.320000e+00
    p[C.kmCTpS] = 2.100000e+00
    p[C.kmCBTpS] = 2.100000e+00
    p[C.kmCMTpS] = 2.100000e+00
    p[C.kmCBMTpS] = 2.100000e+00
    p[C.kTpIKK1] = 1.670000e+01
    p[C.kmTpIKK1] = 7.800000e-01
    p[C.kTpu] = 3.400000e+00
    p[C.kmTpu] = 1.790000e+00
    p[C.kTpC0] = 3.300000e-01
    p[C.kCTpCS] = 5.500000e+00
    p[C.kCBTpCS] = 1.750000e+01
    p[C.kCMTpCS] = 1.650000e+01
    p[C.kCBMTpCS] = 2.450000e+01
    p[C.kmCTpCS] = 5.200000e-01
    p[C.kmCBTpCS] = 5.200000e-01
    p[C.kmCMTpCS] = 5.200000e-01
    p[C.kmCBMTpCS] = 5.200000e-01
    p[C.kTpCIKK1] = 1.000000e+01
    p[C.kmTpCIKK1] = 7.800000e-01
    p[C.kTpCIKK2] = 1.710000e+01
    p[C.kmTpCIKK2] = 1.300000e-01
    p[C.kTpCIKK3] = 8.000000e+01
    p[C.kmTpCIKK3] = 7.000000e-01
    p[C.kTpCu] = 1.533000e+01
    p[C.kmTpCu] = 8.000000e-01
    p[C.kTaB] = 1.790000e+02
    p[C.kTaM] = 1.690000e+02
    p[C.kTaBM] = 3.690000e+02
    p[C.kTpaB] = 3.790000e+00
    p[C.kTpaM] = 3.790000e+00
    p[C.kTpaBM] = 3.790000e+00
    p[C.kTCd] = 4.000000e+02
    p[C.kTpCd] = 3.210000e+00

    p[C.kIp0] = 5.500000e-02
    p[C.kIpTAKp] = 1.930000e+00
    p[C.kmIpTAKp] = 4.000000e-01
    p[C.kIpTAKpC] = 3.800000e+00
    p[C.kmIpTAKpC] = 5.800000e-01
    p[C.kIpu] = 1.500000e+01
    p[C.kmIpu] = 6.700000e-01
    p[C.kICp0] = 5.100000e-02
    p[C.kICpTAKp] = 8.640000e+00
    p[C.kmICpTAKp] = 7.100000e+00
    p[C.kICpTAKpC] = 1.230000e+01
    p[C.kmICpTAKpC] = 7.177000e+00
    p[C.kIpCu] = 1.250000e+01
    p[C.kmIpCu] = 5.410000e+00
    p[C.kIaB] = 3.030000e+01
    p[C.kIaM] = 2.230000e+01
    p[C.kIaBM] = 5.030000e+01
    p[C.kIpaB] = 3.130000e+01
    p[C.kIpaM] = 2.330000e+01
    p[C.kIpaBM] = 6.030000e+01
    p[C.kICd] = 3.280000e+01
    p[C.kIpCd] = 1.320000e+00
    p[C.kIpCfaIKKpC] = 1.000000e-05
    p[C.kmIpCfaIKKpC] = 2.560000e+00
    p[C.kIpCfaIKKppC] = 1.085700e+02
    p[C.kmIpCfaIKKppC] = 2.010000e+00
    p[C.kIppCu] = 1.530000e+00
    p[C.kmIppCu] = 2.000000e-01
    p[C.kIppCd] = 4.500000e-01
    p[C.kIppChf] = 1.000000e-01
    p[C.kmIppChf] = 5.000000e-01
    p[C.kIpphf] = 9.000000e-03
    p[C.kmIpphf] = 1.520000e+00
    p[C.kIir] = 1.440000e+00
    p[C.kmIir] = 3.440000e+00
    p[C.kIpA20] = 2.000000e+01
    p[C.kmIpA20] = 1.000000e-03
    p[C.kIpCA20] = 2.000000e+01
    p[C.kmIpCA20] = 1.000000e-03

    p[C.kassaikk] = 1.631051e-03
    p[C.kassaikknfkb] = 1.305073e+01
    p[C.kassanfkbikk] = 4.812297e-02
    p[C.kdisaikk] = 1.375438e-01
    p[C.kdisaikknfkb] = 2.296074e-05
    p[C.kdisanfkbikk] = 3.030804e-02
    p[C.kassa] = 3.655315e+01
    p[C.kdisa] = 3.495116e-04
    p[C.kdegbounda] = 2.089360e-03
    p[C.kdegboundaIKK] = 1.774837e-01
    p[C.kdegfreea] = 2.010515e-02
    p[C.kdegfreeaIKK] = 2.265038e-03
    p[C.kshutboundikbain] = 4.102701e-01
    p[C.kshutboundikbaout] = 1.632918e-01
    p[C.kshutfreeikbain] = 1.011242e-01
    p[C.kshutfreeikbaout] = 1.158262e-03
    p[C.k0mrnaikba] = 3.450000e-05
    p[C.kdegmrnaikba] = 8.326217e-02
    p[C.khillprodmrnaikba] = 1.620503e+00
    p[C.kpikba] = 8.484488e+00
    p[C.kprodmrnaikba] = 1.290500e-01

    p[C.kassbikk] = 5.609029e-03
    p[C.kassbikknfkb] = 1.423728e+02
    p[C.kassbnfkbikk] = 1.541680e-02
    p[C.kdisbikk] = 2.771705e-02
    p[C.kdisbikknfkb] = 1.742842e-03
    p[C.kdisbnfkbikk] = 1.068148e-01
    p[C.kassb] = 2.437209e+01
    p[C.kdegboundb] = 2.804468e-04
    p[C.kdegboundbIKK] = 3.645044e-01
    p[C.kdegfreeb] = 1.486087e-01
    p[C.kdegfreebIKK] = 2.818781e-04
    p[C.kdisb] = 1.856372e-04
    p[C.kshutboundikbbin] = 2.155431e-01
    p[C.kshutboundikbbout] = 7.783843e-01
    p[C.kshutfreeikbbin] = 1.505364e-03
    p[C.kshutfreeikbbout] = 8.350086e-04
    p[C.k0mrnaikbb] = 9.336099e-05
    p[C.kdegmrnaikbb] = 4.101465e-03
    p[C.khillprodmrnaikbb] = 0.000000e+00
    p[C.kpikbb] = 4.171048e-02
    p[C.kprodmrnaikbb] = 0.000000e+00
    p[C.kasseikk] = 4.423523e-03
    p[C.kasseikknfkb] = 1.485352e+02
    p[C.kassenfkbikk] = 1.531672e-02
    p[C.kdiseikk] = 3.827432e-01
    p[C.kdiseikknfkb] = 2.124084e-04
    p[C.kdisenfkbikk] = 2.691468e-01

    p[C.kasse] = 2.413478e+02
    p[C.kdegbounde] = 1.712910e-02
    p[C.kdegboundeIKK] = 6.934859e-02
    p[C.kdegfreee] = 2.201754e-01
    p[C.kdegfreeeIKK] = 4.819481e-05
    p[C.kdise] = 2.215707e-02
    p[C.kshutboundikbein] = 8.666062e-03
    p[C.kshutboundikbeout] = 4.006261e-02
    p[C.kshutfreeikbein] = 1.364442e-03
    p[C.kshutfreeikbeout] = 1.014252e-04
    p[C.k0mrnaikbe] = 1.232627e-04
    p[C.kdegmrnaikbe] = 9.773320e-03
    p[C.khillprodmrnaikbe] = 1.944545e+00
    p[C.kpikbe] = 2.384772e-01
    p[C.kprodmrnaikbe] = 1.052365e-01

    p[C.kshutnfkbin] = 2.072176e-01
    p[C.kshutnfkbout] = 5.537618e-04

    p[C.k0mrnaa20] = 1.118559e-03
    p[C.kdegmrnaa20] = 2.809754e-01
    p[C.kprodmrnaa20] = 1.788303e+00
    p[C.kpa20] = 2.445338e-01
    p[C.kdega20] = 6.217866e+00
    p[C.khillprodmrnaa20] = 1.243583e+00
    p[C.delayrnae] = 4.500000e+01

    p[C.X] = 0.0
    p[C.term] = 1.0

    return p
end

function initial_values()

    u0 = zeros(V.NUM)

    u0[V.B] = 1.0 # B
    u0[V.M] = 1.0 # M
    u0[V.C] = 1.0 # C
    u0[V.TAK1] = 1.0 # TAK1
    u0[V.IKK] = 1.0 # IKK
    u0[V.NFKBc] = 1.0 # NFKBc

    return u0
end


text:
function get_initial_population(
        n_population::Int64,
        n_gene::Int64)::Matrix{Float64}
    println(
        "Generating the initial population...\n"
    )
    flush(stdout)
    population::Matrix{Float64} = fill(
        Inf, (n_population, n_gene + 1)
    )
    for i = 1:n_population
        while !isfinite(population[i,end])
            for j = 1:n_gene
                population[i,j] = rand()
            end
            population[i,end] = objective(
                population[i,1:n_gene]
            )
        end
        println(
            @sprintf(
                "%d / %d", i, n_population
            )
        )
        flush(stdout)
    end
    println(
        "\n----------------------------------------\n"
    )
    population = sortslices(population, dims = 1, by = x->x[end])

    return population
end

function get_initial_population_continue(
        nth_param_set::Int64,
        n_population::Int64,
        n_gene::Int64,
        p0_bounds::Vector{Float64})::Matrix{Float64}
    generation::Int64 = readdlm(
        "./fitparam/$nth_param_set/generation.dat"
    )[1,1]
    best_indiv::Vector{Float64} = readdlm(
        @sprintf(
            "./fitparam/%d/fit_param%d.dat", nth_param_set, generation
        )
    )[:,1]
    println(
        "\n----------------------------------------\n"*
        "Generating the initial population...\n"
    )
    flush(stdout)
    population::Matrix{Float64} = fill(
        Inf, (n_population, n_gene + 1)
    )
    for i = 1:n_population
        while !isfinite(population[i,end])
            for j = 1:n_gene
                population[i,j] = encode_bestIndivVal2randGene(
                    j, best_indiv, p0_bounds
                )
                if population[i,j] > 1.0
                    population[i,j] = 1.0
                elseif population[i,j] < 0.0
                    population[i,j] = 0.0
                end
            end
            population[i,end] = objective(
                population[i,1:n_gene]
            )
        end
        println(
            @sprintf(
                "%d / %d", i, n_population
            )
        )
        flush(stdout)
    end
    println(
        "\n----------------------------------------\n"
    )
    population = sortslices(population, dims = 1, by = x->x[end])

    return population
end


text:
@testset "Durou" begin
    p, q = synthetic_gradient(SynthSphere(50))
    Z = Durou()(p, q)
    @test Z[38,38] ≈ 1.5742826661339246
    @test Z[1,1] ≈ 0.44970037692631165
    @test Z[76,76] ≈ 47.228212708897075
end


text:
function GridPolygon=PolygonToGrid(Polygon,Grid,OrientFace)
GridPolygon.nBar=[ 0  1   0   -1
                   1  0  -1    0];
GridPolygon.Dim=3;
GridPolygon.Type='Quad';
NumNodesFace=size(Polygon.N,2);
NumNodes=2*NumNodesFace+1;
Nodes(1:NumNodes)=Node([0 0 0],0);
NodeNumber=1;
BoundaryNodeNumber=1;
BoundaryNode=zeros(1,2*NumNodesFace);
%BoundaryNode=zeros(1,1);
MidPoint=Grid.Nodes(Polygon.N(end)).P;
for i=1:NumNodesFace-1
  Nodes(NodeNumber)=Node(Grid.Nodes(Polygon.N(i)).P',NodeNumber);
  MidPoint=MidPoint+Grid.Nodes(Polygon.N(i)).P;
  BoundaryNode(BoundaryNodeNumber)=NodeNumber;
  BoundaryNodeNumber=BoundaryNodeNumber+1;
  NodeNumber=NodeNumber+1;
  Nodes(NodeNumber)=Node(0.5*(Grid.Nodes(Polygon.N(i)).P'+Grid.Nodes(Polygon.N(i+1)).P'),NodeNumber);
  BoundaryNode(BoundaryNodeNumber)=NodeNumber;
  BoundaryNodeNumber=BoundaryNodeNumber+1;
  NodeNumber=NodeNumber+1;
end
Nodes(NodeNumber)=Node(Grid.Nodes(Polygon.N(NumNodesFace)).P',NodeNumber);
BoundaryNode(BoundaryNodeNumber)=NodeNumber;
BoundaryNodeNumber=BoundaryNodeNumber+1;
NodeNumber=NodeNumber+1;
Nodes(NodeNumber)=Node(0.5*(Grid.Nodes(Polygon.N(NumNodesFace)).P'+Grid.Nodes(Polygon.N(1)).P'),NodeNumber);
BoundaryNode(BoundaryNodeNumber)=NodeNumber;
NodeNumber=NodeNumber+1;
MidPoint=(1/NumNodesFace)*MidPoint;
Nodes(NodeNumber)=Node(MidPoint',NodeNumber);

GridPolygon.Nodes=Nodes;
GridPolygon.NumNodes=NumNodes;
GridPolygon.BoundaryNode=BoundaryNode;
InteriorNode(1)=NodeNumber;
GridPolygon.InteriorNode=InteriorNode;

NumEdges=3*NumNodesFace;
Edges(1:NumEdges)=Edge([1 2],GridPolygon,0,0);
EdgeNumber=1;
BoundaryEdgeNumber=1;
BoundaryEdge=zeros(1,2*NumNodesFace);
InteriorEdgeNumber=1;
InteriorEdge=zeros(1,NumNodesFace);
N1=1;
N2=2;
EdgeNumberI=1;
EdgeNumberB=1;
for i=1:NumNodesFace-1
  Edges(EdgeNumber)=Edge([N1 N2],GridPolygon,EdgeNumber,EdgeNumberB+NumNodesFace);
  BoundaryEdge(BoundaryEdgeNumber)=EdgeNumberB+NumNodesFace;
  BoundaryEdgeNumber=BoundaryEdgeNumber+1;
  EdgeNumber=EdgeNumber+1;
  EdgeNumberB=EdgeNumberB+1;
  N1=N1+1;
  N2=N2+1;
  Edges(EdgeNumber)=Edge([N1 N2],GridPolygon,EdgeNumber,EdgeNumberB+NumNodesFace);
  BoundaryEdge(BoundaryEdgeNumber)=EdgeNumberB+NumNodesFace;
  BoundaryEdgeNumber=BoundaryEdgeNumber+1;
  EdgeNumber=EdgeNumber+1;
  EdgeNumberB=EdgeNumberB+1;
  N1=N1+1;
  N2=N2+1;
end
Edges(EdgeNumber)=Edge([N1 N2],GridPolygon,EdgeNumber,EdgeNumberB+NumNodesFace);
BoundaryEdge(BoundaryEdgeNumber)=EdgeNumberB+NumNodesFace;
BoundaryEdgeNumber=BoundaryEdgeNumber+1;
EdgeNumber=EdgeNumber+1;
EdgeNumberB=EdgeNumberB+1;
N1=N1+1;
N2=1;
Edges(EdgeNumber)=Edge([N1 N2],GridPolygon,EdgeNumber,EdgeNumberB+NumNodesFace);
BoundaryEdge(BoundaryEdgeNumber)=EdgeNumberB+NumNodesFace;
EdgeNumber=EdgeNumber+1;
EdgeNumberB=EdgeNumberB+1;

N1=2;
N2=2*NumNodesFace+1;
for i=1:NumNodesFace
  Edges(EdgeNumber)=Edge([N1 N2],GridPolygon,EdgeNumber,EdgeNumberI);
  InteriorEdge(InteriorEdgeNumber)=EdgeNumberI;
  InteriorEdgeNumber=InteriorEdgeNumber+1;
  EdgeNumber=EdgeNumber+1;
  EdgeNumberI=EdgeNumberI+1;
  N1=N1+2;
end
GridPolygon.Edges=Edges;
GridPolygon.NumEdges=NumEdges;
GridPolygon.BoundaryEdge=BoundaryEdge;
GridPolygon.InteriorEdge=InteriorEdge;
GridPolygon.NumEdgesI=EdgeNumberI-1;
GridPolygon.NumEdgesB=EdgeNumberB-1;

NumFaces=NumNodesFace;
Faces(1:NumFaces)=Face([0 0 0 0],GridPolygon,0,'x');
E1=2*NumNodesFace;
E2=3*NumNodesFace;
E3=2*NumNodesFace+1;
E4=1;
FaceNumber=1;
Type='o';
[Faces(FaceNumber),GridPolygon]=Face([E1 E2 E3 E4],GridPolygon,FaceNumber,Type,OrientFace);
FaceNumber=FaceNumber+1;
E1=2;
E2=2*NumNodesFace+1;
E3=2*NumNodesFace+2;
E4=3;
for i=2:NumNodesFace
  [Faces(FaceNumber),GridPolygon]=Face([E1 E2 E3 E4],GridPolygon,FaceNumber,Type,OrientFace);
  FaceNumber=FaceNumber+1;
  E1=E1+2;
  E2=E2+1;
  E3=E3+1;
  E4=E4+2;
end
GridPolygon.Faces=Faces;
GridPolygon.NumFaces=NumFaces;
GridPolygon=Orientation(GridPolygon);
GridPolygon=Renumbering(GridPolygon);
end


text:
# * Split Hamiltonian

"""
    OrbitalHamiltonianTerm(i, j, coeff, A, integrals)

Represents a term in the orbital Hamiltonian arising from a variation
of the energy expressions between configurations `i` and `j` in the
multi-configurational expansion. `coeff` is the numeric coefficient,
`A` is the operator acting on the orbital, and `integrals` is a vector
of [`OrbitalIntegral`](@ref)s arising from the presence of
non-orthogonal orbitals and whose values should be multiplied to form
the overall coefficient.
"""
struct OrbitalHamiltonianTerm{aO,bO,T,QO,VOint}
    i::Int
    j::Int
    coeff::T
    A::QO
    integrals::VOint
end

function Base.show(io::IO, t::OrbitalHamiltonianTerm)
    write(io, "OrbitalHamiltonianTerm at [$(t.i),$(t.j)]: ")
    show(io, t.A)
    write(io, " × ($(t.coeff)) (#∫: $(length(t.integrals)))")
end

Base.show(io::IO, ::MIME"text/plain", t::OrbitalHamiltonianTerm) = show(io, t)

"""
    coefficient(term::OrbitalHamiltonianTerm)

Return the multiplicative coefficient pertaining to `term`,
_excluding_ the `conj(c_i)*c_j` mixing coefficients, due to the
configuration-interaction.
"""
function coefficient(term::OrbitalHamiltonianTerm)
    isempty(term.integrals) && return term.coeff
    term.coeff*prod(integral_value, term.integrals)
end

"""
    coefficient(term::OrbitalHamiltonianTerm, c::Vector)

Return the multiplicative coefficient pertaining to `term`,
_including_ the `conj(c_i)*c_j` mixing coefficients, due to the
configuration-interaction.
"""
function coefficient(term::OrbitalHamiltonianTerm, c::AbstractVector)
    coeff = coefficient(term)
    # We use zero-valued i,j coordinates to designate terms which do
    # not originate from the configuration expansion, but rather added
    # in manually, such as energy shifts.
    (iszero(term.i) || iszero(term.j)) && return coeff

    conj(c[term.i])*c[term.j]*coeff
end

update!(::AbstractMatrix) = nothing
update!(::AbstractMatrix, ::Atom) = nothing
update!(t::OrbitalHamiltonianTerm) = update!(t.A)
update!(t::OrbitalHamiltonianTerm, atom::Atom) = update!(t.A, atom)

"""
    OrbitalHamiltonian(R, terms, mix_coeffs, projector, orbital)

The Hamiltonian for `orbital` is constructed from a radial basis `R`,
a set of [`OrbitalHamiltonianTerm`](@ref) `terms` that describe the
various interactions between orbitals, `mix_coeffs` which are the
mixing coefficents for the multi-configurational expansion. The
`projector` ensures orthogonality between orbital pairs which have
Lagrange multipliers associated with them, by projecting out
components of other orbitals every time the `OrbitalHamiltonian`
action on `orbital` is computed.
"""
mutable struct OrbitalHamiltonian{aO,bO,O,T,Proj,RT}
    R::RT
    terms::Vector{OrbitalHamiltonianTerm{aO,bO,T}}
    mix_coeffs::Vector{T}
    projector::Proj
    orbital::O
end

Base.axes(hamiltonian::OrbitalHamiltonian) =
    (axes(hamiltonian.R,1),axes(hamiltonian.R,1))

Base.axes(hamiltonian::OrbitalHamiltonian, i) =
    axes(hamiltonian)[i]

Base.eltype(hamiltonian::OrbitalHamiltonian{aO,bO,O,T,Proj}) where {aO,bO,O,T,Proj} = T

update!(h::OrbitalHamiltonian) = foreach(t -> update!(t), h.terms)
function update!(h::OrbitalHamiltonian, atom::Atom)
    foreach(t -> update!(t, atom), h.terms)
    h.mix_coeffs = atom.mix_coeffs
end

"""
    energy_matrix!(H, hamiltonian, ϕ)

Compute the contribution of `hamiltonian` to the Hamiltonian matrix
`H` by repeatedly acting on the associated radial orbital `ϕ` with the
different multi-configurational [`OrbitalHamiltonianTerm`](@ref)s of
`hamiltonian`.
"""
function SCF.energy_matrix!(H::HM, hamiltonian::OrbitalHamiltonian{aO,bO,O,T},
                            ϕ::RadialOrbital{T}) where {HM<:AbstractMatrix,aO,bO,O,T}
    tmp = similar(ϕ)
    for term in hamiltonian.terms
        materialize!(MulAdd(coefficient(term), term.A, ϕ, zero(T), tmp))
        H[term.i,term.j] += (ϕ'tmp)[1]
    end
    H
end

"""
    filter(fun::Function, H::OrbitalHamiltonian)

Filter the [`OrbitalHamiltonianTerm`](@ref)s of `H` according to the
predicate `fun`.
"""
Base.filter(fun::Function, H::OrbitalHamiltonian) =
    OrbitalHamiltonian(H.R, filter(fun, H.terms),
                       H.mix_coeffs, H.projector, H.orbital)

function Base.getindex(H::OrbitalHamiltonian, term::Symbol)
    term == :total && return H
    operator_type = if term == :onebody
        AtomicOneBodyHamiltonian
    elseif term == :direct
        DirectPotential
    elseif term == :exchange
        ExchangePotential
    elseif term == :source
        SourceTerm
    elseif term == :shift
        ShiftTerm
    else
        throw(ArgumentError("Unknown Hamiltonian term $term"))
    end

    filter(t -> t.A isa operator_type, H)
end

# ** Materialization

const OrbitalHamiltonianMatrixElement{aO,bO,O,T} =
    Mul{<:Any,<:Tuple{<:AdjointRadialOrbital{T},
                      <:OrbitalHamiltonian{aO,bO,O,T},
                      <:RadialOrbital{T}}}

const OrbitalHamiltonianMatrixVectorProduct{aO,bO,O,T} =
    Mul{<:Any,<:Tuple{<:OrbitalHamiltonian{aO,bO,O,T},<:RadialOrbital{T}}}

const OrbitalHamiltonianMatrixMatrixProduct{aO,bO,O,T} =
    Mul{<:Any,<:Tuple{<:OrbitalHamiltonian{aO,bO,O,T},<:RadialOrbitals{T}}}

Base.eltype(::OrbitalHamiltonianMatrixVectorProduct{aO,bO,O,T}) where {aO,bO,O,T} = T

function Base.copyto!(dest::RadialOrbital{T},
                      matvec::OrbitalHamiltonianMatrixVectorProduct{aO,bO,O,T}) where {aO,bO,O,T}
    axes(dest) == axes(matvec) || throw(DimensionMismatch("axes must be the same"))
    R′,v = dest.args
    hamiltonian,b = matvec.args

    c = hamiltonian.mix_coeffs

    v .= zero(T)

    # Compute the additive action of each of the terms in the
    # Hamiltonian, weighting them by their coefficients which are
    # constructed from
    # 1) the configuration-interaction mixing coefficients,
    # 2) coefficients due to the energy expression, and
    # 3) values of various orbitals integrals resulting from
    #    non-orthogonal orbitals.
    for term in hamiltonian.terms
        coeff = coefficient(term, c)
        materialize!(MulAdd(coeff, term.A, b, one(T), dest))
    end

    dest
end

function Base.copyto!(dest::RadialOrbitals,
                      matvec::OrbitalHamiltonianMatrixMatrixProduct)
    axes(dest) == axes(matvec) || throw(DimensionMismatch("axes must be the same"))
    R′,dv = dest.args
    hamiltonian,b = matvec.args
    bv = b.args[2]

    n = size(dv,2)
    for j = 1:n
        copyto!(applied(*, hamiltonian.R, view(dv, :, j)),
                hamiltonian ⋆ applied(*, hamiltonian.R, view(bv, :, j)))
    end
    dest
end

Base.similar(matvec::OrbitalHamiltonianMatrixVectorProduct) =
    similar(matvec.args[2])

LazyArrays.materialize(matvec::OrbitalHamiltonianMatrixVectorProduct) =
    copyto!(similar(matvec), matvec)

function LazyArrays.materialize(matel::OrbitalHamiltonianMatrixElement)
    a,op,b = matel.args
    materialize(applied(*, a, materialize(op⋆b)))
end

# *** Materialization into a matrix

get_operator_matrix(A::AtomicOneBodyHamiltonian) = A.op.args[2]
get_operator_matrix(A::DirectPotential) = A.V̂.args[2]
get_operator_matrix(A::ShiftTerm) = A.shift
get_operator_matrix(A::RadialOperator) = A.args[2]
get_operator_matrix(A::Diagonal) = A
get_operator_matrix(A::IdentityOperator) = I
get_operator_matrix(A::SourceTerm) = get_operator_matrix(A.operator)
get_operator_matrix(A) =
    throw(ArgumentError("Don't know how to materialize a $(typeof(A)) as a matrix"))

"""
    copyto!(dest::AbstractMatix, hamiltonian::OrbitalHamiltonian)

Materialize the orbital `hamiltonian` into matrix form and store it in
`dest`, using the current values of all other orbitals. This is only
possible if the orbital `hamiltonian` does *not* contain any
[`ExchangePotential`](@ref)s or [`SourceTerm`](@ref)s (which are not
diagonal in orbital space), since the former is non-local (and thus
not representable as a matrix) and the latter is not a linear operator
(but an affine one).

Typical usage is to compute an easily factorizable matrix that can be
used for preconditioning the solution of the full equation.
"""
function Base.copyto!(dest::M, hamiltonian::OrbitalHamiltonian) where {T,M<:AbstractMatrix{T}}
    m = size(hamiltonian.R,2)
    size(dest) == (m,m) || throw(DimensionMismatch("axes must be the same"))

    c = hamiltonian.mix_coeffs

    dest .= zero(T)

    for term in hamiltonian.terms
        term.A isa SourceTerm && term.A.source_orbital ≠ hamiltonian.orbital &&
            throw(ArgumentError("It is not possible to materialize an orbitally off-diagonal $(typeof(term.A)) as a matrix"))
        term.A isa ExchangePotential &&
            throw(ArgumentError("It is not possible to materialize a $(typeof(term.A)) as a $M"))
        coeff = coefficient(term, c)

        dest += coeff*get_operator_matrix(term.A)
    end

    dest
end

function Base.similar(h::OrbitalHamiltonian{aO,bO,O,T,Proj,RT}, ::Type{T}) where {aO,bO,O,T,Proj,RT<:AbstractFiniteDifferences}
    R = h.R
    m = size(R,2)
    # TODO: This is only valid for RadialDifferences of
    # FiniteDifferencesQuasi.jl
    o = ones(T, m)
    SymTridiagonal(o,0*o[2:end])
end

Base.similar(h::OrbitalHamiltonian{aO,bO,O,T,Proj,RT}, ::Type{T}) where {aO,bO,O,T,Proj,RT<:BasisOrRestricted{<:FEDVR}} =
    Matrix(undef, h.R)

LazyArrays.materialize(h::OrbitalHamiltonian) =
    copyto!(similar(h, eltype(h)), h)

# ** Arithmetic

"""
    h::OrbitalHamiltonian + λ::UniformScaling

Shift the [`OrbitalHamiltonian`](@ref) `h` by `λ`.
"""
function Base.:(+)(h::OrbitalHamiltonian{aO,bO,O,T,Proj,RT}, λ::UniformScaling) where {aO,bO,O,T,Proj,RT}
    # The zeros designate that the shift is not to be weighted by the
    # mixing coefficients
    shift_term = OrbitalHamiltonianTerm(0, 0, one(T), ShiftTerm(λ),
                                        Vector{OrbitalIntegral{<:Any,aO,bO,T}}())
    OrbitalHamiltonian{aO,bO,O,T,Proj,RT}(h.R, vcat(h.terms, shift_term),
                                          h.mix_coeffs, h.projector, h.orbital)
end

"""
    h::OrbitalHamiltonian - λ::UniformScaling

Shift the [`OrbitalHamiltonian`](@ref) `h` by `-λ`.
"""
Base.:(-)(h::OrbitalHamiltonian, λ::UniformScaling) = h + (-λ)

# ** Krylov wrapper

"""
    SCF.KrylovWrapper(hamiltonian::OrbitalHamiltonian)

Construct a `KrylovWrapper` such that `hamiltonian`, that acts on
function spaces, can be used in a Krylov solver, which works with
linear algebra vector spaces.
"""
SCF.KrylovWrapper(hamiltonian::OrbitalHamiltonian{aO,bO,O,T,Proj}) where {aO,bO,O,T,Proj} =
    KrylovWrapper{T,OrbitalHamiltonian{aO,bO,O,T,Proj}}(hamiltonian)

Base.size(hamiltonian::OrbitalHamiltonian, ::SCF.KrylovWrapper) =
    (size(hamiltonian.R,2),size(hamiltonian.R,2))

"""
    mul!(y, A::KrylovWrapper{T,<:OrbitalHamiltonian}, x)

Materialize the action of the [`OrbitalHamiltonian`](@ref) on the
linear algebra vector `x` and store the result in `y`, by wrapping
them both with the `QuasiMatrix` necessary to transform `x` and `y` to
the function space of the Hamiltonian.
"""
LinearAlgebra.mul!(y::V₁, A::KrylovWrapper{T,Hamiltonian}, x::V₂) where {V₁,V₂,T,Hamiltonian<:OrbitalHamiltonian} =
    copyto!(applied(*, A.hamiltonian.R, y),
            A.hamiltonian⋆(applied(*, A.hamiltonian.R, x)))

function SCF.OrthogonalKrylovWrapper(hamiltonian::OrbitalHamiltonian)
    R = hamiltonian.R
    S = R'R
    SCF.OrthogonalKrylovWrapper(hamiltonian, KrylovWrapper(hamiltonian),
                                length(hamiltonian.projector.orbitals), S)
end

function update!(okw::SCF.OrthogonalKrylovWrapper{<:OrbitalHamiltonian})
    for (j,ϕ) in enumerate(okw.h.projector.ϕs)
        copyto!(view(okw.C, :, j), ϕ.args[2])
    end
    okw
end

# ** Preconditioner

"""
    IterativeFactorizations.preconditioner(hamiltonian::OrbitalHamiltonian)

Return a factorization of the matrix corresponding to `hamiltonian`,
where all terms arising from exchange and configuration interaction
have been removes, since they cannot be represented by a matrix.
"""
function IterativeFactorizations.preconditioner(hamiltonian::OrbitalHamiltonian)
    # To form the preconditioner, we select all terms of the shifted
    # Hamiltonian, except the exchange potentials and source terms,
    # which are not factorizable.
    Ph = filter(t -> !(t.A isa ExchangePotential || t.A isa SourceTerm), hamiltonian)
    hm = materialize(Ph)
    factorize(hm)
end

# ** Factorization
IterativeFactorizations.factorization(hamiltonian::OrbitalHamiltonian, args...; kwargs...) =
    IterativeFactorization(KrylovWrapper(hamiltonian), args...;
                           isposdefA=false, verbosity=0, kwargs...)

LinearAlgebra.factorize(hamiltonian::OrbitalHamiltonian) =
    factorization(hamiltonian)


text: Dict{String, Any}("tbl_tbl_empty" => Dict{String, Any}("tbl_0" => Dict{String, Any}()), "tbl_tbl_val" => Dict{String, Any}("tbl_1" => Dict{String, Any}("one" => Dict{String, Any}("value" => "1", "type" => "integer"))), "arr_arr_tbls" => Any[Any[Dict{String, Any}("one" => Dict{String, Any}("value" => "1", "type" => "integer")), Dict{String, Any}("two" => Dict{String, Any}("value" => "2", "type" => "integer"))]], "tbl_arr_tbl" => Dict{String, Any}("arr_tbl" => Any[Dict{String, Any}("one" => Dict{String, Any}("value" => "1", "type" => "integer"))]), "arr_arr_tbl_val" => Any[Any[Dict{String, Any}("one" => Dict{String, Any}("value" => "1", "type" => "integer"))]], "arr_tbl_tbl" => Any[Dict{String, Any}("tbl" => Dict{String, Any}("one" => Dict{String, Any}("value" => "1", "type" => "integer")))], "arr_arr_tbl_empty" => Any[Any[Dict{String, Any}()]])


text:
module Datasets

include("Misc.jl")
using .Misc

import NPZ, CSV
import DataFrames: dropmissing!, names!

export DataDesc, ngram, load

struct DataDesc
    path::String
    features::Array{Symbol}
    dims::Array{Int}
    types::Array{DataType}
end

descriptions = Dict(
    "Abalone"       => DataDesc("data/abalone/abalone.data",
                         [:sex, :length, :diameter, :height, :whole, :shucked, :viscera, :shell, :rings],
                         [3,10,10,10,10,10,10,10,10],
                         [String,Float64,Float64,Float64,Float64,Float64,Float64,Float64,Float64]),

    "Contraception" => DataDesc("data/contraception/cmc.data",
                         [:age, :education, :heducation, :children, :religion, :working, :hoccupation, :living, :exposure, :method],
                         [4,4,4,4,2,2,4,4,2,3],
                         [Float64,Int,Int,Float64,Int,Int,Int,Int,Int,Int]),

    "Citation"      => DataDesc("data/citation/citation.data",
                         [:source, :target],
                         [2555,2555],
                         [Int,Int]),

    "LastFM"        => DataDesc("data/lastfm/user_artists.dat",
                         [:user, :artist],
                         [1892,17632],
                         [String,String])
    )

function load(name::String)
    desc = descriptions[name]
    df = CSV.read(desc.path; delim=',', datarow=1, types=desc.types)
    dropmissing!(df; disallowmissing=true)

    T,N = size(df)
    tab = Array{Int,2}(undef, T, N)

    for n in 1:N
        tab[:,n] .= categorize(df[n],desc.dims[n])
    end

    return cooccurrences(tab,desc.dims)
end

function load(name::String,features::Array{Symbol})
    desc = descriptions[name]
    df = CSV.read(desc.path; delim=',', datarow=1, types=desc.types)
    dropmissing!(df; disallowmissing=true)
    names!(df,desc.features)

    T,N = size(df)
    tab = Array{Int,2}(undef, T, length(features))
    dims = Array{Int,1}(undef, length(features))

    for n in 1:length(features)
        i = findfirst(desc.features .== features[n])
        tab[:,n] .= categorize(df[i],desc.dims[i])
        dims[n] = desc.dims[i]
    end

    return cooccurrences(tab,dims)
end

function load(name::String,features::Array{String})
    return load(name, Symbol.(features))
end

function load(name::String,features::AbstractArray{Int})
    return load(name, descriptions[name].features[features])
end

function ngram(n::Int, T::Int)
    freq = NPZ.npzread("data/letters/transitions_$(n)L.npy");
    θ = freq ./ sum(freq)
    return finucan(Float64(T),θ)
end

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
# the functions are somewhat trivial so the tests are also a bit trivial

Random.seed!(1240)

d = 10
normal = randn(d)
v1 = randn(d)
v2 = copy(v1)
II = diagm(0=>ones(d))

@test norm(reflect_bps!(normal,v1)-reflect_bps!(normal,v2,II)) <= 1e-12

v3 = v1-2.0dot(normal,v1)*normal/norm(normal)^2

@test norm(reflect_bps!(normal,v1)-v3) <= 1e-12

mask     = rand(1:d, 3)
v4       = copy(v1)
v4[mask] = v4[mask]*(-1.0)

@test norm(reflect_zz!(mask,v1) - v4) <= 1e-12

v5          = copy(v1)
v5[mask[1]] = -v5[mask[1]]

@test norm(reflect_zz!(mask[1], v1) - v5) <= 1e-12

Random.seed!(12); v  = randn(d)
Random.seed!(12); v1 = refresh_global!(v1)

@test norm(v-v1) <= 1e-12

Random.seed!(31); v  = randn(d); v /= norm(v)
Random.seed!(31); v1 = refresh_restricted!(v1)

@test norm(v-v1) <= 1e-12

beta  = Beta(abs(randn()))
Random.seed!(53)
w  = randn(d)
w /= norm(w)
w *= tan( rand(beta) * 2 * pi )
v  = v1 + w
v /= norm(v)
Random.seed!(53)
v1 = refresh_partial!(beta)(v1)

@test norm(v-v1) <= 1e-12


text:
# --------------------------------------------------------------------------- #
# Abstract type of a synthetic Cell

abstract type Cell end

# --------------------------------------------------------------------------- #
# Simple implementation of a round synthetic Cell that has sharp boundaries

struct SharpCircleCell <: Cell
  radius :: Tuple{Float64, Float64}
  mag :: Tuple{Float64, Float64}
  SharpCircleCell(r, m) = new((r...), (m...))
end

SharpCircleCell(r::Real, m::Real = 1.) = SharpCircleCell((r, r), (m, m))
SharpCircleCell(r, m::Real = 1.) = SharpCircleCell((r...), (m, m))
SharpCircleCell(r::Real, m) = SharpCircleCell((r, r), (m...))

function (c::SharpCircleCell)()
  s = rand(2)
  r = s[1] * c.radius[1] + (1 - s[1]) * c.radius[2]
  m = s[2] * c.mag[1]    + (1 - s[2]) * c.mag[2]
  l = floor(Int, ceil(2 * r) / 2) * 2 + 1
  k = ceil(Int, l / 2)
  patch = zeros(l, l)

  for i in 1:l, j in 1:l
    patch[i, j] = (i - k)^2 + (j - k)^2 <= r^2 ? m : 0
  end

  return patch
end

# --------------------------------------------------------------------------- #
# Simple implementation of an elliptic Cell that has sharp boundaries

struct SharpEllipticCell <: Cell
  axes :: Tuple{Float64, Float64}
  mag :: Tuple{Float64, Float64}
  SharpEllipticCell(r, m) = new((r...), (m...))
end

SharpEllipticCell(r::Real, m::Real = 1.) = SharpEllipticCell((r, r), (m, m))
SharpEllipticCell(r, m::Real = 1.) = SharpEllipticCell((r...), (m, m))
SharpEllipticCell(r::Real, m) = SharpEllipticCell((r, r), (m...))

function (c::SharpEllipticCell)()
  s = rand(4)
  a = s[1] * c.axes[1] + (1 - s[1]) * c.axes[2]
  b = s[2] * c.axes[1] + (1 - s[2]) * c.axes[2]
  m = s[3] * c.mag[1]  + (1 - s[3]) * c.mag[2]
  t = s[4] * 2pi

  l = floor(Int, ceil(2 * max(a, b)) / 2) * 2 + 1
  k = ceil(Int, l / 2)
  patch = zeros(l, l)

  for i in 1:l, j in 1:l
    v = [ cos(t) (-sin(t)); sin(t) cos(t) ] * [ i - k, j - k ]
    patch[i, j] = v[1]^2 / a^2 + v[2]^2 / b^2 <= 1 ? m : 0
  end

  return patch
end

# --------------------------------------------------------------------------- #
# Function that creates synthesized labeled images for provided Cell types

function synthesize( width, height, n :: Tuple{Integer, Integer};
                     cell = SharpCircleCell((5,10)), op = Id() )

  # Initialize image with background and prepare label
  image = zeros(Float32, width, height)

  n = rand(n[1]:n[2])
  label = zeros(Int, 2, n)

  # Iterate over all cell-realizations
  for i in 1:n
    c = cell()
    s = size(c)
    lx = floor(Int, s[1] / 2) + 1
    ly = floor(Int, s[2] / 2) + 1
    x, y = rand(lx:height-lx), rand(ly:width-ly)
    label[:, i] = [x, y]
    xreg = x-lx+1 : x+lx-1
    yreg = y-ly+1 : y+ly-1
    image[yreg, xreg] += c
  end

  lbl = Label(label)
  img = GreyscaleImage(image)

  # Apply post-processing steps and return
  return apply(op, img, lbl)
end

synthesize( w, h, n :: Integer; kwargs...) = synthesize(w, h, (n, n); kwargs...)

max_stars_repo_name: tscode/DCellC

max_issues_repo_name: tscode/DCellC


text:
struct Fiber
    D::Union{Float32} # dispersion factor
    α # exponential power loss
    γ # nonlinear Kerr factor
    L # fiber length
end

struct Bundle
    elems::Vector{Fiber}
    Ns::Union{Int, Vector{Int}}
end

struct SimulationConfig
    SS::Bool
    SRS::Bool
    SBS::Bool
end

default_config = SimulationConfig(false, false, false)

struct Simulation # Simulation struct containing all needed data
    Ψ::Matrix{Complex{Float32}}
    dz
    dt
    Nₜ::Int
    Nₗ::Int
    D
    γ
    α
    L
    T
    λ
    β2
    t::Vector{Float32}
    conf::SimulationConfig
end


text:
@safetestset "HTML+Julia rendering" begin

  @safetestset "Simple tag rendering" begin
    using Genie
    using Genie.Renderer.Html

    copy = "Hello Genie"

    @test Html.p(copy) == "
$copy
"
    @test Html.div() == "

"
    @test Html.br() ==  "
"

    message = "Important message"
    @test Html.span(message, class = "focus") == "$message"
    @test Html.span(message, class = "focus"; NamedTuple{(Symbol("data-process"),)}(("pre-process",))...) == "Important message"
    @test Html.span("Important message", class = "focus"; NamedTuple{(Symbol("data-process"),)}(("pre-process",))...) do
            Html.a("Click here to read message")
          end == "Click here to read message"
  end;

end;


text:
using Test

include("../src/types.jl")
include("../src/meshes.jl")
include("../src/intersecttest.jl")
include("../src/hyperslicesets.jl")
include("../src/hyperslice.jl")

# using MdSlicing

QHypersliceSeg = Tuple{Array{Float64,1},Int64,Int64,Tuple{Float64,Float64},Tuple{Float64,Float64}}

function Base.isapprox(x::QHypersliceSeg, y::QHypersliceSeg)
  all([
    x[1] ≈ y[1],
    x[2] ≈ y[2],
    x[3] ≈ y[3],
    x[4][1] ≈ y[4][1],
    x[4][2] ≈ y[4][2],
    x[5][1] ≈ y[5][1],
    x[5][2] ≈ y[5][2]
  ])
end

function approxin(x, y::Set)
  for e in y
    if x ≈ e return true end
  end
  return false
end

function Base.isapprox(x::Set, y::Set)
  if length(x) != length(y) return false end

  # subset
  all([approxin(e, y) for e in x])
end

function Base.isapprox(x::Tuple{Float64,Float64,Float64,Float64}, y::Tuple{Float64,Float64,Float64,Float64})
  all([
    x[1] ≈ y[1]
    x[2] ≈ y[2]
    x[3] ≈ y[3]
    x[4] ≈ y[4]
  ])
end

function simplifyseg(s::HypersliceSegment)
  (s.fp, Int64(s.d1), Int64(s.d2), (s.p1d1, s.p1d2), (s.p2d1, s.p2d2))
end

function checksegs(v::Vector{Intersect2D}, exp::Vector{QHypersliceSeg})
  v2 = Set([(x.p1d1, x.p1d2, x.p2d1, x.p2d2) for x in v])
  e2 = Set([(y[4][1], y[4][2], y[5][1], y[5][2]) for y in exp])
  @test v2 ≈ e2
end

function checksegs(v::SliceSet2D, exp::Vector{QHypersliceSeg})
  cv = Set(map(simplifyseg, v.slices))
  @test cv ≈ Set(exp)
end

function checksegs(v::Vector{Tuple{Simplex,Intersect2D}}, exp::Vector{QHypersliceSeg})
  i2ds = [x[2] for x in v]
  checksegs(i2ds, exp)
end

@testset "Type sanity checks" begin
  trimesh = ConvexMesh(
    ProblemSpec("x1" => (0., 1.), "x2" => (0., 1.), "x3" => (0., 1.)),
    [ 0. 0. 0. ;
      0. 1. 0. ;
      1. 0. 1.
    ],
    [ 1 2 3 ]  )

  @inferred Vector{Tuple{Simplex,Intersect2D}} slice(trimesh, [0.5, 0.5, 0.5], Dim(1), Dim(2))
end

@testset "One triangle mesh" begin
  trimesh = ConvexMesh(
    ProblemSpec("x1" => (0., 1.), "x2" => (0., 1.), "x3" => (0., 1.)),
    [ 0. 0. 0. ;
      0. 1. 0. ;
      1. 0. 1.
    ],
    [ 1 2 3 ]  )

  checksegs(slice2d(trimesh, [[1.5, 1.5, 1.5]]), QHypersliceSeg[])
  # TODO: define macro to compare slices from expected slices
  checksegs(slice2d(trimesh, [[0.5, 0.5, 0.5]]),
    [ ([0.5, 0.5, 0.5], 1, 2, (0.5, 0.0), (0.5, 0.5)),
      ([0.5, 0.5, 0.5], 1, 3, (0.0, 0.0), (0.5, 0.5)),
      ([0.5, 0.5, 0.5], 2, 3, (0.0, 0.5), (0.5, 0.5)) ])
  checksegs(slice2d(trimesh, [[0.5, 0.5, 0.8]]),
    [ ([0.5, 0.5, 0.8], 1, 2, (0.8, 0.0), (0.8, 0.2)),
      ([0.5, 0.5, 0.8], 1, 3, (0.0, 0.0), (0.5, 0.5)),
      ([0.5, 0.5, 0.8], 2, 3, (0.0, 0.5), (0.5, 0.5)) ])

  checksegs(slice(trimesh, [1.5, 1.5, 1.5], Dim(1), Dim(2)), QHypersliceSeg[])
  checksegs(slice(trimesh, [1.5, 1.5, 1.5], Dim(1), Dim(3)), QHypersliceSeg[])
  checksegs(slice(trimesh, [1.5, 1.5, 1.5], Dim(2), Dim(3)), QHypersliceSeg[])
end

@testset "3D cube" begin
  cube = ConvexMesh(
    ProblemSpec("x1" => (0., 1.), "x2" => (0., 1.), "x3" => (0., 1.)),
    [ 0. 0. 0. ;
      0. 0. 1. ;
      0. 1. 0. ;
      0. 1. 1. ;
      1. 0. 0. ;
      1. 0. 1. ;
      1. 1. 0. ;
      1. 1. 1.
    ],
    [ 4  3  1 ;
      4  2  1 ;
      6  2  1 ;
      6  5  1 ;
      6  4  2 ;
      6  4  8 ;
      7  3  1 ;
      7  5  1 ;
      7  4  3 ;
      7  4  8 ;
      7  6  5 ;
      7  6  8 ]
  )

  checksegs(slice(cube, [0.5, 0.5, 0.5], Dim(1), Dim(2)),
    [ ([0.5, 0.5, 0.5], 1, 2, (0.0, 0.0), (0.5, 0.0)),
      ([0.5, 0.5, 0.5], 1, 2, (0.5, 0.0), (1.0, 0.0)),
      ([0.5, 0.5, 0.5], 1, 2, (1.0, 0.0), (1.0, 0.5)),
      ([0.5, 0.5, 0.5], 1, 2, (1.0, 0.5), (1.0, 1.0)),
      ([0.5, 0.5, 0.5], 1, 2, (0.0, 1.0), (0.5, 1.0)),
      ([0.5, 0.5, 0.5], 1, 2, (0.5, 1.0), (1.0, 1.0)),
      ([0.5, 0.5, 0.5], 1, 2, (0.0, 0.0), (0.0, 0.5)),
      ([0.5, 0.5, 0.5], 1, 2, (0.0, 0.5), (0.0, 1.0)) ]
  )
  checksegs(slice2d(cube, [[0.5, 0.5, 0.5]]),
    [ ([0.5, 0.5, 0.5], 1, 2, (0.0, 0.0), (0.5, 0.0)),
      ([0.5, 0.5, 0.5], 1, 2, (0.5, 0.0), (1.0, 0.0)),
      ([0.5, 0.5, 0.5], 1, 2, (1.0, 0.0), (1.0, 0.5)),
      ([0.5, 0.5, 0.5], 1, 2, (1.0, 0.5), (1.0, 1.0)),
      ([0.5, 0.5, 0.5], 1, 2, (0.0, 1.0), (0.5, 1.0)),
      ([0.5, 0.5, 0.5], 1, 2, (0.5, 1.0), (1.0, 1.0)),
      ([0.5, 0.5, 0.5], 1, 2, (0.0, 0.0), (0.0, 0.5)),
      ([0.5, 0.5, 0.5], 1, 2, (0.0, 0.5), (0.0, 1.0)),
      ([0.5, 0.5, 0.5], 1, 3, (0.0, 0.0), (0.5, 0.0)),
      ([0.5, 0.5, 0.5], 1, 3, (0.5, 0.0), (1.0, 0.0)),
      ([0.5, 0.5, 0.5], 1, 3, (1.0, 0.0), (1.0, 0.5)),
      ([0.5, 0.5, 0.5], 1, 3, (1.0, 0.5), (1.0, 1.0)),
      ([0.5, 0.5, 0.5], 1, 3, (0.0, 1.0), (0.5, 1.0)),
      ([0.5, 0.5, 0.5], 1, 3, (0.5, 1.0), (1.0, 1.0)),
      ([0.5, 0.5, 0.5], 1, 3, (0.0, 0.0), (0.0, 0.5)),
      ([0.5, 0.5, 0.5], 1, 3, (0.0, 0.5), (0.0, 1.0)),
      ([0.5, 0.5, 0.5], 2, 3, (0.0, 0.0), (0.5, 0.0)),
      ([0.5, 0.5, 0.5], 2, 3, (0.5, 0.0), (1.0, 0.0)),
      ([0.5, 0.5, 0.5], 2, 3, (1.0, 0.0), (1.0, 0.5)),
      ([0.5, 0.5, 0.5], 2, 3, (1.0, 0.5), (1.0, 1.0)),
      ([0.5, 0.5, 0.5], 2, 3, (0.0, 1.0), (0.5, 1.0)),
      ([0.5, 0.5, 0.5], 2, 3, (0.5, 1.0), (1.0, 1.0)),
      ([0.5, 0.5, 0.5], 2, 3, (0.0, 0.0), (0.0, 0.5)),
      ([0.5, 0.5, 0.5], 2, 3, (0.0, 0.5), (0.0, 1.0)) ]
  )

end


text:
module node_partitioning_examples
using FinEtools
using UnicodePlots
using Test
import LinearAlgebra: norm

colors = [:magenta, :white, :red, :green, :blue, :yellow, :orange, :brown,  ]

function node_partitioning_1()

    H = 100. # strip width
    a = 10. # radius of the hole
    L = 200. # length of the strip
    nL = 15
    nH = 10
    nR = 50
    fens,fes = Q4elliphole(a, a, L/2, H/2, nL, nH, nR)

    npartitions = 4
    partitioning = nodepartitioning(fens, npartitions)
    partitionnumbers = unique(partitioning)

    # set(axis="normal", plotstyle="points", linewidth=2, pointsize = 1, color = "black", xlabel = "x", ylabel = "y ", grid="on", title = "")
    ixxxx = 1
    i1 = findall(x -> x == partitionnumbers[ixxxx], partitioning)
    pl = scatterplot(vec(fens.xyz[i1, 1]), vec(fens.xyz[i1, 2]), color = colors[ixxxx], xlim = [0 L/2], ylim = [0 H/2])
    for ixxxx = 2:length(partitionnumbers)
        i1 = findall(x -> x == partitionnumbers[ixxxx], partitioning)
        scatterplot!(pl, vec(fens.xyz[i1, 1]), vec(fens.xyz[i1, 2]), color = colors[ixxxx])
    end
    display(pl)

    # Render the  points as dots ("Points" in paraview)
    Points = FESetP1(reshape(collect(1:count(fens)), count(fens), 1))
    File = "node_partitioning_1.vtk"
    vtkexportmesh(File, fens, Points;
    scalars=[("Partition", vec(partitioning))])
    @async run(`"paraview.exe" $File`)

end # node_partitioning_1

function node_partitioning_2()
    a = 10. # radius of the hole
    nC = 20
    nR = 4
    fens,fes = Q4annulus(a, 1.5*a, nR, nC, 2*pi)

    npartitions = 8
    partitioning = nodepartitioning(fens, npartitions)
    partitionnumbers = unique(partitioning)

    # set(axis="normal", plotstyle="points", linewidth=2, pointsize = 2, color = "black", xlabel = "x", ylabel = "y ", grid="on", title = "")

    ixxxx = 1
    i1 = findall(x -> x == partitionnumbers[ixxxx], partitioning)
    pl = scatterplot(vec(fens.xyz[i1, 1]), vec(fens.xyz[i1, 2]), color = colors[ixxxx])
    for ixxxx = 2:length(partitionnumbers)
    	i1 = findall(x -> x == partitionnumbers[ixxxx], partitioning)
    	scatterplot!(pl, vec(fens.xyz[i1, 1]), vec(fens.xyz[i1, 2]), color = colors[ixxxx])
    end
    display(pl)

end # node_partitioning_2

function node_partitioning_3()
    H = 100. # strip width
    a = 10. # radius of the hole
    L = 200. # length of the strip
    nL = 15
    nH = 10
    nR = 50
    fens,fes = Q4elliphole(a, a, L/2, H/2, nL, nH, nR)

    npartitions = 4
    partitioning = nodepartitioning(fens, npartitions)
    partitionnumbers = unique(partitioning)

    # Render the  points as dots ("Points" in paraview)
    Points = FESetP1(reshape(collect(1:count(fens)), count(fens), 1))
    File = "Q4elliphole_mesh.vtk"
    vtkexportmesh(File, fens, Points;
    scalars=[("Partition", vec(partitioning))])
    @async run(`"paraview.exe" $File`)

end # node_partitioning_3

function node_partitioning_4()
    H = 30. # strip width
    R = 10. # radius of the hole
    L = 20. # length of the strip
    nL = 15
    nH = 10
    nR = 5
    fens,fes = H8block(L, H, R, nL, nH, nR)

    npartitions = 16
    partitioning = nodepartitioning(fens, npartitions)
    partitionnumbers = unique(partitioning)
    @test norm(sort(partitionnumbers) - sort(1:npartitions)) < 1.e-5

    # Visualize partitioning
    for gp = partitionnumbers
      groupnodes = findall(k -> k == gp, partitioning)
      File =  "partition-nodes-$(gp).vtk"
      vtkexportmesh(File, fens, FESetP1(reshape(groupnodes, length(groupnodes), 1)))
    end
    File =  "partition-mesh.vtk"
    vtkexportmesh(File, fens, fes)
    @async run(`"paraview.exe" $File`)
end

function node_partitioning_5()
    H = 30. # strip width
    L = 20. # length of the strip
    nL = 15
    nH = 10
    fens,fes = Q4block(L, H, nL, nH)
    reg1 = selectelem(fens, fes; box = [0.0 L/3 0.0 H], inflate = L/nL/1000)
    fes1 = subset(fes, reg1)
    fes2 = subset(fes, setdiff(1:count(fes), reg1))

    # Partitioning of all the nodes
    partitioning = fill(0, count(fens))

    # Find the partitioning of the nodes in region 1
    nincluded1 = fill(false, count(fens))
    for i = connectednodes(fes1)
        nincluded1[i] = true
    end
    npartitions = 8
    partitioning1 = nodepartitioning(fens, nincluded1, npartitions)
    npartitions1 = maximum(unique(partitioning1))

    # Transfer the partitioning of region 1 into the overall partitioning
    partitioning[nincluded1] = partitioning1[nincluded1]

    # Find the partitioning of the nodes connected to region 2, but not region 1
    nincluded2 = fill(false, count(fens))
    for i = connectednodes(fes2)
        nincluded2[i] = (!nincluded1[i]) && true
    end
    npartitions = 4
    partitioning2 = nodepartitioning(fens, nincluded2, npartitions)
    partitioning2 = partitioning2 .+ npartitions1 # shift by the number of partitions in region 1
    partitionnumbers = unique(vcat(partitioning1, partitioning2))
    npartitions = maximum(partitionnumbers)

    # Transfer the partitioning of region 2 into the overall partitioning
    partitioning[nincluded2] = partitioning2[nincluded2]

    # Visualize partitioning
    for gp = partitionnumbers
      groupnodes = findall(k -> k == gp, partitioning)
      File =  "partition-nodes-$(gp).vtk"
      vtkexportmesh(File, fens, FESetP1(reshape(groupnodes, length(groupnodes), 1)))
    end
    File =  "mesh-1.vtk"
    vtkexportmesh(File, fens, fes1)
    File =  "mesh-2.vtk"
    vtkexportmesh(File, fens, fes2)
    @async run(`"paraview.exe" $File`)
end

function node_partitioning_6()
    H = 30. # strip width
    L = 20. # length of the strip
    nL = 15
    nH = 10
    fens,fes = Q4block(L, H, nL, nH)
    reg1 = selectelem(fens, fes; box = [0.0 L/3 0.0 H], inflate = L/nL/1000)
    reg2 = selectelem(fens, fes; box = [L/3 L 0.0 H/2], inflate = L/nL/1000)
    reg3 = selectelem(fens, fes; box = [L/3 L H/2 H], inflate = L/nL/1000)
    fesarr = vec([subset(fes, reg1) subset(fes, reg2) subset(fes, reg3)])

    # Partitioning of all the nodes
    partitioning = nodepartitioning(fens, fesarr, vec([2 4 2]))
    partitionnumbers = unique(partitioning)

    # Visualize partitioning
    for gp = partitionnumbers
      groupnodes = findall(k -> k == gp, partitioning)
      File =  "partition-nodes-$(gp).vtk"
      vtkexportmesh(File, fens, FESetP1(reshape(groupnodes, length(groupnodes), 1)))
    end
    for i = 1:length(fesarr)
        File =  "mesh-$(i).vtk"
        vtkexportmesh(File, fens, fesarr[i])
    end
    File =  "mesh-1.vtk"
    @async run(`"paraview.exe" $File`)
end

function allrun()
    println("#####################################################")
    println("# node_partitioning_1 ")
    node_partitioning_1()
    println("#####################################################")
    println("# node_partitioning_2 ")
    node_partitioning_2()
    println("#####################################################")
    println("# node_partitioning_3 ")
    node_partitioning_3()
    println("#####################################################")
    println("# node_partitioning_4 ")
    node_partitioning_4()
    println("#####################################################")
    println("# node_partitioning_5 ")
    node_partitioning_5()
    println("#####################################################")
    println("# node_partitioning_6 ")
    node_partitioning_6()
    return true
end # function allrun

end # module node_partitioning_examples


text:
import UncertainData:
    AbstractBinnedResampling,
    AbstractBinnedUncertainValueResampling,
    AbstractBinnedSummarisedResampling,
    AbstractUncertainIndexValueDataset,
    resample

export BinnedDataCausalityTest, causality

"""
    BinnedDataCausalityTest(test::CT, binning::AbstractBinnedUncertainValueResampling, n_realizations::Int)
    BinnedDataCausalityTest(test::CT, binning::AbstractBinnedSummarisedResampling)

A causality test where the data are binned before applying the test. If `binning` is some
binning scheme that returns an uncertain value for each
bin (e.g. `BinnedResampling` or `BinnedWeightedResampling`), then the `test` is
applied `n_realizations` times. If `binning` returns a single value for each bin
(e.g. `BinnedMeanResampling` or `BinnedMeanWeightedResampling`, the `test` is applied
only once.

!!! note

    This method uses [`UncertainData.jl`](https://github.com/kahaaga/UncertainData.jl) for
    uncertainty handling. To use it, first run `using UncertainData` in the Julia console.

## Fields

- **`test::CausalityTest`**. An instance of a causality test, e.g. `VisitationFrequencyTest`,
    `PredictiveAsymmetryTest`, `CrossMappingTest` or `JointDistanceDistributionTest`.

- **`binning::AbstractBinnedResampling`**. An instance of a resampling scheme
    indicating the type of binning, e.g `BinnedResampling`, `BinnedWeightedResampling`,
    `BinnedMeanResampling` or `BinnedMeanWeightedResampling`.

- **`n_realizations::Int`**: The number of independent draws of the binned data set over
    which to compute the causality `test`. Only taken into consideration of `binning` is
    a scheme resulting in each bin being represented by an uncertain value. If bins are
    summarized, then the `test` is always applied only once.

## Examples

### Binned convergent cross mapping test

Let's say we want to bin the data by drawing n_draws = 5000 realisations of
each uncertain data point, then assign the draws to the correct bins, and
finally get a kernel density estimate to the distribution of values in each bin.

```julia
using UncertainData, CausalityTools

grid = 0:10:1000 # left bin edges
n_draw = 5000 # sample each point 5000 times and distribute among bins
binning = BinnedResampling(grid, 5000)

ccm_test = ConvergentCrossMappingTest(timeseries_lengths = 15:5:100 |> collect,
    n_reps = 100, libsize = 100, τ = 1, replace = true)

test = BinnedDataCausalityTest(ccm_test, binning, n_realizations)
```

### Binned transfer entropy test

Let's say we want to bin the data by drawing n_draws = 7500 realisations of
each uncertain data point, then assign the draws to the correct bins, and
finally get a kernel density estimate to the distribution of values in each bin.

```julia
grid = 0:5:1000 # left bin edges
n_draws = 7500
binning = BinnedResampling(grid, n_draws)

# A transfer entropy causality test using the visitation frequency estimator
state_space_binning = RectangularBinning(5)
te_test = VisitationFrequencyTest(binning = state_space_binning, ηs = 1)

test = BinnedDataCausalityTest(te_test, binning, n_realizations)
```
"""
struct BinnedDataCausalityTest{
    CT <: CausalityTest,
    BR <: AbstractBinnedResampling} <: PreprocessedDataCausalityTest{CT}

    test::CT
    binning::BR
    n_realizations::Int
end

function BinnedDataCausalityTest(test::CR, binning::BR) where {
    CR <: CausalityTest, BR <: AbstractBinnedSummarisedResampling}

    BinnedDataCausalityTest(test, binning, 1)
end

export BinnedDataCausalityTest