text:
"""
    tr(r)

Transpose an r-matrix.

### Given ###

- `r`: R-matrix

### Returned ###

- `rt`: Transpose

### Note ###

   It is permissible for r and rt to be the same array.

### Called ###

- `eraCr`: copy r-matrix

"""
function tr(r)
    rt = zeros((3, 3))
    ccall((:eraTr, liberfa), Cvoid,
          (Ptr{Cdouble}, Ptr{Cdouble}),
          r, rt)
    rt
end

"""
    trxpv(r, pv)

Multiply a pv-vector by the transpose of an r-matrix.

### Given ###

- `r`: R-matrix
- `pv`: Pv-vector

### Returned ###

- `trpv`: R * pv

### Note ###

   It is permissible for pv and trpv to be the same array.

### Called ###

- `eraTr`: transpose r-matrix
- `eraRxpv`: product of r-matrix and pv-vector

"""
function trxpv(r, p)
    rp = zeros((2, 3))
    ccall((:eraTrxpv, liberfa), Cvoid,
            (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
            r, p, rp)
    rp
end

"""
    trxp(r, p)

Multiply a p-vector by the transpose of an r-matrix.

### Given ###

- `r`: R-matrix
- `p`: P-vector

### Returned ###

- `trp`: R * p

### Note ###

   It is permissible for p and trp to be the same array.

### Called ###

- `eraTr`: transpose r-matrix
- `eraRxp`: product of r-matrix and p-vector

"""
function trxp(r, p)
    rp = zeros(3)
    ccall((:eraTrxp, liberfa), Cvoid,
            (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
            r, p, rp)
    rp
end

"""
    taiut1(tai1, tai2, dta)

Time scale transformation:  International Atomic Time, TAI, to
Universal Time, UT1.

### Given ###

- `tai1`, `tai2`: TAI as a 2-part Julian Date
- `dta`: UT1-TAI in seconds

### Returned ###

- `ut11`, `ut12`: UT1 as a 2-part Julian Date

### Notes ###

1. tai1+tai2 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where tai1 is the Julian
   Day Number and tai2 is the fraction of a day.  The returned
   UT11,UT12 follow suit.

2. The argument dta, i.e. UT1-TAI, is an observed quantity, and is
   available from IERS tabulations.

### Reference ###

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992)

"""
taiut1

"""
    tdbtt(tdb1, tdb2, dtr)

Time scale transformation:  Barycentric Dynamical Time, TDB, to
Terrestrial Time, TT.

### Given ###

- `tdb1`, `tdb2`: TDB as a 2-part Julian Date
- `dtr`: TDB-TT in seconds

### Returned ###

- `tt1`, `tt2`: TT as a 2-part Julian Date

### Notes ###

1. tdb1+tdb2 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where tdb1 is the Julian
   Day Number and tdb2 is the fraction of a day.  The returned
   tt1,tt2 follow suit.

2. The argument dtr represents the quasi-periodic component of the
   GR transformation between TT and TCB.  It is dependent upon the
   adopted solar-system ephemeris, and can be obtained by numerical
   integration, by interrogating a precomputed time ephemeris or by
   evaluating a model such as that implemented in the ERFA function
   eraDtdb.   The quantity is dominated by an annual term of 1.7 ms
   amplitude.

3. TDB is essentially the same as Teph, the time argument for the
   JPL solar system ephemerides.

### References ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

- IAU 2006 Resolution 3

"""
tdbtt

"""
    tttdb(tt1, tt2, dtr)

Time scale transformation:  Terrestrial Time, TT, to Barycentric
Dynamical Time, TDB.

### Given ###

- `tt1`, `tt2`: TT as a 2-part Julian Date
- `dtr`: TDB-TT in seconds

### Returned ###

- `tdb1`, `tdb2`: TDB as a 2-part Julian Date

### Notes ###

1. tt1+tt2 is Julian Date, apportioned in any convenient way between
   the two arguments, for example where tt1 is the Julian Day Number
   and tt2 is the fraction of a day.  The returned tdb1,tdb2 follow
   suit.

2. The argument dtr represents the quasi-periodic component of the
   GR transformation between TT and TCB.  It is dependent upon the
   adopted solar-system ephemeris, and can be obtained by numerical
   integration, by interrogating a precomputed time ephemeris or by
   evaluating a model such as that implemented in the ERFA function
   eraDtdb.   The quantity is dominated by an annual term of 1.7 ms
   amplitude.

3. TDB is essentially the same as Teph, the time argument for the JPL
   solar system ephemerides.

### References ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

- IAU 2006 Resolution 3

"""
tttdb

"""
    ttut1(tt1, tt2, dt)

Time scale transformation:  Terrestrial Time, TT, to Universal Time,
UT1.

### Given ###

- `tt1`, `tt2`: TT as a 2-part Julian Date
- `dt`: TT-UT1 in seconds

### Returned ###

- `ut11`, `ut12`: UT1 as a 2-part Julian Date

### Notes ###

1. tt1+tt2 is Julian Date, apportioned in any convenient way between
   the two arguments, for example where tt1 is the Julian Day Number
   and tt2 is the fraction of a day.  The returned ut11,ut12 follow
   suit.

2. The argument dt is classical Delta T.

### Reference ###

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992)

"""
ttut1

for name in ("taiut1",
             "tdbtt",
             "tttdb",
             "ttut1")
    f = Symbol(name)
    fc = "era" * uppercasefirst(name)
    @eval begin
        function ($f)(a, b, c)
            r1 = Ref(0.0)
            r2 = Ref(0.0)
            i = ccall(($fc, liberfa), Cint,
                      (Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
                      a, b, c, r1, r2)
            @assert i == 0
            r1[], r2[]
        end
    end
end

"""
    taitt(tai1, tai2)

Time scale transformation:  International Atomic Time, TAI, to
Terrestrial Time, TT.

### Given ###

- `tai1`, `tai2`: TAI as a 2-part Julian Date

### Returned ###

- `tt1`, `tt2`: TT as a 2-part Julian Date

### Note ###

   tai1+tai2 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where tai1 is the Julian
   Day Number and tai2 is the fraction of a day.  The returned
   tt1,tt2 follow suit.

### References ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992)

"""
taitt

"""
    taiutc(tai1, tai2)

Time scale transformation:  International Atomic Time, TAI, to
Coordinated Universal Time, UTC.

### Given ###

- `tai1`, `tai2`: TAI as a 2-part Julian Date (Note 1)

### Returned ###

- `utc1`, `utc2`: UTC as a 2-part quasi Julian Date (Notes 1-3)

### Notes ###

1. tai1+tai2 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where tai1 is the Julian
   Day Number and tai2 is the fraction of a day.  The returned utc1
   and utc2 form an analogous pair, except that a special convention
   is used, to deal with the problem of leap seconds - see the next
   note.

2. JD cannot unambiguously represent UTC during a leap second unless
   special measures are taken.  The convention in the present
   function is that the JD day represents UTC days whether the
   length is 86399, 86400 or 86401 SI seconds.  In the 1960-1972 era
   there were smaller jumps (in either direction) each time the
   linear UTC(TAI) expression was changed, and these "mini-leaps"
   are also included in the ERFA convention.

3. The function eraD2dtf can be used to transform the UTC quasi-JD
   into calendar date and clock time, including UTC leap second
   handling.

4. The warning status "dubious year" flags UTCs that predate the
   introduction of the time scale or that are too far in the future
   to be trusted.  See eraDat for further details.

### Called ###

- `eraUtctai`: UTC to TAI

### References ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992)

"""
taiutc

"""
    tcbtdb(tcb1, tcb2)

Time scale transformation:  Barycentric Coordinate Time, TCB, to
Barycentric Dynamical Time, TDB.

### Given ###

- `tcb1`, `tcb2`: TCB as a 2-part Julian Date

### Returned ###

- `tdb1`, `tdb2`: TDB as a 2-part Julian Date

### Notes ###

1. tcb1+tcb2 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where tcb1 is the Julian
   Day Number and tcb2 is the fraction of a day.  The returned
   tdb1,tdb2 follow suit.

2. The 2006 IAU General Assembly introduced a conventional linear
   transformation between TDB and TCB.  This transformation
   compensates for the drift between TCB and terrestrial time TT,
   and keeps TDB approximately centered on TT.  Because the
   relationship between TT and TCB depends on the adopted solar
   system ephemeris, the degree of alignment between TDB and TT over
   long intervals will vary according to which ephemeris is used.
   Former definitions of TDB attempted to avoid this problem by
   stipulating that TDB and TT should differ only by periodic
   effects.  This is a good description of the nature of the
   relationship but eluded precise mathematical formulation.  The
   conventional linear relationship adopted in 2006 sidestepped
   these difficulties whilst delivering a TDB that in practice was
   consistent with values before that date.

3. TDB is essentially the same as Teph, the time argument for the
   JPL solar system ephemerides.

### Reference ###

- IAU 2006 Resolution B3

"""
tcbtdb

"""
    tcgtt(tcg1, tcg2)

Time scale transformation:  Geocentric Coordinate Time, TCG, to
Terrestrial Time, TT.

### Given ###

- `tcg1`, `tcg2`: TCG as a 2-part Julian Date

### Returned ###

- `tt1`, `tt2`: TT as a 2-part Julian Date

### Note ###

   tcg1+tcg2 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where tcg1 is the Julian
   Day Number and tcg22 is the fraction of a day.  The returned
   tt1,tt2 follow suit.

### References ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),.
    IERS Technical Note No. 32, BKG (2004)

- IAU 2000 Resolution B1.9

"""
tcgtt

"""
    tdbtcb(tdb1, tdb2)

Time scale transformation:  Barycentric Dynamical Time, TDB, to
Barycentric Coordinate Time, TCB.

### Given ###

- `tdb1`, `tdb2`: TDB as a 2-part Julian Date

### Returned ###

- `tcb1`, `tcb2`: TCB as a 2-part Julian Date

### Notes ###

1. tdb1+tdb2 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where tdb1 is the Julian
   Day Number and tdb2 is the fraction of a day.  The returned
   tcb1,tcb2 follow suit.

2. The 2006 IAU General Assembly introduced a conventional linear
   transformation between TDB and TCB.  This transformation
   compensates for the drift between TCB and terrestrial time TT,
   and keeps TDB approximately centered on TT.  Because the
   relationship between TT and TCB depends on the adopted solar
   system ephemeris, the degree of alignment between TDB and TT over
   long intervals will vary according to which ephemeris is used.
   Former definitions of TDB attempted to avoid this problem by
   stipulating that TDB and TT should differ only by periodic
   effects.  This is a good description of the nature of the
   relationship but eluded precise mathematical formulation.  The
   conventional linear relationship adopted in 2006 sidestepped
   these difficulties whilst delivering a TDB that in practice was
   consistent with values before that date.

3. TDB is essentially the same as Teph, the time argument for the
   JPL solar system ephemerides.

### Reference ###

- IAU 2006 Resolution B3

"""
tdbtcb

"""
    tttai(tt1, tt2)

Time scale transformation:  Terrestrial Time, TT, to International
Atomic Time, TAI.

### Given ###

- `tt1`, `tt2`: TT as a 2-part Julian Date

### Returned ###

- `tai1`, `tai2`: TAI as a 2-part Julian Date

### Note ###

   tt1+tt2 is Julian Date, apportioned in any convenient way between
   the two arguments, for example where tt1 is the Julian Day Number
   and tt2 is the fraction of a day.  The returned tai1,tai2 follow
   suit.

### References ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

- Explanatory Supplement to the Astronomical Almanac,
    P. Kenneth Seidelmann (ed), University Science Books (1992)

"""
tttai

"""
    tttcg(tt1, tt2)

Time scale transformation:  Terrestrial Time, TT, to Geocentric
Coordinate Time, TCG.

### Given ###

- `tt1`, `tt2`: TT as a 2-part Julian Date

### Returned ###

- `tcg1`, `tcg2`: TCG as a 2-part Julian Date

### Note ###

   tt1+tt2 is Julian Date, apportioned in any convenient way between
   the two arguments, for example where tt1 is the Julian Day Number
   and tt2 is the fraction of a day.  The returned tcg1,tcg2 follow
   suit.

### References ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

- IAU 2000 Resolution B1.9

"""
tttcg

for name in ("taitt",
             "taiutc",
             "tcbtdb",
             "tcgtt",
             "tdbtcb",
             "tttai",
             "tttcg")
    f = Symbol(name)
    fc = "era" * uppercasefirst(name)
    @eval begin
        function ($f)(a, b)
            r1 = Ref(0.0)
            r2 = Ref(0.0)
            i = ccall(($fc, liberfa), Cint,
                      (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
                      a, b, r1, r2)
            @assert i == 0
            r1[], r2[]
        end
    end
end

"""
    tf2a(s, ihour, imin, sec)

Convert hours, minutes, seconds to radians.

### Given ###

- `s`: Sign:  '-' = negative, otherwise positive
- `ihour`: Hours
- `imin`: Minutes
- `sec`: Seconds

### Returned ###

- `rad`: Angle in radians

### Notes ###

1.  The result is computed even if any of the range checks fail.

2.  Negative ihour, imin and/or sec produce a warning status, but
    the absolute value is used in the conversion.

3.  If there are multiple errors, the status value reflects only the
    first, the smallest taking precedence.

"""
tf2a

"""
    tf2d(s, ihour, imin, sec)

Convert hours, minutes, seconds to days.

### Given ###

- `s`: Sign:  '-' = negative, otherwise positive
- `ihour`: Hours
- `imin`: Minutes
- `sec`: Seconds

### Returned ###

- `days`: Interval in days

### Notes ###

1.  The result is computed even if any of the range checks fail.

2.  Negative ihour, imin and/or sec produce a warning status, but
    the absolute value is used in the conversion.

3.  If there are multiple errors, the status value reflects only the
    first, the smallest taking precedence.

"""
tf2d

for name in ("tf2a",
             "tf2d")
    f = Symbol(name)
    fc = "era" * uppercasefirst(name)
    @eval begin
        function ($f)(s, ideg, iamin, asec)
            rad = Ref(0.0)
            i = ccall(($fc, liberfa), Cint,
                      (Cchar, Cint, Cint, Cdouble, Ref{Cdouble}),
                       s, ideg, iamin, asec, rad)
            if i == 1
                throw(ERFAException("ihour outside range 0-23"))
            elseif i == 2
                throw(ERFAException("imin outside range 0-59"))
            elseif i == 3
                throw(ERFAException("sec outside range 0-59.999..."))
            end
            rad[]
        end
    end
end


text:
export ConstantWidthKernel

"SPH smoothing kernel with constant width."
struct ConstantWidthKernel{N,K} <: Kernel{N}
    W::K
    h::Float
    h⁻¹::Float
    h⁻ᴺ::Float
    h⁻ᴺ⁻¹::Float

    ConstantWidthKernel(W::K, h::Number) where {N,K<:Kernel{N}} =
        @check(h > 0) && new{N,K}(W, h, 1 / h, 1 / h^N, 1 / h^(N + 1))
end

"""
    (W::ConstantWidthKernel)(q::Number) -> Float

Evaluate the kernel with unit width at relative distance `q`.
"""
(wrapper::ConstantWidthKernel)(q::Number) = wrapper.W(q)

"""
    ddr(W::ConstantWidthKernel, q::Number) -> Float

Evaluate the derivative of the kernel with unit width at relative distance `q`.
The derivative is with respect to `r`.
"""
ddr(wrapper::ConstantWidthKernel, q::Number) = ddr(wrapper.W, q)

"""
    ddh(W::ConstantWidthKernel, q::Number) -> Float

Evaluate the derivative of the kernel with unit width at relative distance `q`.
The derivative is with respect to `h`.
"""
ddh(wrapper::ConstantWidthKernel, q::Number) = ddh(wrapper.W, q)

"""
    (W::ConstantWidthKernel)(q::Number, ::Number) -> Float

Evaluate the kernel with constant width `h` at relative distance `q = r/h`.
"""
(W::ConstantWidthKernel)(q::Number, ::Number) = W(q) * W.h⁻ᴺ

"""
    ddr(W::ConstantWidthKernel, q::Number, ::Number) -> Float

Evaluate the derivative of the kernel with constant width `h` at relative
distance `q = r/h`. The derivative is with respect to `r`.
"""
ddr(W::ConstantWidthKernel, q::Number, ::Number) = ddr(W, q) * W.h⁻ᴺ⁻¹

"""
    ddh(W::ConstantWidthKernel, q::Number, ::Number) -> Float

Evaluate the derivative of the kernel with constant width `h` at relative
distance `q = r/h`. The derivative is with respect to `h`.
"""
ddh(W::ConstantWidthKernel, q::Number, ::Number) = ddh(W, q) * W.h⁻ᴺ⁻¹

get_support(wrapper::ConstantWidthKernel) = get_support(wrapper.W)

get_kernel_width(wrapper::ConstantWidthKernel) = wrapper.h

get_inverse_kernel_width(wrapper::ConstantWidthKernel) = wrapper.h⁻¹

ext: jl

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
using DiffEqBase, Catalyst
# Jump Example Problems
export prob_jump_dnarepressor, prob_jump_constproduct, prob_jump_nonlinrxs,
# examples mixing mass action and constant rate jumps
       prob_jump_osc_mixed_jumptypes,
# examples used in published benchmarks / comparisions
       prob_jump_multistate, prob_jump_twentygenes, prob_jump_dnadimer_repressor,
# examples approximating diffusion by continuous time random walks
       prob_jump_diffnetwork

"""
    General structure to hold JumpProblem info. Needed since
    the JumpProblem constructor requires the algorithm, so we
    don't create the JumpProblem here.
"""

struct JumpProblemNetwork
    network         # Catalyst network
    rates           # vector of rate constants or nothing
    tstop           # time to end simulation
    u0              # initial values
    discrete_prob   # initialized discrete problem
    prob_data       # additional problem data, stored as a Dict
end

dna_rs = @reaction_network begin
    k1, DNA --> mRNA + DNA
    k2, mRNA --> mRNA + P
    k3, mRNA --> 0
    k4, P --> 0
    k5, DNA + P --> DNAR
    k6, DNAR --> DNA + P
end k1 k2 k3 k4 k5 k6
rates = [.5, (20*log(2.)/120.), (log(2.)/120.), (log(2.)/600.), .025, 1.]
tf = 1000.0
u0 = [1,0,0,0]
prob = DiscreteProblem(dna_rs, u0, (0.0, tf), rates)
Nsims = 8000
expected_avg = 5.926553750000000e+02
prob_data = Dict("num_sims_for_mean" => Nsims, "expected_mean" => expected_avg)
"""
    DNA negative feedback autoregulatory model. Protein acts as repressor.
"""
prob_jump_dnarepressor = JumpProblemNetwork(dna_rs, rates, tf, u0, prob, prob_data)

bd_rs = @reaction_network begin
    k1, 0 --> A
    k2, A --> 0
end k1 k2
rates = [1000., 10.]
tf = 1.0
u0 = [0]
prob = DiscreteProblem(bd_rs, u0, (0., tf), rates)
Nsims = 16000
expected_avg = t -> rates[1] / rates[2] .* ( 1. - exp.(-rates[2] * t))
prob_data = Dict("num_sims_for_mean" => Nsims, "expected_mean_at_t" => expected_avg)
"""
    Simple birth-death process with constant production and degradation.
"""
prob_jump_constproduct = JumpProblemNetwork(bd_rs, rates, tf, u0, prob, prob_data)

nonlin_rs = @reaction_network begin
    k1, 2A --> B
    k2, B --> 2A
    k3, A + B --> C
    k4, C --> A + B
    k5, 3C --> 3A
end k1 k2 k3 k4 k5
rates = [1., 2., .5, .75, .25]
tf = .01
u0 = [200, 100, 150]
prob = DiscreteProblem(nonlin_rs, u0, (0., tf), rates)
Nsims = 32000
expected_avg = 84.876015624999994
prob_data = Dict("num_sims_for_mean" => Nsims, "expected_mean" => expected_avg)
"""
    Example with a mix of nonlinear reactions, including third order
"""
prob_jump_nonlinrxs = JumpProblemNetwork(nonlin_rs, rates, tf, u0, prob, prob_data)

oscil_rs = @reaction_network begin
    0.01, (X,Y,Z) --> 0
    hill(X,3.,100.,-4), 0 --> Y
    hill(Y,3.,100.,-4), 0 --> Z
    hill(Z,4.5,100.,-4), 0 --> X
    hill(X,2.,100.,6), 0 --> R
    hill(Y,15.,100.,4)*0.002, R --> 0
    20, 0 --> S
    R*0.005, S --> SP
    0.01, SP + SP --> SP2
    0.05, SP2 --> 0
end
u0 = [200.,60.,120.,100.,50.,50.,50.]  # Hill equations force use of floats!
tf = 4000.
prob = DiscreteProblem(oscil_rs, u0, (0.,tf))
"""
    Oscillatory system, uses a mixture of jump types.
"""
prob_jump_osc_mixed_jumptypes = JumpProblemNetwork(oscil_rs, nothing, tf, u0, prob, nothing)

specs_sym_to_name = Dict(
    :S1 => "R(a,l)",
    :S2 => "L(r)",
    :S3 => "A(Y~U,r)",
    :S4 => "L(r!1).R(a,l!1)",
    :S5 => "A(Y~U,r!1).R(a!1,l)",
    :S6 => "A(Y~U,r!1).L(r!2).R(a!1,l!2)",
    :S7 => "A(Y~P,r!1).L(r!2).R(a!1,l!2)",
    :S8 => "A(Y~P,r!1).R(a!1,l)",
    :S9 => "A(Y~P,r)")
rates_sym_to_idx = Dict(
    :R0 => 1,   :L0 => 2,    :A0 => 3,  :kon => 4, :koff => 5,
    :kAon => 6, :kAoff => 7, :kAp => 8, :kAdp => 9)
params = [5360, 1160, 5360, 0.01, 0.1, 0.01, 0.1, 0.01, 0.1]
rs = @reaction_network begin
    kon, S1 + S2 --> S4
    kAon, S1 + S3 --> S5
    kon, S2 + S5 --> S6
    koff, S4 --> S1 + S2
    kAon, S3 + S4 --> S6
    kAoff, S5 --> S1 + S3
    koff, S6 --> S2 + S5
    kAoff, S6 --> S3 + S4
    kAp, S6 --> S7
    koff, S7 --> S2 + S8
    kAoff, S7 --> S4 + S9
    kAdp, S7 --> S6
    kon, S2 + S8 --> S7
    kAon, S1 + S9 --> S8
    kAon, S4 + S9 --> S7
    kAoff, S8 --> S1 + S9
    kAdp, S8 --> S5
    kAdp, S9 --> S3
end kon kAon koff kAoff kAp kAdp
rsi   = rates_sym_to_idx
rates = params[[rsi[:kon], rsi[:kAon], rsi[:koff], rsi[:kAoff], rsi[:kAp], rsi[:kAdp]]]
u0    = zeros(Int,9)
statesyms = ModelingToolkit.tosymbol.(ModelingToolkit.operation.(rs.states))
u0[ findfirst(isequal(:S1), statesyms)] = params[1]
u0[ findfirst(isequal(:S2), statesyms)] = params[2]
u0[ findfirst(isequal(:S3), statesyms)] = params[3]
tf    = 100.
prob = DiscreteProblem(rs, u0, (0., tf), rates)
"""
    Multistate model from Gupta and Mendes,
    "An Overview of Network-Based and -Free Approaches for Stochastic Simulation of Biochemical Systems",
    Computation 2018, 6, 9; doi:10.3390/computation6010009
    Translated from supplementary data file: Models/Multi-state/fixed_multistate.xml
"""
prob_jump_multistate = JumpProblemNetwork(rs, rates, tf, u0, prob,
                Dict("specs_to_sym_name" => specs_sym_to_name, "rates_sym_to_idx" => rates_sym_to_idx, "params" => params))

# generate the network
N = 10  # number of genes
@parameters t
@variables G[1:2N](t) M[1:2N](t) P[1:2N](t) G_ind[1:2N](t)

function construct_genenetwork(N)
    genenetwork = make_empty_network()
    for i in 1:N
        addspecies!(genenetwork,G[2*i-1])
        addspecies!(genenetwork,M[2*i-i])
        addspecies!(genenetwork,P[2*i-i])
        addreaction!(genenetwork, Reaction(10.0, [G[2*i-i]], [G[2*i-i],M[2*i-i]]))
        addreaction!(genenetwork, Reaction(10.0, [M[2*i-i]], [M[2*i-i],P[2*i-i]]))
        addreaction!(genenetwork, Reaction(1.0, [M[2*i-i]], nothing))
        addreaction!(genenetwork, Reaction(1.0, [P[2*i-i]], nothing))
        # genenetwork *= "\t 10.0, G$(2*i-1) --> G$(2*i-1) + M$(2*i-1)\n"
        # genenetwork *= "\t 10.0, M$(2*i-1) --> M$(2*i-1) + P$(2*i-1)\n"
        # genenetwork *= "\t 1.0,  M$(2*i-1) --> 0\n"
        # genenetwork *= "\t 1.0,  P$(2*i-1) --> 0\n"

        addspecies!(genenetwork,G[2*i])
        addspecies!(genenetwork,M[2*i])
        addspecies!(genenetwork,P[2*i])
        addreaction!(genenetwork, Reaction(5.0, [G[2*i]], [G[2*i],M[2*i]]))
        addreaction!(genenetwork, Reaction(5.0, [M[2*i]], [M[2*i],P[2*i]]))
        addreaction!(genenetwork, Reaction(1.0, [M[2*i]], nothing))
        addreaction!(genenetwork, Reaction(1.0, [P[2*i]], nothing))
        # genenetwork *= "\t 5.0, G$(2*i) --> G$(2*i) + M$(2*i)\n"
        # genenetwork *= "\t 5.0, M$(2*i) --> M$(2*i) + P$(2*i)\n"
        # genenetwork *= "\t 1.0,  M$(2*i) --> 0\n"
        # genenetwork *= "\t 1.0,  P$(2*i) --> 0\n"

        addspecies!(genenetwork, G_ind[2*i])
        addreaction!(genenetwork, Reaction(0.0001, [G[2*i],P[2*i-i]], [G_ind[2*i]]))
        addreaction!(genenetwork, Reaction(100.0, [G_ind[2*i]], [G_ind[2*i],M[2*i]]))
        # genenetwork *= "\t 0.0001, G$(2*i) + P$(2*i-1) --> G$(2*i)_ind \n"
        # genenetwork *= "\t 100., G$(2*i)_ind --> G$(2*i)_ind + M$(2*i)\n"
    end
    genenetwork
end
rs = construct_genenetwork(N)
u0 = zeros(Int, length(rs.states))
statesyms = ModelingToolkit.tosymbol.(ModelingToolkit.operation.(rs.states))
for i = 1:(2*N)
    u0[findfirst(isequal(G[i]),rs.states)] = 1
end
tf = 2000.0
prob = DiscreteProblem(rs, u0, (0.0, tf))
"""
    Twenty-gene model from McCollum et al,
    "The sorting direct method for stochastic simulation of biochemical systems with varying reaction execution behavior"
    Comp. Bio. and Chem., 30, pg. 39-49 (2006).
"""
prob_jump_twentygenes = JumpProblemNetwork(rs, nothing, tf, u0, prob, nothing)

rn = @reaction_network begin
    c1, G --> G + M
    c2, M --> M + P
    c3, M --> 0
    c4, P --> 0
    c5, 2P --> P2
    c6, P2 --> 2P
    c7, P2 + G --> P2G
    c8, P2G --> P2 + G
end c1 c2 c3 c4 c5 c6 c7 c8
rnpar = [.09, .05, .001, .0009, .00001, .0005, .005, .9]
varlabels = ["G", "M", "P", "P2","P2G"]
u0 = [1000, 0, 0, 0,0]
tf = 4000.
prob = DiscreteProblem(rn, u0, (0.0, tf), rnpar)
"""
    Negative feedback autoregulatory gene expression model. Dimer is the repressor.
    Taken from Marchetti, Priami and Thanh,
    "Simulation Algorithms for Comptuational Systems Biology",
    Springer (2017).
"""
prob_jump_dnadimer_repressor = JumpProblemNetwork(rn, rnpar, tf, u0, prob,
                                                Dict("specs_names" => varlabels))

# diffusion model
function getDiffNetwork(N)
    diffnetwork = make_empty_network()
    @parameters t K
    @variables X[1:N](t)
    for i = 1:N
        addspecies!(diffnetwork, X[i])
    end
    for i in 1:(N-1)
        addreaction!(diffnetwork, Reaction(K, [X[i]], [X[i+1]]))
        addreaction!(diffnetwork, Reaction(K, [X[i+1]], [X[i]]))
    end
    diffnetwork
end
params = (1.,)
function getDiffu0(N)
    10*ones(Int64, N)
end
tf = 10.
"""
    Continuous time random walk (i.e. diffusion approximation) example.
    Here the network in the JumpProblemNetwork is a function that returns a
    network given the number of lattice sites.
    u0 is a similar function that returns the initial condition vector.
"""
prob_jump_diffnetwork = JumpProblemNetwork(getDiffNetwork, params, tf, getDiffu0, nothing, nothing)


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# src/forward_operators/concave_increasing.jl
# Contains definitions of commonly used activation functions:
# relu, parametric relu, leaky relu, maxsig, elu, selu, gelu, sigmoid,
# bisigmoid, swish, maxtanh, pentanh, softsign, softplus.
#############################################################################

# RELU DEFINITION
#=
relu

The Rectified Linear Unit (ReLU) activation function `relu(x) = max(x, 0.0)`.
=#
relu_deriv(x) = x > 0.0 ? 1.0 : 0.0
relu_deriv2(x) =  0.0
function relu_kernel(x::MC{N,T}, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
    max_kernel(x, 0.0, z)
end
function relu(x::MC{N,T}) where {N, T<:Union{NS,MV}}
    relu_kernel(x, max(x.Intv, 0.0))
end

# PARAMETRIC RELU DEFINITION
"""
param_relu

The parametric Rectified Linear Unit activation function `param_relu(x, α) = (max(x, αx)` with α in [0,1].
"""
param_relu(x, α) = x > 0.0 ? x : α*x
function param_relu(x::Interval{Float64}, α::Float64)
    xL = x.lo
    xU = x.hi
    (xL < 0.0) && (xL *= α)
    (xU < 0.0) && (xU *= α)
    return Interval{Float64}(xL, xU)
end
param_relu(x::Float64, α::Float64) = x > 0.0 ? x : α*x
param_relu_deriv(x::Float64, α::Float64) = x > 0.0 ? 1.0 : α

function param_relu_grad(g, x::Float64, α::Float64)
    if x > 0.0
        g[1] = 1.0
        g[2] = 0.0
    end
    g[1] = α
    g[2] = x
    nothing
end
function param_relu_kernel(x::MC{N,T}, α::Float64, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
    @assert α >= 0.0
    xLc = z.lo
    xUc = z.hi
    xL = x.Intv.lo
    xU = x.Intv.hi
    midcv, cv_id = mid3(x.cc, x.cv, xL)
    midcc, cc_id = mid3(x.cc, x.cv, xU)
    dcc = (xUc > xLc) ? (xUc - xLc)/(xU - xL) : 0.0
    convex = param_relu(midcv, α)
    concave = dcc*(midcc - xL) + xLc
    concave_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    convex_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*param_relu_deriv(midcv, α)
    convex, concave, convex_grad, concave_grad = cut(xLc, xUc, convex, concave, convex_grad, concave_grad)
    return MC{N, T}(convex, concave, z, convex_grad, concave_grad, x.cnst)
end
@inline param_relu(x::MC, α::Float64) = param_relu_kernel(x, α, param_relu(x.Intv, α))

# LEAKY RELU DEFINITION
"""
leaky_relu

The leaky Rectified Linear Unit activation function `leaky_relu(x) = max(x, 0.01x)`.
"""
@inline leaky_relu(x) = param_relu(x, 0.01)
@inline leaky_relu_kernel(x::MC, z::Interval{Float64}) = param_relu_kernel(x, 0.01, z)
@inline leaky_relu(x::MC) = leaky_relu_kernel(x, param_relu(x.Intv, 0.01))
@inline leaky_relu_deriv(x) = x > 0.0 ? 1.0 : 0.01
@inline leaky_relu_deriv2(x) = 0.0

# DEFINE MAXSIG
"""
maxsig

The `maxsig` activation function  `maxsig(x) = max(x, 1.0/(1.0 + exp(-x)))`.
"""
@inline maxsig(x) = max(x, 1.0/(1.0 + exp(-x)))
@inline maxsig(x::Float64) = max(x, 1.0/(1.0 + exp(-x)))
@inline maxsig(x::Interval{Float64}) = max(x, 1.0/(1.0 + exp(-x)))
@inline function maxsig_deriv(x::Float64)
    if x > 1.0/(exp(-x) + 1.0)
        return 1.0
    end
    return exp(-x)/(exp(-x) + 1.0)^2
end
@inline function maxsig_deriv2(x::Float64)
    if x > 1.0/(exp(-x) + 1.0)
        return 0.0
    end
    return exp(x)*(1.0 - exp(x))/(exp(x) + 1.0)^3 #
end
function maxsig_kernel(x::MC{N,T}, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
    xLc = z.lo
    xUc = z.hi
    xL = x.Intv.lo
    xU = x.Intv.hi
    midcv, cv_id = mid3(x.cc, x.cv, xL)
    midcc, cc_id = mid3(x.cc, x.cv, xU)
    dcc = (xUc > xLc) ? (xUc - xLc)/(xU - xL) : 0.0
    convex = maxsig(midcv)
    concave = dcc*(midcc - xL) + xLc
    concave_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    convex_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*maxsig_deriv(midcv)
    convex, concave, convex_grad, concave_grad = cut(xLc, xUc, convex, concave, convex_grad, concave_grad)
    return MC{N, T}(convex, concave, z, convex_grad, concave_grad, x.cnst)
end
maxsig(x::MC{N,T}) where {N, T<:Union{NS,MV}} = maxsig_kernel(x, maxsig(x.Intv))

# DEFINE ELU
#=
elu

The Exponential Linear Unit (ELU) activation function  `elu(x, α) = x > 0 ? x : α*(exp(x) - 1.0)`.
=#
function elu(x::Interval{Float64}, α::Float64)
    xL = x.lo
    xU = x.hi
    if xU < 0.0
        return α*(exp(x) - 1.0)
    elseif xL > 0.0
        return x
    end
    xLIntv = α*(exp(x) - 1.0)
    Interval(xLIntv.lo, x.hi)
end
@inline elu_deriv(x::Float64, α::Float64) = x > 0.0 ? 1.0 : α*exp(x)
function elu_grad(g, x::Float64, α::Float64)
    g[1] = elu_deriv(x, α)
    g[2] = x > 0.0 ? 0.0 : exp(x)
    nothing
end
function elu_kernel(x::MC{N,T}, α::Float64, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
    xLc = z.lo
    xUc = z.hi
    xL = x.Intv.lo
    xU = x.Intv.hi
    midcv, cv_id = mid3(x.cc, x.cv, xL)
    midcc, cc_id = mid3(x.cc, x.cv, xU)
    dcc = (xUc > xLc) ? (xUc - xLc)/(xU - xL) : 0.0
    convex = elu(midcv, α)
    concave = dcc*(midcc - xL) + xLc
    concave_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    convex_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*elu_deriv(midcv, α)
    convex, concave, convex_grad, concave_grad = cut(xLc, xUc, convex, concave, convex_grad, concave_grad)
    return MC{N,T}(convex, concave, z, convex_grad, concave_grad, x.cnst)
end
elu(x::MC{N,T}, α::Float64) where {N, T<:Union{NS,MV}} = elu_kernel(x, α, elu(x.Intv, α))

"""
selu

The Scaled Exponential Linear Unit (SELU) activation function  `selu(x, α, λ) = λ*elu(x, α)`.
"""
selu(x, α, λ) = λ*elu(x, α)
function selu_grad(g, x::Float64, α::Float64, λ::Float64)
    g[1] = λ*elu_deriv(x, α)
    g[2] = λ*(x > 0.0 ? 0.0 : exp(x))
    g[3] = elu(x, α)
    nothing
end

"""
maxtanh

The `maxtanh` activation function  `maxtanh(x) = max(x, tanh(x))`.
"""
@inline maxtanh(x) = max(x, tanh(x))
@inline maxtanh(x::Float64) = max(x, tanh(x))
@inline function maxtanh(x::Interval{Float64})
    xLintv = Interval(x.lo)
    xUintv = Interval(x.hi)
    xLc = max(xLintv, tanh(xLintv))
    xUc = max(xUintv, tanh(xUintv))
    Interval(xLc.lo, xUc.hi)
end
@inline function maxtanh_deriv(x::Float64)
    if x > tanh(x)
        return 1.0
    end
    return sech(x)^2
end
@inline function maxtanh_deriv2(x::Float64)
    if x > tanh(x)
        return 0.0
    end
    return -2.0*tanh(x)*sech(x)^2
end
function maxtanh_kernel(x::MC{N,T}, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
    xLc = z.lo
    xUc = z.hi
    xL = x.Intv.lo
    xU = x.Intv.hi
    midcv, cv_id = mid3(x.cc, x.cv, xL)
    midcc, cc_id = mid3(x.cc, x.cv, xU)
    dcc = (xUc > xLc) ? (xUc - xLc)/(xU - xL) : 0.0
    convex = maxtanh(midcv)
    concave = dcc*(midcc - xL) + xLc
    concave_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    convex_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*maxtanh_deriv(midcv)
    convex, concave, convex_grad, concave_grad = cut(xLc, xUc, convex, concave, convex_grad, concave_grad)
    return MC{N, T}(convex, concave, z, convex_grad, concave_grad, x.cnst)
end
maxtanh(x::MC{N,T}) where {N, T<:Union{NS,MV}} = maxtanh_kernel(x, maxtanh(x.Intv))

#=
softplus

The `softplus` activation function  `softplus(x) = log(1.0 + exp(x))`.
=#
@inline softplus(x::Interval{Float64}) = log(1.0 + exp(x))
@inline softplus_deriv(x::Float64) = 1.0/(exp(-x) + 1.0)
@inline softplus_deriv2(x::Float64) = exp(-x)/(exp(-x) + 1.0)^2
function softplus_kernel(x::MC{N,T}, z::Interval{Float64}) where {N, T<:Union{NS,MV}}
    xLc = z.lo
    xUc = z.hi
    xL = x.Intv.lo
    xU = x.Intv.hi
    midcv, cv_id = mid3(x.cc, x.cv, xL)
    midcc, cc_id = mid3(x.cc, x.cv, xU)
    dcc = (xUc > xLc) ? (xUc - xLc)/(xU - xL) : 0.0
    convex = softplus(midcv)
    concave = dcc*(midcc - xL) + xLc
    concave_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    convex_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*softplus_deriv(midcv)
    convex, concave, convex_grad, concave_grad = cut(xLc, xUc, convex, concave, convex_grad, concave_grad)
    return MC{N, T}(convex, concave, z, convex_grad, concave_grad, x.cnst)
end
softplus(x::MC{N,T}) where {N, T<:Union{NS,MV}} = softplus_kernel(x, softplus(x.Intv))

"""
pentanh

The `pentanh` activation function `pentanh(x) = x > 0.0 ? tanh(x) : tanh(0.25*x)`.
"""
@inline pentanh(x) = x > 0.0 ? tanh(x) : tanh(0.25*x)
@inline pentanh(x::Float64) = x > 0.0 ? tanh(x) : tanh(0.25*x)
function pentanh(x::Interval{Float64})
    (x.lo >= 0.0) && return tanh(x)
    (x.hi <= 0.0) && return tanh(0.25*x)
    lo_part = tanh(0.25*x)
    hi_part = tanh(x)
    Interval(lo_part.lo, hi_part.hi)
end
function pentanh_deriv(x::Float64)
    if x > 0.0
        return 1.0 - tanh(x)^2
    end
    0.25 - 0.25*tanh(0.25*x)^2
end
function pentanh_deriv2(x::Float64)
    if x > 0.0
        return -2.0*tanh(x)*sech(2)^2
    end
    -0.125*tanh(0.25*x)*sech(0.25*x)^2
end
@inline function pentanh_env(x::Float64, y::Float64, z::Float64)
    (x - y) - (pentanh(x) - pentanh(y))/pentanh_deriv(x)
end
@inline function cv_pentanh(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return dline_seg(pentanh, pentanh_deriv, x, xL, xU)..., p)
    (xU <= 0.0) && (return pentanh(x), pentanh_deriv(x), p)
    if p == Inf
        p, flag = secant(xL, 0.0, xL, 0.0, pentanh_env, xU, 0.0)
        flag && (p = golden_section(xL, 0.0, pentanh_env, xU, 0.0))
    end
    (x <= p) && (return pentanh(x), pentanh_deriv(x), p)
    return dline_seg(pentanh, pentanh_deriv, x, p, xU)..., p
end
@inline function cc_pentanh(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return pentanh(x), pentanh_deriv(x), p)
    (xU <= 0.0) && (return dline_seg(pentanh, pentanh_deriv, x, xL, xU)..., p)
    if p == Inf
        p, flag = secant(0.0, xU, 0.0, xU, pentanh_env, xL, 0.0)
        flag && (p = golden_section(0.0, xU, pentanh_env, xL, 0.0))
    end
    (x <= p) && (return dline_seg(pentanh, pentanh_deriv, x, xL, p)..., p)
    return pentanh(x), pentanh_deriv(x), p
end

#=
sigmoid

The `sigmoid` activation function `sigmoid(x) = 1.0/(1.0 + exp(-x))`.
=#
@inline sigmoid(x::Interval{Float64}) = 1.0/(1.0 + exp(-x))
@inline sigmoid_deriv(x::Float64) = exp(-x)/(1.0 + exp(-x))^2
@inline sigmoid_deriv2(x::Float64) = -exp(x)*(exp(x) - 1.0)/(exp(x) + 1.0)^3
@inline function sigmoid_env(x::Float64, y::Float64, z::Float64)
    (x - y)*sigmoid_deriv(x) - (sigmoid(x) - sigmoid(y))
end
@inline function cv_sigmoid(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return dline_seg(sigmoid, sigmoid_deriv, x, xL, xU)..., p)
    (xU <= 0.0) && (return sigmoid(x), sigmoid_deriv(x), p)
    if p == Inf
        p, flag = secant(xL, 0.0, xL, 0.0, sigmoid_env, xU, 0.0)
        flag && (p = golden_section(xL, 0.0, sigmoid_env, xU, 0.0))
    end
    (x <= p) && (return sigmoid(x), sigmoid_deriv(x), p)
    return dline_seg(sigmoid, sigmoid_deriv, x, p, xU)..., p
end
@inline function cc_sigmoid(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return sigmoid(x), sigmoid_deriv(x), p)
    (xU <= 0.0) && (return dline_seg(sigmoid, sigmoid_deriv, x, xL, xU)..., p)
    if p == Inf
        p, flag = secant(0.0, xU, 0.0, xU, sigmoid_env, xL, 0.0)
        flag && (p = golden_section(0.0, xU, sigmoid_env, xL, 0.0))
    end
    (x <= p) && (return dline_seg(sigmoid, sigmoid_deriv, x, xL, p)..., p)
    return sigmoid(x), sigmoid_deriv(x), p
end

"""
bisigmoid

The `bisigmoid` activation function `bisigmoid(x) = (1.0 - exp(-x))/(1.0 + exp(-x))`.
"""
@inline bisigmoid(x) = (1.0 - exp(-x))/(1.0 + exp(-x))
@inline bisigmoid(x::Float64) = (1.0 - exp(-x))/(1.0 + exp(-x))
@inline function bisigmoid(x::Interval{Float64})
    xLintv = Interval(x.lo)
    xUintv = Interval(x.hi)
    xLc = (1.0 - exp(-xLintv))/(1.0 + exp(-xLintv))
    xUc = (1.0 - exp(-xUintv))/(1.0 + exp(-xUintv))
    return Interval(xLc.hi, xUc.hi)
end
@inline bisigmoid_deriv(x::Float64) = 2.0*exp(x)/(exp(x) + 1.0)^2
@inline function bisigmoid_deriv2(x::Float64)
    term1 = exp(-x)/(exp(-x) + 1.0)
    term2 = 2.0*exp(-2.0*x)/(exp(-x) + 1.0)^2
    term3 = (1.0 - exp(-x))*(term2 - term1)/(exp(-x) + 1.0)
    return term3 - term1 + term2
end
@inline function bisigmoid_env(x::Float64, y::Float64, z::Float64)
    bisigmoid(y) - bisigmoid(x) - bisigmoid_deriv(x)*(y - x)
end
@inline function cv_bisigmoid(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return dline_seg(bisigmoid, bisigmoid_deriv, x, xL, xU)..., p)
    (xU <= 0.0) && (return bisigmoid(x), bisigmoid_deriv(x), p)
    if p === Inf
        p, flag = secant(xL, 0.0, xL, 0.0, bisigmoid_env, xU, 0.0)
        flag && (p = golden_section(xL, 0.0, bisigmoid_env, xU, 0.0))
    end
    (x <= p) && (return bisigmoid(x), bisigmoid_deriv(x), p)
    return dline_seg(bisigmoid, bisigmoid_deriv, x, p, xU)..., p
end
@inline function cc_bisigmoid(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return bisigmoid(x), bisigmoid_deriv(x), p)
    (xU <= 0.0) && (return dline_seg(bisigmoid, bisigmoid_deriv, x, xL, xU)..., p)
    if p === Inf
        p, flag = secant(0.0, xU, 0.0, xU, bisigmoid_env, xL, 0.0)
        flag && (p = golden_section(0.0, xU, bisigmoid_env, xL, 0.0))
    end
    (x <= p) && (return dline_seg(bisigmoid, bisigmoid_deriv, x, xL, p)..., p)
    return bisigmoid(x), bisigmoid_deriv(x), p
end

#=
softsign

The `softsign` activation function `softsign(x) = x/(1.0 + abs(x))`.
=#
@inline function softsign(x::Interval{Float64})
    xLintv = Interval(x.lo)
    xUintv = Interval(x.hi)
    xLc = xLintv/(1.0 + abs(xLintv))
    xUc = xUintv/(1.0 + abs(xUintv))
    return Interval(xLc.hi, xUc.hi)
end
@inline softsign_deriv(x::Float64) = 1.0/(1.0 + abs(x))^2
@inline softsign_deriv2(x::Float64) = -2.0*x/(abs(x)*(1.0 + abs(x))^(3))
@inline function softsign_env(x::Float64, y::Float64, z::Float64)
    (x - y)*softsign_deriv(x) - (softsign(x) - softsign(y))
end
@inline function cv_softsign(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return dline_seg(softsign, softsign_deriv, x, xL, xU)..., p)
    (xU <= 0.0) && (return softsign(x), softsign_deriv(x), p)
    if p === Inf
        p, flag = secant(xL, 0.0, xL, 0.0, softsign_env, xU, 0.0)
        flag && (p = golden_section(xL, 0.0, softsign_env, xU, 0.0))
    end
    (x <= p) && (return softsign(x), softsign_deriv(x), p)
    return dline_seg(softsign, softsign_deriv, x, p, xU)..., p
end
@inline function cc_softsign(x::Float64, xL::Float64, xU::Float64, p::Float64)
    (xL >= 0.0) && (return softsign(x), softsign_deriv(x), p)
    (xU <= 0.0) && (return dline_seg(softsign, softsign_deriv, x, xL, xU)..., p)
    if p === Inf
        p, flag = secant(0.0, xU, 0.0, xU, softsign_env, xL, 0.0)
        flag && (p = golden_section(0.0, xU, softsign_env, xL, 0.0))
    end
    (x <= p) && (return dline_seg(softsign, softsign_deriv, x, xL, p)..., p)
    return softsign(x), softsign_deriv(x), p
end

const GELU_MIN = -0.751791524693564457457904946779522
const GELU_2D_ROOT1 = -sqrt(2)
const GELU_2D_ROOT2 = sqrt(2)

const SWISH1_MIN = -1.27846454276107379510935873902298015543947748
const SWISH1_2D_ROOT1 = -2.399357280515467667832739697282283888523
const SWISH1_2D_ROOT2 = 2.399357280515467667832739697282283888523

"""
gelu

The Gaussian Error Linear Unit `gelu` activation function `gelu(x) = x/(1.0 + abs(x))`.
"""
@inline function gelu(x::Interval{Float64})
    xLintv = Interval(x.lo)
    xUintv = Interval(x.hi)
    xLc = xLintv*(1.0 + erf(xLintv/sqrt(2)))/2.0
    xUc = xUintv*(1.0 + erf(xUintv/sqrt(2)))/2.0
    if x.hi < GELU_MIN
        xLcv = xUc.lo
        xUcv = xLc.hi
    elseif GELU_MIN < x.lo
        xLcv = xLc.lo
        xUcv = xUc.hi
    else
        xLcv = GELU_MIN
        xUcv = max(xLc.hi, xUc.hi)
    end
    return Interval(xLcv, xUcv)
end
function gelu_deriv(x::Float64)
    0.5*(1.0 + erf(x/sqrt(2))) + (x/sqrt(2*pi))*exp((-x^2)/2.0)
end
function gelu_deriv2(x::Float64)
    (0.797885 - 0.398942*x^2)*exp((-x^2)/2.0)
end
function gelu_env(x::Float64, y::Float64, z::Float64)
    (y - x)*gelu_deriv(x) + gelu(x) - gelu(y)
end
function gelu_denv(x::Float64, y::Float64, z::Float64)
    (y - x)*gelu_deriv2(x)
end

@inline function gelu_envm(x::Float64, y::Float64, z::Float64)
    gelu_deriv(y)*(x - y) - (gelu(x) - gelu(y))
end

@inline function gelu_rt1(x::Float64, y::Float64, z::Float64)
    gelu_deriv(y)*(x - y) + gelu(y) -gelu(x)
end

@inline function gelu_rt1_deriv(x::Float64, y::Float64, z::Float64)
    gelu_deriv(y) - gelu_deriv(x)
end

@inline function cc_gelu(x::Float64, xL::Float64, xU::Float64, p1::Float64, p2::Float64)
    # Single convexity regions
    if xL >= GELU_2D_ROOT2
        return gelu(x), gelu_deriv(x), p1, p2
    elseif xU <= GELU_2D_ROOT1
        return gelu(x), gelu_deriv(x), p1, p2
    elseif (GELU_2D_ROOT1 <= xL) && (xU <= GELU_2D_ROOT2)
        return dline_seg(gelu, gelu_deriv, x, xL, xU)..., p1, p2
    end

    if xL < GELU_2D_ROOT1
        p2, flag = newton(0.0, GELU_MIN, 0.0, gelu_rt1, gelu_rt1_deriv, xL, 0.0)
        flag && (p2 = golden_section(GELU_2D_ROOT1, 0.0, gelu_rt1, xL, 0.0))
        if p2 > xU
            if p1 === Inf
                p1, flag = secant(xL, GELU_2D_ROOT1, xL, GELU_2D_ROOT1, gelu_env, xU, 0.0)
                flag && (p1 = golden_section(xL, GELU_2D_ROOT1, gelu_env, xU, 0.0))
            end
            (x >= p1) && (return dline_seg(gelu, gelu_deriv, x, p1, xU)..., p1, p2)
            return gelu(x), gelu_deriv(x), p1, p2
        else
            return dline_seg(gelu, gelu_deriv, x, xL, xU)..., p1, p2
        end
    end

    if xL > GELU_2D_ROOT1
        p2, flag = newton(0.5*(GELU_MIN + GELU_2D_ROOT2), GELU_MIN, GELU_2D_ROOT2, gelu_rt1, gelu_rt1_deriv, xU, 0.0)
        flag && (p2 = golden_section(GELU_MIN, GELU_2D_ROOT2, gelu_rt1, xU, 0.0))
        if p2 < xL
            if p1 === Inf
                p1, flag = secant(GELU_2D_ROOT2, xU, GELU_2D_ROOT2, xU, gelu_env, xL, 0.0)
                flag && (p1 = golden_section(GELU_2D_ROOT2, xU, gelu_env, xL, 0.0))
            end
            (x >= p1) && (return dline_seg(gelu, gelu_deriv, x, p1, xU)..., p1, p2)
            return gelu(x), gelu_deriv(x), p1, p2
        else
            return dline_seg(gelu, gelu_deriv, x, xL, xU)..., p1, p2
        end
    end
end
@inline function cv_gelu(x::Float64, xL::Float64, xU::Float64, p1::Float64, p2::Float64)

    # Single convexity regions
    if xL >= GELU_2D_ROOT2
        return dline_seg(gelu, gelu_deriv, x, xL, xU)..., p1, p2
    elseif xU <= GELU_2D_ROOT1
        return dline_seg(gelu, gelu_deriv, x, xL, xU)..., p1, p2
    elseif (GELU_2D_ROOT1 <= xL) && (xU <= GELU_2D_ROOT2)
        return gelu(x), gelu_deriv(x), p1, p2
    end

    if xL < GELU_2D_ROOT1
        if p1 === Inf
            p1, flag = newton(0.5*(GELU_2D_ROOT1 + GELU_MIN), GELU_2D_ROOT1, GELU_MIN, gelu_env, gelu_denv, xL, 0.0)
            flag && (p1 = golden_section(GELU_2D_ROOT1, GELU_MIN, gelu_env, xL, 0.0))
        end
    else
        p1 = -Inf
    end

    if xU > GELU_2D_ROOT2
        if p2 === Inf
            p2, flag = newton(0.5*(GELU_MIN + GELU_2D_ROOT2), GELU_MIN, GELU_2D_ROOT2, gelu_env, gelu_denv, xU, 0.0)
            flag && (p2 = golden_section(GELU_MIN, GELU_2D_ROOT2, gelu_env, xU, 0.0))
        end
    else
        p2 = Inf
    end

    if x < p1
        return dline_seg(gelu, gelu_deriv, x, xL, p1)..., p1, p2
    elseif x > p2
        return dline_seg(gelu, gelu_deriv, x, p2, xU)..., p1, p2
    end
    return gelu(x), gelu_deriv(x), p1, p2
end

"""
swish

The Swish-1 activation function `swish(x) = x/(1.0 + exp(-x))`.
"""
@inline function swish(x::Interval{Float64})
    xLintv = Interval(x.lo)
    xUintv = Interval(x.hi)
    xLc = xLintv/(1.0 + exp(-xLintv))
    xUc = xUintv/(1.0 + exp(-xUintv))
    if x.hi < SWISH1_MIN
        xLcv = xUc.lo
        xUcv = xLc.hi
    elseif SWISH1_MIN < x.lo
        xLcv = xLc.lo
        xUcv = xUc.hi
    else
        xmin = Interval(SWISH1_MIN)
        xminI = xmin/(1.0 + exp(-xmin))
        xLcv = xminI.lo
        xUcv = max(xLc.hi, xUc.hi)
    end
    return Interval(xLcv, xUcv)
end
@inline function swish_deriv(x::Float64)
    z = 1.0/(1.0 + exp(-x))
    swish(x) + z*(1-swish(x))
end
@inline function swish_deriv2(x::Float64)
    z = 1.0/(1.0 + exp(-x))
    dz = exp(-x)/(1.0 + exp(-x))^2
    swish_deriv(x) + dz*(1.0 - swish(x)) - z*swish_deriv(x)
end
@inline function swish_env(x::Float64, y::Float64, z::Float64)
    swish_deriv(x)*(y - x) + swish(x) - swish(y)
end
@inline function swish_denv(x::Float64, y::Float64, z::Float64)
    swish_deriv2(x)*(y - x)
end

@inline function swish_envm(x::Float64, y::Float64, z::Float64)
    swish_deriv(y)*(x - y) - (swish(x) - swish(y))
end

@inline function swish_rt1(x::Float64, y::Float64, z::Float64)
    swish_deriv(y)*(x - y) + swish(y) - swish(x)
end
@inline function swish_rt1_deriv(x::Float64, y::Float64, z::Float64)
    swish_deriv(y) - swish_deriv(x)
end

@inline function cc_swish(x::Float64, xL::Float64, xU::Float64, p1::Float64, p2::Float64)
    # Single convexity regions
    if xL >= SWISH1_2D_ROOT2
        return swish(x), swish_deriv(x), p1, p2
    elseif xU <= SWISH1_2D_ROOT1
        return swish(x), swish_deriv(x), p1, p2
    elseif (SWISH1_2D_ROOT1 <= xL) && (xU <= SWISH1_2D_ROOT2)
        return dline_seg(swish, swish_deriv, x, xL, xU)..., p1, p2
    end

    if xL < SWISH1_2D_ROOT1
        p2, flag = newton(0.0, SWISH1_MIN, 0.0, swish_rt1, swish_rt1_deriv, xL, 0.0)
        flag && (p2 = golden_section(SWISH1_2D_ROOT1, 0.0, swish_rt1, xL, 0.0))
        if p2 > xU
            if p1 === Inf
                p1, flag = secant(xL, SWISH1_2D_ROOT1, xL, SWISH1_2D_ROOT1, swish_env, xU, 0.0)
                flag && (p1 = golden_section(xL, SWISH1_2D_ROOT1, swish_env, xU, 0.0))
            end
            (x >= p1) && (return dline_seg(swish, swish_deriv, x, p1, xU)..., p1, p2)
            return swish(x), swish_deriv(x), p1, p2
        else
            return dline_seg(swish, swish_deriv, x, xL, xU)..., p1, p2
        end
    end

    if xL > SWISH1_2D_ROOT1
        p2, flag = newton(0.5*(SWISH1_MIN + SWISH1_2D_ROOT2), SWISH1_MIN, SWISH1_2D_ROOT2, swish_rt1, swish_rt1_deriv, xU, 0.0)
        flag && (p2 = golden_section(SWISH1_MIN, SWISH1_2D_ROOT2, swish_rt1, xU, 0.0))
        if p2 < xL
            if p1 === Inf
                p1, flag = secant(SWISH1_2D_ROOT2, xU, SWISH1_2D_ROOT2, xU, swish_env, xL, 0.0)
                flag && (p1 = golden_section(SWISH1_2D_ROOT2, xU, swish_env, xL, 0.0))
            end
            (x >= p1) && (return dline_seg(swish, swish_deriv, x, p1, xU)..., p1, p2)
            return swish(x), swish_deriv(x), p1, p2
        else
            return dline_seg(swish, swish_deriv, x, xL, xU)..., p1, p2
        end
    end
end
@inline function cv_swish(x::Float64, xL::Float64, xU::Float64, p1::Float64, p2::Float64)

    # Single convexity regions
    if xL >= SWISH1_2D_ROOT2
        return dline_seg(swish, swish_deriv, x, xL, xU)..., p1, p2
    elseif xU <= SWISH1_2D_ROOT1
        return dline_seg(swish, swish_deriv, x, xL, xU)..., p1, p2
    elseif (SWISH1_2D_ROOT1 <= xL) && (xU <= SWISH1_2D_ROOT2)
        return swish(x), swish_deriv(x), p1, p2
    end

    if xL < SWISH1_2D_ROOT1
        if p1 === Inf
            p1, flag = newton(0.5*(SWISH1_2D_ROOT1 + SWISH1_MIN), SWISH1_2D_ROOT1, SWISH1_MIN, swish_env, swish_denv, xL, 0.0)
            flag && (p1 = golden_section(SWISH1_2D_ROOT1, SWISH1_MIN, swish_env, xL, 0.0))
        end
    else
        p1 = -Inf
    end

    if xU > SWISH1_2D_ROOT2
        if p2 === Inf
            p2, flag = newton(0.5*(SWISH1_MIN + SWISH1_2D_ROOT2), SWISH1_MIN, SWISH1_2D_ROOT2, swish_env, swish_denv, xU, 0.0)
            flag && (p2 = golden_section(SWISH1_MIN, SWISH1_2D_ROOT2, swish_env, xU, 0.0))
        end
    else
        p2 = Inf
    end

    if x < p1
        return dline_seg(swish, swish_deriv, x, xL, p1)..., p1, p2
    elseif x > p2
        return dline_seg(swish, swish_deriv, x, p2, xU)..., p1, p2
    end
    return swish(x), swish_deriv(x), p1, p2
end

# define kernel and operator for sigmoid, bisigmoid, softsign, gelu
for expri in (:pentanh, :sigmoid, :bisigmoid, :softsign)
    expri_cv = Symbol("cv_"*String(expri))
    expri_cc = Symbol("cc_"*String(expri))
    expri_kernel = Symbol(String(expri)*"_kernel")
    eps_min = :xL
    eps_max = :xU
    @eval @inline function ($expri_kernel)(x::MC{N, T}, y::Interval{Float64},
                            cv_p::Float64, cc_p::Float64) where {N,T<:Union{NS,MV}}
        xL = x.Intv.lo
        xU = x.Intv.hi
        midcv, cv_id = mid3(x.cc, x.cv, $eps_min)
        midcc, cc_id = mid3(x.cc, x.cv, $eps_max)
        cv, dcv, cv_p = $(expri_cv)(midcv, xL, xU, cv_p)
        cc, dcc, cc_p = $(expri_cc)(midcc, xL, xU, cc_p)
        cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
        cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
        cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad)
        return MC{N, T}(cv, cc, y, cv_grad, cc_grad, x.cnst), cv_p, cc_p
    end
    @eval @inline function ($expri)(x::MC{N,T}) where {N, T<:RelaxTag}
        z, tp1, tp2 = ($expri_kernel)(x, ($expri)(x.Intv), Inf, Inf)
        return z
    end
end

for expri in (:swish, :gelu)
    expri_cv = Symbol("cv_"*String(expri))
    expri_cc = Symbol("cc_"*String(expri))
    expri_kernel = Symbol(String(expri)*"_kernel")
    if expri == swish
        eps_min = :(SWISH1_MIN < xL ? xL : (SWISH1_MIN > xU ? xU : SWISH1_MIN))
        eps_max = :(swish(xL) < swish(xU) ? xU : xL)
    else
        eps_min = :(GELU_MIN > xU ? xU : (GELU_MIN < xL ? xL : GELU_MIN))
        eps_max = :(gelu(xL) < gelu(xU) ? xU : xL)
    end
    @eval @inline function ($expri_kernel)(x::MC{N, T}, y::Interval{Float64},
                            cv_p::Float64, cc_p::Float64, cv_p2::Float64,
                            cc_p2::Float64) where {N,T<:Union{NS,MV}}
        xL = x.Intv.lo
        xU = x.Intv.hi
        midcv, cv_id = mid3(x.cc, x.cv, $eps_min)
        midcc, cc_id = mid3(x.cc, x.cv, $eps_max)
        cv, dcv, cv_p, cv_p2 = $(expri_cv)(midcv, xL, xU, cv_p, cv_p2)
        cc, dcc, cc_p, cc_p2 = $(expri_cc)(midcc, xL, xU, cc_p, cc_p2)
        cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
        cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
        cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad)
        return MC{N, T}(cv, cc, y, cv_grad, cc_grad, x.cnst), cv_p, cc_p, cv_p2, cc_p2
    end
    @eval @inline function ($expri)(x::MC{N,T}) where {N, T<:RelaxTag}
        z, tp1, tp2, tp3, tp4 = ($expri_kernel)(x, ($expri)(x.Intv), Inf, Inf, Inf, Inf)
        return z
    end
end

"""
logcosh

The function `logcosh` is defined as `logcosh(x) = log(cosh(x))`.
"""
logcosh_deriv(x::Float64) = tanh(x)
logcosh_deriv2(x::Float64) = sech(x)^2
cc_logcosh(x::Float64, xL::Float64, xU::Float64) = dline_seg(logcosh, logcosh_deriv, x, xL, xU)
cv_logcosh(x::Float64, xL::Float64, xU::Float64) = logcosh(x), logcosh_deriv(x)
function logcosh(x::Interval{Float64})
    (abs(x.lo) < abs(x.hi)) ? (logcosh_xU = log(cosh(Interval(x.hi)))) : (logcosh_xU = log(cosh(Interval(x.lo))))
    logcosh_xL = Interval(0.0)
    (x.lo >= 0.0) && (logcosh_xL = log(cosh(Interval(x.lo))))
    (x.hi <= 0.0) && (logcosh_xL = log(cosh(Interval(x.hi))))
    Interval{Float64}(logcosh_xL.lo, logcosh_xU.hi)
end
@inline function logcosh_kernel(x::MC{N, T}, y::Interval{Float64}) where {N,T<:Union{NS,MV}}
    xL = x.Intv.lo
    xU = x.Intv.hi
    eps_min = in(0.0, x) ? 0.0 : (xU <= 0.0 ? xU : xL)
    eps_max = (abs(xL) < abs(xU)) ? xU : xL
    midcv, cv_id = mid3(x.cc, x.cv, eps_min)
    midcc, cc_id = mid3(x.cc, x.cv, eps_max)
    cv, dcv = cv_logcosh(midcv, xL, xU)
    cc, dcc = cc_logcosh(midcc, xL, xU)
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv, cc, cv_grad, cc_grad = cut(y.lo, y.hi, cv, cc, cv_grad, cc_grad)
    return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end
@inline function logcosh_kernel(x::MC{N, T}, y::Interval{Float64}) where {N,T<:Diff}
    xL = x.Intv.lo
    xU = x.Intv.hi
    eps_min = in(0.0, x) ? 0.0 : (xU <= 0.0 ? xU : xL)
    eps_max = (abs(xL) < abs(xU)) ? xU : xL
    midcv, cv_id = mid3(x.cc, x.cv, eps_min)
    midcc, cc_id = mid3(x.cc, x.cv, eps_max)
    cv, dcv = cv_logcosh(midcv, xL, xU)
    cc, dcc = cc_logcosh(midcc, xL, xU)
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    return MC{N,T}(cv, cc, y, cv_grad, cc_grad, x.cnst)
end
function logcosh(x::MC{N,T}) where {N, T <: RelaxTag}
	logcosh_kernel(x, logcosh(x.Intv))
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function Annotate(x::Real, y::Real, pix_offset::ImVec2, fmt::String...)

    ccall((:ImPlot_AnnotateStr, libcimplot), Cvoid, (Cdouble,Cdouble,ImVec2,Cstring...),
          x, y, pix_offset, fmt...)
end

function AnnotateClamped(x::Real, y::Real, pix_offset::ImVec2, fmt::String...)

    ccall((:ImPlot_AnnotateClampedStr, libcimplot), Cvoid, (Cdouble,Cdouble,ImVec2,Cstring...),
          x, y, pix_offset, fmt...)
end

function Annotate(x::Real, y::Real, pix_offset::ImVec2, color::ImVec4, fmt::String...)

    ccall((:ImPlot_AnnotateVec4, libcimplot), Cvoid, (Cdouble,Cdouble,ImVec2,ImVec4,Cstring...),
          x, y, pix_offset, color, fmt...)
end

function AnnotateClamped(x::Real, y::Real, pix_offset::ImVec2, color::ImVec4, fmt::String...)

    ccall((:ImPlot_AnnotateClampedVec4, libcimplot), Cvoid, (Cdouble,Cdouble,ImVec2,ImVec4,Cstring...),
          x, y, pix_offset, color, fmt...)
end

BeginLegendDragDropSource(label_id) = BeginLegendDragDropSource(label_id, 0)
BeginLegendPopup(label_id) = BeginLegendPopup(label_id, 1)

export ImPlotPoint, ImPlotRange, ImPlotLimits, ImPlotStyle, ImPlotInputMap, ImPlotContext
export CreateContext, DestroyContext, GetCurrentContext, SetCurrentContext, SetImGuiContext
export EndPlot
export SetNextPlotLimits, SetNextPlotLimitsX, SetNextPlotLimitsY, LinkNextPlotLimits, FitNextPlotAxes
export SetPlotYAxis, HideNextItem, IsPlotHovered, IsPlotXAxisHovered, IsPlotYAxisHovered, IsPlotQueried
export SetLegendLocation, SetMousePosLocation, IsLegendEntryHovered
export BeginLegendDragDropSource, EndLegendDragDropSource, BeginLegendPopup, EndLegendPopup
export GetStyle
export GetStyleColorName, GetMarkerName, PopColormap, GetColormapSize, ShowColormapScale
export GetColormapName, PushPlotClipRect, PopPlotClipRect
export ShowStyleSelector, ShowColormapSelector, ShowUserGuide, ShowDemoWindow, ShowMetricsWindow
export Annotate, AnnotateClamped, PushPlotClipRect, PopPlotClipRect, GetPlotDrawList


text:
server = LanguageServerInstance(IOBuffer(), IOBuffer(), dirname(Pkg.Types.Context().env.project_file), first(Base.DEPOT_PATH))
server.runlinter = true
server.jr_endpoint = nothing

LanguageServer.initialize_request(init_request, server, nothing)

testtext = """
decode_overlong
"""
LanguageServer.textDocument_didOpen_notification(LanguageServer.DidOpenTextDocumentParams(LanguageServer.TextDocumentItem("testdoc", "julia", 0, testtext)), server, nothing)

doc = LanguageServer.getdocument(server, LanguageServer.URI2("testdoc"))

@test !isempty(LanguageServer.textDocument_codeAction_request(LanguageServer.CodeActionParams(LanguageServer.TextDocumentIdentifier("testdoc"), LanguageServer.Range(LanguageServer.Position(0, 1), LanguageServer.Position(0, 1)), LanguageServer.CodeActionContext(LanguageServer.Diagnostic[doc.diagnostics[1]], missing)), server, nothing))


text:
"""
Google Cloud Storage API
"""
module _storage

export storage

using ..api
using ...root

"""
Google Cloud Storage API root.
"""
storage = APIRoot(
    "$API_ROOT/storage/v1",
    Dict(
        "devstorage.full_control" => "Read/write and ACL management access to Google Cloud Storage",
        "devstorage.read_write" => "Read/write access to Google Cloud Storage",
        "devstorage.read_only" => "Read-only access to Google Cloud Storage",
    );
    Bucket=APIResource("b";
        delete=APIMethod(:DELETE, "{bucket}", "Permanently deletes an empty bucket."),
        get=APIMethod(:GET, "{bucket}", "Returns metadata for the specified bucket."),
        insert=APIMethod(:POST, "", "Creates a new bucket.",
            Dict(:project => :project_id)
        ),
        list=APIMethod(:GET, "", "Retrieves a list of buckets for a given project.",
            Dict(:project => :project_id);
            transform=(x, t) -> map(t, get(x, :items, []))
        ),
        patch=APIMethod(:PATCH, "{bucket}", "Updates a bucket. This method supports patch semantics."),
        update=APIMethod(:PUT, "{bucket}", "Updates a bucket."),
    ),
    BucketAccessControl=APIResource("b/{bucket}/acl";
        delete=APIMethod(:DELETE, "{entity}", "Permanently deletes the ACL entry for the specified entity on the specified bucket."),
        get=APIMethod(:GET, "{entity}", "Returns the ACL entry for the specified entity on the specified bucket."),
        insert=APIMethod(:POST, "", "Creates a new ACL entry on the specified bucket."),
        list=APIMethod(:GET, "", "Retrieves ACL entries on a specified bucket."; transform=(x, t) -> map(t, get(x, :items, []))),
        patch=APIMethod(:PATCH, "{entity}", "Updates an ACL entry on the specified bucket. This method supports patch semantics."),
        update=APIMethod(:PUT, "{entity}", "Updates an ACL entry on the specified bucket."),
    ),
    Object=APIResource("b/{bucket}/o";
        compose=APIMethod(:POST, "{destinationObject}/compose", "Concatenates a list of existing objects into a new object in the same bucket."),
        copy=APIMethod(:POST, "{sourceObject}/copyTo/b/{destinationBucket}/o/{destinationObject}", "Copies a source object to a destination object. Optionally overrides metadata."),
        delete=APIMethod(:DELETE, "{object}", "Deletes an object and its metadata. Deletions are permanent if versioning is not enabled for the bucket, or if the generation parameter is used."),
        get=APIMethod(:GET, "{object}", "Retrieves an object or its metadata.",
            Dict(:alt => "media")
        ),
        insert=APIMethod(:POST, "$API_ROOT/upload/storage/v1/b/{bucket}/o", "Stores a new object and metadata.",
            Dict(:uploadType => "media")
        ),
        list=APIMethod(:GET, "", "Retrieves a list of objects matching the criteria."; transform=(x, t) -> map(t, get(x, :items, []))),
        patch=APIMethod(:PATCH, "{object}", "Updates a data blob's associated metadata. This method supports patch semantics."),
        rewrite=APIMethod(:POST, "{sourceObject}/rewriteTo/b/{destinationBucket}/o/{destinationObject}", "Rewrites a source object to a destination object. Optionally overrides metadata."),
        update=APIMethod(:PUT, "{object}", "Updates an object's metadata."),
        watchAll=APIMethod(:POST, "watch", "Watch for changes on all objects in a bucket."),
    ),
    ObjectAccessControl=APIResource("b/{bucket}/o/{object}/acl";
        delete=APIMethod(:DELETE, "{entity}", "Permanently deletes the ACL entry for the specified entity on the specified object."),
        get=APIMethod(:GET, "{entity}", "Returns the ACL entry for the specified entity on the specified object."),
        insert=APIMethod(:POST, "", "Creates a new ACL entry on the specified object."),
        list=APIMethod(:GET, "", "Retrieves ACL entries on the specified object."; transform=(x, t) -> map(t, get(x, :items, []))),
        patch=APIMethod(:PATCH, "{entity}", "Updates an ACL entry on the specified object. This method supports patch semantics."),
        update=APIMethod(:PUT, "{entity}", "Updates an ACL entry on the specified object."),
    ),
    DefaultObjectAccessControl=APIResource("b/{bucket}/defaultObjectAcl";
        delete=APIMethod(:DELETE, "entity", "Permanently deletes the default object ACL entry for the specified entity on the specified bucket."),
        get=APIMethod(:GET, "entity", "Returns the default object ACL entry for the specified entity on the specified bucket."),
        insert=APIMethod(:POST, "", "Creates a new default object ACL entry on the specified bucket."),
        list=APIMethod(:GET, "", "Retrieves default object ACL entries on the specified bucket."; transform=(x, t) -> map(t, get(x, :items, []))),
        patch=APIMethod(:PATCH, "entity", "Updates a default object ACL entry on the specified bucket. This method supports patch semantics."),
        update=APIMethod(:PUT, "entity", "Updates a default object ACL entry on the specified bucket."),
    ),
    Channel=APIResource("channels";
        stop=APIMethod(:POST, "stop", "Stop receiving object change notifications through this channel."),
    ),
)

end


text:
using Test
using CGP
CGP.Config.init("cfg/test.yaml")

CTYPES = [CGPChromo, PCGPChromo]

function test_crossover(p1::Chromosome, p2::Chromosome, child::Chromosome)
      @test child != p1
      @test child != p2
      @test child.genes != p1.genes
      @test child.genes != p2.genes
      @test child.nin == p1.nin
      @test child.nin == p2.nin
      @test child.nout == p1.nout
      @test child.nout == p2.nout
      @test distance(p1, child) > 0
      @test distance(p2, child) > 0
      @test length(intersect([n.p for n in p1.nodes], [n.p for n in child.nodes])) > 0
      @test length(intersect([n.p for n in p2.nodes], [n.p for n in child.nodes])) > 0
      @test length(intersect([n.f for n in p1.nodes], [n.f for n in child.nodes])) > 0
      @test length(intersect([n.f for n in p2.nodes], [n.f for n in child.nodes])) > 0
      # @test length(intersect(forward_connections(p1), forward_connections(child))) > 0
      # @test length(intersect(forward_connections(p2), forward_connections(child))) > 0
      # @test length(intersect(get_output_trace(p1), get_output_trace(child))) > 0
      # @test length(intersect(get_output_trace(p2), get_output_trace(child))) > 0
end

@testset "Crossover tests" begin
    for ct in CTYPES
        println(ct)
        nin = rand(1:100); nout = rand(1:100)
        p1 = ct(nin, nout)
        p2 = ct(nin, nout)
        @test p1.genes != p2.genes
        @testset "Single point $ct" begin
            child = single_point_crossover(p1, p2)
            test_crossover(p1, p2, child)
        end
        @testset "Random node $ct" begin
            child = random_node_crossover(p1, p2)
            test_crossover(p1, p2, child)
        end
        @testset "Aligned node $ct" begin
            child = aligned_node_crossover(p1, p2)
            test_crossover(p1, p2, child)
        end
        @testset "Proportional $ct" begin
            child = proportional_crossover(p1, p2)
            test_crossover(p1, p2, child)
        end
        @testset "Output graph $ct" begin
            child = output_graph_crossover(p1, p2)
            test_crossover(p1, p2, child)
        end
        @testset "Subgraph $ct" begin
            child = subgraph_crossover(p1, p2)
            test_crossover(p1, p2, child)
        end
    end
end

nothing


text:
## include("lptest.jl")
## lptest("problem-00003-00003-00000a.txt")

# using Printf
using JuMP
using Clp

function lptest(filename)
    lines = open(filename, "r") do fp
        readlines(fp)
    end

    m = parse(Int, lines[1])
    n = parse(Int, lines[2])

    b = parse.(Float64, (split(lines[3],[','])))
    c = parse.(Float64, (split(lines[4],[','])))
    A = zeros(Float64, m, n)
    for i=1:m
        A[i,:] = parse.(Float64, (split(lines[4+i],[','])))
    end

    # @show b
    # @show c
    # @show A

    before_solve = time()
    model1 = Model(Clp.Optimizer)
    @variable(model1, x[1:n]>=0)
    @objective(model1, Max, dot(c, x))
    @constraint(model1, A*x .<= b)
    optimize!(model1)
    time_solve = time() - before_solve
    status = termination_status(model1)
    x_value = value.(x)

    @show x_value
    @show status
    @show time_solve
end

max_stars_repo_name: makoto-yamashita/lp-generator

max_issues_repo_name: makoto-yamashita/lp-generator

max_forks_repo_name: makoto-yamashita/lp-generator


text:
# this file tests various aspects of the advection routines
using MagmaThermoKinematics
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
using Plots
using LinearAlgebra
using SpecialFunctions
using Test

const CreatePlots = false      # easy way to deactivate plotting throughout

# Initialize for multiple threads (GPU is not tested here)
@init_parallel_stencil(Threads, Float64, 3);    # initialize parallel stencil in 3D

function test_Interpolation(Dimension="2D", InterpolationMethod="Linear")
  # test interpolation methods in 2D and 3D

  if Dimension=="2D"
    # Model parameters
    W,H                     =   1.,  1.;                                    # Width, Length, Height

    # Define coarse grid
    Nx, Nz                  =   33, 33;                                     # resolution of coarse grid
    dx,dz                   =   W/(Nx-1), H/(Nz-1);                         # grid size [m]
    x,z                     =   0:dx:((Nx-1)*dx), 0:dz:((Nz-1)*dz);         # 1D coordinate arrays
    coords                  =   collect(Iterators.product(x,z))             # generate coordinates from 1D coordinate vectors
    X,Z                     =   (x->x[1]).(coords), (x->x[2]).(coords);     # transfer coords to 3D arrays
    Grid, Spacing           =   (x,z), (dx,dz);

    # Define fine grid
    Nx_f, Nz_f              =   65, 65;                                     # resolution of fine grid
    dx_f,dz_f               =   W/(Nx_f-1), H/(Nz_f-1);                     # grid size [m]
    x_f,z_f                 =   0:dx_f:H, 0:dz_f:H;                         # 1D coordinate arrays
    coords_f                =   collect(Iterators.product(x_f,z_f))         # generate coordinates from 1D coordinate vectors
    X_f,Z_f                 =   (x->x[1]).(coords_f), (x->x[2]).(coords_f); # transfer coords to 3D arrays
    Grid_f                  =   (x_f, z_f);     # note
    Y_f = [];
    # Define function on coarse grid
    T                       =   cos.(pi.*X).*sin.(2*pi.*Z)

    Data_coarse1            =   tuple(T);         # 1 dataset on grid
    Data_coarse2            =   (T, T);           # 2 datasets on grid

    Grid_fine               =   (X_f, Z_f)

    Data_fine1              =   (X_f*0,);
    Data_fine2              =   (X_f.*0, Z_f.*0);

  elseif Dimension=="3D"
      # Model parameters
      W,L,H                 =   1., 1., 1.;                                    # Width, Length, Height

      # Define coarse grid
      Nx, Ny, Nz              =   33,33,33;                                                    # resolution of coarse grid

      Nx, Ny, Nz              =   150,150,150;                                                  # resolution of coarse grid

      dx,dy,dz                =   W/(Nx-1), L/(Ny-1), H/(Nz-1);                                 # grid size [m]
      x,y,z                   =   0:dx:((Nx-1)*dx),  0:dy:((Ny-1)*dy), 0:dz:((Nz-1)*dz);        # 1D coordinate arrays
      coords                  =   collect(Iterators.product(x,y,z))                             # generate coordinates from 1D coordinate vectors
      X,Y,Z                   =   (x->x[1]).(coords), (x->x[2]).(coords), (x->x[3]).(coords);   # transfer coords to 3D arrays
      Grid, Spacing           =   (x,y,z), (dx,dy,dz);

      # Define fine grid
      Nx_f, Ny_f, Nz_f        =   65, 65, 65;                                                       # resolution of fine grid
      Nx_f, Ny_f, Nz_f        =   165, 165, 165;                                                       # resolution of fine grid
      dx_f,dy_f,dz_f          =   W/(Nx_f-1), L/(Ny_f-1), H/(Nz_f-1);                               # grid size [m]
      x_f,y_f,z_f             =   0:dx_f:H, 0:dy_f:L,0:dz_f:H;                                      # 1D coordinate arrays
      coords_f                =   collect(Iterators.product(x_f,y_f, z_f))                          # generate coordinates from 1D coordinate vectors

      X_f,Y_f, Z_f            =   (x->x[1]).(coords_f), (x->x[2]).(coords_f), (x->x[3]).(coords_f); # transfer coords to 3D arrays
      Grid_f                  =   (x_f, y_f, z_f);

      # Define function on coarse grid
      T                       =   cos.(pi.*X).*sin.(2*pi.*Z).*sin.(2*pi.*Y)

      Data_coarse1            =   tuple(T);         # 1 dataset on grid
      Data_coarse2            =   (T, T, T);        # 3 datasets on grid

      Grid_fine               =   (X_f, Y_f, Z_f)
      Data_fine1              =   (X_f.*0,);
      Data_fine2              =   (X_f.*0, Y_f.*0, Z_f.*0);

  end

  # Perform interpolation from coarse-fine grid
  Interpolate!(Data_fine1, Grid, Data_coarse1, Grid_fine,   InterpolationMethod);       # one field
  Interpolate!(Data_fine2, Grid, Data_coarse2, Grid_fine,   InterpolationMethod);       # several fields

  # Compute error
  if Dimension=="2D"
    Tanal       =   cos.(pi.*X_f).*sin.(2*pi.*Z_f);
    Terror      =   Tanal - Data_fine2[2];

    if CreatePlots
      p1          =   contourf(x, z,      Data_coarse1[1]',       aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Coarse ",  dpi=150, levels=10)
      p2          =   contourf(x_f, z_f,  Data_fine2[2]',         aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Fine ",    dpi=150, levels=10)
      p3          =   contourf(x_f, z_f,  Terror',                aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="error ",    dpi=150, levels=10)

      plot(p1,p2,p3);

      png("Interpolation_2D_$InterpolationMethod")
    end

  elseif Dimension=="3D"
    Tanal       =   cos.(pi.*X_f).*sin.(2*pi.*Z_f).*sin.(2*pi.*Y_f);
    Terror      =   Tanal - Data_fine2[2];

    if CreatePlots
      p1          =   contourf(x, z,      Data_coarse1[1][:,10,:]', aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Coarse ",  dpi=150, levels=10)
      p2          =   contourf(x_f, z_f,  Data_fine2[2][:,10,:]',   aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Fine ",    dpi=150, levels=10)
      p3          =   contourf(x_f, z_f,  Terror[:,10,:]',          aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="error ",    dpi=150, levels=10)

      plot(p1,p2,p3);
      png("Interpolation_3D_$InterpolationMethod")
    end

  end

  error = norm(Terror[:],2)/length(Terror[:]);
  return error;        # return error
end

function test_SemiLagrangian2D(Method="ConstantZ",  InterpolationMethod="Linear", AdvectionMethod="RK2")
  # SemiLagrangian advection test in 2D

  # Model parameters
  W,H                     =   1.,  1.;                      # Width, Length, Height
  Nx, Nz                  =   129, 129;                      # resolution
  σ                       =   0.1;                          # halfwidth of gaussian
  Tmax                    =   2;                            # max. of gaussian
  dx,dz                   =   W/(Nx-1), H/(Nz-1);           # grid size [m]

  # Array initializations (1 - main arrays on which we can initialize properties)
  T                       =   @zeros(Nx,Nz);
  Tnew                    =   @zeros(Nx,Nz);

  # Set up model geometry & initial T structure
  x,z                     =   0:dx:((Nx-1)*dx), 0:dz:((Nz-1)*dz);
  coords                  =   collect(Iterators.product(x,z))             # generate coordinates from 1D coordinate vectors
  X,Z                     =   (x->x[1]).(coords), (x->x[2]).(coords);     # transfer coords to 2D arrays
  Grid,GridFull,Spacing   =   (x,z), (X,Z), (dx,dz);

  Xc                      =   0.5;
  Zc                      =   0.75;
  if      Method=="ConstantZ"
    # Simple advection in z-direction
    Vx        =   X.*0;
    Vz        =   Z.*0 .- 1;
    TotalTime =   0.5;                         # thermal cooling age

  elseif  Method=="Rotation"
    # Rigid body rotation:
    Vx        =    (Z .- 0.5);
    Vz        =   -(X .- 0.5);
    TotalTime =   2*pi;

  elseif  Method=="Shear"
    # Shear velocity field:
    Vx        =  -sin.(pi.*X).*cos.(pi.*Z);
    Vz        =   cos.(pi.*X).*sin.(pi.*Z);

  end

  T                      .=   Tmax.*exp.(  -( (X .- Xc).^2 .+  (Z .- Zc).^2 )./σ^2);                 # initial gaussian profile
  Tnew                   .=   T;

  maxVel      =   max( maximum(abs.(Vx)), maximum(abs.(Vz)) );
  dt          =   min(dx,dz)/2.0/maxVel;   # stable timestep (to ensure we move no more than 0.5*(dx|dz) per timesstep )
  numTime     =   ceil(TotalTime/dt);
  nt          =   Int(numTime);
  dt          =   TotalTime/nt;

  #ENV["GKSwstype"]="nul"; if isdir("viz2D_out")==false mkdir("viz2D_out") end; loadpath = "./viz2D_out/"; anim = Animation(loadpath,String[])
  #println("Animation directory: $(anim.dir)")

  time,time_kyrs          = 0.0, 0.0;
  err = 100;
  it = 0;
  #nt = 100
  for it=1:nt

      # Perform an advection step for temperature
      Tnew = AdvectTemperature(T,    Grid,  GridFull, (Vx,Vz),    dt, AdvectionMethod, InterpolationMethod);

      T, Tnew     =   Tnew, T;                                                                # Update temperature
      time        =   time + dt;                                            # Keep track of evolved time

      if mod(it,1000)==0  # print progress
          println(" Timestep $it = $((time)) ")
          #p1          =   heatmap(x_km, z_km, T[:,Int(Ny/2),:]',         aspect_ratio=1, xlims=(x_km[1],x_km[end]), ylims=(z_km[1],z_km[end]),   c=:inferno, title="Temperature, $(round(time_kyrs, digits=2)) kyrs",  dpi=150)
          p1          =   contourf(x, z, T[:,:]',         aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Temperature, $(round(time, digits=2)) ",  dpi=150, levels=10)

          plot(p1); frame(anim)
      end

  end

  x_km, z_km  =   x, z;

  # compute analytical solution
  if      Method=="ConstantZ"
    # Simple advection in z-direction
    Xc    =   0.5;
    Zc    =   0.25;

  elseif      Method=="Rotation"
      # Should arrive @ same point
  end
  Tanal =   Tmax.*exp.(  -( (X .- Xc).^2 .+  (Z .- Zc).^2 )./σ^2);                 # initial gaussian profile
  fname       = "Advection_2D_$(Method)_$(AdvectionMethod)_$(InterpolationMethod)";

  Terror      =  T - Tanal;   # error

  if CreatePlots
    # create plot
    p1          =   contourf(x, z, T',      aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Tnumeric, $(round(time, digits=2)) ",  dpi=300, levels=10)
    p2          =   contourf(x, z, Tanal',  aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Tanal, $(round(time, digits=2)) ",     dpi=300, levels=10)
    p3          =   contourf(x, z, Terror', aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Terror, $(round(time, digits=2)) ",    dpi=300, levels=10)
    plot(p1,p2,p3);
    png(fname)
  end

  error = norm(Terror[:],2)/length(Terror[:]);
  return error;        # return error
end

function test_SemiLagrangian3D(Method="ConstantZ",  InterpolationMethod="Linear", AdvectionMethod="RK2")
    # SemiLagrangian advection test in 3D

    # Model parameters
    W,L,H                   =   1., 1., 1.;                         # Width, Length, Height
    Nx, Ny, Nz              =   50, 50, 50;                         # resolution
    σ                       =   0.1;                                # halfwidth of gaussian
    Tmax                    =   2;                                  # max. of gaussian
    TotalTime               =   1;                                  # thermal cooling age
    dx,dy,dz                =   W/(Nx-1), L/(Ny-1), H/(Nz-1);       # grid size [m]
    dt                      =   min(dx,dy,dz)./2;                   # stable timestep (required for explicit FD)

    numTime                 =   ceil(TotalTime/dt);
    nt                      =   Int(numTime);

    # Array initializations (1 - main arrays on which we can initialize properties)
    T                       =   @zeros(Nx,Ny,Nz);
    Tnew                    =   @zeros(Nx,Ny,Nz);

    # Set up model geometry & initial T structure
    x,y,z                   =   0:dx:((Nx-1)*dx), 0:dy:((Ny-1)*dy), 0:dz:((Nz-1)*dz);
    coords                  =   collect(Iterators.product(x,y,z))                               # generate coordinates from 1D coordinate vectors
    X,Y,Z                   =   (x->x[1]).(coords), (x->x[2]).(coords), (x->x[3]).(coords);     # transfer coords to 3D arrays
    Grid, GridFull, Spacing =   (x,y,z), (X,Y,Z), (dx,dy,dz);
    Vx                      =  -sin.(pi.*X).*cos.(pi.*Z);
    Vy                      =   zeros(size(X));
    Vz                      =   cos.(pi.*X).*sin.(pi.*Z);

    Xc                      =   0.5;
    Yc                      =   0.5;
    Zc                      =   0.75;
    if      Method=="ConstantZ"
      # Simple advection in z-direction
      Vx        =   X.*0;
      Vy        =   Y.*0;
      Vz        =   Z.*0 .- 1;
      TotalTime =   0.5;                         # thermal cooling age

    elseif  Method=="Rotation_alongY"
      # Rigid body rotation:
      Vx        =    (Z .- 0.5);
      Vy        =   X.*0;
      Vz        =   -(X .- 0.5);
      TotalTime =   2*pi;

    elseif  Method=="Rotation_alongX"
      # Rigid body rotation:
      Vy        =    (Z .- 0.5);
      Vx        =   Y.*0;
      Vz        =   -(Y .- 0.5);
      TotalTime =   2*pi;

    elseif  Method=="Rotation_alongZ"
      # Rigid body rotation:
      Vy        =   -(X .- 0.5);
      Vx        =    (Y .- 0.5);
      Vz        =   Z.*0.0;
      TotalTime =   2*pi;
      Xc                      =   0.5;
      Yc                      =   0.75;
      Zc                      =   0.5;
   end

    T           .=   Tmax.*exp.(  -(((X .- Xc).^2 .+ (Y .- Yc).^2 .+ (Z .- Zc).^2)./(σ^2)) );                 # initial gaussian profile
    Tnew        .=   T;

    maxVel      =   max( maximum(abs.(Vx)), maximum(abs.(Vy)), maximum(abs.(Vz)) );
    dt          =   min(dx,dy,dz)/2.0/maxVel;   # stable timestep (to ensure we move no more than 0.5*(dx|dz) per timesstep )
    numTime     =   ceil(TotalTime/dt);
    nt          =   Int(numTime);
    dt          =   TotalTime/nt;

    #ENV["GKSwstype"]="nul"; if isdir("viz2D_out")==false mkdir("viz2D_out") end; loadpath = "./viz2D_out/"; anim = Animation(loadpath,String[])
    #println("Animation directory: $(anim.dir)")

    time          = 0.0;
    for it=1:nt

        # Perform an advection step for temperature

        Tnew = AdvectTemperature(T,    Grid,  GridFull, (Vx,Vy,Vz),    dt, AdvectionMethod, InterpolationMethod);

        T, Tnew         =   Tnew, T;                                                                # Update temperature
        time            =   time + dt;                                            # Keep track of evolved time

        if mod(it,1000)==0  # print progress
            println(" Timestep $it = $((time)) ")
            #p1          =   heatmap(x_km, z_km, T[:,Int(Ny/2),:]',         aspect_ratio=1, xlims=(x_km[1],x_km[end]), ylims=(z_km[1],z_km[end]),   c=:inferno, title="Temperature, $(round(time_kyrs, digits=2)) kyrs",  dpi=150)
            #p1          =   heatmap(x, z, T[:,Int(Ny/2),:]',         aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Temperature, $time",  dpi=150)
         #   p1          =   heatmap(x, y, T[:,:,Int(Nz/2)]',         aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Temperature, $time",  dpi=150)

           # p1          =   heatmap(x, z, T[Int(Nx/2),:,:]',         aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Temperature, $time",  dpi=150)

          #  plot(p1); frame(anim)
        end

    end

    x_km, z_km  =   x, z;

    # compute analytical solution
    if      Method=="ConstantZ"
      # Simple advection in z-direction
      Xc    =   0.5;
      Yc    =   0.5;
      Zc    =   0.25;

    elseif  (Method=="Rotation_alongY") || (Method=="Rotation_alongX")
      # Should arrive @ same point

    end
    Tanal =   Tmax.*exp.(  -( (X .- Xc).^2 .+  (Y .- Yc).^2 .+  (Z .- Zc).^2 )./σ^2);                 # gaussian profile
    fname = "Advection_3D_$(Method)_$(AdvectionMethod)_$(InterpolationMethod)";

    Terror      =  T - Tanal;   # error

    Tslice  = T[:,Int(Ny/2),:];
    Tanal1  = Tanal[:,Int(Ny/2),:];
    Terror1 = Terror[:,Int(Ny/2),:];

    Vx1 = Vx[:,Int(Ny/2),:];
    Vy1 = Vy[:,Int(Ny/2),:];
    Vz1 = Vz[:,Int(Ny/2),:];
    X1  = X[:,Int(Ny/2),:];
    Y1  = Y[:,Int(Ny/2),:];
    Z1  = Z[:,Int(Ny/2),:];

    if CreatePlots
      # create plot
      p1          =   contourf(x, z, Tslice',      aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Tnumeric, $(round(time, digits=2)) ",  dpi=300, levels=10)
      p2          =   contourf(x, z, Tanal1',  aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Tanal 3D, $(round(time, digits=2)) ",     dpi=300, levels=10)
      p3          =   contourf(x, z, Terror1', aspect_ratio=1, xlims=(x[1],x[end]), ylims=(z[1],z[end]),   c=:inferno, title="Terror, $(round(time, digits=2)) ",    dpi=300, levels=10)
      plot(p1,p2,p3);

      # st=20;
      # quiver!(X1[1:st:end], Z1[1:st:end], gradient=(Vx1[1:st:end]*0.1,Vz1[1:st:end]*0.1), arrow = :arrow, color = :white)

      png(fname)
    end

    error = norm(Terror[:],2)/length(Terror[:]);
    return error;        # return error
end

function test_AdvectTracers2D(Method="ConstantZ",  AdvectionMethod="RK2")
  # Test tracer advection in 2D

  # Model parameters
  W,H                     =   1.,  1.;                      # Width, Length, Height
  Nx, Nz                  =   129, 129;                     # resolution
  dx,dz                   =   W/(Nx-1), H/(Nz-1);           # grid size [m]

  # Set up model geometry
  x,z                     =   0:dx:((Nx-1)*dx), 0:dz:((Nz-1)*dz);
  coords                  =   collect(Iterators.product(x,z))             # generate coordinates from 1D coordinate vectors
  X,Z                     =   (x->x[1]).(coords), (x->x[2]).(coords);     # transfer coords to 2D arrays
  Grid, Spacing           =   (x,z), (dx,dz);

  # generate tracer coordinates
  Nxt, Nzt                =   50, 50          # num of tracers in x,z
  Wt,Ht                   =   0.1, 0.1;
  dxt,dzt                 =   Wt/(Nxt-1), Ht/(Nzt-1);
  xt,zt                   =   -Wt/2:dxt:Wt/2, -Ht/2:dz:Ht/2;
  coordst                 =   collect(Iterators.product(xt,zt))             # generate coordinates from 1D coordinate vectors
  Xt,Zt                   =   (x->x[1]).(coordst), (x->x[2]).(coordst);     # transfer coords to 2D arrays

  Xc                      =   0.5;
  Zc                      =   0.75;
  if      Method=="ConstantZ"
    # Simple advection in z-direction
    Vx        =   X.*0;
    Vz        =   Z.*0 .- 1;
    TotalTime =   0.5;                         # thermal cooling age

  elseif  Method=="Rotation"
    # Rigid body rotation:
    Vx        =    (Z .- 0.5);
    Vz        =   -(X .- 0.5);
    TotalTime =   2*pi;

  elseif  Method=="Shear"
    # Shear velocity field:
    Vx        =  -sin.(pi.*X).*cos.(pi.*Z);
    Vz        =   cos.(pi.*X).*sin.(pi.*Z);

  end

  maxVel      =   max( maximum(abs.(Vx)), maximum(abs.(Vz)) );
  dt          =   min(dx,dz)/2.0/maxVel;   # stable timestep (to ensure we move no more than 0.5*(dx|dz) per timesstep )
  numTime     =   ceil(TotalTime/dt);
  nt          =   Int(numTime);
  dt          =   TotalTime/nt;

  # Create tracer structure
  Xt .=  Xt .+ Xc; Zt .=  Zt .+ Zc;   # shift coordinates of tracers
  Xt0, Zt0    =   Xt, Zt;

  new_tracer  =   Tracer(num=1, coord=[Xt[1]; Zt[1]], T=0.);            # Create new tracer
  Tracers     =   StructArray([new_tracer]);                            # Create tracer array
  for i=firstindex(Xt)+1:lastindex(Xt)
    new_tracer  =   Tracer(num=i, coord=[Xt[i]; Zt[i]], T=0.);          # Create new tracer
    push!(Tracers, new_tracer);                             # Add new point to existing array
  end
  Tracers0 = copy(Tracers);         # create a copy of the original tracers

  #  ENV["GKSwstype"]="nul"; if isdir("viz2D_out")==false mkdir("viz2D_out") end; loadpath = "./viz2D_out/"; anim = Animation(loadpath,String[])
  #  println("Animation directory: $(anim.dir)")

  time,time_kyrs          = 0.0, 0.0;
  err = 100;
  it = 0;
  #nt = 100
  for it=1:nt

      # Perform an advection step for temperature
      AdvectTracers!(Tracers, Grid, (Vx,Vz), dt, AdvectionMethod);

      time        =   time + dt;                                            # Keep track of evolved time

      if mod(it,1000)==0  # print progress
          println(" Timestep $it = $((time)) ")

          Tr_coord = Tracers.coord; Tr_coord = hcat(Tr_coord...)';       # extract array with coordinates of tracers
          p1 = scatter(Tr_coord[:,1], Tr_coord[:,2], zcolor = Tracers.T, m = (:inferno , 0.8, Plots.stroke(0.01, :black)), markersize=1.0, xlims=(0,1), ylims=(0,1))
          plot(p1); frame(anim)
      end

  end

  x_km, z_km  =   x, z;

  # compute analytical solution
  if      Method=="ConstantZ"
    # Analytical location of tracers
    Tr_coord_anal = Tracers0.coord; Tr_coord_anal = hcat(Tr_coord_anal...)';       # extract array with coordinates of tracers
    Tr_coord_anal[:,2] .= Tr_coord_anal[:,2] .- 0.5;

  elseif      Method=="Rotation"
      # Should arrive @ same point
      Tr_coord_anal = Tracers0.coord; Tr_coord_anal = hcat(Tr_coord_anal...)';       # extract array with coordinates of tracers

  end
  Tr_coord = Tracers.coord; Tr_coord = hcat(Tr_coord...)';

  fname       = "TracersAdvection_2D_$(Method)_$(AdvectionMethod)";

  Terror      =  Tr_coord - Tr_coord_anal;   # error

  if CreatePlots
    # create plot
    p1 = plot(Tr_coord_anal[:,1], Tr_coord_anal[:,2], seriestype = :scatter, markersize=5.0, markershape=:circle, markerstrokecolor=:red, markercolor=:white, linewidth=0, label="analytics")
        plot!(Tr_coord[:,1],      Tr_coord[:,2],      seriestype = :scatter, markersize=2.0, markershape=:circle, markerstrokecolor=:black, linewidth=0, label="numerics")
    plot(p1);
    png(fname)
  end

  error = norm(Terror[:],2)/length(Terror[:]);
  return error;        # return error
end

function test_AdvectTracers3D(Method="ConstantZ",  AdvectionMethod="RK2")
  # Test tracer advection in 3D

  # Model parameters
  W,L,H                   =   1.,  1., 1.;                      # Width, Length, Height
  Nx,Ny, Nz               =   65, 65, 65;                     # resolution
  dx,dy,dz                =   W/(Nx-1), L/(Ny-1),H/(Nz-1);           # grid size [m]

  # Set up model geometry
  x,y,z                   =   0:dx:((Nx-1)*dx), 0:dy:((Ny-1)*dy), 0:dz:((Nz-1)*dz);
  coords                  =   collect(Iterators.product(x,y,z))                               # generate coordinates from 1D coordinate vectors
  X,Y,Z                   =   (x->x[1]).(coords), (x->x[2]).(coords), (x->x[3]).(coords);     # transfer coords to 2D arrays
  Grid, Spacing           =   (x,y,z), (dx,dy,dz);

  # generate tracer coordinates
  Nxt, Nyt, Nzt           =   50, 50, 50          # num of tracers in x,z
  Wt,Lt,Ht                =   0.1, 0.1, 0.1;
  dxt,dyt,dzt             =   Wt/(Nxt-1), Lt/(Nyt-1), Ht/(Nzt-1);
  xt,yt,zt                =   -Wt/2:dxt:Wt/2, -Lt/2:dyt:Lt/2, -Ht/2:dz:Ht/2;
  coordst                 =   collect(Iterators.product(xt,yt,zt))             # generate coordinates from 1D coordinate vectors
  Xt,Yt,Zt                =   (x->x[1]).(coordst), (x->x[2]).(coordst),(x->x[3]).(coordst);     # transfer coords to 2D arrays

  Xc                      =   0.5;
  Yc                      =   0.5;
  Zc                      =   0.75;
  if      Method=="ConstantZ"
    # Simple advection in z-direction
    Vx        =   X.*0;
    Vy        =   Y.*0;
    Vz        =   Z.*0 .- 1;
    TotalTime =   0.5;                         # thermal cooling age

  elseif  Method=="Rotation_alongY"
    # Rigid body rotation:
    Vx        =    (Z .- 0.5);
    Vy        =   X.*0;
    Vz        =   -(X .- 0.5);
    TotalTime =   2*pi;

  elseif  Method=="Rotation_alongX"
    # Rigid body rotation:
    Vy        =    (Z .- 0.5);
    Vx        =   Y.*0;
    Vz        =   -(Y .- 0.5);
    TotalTime =   2*pi;

  elseif  Method=="Rotation_alongZ"
    # Rigid body rotation:
    Vy        =   -(X .- 0.5);
    Vx        =    (Y .- 0.5);
    Vz        =   Z.*0.0;
    TotalTime =   2*pi;
    Xc                      =   0.5;
    Yc                      =   0.75;
    Zc                      =   0.5;
  end

  maxVel      =   max( maximum(abs.(Vx)), maximum(abs.(Vy)), maximum(abs.(Vz)) );
  dt          =   min(dx,dy,dz)/2.0/maxVel;   # stable timestep (to ensure we move no more than 0.5*(dx|dz) per timesstep )
  numTime     =   ceil(TotalTime/dt);
  nt          =   Int(numTime);
  dt          =   TotalTime/nt;

  # Create tracer structure
  Xt .=  Xt .+ Xc; Yt .=  Yt .+ Yc; Zt .=  Zt .+ Zc;   # shift coordinates of tracers

  new_tracer  =   Tracer(num=1, coord=[Xt[1]; Yt[1]; Zt[1]], T=0.);            # Create new tracer
  Tracers     =   StructArray([new_tracer]);                # Create tracer array
  for i=firstindex(Xt)+1:lastindex(Xt)
    new_tracer  =   Tracer(num=i, coord=[Xt[i]; Yt[i]; Zt[i]], T=0.);          # Create new tracer
    push!(Tracers, new_tracer);                             # Add new point to existing array
  end
  Tracers0 = copy(Tracers);         # create a copy of the original tracers

  time    = 0.0;
  err     = 100;
  it = 0;
  #nt = 100
  for it=1:nt

      # Perform an advection step for temperature
     AdvectTracers!(Tracers, Grid, (Vx,Vy,Vz), dt, AdvectionMethod);

      time        =   time + dt;                                            # Keep track of evolved time

      if mod(it,1000)==0  # print progress
          println(" Timestep $it = $((time)) ")

          Tr_coord = Tracers.coord; Tr_coord = hcat(Tr_coord...)';       # extract array with coordinates of tracers
          p1 = scatter(Tr_coord[:,1], Tr_coord[:,3], zcolor = Tracers.T, m = (:inferno , 0.8, Plots.stroke(0.01, :black)), markersize=1.0, xlims=(0,1), ylims=(0,1))
          plot(p1);
      end

  end

  x_km, z_km  =   x, z;

  # compute analytical solution
  if      Method=="ConstantZ"
    # Analytical location of tracers
    Tr_coord_anal = Tracers0.coord; Tr_coord_anal = hcat(Tr_coord_anal...)';       # extract array with coordinates of tracers
    Tr_coord_anal[:,3] .= Tr_coord_anal[:,3] .- 0.5;

  else

      # Should arrive @ same point
      Tr_coord_anal = Tracers0.coord; Tr_coord_anal = hcat(Tr_coord_anal...)';       # extract array with coordinates of tracers

  end
  Tr_coord = Tracers.coord; Tr_coord = hcat(Tr_coord...)';

  fname       = "TracersAdvection_3D_$(Method)_$(AdvectionMethod)";

  Terror      =  Tr_coord - Tr_coord_anal;   # error

  if CreatePlots
    # create plot
    p1 = plot(Tr_coord_anal[:,1], Tr_coord_anal[:,3], seriestype = :scatter, markersize=5.0, markershape=:circle, markerstrokecolor=:red, markercolor=:white, linewidth=0, label="analytics")
        plot!(Tr_coord[:,1],      Tr_coord[:,3],      seriestype = :scatter, markersize=2.0, markershape=:circle, markerstrokecolor=:black, linewidth=0, label="numerics")
    plot(p1);
    png(fname)
  end

  error = norm(Terror[:],2)/length(Terror[:]);
  return error;        # return error
end

# ===================================================================================================
if 1==1

@testset "Interpolation" begin
  @test test_Interpolation("2D", "Linear")      ≈  2.9880072526933544e-5    atol=1e-8;
  @test test_Interpolation("2D", "Quadratic")   ≈  9.72101621785853e-7      atol=1e-8;
  @test test_Interpolation("2D", "Cubic")       ≈  5.752004968381702e-7     atol=1e-8;
  @test test_Interpolation("3D", "Linear")      ≈  5.758940704356386e-8     atol=1e-8;
  @test test_Interpolation("3D", "Cubic")       ≈  1.0924731223751025e-10   atol=1e-8;
end;

@testset "2D semi-lagrangian advection" begin
  @test test_SemiLagrangian2D("ConstantZ","Cubic", "RK4"    )   ≈ 7.257090264260951e-6    atol=1e-8;
  @test test_SemiLagrangian2D("ConstantZ","Linear","Euler"  )   ≈ 0.00026385729240145655  atol=1e-6;

  @test test_SemiLagrangian2D("Rotation", "Linear", "Euler" )   ≈ 0.0009061021221256843   atol=1e-8;
  @test test_SemiLagrangian2D("Rotation", "Linear", "RK2"   )   ≈ 0.000892833826724465    atol=1e-8;
  @test test_SemiLagrangian2D("Rotation", "Linear", "RK4"   )   ≈ 0.0008928336158448796   atol=1e-8;

  @test test_SemiLagrangian2D("Rotation", "Quadratic", "RK2")   ≈ 2.2477968913260075e-5   atol=1e-8;

  @test test_SemiLagrangian2D("Rotation", "Cubic", "Euler"  )   ≈ 0.00013400741577631413  atol=1e-8;
  @test test_SemiLagrangian2D("Rotation", "Cubic", "RK2"    )   ≈ 1.0762102257417536e-6   atol=1e-8;
  @test test_SemiLagrangian2D("Rotation", "Cubic", "RK4"    )   ≈ 1.069882664228221e-6    atol=1e-8;
end;

@testset "3D semi-lagrangian advection" begin
  @test test_SemiLagrangian3D("ConstantZ","Linear","RK2"   )          ≈ 6.822115778040405e-5    atol=1e-8;
  @test test_SemiLagrangian3D("ConstantZ","Cubic", "RK2"   )          ≈ 1.1038374380141209e-6   atol=1e-8;

  @test test_SemiLagrangian3D("Rotation_alongY","Linear", "Euler")    ≈ 0.00017064553949646487  atol=1e-8;
  @test test_SemiLagrangian3D("Rotation_alongY","Cubic", "RK2"   )    ≈ 1.899516467222456e-6    atol=1e-8;
  @test test_SemiLagrangian3D("Rotation_alongY","Cubic", "RK2"   )    ≈ 1.899516467222456e-6    atol=1e-8;

  @test test_SemiLagrangian3D("Rotation_alongX","Cubic", "RK2"   )    ≈ 1.8995164672224572e-6    atol=1e-8;
  @test test_SemiLagrangian3D("Rotation_alongY","Cubic", "RK2"   )    ≈ 1.899516467222456e-6    atol=1e-8;
  @test test_SemiLagrangian3D("Rotation_alongZ","Cubic", "RK2"   )    ≈ 1.8995164672224786e-6   atol=1e-8;
end;

@testset "2D tracer advection" begin
  @test test_AdvectTracers2D("Rotation","Euler")        ≈ 0.00012187143210808163  atol=1e-8;
  @test test_AdvectTracers2D("Rotation","RK2")          ≈ 3.132104548068979e-7    atol=1e-8;
  @test test_AdvectTracers2D("Rotation","RK4")          ≈ 1.5661110135721713e-7   atol=1e-8;

  @test test_AdvectTracers2D("ConstantZ","Euler")       ≈ 0.0   atol=1e-8;
end;

@testset "3D tracer advection" begin
  @test test_AdvectTracers3D("ConstantZ","Euler")       ≈ 0.0   atol=1e-8;

  @test test_AdvectTracers3D("Rotation_alongX","RK2")   ≈ 1.60725931984687e-7    atol=1e-8;
  @test test_AdvectTracers3D("Rotation_alongY","RK2")   ≈ 1.607259319846869e-7   atol=1e-8;
  @test test_AdvectTracers3D("Rotation_alongZ","RK2")   ≈ 1.6256291449684137e-7  atol=1e-8;

  @test test_AdvectTracers3D("Rotation_alongY","Euler") ≈ 3.169439899086836e-5   atol=1e-8;
  @test test_AdvectTracers3D("Rotation_alongY","RK4")   ≈ 8.037493317907375e-8   atol=1e-8;
end;

end

#@time test_SemiLagrangian3D("ConstantZ","Cubic", "RK4"    )
#test_SemiLagrangian3D("ConstantZ","Linear", "Euler"    )
#test_SemiLagrangian2D("ConstantZ","Linear", "Euler"    )
#test_AdvectTracers2D("Rotation","Euler")
#test_AdvectTracers2D("ConstantZ","Euler")

#test_AdvectTracers3D("Rotation_alongX","RK2")

#test_SemiLagrangian2D("Rotation", "Linear", "RK2"   )