text:
using Flux, PerceptualColourMaps, ProgressMeter, Distributions, ArgParse

# change to CairoMakie if without gpu
# tried an if-else but that borked my Julia for some reason...
# (ERROR: LoadError: InitError: could not load library libharfbuzz-gobject.so")
using GLMakie

s = ArgParseSettings()

@add_arg_table! s begin
    "--resolution-upscaling"
    help = "Modifier to 1920x1080:\n -2 = 480x270,\n -1 = 960x540,\n 1 = 1920x1080 (DEFAULT),\n 2 = 3840x2160"
    arg_type = Int
    default = 1
end

Arguments = parse_args(s)

@info Arguments

scale = Arguments["resolution-upscaling"]
res = (1920,1080)
if res == -2
    res = (480,270)
elseif res == -1
    res = (960,540)
elseif res == 2
    res = (3840,2160)
end

include("visualisation.jl")

# comment out these two lines to watch live, it's much faster without though
set_window_config!(framerate = Inf, vsync = false)
inline!(true)

fps = 5
# 1D

# no noise 1d linear

model = Chain(Dense(1,3),Dense(3,3),Dense(3,1))

for l in model.layers
    l.W .= .1
    l.b .= 0.0
end

f(x) = 5.0*x

x = collect(0.0:0.01:1.0)
y = f.(x);
pred = [model([i]) for i in x]

params_ = Flux.params(model);

Data = [([x[i]],[y[i]]) for i in 1:length(x)];

loss(model,x,y) = sum((model(x).-y).^2.0)

WatchNetwork(model,Data,loss,y,fps=fps,Epochs=1000,
res=(1920*2,1080*2),
filename="linear.mp4",opt=Descent(0.0001))

# Cauchy + normal noise in 1D

model = Chain(Dense(1,3),Dense(3,3),Dense(3,1))

for l in model.layers
    l.W .= .1
    l.b .= 0.
end
η=0.25
f(x) = 5.0*x+(rand(Cauchy(0,η))%5)+randn()*η

x = collect(0.0:0.01:1.0)
y = f.(x);
pred = [model([i]) for i in x]

loss(x,y) = sum((model(x).-y).^2.0)
params_ = Flux.params(model);

Data = [([x[i]],[y[i]]) for i in 1:length(x)];
cb = () -> @show sum([loss(x,y) for (x,y) in Data])/length(Data)

loss(model,x,y) = sum((model(x).-y).^2.0)

WatchNetwork(model,Data,loss,y,fps=fps,Epochs=600,res=(1920*2,1080*2), opt=Descent(0.001))

# 2D

# linear problem

model = Chain(Dense(2,2),Dense(2,2))
f(x,y)= 0.4*x < y

N = 32
x = range(0,stop=1.,length=N)
y = copy(x)
t = [f(i,j)==1 ? [0,1] : [1,0] for i in x for j in y];
pred = [model([i,j]) for i in x for j in y]

loss(x,y) = Flux.Losses.logitbinarycrossentropy(softmax(model(x)),y)
params_ = Flux.params(model);

Data = [([i,j],f(i,j)==1 ? [0,1] : [1,0]) for i in x for j in y];

n = [d==[0,1] for d in map(x->x[2],Data)]
m = [d==[1,0] for d in map(x->x[2],Data)]

Data_old = copy(Data)

if sum(n) < sum(m)
    for i in 1:abs(sum(n) - sum(m))
        push!(Data,rand(Data_old[n]))
    end
else
    for i in 1:abs(sum(n) - sum(m))
        push!(Data,rand(Data_old[m]))
    end
end
n = [d==[0,1] for d in map(x->x[2],Data)]
m = [d==[1,0] for d in map(x->x[2],Data)]
sum(n),sum(m)

loss(model,x,y) = Flux.Losses.logitbinarycrossentropy(softmax(model(x)),y)

WatchNetwork(model,Data,loss,t,fps=5,Epochs=2000,res=(1920*2,1080*2),
filename="2d-linear.mp4",opt=Descent(0.005))

# nonlinear with linear network

model = Chain(Dense(2,2),Dense(2,2),Dense(2,2))

f(x,y)= x<y || x^2+y^2 < 0.5
N = 32
x = range(0,stop=1.,length=N)
y = copy(x)
t = [f(i,j)==1 ? [0,1] : [1,0] for i in x for j in y];
pred = [model([i,j]) for i in x for j in y]

loss(x,y) = Flux.Losses.logitbinarycrossentropy(softmax(model(x)),y)
params_ = Flux.params(model);

Data = [([i,j],f(i,j)==1 ? [0,1] : [1,0]) for i in x for j in y];

n = [d==[0,1] for d in map(x->x[2],Data)]
m = [d==[1,0] for d in map(x->x[2],Data)]

Data_old = copy(Data)

if sum(n) < sum(m)
    for i in 1:abs(sum(n) - sum(m))
        push!(Data,rand(Data_old[n]))
    end
else
    for i in 1:abs(sum(n) - sum(m))
        push!(Data,rand(Data_old[m]))
    end
end
n = [d==[0,1] for d in map(x->x[2],Data)]
m = [d==[1,0] for d in map(x->x[2],Data)]
sum(n),sum(m)

loss(model,x,y) = Flux.Losses.logitbinarycrossentropy(softmax(model(x)),y)

WatchNetwork(model,Data,loss,t,fps=5,Epochs=1000,res=(1920*2,1080*2),
filename="2d-nonlinear-bad.mp4",opt=Descent(0.001))

# non-linear problem with activations

model = Chain(Dense(2,2,tanh),Dense(2,2,tanh),Dense(2,2,tanh))

f(x,y)= x<y || x^2+y^2 < 0.5
N = 32
x = range(0,stop=1.,length=N)
y = copy(x)
t = [f(i,j)==1 ? [0,1] : [1,0] for i in x for j in y];
pred = [model([i,j]) for i in x for j in y]

loss(x,y) = Flux.Losses.logitbinarycrossentropy(softmax(model(x)),y)
params_ = Flux.params(model);

Data = [([i,j],f(i,j)==1 ? [0,1] : [1,0]) for i in x for j in y];

n = [d==[0,1] for d in map(x->x[2],Data)]
m = [d==[1,0] for d in map(x->x[2],Data)]

Data_old = copy(Data)

if sum(n) < sum(m)
    for i in 1:abs(sum(n) - sum(m))
        push!(Data,rand(Data_old[n]))
    end
else
    for i in 1:abs(sum(n) - sum(m))
        push!(Data,rand(Data_old[m]))
    end
end
n = [d==[0,1] for d in map(x->x[2],Data)]
m = [d==[1,0] for d in map(x->x[2],Data)]
sum(n),sum(m)

loss(model,x,y) = Flux.Losses.logitbinarycrossentropy(softmax(model(x)),y)

WatchNetwork(model,Data,loss,t,fps=5,Epochs=3000,res=(1920*2,1080*2),
filename="2d-nonlinear.mp4",opt=Descent(0.001),neuron_colour=RGBAf0(0.6,0.6,0,1))

# Polkadots!

model = Chain(Dense(2,4,tanh),Dense(4,8,tanh),Dense(8,8,tanh),Dense(8,4,tanh),Dense(4,2,tanh))

f(x,y)= x^2+y^2 < 0.1 || (x-0.6)^2+(y-0.5)^2 < 0.1
N = 32
x = range(0,stop=1.,length=N)
y = copy(x)
t = [f(i,j)==1 ? [0,1] : [1,0] for i in x for j in y];
pred = [model([i,j]) for i in x for j in y]

loss(x,y) = Flux.Losses.logitbinarycrossentropy(softmax(model(x)),y)
params_ = Flux.params(model);

Data = [([i,j],f(i,j)==1 ? [0,1] : [1,0]) for i in x for j in y];

n = [d==[0,1] for d in map(x->x[2],Data)]
m = [d==[1,0] for d in map(x->x[2],Data)]

Data_old = copy(Data)

if sum(n) < sum(m)
    for i in 1:abs(sum(n) - sum(m))
        push!(Data,rand(Data_old[n]))
    end
else
    for i in 1:abs(sum(n) - sum(m))
        push!(Data,rand(Data_old[m]))
    end
end
n = [d==[0,1] for d in map(x->x[2],Data)]
m = [d==[1,0] for d in map(x->x[2],Data)]
sum(n),sum(m)

loss(model,x,y) = Flux.Losses.logitbinarycrossentropy(softmax(model(x)),y)

WatchNetwork(model,Data,loss,t,fps=5,Epochs=6000,res=(1920*2,1080*2),
filename="2d-polkadot.mp4",opt=ADAGrad(0.01),neuron_colour=RGBAf0(0.6,0.6,0,1))

# Some sort of triangle thing

model = Chain(Dense(2,4,tanh),Dense(4,8,tanh),Dense(8,8,tanh),Dense(8,4,tanh),Dense(4,2,tanh))

f(x,y)= 0.4*x < y && 0.4*y < x
N = 32
x = range(0,stop=1.,length=N)
y = copy(x)
t = [f(i,j)==1 ? [0,1] : [1,0] for i in x for j in y];
pred = [model([i,j]) for i in x for j in y]

loss(x,y) = Flux.Losses.logitbinarycrossentropy(softmax(model(x)),y)
params_ = Flux.params(model);

Data = [([i,j],f(i,j)==1 ? [0,1] : [1,0]) for i in x for j in y];

n = [d==[0,1] for d in map(x->x[2],Data)]
m = [d==[1,0] for d in map(x->x[2],Data)]

Data_old = copy(Data)

if sum(n) < sum(m)
    for i in 1:abs(sum(n) - sum(m))
        push!(Data,rand(Data_old[n]))
    end
else
    for i in 1:abs(sum(n) - sum(m))
        push!(Data,rand(Data_old[m]))
    end
end
n = [d==[0,1] for d in map(x->x[2],Data)]
m = [d==[1,0] for d in map(x->x[2],Data)]
sum(n),sum(m)

loss(model,x,y) = Flux.Losses.logitbinarycrossentropy(softmax(model(x)),y)

WatchNetwork(model,Data,loss,t,fps=5,Epochs=6000,res=(1920*2,1080*2),
filename="2d-triangle.mp4",opt=ADAGrad(0.001),neuron_colour=RGBAf0(0.6,0.6,0,1))

max_stars_repo_name: Jerboa-app/MLPVis

max_issues_repo_name: Jerboa-app/MLPVis

max_forks_repo_name: Jerboa-app/MLPVis


text:
using FFTW
using PyPlot

function dsdemo6()
    println("Half-band filter design")
    println("ALPHA VERSION")

    if false
        # Use one of the canned examples
        f1, f2 = exampleHBF(2)
    else
        fp = 0.9 * 0.25
        delta = undbv(-100)
        fig = figure(1)
        fig.clf()
        fig.canvas.set_window_title("designHBF Iterations")
        f1, f2, = designHBF(fp, delta, true)
    end

    # interleave the even and odd decimated impulse responses
    Nimp = 2^11
    imp = simulateHBF([1; zeros(Nimp-1)], f1[1], f2[1])

    mag = abs.(fft(imp))
    mag = mag[1:end÷2+1]

    fig = figure(2)
    fig.clf()
    fig.canvas.set_window_title("designHBF Result")

    plot(range(0, 0.5, length=length(mag)), dbv(mag))
    figureMagic([0, 0.5], 0.05, 2, [-150, 3], 10, 5)

    return nothing
end

ext: jl

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
"""
$(TYPEDEF)

Defines a jump process with a rate (i.e. hazard, intensity, or propensity) that
does not *explicitly* depend on time. More precisely, one where the rate
function is constant *between* the occurrence of jumps. For detailed examples
and usage information see the
- [Tutorial](https://diffeq.sciml.ai/stable/tutorials/discrete_stochastic_example/)

## Fields

$(FIELDS)

## Examples
Suppose `u[1]` gives the amount of particles and `p[1]` the probability per time
each particle can decay away. A corresponding `ConstantRateJump` for this jump
process is
```julia
rate(u,p,t) = p[1]*u[1]
affect!(integrator) = integrator.u[1] -= 1
crj = ConstantRateJump(rate, affect!)
```
Notice, here that `rate` changes in time, but is constant between the occurrence
of jumps (when `u[1]` will decrease).
"""
struct ConstantRateJump{F1,F2} <: AbstractJump
  """Function `rate(u,p,t)` that returns the jump's current rate."""
  rate::F1
  """Function `affect(integrator)` that updates the state for one occurrence of the jump."""
  affect!::F2
end

"""
$(TYPEDEF)

Defines a jump process with a rate (i.e. hazard, intensity, or propensity) that
may explicitly depend on time. More precisely, one where the rate function is
allowed to change *between* the occurrence of jumps. For detailed examples and
usage information see the
- [Tutorial](https://diffeq.sciml.ai/stable/tutorials/discrete_stochastic_example/)

## Fields

$(FIELDS)

## Examples
Suppose `u[1]` gives the amount of particles and `t*p[1]` the probability per
time each particle can decay away. A corresponding `VariableRateJump` for this
jump process is
```julia
rate(u,p,t) = t*p[1]*u[1]
affect!(integrator) = integrator.u[1] -= 1
crj = VariableRateJump(rate, affect!)
```

## Notes
- **`VariableRateJump`s result in `integrator`s storing an effective state type
  that wraps the main state vector.** See [`ExtendedJumpArray`](@ref) for
  details on using this object. Note that the presence of *any*
  `VariableRateJump`s will result in all `ConstantRateJump`, `VariableRateJump`
  and callback `affect!` functions receiving an integrator with `integrator.u`
  an [`ExtendedJumpArray`](@ref).
- Must be used with `ODEProblem`s or `SDEProblem`s to be correctly simulated
  (i.e. can not currently be used with `DiscreteProblem`s).
- Salis H., Kaznessis Y.,  Accurate hybrid stochastic simulation of a system of
  coupled chemical or biochemical reactions, Journal of Chemical Physics, 122
  (5), DOI:10.1063/1.1835951 is used for calculating jump times with
  `VariableRateJump`s within ODE/SDE integrators.
"""
struct VariableRateJump{R,F,I,T,T2} <: AbstractJump
  """Function `rate(u,p,t)` that returns the jump's current rate."""
  rate::R
  """Function `affect(integrator)` that updates the state for one occurrence of the jump."""
  affect!::F
  idxs::I
  rootfind::Bool
  interp_points::Int
  save_positions::Tuple{Bool,Bool}
  abstol::T
  reltol::T2
end

VariableRateJump(rate,affect!;
                   idxs = nothing,
                   rootfind=true,
                   save_positions=(true,true),
                   interp_points=10,
                   abstol=1e-12,reltol=0) = VariableRateJump(
                              rate,affect!,idxs,
                              rootfind,interp_points,
                              save_positions,abstol,reltol)

struct RegularJump{iip,R,C,MD}
    rate::R
    c::C
    numjumps::Int
    mark_dist::MD
    function RegularJump{iip}(rate,c,numjumps::Int; mark_dist = nothing) where iip
      new{iip,typeof(rate),typeof(c),typeof(mark_dist)}(rate,c,numjumps,mark_dist)
    end
end

DiffEqBase.isinplace(::RegularJump{iip,R,C,MD}) where {iip,R,C,MD} = iip

RegularJump(rate,c,numjumps::Int; kwargs...) = RegularJump{DiffEqBase.isinplace(rate,4)}(rate,c,numjumps;kwargs...)

# deprecate old call
function RegularJump(rate,c,dc::AbstractMatrix; constant_c=false, mark_dist = nothing)
  @warn("The RegularJump interface has changed to be matrix-free. See the documentation for more details.")
  function _c(du,u,p,t,counts,mark)
    c(dc,u,p,t,mark)
    mul!(du,dc,counts)
  end
  RegularJump{true}(rate,_c,size(dc,2);mark_dist=mark_dist)
end

"""
$(TYPEDEF)

Optimized representation for `ConstantRateJump`s that can be represented in mass
action form, offering improved performance within jump algorithms compared to
`ConstantRateJump`. For detailed examples and usage information see the
- [Main
  Docs](https://diffeq.sciml.ai/stable/types/jump_types/#Defining-a-Mass-Action-Jump)
- [Tutorial](https://diffeq.sciml.ai/stable/tutorials/discrete_stochastic_example/)

### Constructors
- `MassActionJump(reactant_stoich, net_stoich; scale_rates=true,
  param_idxs=nothing)`

Here `reactant_stoich` denotes the reactant stoichiometry for each reaction and
`net_stoich` the net stoichiometry for each reaction.

## Fields

$(FIELDS)

## Keyword Arguments
- `scale_rates=true`, whether to rescale the reaction rate constants according
  to the stoichiometry.
- `nocopy=false`, whether the `MassActionJump` can alias the `scaled_rates` and
  `reactant_stoch` from the input. Note, if `scale_rates=true` this will
  potentially modify both of these.
- `param_idxs=nothing`, indexes in the parameter vector, `JumpProblem.prob.p`,
  that correspond to each reaction's rate.

See the tutorial and main docs for details.

## Examples
An SIR model with `S + I --> 2I` at rate β as the first reaction and `I --> R`
at rate ν as the second reaction can be encoded by
```julia
p        = (β=1e-4, ν=.01)
u0       = [999, 1, 0]       # (S,I,R)
tspan    = (0.0, 250.0)
rateidxs = [1, 2]           # i.e. [β,ν]
reactant_stoich =
[
  [1 => 1, 2 => 1],         # 1*S and 1*I
  [2 => 1]                  # 1*I
]
net_stoich =
[
  [1 => -1, 2 => 1],        # -1*S and 1*I
  [2 => -1, 3 => 1]         # -1*I and 1*R
]
maj = MassActionJump(reactant_stoich, net_stoich; param_idxs=rateidxs)
prob = DiscreteProblem(u0, tspan, p)
jprob = JumpProblem(prob, Direct(), maj)
```

## Notes
- By default reaction rates are rescaled when constructing the `MassActionJump`
  as explained in the [main
  docs](https://diffeq.sciml.ai/stable/types/jump_types/#Defining-a-Mass-Action-Jump).
  Disable this with the kwarg `scale_rates=false`.
- Also see the [main
  docs](https://diffeq.sciml.ai/stable/types/jump_types/#Defining-a-Mass-Action-Jump)
  for how to specify reactions with no products or no reactants.

"""
struct MassActionJump{T,S,U,V} <: AbstractMassActionJump
  """The (scaled) reaction rate constants."""
  scaled_rates::T
  """The reactant stoichiometry vectors."""
  reactant_stoch::S
  """The net stoichiometry vectors."""
  net_stoch::U
  """Parameter mapping functor to identify reaction rate constants with parameters in `p` vectors."""
  param_mapper::V

  function MassActionJump{T,S,U,V}(rates::T, rs_in::S, ns::U, pmapper::V, scale_rates::Bool, useiszero::Bool, nocopy::Bool) where {T <: AbstractVector, S, U, V}
    sr  = nocopy ? rates : copy(rates)
    rs = nocopy ? rs_in : copy(rs_in)
    for i in eachindex(rs)
      if useiszero && (length(rs[i]) == 1) && iszero(rs[i][1][1])
        rs[i] = typeof(rs[i])()
      end
    end

    if scale_rates && !isempty(sr)
      scalerates!(sr, rs)
    end
    new(sr, rs, ns, pmapper)
  end
  function MassActionJump{Nothing,Vector{S},Vector{U},V}(::Nothing, rs_in::Vector{S}, ns::Vector{U}, pmapper::V, scale_rates::Bool, useiszero::Bool, nocopy::Bool) where {S<:AbstractVector, U<:AbstractVector, V}
    rs = nocopy ? rs_in : copy(rs_in)
    for i in eachindex(rs)
      if useiszero && (length(rs[i]) == 1) && iszero(rs[i][1][1])
        rs[i] = typeof(rs[i])()
      end
    end
    new(nothing, rs, ns, pmapper)
  end
  function MassActionJump{T,S,U,V}(rate::T, rs_in::S, ns::U, pmapper::V, scale_rates::Bool, useiszero::Bool, nocopy::Bool) where {T <: Number, S, U, V}
    rs = rs_in
    if useiszero && (length(rs) == 1) && iszero(rs[1][1])
      rs = typeof(rs)()
    end
    sr = scale_rates ? scalerate(rate, rs) : rate
    new(sr, rs, ns, pmapper)
  end
  function MassActionJump{Nothing,S,U,V}(::Nothing, rs_in::S, ns::U, pmapper::V, scale_rates::Bool, useiszero::Bool, nocopy::Bool) where {S, U, V}
    rs = rs_in
    if useiszero && (length(rs) == 1) && iszero(rs[1][1])
      rs = typeof(rs)()
    end
    new(nothing, rs, ns, pmapper)
  end

end
MassActionJump(usr::T, rs::S, ns::U, pmapper::V; scale_rates = true, useiszero = true, nocopy=false) where {T,S,U,V} = MassActionJump{T,S,U,V}(usr, rs, ns, pmapper, scale_rates, useiszero, nocopy)
MassActionJump(usr::T, rs, ns; scale_rates = true, useiszero = true, nocopy=false) where {T <: AbstractVector,S,U} = MassActionJump(usr, rs, ns, nothing; scale_rates=scale_rates, useiszero=useiszero, nocopy=nocopy)
MassActionJump(usr::T, rs, ns; scale_rates = true, useiszero = true, nocopy=false) where {T <: Number,S,U} = MassActionJump(usr, rs, ns, nothing; scale_rates=scale_rates, useiszero=useiszero, nocopy=nocopy)

# with parameter indices or mapping, multiple jump case
function MassActionJump(rs, ns; param_idxs=nothing, param_mapper=nothing, scale_rates = true, useiszero = true, nocopy=false)
  if param_mapper === nothing
    (param_idxs === nothing) && error("If no parameter indices are given via param_idxs, an explicit parameter mapping must be passed in via param_mapper.")
    pmapper = MassActionJumpParamMapper(param_idxs)
  else
    (param_idxs !== nothing) && error("Only one of param_idxs and param_mapper should be passed.")
    pmapper = param_mapper
  end

  MassActionJump(nothing, nocopy ? rs : copy(rs), ns, pmapper; scale_rates=scale_rates,
                 useiszero=useiszero, nocopy=true)
end

using_params(maj::MassActionJump{T,S,U,Nothing}) where {T,S,U} = false
using_params(maj::MassActionJump) = true
using_params(maj::Nothing) = false
@inline get_num_majumps(maj::MassActionJump) = length(maj.scaled_rates)
@inline get_num_majumps(maj::Nothing) = 0

struct MassActionJumpParamMapper{U}
  param_idxs::U
end

# create the initial parameter vector for use in a MassActionJump
# Note these are unscaled
function (ratemap::MassActionJumpParamMapper{U})(params) where {U <: AbstractArray}
  [params[pidx] for pidx in ratemap.param_idxs]
end

# Note this is unscaled
function (ratemap::MassActionJumpParamMapper{U})(params) where {U <: Int}
  params[ratemap.param_idxs]
end

# update a maj with parameter vectors
function (ratemap::MassActionJumpParamMapper{U})(maj::MassActionJump, newparams; scale_rates, kwargs...) where {U <: AbstractArray}
  for i in 1:get_num_majumps(maj)
    maj.scaled_rates[i] = newparams[ratemap.param_idxs[i]]
  end
  scale_rates && scalerates!(maj.scaled_rates, maj.reactant_stoch)
  nothing
end

function to_collection(ratemap::MassActionJumpParamMapper{Int})
  MassActionJumpParamMapper([ratemap.param_idxs])
end

function Base.merge!(pmap1::MassActionJumpParamMapper{U}, pmap2::MassActionJumpParamMapper{U}) where {U <: AbstractVector}
  append!(pmap1.param_idxs, pmap2.param_idxs)
end

function Base.merge!(pmap1::MassActionJumpParamMapper{U}, pmap2::MassActionJumpParamMapper{V}) where {U <: AbstractVector, V <: Int}
  push!(pmap1.param_idxs, pmap2.param_idxs)
end

function Base.merge(pmap1::MassActionJumpParamMapper{Int}, pmap2::MassActionJumpParamMapper{Int})
  MassActionJumpParamMapper([pmap1.param_idxs, pmap2.param_idxs])
end

"""
  update_parameters!(maj::MassActionJump, newparams; scale_rates=true)

Updates the passed in MassActionJump with the parameter values in `newparams`.

Notes:
  - Requires the jump to have been constructed with a user-passed `param_idxs` or `param_mapper`.
  - `scale_rates=true` will scale the parameter representing the jump rate by an
    appropriate combinatoric factor. i.e for 3A --> B at rate k it will scale
    k --> k/3!.
"""
function update_parameters!(maj::MassActionJump, newparams; scale_rates=true, kwargs...)
  (maj.param_mapper === nothing) && error("MassActionJumps must be constructed with param_idxs or a param_mapper to be updateable.")
  maj.param_mapper(maj, newparams; scale_rates, kwargs)
end

struct JumpSet{T1,T2,T3,T4} <: AbstractJump
  variable_jumps::T1
  constant_jumps::T2
  regular_jump::T3
  massaction_jump::T4
end
JumpSet(vj, cj, rj, maj::MassActionJump{S,T,U,V}) where {S <: Number,T,U,V} = JumpSet(vj, cj, rj, check_majump_type(maj))

JumpSet(jump::ConstantRateJump) = JumpSet((),(jump,),nothing,nothing)
JumpSet(jump::VariableRateJump) = JumpSet((jump,),(),nothing,nothing)
JumpSet(jump::RegularJump)      = JumpSet((),(),jump,nothing)
JumpSet(jump::AbstractMassActionJump)   = JumpSet((),(),nothing,jump)
function JumpSet(; variable_jumps=(), constant_jumps=(),
                   regular_jumps=nothing, massaction_jumps=nothing)
  JumpSet(variable_jumps, constant_jumps, regular_jumps, massaction_jumps)
end
JumpSet(jb::Nothing) = JumpSet()

# For Varargs, use recursion to make it type-stable
JumpSet(jumps::AbstractJump...) = JumpSet(split_jumps((), (), nothing, nothing, jumps...)...)

# handle vector of mass action jumps
function JumpSet(vjs, cjs, rj, majv::Vector{T}) where {T <: MassActionJump}
  if isempty(majv)
    error("JumpSets do not accept empty mass action jump collections; use \"nothing\" instead.")
  end

  maj = setup_majump_to_merge(majv[1].scaled_rates, majv[1].reactant_stoch, majv[1].net_stoch, majv[1].param_mapper)
  for i = 2:length(majv)
    massaction_jump_combine(maj, majv[i])
  end

  JumpSet(vjs, cjs, rj, maj)
end

@inline get_num_majumps(jset::JumpSet) = get_num_majumps(jset.massaction_jump)

@inline split_jumps(vj, cj, rj, maj) = vj, cj, rj, maj
@inline split_jumps(vj, cj, rj, maj, v::VariableRateJump, args...) = split_jumps((vj..., v), cj, rj, maj, args...)
@inline split_jumps(vj, cj, rj, maj, c::ConstantRateJump, args...) = split_jumps(vj, (cj..., c), rj, maj, args...)
@inline split_jumps(vj, cj, rj, maj, c::RegularJump, args...) = split_jumps(vj, cj, regular_jump_combine(rj,c), maj, args...)
@inline split_jumps(vj, cj, rj, maj, c::MassActionJump, args...) = split_jumps(vj, cj, rj, massaction_jump_combine(maj,c), args...)
@inline split_jumps(vj, cj, rj, maj, j::JumpSet, args...) = split_jumps((vj...,j.variable_jumps...),
                                                                        (cj..., j.constant_jumps...),
                                                                        regular_jump_combine(rj,j.regular_jump),
                                                                        massaction_jump_combine(maj,j.massaction_jump), args...)

regular_jump_combine(rj1::RegularJump,rj2::Nothing) = rj1
regular_jump_combine(rj1::Nothing,rj2::RegularJump) = rj2
regular_jump_combine(rj1::Nothing,rj2::Nothing) = rj1
regular_jump_combine(rj1::RegularJump,rj2::RegularJump) = error("Only one regular jump is allowed in a JumpSet")

# functionality to merge two mass action jumps together
check_majump_type(maj::MassActionJump{S,T,U,V}) where {S <: Number,T,U,V} = setup_majump_to_merge(maj.scaled_rates, maj.reactant_stoch, maj.net_stoch, maj.param_mapper)
check_majump_type(maj::MassActionJump{Nothing,T,U,V}) where {T,U,V} = setup_majump_to_merge(maj.scaled_rates, maj.reactant_stoch, maj.net_stoch, maj.param_mapper)

# if given containers of rates and stoichiometry directly create a jump
function setup_majump_to_merge(sr::T, rs::AbstractVector{S}, ns::AbstractVector{U}, pmapper) where {T <: AbstractVector, S <: AbstractArray, U <: AbstractArray}
   MassActionJump(sr, rs, ns, pmapper; scale_rates=false)
end

# if just given the data for one jump (and not in a container) wrap in a vector
function setup_majump_to_merge(sr::S, rs::T, ns::U, pmapper) where {S <: Number, T <: AbstractArray, U <: AbstractArray}
  MassActionJump([sr], [rs], [ns], (pmapper === nothing) ? pmapper : to_collection(pmapper); scale_rates=false)
end

# if no rate field setup yet
function setup_majump_to_merge(::Nothing, rs::T, ns::U, pmapper) where {T <: AbstractArray, U <: AbstractArray}
  MassActionJump(nothing, [rs], [ns], (pmapper === nothing) ? pmapper : to_collection(pmapper); scale_rates=false)
end

# when given a collection of reactions to add to maj
function majump_merge!(maj::MassActionJump{U,<:AbstractVector{V},<:AbstractVector{W},X}, sr::U, rs::AbstractVector{V}, ns::AbstractVector{W}, param_mapper) where {U <: Union{AbstractVector,Nothing}, V <: AbstractVector, W <: AbstractVector, X}
  (U <: AbstractVector) && append!(maj.scaled_rates, sr)
  append!(maj.reactant_stoch, rs)
  append!(maj.net_stoch, ns)
  if maj.param_mapper === nothing
    (param_mapper === nothing) || error("Error, trying to merge a MassActionJump with a parameter mapping to one without a parameter mapping.")
  else
    merge!(maj.param_mapper, param_mapper)
  end
  maj
end

# when given a single jump's worth of data to add to maj
function majump_merge!(maj::MassActionJump{U,V,W,X}, sr::T, rs::S1, ns::S2, param_mapper) where {T <: Union{Number,Nothing}, S1 <: AbstractArray, S2 <: AbstractArray, U <: Union{AbstractVector{T},Nothing}, V <: AbstractVector{S1}, W <: AbstractVector{S2}, X}
  (T <: Number) && push!(maj.scaled_rates, sr)
  push!(maj.reactant_stoch, rs)
  push!(maj.net_stoch, ns)
  if maj.param_mapper === nothing
    (param_mapper === nothing) || error("Error, trying to merge a MassActionJump with a parameter mapping to one without a parameter mapping.")
  else
    merge!(maj.param_mapper, param_mapper)
  end

  maj
end

# when maj only stores a single jump's worth of data (and not in a collection)
# create a new jump with the merged data stored in vectors
function majump_merge!(maj::MassActionJump{T,S,U,V}, sr::T, rs::S, ns::U, param_mapper::V) where {T <: Union{Number,Nothing}, S <: AbstractArray{<:Pair}, U <: AbstractArray{<:Pair}, V}
  rates = (T <: Nothing) ? nothing : [maj.scaled_rates, sr]
  if maj.param_mapper === nothing
    (param_mapper === nothing) || error("Error, trying to merge a MassActionJump with a parameter mapping to one without a parameter mapping.")
    return MassActionJump(rates, [maj.reactant_stoch, rs], [maj.net_stoch, ns], param_mapper; scale_rates=false)
  else
    return MassActionJump(rates, [maj.reactant_stoch, rs], [maj.net_stoch, ns], merge(maj.param_mapper,param_mapper); scale_rates=false)
  end
end

massaction_jump_combine(maj1::MassActionJump, maj2::Nothing) = maj1
massaction_jump_combine(maj1::Nothing, maj2::MassActionJump) = maj2
massaction_jump_combine(maj1::Nothing, maj2::Nothing) = maj1
massaction_jump_combine(maj1::MassActionJump, maj2::MassActionJump) = majump_merge!(maj1, maj2.scaled_rates, maj2.reactant_stoch, maj2.net_stoch, maj2.param_mapper)

##### helper methods for unpacking rates and affects! from constant jumps #####
function get_jump_info_tuples(constant_jumps)
  if (constant_jumps !== nothing) && !isempty(constant_jumps)
    rates    = ((c.rate for c in constant_jumps)...,)
    affects! = ((c.affect! for c in constant_jumps)...,)
  else
    rates    = ()
    affects! = ()
  end

  rates, affects!
end

function get_jump_info_fwrappers(u, p, t, constant_jumps)
  RateWrapper   = FunctionWrappers.FunctionWrapper{typeof(t),Tuple{typeof(u), typeof(p), typeof(t)}}
  AffectWrapper = FunctionWrappers.FunctionWrapper{Nothing,Tuple{Any}}

  if (constant_jumps !== nothing) && !isempty(constant_jumps)
    rates    = [RateWrapper(c.rate) for c in constant_jumps]
    affects! = [AffectWrapper(x->(c.affect!(x);nothing)) for c in constant_jumps]
  else
    rates    = Vector{RateWrapper}()
    affects! = Vector{AffectWrapper}()
  end

  rates, affects!
end


text:
#!/usr/bin/env julia

using SEM
using LinearAlgebra,Plots, UnPack
#----------------------------------#
kx=2.0
ky=2.0
kt=2.0
ν =1.0

function utrue(x,y,t)
    ut = @. sin(kx*pi*x)*sin.(ky*pi*y)*cos(kt*pi*t)
    return ut
end

function setIC!(u,x,y,t)
    u .= utrue(x,y,t)
    return
end

function setBC!(ub,x,y,t)
    ub .= @. 0+0*x
    return
end

function setForcing!(f,x,y,t)
    ut = utrue(x,y,t)
    f  .= @. ut*((kx^2+ky^2)*pi^2*ν)
    f .+= @. - sin(kx*pi*x)*sin(ky*pi*y)*sin(kt*pi*t)*(kt*pi)
    return
end

function setVisc!(ν,x,y,t)
    ν .= @. 1+0*x
    return
end

function callback!(dfn::Diffusion)
    @unpack fld,msh = dfn
    @unpack time, istep = dfn.tstep

    ut = utrue(msh.x,msh.y,time[1])
    u  = fld.u
    er = norm(ut-u,Inf)
    println("Step $(istep[1]), Time=$(time[1]), er=$er")
    plt = meshplt(u,msh)
    plt = plot!(zlims=(-1,1))
    display(plt)
    #frame(anim)
    return
end

#----------------------------------#
Ex = 5; nr1 = 8;
Ey = 5; ns1 = 8;

m1 = Mesh(nr1,ns1,Ex,Ey)
bc = ['D','D','D','D']
T = Diffusion(bc,m1,Tf=1.0,dt=0.01) # temperature
#----------------------------------#
function dudt!(dudt,u,msh::Mesh)

    @unpack x,y = msh
    # can play directly with full solution (uh + ub)
    # since we're no longer solving a BVP (explicit time method)

    dudt = lapl(u,msh) + f(msh)
    return dudt
end

prob = ODEProblem(dudt!,)
sol  = solve(prob,SSPRK33(),dt=0.01,saveat=0.25)
#----------------------------------#


text:
@testset "sumlessthan.jl" begin
    @testset "SumLessThan()" begin
        trailer = SeaPearl.Trailer()

        x = SeaPearl.IntVar(2, 6, "x", trailer)
        y = SeaPearl.IntVar(2, 3, "y", trailer)
        ax = SeaPearl.IntVarViewMul(x, 3, "3x")
        minusY = SeaPearl.IntVarViewOpposite(y, "-y")
        vars = [x, y, ax, minusY]
        constraint = SeaPearl.SumLessThan(vars, 10, trailer)

        @test constraint in x.onDomainChange
        @test constraint in y.onDomainChange
        @test constraint.active.value
        @test constraint.numberOfFreeVars.value == 4
        @test constraint.sumOfFixedVars.value == 0
        @test constraint.freeIds == [1, 2, 3, 4]
    end
    @testset "propagate!(::SumLessThan)" begin
        trailer = SeaPearl.Trailer()
        vars = SeaPearl.AbstractIntVar[]

        x = SeaPearl.IntVar(2, 3, "x", trailer)
        push!(vars, x)
        y = SeaPearl.IntVar(5, 15, "y", trailer)
        ax = SeaPearl.IntVarViewMul(x, 3, "3x")
        push!(vars, ax)
        minusY = SeaPearl.IntVarViewOpposite(y, "-y")
        push!(vars, minusY)
        constraint = SeaPearl.SumLessThan(vars, -5, trailer)

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test length(y.domain) == 3
        @test 13 in y.domain
        @test !(12 in y.domain)
        @test !(3 in ax.domain)

    end
end


text:
"""
Direct with rejection sampling
"""

mutable struct RDirectJumpAggregation{T,S,F1,F2,RNG,DEPGR} <: AbstractSSAJumpAggregator
  next_jump::Int
  prev_jump::Int
  next_jump_time::T
  end_time::T
  cur_rates::Vector{T}
  sum_rate::T
  ma_jumps::S
  rates::F1
  affects!::F2
  save_positions::Tuple{Bool,Bool}
  rng::RNG
  dep_gr::DEPGR
  max_rate::T
  counter::Int
  counter_threshold
end

function RDirectJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, maj::S,
                                rs::F1, affs!::F2, sps::Tuple{Bool,Bool}, rng::RNG;
                                num_specs, counter_threshold=length(crs), dep_graph=nothing, kwargs...) where {T,S,F1,F2,RNG,DEPGR}
    # a dependency graph is needed and must be provided if there are constant rate jumps
    if dep_graph === nothing
        if (get_num_majumps(maj) == 0) || !isempty(rs)
            error("To use ConstantRateJumps with the Rejection Direct (RDirect) algorithm a dependency graph must be supplied.")
        else
            dg = make_dependency_graph(num_specs, maj)
        end
    else
        dg = dep_graph

        # make sure each jump depends on itself
        add_self_dependencies!(dg)
    end

    max_rate = maximum(crs)
    return RDirectJumpAggregation{T,S,F1,F2,RNG,typeof(dg)}(nj, nj, njt, et, crs, sr, maj, rs, affs!, sps, rng,
        dg, max_rate, 0, counter_threshold)
end

############################# Required Functions #############################

# creating the JumpAggregation structure (tuple-based constant jumps)
function aggregate(aggregator::RDirect, u, p, t, end_time, constant_jumps,
                    ma_jumps, save_positions, rng; kwargs...)

    # handle constant jumps using function wrappers
    rates, affects! = get_jump_info_fwrappers(u, p, t, constant_jumps)

    build_jump_aggregation(RDirectJumpAggregation, u, p, t, end_time, ma_jumps,
                          rates, affects!, save_positions, rng; num_specs = length(u), kwargs...)
end

# set up a new simulation and calculate the first jump / jump time
function initialize!(p::RDirectJumpAggregation, integrator, u, params, t)
    p.end_time = integrator.sol.prob.tspan[2]
    fill_rates_and_sum!(p, u, params, t)
    p.max_rate = maximum(p.cur_rates)
    generate_jumps!(p, integrator, u, params, t)
    nothing
end

"execute one jump, changing the system state and updating rates"
function execute_jumps!(p::RDirectJumpAggregation, integrator, u, params, t)
    # execute jump
    u = update_state!(p, integrator, u)

    # update rates
    update_dependent_rates!(p, u, params, t)
    nothing
end

"calculate the next jump / jump time"
function generate_jumps!(p::RDirectJumpAggregation, integrator, u, params, t)
    # if no more events possible there is nothing to do
    sum_rate = p.sum_rate
    if nomorejumps!(p, sum_rate)
        return nothing
    end
    @unpack rng, cur_rates, max_rate = p

    num_rxs = length(cur_rates)
    counter = 0
    rx = trunc(Int, rand(rng) * num_rxs)+1
    @inbounds while cur_rates[rx] < rand(rng) * max_rate
        rx = trunc(Int, rand(rng) * num_rxs)+1
        counter += 1
    end

    p.counter = counter
    p.next_jump = rx
    p.next_jump_time = t + randexp(p.rng) / sum_rate
    nothing
end

######################## SSA specific helper routines #########################
function update_dependent_rates!(p::RDirectJumpAggregation, u, params, t)
    @inbounds dep_rxs = p.dep_gr[p.next_jump]
    @unpack ma_jumps, rates, cur_rates, sum_rate = p
    num_majumps = get_num_majumps(ma_jumps)
    max_rate_increased = false
    @inbounds for rx in dep_rxs
        @inbounds new_rate = calculate_jump_rate(ma_jumps, num_majumps, rates, u,params,t,rx)
        sum_rate += new_rate - cur_rates[rx]
        if new_rate > p.max_rate
            p.max_rate = new_rate
            max_rate_increased = true
        end
        cur_rates[rx] = new_rate
    end
    if !max_rate_increased && p.counter > p.counter_threshold
        p.max_rate = maximum(p.cur_rates)
    end

    p.sum_rate = sum_rate
    nothing
end


text:
module CompositeStructs

using Core: apply_type
using Base: datatype_fieldtypes, unwrap_unionall
using MacroTools

export @composite

# given an expression like :(Complex{T}) and a module in which the
# relevant symbols are defined, reconstruct a type object without
# calling @eval
reconstruct_type(__module__, s::Symbol) = getfield(__module__, s)
function reconstruct_type(__module__, ex)
    isexpr(ex, :curly) || error("Invalid @composite syntax.")
    foldl((t,x) -> apply_type(t, TypeVar(x)), ex.args[2:end], init=reconstruct_type(__module__,ex.args[1]))
end

# convert a type or UnionAll back to an expression
to_expr(t::DataType) = isempty(t.parameters) ? t.name.name : :($(t.name.name){$(map(to_expr,t.parameters)...)})
to_expr(t::TypeVar) = t.name
to_expr(t::Symbol) = t
to_expr(t::UnionAll) = :($(to_expr(t.body)) where {$(t.var.lb) <: $(t.var.name) <: $(t.var.ub)})

# given an expression like :(Complex{T}) and a module in which the
# relevant symbols are defined, return an array with the struct's
# fields and their type signatures, in this case: [:(re::T), (im::T)]
function reconstruct_fields(__module__, ex)
    t = reconstruct_type(__module__, ex)
    (t isa UnionAll) && error("Spliced type $ex must not have any free type parameters.")
    map(zip(fieldnames(t),datatype_fieldtypes(t))) do (x,T)
        :($x::$(to_expr(T)))
    end
end

@doc join(readlines(joinpath(@__DIR__, "../README.md"))[6:end], "\n")
macro composite(ex)

    if !(
        ((iskwdef = @capture(ex, @kwdef(structdef_) | Base.@kwdef(structdef_))) || @capture(ex, structdef_)) &&
        @capture(structdef, struct ParentTypeDecl_ parent_body__ end | mutable struct ParentTypeDecl_ parent_body__ end) &&
        @capture(ParentTypeDecl, (ParentType_ <: _) | ParentType_) &&
        @capture(ParentType, ParentName_{ParentTypeArgs__} | ParentName_)
    )
        error("Invalid @composite syntax.")
    end

    generic_child_constructors = []
    concrete_child_constructors = []
    parent_body′ = []
    explicit_parent_fields = []
    constructor_args = []

    # seem to be a bug if you put these in a single @capture
    _field_name(ex) = @capture(ex,name_Symbol::T_=val_) || @capture(ex,name_Symbol::T_) || @capture(ex,name_Symbol=val_) || @capture(ex,name_Symbol) ? name : nothing
    _field_kw(ex)   = @capture(ex,((name_::T_=val_) | (name_=val_))) ? Expr(:kw, name, val) : _field_name(ex)
    _field_decl(ex) = iskwdef && @capture(ex,(decl_ = val_)) ? decl : ex
    _strip_type_bound(ex) = @capture(ex, T_ <: _) ? T : ex

    for x in parent_body
        if @capture(x, ChildType_...) && @capture(ChildType, ChildName_{__} | ChildName_)
            child_fields = (reconstruct_fields(__module__, ChildType)...,)
            child_field_names = _field_name.(child_fields)
            append!(parent_body′, child_fields)
            child_instance = gensym()
            push!(generic_child_constructors,  :($child_instance = $ChildName(; filter(((k,_),)->(k in $child_field_names), kw)...)))
            push!(concrete_child_constructors, :($child_instance = $ChildType(; filter(((k,_),)->(k in $child_field_names), kw)...)))
            for f in child_field_names
                push!(constructor_args, :($child_instance.$f))
            end
        elseif _field_name(x) != nothing
            push!(explicit_parent_fields, x)
            push!(constructor_args, _field_name(x))
            push!(parent_body′, _field_decl(x))
        else
            push!(parent_body′, x)
        end
    end

    structdef.args[3] = :(begin $(parent_body′...) end)

    if !iskwdef
        esc(structdef)
    else
        ret = quote Core.@__doc__ $structdef end
        push!(ret.args, quote
            function $ParentName(;$(_field_kw.(explicit_parent_fields)...), kw...)
                $(generic_child_constructors...)
                $ParentName($(constructor_args...))
            end
        end)
        if ParentTypeArgs!=nothing
            ParentTypeArgsStripped = map(_strip_type_bound, ParentTypeArgs)
            push!(ret.args, quote
                function $ParentName{$(ParentTypeArgsStripped...)}(;$(_field_kw.(explicit_parent_fields)...), kw...) where {$(ParentTypeArgs...)}
                    $(concrete_child_constructors...)
                    $ParentName{$(ParentTypeArgsStripped...)}($(constructor_args...))
                end
            end)
        end
        esc(ret)
    end

end

end


text:
__precompile__()

module MathExpressions

import Lazy: @>
using DataStructures

"""

"""
abstract type Expression <: Number end

"""
Atoms are the not quite indivisible building blocks of mathematical expressions. The expression
    D(x, t) ≅ m * x(t)^2 + C
has atoms `D`, `x`, `t`, `≅`, `m`, `*` , `^`, `2`, `+` and `C`.

All `T` where `{T <: Atom}` are assumed by default to have the fields:

*) `name::Symbol`

*) `value` which may contain a value associated with that atom or `nothing` otherwise
"""
abstract type Atom <: Expression end

abstract type AbstractOperation <: Expression end
abstract type AbstractDomain end
abstract type AbstractOperation <: Expression end

include("Domains.jl")
include("Atoms.jl")
include("Operations.jl")

export

end


text:
struct Model
    x0::Vector{Int}
    F::Function
    H::Function
    ν::Matrix{Int}
    bounds::Vector{Vector{Float64}}
    mapping::Vector{Vector{Int}}
    obs::Vector{Int}
    obsname::Vector{Symbol}
    hor::Vector{Int}
    ns::Int
    nr::Int
    ps::Vector{Int}
    Model(x0::Vector{Int}, F::Function, H::Function, ν::Matrix{Int}, bounds::Vector{Vector{Float64}},mapping::Vector{Vector{Int}}, obs::Vector{Int}, obsname::Vector{Symbol}, hor::Vector{Int}, ns::Int, nr::Int, ps::Vector{Int}) = new(copy(x0), F, H, copy(ν), bounds, copy(mapping), copy(obs), obsname, copy(hor), ns, nr, ps)
end

Model(x0::Vector{Int}, F::Function, H::Function, ν::Matrix{Int}, bounds::Vector{Vector{Float64}}; obs::Vector{Int} = collect(1:length(x0)), obsname::Vector{Symbol} = Symbol[], hor::Vector{Int} = fill(2, length(x0)), ns::Int = length(x0), nr::Int = size(ν)[1], mapping::Vector{Vector{Int}} = [Int[i] for i in 1:length(bounds[1])], ps::Vector{Int} = collect(1:size(ν)[1])) = Model(convert(Vector{Int}, x0), F, H, convert(Matrix{Int}, ν), bounds, mapping, obs, obsname, hor, ns, nr, ps)

function Base.show(io::IO, m::Model)
    print("Model:\n\t\tx0: ", m.x0, ",\n\t\tObservables: ", m.obs, ",\n\t\tFunction: ", m.F, ",\n\t\tState Change Matrix: ", m.ν, ",\n\t\tHOR: ", m.hor, ",\n\t\t", m.ns, " species,\n\t\t", m.nr, " reactions.")
end

max_stars_repo_name: pzuliani/SPICE

max_issues_repo_name: pzuliani/SPICE

max_forks_repo_name: pzuliani/SPICE


text:
@testset "SingleField" begin
    @test_throws ErrorException("Unknown SingleField aggregation") SingleField("weird_aggregation", "field", "name")
end

@testset "StringAgg" begin
    @test_throws ErrorException("Unknown StringAgg aggregation") StringAgg("second", "field", "name")
    @test_throws ErrorException("maxStringBytes must be a positive integer") StringAgg("first", "field", "name", maxStringBytes="a")
    @test_throws ErrorException("maxStringBytes must be a positive integer") StringAgg("first", "field", "name", maxStringBytes=-1)
end


text:
using OrdinaryDiffEq, Flux, DiffEqFlux, Test

@testset "NeuralODE struct" begin
    x = Float32[2.; 0.]
    tspan = (0.0f0,25.0f0)
    dudt = Chain(Dense(2,50,tanh),Dense(50,2))
    slver = Tsit5()
    slver_kwargs = Dict(:save_everystep=>false, :save_start=>false)

    funcsolv = neural_ode(dudt,x,tspan,slver;slver_kwargs...)

    node = NeuralODE(dudt,tspan,slver,slver_kwargs)

    @test node(x) ≈ funcsolv
end


text:
## Load the packages.
using Tullio
using Statistics
using BenchmarkTools

## Create the coarsening function.
function coarsen_field!(lo, hi)
    filter = 1//8 * ones(2, 2, 2)
    @tullio lo[i, j, k] = hi[2(i-1)+l, 2(j-1)+m, 2(k-1)+n] * filter[l, m, n]
end

## Set up the grids.
n_hi = 128
n_lo = 64

a_hi = rand(n_hi, n_hi, n_hi)
a_lo = zeros(n_lo, n_lo, n_lo)

## Compute.
@btime coarsen_field!(a_lo, a_hi)

## Check the output.
println(mean(a_lo) ≈ mean(a_hi))

## Create the refine function.
function refine_field!(hi, lo)
    @tullio hi[2(i-1)+l, 2(j-1)+m, 2(k-1)+n] = lo[i, j, k] (l in 1:2, m in 1:2, n in 1:2)
    return
end

## Set up the grids.
n_hi = 128
n_lo = 64

a_lo = rand(n_lo, n_lo, n_lo)
a_hi = zeros(n_hi, n_hi, n_hi)

## Compute.
@btime refine_field!(a_hi, a_lo)

## Check the output.
println(mean(a_lo) ≈ mean(a_hi))


text:
import Base: LineEdit, REPL

function idl_repl()

    # Setup idl prompt
    prompt = LineEdit.Prompt("IDL> ";
                             prompt_prefix=Base.text_colors[:blue],
                             prompt_suffix=Base.text_colors[:white])

    repl = Base.active_repl

    prompt.on_done = REPL.respond(repl,prompt) do line
        ok2, line, msg = convert_continuations(line)
        if !ok2
            println(msg)
            println()
            return
        end
        ok2, line, msg = replace_interpolated_vars(line)
        if !ok2
            println(msg)
            println()
            return
        end
        ec = execute_converted(line)
        nothing
    end

    main_mode = repl.interface.modes[1]

    # replace existing IDL REPL if present
    i_mode = find_prompt_in_modes(repl.interface.modes, "IDL> ")
    if i_mode < 1
        push!(repl.interface.modes,prompt)
    else
        repl.interface.modes[i_mode] = prompt
    end

    hp = main_mode.hist
    hp.mode_mapping[:idl] = prompt
    prompt.hist = hp

    const idl_keymap = Dict{Any,Any}(
            '>' => function (s,args...)
                if isempty(s)
                    if !haskey(s.mode_state,prompt)
                        s.mode_state[prompt] = LineEdit.init_state(repl.t,prompt)
                    end
                    LineEdit.transition(s,prompt)
                else
                    LineEdit.edit_insert(s,'>')
                end
            end)

    search_prompt, skeymap = LineEdit.setup_search_keymap(hp)
    mk = REPL.mode_keymap(main_mode)

    b = Dict{Any,Any}[skeymap, mk, LineEdit.history_keymap, LineEdit.default_keymap,
                      LineEdit.escape_defaults]
    prompt.keymap_dict = LineEdit.keymap(b)

    main_mode.keymap_dict = LineEdit.keymap_merge(main_mode.keymap_dict, idl_keymap);
    nothing
end

function find_prompt_in_modes(modes, name)
    j = -1
    for (i,mode) in enumerate(modes)
        if :prompt in fieldnames(mode) && mode.prompt == name
            j = i
            break
        end
    end
    return j
end


text:
struct Polynomial{R}
    coeffs::Vector{R}
end

function (p::Polynomial)(x)
    v = p.coeffs[end]
    for i = (length(p.coeffs)-1):-1:1
        v = v*x + p.coeffs[i]
    end
    return v
end

(p::Polynomial)() = p(9)

p = Polynomial([1,9,250])

println(p(3))

println(p())

max_stars_repo_name: JulialangOrgCN/howtojulia

max_issues_repo_name: JulialangOrgCN/howtojulia

max_forks_repo_name: JulialangOrgCN/howtojulia


text:
using Test
using MPI
using Random

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    import CUDA
    ArrayType = CUDA.CuArray
else
    ArrayType = Array
end

MPI.Init()

comm = MPI.COMM_WORLD

root = 0
matsize = (17,17)

for T in Base.uniontypes(MPI.MPIDatatype)
    # This test depends on the stability of the rng and we have observed with
    # CUDA.jl that it is not gurantueed that the same number of rand calls will
    # occur on each rank. (This is a hypothesis). To be sure we shall seed the rng
    # just before we call rand.
    Random.seed!(17)
    A = ArrayType(rand(T, matsize))
    B = MPI.Comm_rank(comm) == root ? A : similar(A)
    MPI.Bcast!(B, comm; root=root)
    @test B == A
end

# Char
A = ['s', 't', 'a', 'r', ' ', 'w', 'a', 'r', 's']
B = MPI.Comm_rank(comm) == root ? A : similar(A)
MPI.Bcast!(B, comm; root=root)
@test B == A

g = x -> x^2 + 2x - 1
if MPI.Comm_rank(comm) == root
    f = g
else
    f = nothing
end
f = MPI.bcast(f, root, comm)
@test f(3) == g(3)
@test f(5) == g(5)
@test f(7) == g(7)

A = Dict("foo" => "bar")
if MPI.Comm_rank(comm) == root
    B = A
else
    B = nothing
end
B = MPI.bcast(B, root, comm)
@test B["foo"] == "bar"

MPI.Finalize()
@test MPI.Finalized()


text:
using BioMedQuery.Entrez

"""
    export_citation(entrez_email, pmid::Int64, citation_type, output_file,verbose)

Export, to an output file, the citation for PubMed article identified by the given pmid

### Arguments
* `citation_type::ASCIIString`: At the moment supported types include: "endnote"
"""
function export_citation(entrez_email, pmid::Int64, citation_type, output_file,
    overwrite = true, verbose=false)
    fetch_dic = Dict("db"=>"pubmed","tool" =>"BioJulia", "email" => entrez_email,
    "retmode" => "xml", "rettype"=>"null")
    efetch_response = efetch(fetch_dic, [pmid])
    if verbose
        xmlASCII2file(efetch_response, "./efetch.xml")
    end
    efetch_dict = eparse(efetch_response)
    config = Dict(:type => citation_type, :output_file => output_file, :overwrite=>overwrite)
    save_article_citations(efetch_dict, config, verbose)
end

"""
    export_citation(entrez_email, pmids::Vector{Int64}, citation_type, output_file,verbose)

Export, to an output file, the citation for collection of PubMed articles identified by the given pmids

### Arguments
* `citation_type::ASCIIString`: At the moment supported types include: "endnote"
"""
function export_citation(entrez_email, pmids::Vector{Int64}, citation_type, output_file,
    overwrite=true, verbose=false)
    fetch_dic = Dict("db"=>"pubmed","tool" =>"BioJulia", "email" => entrez_email,
    "retmode" => "xml", "rettype"=>"null")
    efetch_response = efetch(fetch_dic, pmids)
    if verbose
        xmlASCII2file(efetch_response, "./efetch.xml")
    end
    efetch_dict = eparse(efetch_response)
    config = Dict(:type => citation_type, :output_file => output_file, :overwrite=>overwrite)
    save_article_citations(efetch_dict, config, verbose)
end


text:
# for each parameter sim_select

# generate

println(now())
using Debug
# Initialize simulation functions
fpath="../src/new_SimulationFunctions.jl" # Loads up DataFile , SolveBalances , Balances, Kinetics, Control
include(fpath)

# Establish simulation parameters
TSTART = 0.0;
Ts = 1.0; # 0.1
TSTOP = 150.0; # 250

# Initialize intial data dictionary
data_dictionary = DataFile(TSTART,TSTOP,Ts)
len_cpa = length(data_dictionary["CONTROL_PARAMETER_ARRAY"])
len_rca = length(data_dictionary["RATE_CONSTANT_ARRAY"])
len_sca = length(data_dictionary["SATURATION_CONSTANT_ARRAY"])
sca_indices,sca_values = [],[]
for i in 1:len_sca
  val = data_dictionary["SATURATION_CONSTANT_ARRAY"][i]
if val !=0
  push!(sca_indices,i)
  push!(sca_values,val)
end
end

sca_pack = zip(sca_indices,sca_values)

 # CPA, let's give access to gain AND power
 cpa_gain = data_dictionary["CONTROL_PARAMETER_ARRAY"][:,1] # Gain
 cpa_order= data_dictionary["CONTROL_PARAMETER_ARRAY"][:,2] # Power

 # get lengths of new arrays
 len_cpa_gain = length(cpa_gain)
 len_cpa_order = length(cpa_order)
 len_sca_values = length(sca_values)

 len_sca_prime = length(sca_values)
   A = len_rca + len_cpa_gain
   B = A + len_cpa_order
   C = B + len_sca_prime

# read parameter arrays resulting from POETs

PC = readdlm("PC.dat")

# PC is an array of parameter arrays
PC_results =[]
fraction = 1
popsize = Int(length(PC[1,:])*fraction)
# popsize = 2 # to debug
println("Model Population Size ", popsize)
for i in 1:popsize
  println("Regen ",i," of ",popsize)
  # Refactor data dictionary for model execution
  # need to specify implicit strucutre of parameter_array
  P = PC[:,i]
  data_dictionary["RATE_CONSTANT_ARRAY"] = P[1:len_rca]
  data_dictionary["CONTROL_PARAMETER_ARRAY"][:,1] = P[len_rca+1:A]
  data_dictionary["CONTROL_PARAMETER_ARRAY"][:,2] = P[A+1:B]
# tune VEGF/NFATc1 parameters
  control_parameter_array=data_dictionary["CONTROL_PARAMETER_ARRAY"]
  control_parameter_array[6,1]=0.001 # GAIN NFATc1 -> Ecadherin_Bcatenin
  control_parameter_array[6,2]=1.5 # ORDER
  data_dictionary["CONTROL_PARAMETER_ARRAY"]=control_parameter_array
  rate_constant_array=data_dictionary["RATE_CONSTANT_ARRAY"]
# degradation rates
  rate_constant_array[34]=25.0
  rate_constant_array[36]=100.0
  rate_constant_array[37]=100.0
# PKC->RAF activation rate
  rate_constant_array[11]=0.0
  data_dictionary["RATE_CONSTANT_ARRAY"]=rate_constant_array

  # Rebuild SCA
  sca_pack = zip(sca_indices,P[B+1:C]) # new values
  for (i,val) in sca_pack
  data_dictionary["SATURATION_CONSTANT_ARRAY"][i]=val
  end

  # Run to Steady State
 #println("running to SS")
  XSS = EstimateSteadyState(data_dictionary)
  initial_condition_array = vec(XSS);
  data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
  xss_data_dict = deepcopy(data_dictionary)

  simulations = ["edit_file_none.jl","edit_file_tgfb.jl","edit_file_mapkin.jl","edit_file_ss.jl","edit_file_ssecad.jl","edit_file_ss_tgfb3ab.jl","edit_file_ss_dnlef1.jl","edit_file_tgfb3.jl","edit_file_dnsmad.jl","edit_file_tgfb3_dnlef1.jl","edit_file_tgfb_tgfb3ab.jl"]
# Need new simulations for TGFB, VEGF, TGFB+VEGF, VEGF+VIVIT, TGFB+VIVIT,
  simulations = ["edit_file_none.jl","edit_file_tgfb.jl","edit_file_vegf.jl","edit_file_tgfb_vegf.jl",
  "edit_file_tgfb_vegf_vi.jl","edit_file_tgfb_vi.jl","edit_file_vegf_vi.jl"]
  sim_data = [] # Keep data in RAM
  for sim in simulations
    ##### What data_dictionary should be loaded here??? clearly the active one
    active_data_dictionary = deepcopy(xss_data_dict)
    (TSIM,X) =  SolveModelGenericPOETS(TSTART,TSTOP,Ts,active_data_dictionary,string("../src/",sim))
    push!(sim_data,X)
    # now reset the data dictionary
  end
  push!(PC_results,sim_data)

end

function trapz{Tx<:Number, Ty<:Number}(x::Vector{Tx}, y::Vector{Ty})
    # Trapezoidal integration rule
    local n = length(x)
    if (length(y) != n)
        error("Vectors 'x', 'y' must be of same length")
    end
    r = zero(zero(Tx) + zero(Ty))
    if n == 1; return r; end
    for i in 2:n
        r += (x[i] - x[i-1]) * (y[i] + y[i-1])
    end
    #= correction -h^2/12 * (f'(b) - f'(a))
    ha = x[2] - x[1]
    he = x[end] - x[end-1]
    ra = (y[2] - y[1]) / ha
    re = (y[end] - y[end-1]) / he
    r/2 - ha*he/12 * (re - ra)
    =#
    return r/2
end
#@show "BYPASS ROBUSTNESS CALC"
function compute_robustness(BASE,PERT,TIME_B,TIME_P,K=0)

#base_auc = trapz(TIME_B,BASE)
#pert_auc = trapz(TIME_P,PERT)
#ROBUSTNESS = (pert_auc)/(K+base_auc)

return PERT[t48]
#return ROBUSTNESS
end

# # ok we will overwrite these indices to plot nfat and sp1

# ecadherin protein (unbound)
index_active_nfatc = 34
#index_active_nfatc = 16 # MRNA!!
# ecad_bcat index
#index_active_nfatc_bcat = 53
# vimentin protein
index_psp1 = 50 #AP1_SP1_P
#index_psp1 = 85 # MRNA
t0  = 1
t24 = 74
t48 = 149
t72 = 224
t96 = 298
ti = t0
tf = t48
TIME = collect(ti:tf)

# sim indices
# 1 edit_file_none.jl
# 2 edit_file_tgfb.jl
# 3 edit_file_vegf.jl
# 4 edit_file_tgfb_vegf.jl
# 5 edit_file_tgfb_vegf_vi.jl
# 6 edit_file_tgfb_vi.jl
# 7 edit_file_vegf_vi.jl

robustness_array = []

const K = 0 # offset for robustness calc
for state in PC_results # iterate over parameter sets

state_robustness_values = []

# calculate tgfb robustness
# sample from 0 to 48 hrs and integrate
# get time series for simulation
base_ = state[1]
pert_ = state[2]

# select species and time range from simulation
base_ecad = base_[ti:tf,index_active_nfatc] #+ base_[ti:tf,index_active_nfatc_bcat]
pert_ecad = pert_[ti:tf,index_active_nfatc] #+ pert_[ti:tf,index_active_nfatc_bcat]
base_vim  = base_[ti:tf,index_psp1]
pert_vim  = pert_[ti:tf,index_psp1]

# compute robustness
robustness_ecad = compute_robustness(base_ecad,pert_ecad,TIME,TIME,K)
robustness_vim  = compute_robustness(base_vim,pert_vim,TIME,TIME,K)
push!(state_robustness_values, robustness_ecad)
push!(state_robustness_values, robustness_vim)

# calculate vegf robustness
# sample from 0 to 48 hrs and integrate
base_ = state[1]
pert_ = state[3]

base_ecad = base_[ti:tf,index_active_nfatc] #+ base_[ti:tf,index_active_nfatc_bcat]
pert_ecad = pert_[ti:tf,index_active_nfatc] #+ pert_[ti:tf,index_active_nfatc_bcat]
base_vim  = base_[ti:tf,index_psp1]
pert_vim  = pert_[ti:tf,index_psp1]

robustness_ecad = compute_robustness(base_ecad,pert_ecad,TIME,TIME,K)
robustness_vim  = compute_robustness(base_vim,pert_vim,TIME,TIME,K)
push!(state_robustness_values, robustness_ecad)
push!(state_robustness_values, robustness_vim)

# calculate tgfb vegf robustness
# sample from 0 to 48 hrs and integrate
base_ = state[1]
pert_ = state[4]

base_ecad = base_[ti:tf,index_active_nfatc] #+ base_[ti:tf,index_active_nfatc_bcat]
pert_ecad = pert_[ti:tf,index_active_nfatc] #+ pert_[ti:tf,index_active_nfatc_bcat]
base_vim  = base_[ti:tf,index_psp1]
pert_vim  = pert_[ti:tf,index_psp1]

robustness_ecad = compute_robustness(base_ecad,pert_ecad,TIME,TIME,K)
robustness_vim  = compute_robustness(base_vim,pert_vim,TIME,TIME,K)
push!(state_robustness_values, robustness_ecad)
push!(state_robustness_values, robustness_vim)

# calculate tgfb vegf vi robustness
# sample from 0 to 48 hrs and integrate
base_ = state[1]
pert_ = state[5]

base_ecad = base_[ti:tf,index_active_nfatc] #+ base_[ti:tf,index_active_nfatc_bcat]
pert_ecad = pert_[ti:tf,index_active_nfatc] #+ pert_[ti:tf,index_active_nfatc_bcat]
base_vim  = base_[ti:tf,index_psp1]
pert_vim  = pert_[ti:tf,index_psp1]

robustness_ecad = compute_robustness(base_ecad,pert_ecad,TIME,TIME,K)
robustness_vim  = compute_robustness(base_vim,pert_vim,TIME,TIME,K)
push!(state_robustness_values, robustness_ecad)
push!(state_robustness_values, robustness_vim)

# calculate tgfb vi robustness
# sample from 0 to 48 hrs and integrate
base_ = state[1]
pert_ = state[6]

base_ecad = base_[ti:tf,index_active_nfatc] #+ base_[ti:tf,index_active_nfatc_bcat]
pert_ecad = pert_[ti:tf,index_active_nfatc] #+ pert_[ti:tf,index_active_nfatc_bcat]
base_vim  = base_[ti:tf,index_psp1]
pert_vim  = pert_[ti:tf,index_psp1]

robustness_ecad = compute_robustness(base_ecad,pert_ecad,TIME,TIME,K)
robustness_vim  = compute_robustness(base_vim,pert_vim,TIME,TIME,K)
push!(state_robustness_values, robustness_ecad)
push!(state_robustness_values, robustness_vim)

# calculate vegf vi robustness
# sample from 0 to 48 hrs and integrate
base_ = state[1]
pert_ = state[7]
# base_ = state[2]# test
# pert_ = state[2]

base_ecad = base_[ti:tf,index_active_nfatc] #+ base_[ti:tf,index_active_nfatc_bcat]
pert_ecad = pert_[ti:tf,index_active_nfatc] #+ pert_[ti:tf,index_active_nfatc_bcat]
base_vim  = base_[ti:tf,index_psp1]
pert_vim  = pert_[ti:tf,index_psp1]

robustness_ecad = compute_robustness(base_ecad,pert_ecad,TIME,TIME,K)
robustness_vim  = compute_robustness(base_vim,pert_vim,TIME,TIME,K)
push!(state_robustness_values, robustness_ecad)
push!(state_robustness_values, robustness_vim)

# let's report the base case
push!(state_robustness_values, base_ecad[t48])
push!(state_robustness_values, base_vim[t48])

push!(robustness_array, state_robustness_values)

end

out_dir = "../poets/results_poets/"

writedlm(string(out_dir,"active_nfat_ap_sp1_p_training_set",".dat"),robustness_array,',')
# execfile('/home/dbassen/Dropbox/pyprgs/in_out_line.py')
# execfile('/home/dbassen/Dropbox/pyprgs/in_out_csv.py')
# b = lin.reader("results_poets/Robustness_training_set.dat")
# csv.nestedToFlt([( str.replace('Any[' ,'').replace(']','')).split(',') for str in b]) # python

### END ###

# # TGFB
# robustness_tgfb_ecad = []
# robustness_tgfb_vim = []
# # VEGF
# robustness_vegf_ecad = []
# robustness_vegf_vim = []
# # VEGF + TGFB
# robustness_tgfbvegf_ecad = []
# robustness_tgfbvegf_vim = []
# # VEGF + TGFB + VIVIT
# robustness_tgfbvegfviv_ecad = []
# robustness_tgfbvegfviv_vim = []
# # TGFB + VI
# robustness_tgfbviv_ecad = []
# robustness_tgfbviv_vim = []
# # VEGF + VI
# robustness_vegfviv_ecad = []
# robustness_vegfviv_vim = []


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Random

@test !("f=a=k=e=n=a=m=e" ∈ keys(ENV))

@testset "issue #10994" begin
    @test_throws ArgumentError ENV["bad\0name"] = "ok"
    @test_throws ArgumentError ENV["okname"] = "bad\0val"
    @test_throws ArgumentError Sys.set_process_title("bad\0title")

    withenv("bad"=>"dog") do
        @test_throws ArgumentError ENV["bad\0cat"]
    end
end
@testset "issue #11170" begin
    withenv("TEST"=>"nonempty") do
        @test ENV["TEST"] == "nonempty"
    end
    withenv("TEST"=>"") do
        @test ENV["TEST"] == ""
    end

    let c = collect(ENV)
        @test isa(c, Vector)
        @test length(ENV) == length(c)
        @test isempty(ENV) || first(ENV) in c
    end
end
@testset "non-existent keys" begin
    key = randstring(25)
    @test !haskey(ENV,key)
    @test_throws KeyError ENV[key]
    @test get(ENV,key,"default") == "default"
    @test get(() -> "default", ENV, key) == "default"

    key = randstring(25)
    @test !haskey(ENV, key)
    @test get!(ENV, key, "default") == "default"
    @test haskey(ENV, key)
    @test ENV[key] == "default"
end
@testset "#17956" begin
    @test length(ENV) > 1
    k1, k2 = "__TEST__", "__TEST1__"
    withenv(k1=>k1, k2=>k2) do
        b_k1, b_k2 = false, false
        for (k, v) in ENV
            if k==k1
                b_k1=true
            elseif k==k2
                b_k2=true
            end
        end
        @test b_k1 && b_k2
        io = IOBuffer()
        show(io, ENV)
        s = String(take!(io))
        @test occursin("$k1=$k1", s)
        @test occursin("$k2=$k2", s)

        @test pop!(ENV, k1) == k1
        @test !haskey(ENV, k1)
        ENV[k1] = k1
        @test pop!(ENV, k1) == k1
        @test pop!(ENV, k1, "not_there") == "not_there"

        ENV[k1] = k1
        @test delete!(ENV, k1) == ENV
        @test !haskey(ENV, k1)
    end
end
# Test for #10853
@test withenv(Dict{Any,Any}()...) do; true; end

# Test for #18141
for (k, v) in ENV
    if length(v) > 0
        @test v[end] != '\0'
    end
end

@testset "push" begin
    @test !haskey(ENV, "testing_envdict")
    push!(ENV, "testing_envdict" => "tested")
    @test haskey(ENV, "testing_envdict")
    @test ENV["testing_envdict"] == "tested"
    delete!(ENV, "testing_envdict")
end

if Sys.iswindows()
    @testset "windows case-insensitivity" begin
        for k in ("testing_envdict", "testing_envdict_\u00ee")
            K = uppercase(k)
            v = "tested $k"
            ENV[k] = v
            @test haskey(ENV, K)
            @test ENV[K] == v
            @test K in keys(ENV)
            @test K in collect(keys(ENV))
            @test k ∉ collect(keys(ENV))
            env = copy(ENV)
            @test haskey(env, K)
            @test env[K] == v
            @test !haskey(env, k)
            delete!(ENV, k)
            @test !haskey(ENV, K)
        end
    end
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
using Documenter, WhereTheWaterFlows

makedocs(;
    modules=[WhereTheWaterFlows],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/mauro3/WhereTheWaterFlows.jl/blob/{commit}{path}#L{line}",
    sitename="WhereTheWaterFlows.jl",
    authors="Mauro Werder <mauro3@runbox.com>",
    assets=String[],
)

deploydocs(;
    repo="github.com/mauro3/WhereTheWaterFlows.jl",
)


text:
module SparseDiffTools

using Compat
using FiniteDiff
using ForwardDiff
using LightGraphs
using LightGraphs: SimpleGraph
using Requires
using VertexSafeGraphs
using Adapt

using LinearAlgebra
using SparseArrays, ArrayInterface

using ForwardDiff: Dual, jacobian, partials, DEFAULT_CHUNK_THRESHOLD
using DataStructures: DisjointSets, find_root, union!

using ArrayInterface: matrix_colors

export  contract_color,
        greedy_d1,
        greedy_star1_coloring,
        greedy_star2_coloring,
        matrix2graph,
        matrix_colors,
        forwarddiff_color_jacobian!,
        forwarddiff_color_jacobian,
        ForwardColorJacCache,
        auto_jacvec,auto_jacvec!,
        num_jacvec,num_jacvec!,
        num_hesvec,num_hesvec!,
        numauto_hesvec,numauto_hesvec!,
        autonum_hesvec,autonum_hesvec!,
        num_hesvecgrad,num_hesvecgrad!,
        auto_hesvecgrad,auto_hesvecgrad!,
        JacVec,HesVec,HesVecGrad

include("coloring/high_level.jl")
include("coloring/backtracking_coloring.jl")
include("coloring/contraction_coloring.jl")
include("coloring/greedy_d1_coloring.jl")
include("coloring/acyclic_coloring.jl")
include("coloring/greedy_star1_coloring.jl")
include("coloring/greedy_star2_coloring.jl")
include("coloring/matrix2graph.jl")
include("differentiation/compute_jacobian_ad.jl")
include("differentiation/jaches_products.jl")

Base.@pure __parameterless_type(T) = Base.typename(T).wrapper
parameterless_type(x) = parameterless_type(typeof(x))
parameterless_type(x::Type) = __parameterless_type(x)

function __init__()
    @require Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" begin
        export numback_hesvec, numback_hesvec!, autoback_hesvec, autoback_hesvec!

        include("differentiation/jaches_products_zygote.jl")
    end
end

end # module


text:
###############################################################################
#
#  Descent methods for rank computations and local solubility of quartics
#
#
#  (C) 2022 Jeroen Hanselman
#
###############################################################################

export  quartic_local_solubility, R_soluble, Qp_soluble, quartic_rational_point_search, rank_2_torsion

###############################################################################
#
#  Quartic local solubility
#
###############################################################################

#TODO: Extend this to hyperelliptic curves
#TODO: Do this over number fields

@doc Markdown.doc"""
    quartic_local_solubility(a::fmpq, b::fmpq, c::fmpq, d::fmpq, e::fmpq) -> Bool

Check if the quartic defined by ax^4+bx^3+cx^2+dx+e is soluble over the real field $\mathbb{R}$ and over all local fields $\mathbb{Q}_p$.
"""
function quartic_local_solubility(a, b, c, d, e)

  if !(R_soluble(a,b,c,d,e))
    return false
  end

  if !(Qp_soluble(a,b,c,d,e,2))
    return false
  end

  R, x = PolynomialRing(QQ)

  delta = discriminant(a*x^4 + b*x^3 + c*x^2 +d*x + e)

  fac = factor(numerator(delta))
  p_list = [i[1] for i in fac]
  p_list = sort(p_list)

  for p in p_list
    if (!Qp_soluble(a,b,c,d,e,p))
      return false
    end
  end

  return true
end

@doc Markdown.doc"""
    R_soluble(a::fmpq, b::fmpq, c::fmpq, d::fmpq, e::fmpq) -> Bool

Check if the quartic defined by $ax^4+bx^3+cx^2+dx+e$ has a solution over $\mathbb{R}$.
"""
function R_soluble(a, b, c, d, e)
  if a>0
    return true
  end

  R,x = PolynomialRing(ZZ,"x")

  if (signature(a*x^4+b*x^3+c*x^2+d*x+e)[1]>0)
    return true
  else
    return false
  end
end

@doc Markdown.doc"""
    Qp_soluble(a::fmpq, b::fmpq, c::fmpq, d::fmpq, e::fmpq,
    p::fmpz) -> Bool

Check if the quartic defined by $ax^4+bx^3+cx^2+dx+e$ has a solution over the local field $\mathbb{Q}_p$.
"""
function Qp_soluble(a, b, c, d, e, p)
  R,x = PolynomialRing(QQ,"x")
  if(Zp_soluble(a,b,c,d,e,0,p,0))
    return true
  end

  if(Zp_soluble(e,d,c,b,a,0,p,1))
    return true
  end
  return false
end

function Zp_soluble(a,b,c,d,e,x_k,p,k)
  if (p==2)
    code = lemma7(a,b,c,d,e,x_k,k)
  else
    code = lemma6(a,b,c,d,e,x_k,p,k)
  end

  if code==1
    return true
  end
  if code==-1
    return false
  end

  for t in (0:p-1)
    if Zp_soluble(a,b,c,d,e,x_k+t*p^k,p,k+1)
      return true
    end
  end

  return false
end

#Z_p lifting for odd p

function lemma6(a,b,c,d,e,x,p,n)
  gx = a*x^4+b*x^3+c*x^2+d*x+e

  if(gx == 0)
    return 1
  end

  # Test if square in Zp
  l = valuation(gx,p)
  if iseven(l)
    if jacobi_symbol(numerator(gx//(p^l)),p)==1
      return 1
    end
  end

  gdx = 4*a*x^3+3*b*x^2+2*c*x+d

  if(gdx == 0)
   return 1
  end

  m = valuation(gdx,p)
  if (l>=m+n) && (n>m)
    return 1
  end
  if (l>=2*n) && (m>=n)
    return 0
  end
  return -1
end

#Z_p lifting for p=2

function lemma7(a,b,c,d,e,x,n)
  gx = a*x^4+b*x^3+c*x^2+d*x+e
  # Test if square in Z2
  if(gx == 0)
    return 1
  end

  l = valuation(gx,2)
  if divides(l,2)[1]
    if (mod(fmpz(numerator(gx//(2^l))),8) == 1)
      return 1
    end
  end

  gdx = 4*a*x^3+3*b*x^2+2*c*x+d

  # In this case the valuation of gdx is -infinity
  if(gdx == 0)
   return 1
  end

  m = valuation(gdx,2)
  gxodd = fmpz(numerator(gx//(2^l)))

  gxodd = mod(gxodd,4)

  if (l>=(m+n)) && (n>m)
    return 1
  end
  if (n>m) && (l==(m+n-1)) && iseven(l)
    return 1
  end
  if (n>m) && (l==(m+n-2)) && (gxodd==1) && iseven(l)
    return 1
  end
  if (m>=n) && (l>=2*n)
    return 0
  end
  if (m>=n) && (l==(2*n-2)) && (gxodd==1)
    return 0
  end

  return -1

end

#TODO: Make sieve-assisted
@doc Markdown.doc"""
    quartic_rational_point_search(a::fmpq, b::fmpq, c::fmpq, d::fmpq, e::fmpq,
    lower_bound::Int, upper_bound::Int) -> Bool

Check if the quartic defined by $ax^4+bx^3+cx^2+dx+e$ has a rational point $u/w$
where $\gcd(u,w) = 1$ with lower_bound <= u+w <= upper_bound.
"""
function quartic_rational_point_search(a, b, c, d, e, lower_bound, upper_bound, tmp1 = ZZ(), tmp2 = ZZ())
  R = Globals.Qx
  for n in lower_bound:upper_bound
    if n==1
      if AbstractAlgebra.issquare(a)
        return true
      end
      if AbstractAlgebra.issquare(e)
        return true
      end
    else
      @assert nbits(n) <= 14 # otherwise u^4 will overflow
      for u in (1:n-1)
        if gcd(u,n)==1
          w = n-u
          u2 = u^2
          u3 = u*u2
          u4 = (u2)^2
          w2 = w^2
          w3 = w*w2
          w4 = (w2)^2
          u3w = u3 * w
          u2w2 = u2 * w2
          uw3 = u * w3
          mul!(tmp1, a, u4)
          add!(tmp1, tmp1, b * u3w)
          mul!(tmp2, c, u2w2)
          add!(tmp1, tmp1, tmp2)
          add!(tmp1, tmp1, d * uw3)
          mul!(tmp2, e, w4)
          add!(tmp1, tmp1, tmp2)
          # we have tmp1 == a*u^4+b*u^3*w+c*u^2*w^2+d*u*w^3+e*w^4
          if AbstractAlgebra.issquare(tmp1)
            return true
          end
          add!(tmp1, tmp1, - 2 * b * u3w)
          add!(tmp1, tmp1, - 2 * d * uw3)
          # now tmp1 == a*u^4-b*u^3*w+c*u^2*w^2-d*u*w^3+e*w^4
          if AbstractAlgebra.issquare(tmp1)
            return true
          end
        end
      end
    end
  end
  return false
end

###############################################################################
#
#  Descent via 2-torsion
#
###############################################################################

# Following Cremona p. 87
@doc Markdown.doc"""
    rank_2_torsion(E::EllCrv, lim1::Int, lim2::Int) -> Int, Int, Int, Int

Compute bounds on rank and Sha using descent by 2-torsion isogeny.

lim1 gives a bound on the initial rational point search, lim2 a bound on
the exhaustive rational bound search in case of local solubility
The output consists of:
r_min: lower bound on the rank
r_max: upper bound on the rank
S: upper bound on #III(E)[phi]
S': upper bound on #III(E')[phi']

Here phi: E -> E' is an isogeny defined by a rational 2-torsion point on E
and phi' is its dual isogeny.
"""
function rank_2_torsion(E::EllCrv, lim1=100, lim2 = 1000)
  a1, a2, a3, a4, a6 = map(numerator,(a_invars(E)))
  if (a1==a3==0)
    s2 = a2
    s4 = a4
    s6 = a6
  else
    s2 = a1^2+4*a2
    s4 = 8*(a1*a3+2*a4)
    s6 = 16*(a3^2+4*a6)
  end

  R,x = PolynomialRing(QQ,"x")
  list = roots(x^3+s2*x^2+s4*x+s6)

  if all(!isintegral, list)
    throw(DomainError(E, "No rational 2-torsion"))
  end

  x0 = numerator(list[findfirst(isintegral, list)])

  c = 3*x0+s2
  d = (c+s2)*x0+s4
  _c = -2*c
  _d = c^2-4*d

  if iszero(d*_d)
    throw(DomainError(E, "Curve is singular"))
  end

  p_list = [p for (p,e) in factor(2*d*_d)]
  #print(p_list)
  n1,n2 = count_global_local(c,d,p_list,lim1,lim2)
  _n1,_n2 = count_global_local(_c,_d,p_list,lim1,lim2)

  e1 = log(2,n1)
  e2 = log(2,n2)
  _e1 = log(2,_n1)
  _e2 = log(2,_n2)

  r_min = Int(e1+_e1-2)
  r_max = Int(e2+_e2-2)
  S = Int(_n2//_n1)
  _S = Int(n2//n1)

  return r_min,r_max,S,_S
end

#TODO: Do better bookkeeping for divisors

function count_global_local(c, d, p_list, lim1, lim2)
  n1 = n2 = 1
  _d = c^2-4*d
  #print(d)
  d1_list = _squarefree_divisors(d)
  #print(d1_list)
  length(d1_list)
  tmp1 = ZZ()
  tmp2 = ZZ()
  for d1 in d1_list
    if quartic_rational_point_search(d1, 0, c, 0, divexact(d, d1), 1, lim1, tmp1, tmp2)
      #println(d1)
      n1 = n1+1
      n2 = n2+1
    else
      if everywhere_locally_soluble(c, d, _d, d1, p_list)
        #println(d1)
        n2 = n2+1
        if quartic_rational_point_search(d1, 0, c, 0, divexact(d, d1), lim1+1, lim2, tmp1, tmp2)
          n1 = n1+1
        end
      end
    end
  end
  return n1,n2
end

function everywhere_locally_soluble(c, d, _d, d1, p_list)
  if _d < 0 && d1 < 0
    return false
  end
  if _d > 0 && d1<0
    if isqrt(_d) < -c
      return false
    end
  end
  for p in p_list
    if !(Qp_soluble(d1, 0, c, 0, divexact(d, d1), p))
      return false
    end
  end
  return true
end

###############################################################################
#
#  Miscellaneous
#
###############################################################################

function _squarefree_divisors(n)
  divs = divisors(n)
  return filter(_issquarefree2, append!(-divs, divs))
end

function _issquarefree2(n)
  if (n==1)
    return false
  end
  return AbstractAlgebra.issquarefree(n)
end


text:
using Test
using CartesianGeneticProgramming

function test_inputs(f::Function, inps::AbstractArray)
    out = copy(f(inps...))
    @test typeof(out) == Float64
    @test all(out == f(inps...)) # functions are idempotent
    @test all(out .>= -1.0)
    @test all(out .<= 1.0)
    @test ~any(isnan.(out))
end

function test_functions(functions::Array{Function})
    for f in functions
        i1 = convert(Int16, 1)
        i2 = convert(Int16, 2)
        test_inputs(f, [[-1.0, -1.0], i1, i2])
        test_inputs(f, [[0.0, 0.0], i1, i2])
        test_inputs(f, [[1e-310, 1e-310], i1, i2])
        test_inputs(f, [[-1e-310, -2e-310], i1, i2])
        test_inputs(f, [[1.0, 1.0], i1, i2])
        for i in 1:5
            test_inputs(f, [[2 * rand() - 1, 2 * rand() - 1], i1, i2])
        end
    end
end

@testset "Mathematical functions" begin
    functions = [
        CGPFunctions.f_add,
        CGPFunctions.f_subtract,
        CGPFunctions.f_mult,
        CGPFunctions.f_div,
        CGPFunctions.f_abs,
        CGPFunctions.f_sqrt,
        CGPFunctions.f_pow,
        CGPFunctions.f_exp,
        CGPFunctions.f_sin,
        CGPFunctions.f_cos,
        CGPFunctions.f_tanh,
        CGPFunctions.f_sqrt_xy,
        CGPFunctions.f_lt,
        CGPFunctions.f_gt
    ]
    test_functions(functions)
end

@testset "Logical functions" begin
    functions = [
        CGPFunctions.f_and,
        CGPFunctions.f_or,
        CGPFunctions.f_xor,
        CGPFunctions.f_not
    ]
    test_functions(functions)
end

@testset "Range functions" begin
    functions = [
        CGPFunctions.f_avg,
        CGPFunctions.f_max,
        CGPFunctions.f_min,
        CGPFunctions.f_med
    ]
    test_functions(functions)
end


text:
using TemporalGPs: are_harmonised

function test_harmonise(a, b; recurse=true)
    h = harmonise(a, b)
    @test h isa Tuple
    @test length(h) == 2
    @test are_harmonised(h[1], h[2])

    recurse && test_harmonise(b, a; recurse=false)
    h′ = harmonise(b, a)
    @test h isa Tuple
    @test length(h) == 2
    @test are_harmonised(h′[1], h′[2])
    @test are_harmonised(h[1], h′[1])
    @test are_harmonised(h[1], h′[2])
end

@testset "harmonise" begin
    test_harmonise(5.0, 4.0)

    @testset "AbstractZero" begin
        test_harmonise(5.0, Zero())
        test_harmonise(Zero(), randn(10))
        test_harmonise(Zero(), Zero())
    end

    @testset "Array" begin
        test_harmonise(randn(5), randn(5))
        test_harmonise(
            [(randn(), randn()) for _ in 1:10],
            [Composite{Any}(randn(), rand()) for _ in 1:10],
        )
    end

    @testset "Tuple / Composite{Tuple}" begin
        test_harmonise((5, 4), (5, 4))
        test_harmonise(Composite{Tuple}(5, 4), (5, 4))
        test_harmonise(Composite{Tuple}(5, 4), Composite{Tuple}(5, 4))

        test_harmonise((5, Composite{Tuple}(randn(5))), (5, (randn(5), )))
        test_harmonise(
            Composite{Any}(Composite{Any}(randn(5))),
            (Composite{Any}(randn(5)), ),
        )
    end

    @testset "NamedTuple / Composite{NamedTuple}" begin
        test_harmonise(Composite{Any}(; m=4, P=5), Composite{Gaussian}(; m=5, P=4))
        test_harmonise(Composite{Any}(; m=4, P=5), Composite{Any}(; m=4))
        test_harmonise(Composite{Any}(; m=5), Composite{Any}(; P=4))

        test_harmonise(Composite{Any}(; m=(5, 4)), Composite{Any}(; P=4))

        test_harmonise(Composite{Any}(; m=5, P=4), Gaussian(5, 4))
        test_harmonise(Composite{Any}(; P=4), Gaussian(4, 5))
    end
end


text:
mutable struct OnlineSVGP{
    T<:Real,TLikelihood<:AbstractLikelihood,TInference<:AbstractInference,N
} <: AbstractGPModel{T,TLikelihood,TInference,N}
    f::NTuple{N,OnlineVarLatent{T}}
    likelihood::TLikelihood
    inference::TInference
    verbose::Int
    atfrequency::Int
    trained::Bool
end

"""
    OnlineSVGP(args...; kwargs...)

Online Sparse Variational Gaussian Process

## Arguments
- `kernel::Kernel` : Covariance function, can be any kernel from KernelFunctions.jl
- `likelihood` : Likelihood of the model. For compatibilities, see [`Likelihood Types`](@ref likelihood_user)
- `inference` : Inference for the model, see the [`Compatibility Table`](@ref compat_table))
- `Zalg` : Algorithm selecting how inducing points are selected

## Keywords arguments
- `verbose` : How much does the model print (0:nothing, 1:very basic, 2:medium, 3:everything)
- `optimiser` : Optimiser used for the kernel parameters. Should be an Optimiser object from the [Flux.jl](https://github.com/FluxML/Flux.jl) library, see list here [Optimisers](https://fluxml.ai/Flux.jl/stable/training/optimisers/) and on [this list](https://github.com/theogf/AugmentedGaussianProcesses.jl/tree/master/src/inference/optimisers.jl). Default is `ADAM(0.001)`
- `atfrequency::Int=1` : Choose how many variational parameters iterations are between hyperparameters optimization
- `mean=ZeroMean()` : PriorMean object, check the documentation on it [`MeanPrior`](@ref meanprior)
- `Zoptimiser` : Optimiser for inducing points locations
- `T::DataType=Float64` : Hint for what the type of the data is going to be.
"""
function OnlineSVGP(
    kernel::Kernel,
    likelihood::AbstractLikelihood,
    inference::AbstractInference,
    Zalg::InducingPoints.OnIPSA=OIPS(0.9);
    verbose::Integer=0,
    optimiser=ADAM(0.01),
    atfrequency::Integer=1,
    mean::Union{<:Real,AbstractVector{<:Real},PriorMean}=ZeroMean(),
    Zoptimiser=nothing,
    T::DataType=Float64,
)
    inference isa AnalyticVI || error("The inference object should be of type `AnalyticVI`")
    implemented(likelihood, inference) ||
        error("The ", likelihood, " is not compatible or implemented with the ", inference)

    if isa(optimiser, Bool)
        optimiser = optimiser ? ADAM(0.01) : nothing
    end

    if typeof(mean) <: Real
        mean = ConstantMean(mean)
    elseif typeof(mean) <: AbstractVector{<:Real}
        mean = EmpiricalMean(mean)
    end

    num_latent = n_latent(likelihood)
    latentf = ntuple(num_latent) do _
        return OnlineVarLatent(T, 0, [], Zalg, kernel, mean, optimiser, Zoptimiser)
    end
    # inference.nIter = 1
    return OnlineSVGP{T,typeof(likelihood),typeof(inference),num_latent}(
        latentf, likelihood, inference, verbose, atfrequency, false
    )
end

function Base.show(io::IO, model::OnlineSVGP) where {T}
    return print(
        io,
        "Online Variational Gaussian Process with a ",
        likelihood(model),
        " infered by ",
        inference(model),
    )
end

@traitimpl IsSparse{OnlineSVGP}

Zviews(m::OnlineSVGP) = Zview.(m.f)
objective(m::OnlineSVGP, state, y) = ELBO(m, state, y)


text:
"""
    LinearEquationOfState{FT} <: AbstractEquationOfState

Linear equation of state for seawater.
"""
struct LinearEquationOfState{FT} <: AbstractEquationOfState
    α :: FT
    β :: FT
end

"""
    LinearEquationOfState([FT=Float64;] α=1.67e-4, β=7.80e-4)

Returns parameters for a linear equation of state for seawater with
thermal expansion coefficient `α` [K⁻¹] and haline contraction coefficient
`β` [ppt⁻¹]. The buoyancy perturbation associated with a linear equation of state is

```math
    b = g (α T - β S)
```

Default constants are taken from Table 1.2 (page 33) of Vallis, "Atmospheric and Oceanic Fluid
Dynamics: Fundamentals and Large-Scale Circulation" (2ed, 2017).
"""
LinearEquationOfState(FT=Float64; α=1.67e-4, β=7.80e-4) = LinearEquationOfState{FT}(α, β)

@inline  thermal_expansion(Θ, sᴬ, D, eos::LinearEquationOfState) = eos.α
@inline haline_contraction(Θ, sᴬ, D, eos::LinearEquationOfState) = eos.β

# Shortcuts
@inline  thermal_expansionᶜᶜᶜ(i, j, k, grid, eos::LinearEquationOfState, C) = eos.α
@inline  thermal_expansionᶠᶜᶜ(i, j, k, grid, eos::LinearEquationOfState, C) = eos.α
@inline  thermal_expansionᶜᶠᶜ(i, j, k, grid, eos::LinearEquationOfState, C) = eos.α
@inline  thermal_expansionᶜᶜᶠ(i, j, k, grid, eos::LinearEquationOfState, C) = eos.α

@inline haline_contractionᶜᶜᶜ(i, j, k, grid, eos::LinearEquationOfState, C) = eos.β
@inline haline_contractionᶠᶜᶜ(i, j, k, grid, eos::LinearEquationOfState, C) = eos.β
@inline haline_contractionᶜᶠᶜ(i, j, k, grid, eos::LinearEquationOfState, C) = eos.β
@inline haline_contractionᶜᶜᶠ(i, j, k, grid, eos::LinearEquationOfState, C) = eos.β

const LinearSeawaterBuoyancy = SeawaterBuoyancy{FT, <:LinearEquationOfState} where FT
const LinearTemperatureSeawaterBuoyancy = SeawaterBuoyancy{FT, <:LinearEquationOfState, <:Nothing, <:Number} where FT
const LinearSalinitySeawaterBuoyancy = SeawaterBuoyancy{FT, <:LinearEquationOfState, <:Number, <:Nothing} where FT

@inline buoyancy_perturbation(i, j, k, grid, b::LinearSeawaterBuoyancy, C) =
    @inbounds b.gravitational_acceleration * (  b.equation_of_state.α * C.T[i, j, k]
                                              - b.equation_of_state.β * C.S[i, j, k])

@inline buoyancy_perturbation(i, j, k, grid, b::LinearTemperatureSeawaterBuoyancy, C) =
    @inbounds b.gravitational_acceleration * b.equation_of_state.α * C.T[i, j, k]

@inline buoyancy_perturbation(i, j, k, grid, b::LinearSalinitySeawaterBuoyancy, C) =
    @inbounds - b.gravitational_acceleration * b.equation_of_state.β * C.S[i, j, k]


text:
using MultilayerPerceptron
using Random: seed!, shuffle

using Makie

# Ref: scikit-learn's make_moons
function make_moons(n_samples::Int; noise=0.1)::Tuple{Matrix{Float64}, BitArray}
    seed!()
    n_concave = n_samples ÷ 2
    n_convex = n_samples - n_concave
    🌛 = begin
        angles = LinRange(0, π, n_concave)
        hcat(cos.(angles), sin.(angles))'
    end
    🌜 = begin
        angles = LinRange(0, π, n_convex)
        hcat(1.0 .- cos.(angles), (1.0 - 0.5) .- sin.(angles))'
    end
    X = hcat(🌛, 🌜)
    Y = vcat(falses(n_concave), trues(n_convex))
    X .+= randn(size(X)) * noise

    rnd_ind = shuffle(1:n_samples)
    X = X[:, rnd_ind]
    Y = reshape(Y[rnd_ind], 1, :)

    (X, Y)
end

function moons_example()
    (train_x, train_y) = make_moons(250)

    N = 500
    hidden_layer1 = Layer((3, 2); activation_fn=logistic, grad_fn=∇logistic)
    hidden_layer2 = Layer((3, 3); activation_fn=logistic, grad_fn=∇logistic)
    output_layer = Layer((1, 3))

    mlp = MLP([hidden_layer1, hidden_layer2, output_layer])
    fit!(mlp, train_x, train_y;
        loss_fn=logit_binary_cross_entropy,
        learning_rate=0.2, batch_size=5, n_epochs=N)

    (test_x, test_y) = make_moons(250)
    pred_y = predict(mlp, test_x) .> 0.5

    test_y, pred_y = vec(test_y), vec(pred_y)

    AbstractPlotting.__init__()
    scene = Scene()
    tp, fn, tn, fp =
        @. test_y & pred_y, test_y & (!pred_y), (!test_y) & (!pred_y), (!test_y) & pred_y
    any(tp) && scatter!(scene, test_x[:, tp]'; color=:red)
    any(tn) && scatter!(scene, test_x[:, tn]'; color=:green)
    any(fn) && scatter!(scene, test_x[:, fn]'; color=:red4)
    any(fp) && scatter!(scene, test_x[:, fp]'; color=:green1)
    display(scene)
    println("TP = $(sum(tp)), FN = $(sum(fn)); TN = $(sum(tn)), FP = $(sum(fp))")
end

moons_example()

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
# This file is part of the TaylorSeries.jl Julia package, MIT license
#
# Parameters for HomogeneousPolynomial and TaylorN

"""
    ParamsTaylor1

DataType holding the current variable name for `Taylor1`.

**Field:**

- `var_name   :: String`  Names of the variables

These parameters can be changed using [`set_taylor1_varname`](@ref)
"""
mutable struct ParamsTaylor1
    var_name   :: String
end

const _params_Taylor1_ = ParamsTaylor1("t")

"""
    set_taylor1_varname(var::String)

Change the displayed variable for `Taylor1` objects.
"""
set_taylor1_varname(var::String) = _params_Taylor1_.var_name = strip(var)

"""
    ParamsTaylorN

DataType holding the current parameters for `TaylorN` and
`HomogeneousPolynomial`.

**Fields:**

- `order            :: Int`  Order (degree) of the polynomials
- `num_vars         :: Int`  Number of variables
- `variable_names   :: Vector{String}`  Names of the variables
- `variable_symbols :: Vector{Symbol}`  Symbols of the variables

These parameters can be changed using [`set_variables`](@ref)
"""
mutable struct ParamsTaylorN
    order            :: Int
    num_vars         :: Int
    variable_names   :: Vector{String}
    variable_symbols :: Vector{Symbol}
end

ParamsTaylorN(order, num_vars, variable_names) = ParamsTaylorN(order, num_vars, variable_names, Symbol.(variable_names))

const _params_TaylorN_ = ParamsTaylorN(6, 2, ["x₁", "x₂"])

## Utilities to get the maximum order, number of variables, their names and symbols
get_order() = _params_TaylorN_.order
get_numvars() = _params_TaylorN_.num_vars
get_variable_names() = _params_TaylorN_.variable_names
get_variable_symbols() = _params_TaylorN_.variable_symbols
function lookupvar(s::Symbol)
    ind = findfirst(x -> x==s, _params_TaylorN_.variable_symbols)
    isa(ind, Nothing) && return 0
    return ind
end

"""
    get_variables(T::Type, [order::Int=get_order()])

Return a `TaylorN{T}` vector with each entry representing an
independent variable. It takes the default `_params_TaylorN_` values
if `set_variables` hasn't been changed with the exception that `order`
can be explicitely established by the user without changing internal values
for `num_vars` or `variable_names`. Ommiting `T` defaults to `Float64`.
"""
get_variables(T::Type, order::Int=get_order()) =
    [TaylorN(T, i, order=order) for i in 1:get_numvars()]
get_variables(order::Int=get_order()) =
    [TaylorN(Float64, i, order=order) for i in 1:get_numvars()]

"""
    set_variables([T::Type], names::String; [order=get_order(), numvars=-1])

Return a `TaylorN{T}` vector with each entry representing an
independent variable. `names` defines the output for each variable
(separated by a space). The default type `T` is `Float64`,
and the default for `order` is the one defined globally.
Changing the `order` or `numvars` resets the hash_tables.

If `numvars` is not specified, it is inferred from `names`. If only
one variable name is defined and `numvars>1`, it uses this name with
subscripts for the different variables.

```julia
julia> set_variables(Int, "x y z", order=4)
3-element Array{TaylorSeries.TaylorN{Int},1}:
  1 x + 𝒪(‖x‖⁵)
  1 y + 𝒪(‖x‖⁵)
  1 z + 𝒪(‖x‖⁵)

julia> set_variables("α", numvars=2)
2-element Array{TaylorSeries.TaylorN{Float64},1}:
  1.0 α₁ + 𝒪(‖x‖⁵)
  1.0 α₂ + 𝒪(‖x‖⁵)

julia> set_variables("x", order=6, numvars=2)
2-element Array{TaylorSeries.TaylorN{Float64},1}:
  1.0 x₁ + 𝒪(‖x‖⁷)
  1.0 x₂ + 𝒪(‖x‖⁷)
```
"""
function set_variables(R::Type, names::Vector{T}; order=get_order()) where
        {T<:AbstractString}

    order ≥ 1 || error("Order must be at least 1")

    num_vars = length(names)
    num_vars ≥ 1 || error("Number of variables must be at least 1")

    _params_TaylorN_.variable_names = names
    _params_TaylorN_.variable_symbols = Symbol.(names)

    if !(order == get_order() && num_vars == get_numvars())
        # if these are unchanged, no need to regenerate tables

        _params_TaylorN_.order = order
        _params_TaylorN_.num_vars = num_vars

        resize!(coeff_table,order+1)
        resize!(index_table,order+1)
        resize!(size_table,order+1)
        resize!(pos_table,order+1)

        coeff_table[:], index_table[:], size_table[:], pos_table[:] =
            generate_tables(num_vars, order)
        GC.gc();
    end

    # return a list of the new variables
    TaylorN{R}[TaylorN(R,i) for i in 1:get_numvars()]
end
set_variables(R::Type, symbs::Vector{T}; order=get_order()) where
    {T<:Symbol} = set_variables(R, string.(symbs), order=order)

set_variables(names::Vector{T}; order=get_order()) where {T<:AbstractString} =
    set_variables(Float64, names, order=order)
set_variables(symbs::Vector{T}; order=get_order()) where {T<:Symbol} =
    set_variables(Float64, symbs, order=order)

function set_variables(R::Type, names::T; order=get_order(), numvars=-1) where
        {T<:AbstractString}

    variable_names = split(names)

    if length(variable_names) == 1 && numvars ≥ 1
        name = variable_names[1]
        variable_names = [string(name, subscriptify(i)) for i in 1:numvars]
    end

    set_variables(R, variable_names, order=order)
end
set_variables(R::Type, symbs::Symbol; order=get_order(), numvars=-1) =
    set_variables(R, string(symbs), order=order, numvars=numvars)

set_variables(names::T; order=get_order(), numvars=-1) where {T<:AbstractString} =
    set_variables(Float64, names, order=order, numvars=numvars)
set_variables(symbs::Symbol; order=get_order(), numvars=-1) =
    set_variables(Float64, string(symbs), order=order, numvars=numvars)

"""
    show_params_TaylorN()

Display the current parameters for `TaylorN` and `HomogeneousPolynomial` types.
"""
function show_params_TaylorN()
    @info( """
    Parameters for `TaylorN` and `HomogeneousPolynomial`:
    Maximum order       = $(get_order())
    Number of variables = $(get_numvars())
    Variable names      = $(get_variable_names())
    Variable symbols    = $(Symbol.(get_variable_names()))
    """)
    nothing
end

# Control the display of the big 𝒪 notation
const bigOnotation = Bool[true]
const _show_default = [false]

"""
    displayBigO(d::Bool) --> nothing

Set/unset displaying of the big 𝒪 notation in  the output
of `Taylor1` and `TaylorN` polynomials. The initial value is
`true`.
"""
displayBigO(d::Bool) = (bigOnotation[end] = d; d)

"""
    use_Base_show(d::Bool) --> nothing

Use `Base.show_default` method (default `show` method
in Base), or a custom display. The initial value is
`false`, so customized display is used.
"""
use_show_default(d::Bool) = (_show_default[end] = d; d)