text:
using Test
using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit.BipartiteGraphs
using ModelingToolkit.StructuralTransformations

function check_contract(sys)
    graph = ModelingToolkit.get_tearing_state(sys).structure.graph
    sys = tearing_substitution(sys)
    state = TearingState(sys)
    fullvars = state.fullvars

    eqs = equations(sys)
    var2idx = Dict(enumerate(fullvars))
    for (i, eq) in enumerate(eqs)
        actual = union(ModelingToolkit.vars(eq.lhs), ModelingToolkit.vars(eq.rhs))
        actual = filter(!ModelingToolkit.isparameter, collect(actual))
        current = Set(fullvars[𝑠neighbors(graph, i)])
        @test isempty(setdiff(actual, current))
    end
end

function check_rc_sol(sol)
    @test sol[rc_model.resistor.p.i] == sol[resistor.p.i] == sol[capacitor.p.i]
    @test sol[rc_model.resistor.n.i] == sol[resistor.n.i] == -sol[capacitor.p.i]
    @test sol[rc_model.capacitor.n.i] == sol[capacitor.n.i] == -sol[capacitor.p.i]
    @test iszero(sol[rc_model.ground.g.i])
    @test iszero(sol[rc_model.ground.g.v])
    @test sol[rc_model.resistor.v] == sol[resistor.v] ==
          sol[source.p.v] - sol[capacitor.p.v]
end

include("../examples/rc_model.jl")

@test ModelingToolkit.n_extra_equations(capacitor) == 2
@test length(equations(structural_simplify(rc_model, allow_parameter = false))) > 1
sys = structural_simplify(rc_model)
check_contract(sys)
@test !isempty(ModelingToolkit.defaults(sys))
u0 = [capacitor.v => 0.0
      capacitor.p.i => 0.0
      resistor.v => 0.0]
prob = ODEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Rodas4())
check_rc_sol(sol)
prob = ODAEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Rodas4())
check_rc_sol(sol)

let
    # 1478
    @named resistor2 = Resistor(R = R)
    rc_eqs2 = [connect(source.p, resistor.p)
               connect(resistor.n, resistor2.p)
               connect(resistor2.n, capacitor.p)
               connect(capacitor.n, source.n)
               connect(capacitor.n, ground.g)]

    @named _rc_model2 = ODESystem(rc_eqs2, t)
    @named rc_model2 = compose(_rc_model2,
                               [resistor, resistor2, capacitor, source, ground])
    sys2 = structural_simplify(rc_model2)
    prob2 = ODAEProblem(sys2, u0, (0, 10.0))
    sol2 = solve(prob2, Tsit5())
    @test sol2[source.p.i] == sol2[rc_model2.source.p.i] == -sol2[capacitor.i]
end

# Outer/inner connections
function rc_component(; name)
    R = 1
    C = 1
    @named p = Pin()
    @named n = Pin()
    @named resistor = Resistor(R = R)
    @named capacitor = Capacitor(C = C)
    eqs = [connect(p, resistor.p);
           connect(resistor.n, capacitor.p);
           connect(capacitor.n, n)]
    @named sys = ODESystem(eqs, t)
    compose(sys, [p, n, resistor, capacitor]; name = name)
end

@named ground = Ground()
@named source = ConstantVoltage(V = 1)
@named rc_comp = rc_component()
eqs = [connect(source.p, rc_comp.p)
       connect(source.n, rc_comp.n)
       connect(source.n, ground.g)]
@named sys′ = ODESystem(eqs, t)
@named sys_inner_outer = compose(sys′, [ground, source, rc_comp])
@test_nowarn show(IOBuffer(), MIME"text/plain"(), sys_inner_outer)
expand_connections(sys_inner_outer, debug = true)
sys_inner_outer = structural_simplify(sys_inner_outer)
@test !isempty(ModelingToolkit.defaults(sys_inner_outer))
u0 = [rc_comp.capacitor.v => 0.0
      rc_comp.capacitor.p.i => 0.0
      rc_comp.resistor.v => 0.0]
prob = ODEProblem(sys_inner_outer, u0, (0, 10.0), sparse = true)
sol_inner_outer = solve(prob, Rodas4())
@test sol[capacitor.v] ≈ sol_inner_outer[rc_comp.capacitor.v]

u0 = [
    capacitor.v => 0.0,
]
prob = ODAEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Tsit5())

@test sol[resistor.p.i] == sol[capacitor.p.i]
@test sol[resistor.n.i] == -sol[capacitor.p.i]
@test sol[capacitor.n.i] == -sol[capacitor.p.i]
@test iszero(sol[ground.g.i])
@test iszero(sol[ground.g.v])
@test sol[resistor.v] == sol[source.p.v] - sol[capacitor.p.v]
#using Plots
#plot(sol)

include("../examples/serial_inductor.jl")
sys = structural_simplify(ll_model)
check_contract(sys)
u0 = [inductor1.i => 0.0
      inductor2.i => 0.0
      inductor2.v => 0.0]
@test_throws Any ODEProblem(sys, u0, (0, 10.0))
@test_throws Any ODAEProblem(sys, u0, (0, 10.0))
prob = DAEProblem(sys, Differential(t).(states(sys)) .=> 0, u0, (0, 0.5))
@test_nowarn sol = solve(prob, DFBDF())

@variables t x1(t) x2(t) x3(t) x4(t)
D = Differential(t)
@named sys1_inner = ODESystem([D(x1) ~ x1], t)
@named sys1_partial = compose(ODESystem([D(x2) ~ x2], t; name = :foo), sys1_inner)
@named sys1 = extend(ODESystem([D(x3) ~ x3], t; name = :foo), sys1_partial)
@named sys2 = compose(ODESystem([D(x4) ~ x4], t; name = :foo), sys1)
@test_nowarn sys2.sys1.sys1_inner.x1 # test the correct nesting

# compose tests
@parameters t

function record_fun(; name)
    pars = @parameters a=10 b=100
    ODESystem(Equation[], t, [], pars; name)
end

function first_model(; name)
    @named foo = record_fun()

    defs = Dict()
    defs[foo.a] = 3
    defs[foo.b] = 300
    pars = @parameters x=2 y=20
    compose(ODESystem(Equation[], t, [], pars; name, defaults = defs), foo)
end
@named goo = first_model()
@unpack foo = goo
@test ModelingToolkit.defaults(goo)[foo.a] == 3
@test ModelingToolkit.defaults(goo)[foo.b] == 300

#=
model Circuit
  Ground ground;
  Load load;
  Resistor resistor;
equation
  connect(load.p , ground.p);
  connect(resistor.p, ground.p);
end Circuit;
model Load
  extends TwoPin;
  Resistor resistor;
equation
  connect(p, resistor.p);
  connect(resistor.n, n);
end Load;
=#

function Load(; name)
    R = 1
    @named p = Pin()
    @named n = Pin()
    @named resistor = Resistor(R = R)
    eqs = [connect(p, resistor.p);
           connect(resistor.n, n)]
    @named sys = ODESystem(eqs, t)
    compose(sys, [p, n, resistor]; name = name)
end

function Circuit(; name)
    R = 1
    @named ground = Ground()
    @named load = Load()
    @named resistor = Resistor(R = R)
    eqs = [connect(load.p, ground.g);
           connect(resistor.p, ground.g)]
    @named sys = ODESystem(eqs, t)
    compose(sys, [ground, resistor, load]; name = name)
end

@named foo = Circuit()
@test structural_simplify(foo) isa ModelingToolkit.AbstractSystem

# BLT tests
using LinearAlgebra
function parallel_rc_model(i; name, source, ground, R, C)
    resistor = HeatingResistor(name = Symbol(:resistor, i), R = R)
    capacitor = Capacitor(name = Symbol(:capacitor, i), C = C)
    heat_capacitor = HeatCapacitor(name = Symbol(:heat_capacitor, i))

    rc_eqs = [connect(source.p, resistor.p)
              connect(resistor.n, capacitor.p)
              connect(capacitor.n, source.n, ground.g)
              connect(resistor.h, heat_capacitor.h)]

    compose(ODESystem(rc_eqs, t, name = Symbol(name, i)),
            [resistor, capacitor, source, ground, heat_capacitor])
end
V = 2.0
@named source = ConstantVoltage(V = V)
@named ground = Ground()
N = 50
Rs = 10 .^ range(0, stop = -4, length = N)
Cs = 10 .^ range(-3, stop = 0, length = N)
rc_systems = map(1:N) do i
    parallel_rc_model(i; name = :rc, source = source, ground = ground, R = Rs[i], C = Cs[i])
end;
@variables E(t) = 0.0
eqs = [
    D(E) ~ sum(((i, sys),) -> getproperty(sys, Symbol(:resistor, i)).h.Q_flow,
               enumerate(rc_systems)),
]
@named _big_rc = ODESystem(eqs, t, [E], [])
@named big_rc = compose(_big_rc, rc_systems)
ts = TearingState(expand_connections(big_rc))
@test istriu(but_ordered_incidence(ts)[1])


text:
export PosteriorFactor

"""
A `PosteriorFactor` specifies the subset of variables that comprise
a joint factor in the posterior factorization. A PosteriorFactor can
be defined by a (selection of) variable(s), or on the complete graph.
"""
mutable struct PosteriorFactor
    id::Symbol
    internal_edges::Set{Edge}

    # Fields set by algorithm construction
    target_variables::Set{Variable} # Target variables for which marginals are required
    target_clusters::Set{Cluster} # Target clusters for which marginals are required
    breaker_interfaces::Set{Interface} # Target interfaces for which messages are required
    ep_sites::Set{Interface} # Sites for expectation propagation update

    # Fields set by algorithm assembler
    algorithm_id::Symbol # Specify the algorithm id for this posterior_factor
    schedule::Schedule # Specify the internal message passing schedule for this posterior factor
    marginal_table::MarginalTable # Specify the marginal updates for internal variables
    optimize::Bool # Indicate the need for an optimization block
    initialize::Bool # Indicate the need for a message initialization block

    # Contruct a posterior factor that encompasses all variables in the graph
    function PosteriorFactor(fg::FactorGraph; pfz=currentPosteriorFactorization(), id=generateId(PosteriorFactor))
        internal_edges = edges(fg) # Include all edges in a single posterior factor
        self = new(id, internal_edges, Set{Variable}(), Set{Cluster}(), Set{Interface}(), Set{Interface}()) # Initialize posterior factor
        pfz.posterior_factors[id] = self # Register self with the algorithm

        return self
    end

    # Construct a posterior factor that includes all variables deterministically linked to the target variables
    function PosteriorFactor(seed_variables::Set{Variable}; pfz=currentPosteriorFactorization(), id=generateId(PosteriorFactor))
        internal_edges = extend(edges(seed_variables)) # Include all deterministically liked variables in a single posterior factor
        self = new(id, internal_edges, Set{Variable}(), Set{Cluster}(), Set{Interface}(), Set{Interface}()) # Initialize posterior factor
        pfz.posterior_factors[id] = self # Register self with the posterior factorization

        return self
    end
end

PosteriorFactor(seed_variable::Variable; pfz=currentPosteriorFactorization(), id=generateId(PosteriorFactor)) = PosteriorFactor(Set([seed_variable]), pfz=pfz, id=id)
PosteriorFactor(seed_variables::Vector{Variable}; pfz=currentPosteriorFactorization(), id=generateId(PosteriorFactor)) = PosteriorFactor(Set(seed_variables), pfz=pfz, id=id)

"""
`messagePassingSchedule()` generates a message passing schedule for the posterior factor
"""
function messagePassingSchedule(pf::PosteriorFactor)
    nodes_connected_to_external_edges = nodesConnectedToExternalEdges(pf)

    # Schedule messages towards posterior distributions and target sites, limited to the internal edges
    breaker_interfaces = sort(collect(pf.breaker_interfaces), rev=true)
    schedule = summaryPropagationSchedule(sort(collect(pf.target_variables), rev=true),
                                          sort(collect(pf.target_clusters), rev=true),
                                          limit_set=pf.internal_edges,
                                          target_sites=breaker_interfaces)
    for entry in schedule
        if entry.interface in pf.ep_sites
            entry.message_update_rule = ExpectationPropagationRule{typeof(entry.interface.node)}
        elseif (entry.interface.node in nodes_connected_to_external_edges) && !isa(entry.interface.node, DeltaFactor)
            local_posterior_factors = localPosteriorFactors(entry.interface.node)
            if allunique(local_posterior_factors) # Local posterior factorization is naive
                entry.message_update_rule = NaiveVariationalRule{typeof(entry.interface.node)}
            else
                entry.message_update_rule = StructuredVariationalRule{typeof(entry.interface.node)}
            end
        else
            entry.message_update_rule = SumProductRule{typeof(entry.interface.node)}
        end
    end

    breaker_types = breakerTypes(collect(pf.breaker_interfaces))
    inferUpdateRules!(schedule, inferred_outbound_types=breaker_types)

    return schedule
end

function draw(pf::PosteriorFactor; schedule=ScheduleEntry[], args...)
    subgraph_nodes = nodes(pf.internal_edges)
    external_edges = setdiff(edges(subgraph_nodes), pf.internal_edges)
    ForneyLab.graphviz(ForneyLab.genDot(subgraph_nodes, pf.internal_edges, schedule=schedule, external_edges=external_edges); args...)
end

"""
Find the smallest legal subgraph that includes the argument edges. Default setting terminates the search at soft factors
and does not constrain the search to a limiting set (as specified by an empty `limit_set` argument).
"""
function extend(edge_set::Set{Edge}; terminate_at_soft_factors=true, limit_set=Set{Edge}())
    extension = Set{Edge}() # Initialize extension
    stack = copy(edge_set) # Initialize stack
    while !isempty(stack) # As long as there are unchecked edges connected through deterministic nodes
        current_edge = pop!(stack) # Pick one
        push!(extension, current_edge) # Add to edge extension

        connected_nodes = [] # Find nodes connected to edge (as a vector)
        (current_edge.a == nothing) || push!(connected_nodes, current_edge.a.node)
        (current_edge.b == nothing) || push!(connected_nodes, current_edge.b.node)

        for node in connected_nodes # Check both head and tail node (if present)
            if !terminate_at_soft_factors || isa(node, DeltaFactor)
                for interface in node.interfaces
                    if (interface.edge !== current_edge) && !(interface.edge in extension) && ( isempty(limit_set) || (interface.edge in limit_set) ) # No backtracking, if edge is not already visited and edge is contained within limit set
                        push!(stack, interface.edge) # Add unseen edges to the stack (to visit sometime in the future)
                    end
                end
            end
        end
    end

    return extension
end

extend(edge::Edge; terminate_at_soft_factors=true, limit_set=Set{Edge}()) = extend(Set{Edge}([edge]), terminate_at_soft_factors=terminate_at_soft_factors, limit_set=limit_set)

"""
Find the `PosteriorFactor` that `edge` belongs to (if available)
"""
function posteriorFactor(edge::Edge)
    dict = current_posterior_factorization.edge_to_posterior_factor
    if haskey(dict, edge)
        pf = dict[edge]
    else # No posterior factor is found or edge is clamped; return the edge itself
        pf = edge
    end

    return pf::Union{PosteriorFactor, Edge}
end

"""
Check if edge or interface is terminated by a Clamp node
"""
isClamped(interface::Interface) = isdefined(interface, :node) && isa(interface.node, Clamp)
isClamped(edge::Edge) = isClamped(edge.a) || isClamped(edge.b)
isClamped(::Nothing) = false

"""
Check if node imposes a point-mass constraint. This enables external overloading
of what it means to be a point-mass constrained.
"""
isDeltaConstraint(node::FactorNode) = false
isDeltaConstraint(node::Clamp) = true #

"""
Check if edge or interface is terminated by a point-mass constraint
"""
isDeltaConstrained(interface::Interface) = isdefined(interface, :node) && isDeltaConstraint(interface.node)
isDeltaConstrained(edge::Edge) = isDeltaConstrained(edge.a) || isDeltaConstrained(edge.b)
isDeltaConstrained(::Nothing) = false

"""
Return the ids of the posterior factors to which edges connected to `node` belong
"""
localPosteriorFactors(node::FactorNode) = Any[posteriorFactor(interface.edge) for interface in node.interfaces]

"""
Find the nodes in `posterior_factor` that are connected to external edges
"""
function nodesConnectedToExternalEdges(pf::PosteriorFactor)
    subgraph_nodes = nodes(pf.internal_edges)
    external_edges = setdiff(edges(subgraph_nodes), pf.internal_edges)

    # nodes_connected_to_external_edges are the nodes connected to external edges that are also connected to internal edges
    nodes_connected_to_external_edges = intersect(nodes(external_edges), subgraph_nodes)

    return nodes_connected_to_external_edges
end

"""
Return a dictionary from posterior factors to regions local to `node`
"""
function localEdgeToRegion(node::FactorNode)
    # For each edge connected to node, collect the respective posterior factors and regions
    local_edge_to_region = Dict{Edge, Region}()
    for interface in node.interfaces
        local_edge_to_region[interface.edge] = region(node, interface.edge)
    end

    return local_edge_to_region
end

function collectAverageEnergyInbounds(node::FactorNode)
    inbounds = Any[]
    local_edge_to_region = localEdgeToRegion(node)

    encountered_regions = Region[] # Keep track of encountered regions
    for node_interface in node.interfaces
        inbound_interface = ultimatePartner(node_interface)
        current_region = region(node_interface.node, node_interface.edge)

        if isClamped(inbound_interface)
            # Hard-code marginal of constant node in schedule
            push!(inbounds, assembleClamp!(copy(inbound_interface.node), ProbabilityDistribution)) # Copy Clamp before assembly to prevent overwriting dist_or_msg field
        elseif !(current_region in encountered_regions)
            # Collect marginal entry from marginal dictionary (if marginal entry is not already accepted)
            target = local_edge_to_region[node_interface.edge]
            push!(inbounds, current_inference_algorithm.target_to_marginal_entry[target])
        end

        push!(encountered_regions, current_region)
    end

    return inbounds
end

ext: jl


text:
# Used to check up on the performance of non-allocating rod functions,
# and quickly check that they only allocate views, if anything.

using DiscreteElasticRods
using BenchmarkTools
using Random

N = 500
suite = BenchmarkGroup()
Random.seed!(42)

# Linear ALgebra Sections
a = rand(N,3)
b = rand(N,3)
c = rand(N,3)

suite["cross3!"] = @benchmarkable DER.cross3!($c,$a,$b)

d = rand(N,1)
suite["dot3!"] = @benchmarkable DER.dot3!($d,$a,$b) samples = 50000
suite["norm3!(a,b)"] = @benchmarkable DER.norm3!($d, $a)
suite["norm3(a)!"] = @benchmarkable DER.norm3!($a)

e = rand(N-1,3)
suite["edges!"] = @benchmarkable DER.edges!($e,$b) samples = 50000
suite["tangents!"] = @benchmarkable DER.tangents!($e,$b) samples = 50000

# Parallel Transport and Rotation
v1 = rand(N,3)
v2 = similar(v1)
t1 = rand(N,3)
t2 = rand(N,3)
DER.norm3!(t1)
DER.norm3!(t2)
DER.norm3!(v1)

cache = zeros(N,6)
suite["ptransport!"] = @benchmarkable DER.ptransport!($v2,$v1,$t1,$t2,$cache)

θ = rand(N,1)
suite["rotate_orthogonal_unit!"] = @benchmarkable DER.rotate_orthogonal_unit!($v2,$t2,$θ,$cache) setup=(v2 .= $v1)

# Update and Kinematics
x0 = rand(N,3)
d0 = rand(N-1,3)
DER.norm3!(d0)
x1 = zeros(N,3)
d1 = zeros(N-1,3)
Δx = rand(N,3)
Δθ = rand(N-1,1)
cache2 = DER.allocate_cache(DER.rod_update!, Float64, N-2)
suite["rod_update!"] =
    @benchmarkable DER.rod_update!($x1,$d1,$x0,$d0,$Δx,$Δθ,$cache2)

l = zeros(N-1)
κ = zeros(N-2,2)
τ = zeros(N-2)
cache3 = DER.allocate_cache(DER.full_kinematics!, Float64, N-2)
suite["full_kinematics!"] =
    @benchmarkable DER.full_kinematics!($l,$κ,$τ,$x0,$d0,$cache3)

l0 = zeros(N-1)
κ0 = zeros(N-2,2)
τ0 = zeros(N-2)
k = rand()
B = rand(1,2)
β = rand()
suite["elastic_energy"] =
    @benchmarkable DER.elastic_energy($l,$κ,$τ,$l0,$κ0,$τ0,$k,$B,$β,$cache3[2])

m = rand(N-1)
g = 10.0
suite["gravitational_energy"] =
    @benchmarkable DER.gravitational_energy($x0,$m,$g,$cache3[1])

# Set: tune, run, and save a benchmark
# tune!(suite)
# BenchmarkTools.save("benchmark/params.json", params(suite));
# BenchmarkTools.save("benchmark/results.json", minimum(run(suite)))

# Check: load params/results, then run and compare
prior_params = BenchmarkTools.load("benchmark/params.json")[1]
loadparams!(suite, prior_params, :evals, :samples)
prior = BenchmarkTools.load("benchmark/results.json")[1]
posterior = minimum(run(suite))
print(judge(posterior, prior))


text:
using MechanicalSketch
import MechanicalSketch: foil_spline_local
import MechanicalSketch: text, circle, Turtle, Pencolor, Penwidth, Forward, Turn
import MechanicalSketch: HueShift, O, sethue, finish, EM, WI, background, empty_figure
let
    empty_figure(filename = joinpath(@__DIR__, "test_1.png"));
    background("midnightblue")

    posx = -WI / 2
    sethue("green")
    # The unicode ∈ may not be supported by the default font and show as a square instead.
    stri = "Da jeg var på vei til kirken ∈ dag morges så kom jeg forbi en liten sjømann. En frisk og hyggelig liten sjømann som hilste meg."
    text(stri, posx, 0)
    sethue("yellow")

    text("1234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890",
        posx, EM)
    circle(O, 20, :stroke)
    t = Turtle()
    Pencolor(t, "cyan")
    Penwidth(t, 1.5)
    n = 5
    for i in 1:400
        Forward(t, n)
        Turn(t, 89.5)
        HueShift(t)
        n += 0.75
    end
    finish()
end


text:
using JuliaInterpreter: root, locals, moduleof, callee
import ..Atom: wsitem, handle

function contexts(state::DebuggerState = STATE)
  state.frame === nothing && return []
  ctx = []
  trace = ""
  frame = root(state.frame)
  active_callee = callee(active_frame(state))
  while frame ≠ nothing && frame ≠ active_callee
    trace = string(trace, "/", frame.framecode.scope isa Method ?
                                frame.framecode.scope.name : "top-level")
    c = Dict(:context => string("Debug: ", trace), :items => localvars(frame))
    pushfirst!(ctx, c)

    frame = frame.callee
  end
  ctx
end

function localvars(frame)
  vars = locals(frame)
  items = []
  mod = moduleof(frame)
  for v in vars
    # ref: https://github.com/JuliaDebug/JuliaInterpreter.jl/blob/master/src/utils.jl#L365-L370
    v.name == Symbol("#self#") && (isa(v.value, Type) || sizeof(v.value) == 0) && continue
    item = wsitem(mod, v.name, v.value)
    # Julia doesn't support "constantness" for local variables
    item[:type] == "constant" && (item[:type] = "variable")
    item[:icon] == "c" && (item[:icon] = "v")
    push!(items, item)
  end
  items
end


text:
using Turing
using Test
using Random

include("../utility.jl")

Random.seed!(123)

alg1 = HMCDA(3000, 1000, 0.65, 0.15)
alg2 = PG(20, 500)
alg3 = Gibbs(500, PG(30, 10, :s), HMCDA(1, 500, 0.65, 0.05, :m))

@model gdemo(x) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0, sqrt(s))
  x[1] ~ Normal(m, sqrt(s))
  x[2] ~ Normal(m, sqrt(s))
  return s, m
end

chn1 = sample(gdemo([1.5, 2.0]), alg1; save_state=true)

check_numerical(chn1, [:s, :m], [49/24, 7/6])

chn1_resumed = resume(chn1, 1000)

check_numerical(chn1_resumed, [:s, :m], [49/24, 7/6])

chn1_contd = sample(gdemo([1.5, 2.0]), alg1; resume_from=chn1)

check_numerical(chn1_contd, [:s, :m], [49/24, 7/6])

chn1_contd2 = sample(gdemo([1.5, 2.0]), alg1; resume_from=chn1, reuse_spl_n=1000)

check_numerical(chn1_contd2, [:s, :m], [49/24, 7/6])

chn2 = sample(gdemo([1.5, 2.0]), alg2; save_state=true)

check_numerical(chn2, [:s, :m], [49/24, 7/6])

chn2_contd = sample(gdemo([1.5, 2.0]), alg2; resume_from=chn2)

check_numerical(chn2_contd, [:s, :m], [49/24, 7/6])

chn3 = sample(gdemo([1.5, 2.0]), alg3; save_state=true)

check_numerical(chn3, [:s, :m], [49/24, 7/6])

chn3_contd = sample(gdemo([1.5, 2.0]), alg3; resume_from=chn3)

check_numerical(chn3_contd, [:s, :m], [49/24, 7/6])


text:
using PSCOPF

using .PSCOPFFixtures

using Test
using JuMP
using Dates
using DataStructures
using Printf

@testset verbose=true "test_bilevel_pilotables" begin

    TS = [DateTime("2015-01-01T11:00:00")]
    ech = DateTime("2015-01-01T07:00:00")
    next_ech = DateTime("2015-01-01T07:30:00")
    function create_instance(load_1, load_2,
                            limit::Float64=35.,
                            logs=nothing)
        network = PSCOPFFixtures.network_2buses(limit=limit)
        # Pilotables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                                0., 40.,
                                                0., 10.,
                                                Dates.Second(0), Dates.Second(0))
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_2", PSCOPF.Networks.PILOTABLE,
                                                0., 200.,
                                                0., 12.,
                                                Dates.Second(0), Dates.Second(0))
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_2", "prod_2_1", PSCOPF.Networks.PILOTABLE,
                                                0., 500.,
                                                0., 50.,
                                                Dates.Second(0), Dates.Second(0))

        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", TS[1], "S1", load_1)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_2", TS[1], "S1", load_2)

        gen_initial_state = SortedDict{String,PSCOPF.GeneratorState}(
                                "prod_1_1" => PSCOPF.ON,
                                "prod_1_2" => PSCOPF.OFF,
                                "prod_2_1" => PSCOPF.ON,
                            )
        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                    gen_initial_state, uncertainties, nothing, logs)

        return context
    end

    #=
    TS: [11h]
    S: [S1]
                        bus 1                   bus 2
                        |                      |
    (limitable) prod_1_1|       "1_2"          |prod_2_1
    Pmin=0, Pmax=40     |                      |Pmin=0, Pmax=500
    Csta=0, Cprop=10    |                      |Csta=0, Cprop=50
                        |----------------------|
                        |         35           |
                        |                      |
                        |                      |
    (limitable) prod_1_2|                      |
    Pmin=0, Pmax=200    |                      |
    Csta=0, Cprop=12    |                      |
                        |                      |
           load(bus_1)  |                      |load(bus_2)
      S1:  10           |                      | S1: 500

      ideally the market would use prod_1_1 and prod_1_2 first
      but the RSO constraint would not be satisfied
      => TSO needs to impose a maximum prod of 45 on bus 1 or a minimum prod of 465 on bus 2
      cost of unit prod_2_1 is high => impositions are high too (FIXME? this behaviour needs changing wrt to maximum imposition)
      TSO will :
        impose a maximum prod of 45 on bus 1:
            force prod_1_1 to 0 (it's the cheapest to impose)
            maximum of 45 on prod_1_2
            free prod_2_1

      The market is then obliged to use prod_1_2 which is more expensive than prod_1_1,
       when both satisfy the RSO constraint

    =#
    @testset "no_reference_market_no_delta" begin

        context = create_instance(10., 500.,
                                35.)

        tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(CONSIDER_DELTAS=false, USE_UNITS_PROP_COST_AS_TSO_BOUNDING_COST=true))
        firmness = PSCOPF.compute_firmness(tso,
                                            ech, next_ech,
                                            TS, context)
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

        # #TSO RSO constraints are OK
        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

        # #Market EOD constraints are OK
        @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

        #TSO sets the bounds for pilotable production respecting units' pmin and pmax
        #prod_1_1 is not bound
        @test value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"]) < 1e-09
        @test value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"]) < 1e-09
        #prod_1_2 is not bound
        @test value(result.upper.pilotable_model.p_imposition_min["prod_1_2",TS[1],"S1"]) < 1e-09
        @test 45. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_2",TS[1],"S1"])
        #prod_2_1 bounding + EOD, guarantees RSO constraints
        @test value(result.upper.pilotable_model.p_imposition_min["prod_2_1",TS[1],"S1"]) < 1e-09
        @test 500. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_2_1",TS[1],"S1"])

        #Market chooses the levels respecting the bounds for pilotable production
        @test value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"]) < 1e-09
        @test 45. ≈ value(result.lower.pilotable_model.p_injected["prod_1_2",TS[1],"S1"])
        @test 465. ≈ value(result.lower.pilotable_model.p_injected["prod_2_1",TS[1],"S1"])

        @test value(PSCOPF.get_upper_obj_expr(result)) ≈ ( 40. *10. + 155. *12. )
        @test value(PSCOPF.get_lower_obj_expr(result)) ≈ ( 45*12 + 465*50. )

    end

    #=
    TS: [11h]
    S: [S1]
                        bus 1                   bus 2
                        |                      |
    (limitable) prod_1_1|       "1_2"          |prod_2_1
    Pmin=0, Pmax=40     |                      |Pmin=0, Pmax=500
    Csta=0, Cprop=10    |                      |Csta=0, Cprop=50
                        |----------------------|
                        |         35           |
                        |                      |
                        |                      |
    (limitable) prod_1_2|                      |
    Pmin=0, Pmax=200    |                      |
    Csta=0, Cprop=12    |                      |
                        |                      |
           load(bus_1)  |                      |load(bus_2)
      S1:  10           |                      | S1: 500

      refenece market schedule is not considered here in terms of deltas
        (it still affects the cost function if there were previous UC decisions).
        => same behaviour as "no_reference_market_no_delta"

      ideally the market would use prod_1_1 and prod_1_2 first
      but the RSO constraint would not be satisfied
      => TSO needs to impose a maximum prod of 45 on bus 1 or a minimum prod of 465 on bus 2
      cost of unit prod_2_1 is high => impositions are high too (FIXME? this behaviour needs changing wrt to maximum imposition)
      TSO will :
        impose a maximum prod of 45 on bus 1:
            force prod_1_1 to 0 (it's the cheapest to impose)
            maximum of 45 on prod_1_2
            free prod_2_1

      The market is then obliged to use prod_1_2 which is more expensive than prod_1_1,
       when both satisfy the RSO constraint
    =#
    @testset "reference_market_no_delta" begin

        context = create_instance(10., 500.,
                                35.)

        context.market_schedule = PSCOPF.Schedule(PSCOPF.Market(), ech, SortedDict(
            "prod_1_1" => PSCOPF.GeneratorSchedule("prod_1_1",
                    SortedDict(),
                    SortedDict(TS[1] => PSCOPF.UncertainValue{PSCOPF.Float64}(missing,
                                                                                SortedDict("S1"=>40.,)))
                    ),
            "prod_1_2" => PSCOPF.GeneratorSchedule("prod_1_2",
                    SortedDict(),
                    SortedDict(TS[1] => PSCOPF.UncertainValue{PSCOPF.Float64}(missing,
                                                                                SortedDict("S1"=>20.,)))
                    ),
            "prod_2_1" => PSCOPF.GeneratorSchedule("prod_2_1",
                    SortedDict(),
                    SortedDict(TS[1] => PSCOPF.UncertainValue{PSCOPF.Float64}(missing,
                                                                                SortedDict("S1"=>450.,)))
                    )
            )
        )

        tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(CONSIDER_DELTAS=false, USE_UNITS_PROP_COST_AS_TSO_BOUNDING_COST=true))
        firmness = PSCOPF.compute_firmness(tso,
                                            ech, next_ech,
                                            TS, context)
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

        # #TSO RSO constraints are OK
        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

        # #Market EOD constraints are OK
        @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

        #TSO sets the bounds for pilotable production respecting units' pmin and pmax
        #prod_1_1 is not bound
        @test value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"]) < 1e-09
        @test value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"]) < 1e-09
        #prod_1_2 is not bound
        @test value(result.upper.pilotable_model.p_imposition_min["prod_1_2",TS[1],"S1"]) < 1e-09
        @test 45. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_2",TS[1],"S1"])
        #prod_2_1 bounding + EOD, guarantees RSO constraints
        @test value(result.upper.pilotable_model.p_imposition_min["prod_2_1",TS[1],"S1"]) < 1e-09
        @test 500. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_2_1",TS[1],"S1"])

        #Market chooses the levels respecting the bounds for pilotable production
        @test value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"]) < 1e-09
        @test 45. ≈ value(result.lower.pilotable_model.p_injected["prod_1_2",TS[1],"S1"])
        @test 465. ≈ value(result.lower.pilotable_model.p_injected["prod_2_1",TS[1],"S1"])

        @test value(PSCOPF.get_upper_obj_expr(result)) ≈ ( 40. *10. + 155. *12. )
        @test value(PSCOPF.get_lower_obj_expr(result)) ≈ ( 45*12 + 465*50. )

    end

    #=
    TS: [11h]
    S: [S1]
                        bus 1                   bus 2
                        |                      |
    (limitable) prod_1_1|       "1_2"          |prod_2_1
    Pmin=0, Pmax=40     |                      |Pmin=0, Pmax=500
    Csta=0, Cprop=10    |                      |Csta=0, Cprop=50
                        |----------------------|
                        |         35           |
                        |                      |
                        |                      |
    (limitable) prod_1_2|                      |
    Pmin=0, Pmax=200    |                      |
    Csta=0, Cprop=12    |                      |
                        |                      |
           load(bus_1)  |                      |load(bus_2)
      S1:  10           |                      | S1: 500

      ideally the market would use prod_1_1 and prod_1_2 first

      an artificial (not optimal) schedule is given here :
        prod_1_1:40, prod_1_2:20, prod_2_1:450
      => production on bus 1 is problematic (causes RSO constraint)

      => TSO needs to impose a maximum prod of 45 on bus 1 or a minimum prod of 465 on bus 2
      The TSO tries to deviate the least from the arket refence schedule
      cost of unit prod_2_1 is high => impositions are high too (FIXME? this behaviour needs changing wrt to maximum imposition)
      TSO will :
        impose a maximum prod of 45 on bus 1:
            constrain prod_1_1 (it's the cheapest to impose) => max at 20 (=> delta of 15)
            maxi at 20 on prod_1_2 (delta of 0)
            free prod_2_1 (production will compensate the reduced 15MW)
    =#
    @testset "reference_market_and_consider_delta" begin

        context = create_instance(10., 500.,
                                35.)

        context.market_schedule = PSCOPF.Schedule(PSCOPF.Market(), ech, SortedDict(
            "prod_1_1" => PSCOPF.GeneratorSchedule("prod_1_1",
                    SortedDict(),
                    SortedDict(TS[1] => PSCOPF.UncertainValue{PSCOPF.Float64}(missing,
                                                                                SortedDict("S1"=>40.,)))
                    ),
            "prod_1_2" => PSCOPF.GeneratorSchedule("prod_1_2",
                    SortedDict(),
                    SortedDict(TS[1] => PSCOPF.UncertainValue{PSCOPF.Float64}(missing,
                                                                                SortedDict("S1"=>20.,)))
                    ),
            "prod_2_1" => PSCOPF.GeneratorSchedule("prod_2_1",
                    SortedDict(),
                    SortedDict(TS[1] => PSCOPF.UncertainValue{PSCOPF.Float64}(missing,
                                                                                SortedDict("S1"=>450.,)))
                    )
            )
        )

        tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(CONSIDER_DELTAS=true, USE_UNITS_PROP_COST_AS_TSO_BOUNDING_COST=true))
        firmness = PSCOPF.compute_firmness(tso,
                                            ech, next_ech,
                                            TS, context)
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

        # #TSO RSO constraints are OK
        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

        # #Market EOD constraints are OK
        @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

        #TSO sets the bounds for pilotable production respecting units' pmin and pmax
        #prod_1_1 is not bound
        @test value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"]) < 1e-09
        @test 25. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"])
        #prod_1_2 is not bound
        @test value(result.upper.pilotable_model.p_imposition_min["prod_1_2",TS[1],"S1"]) < 1e-09
        @test 20. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_2",TS[1],"S1"])
        #prod_2_1 bounding + EOD, guarantees RSO constraints
        @test value(result.upper.pilotable_model.p_imposition_min["prod_2_1",TS[1],"S1"]) < 1e-09
        @test 500. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_2_1",TS[1],"S1"])

        #Market chooses the levels respecting the bounds for pilotable production
        @test 25. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"])
        @test 20. ≈ value(result.lower.pilotable_model.p_injected["prod_1_2",TS[1],"S1"])
        @test 465. ≈ value(result.lower.pilotable_model.p_injected["prod_2_1",TS[1],"S1"])

        @test value(PSCOPF.get_upper_obj_expr(result)) ≈ ( 15. *10. + 180. *12. )
        @test value(PSCOPF.get_lower_obj_expr(result)) ≈ ( 25. *10 + 20. *12. + 465*50. )

    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    filter(pred, sdata)

Retain all locations in spatial data `sdata` according to
a predicate function `pred`. A predicate function takes
table rows as input, e.g. `pred(r) = r.state == "CA"`.
"""
function filter(pred, sdata::AbstractData)
  𝒯 = values(sdata)
  𝒟 = domain(sdata)

  # row table view
  ctor = Tables.materializer(𝒯)
  rows = Tables.rows(𝒯)

  # indices to retain
  inds = findall(pred, rows)

  # return point set
  table = ctor(rows[inds])
  coord = coordinates(𝒟, inds)

  georef(table, coord)
end


text:
"""

	F = G ∩ H

"""

struct Intersection <: Constraint

	F::Variable
	G::Variable
	H::Variable

end

import Base.show
function Base.show(io::IO, inter::Intersection)
	print(io, inter.F.name * " = " * inter.G.name * " ∩ " * inter.H.name)
end

function Intersection(model::Model, F, G, H)
	constraint = Intersection(F, G, H)
	push!(F.linkedConstraint, constraint)
	push!(G.linkedConstraint, constraint)
	push!(H.linkedConstraint, constraint)

	push!(model.constraints, constraint)

	constraint
end

function filtrage!(inter::Intersection)

	stop = !(isempty(setdiff(intersect(inter.G.lowerBound, inter.H.lowerBound), inter.F.upperBound)))
	stop = stop || !(isempty(setdiff(inter.F.lowerBound, intersect(inter.G.upperBound, inter.H.upperBound))))

	changeVariable = Vector{Variable}(undef, 3)
	nbChange = 0

	changeF = Change(inter.F)
	changeG = Change(inter.G)
	changeH = Change(inter.H)

	if stop
		return [], (changeF, changeG, changeH), true

	#Set
	else
		if !inter.F.isFixed

			rem = setdiff(inter.F.upperBound, intersect(inter.G.upperBound, inter.H.upperBound))
			add = setdiff(intersect(inter.G.lowerBound, inter.H.lowerBound), inter.F.lowerBound)

			setdiff!(inter.F.upperBound, rem)
			union!(inter.F.lowerBound, add)

			cardInf = max(inter.F.cardinalInf, length(inter.F.lowerBound))
			cardSup = min(inter.F.cardinalSup, length(inter.F.upperBound), inter.G.cardinalSup, inter.H.cardinalSup)

			changeF.cardAdded = cardInf - inter.F.cardinalInf
			changeF.cardRemoved = cardSup - inter.F.cardinalSup

			inter.F.cardinalInf += changeF.cardAdded
			inter.F.cardinalSup += changeF.cardRemoved

			if !(isempty(rem) && isempty(add))
				nbChange += 1
				changeVariable[nbChange] = inter.F

				if isFixed(inter.F)
					fix!(inter.F, changeF)
				end
			end
			for elt in add
				push!(changeF.added, elt)
			end
			for elt in rem
				push!(changeF.removed, elt)
			end
		end

		if !inter.G.isFixed

			rem = intersect(inter.G.upperBound, setdiff(inter.H.lowerBound, inter.F.upperBound))
			add = setdiff(inter.F.lowerBound, inter.G.lowerBound)

			setdiff!(inter.G.upperBound, rem)
			union!(inter.G.lowerBound, add)

			cardInf = max(inter.G.cardinalInf, length(inter.G.lowerBound), inter.F.cardinalInf)
			cardSup = min(inter.G.cardinalSup, length(inter.G.upperBound))

			changeG.cardAdded = cardInf - inter.G.cardinalInf
			changeG.cardRemoved = cardSup - inter.G.cardinalSup

			inter.G.cardinalInf += changeG.cardAdded
			inter.G.cardinalSup += changeG.cardRemoved

			if !(isempty(rem) && isempty(add))
				nbChange += 1
				changeVariable[nbChange] = inter.G

				if isFixed(inter.G)
					fix!(inter.G, changeG)
				end
			end

			for elt in add
				push!(changeG.added, elt)
			end
			for elt in rem
				push!(changeG.removed, elt)
			end
		end

		if !inter.H.isFixed
			rem = intersect(inter.H.upperBound, setdiff(inter.G.lowerBound, inter.F.upperBound))
			add = setdiff(inter.F.lowerBound, inter.H.lowerBound)

			setdiff!(inter.H.upperBound, rem)
			union!(inter.H.lowerBound, add)

			cardInf = max(inter.H.cardinalInf, length(inter.H.lowerBound), inter.F.cardinalInf)
			cardSup = min(inter.H.cardinalSup, length(inter.H.upperBound))

			changeH.cardAdded = cardInf - inter.H.cardinalInf
			changeH.cardRemoved = cardSup - inter.H.cardinalSup

			inter.H.cardinalInf += changeH.cardAdded
			inter.H.cardinalSup += changeH.cardRemoved

			if !(isempty(rem) && isempty(add))
				nbChange += 1
				changeVariable[nbChange] = inter.H

				if isFixed(inter.H)
					fix!(inter.H, changeH)
				end
			end

			for elt in add
				push!(changeH.added, elt)
			end
			for elt in rem
				push!(changeH.removed, elt)
			end
		end

		#Error
		stop = stop || !(inter.F.cardinalInf <= inter.F.cardinalSup)
		stop = stop || !(inter.G.cardinalInf <= inter.G.cardinalSup)
		stop = stop || !(inter.H.cardinalInf <= inter.H.cardinalSup)

		changes = (changeF, changeG, changeH)
		for c in changes
			if !(isempty(c.added))
				#println("$inter ajoute $(c.added) à $(c.var)")
			end
			if !(isempty(c.removed))
				#println("$inter enlève $(c.removed) à $(c.var)")
			end
		end

		return changeVariable[1:nbChange], (changeF, changeG, changeH), stop
	end

end


text:
module ReduceNoise

using Base.Iterators
using ImageCore
using ImageCore.MappedArrays
using ImageCore.OffsetArrays
using ImageCore.PaddedViews
using ImageCore: NumberLike, GenericGrayImage, GenericImage
import ..NoiseAPI: AbstractImageDenoiseAlgorithm, reduce_noise, reduce_noise!

include("compat.jl")
include("BM3DDenoise.jl")
include("NonlocalMean.jl")

export
    reduce_noise, reduce_noise!,

    # BM3D
    BM3D,
    # Non-local mean filter for gaussian noise
    NonlocalMean, get_NonlocalMean_rp

end # module


text:
####################################################################
# Hyper rectangles
####################################################################
# Hyper rectangles are used to bound points in space.
# For an inner node all it's children are bounded by the
# inner nodes hyper rectangle.
immutable HyperRectangle{T <: AbstractFloat}
    mins::Vector{T}
    maxes::Vector{T}
end

# Splits a hyper rectangle into two rectangles by dividing the
# rectangle at a specific value in a given dimension.
function split_hyper_rec{T <: AbstractFloat}(hyper_rec::HyperRectangle{T},
                                             dim::Int,
                                             value::T)
    new_max = copy(hyper_rec.maxes)
    new_max[dim] = value

    new_min = copy(hyper_rec.mins)
    new_min[dim] = value

    return HyperRectangle(hyper_rec.mins, new_max),
           HyperRectangle(new_min, hyper_rec.maxes)
end

############################################
# Rectangle - Point functions
############################################
@inline function get_min_dim{T <: AbstractFloat}(rec::HyperRectangle{T}, point::Vector{T}, dim::Int)
    @inbounds d = abs2(max(0, max(rec.mins[dim] - point[dim], point[dim] - rec.maxes[dim])))
    d
end

@inline function get_max_dim{T <: AbstractFloat}(rec::HyperRectangle{T}, point::Vector{T}, dim::Int)
    @inbounds d = abs2(max(rec.maxes[dim] - point[dim], point[dim] - rec.mins[dim]))
    d
end

# Max distance between rectangle and point
@inline function get_max_distance{T <: AbstractFloat}(rec::HyperRectangle{T}, point::Vector{T})
    max_dist = zero(T)
    @inbounds @simd for dim in eachindex(point)
        max_dist += abs2(max(rec.maxes[dim] - point[dim], point[dim] - rec.mins[dim]))
    end
    return max_dist
end

# Min distance between rectangle and point
@inline function get_min_distance{T <: AbstractFloat}(rec::HyperRectangle{T}, point::Vector{T})
    min_dist = zero(T)
    @inbounds @simd for dim in eachindex(point)
        min_dist += abs2(max(0, max(rec.mins[dim] - point[dim], point[dim] - rec.maxes[dim])))
    end
    return min_dist
end

# (Min, Max) distance between rectangle and point
@inline function get_min_max_distance{T <: AbstractFloat}(rec::HyperRectangle{T}, point::Vector{T})
    min_dist = get_min_distance(rec, point)
    max_dist = get_max_distance(rec, point)
    return min_dist, max_dist
end

# (Min, Max) distance between rectangle and point for a certain dim
@inline function get_min_max_dim{T <: AbstractFloat}(rec::HyperRectangle{T}, point::Vector{T}, dim::Int)
    min_dist_dim = get_min_dim(rec, point, dim)
    max_dist_dim = get_max_dim(rec, point, dim)
    return min_dist_dim, max_dist_dim
end

############################################
# Rectangle-rectangle functions
############################################
@inline function get_min_dim{T <: AbstractFloat}(rec1::HyperRectangle{T},  rec2::HyperRectangle{T}, dim::Int)
    return abs2(max(0, max(rec1.mins[dim] - rec2.maxes[dim], rec2.mins[dim] - rec1.maxes[dim])))
end

@inline function get_max_dim{T <: AbstractFloat}(rec1::HyperRectangle{T},  rec2::HyperRectangle{T}, dim::Int)
    return abs2(max(rec1.maxes[dim] - rec2.mins[dim],  rec2.maxes[dim] - rec1.mins[dim]))
end

@inline function get_min_max_distance{T <: AbstractFloat}(rec1::HyperRectangle{T}, rec2::HyperRectangle{T})
    min_dist = get_min_distance(rec1, rec2)
    max_dist = get_max_distance(rec1, rec2)
    return min_dist, max_dist
end

@inline function get_min_distance{T <: AbstractFloat}(rec1::HyperRectangle{T}, rec2::HyperRectangle{T})
    min_dist = zero(T)
    @inbounds for dim in 1:length(rec1.maxes)
        min_dist += get_min_dim(rec1, rec2, dim)
    end
    return min_dist
end

@inline function get_max_distance{T <: AbstractFloat}(rec1::HyperRectangle{T}, rec2::HyperRectangle{T})
    max_dist = zero(T)
    @inbounds for dim in 1:length(rec1.maxes)
        max_dist += get_max_dim(rec1, rec2, dim)
    end
    return max_dist
end

@inline function get_min_max_dim{T <: AbstractFloat}(rec1::HyperRectangle{T},  rec2::HyperRectangle{T}, dim::Int)
    min_dist_dim = get_min_distance(rec1, rec2, dim)
    max_dist_dim = get_max_distance(rec1, rec2, dim)
    return min_dist_dim, max_dist_dim
end


text:
"""
# Download data from FRED

    get_symbols(symbol)

## Arguments

    - `symbol`: symbol, e.g. "GDPC1" or "FEDFUNDS"

## Examples

```jldoctest
julia> get_symbols("GDPC1")
julia> get_symbols("FEDFUNDS")
julia> get_symbols("T10Y2Y")
```
"""
function get_symbols(symbol::String)

    url = "https://api.stlouisfed.org/fred/series/observations"
    parameters = Dict("api_key" => ENV["API_FRED"],
                    "file_type" => "json",
                    "series_id" => symbol)

    response = HTTP.request("GET", url; query = parameters)
    observations = JSON.parse(String(response.body))["observations"]

    time = Dates.Date.([observations[i]["date"] for i ∈ eachindex(observations)])
    observations = [observations[i]["value"] for i ∈ eachindex(observations)]

    index_missing = findall(x -> x == ".", observations)
    deleteat!(observations, index_missing)
    deleteat!(time, index_missing)

    x = DataFrame(
        Value = [parse(Float64, observations[i]) for i ∈ eachindex(observations)],
        Time = time
    )

    x = TimeArray(x, timestamp = :Time)
    return(x)
end

"""
# Download data from FRED

    get_symbols(symbol, start_interval, end_interval)

## Arguments

    - `symbol`: symbol, e.g. "GDPC1" or "FEDFUNDS"
    - `f` : Date character, e.g. "2018-12-13" in the form of "YYYY-MM-DD"
    - `t` : Date character, e.g. "2019-12-13" in the form of "YYYY-MM-DD"

## Examples

```jldoctest
julia> get_symbols("GDPC1", "2010-01-01", "2019-01-01")
julia> get_symbols("FEDFUNDS", "2000-01-01", "2019-01-01")
julia> get_symbols("T10Y2Y", "2000-01-01", "2019-01-01")
```
"""
function get_symbols(symbol::String, start_interval::String, end_interval::String)

    x = get_symbols(symbol::String)

    x = from(x, Date(start_interval))
    x = to(x, Date(end_interval))

    return(x)
end


text:
using ADCME
using LinearAlgebra
using MAT
using AdFem
using PyPlot
using SparseArrays

function Q_exact(x, y, t)
    (2*x^2*y*(x - 1)*(2*y - 1)*(x + y + 1) + x*y*(y - 1)*(y^2 + 1)^2*(-x - y*(2*x - 1)*(y - 1) + 1)*(x + y + 1) + x*(x - 1)*(2*y - 1)*(y^2 + 1)^2 - (y^2 + 1)*(x + y + 1)*(2*x*(x - 1)*(x + (x^2 + 1)*(y^2 + 1)) + 2*y*(x + (x^2 + 1)*(y^2 + 1))*(y - 1) + y*(2*x - 1)*(y - 1)*(2*x*(y^2 + 1) + 1)))*exp(-t)/((y^2 + 1)^2*(x + y + 1))
end

function T_exact(x, y, t)
    x * (1-x) * y * (1-y) * exp(-t)
end

function u_exact(x, y)
    x * y * (1-y)
end

function v_exact(x, y)
    1 / (1 + x + y)
end

function k_exact(x, y)
    1 + x^2 + x / (1 + y^2)
end

function k_nn(x, y)
    out = fc([x y], [20,20,20,1])^2 + 0.1 # N x 1
    squeeze(out)
end

#---------------------------------------------
# grid setup
m = 20
n = 20
h = 1/n
NT = 10
t_final = 1.0
dt = t_final/NT
#---------------------------------------------
# discretized governing equation
#  J * u^{n+1} - J * u^n = - v ⋅ grad u + K * u^{n+1} + F^{n+1}
#---------------------------------------------
xy = fem_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
T0 = @. T_exact(x, y, 0.0)
u = @. u_exact(x,y)
v = @. v_exact(x,y)
k = @. k_nn(x, y); k=stack(k)

bd = bcnode("all", m, n, h)

ugauss = fem_to_gauss_points(u, m, n, h)
vgauss = fem_to_gauss_points(v, m, n, h)
kgauss = fem_to_gauss_points(k, m, n, h)

# precompute source term
Q = zeros(NT, (m+1)*(n+1))
for i = 1:NT
    t = i*dt
    f_value = eval_f_on_gauss_pts((x,y)->Q_exact(x, y, t), m, n, h)
    Q[i,:] = compute_fem_source_term1(f_value, m, n, h)
end
Q = constant(Q)

Mass = constant(compute_fem_mass_matrix1(m, n, h))
Advection = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h))
Laplace = constant(compute_fem_laplace_matrix1(kgauss, m, n, h))
A = 1/dt * Mass + Advection + Laplace
A, _ = fem_impose_Dirichlet_boundary_condition1(A, bd, m, n, h)

function solve_chtenergy_one_step(T, source)
    Tgauss = fem_to_gauss_points(T, m, n, h)
    b = constant(compute_fem_source_term1(constant(Tgauss / dt), m, n, h)) + source
    b = scatter_update(b, bd, zeros(length(bd)))
    sol = A\b
end

# ---------------------------------------------
function condition(i, T_arr)
    i<=NT+1
end

function body(i, T_arr)
    T = read(T_arr, i-1)
    T_new = solve_chtenergy_one_step(T, Q[i-1])
    T_arr = write(T_arr, i, T_new)
    i+1, T_arr
end

i = constant(2, dtype=Int32)

T_arr = TensorArray(NT+1)
T_arr = write(T_arr, 1, T0)
_, T_arr = while_loop(condition, body, [i, T_arr])
T_arr = set_shape(stack(T_arr), (NT+1, (m+1)*(n+1)))

T_data = matread("chtenergy_transient_data.mat")["T"]
# loss = mean((T_data .- T_arr)^2)

sample_size = 40
idx = rand(1:(m+1)*(n+1), sample_size)
observed_data = T_data[:, idx]
loss = mean((T_arr[:, idx].- observed_data)^2)

loss = loss * 1e10

# ---------------------------------------------------
# create a session and run
max_iter = 100
sess = Session(); init(sess)
loss_ = BFGS!(sess, loss, max_iter)
figure(); semilogy(loss_); savefig("chtenergy_transient_loss.png")

figure(figsize=(14,4));
subplot(131)
visualize_scalar_on_fem_points(k_exact.(x,y), m, n, h); title("conductivity exact")
subplot(132)
visualize_scalar_on_fem_points(run(sess, k), m, n, h); title("conductivity prediction")
subplot(133)
visualize_scalar_on_fem_points(k_exact.(x,y).-run(sess, k), m, n, h); title("conductivity difference")
savefig("chtenergy_transient_k.png")


text:
abstract type Runner; end

function target_nbits(target::AbstractString)
    if startswith(target, "i686-") || startswith(target, "arm-")
        return "32"
    else
        return "64"
    end
end

function target_proc_family(target::AbstractString)
    if startswith(target, "arm") || startswith(target, "aarch")
        return "arm"
    elseif startswith(target, "power")
        return "power"
    else
        return "intel"
    end
end

function target_dlext(target::AbstractString)
    if endswith(target, "-mingw32")
        return "dll"
    elseif occursin("-apple-", target)
        return "dylib"
    else
        return "so"
    end
end

"""
    platform_envs(platform::Platform)

Given a `platform`, generate a `Dict` mapping representing all the environment
variables to be set within the build environment to force compiles toward the
defined target architecture.  Examples of things set are `PATH`, `CC`,
`RANLIB`, as well as nonstandard things like `target`.
"""
function platform_envs(platform::Platform, host_target="x86_64-linux-gnu")
    global use_ccache

    # Convert platform to a triplet, but strip out the ABI parts
    target = triplet(abi_agnostic(platform))

    # Helper function to generate paths such as /opt/x86_64-apple-darwin14/bin/llvm-ar
    tool_path(n, t = target) = "/opt/$(t)/bin/$(n)"
    # Helper functions to generate paths such as /opt/x86_64-linux-gnu/bin/x86_64-linux-gnu-gcc
    target_tool_path(n, t = target) = tool_path("$(t)-$(n)", t)
    target_lib_dir(t) = "/opt/$(t)/lib64:/opt/$(t)/lib:/opt/$(t)/$(t)/lib64:/opt/$(t)/$(t)/lib"
    host_tool_path(n, t = host_target) = tool_path("$(t)-$(n)", t)
    llvm_tool_path(n) = "/opt/$(host_target)/tools/$(n)"

    # Start with the default musl ld path
    lib_path = "/usr/local/lib64:/usr/local/lib:/usr/local/lib:/usr/lib"

    # Add our glibc directory (this coupled with our Glibc patches to reject musl
    # binaries gives us a dual-libc environment)
    lib_path *= ":/lib64:/lib"

    # Then add on our target-specific library directories and our host-target library directories
    lib_path *= ":" * target_lib_dir(host_target)
    lib_path *= ":" * target_lib_dir(target)

    # Finally add on our destination location
    lib_path *= ":/workspace/destdir/lib64:/workspace/destdir/lib"

    # Start with the standard PATH:
    path = "/usr/local/bin:/usr/local/sbin:/usr/bin:/usr/sbin:/bin:/sbin"

    # Slip our tools into the front, followed by host tools
    path = "/opt/$(target)/bin:/opt/$(host_target)/bin:/opt/$(host_target)/tools:" * path

    # Then slip $prefix/bin onto the end, so that dependencies naturally show up
    path = path * ":/workspace/destdir/bin"

    mapping = Dict(
        # Activate the given target via `PATH` and `LD_LIBRARY_PATH`
        "PATH" => path,
        "LD_LIBRARY_PATH" => lib_path,

        # We conditionally add on some compiler flags; we'll cull empty ones at the end
        "CFLAGS" => "",
        "CPPFLAGS" => "",
        "LDFLAGS" => "",

        # binutils/toolchain envvars
        "DSYMUTIL" => llvm_tool_path("llvm-dsymutil"),
        "LIBTOOL" => target_tool_path("libtool"),
        "LIPO" => target_tool_path("lipo"),
        "OTOOL" => target_tool_path("otool"),
        "INSTALL_NAME_TOOL" => target_tool_path("install_name_tool"),
        "OBJCOPY" => target_tool_path("objcopy"),
        "READELF" => target_tool_path("readelf"),
        "RANLIB" => target_tool_path("ranlib"),
        "STRIP" => target_tool_path("strip"),
        "WINDRES" => target_tool_path("windres"),
        "WINMC" => target_tool_path("winmc"),
        "LLVM_TARGET" => target,
        "LLVM_HOST_TARGET" => host_target,

        # Useful tools for our buildscripts
        "target" => target,
        "nproc" => "$(Sys.CPU_THREADS)",
        "nbits" => target_nbits(target),
        "proc_family" => target_proc_family(target),
        "dlext" => target_dlext(target),
        "TERM" => "screen",

        # We should always be looking for packages already in the prefix
        "PKG_CONFIG_PATH" => "/workspace/destdir/lib/pkgconfig:/workspace/destdir/share/pkgconfig",
        "PKG_CONFIG_SYSROOT_DIR" => "/workspace/destdir",

        # We like to be able to get at our .bash_history afterwards. :)
        "HISTFILE"=>"/meta/.bash_history",
    )

    # If we're on MacOS or FreeBSD, we default to LLVM tools instead of GCC.
    # On all of our clangy platforms we actually have GCC tools available as well,
    # they're just not used by default.  Override these environment variables in
    # your scripts if you want to use them.
    if occursin("-apple-", target) || occursin("-freebsd", target)
        mapping["AR"] = llvm_tool_path("llvm-ar")
        mapping["AS"] = llvm_tool_path("llvm-as")
        # Because there's only a single `clang` binary, we store it in `x86_64-linux-gnu`,
        # but LLVMBuilder puts the `clang` binary into `tools`, not `bin`.
        mapping["CC"] = llvm_tool_path("clang -target $(target) --sysroot /opt/$(target)/$(target)/sys-root")
        mapping["CXX"] = llvm_tool_path("clang++ -target $(target) --sysroot /opt/$(target)/$(target)/sys-root")

        # FreeBSD is special-cased within the LLVM source tree to not allow for
        # things like the -gcc-toolchain option, which means that we have to manually add
        # the location of libgcc_s.  LE SIGH.
        # https://github.com/llvm-mirror/clang/blob/f3b7928366f63b51ffc97e74f8afcff497c57e8d/lib/Driver/ToolChains/FreeBSD.cpp
        if occursin("-freebsd", target)
            mapping["LDFLAGS"] *= " -L/opt/$(target)/$(target)/lib"
        end

        # OSX must be linked against libc++.
        if occursin("-apple-", target)
            mapping["CXX"] *= " -stdlib=libc++"
        end
        # flang isn't a realistic option yet, so we still use gfortran here
        mapping["FC"] = target_tool_path("gfortran")
        mapping["LD"] = llvm_tool_path("llvm-ld")
        mapping["NM"] = llvm_tool_path("llvm-nm")
        mapping["OBJDUMP"] = llvm_tool_path("llvm-objdump")
    else
        mapping["AR"] = target_tool_path("ar")
        mapping["AS"] = target_tool_path("as")
        mapping["CC"] = target_tool_path("gcc")
        mapping["CXX"] = target_tool_path("g++")
        mapping["FC"] = target_tool_path("gfortran")
        mapping["LD"] = target_tool_path("ld")
        mapping["NM"] = target_tool_path("nm")
        mapping["OBJDUMP"] = target_tool_path("objdump")
    end

    # On OSX, we need to do a little more work.
    if typeof(platform) <: MacOS
        # First, tell CMake what our deployment target is, so that it tries to
        # set -mmacosx-version-min and such appropriately
        mapping["MACOSX_DEPLOYMENT_TARGET"] = "10.8"

        # Also put this into LDFLAGS because some packages are hard of hearing
        mapping["LDFLAGS"] *= " -mmacosx-version-min=10.8"

        # If we're on an old GCC, the -syslibroot that gets passed in isn't
        # calculated correctly, so we have to manually add it in.
        if compiler_abi(platform).gcc_version in (:gcc4, :gcc_any)
            mapping["LDFLAGS"] *= " -Wl,-syslibroot,/opt/$(target)/$(target)/sys-root"
        end
    end

    # There is no broad agreement on what host compilers should be called,
    # so we set all the environment variables that we've seen them called
    # and hope for the best.
    for host_map in (tool -> "HOST$(tool)", tool -> "$(tool)_FOR_BUILD", tool -> "BUILD_$(tool)")
        mapping[host_map("AR")] = host_tool_path("ar")
        mapping[host_map("AS")] = host_tool_path("as")
        mapping[host_map("CC")] = host_tool_path("gcc")
        mapping[host_map("CXX")] = host_tool_path("g++")
        mapping[host_map("DSYMUTIL")] = llvm_tool_path("llvm-dsymutil")
        mapping[host_map("FC")] = host_tool_path("gfortran")
        mapping[host_map("LIPO")] = host_tool_path("lipo")
        mapping[host_map("LD")] = host_tool_path("ld")
        mapping[host_map("NM")] = host_tool_path("nm")
        mapping[host_map("RANLIB")] = host_tool_path("ranlib")
        mapping[host_map("READELF")] = host_tool_path("readelf")
        mapping[host_map("OBJCOPY")] = host_tool_path("objcopy")
        mapping[host_map("OBJDUMP")] = host_tool_path("objdump")
        mapping[host_map("STRIP")] = host_tool_path("strip")
    end

    # If we're using `ccache`, prepend it to `CC`, `CXX`, `FC`, `HOSTCC`, etc....
    if use_ccache
        for tool in ("CC", "CXX", "FC")
            for m in (tool, "BUILD_$(tool)", "HOST$(tool)", "$(tool)_FOR_BUILD")
                mapping[m] = string("ccache ", mapping[m])
            end
        end
    end

    return mapping
end

runner_override = ""
function preferred_runner()
    global runner_override
    if runner_override != ""
        if runner_override in ["userns", "privileged"]
            return UserNSRunner
        elseif runner_override in ["qemu"]
            return QemuRunner
        elseif runner_override in ["docker"]
            return DockerRunner
        end
    end

    @static if Sys.islinux()
        return UserNSRunner
    else
        return QemuRunner
    end
end

"""
    runshell(platform::Platform = platform_key_abi())

Launch an interactive shell session within the user namespace, with environment
setup to target the given `platform`.
"""
function runshell(platform::Platform = platform_key_abi(); verbose::Bool = false)
    runshell(preferred_runner(), platform; verbose=verbose)
end

function runshell(r::Runner, args...; kwargs...)
    run_interactive(r, `/bin/bash`, args...; kwargs...)
end

function runshell(::Type{R}, platform::Platform = platform_key_abi(); verbose::Bool = false) where {R <: Runner}
    return runshell(R(pwd(); cwd="/workspace/", platform=platform, verbose=verbose))
end


text:
# ------------ Identify blocks which are viable as kernels ------------ #

struct KernelHint <: ProgramStructureHint
    loops::Set{NaturalLoop}
    mapped_ir::Dict{NaturalLoop, Vector{Block}}
end
function Base.display(kh::KernelHint)
    if isempty(kh.loops)
        println("\u001b[32m ✓ (KernelHint): Compiler detected no kernels in your model code.\u001b[0m")
    else
        println("\u001b[33m (KernelHint): Compiler detected the following kernels in your model code.\u001b[0m")
        println("________________________\n")
        for l in kh.loops
            for b in kh.mapped_ir[l]
                display(b)
            end
            println("________________________\n")
        end
        println("\u001b[32m (Recommendation): extract kernels into a Fold combinator for easier analysis and optimization.")
    end
end

function detect_kernels(ir::IR)
    loops = detectloops(ir)
    kh = KernelHint(Set{NaturalLoop}([]), Dict{NaturalLoop, Set{Block}}())
    for l in loops
        length(l.backedges) > 1 && continue
        backedge = collect(l.backedges)[1]
        back_blk = block(ir, backedge)
        header = l.header
        head_blk = block(ir, header)
        any(map(branches(back_blk)) do br
                br.block == header && length(br.args) == length(arguments(head_blk))
            end) || continue
        push!(kh.loops, l)
        kh.mapped_ir[l] = map(sort(collect(l.body))) do ind
            block(ir, ind)
        end
    end
    kh
end

function detect_kernels(fn, arg_types...)
    ir = lower_to_ir(fn, arg_types...)
    display(detect_kernels(ir))
end


text:
using OrdinaryDiffEq,DiffEqProblemLibrary, DiffEqDevTools, Test, ODEInterfaceDiffEq

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_linear

prob = prob_ode_linear
sol = solve(prob,DP5())

sol2 = solve(prob,dopri5())
@test sol.t[2] ≈ sol2.t[2]

prob = prob_ode_2Dlinear
sol = solve(prob,DP5(),internalnorm=(u)->sqrt(sum(abs2,u)))

# Change the norm due to error in dopri5.f
sol2 = solve(prob,dopri5())
@test sol.t[2] ≈ sol2.t[2]

prob = deepcopy(prob_ode_linear)
prob2 = ODEProblem(prob.f, prob.u0, (1.0,0.0),1.01)
sol = solve(prob2,DP5())

sol2 = solve(prob2,dopri5())
@test sol.t[2] ≈ sol2.t[2]

prob = deepcopy(prob_ode_2Dlinear)
prob2 = ODEProblem(prob.f, prob.u0, (1.0,0.0), 1.01)
sol = solve(prob2,DP5(),internalnorm=(u)->sqrt(sum(abs2,u)))

# Change the norm due to error in dopri5.f
sol2 = solve(prob2,dopri5())
@test sol.t[2] ≈ sol2.t[2]


text:
type Node{T}
    data::T
    next::Nullable{Node{T}}

    function Node(data::T)
        n = new()
        n.data = data
        # To mark the end of the list we use the Nullable{T}() function .
        n.next = Nullable{Node{T}}()
        n
    end
end

# convenience. Let use write Node(10) or Node(10.0) instead of Node{Int64}(10), Node{Float64}(10.0)
function Node(data)
    return Node{typeof(data)}(data)
end

islast(n::Node) = (isnull(n.next))

function append{T}(n::Node{T}, data::T)
    tmp = Node(data)
    if !islast(n)
        tmp.next = n.next
    end
    n.next = tmp
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function __init__()
    global HANDLE_GIVEN_ID = Dict{UInt64,Symbol}()
    global DEFAULT_RESOURCE = Ref{AbstractResource}(CPU1())
    global DEFAULT_SCITYPE_CHECK_LEVEL = Ref{Int}(1)
    global SHOW_COLOR = Ref{Bool}(true)

    # for testing asynchronous training of learning networks:
    global TESTING = parse(Bool, get(ENV, "TEST_MLJBASE", "false"))
    if TESTING
        global MACHINE_CHANNEL =
            RemoteChannel(() -> Channel(100), myid())
    end

    MLJModelInterface.set_interface_mode(MLJModelInterface.FullInterface())
end


text:
export amoeba

"""
    amoeba(f; alg=Polygonal(), options...)

Compute the amoeba of `f` with the given algorithm `alg` which can be
* `Polygonal()`
* `Greedy()`
* `ArchTrop()`
* `Simple()`

but you probably only want to use `Polygonal()` or `Greedy()`.

Depending on the given algorithm the function takes different keyword arguments.

## Examples

```julia
@polyvar x y

# This uses the `Polygonal()` algorithm
amoeba(x^2 + y^2 + 1)

# We only want a crude approximation
amoeba(x^2 + y^2 + 1, accuracy=0.1)

# Use the `Greedy()` algorithm.
amoeba(x^2 + y^2 + 1, alg=Greedy())

# Use the `Greedy()` algorithm with a custom grid.
grid = Grid2D(xlims=(-5, 5), ylims=(-4, 4), res=(500, 400))
amoeba(x^2 + y^2 + 1, alg=Greedy(), grid=grid)

# Use the `Greedy()` algorithm with the default domain but a higher resolution
amoeba(x^2 + y^2 + 1, alg=Greedy(), resolution=800)
```

## `Polygonal()`
This algorithm computes an approximation of the amoeba ``\\mathcal{A}(f)`` in the provided `domain` Ω
by computing a set of polygons 𝒫 with |𝒫| ⊂ Ω such that the union of these polygons approximates ``\\mathcal{A}(f)`` ∩ Ω
from the outside.

The possible (optional) arguments are
* `domain`: A tuple in the form `(xmin, xmax, ymin, ymax)` which defines a section Ω
for which the amoeba ``\\mathcal{A}(f)`` is computed. This domain has to be such that the intersection
Ω ∩ ``\\mathcal{A}(f)`` still captures the correct topology of ``\\mathcal{A}(f)``.
* `accuracy=0.01`: The maximal allowed error |𝒫 - ``\\mathcal{A}(f)`` ∩ Ω|. Note that we only compute an upper limit of the error.
The algorithms stops if the given accuracy is reached.
* `spine`: This algorithm needs the spine of ``\\mathcal{A}(f)``.
* `minimal_component_size=0.01`: The minimal size of the components of the complement. This is only used if no spine
is passed explicitly.
* `iterations=2000`: The maximal number of iterations.
* `vertices_accuracy=accuracy*1e-3`: During the algorithm we approximate points on the boundary of ``\\mathcal{A}(f)``. This
is the accuracy with which we compute them. Note hat this influences the minimal error of the approximation
and it should always be some magnitudes smaller than `accuracy`.
* `membership_options=[MembershipTestOptions()](@ref)`: As a subroutine a membership test is used.

## `Greedy()`, `Simple()`, `ArchTrop()`

These algorithms are all approximations of the amoeba ``\\mathcal{A}(f)`` based on a grid. This basically
applies the membership test for different grid points. `Greedy()` is the fastest and
`Simple()` the slowest.
The grid can be passed explicitly, otherwise it will be computed based on a heuristic.

* `resolution=600`: The resolution of the grid if not passed explicitly.
* `grid`: If passed explicitly this grid is used, otherwise it will be computed based on a heuristic.
* `membership_options=[MembershipTestOptions()](@ref)`: The options for the membership test.
"""
function amoeba end

function amoeba(f::MP.AbstractPolynomial; kwargs...)
    if MP.nvariables(f) == 2
        F = AmoebaFiber2D(f)
    elseif MP.nvariables(f) == 3
        F = AmoebaFiber3D(f)
    else
        throw(error("Currently only 2 and 3-dimensional PolynomialAmoebas are supported."))
    end
    amoeba(F, f; kwargs...)
end

function amoeba(
        F::AbstractAmoebaFiber,
        f::MP.AbstractPolynomial;
        alg::AbstractAlgorithm=default_algorithm(F),
        kwargs...)
    amoeba(F, f, alg; kwargs...)
end

default_algorithm(F::AbstractAmoebaFiber{T, 2}) where T = Polygonal()
default_algorithm(F::AbstractAmoebaFiber{T, 3}) where T = Greedy()

function amoeba(F::AmoebaFiber3D, f, alg::Polygonal; kwargs...)
    throw(error("The `Polygonal` algorithm is only supported for bivariate polynomials."))
end
function amoeba(
        F::AmoebaFiber2D{T},
        f::MP.AbstractPolynomial,
        alg::Polygonal;
        domain=default_amoeba_domain(f),
        accuracy=1e-2,
        vertices_accuracy = accuracy * 1e-3,
        iterations = 2000,
        membership_options = MembershipTestOptions(),
        minimal_component_size = 0.01,
        spine = Spine2D(f, minimal_component_size=minimal_component_size, domain=domain),
        report_progress=false
        ) where {T}

    C = Covering(spine, domain...)
    M = CoveringFitter(C, F;
        desired_maximal_error=accuracy, vertices_accuracy=vertices_accuracy,
        maximal_iterations=iterations, options=membership_options)
    fit!(M, report_progress=report_progress)

    polygons = cartesian_polygons(M.covering)
    error_bound = M.approximated_global_error

    PolygonalAmoeba(spine, polygons, float.(domain), error_bound, M.iteration)
end

function amoeba(F::AbstractAmoebaFiber{T, N}, f::MP.AbstractPolynomial,
    alg::AbstractGridAlgorithm;
    resolution=600,
    grid=default_grid(F, f, resolution),
    callback=_do_nothing,
    membership_options=MembershipTestOptions()
    ) where {T, N}
    if !(typeof(grid) <: AbstractGrid{N})
        return throw(error("The grid is not of dimension $N"))
    end

    _amoeba(alg, F, f, grid, TorusStartValueGenerator{N}(), membership_options, callback)
end

@inline default_grid(F::AbstractAmoebaFiber{T, 2}, f::MP.AbstractPolynomial, r) where T = Grid2D(f, r)
@inline default_grid(F::AbstractAmoebaFiber{T, 3}, f::MP.AbstractPolynomial, r) where T = throw(error("You have to provide a `Grid3D` via grid=..."))

function default_amoeba_domain(f::MP.AbstractPolynomial; factor=1.5, aspect_ratio=:default)
    try
        return amoeba_carcase_domain_heuristic(f, factor=factor, aspect_ratio=aspect_ratio)
    catch err
        throw(error("You have to pass an explicit domain."))

    end
end

_amoeba(::Simple, F, f, grid, gen, options, callback) = simple_grid(F, grid, gen, options=options)
function _amoeba(::ArchTrop, F::AmoebaFiber2D, f, grid, gen, options, callback)
    nb = archimedean_neighbourhood(f, grid).data
    avoid_check = (B, k) -> (!nb[k], false)
    simple_grid(F, grid, gen, avoid_check, options=options)
end

function _amoeba(::Greedy, F::AmoebaFiber2D, f, grid, gen, options, callback)
    initial_queue = findall(grid_contour(f, grid).data)

    greedy_grid(F, grid, gen, initial_queue, options=options, callback=callback)
end

function _amoeba(::Greedy, F::AmoebaFiber3D, f, grid, gen, options, callback)
    n1, n2, n3 = size(grid)

    z = MP.variables(f)[end]
    k0 = div(n3, 3)
    f_z0 = MP.subs(f, z => exp(grid.zrange[k0]))
    grid_k0 = Grid2D(grid.xrange, grid.yrange)
    contour_f_z0 = grid_contour(f_z0, grid_k0).data
    initial_queue = map(findall(contour_f_z0)) do ij
        CartesianIndex(Tuple(ij)..., k0)
    end

    greedy_grid(F, grid, gen, initial_queue, options=options, callback=callback)
end


text:
@testset "simulate to maxtime" begin

    rng = MersenneTwister(12)
    input = MultilevelInput(
        modulesize=4,
        fixedmu=true,
        b=0.1,
        d=0.1,
        bdrate=0.1,
        clonalmutations=0,
        maxtime=400,
        maxmodules=1000,
        branchrate=3/365,
        branchfraction=0.2,
        μ=1
    )
    population = multilevel_simulation(input, rng, :fixedtime)
    @test age(population) < 400

    b, d = 0.1, 0.0
    modulesize = 6
    branchinitsize = 2
    branchrate = 1/5
    maxtime = 100

    moduletracker = SomaticEvolution.initialize_population(
        modulesize,
        clonalmutations=0
    )[1]

    rng = MersenneTwister(12)

    moduletracker, newmoduletracker =
        SomaticEvolution.module_simulate_to_branching!(
            moduletracker,
            maxtime,
            b,
            d,
            b,
            branchrate,
            modulesize,
            branchinitsize,
            2,
            rng
    )
    @test moduletracker.tvec[end] <= maxtime
    @test moduletracker.Nvec[end] == modulesize - branchinitsize
    @test newmoduletracker.Nvec[end] == branchinitsize
    @test newmoduletracker.tvec[end] == moduletracker.tvec[end]

    newmaxtime = moduletracker.tvec[end] + 0.1
    moduletracker, newmoduletracker =
        SomaticEvolution.module_simulate_to_branching!(
            moduletracker,
            newmaxtime,
            b,
            d,
            b,
            branchrate,
            modulesize,
            branchinitsize,
            2,
            rng
    )
    @test newmoduletracker === nothing
end

input = MultilevelInput(
    modulesize=4,
    fixedmu=true,
    b=0.1,
    d=0,
    bdrate=0.01,
    clonalmutations=0,
    maxtime=1*365,
    branchrate=3/365,
    branchfraction=0.2,
    μ=1
)

mt1 = SomaticEvolution.ModuleTracker(
    [1, 2, 3, 4, 4, 4, 4, 4, 4, 4, 3, 4, 4, 4, 4],
    [0.0, 14.437821887595856, 38.636412371018324, 38.80067279227201, 39.07905423063021, 56.24183190197407, 61.57555456465887, 93.99343227770015, 174.7383395589789, 177.87796865140686, 187.1205185574646, 190.3156218593258, 222.64188808574013, 235.8640994797343, 242.08245554076214],
    Cell[Cell([3, 4, 20], 1), Cell([3, 5, 6, 8, 10, 17, 18, 19], 1), Cell([3, 4, 21, 22, 23], 1), Cell([1, 24], 1)],
    SomaticEvolution.CloneTracker[], 1, 0)
mt2 = SomaticEvolution.ModuleTracker(
    [1, 2, 3, 4, 4, 4, 3],
    [187.1205185574646, 207.5419622150601, 209.92705453342273, 216.2145518957684, 243.24240821372297, 252.31115142542512, 257.22794191422554],
    Cell[Cell([1, 13, 25, 26, 27], 1), Cell([1, 13], 1), Cell([1, 13, 25, 28, 29], 1)],
    SomaticEvolution.CloneTracker[], 2, 1)
mt3 = SomaticEvolution.ModuleTracker(
    [1],
    [257.22794191422554],
    Cell[Cell([1], 1)], SomaticEvolution.CloneTracker[], 3, 2)

@testset "module sampling" begin
    rng = MersenneTwister(12)
    moduletracker = deepcopy(mt1)
    moduletracker, newmoduletracker =
        SomaticEvolution.sample_new_module!(moduletracker, 2, 1, 1.0, rng)
    @test length(moduletracker) == 3
    @test length(newmoduletracker) == 1
    @test newmoduletracker.id == 2
    @test newmoduletracker.Nvec == [1]
    @test newmoduletracker.tvec == [1.0]
end

population = SomaticEvolution.Population(input, [mt1, mt2, mt3])

#check that statistics are calculated correctly
@testset "mutation statistics" begin
    @test mutations_per_cell(mt1) == [3, 8, 5, 2]
    @test average_mutations_per_module(population) ≈ [4.5, 4, 1]
    @test all(average_mutations(population, true) .≈ (3.875, 5.267857142857143))
    @test clonal_mutations(population) == [0, 2, 1]
    @test clonal_mutation_ids(population) == [[], [1, 13], [1]]
    @test pairwise_fixed_differences_matrix(population, diagonals=true) == [0 0 0; 2 2 0; 1 1 1]
    @test pairwise_fixed_differences_matrix(population, [2,3]) == [0 0; 1 0]
    @test pairwise_fixed_differences(population) == Dict{Int64, Int64}(1=>2, 2=>1)
    @test pairwise_fixed_differences(population, [2,3]) == Dict{Int64, Int64}(1=>1)
    @test shared_fixed_mutations(population) == Dict{Int64, Int64}(1=>1, 2=>1)
    @test shared_fixed_mutations(population, [2,3]) == Dict{Int64, Int64}(1=>1, 2=>1)
    @test all(pairwise_fixed_differences_statistics(population, clonal=true) .≈ (1.3333333333333333,0.33333333333333333,1,1))
    @test newmoduletimes(population) ≈ [0.0, 187.1205185574646, 257.22794191422554]
    @test all(cellpopulationsize(population, 50) .≈ ([0,50,100,150,200,250],[1, 4, 4, 4, 5, 8]))
end

@testset "updates" begin
    rng = MersenneTwister(12)
    populationtracker = [deepcopy(mt1), deepcopy(mt1), deepcopy(mt3)]

    #kills only cell in module 3 => population size goes to 2
    SomaticEvolution.deathupdate!(populationtracker, 4, 3, rng)
    @test length(populationtracker) == 2

    #one module splits into two modules of length two
    SomaticEvolution.modulebranchingupdate!(populationtracker, 4, 2, 1.0, rng)
    @test length(populationtracker) == 3
    @test sum(map(x -> length(x.cells), populationtracker)) == 8
    @test sum(map(x -> x.Nvec[end], populationtracker)) == 8

    #birth update one module gets an extra cell
    SomaticEvolution.birthupdate!(populationtracker, 4, 2, 1.0, 1, rng)
    @test sum(map(x -> length(x.cells), populationtracker)) == 9
    @test sum(map(x -> x.Nvec[end], populationtracker)) == 9

    #moran update population size stays the same
    SomaticEvolution.moranupdate!(populationtracker, 4, 3, 1.0, 1, rng)
    @test sum(map(x -> length(x.cells), populationtracker)) == 9
    @test sum(map(x -> x.Nvec[end], populationtracker)) == 9

    N0 = populationtracker[1].Nvec[end]
    @test N0 == length(populationtracker[1].cells)
    L = length(populationtracker[1].Nvec)
    @test L == length(populationtracker[1].tvec)
    SomaticEvolution.updatemodulehistory!(populationtracker[1], 1, 20.0)
    @test populationtracker[1].tvec[end] == 20.0
    @test populationtracker[1].Nvec[end] == 1 + N0
    @test length(populationtracker[1].Nvec) == length(populationtracker[1].tvec) == L + 1
    @test_throws ErrorException SomaticEvolution.check_module_number(populationtracker, 1)

end

@testset "simulate to fixed size" begin
    b = 0.1
    d = 0.0
    bdrate = 0.01
    modulesize = 4
    branchinitsize = 1
    branchrate = 3/365
    maxtime = 1000
    maxmodules = 5
    populationtracker = SomaticEvolution.initialize_population(modulesize, clonalmutations=0)
    mutID = SomaticEvolution.getmutID(populationtracker[1].cells) #get id of next mutation
    t = 0
    rng = MersenneTwister(1)

    #check first update is correct
    transitionrates = SomaticEvolution.get_transitionrates(
        populationtracker,
        b,
        d,
        bdrate,
        branchrate,
        modulesize
    )
    #only tranisition with non-zero rate is birth which has rate b=0.1
    @test transitionrates == [0.0, 0.1, 0.0, 0.0]
    populationtracker, mutID, t =
        SomaticEvolution.update_population!(
            populationtracker,
            b,
            d,
            bdrate,
            branchrate,
            modulesize,
            branchinitsize,
            mutID,
            0,
            maxtime,
            rng
    )
    #popsize is now 2
    @test populationtracker[1].Nvec[end] == length(populationtracker[1].cells)
    @test length(populationtracker[1].Nvec) === length(populationtracker[1].tvec) == 2
    #each cell has one mutation
    @test length(populationtracker[1].cells[1].mutations) == length(populationtracker[1].cells[2].mutations) == 1

    #check transition rates for population of modules m1, m2, m3
    transitionrates = SomaticEvolution.get_transitionrates(
        [mt1, mt2, mt3],
        b,
        d,
        bdrate,
        branchrate,
        modulesize
    )
    @test all(transitionrates .≈ [0.04, 0.4, 0.0, 3/365])

    #check transitionrates for non-zero death rate
    transitionrates = SomaticEvolution.get_transitionrates(
        [mt1, mt2, mt3],
        b,
        0.01,
        bdrate,
        branchrate,
        modulesize
    )
    @test all(transitionrates .≈ [0.04, 0.4, 0.04, 3/365])

    #check only homeostatic modules are chosen for Moran and module branching and that prob of choosing cells is equal
    @test all(SomaticEvolution.choose_homeostaticmodule([mt1,mt2,mt3], 4, rng)==mt1 for i in 1:10)
    @test isapprox(mean(reduce(hcat, [[v for v in SomaticEvolution.choose_homeostaticmodule_cells([mt1,mt2,mt3], 4, rng)[2:3]] for i in 1:10000]), dims=2), [2.5, 2.5], atol=0.1)
    #check prob of choosing non-homeostatic modules is proportional to number of cells
    @test isapprox(mean(length(SomaticEvolution.choose_growingmodule_cell([mt1,mt2,mt3], 4, rng)[1].cells) for i in 1:100000), 2.5, atol=0.1)

    #check multilevel_simulation function
    rng = MersenneTwister(12)
    input = MultilevelInput(
            modulesize=4,
            fixedmu=true,
            b=0.1,
            d=0.01,
            bdrate=0.01,
            clonalmutations=0,
            maxtime=20*365,
            maxmodules=10,
            branchrate=3/365,
            branchfraction=0.2,
            μ=1
    )
    population = multilevel_simulation(input, rng, :normal)
    @test length(population) == 10
    input = MultilevelInput(
            modulesize=4,
            fixedmu=true,
            b=0.1,
            d=0.1,
            bdrate=0.1,
            clonalmutations=0,
            maxtime=365,
            maxmodules=1000,
            branchrate=3/365,
            branchfraction=0.2,
            μ=1
    )
    population = multilevel_simulation(input, rng, :normal)
    @test age(population) <= 365

end


text:
module NonlinearSamplingTest

using Test
using Random
using ForneyLab
using ForneyLab: outboundType, isApplicable, prod!, logPdf, unsafeMean, unsafeVar, Sampling, requiresBreaker, breakerParameters
using ForneyLab: SPNonlinearSInGX, SPNonlinearSOutNGX, SPNonlinearSOutNM, SPNonlinearSIn1MN, SPNonlinearSOutNMX, SPNonlinearSInMX, MNonlinearSInMGX, gradientOptimization
using ForwardDiff

Random.seed!(1234)

f(x) = x
g(x) = x^2
h(x, y) = x + y

@testset "gradientOptimization" begin
    f(x) = -x[1]^2 - x[2]^2
    grad(x) = ForwardDiff.gradient(f, x)
    res = gradientOptimization(f, grad, [4.0, 1.0], 0.001)

    # Optimum is [0.0, 0.0]
    @test isapprox(res, [0.0, 0.0], atol=1e-16)
end

@testset "requiresBreaker and breakerParameters" begin
    fg = FactorGraph()
    x = Variable()
    y = Variable()
    nd = GaussianMeanVariance(x, 0.0, 1.0)
    Nonlinear{Sampling}(y, x, g=g)

    @test !requiresBreaker(nd.i[:out]) # Single-input Nonlinear{Sampling} does not require breaker

    fg = FactorGraph()
    x = Variable()
    z = Variable()
    y = Variable()
    GaussianMeanVariance(z, 0.0, 1.0)
    nd = GaussianMeanVariance(x, 0.0, 1.0)
    Nonlinear{Sampling}(y, x, z, g=h)

    @test requiresBreaker(nd.i[:out]) # Multi-input Nonlinear{Sampling} requires breaker
    @test_throws Exception breakerParameters(nd.i[:out].partner)
    @test breakerParameters(nd.i[:out]) == (Message{GaussianMeanVariance, Univariate}, ())
end

#-------------
# Update rules
#-------------

@testset "SPNonlinearSOutNM" begin
    @test SPNonlinearSOutNM <: SumProductRule{Nonlinear{Sampling}}
    @test outboundType(SPNonlinearSOutNM) == Message{SampleList}
    @test isApplicable(SPNonlinearSOutNM, [Nothing, Message{Gaussian}])
    @test isApplicable(SPNonlinearSOutNM, [Nothing, Message{Bernoulli}])

    msg = ruleSPNonlinearSOutNM(f, nothing, Message(Univariate, GaussianMeanVariance, m=2.0, v=tiny), n_samples=1)
    @test isapprox(msg.dist.params[:s][1], 2.0, atol=1e-4)
    @test msg.dist.params[:w] == [1.0]
    msg = ruleSPNonlinearSOutNM(f, nothing, Message(Multivariate, GaussianMeanVariance, m=[2.0], v=mat(tiny)), n_samples=1)
    @test isapprox(msg.dist.params[:s][1][1], 2.0, atol=1e-4)
    @test msg.dist.params[:w] == [1.0]
end

@testset "SPNonlinearSIn1MN" begin
    @test SPNonlinearSIn1MN <: SumProductRule{Nonlinear{Sampling}}
    @test outboundType(SPNonlinearSIn1MN) == Message{Function}
    @test isApplicable(SPNonlinearSIn1MN, [Message{Union{Bernoulli, Beta, Categorical, Dirichlet, Gaussian, Gamma, LogNormal, Poisson, Wishart}}, Nothing])

    log_pdf(x) = ruleSPNonlinearSIn1MN(f, Message(Univariate, GaussianMeanVariance, m=2.0, v=1.0), nothing, n_samples=0).dist.params[:log_pdf](x)
    @test log_pdf(1.5) == logPdf(ProbabilityDistribution(Univariate, GaussianMeanVariance, m=2.0, v=1.0), 1.5)
end

@testset "SPNonlinearSOutNGX" begin
    @test SPNonlinearSOutNGX <: SumProductRule{Nonlinear{Sampling}}
    @test outboundType(SPNonlinearSOutNGX) == Message{SampleList}
    @test !isApplicable(SPNonlinearSOutNGX, [Nothing, Message{Poisson}])
    @test isApplicable(SPNonlinearSOutNGX, [Nothing, Message{Gaussian}, Message{Gaussian}])

    msg = ruleSPNonlinearSOutNGX(h, nothing, Message(GaussianMeanVariance, m=1.0, v=tiny), Message(GaussianMeanVariance, m=2.0, v=tiny), n_samples=1)
    @test isapprox(msg.dist.params[:s][1], 3.0, atol=1e-4)
    @test msg.dist.params[:w] == [1.0]
end

@testset "SPNonlinearSInGX" begin
    @test SPNonlinearSInGX <: SumProductRule{Nonlinear{Sampling}}
    @test outboundType(SPNonlinearSInGX) == Message{GaussianWeightedMeanPrecision}
    @test !isApplicable(SPNonlinearSInGX, [Nothing, Message{Gamma}])
    @test isApplicable(SPNonlinearSInGX, [Message{Gaussian}, Message{Gaussian}, Nothing])

    msg_out = Message(Univariate, GaussianMeanVariance, m=2.0, v=1.0)
    msg_in1 = Message(Univariate, GaussianMeanVariance, m=0.0, v=1.0)
    msg_in2 = Message(Univariate, GaussianMeanVariance, m=1.0, v=1.0)

    res = ruleSPNonlinearSInGX(h, 1, msg_out, msg_in1, msg_in2, n_samples=1000)

    @test isapprox(mean(res.dist), 0.9998084187686186, atol=0.1)
    @test isapprox(var(res.dist), 1.9999999999999982, atol=0.1)
end

@testset "prod!" begin
    d = prod!(ProbabilityDistribution(Multivariate, Function, log_pdf=(s)->s), ProbabilityDistribution(Multivariate, Function, log_pdf=(s)->s))
    @test isa(d, ProbabilityDistribution{Multivariate,Function})
    @test d.params[:log_pdf](1) == 2
end

@testset "SPNonlinearSOutNMX" begin
    @test SPNonlinearSOutNMX <: SumProductRule{Nonlinear{Sampling}}
    @test outboundType(SPNonlinearSOutNMX) == Message{SampleList}
    @test !isApplicable(SPNonlinearSOutNMX, [Nothing, Message{Gamma}])
    @test !isApplicable(SPNonlinearSOutNMX, [Message{Gaussian}, Message{Gaussian}, Nothing])
    @test !isApplicable(SPNonlinearSOutNMX, [Message{Gaussian}, Message{Gamma}, Nothing])
    @test isApplicable(SPNonlinearSOutNMX, [Nothing, Message{Gaussian}, Message{Gamma}])
end

@testset "SPNonlinearSInMX" begin
    @test SPNonlinearSInMX <: SumProductRule{Nonlinear{Sampling}}
    @test outboundType(SPNonlinearSInMX) == Message{Function}
    @test !isApplicable(SPNonlinearSInMX, [Message{FactorFunction}, Message{Gamma}])
    @test !isApplicable(SPNonlinearSInMX, [Message{Gaussian}, Message{Gaussian}, Nothing])
    @test !isApplicable(SPNonlinearSInMX, [Nothing, Message{Gaussian}, Message{Gamma}])
    @test isApplicable(SPNonlinearSInMX, [Message{FactorFunction}, Nothing, Message{Gaussian}, Message{Gamma}])
end

@testset "MNonlinearSInMGX" begin
    @test MNonlinearSInMGX <: MarginalRule{Nonlinear{Sampling}}
    @test isApplicable(MNonlinearSInMGX, [Nothing, Message{Gaussian}, Message{Gaussian}])
    @test !isApplicable(MNonlinearSInMGX, [Nothing, Message{Gaussian}])
    @test !isApplicable(MNonlinearSInMGX, [Message{Gaussian}, Message{Gaussian}, ProbabilityDistribution])
end

#------------
# Integration
#------------

@testset "Nonlinear integration via sampling" begin
    fg = FactorGraph()

    @RV x ~ GaussianMeanVariance(2.0, 1.0)
    @RV y ~ GaussianMeanVariance(2.0, 3.0)
    n = Nonlinear{Sampling}(y, x, g=g, n_samples=2000)

    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(y)
    code = algorithmSourceCode(algo)
    @test occursin("ruleSPNonlinearSOutNM(g, nothing, messages[2], n_samples=2000)", code)
end

@testset "Nonlinear integration via sampling" begin
    fg = FactorGraph()

    @RV x ~ GaussianMeanVariance(2.0, 1.0)
    @RV y ~ GaussianMeanVariance(2.0, 3.0)
    @RV z ~ GaussianMeanVariance(2.0, 3.0)
    n = Nonlinear{Sampling}(z, x, y, g=g)

    # Forward; g_inv should not be present in call
    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(y)
    code = algorithmSourceCode(algo)

    @test occursin("ruleSPNonlinearSInGX(g, 2, messages[3], messages[2], messages[1])", code)
end

@testset "Nonlinear integration via sampling with specified variate types" begin
    fg = FactorGraph()

    @RV x ~ GaussianMeanVariance([2.0], mat(1.0))
    @RV y ~ GaussianMeanVariance([2.0], mat(3.0))
    @RV z ~ GaussianMeanVariance([2.0], mat(3.0))
    n = Nonlinear{Sampling}(z, x, y, g=g, dims=[(1,), (1,), (1,)])

    # Forward; g_inv should not be present in call
    pfz = PosteriorFactorization(fg)
    algo = messagePassingAlgorithm(y)
    code = algorithmSourceCode(algo)

    @test occursin("ruleSPNonlinearSInGX(g, 2, messages[3], messages[2], messages[1], dims=(1,))", code)
end

end # module