text:
export smooth, forget, propagate_constants, deepcopy, conjoin, replace_node,
    split, clone, merge, split_candidates, random_split, split_step, struct_learn,
    clone_candidates, standardize_circuit, make_vars_contiguous

"""
    smooth(root::StructLogicCircuit)::StructLogicCircuit

Create an equivalent smooth circuit from the given circuit.
"""
function smooth(root::StructLogicCircuit)::StructLogicCircuit
    (false_node, true_node) = canonical_constants(root)
    @assert (false_node === nothing && true_node === nothing) "You should propagate constants before smoothing a structured circuit!"

    lit_nodes = canonical_literals(root)
    f_con(n) = (n, BitSet())
    f_lit(n) = (n, BitSet(variable(n)))
    f_a(n, call) = begin
        (prime, pscope) = call(n.prime)
        prime = fill_missing_vtree(prime, vtree(n).left, vtree(prime), lit_nodes)
        pscope = variables(vtree(n).left)
        (sub, sscope) = call(n.sub)
        sub = fill_missing_vtree(sub, vtree(n).right, vtree(sub), lit_nodes)
        sscope = variables(vtree(n).right)

        parent_scope = pscope ∪ sscope
        smoothed = conjoin([prime, sub])
        @assert variables(vtree(smoothed)) == parent_scope
        (smoothed, parent_scope)
    end
    f_o(n, call) = begin
        parent_scope = mapreduce(c -> call(c)[2], union, children(n))
        smooth_children = Vector{LogicCircuit}(undef, num_children(n))
        map!(smooth_children, children(n)) do child
            (smooth_child, scope) = call(child)
            smooth_node(smooth_child, parent_scope, scope, lit_nodes)
        end
        return (disjoin([smooth_children...]; reuse=n), parent_scope)
    end
    foldup(root, f_con, f_lit, f_a, f_o, Tuple{LogicCircuit, BitSet})[1]
end
"""
    smooth(root::LogicCircuit)

Create an equivalent smooth circuit from the given circuit.
"""
function smooth(root::LogicCircuit)
    lit_nodes = canonical_literals(root)
    f_con(n) = (n, BitSet())
    f_lit(n) = (n, BitSet(variable(n)))
    f_a(n, call) = begin
        parent_scope = BitSet()
        smooth_children = Vector{LogicCircuit}(undef, num_children(n))
        map!(smooth_children, children(n)) do child
            (smooth_child, scope) = call(child)
            union!(parent_scope, scope)
            smooth_child
        end
        (conjoin([smooth_children...]; reuse=n), parent_scope)
    end
    f_o(n, call) = begin
        parent_scope = mapreduce(c -> call(c)[2], union, children(n))
        smooth_children = Vector{LogicCircuit}(undef, num_children(n))
        map!(smooth_children, children(n)) do child
            (smooth_child, scope) = call(child)
            smooth_node(smooth_child, setdiff(parent_scope, scope), lit_nodes)
        end
        return (disjoin([smooth_children...]; reuse=n), parent_scope)
    end
    foldup(root, f_con, f_lit, f_a, f_o, Tuple{LogicCircuit, BitSet})[1]
end

"""
    smooth_node(node::StructLogicCircuit, parent_scope, scope, lit_nodes)

Return a smooth version of the node where
the are added to the scope by filling the gap in vtrees, using literals from `lit_nodes`
"""
function smooth_node(node::StructLogicCircuit, parent_scope, scope, lit_nodes)
    # Compute the vtrees based on variables used
    target_vtr = lca(map(x -> vtree(lit_nodes[x]), collect(parent_scope))...)
    curr_vtr = lca(map(x -> vtree(lit_nodes[x]), collect(scope))...)
    if curr_vtr == target_vtr
        @assert isempty(setdiff(parent_scope, scope))
        node # If the node is where it should be on the vtree, we're done
    else
        fill_missing_vtree(node, target_vtr, curr_vtr, lit_nodes)
    end
end

"Construct a smoothed node from start_vtr, when you get to end_vtr insert the original node"
function fill_missing_vtree(node::StructLogicCircuit, start_vtr, end_vtr, lit_nodes)
    # If we're at the end just return
    if start_vtr == end_vtr
        node
    elseif isleaf(start_vtr)
        # If we're at a leaf node, just return the disjunction of literals
        get_lit(l,vtr) = get!(() -> compile(StructLogicCircuit, vtr, l), lit_nodes, l)
        lit = var2lit(start_vtr.var)
        lit_node = get_lit(lit, start_vtr)
        not_lit_node = get_lit(-lit, start_vtr)
        disjoin([lit_node, not_lit_node]; use_vtree=start_vtr)
    else
        # We're at an inner node, so recurse and conjoin the results?
        left = fill_missing_vtree(node, start_vtr.left, end_vtr, lit_nodes)
        right = fill_missing_vtree(node, start_vtr.right, end_vtr, lit_nodes)
        conjoin(left, right, use_vtree=start_vtr)
    end
end

"""
    smooth_node(node::LogicCircuit, missing_scope, lit_nodes)

Return a smooth version of the node where
the `missing_scope` variables are added to the scope, using literals from `lit_nodes`
"""
function smooth_node(node::LogicCircuit, missing_scope, lit_nodes)
    if isempty(missing_scope)
        return node # avoid adding unnecessary nodes
    else
        get_lit(l) = get!(() -> compile(node, l), lit_nodes, l)
        ors = map(collect(missing_scope)) do v
            lit = var2lit(Var(v))
            lit_node = get_lit(lit)
            not_lit_node = get_lit(-lit)
            disjoin([lit_node, not_lit_node]; reuse=node)
        end
        conjoin([node, ors...]; reuse=node)
    end
end

"""
    forget(root::LogicCircuit, is_forgotten::Function)

Forget variables from the circuit.
Warning: this may or may not destroy the determinism property.
"""
function forget(root::LogicCircuit, is_forgotten::Function, )
    (_, true_node) = canonical_constants(root)
    if isnothing(true_node)
        true_node = compile(typeof(root), true)
    end
    f_con(n) = n
    f_lit(n) = is_forgotten(variable(n)) ? true_node : n
    f_a(n, cn) = conjoin([cn...]; reuse=n) # convert type of cn
    f_o(n, cn) = disjoin([cn...]; reuse=n)
    foldup_aggregate(root, f_con, f_lit, f_a, f_o, LogicCircuit)
end

"""
    propagate_constants(root::LogicCircuit)

Remove all constant leafs from the circuit
"""
function propagate_constants(root::LogicCircuit; remove_unary=false)
    (false_node, true_node) = canonical_constants(root)
    if isnothing(true_node)
        true_node = compile(typeof(root), true)
    end
    if isnothing(false_node)
        false_node = compile(typeof(root), false)
    end
    f_con(n) = n
    f_lit(n) = n
    f_a(n, cn) = begin
        if any(c -> isfalse(c), cn)
            false_node
        else
            T = promote_type(typeof(n), typeof.(cn)...)
            proped_children = convert(Vector{T}, filter(c -> !istrue(c), cn))
            if isempty(proped_children)
                return true_node
            elseif remove_unary && (length(proped_children) == 1)
                return proped_children[1]
            else
                return conjoin(proped_children; reuse=n)
            end
        end
    end
    f_o(n, cn) = begin
        if any(c -> istrue(c), cn)
            true_node
        else
            T = promote_type(typeof(n), typeof.(cn)...)
            proped_children = convert(Vector{T}, filter(c -> !isfalse(c), cn))
            if isempty(proped_children)
                return false_node
            elseif remove_unary && (length(proped_children) == 1)
                return proped_children[1]
            else
                return disjoin(proped_children; reuse=n)
            end
        end
    end
    foldup_aggregate(root, f_con, f_lit, f_a, f_o, LogicCircuit)
end

import Base.deepcopy

"""
    deepcopy(n::LogicCircuit, depth::Int64)

Recursively create a copy circuit rooted at `n` to a certain depth `depth`
"""
function deepcopy(n::LogicCircuit, depth::Int64, old2new::Dict{LogicCircuit, LogicCircuit} = Dict{LogicCircuit, LogicCircuit}(); cache=true)
    # TODO: can we remove the Dict and use the data field instead?
    if depth == 0 || isliteralgate(n) || isconstantgate(n)
        return n
    elseif cache && haskey(old2new, n)
        return old2new[n]
    else
        cns = map(children(n)) do c
            deepcopy(c, depth - 1, old2new; cache=cache)
        end
        if is⋀gate(n)
            new_n = conjoin(cns)
        else
            @assert is⋁gate(n)
            new_n = disjoin(cns)
        end
        if cache
            old2new[n] = new_n
        end
        return new_n
    end
end

"""
    conjoin(root::LogicCircuit, lit::Lit; callback::Function)

Return the circuit conjoined with th given literal constrains
`callback` is called after modifying conjunction node
"""
function conjoin(root::LogicCircuit, lit::Lit; callback=noop)
    literals = canonical_literals(root)
    (false_node, ) = canonical_constants(root) # reuse constants when possible
    if isnothing(false_node)
        try
            false_node = compile(typeof(root), false)
        catch e
            false_node = nothing
        end
    end

    if !haskey(literals, -lit)
        root
    elseif haskey(literals, -lit) && !haskey(literals, lit)
        false_node
    else
        f_con(n) = (n, true)
        f_lit(n) = begin
            if variable(n) == lit2var(lit) && literal(n) != lit
                (false_node, false)
            else
                (n, true)
            end
        end
        f_a(n, cv) = begin
            kept = last.(cv)
            new_children = first.(cv)[kept]
            if all(kept)
                (conjoin([new_children...]; reuse=n), true)
            else
                (false_node, false)
            end
        end
        f_o(n, cv) = begin
            kept = last.(cv)
            new_children = first.(cv)[kept]
            if any(kept) && (length(new_children) == 1) && isliteralgate(new_children[1])
                (new_children[1], true)
            elseif any(kept)
                new_n = disjoin([new_children...]; reuse=n)
                callback(new_n, n, kept)
                (new_n, true)
            else
                (false_node, false)
            end
        end
        foldup_aggregate(root, f_con, f_lit, f_a, f_o, Tuple{Union{Nothing,LogicCircuit}, Bool})[1]
    end
end

import Base.split

"""
    split(root::LogicCircuit, (or, and)::Tuple{LogicCircuit, LogicCircuit}, var::Var; depth=0, sanity_check=true)

Return the circuit after spliting on edge `edge` and variable `var`
"""
function split(root::LogicCircuit, (or, and)::Tuple{LogicCircuit, LogicCircuit}, var::Var; depth=0, sanity_check=true)
    # sanity check
    if sanity_check
        @assert depth >= 0
        @assert is⋁gate(or) && is⋀gate(and) && and in children(or)
        @assert or in or_nodes(root) && and in and_nodes(root)
        literals = canonical_literals(and)
        @assert haskey(literals, var2lit(var)) && haskey(literals, - var2lit(var))
    end

    # split
    new_children1 = map(children(and)) do c
        conjoin(c, var2lit(var))
    end

    new_children2 = map(children(and)) do c
        conjoin(c, - var2lit(var))
    end

    new_and1 = deepcopy(conjoin(new_children1), depth; cache=false)
    new_and2 = deepcopy(conjoin(new_children2), depth; cache=false)
    new_or = disjoin([[new_and1, new_and2]; filter(c -> c != and, children(or))])

    replace_node(root, or, new_or), new_or
end

"""
    split_candidates(circuit::LogicCircuit)::Tuple{Vector{Tuple{LogicCircuit, LogicCircuit}}, Dict{LogicCircuit, BitSet}}

Return the edges and variables which can be splited on
"""
function split_candidates(circuit::LogicCircuit)::Tuple{Vector{Tuple{LogicCircuit, LogicCircuit}}, Dict{LogicCircuit, BitSet}}
    candidates::Vector{Tuple{LogicCircuit, LogicCircuit}} = Vector{Tuple{LogicCircuit, LogicCircuit}}(undef, 0)
    variable_scope::Dict{LogicCircuit, BitSet} = Dict{LogicCircuit, BitSet}()

    f_con(n) = begin
        variable_scope[n] = BitSet()
        (false, BitSet(), BitSet())
    end
    f_lit(n) = begin
        lit = literal(n)
        var = variable(n)
        variable_scope[n] = BitSet(var)

        lit > 0 ? (false, variable_scope[n], BitSet()) : (false, BitSet(), variable_scope[n])
    end
    f_a(n, cv) = begin
        pos_scope = reduce(union, [item[2] for item in cv])
        neg_scope = reduce(union, [item[3] for item in cv])
        variable_scope[n] = intersect(pos_scope, neg_scope)

        (!isdisjoint(pos_scope, neg_scope), pos_scope, neg_scope)
    end
    f_o(n, cv) = begin
        pos_scope = reduce(union, [item[2] for item in cv])
        neg_scope = reduce(union, [item[3] for item in cv])
        variable_scope[n] = intersect(pos_scope, neg_scope)

        map(zip(first.(cv), children(n))) do (splitable, c)
            if splitable
                push!(candidates, (n, c))
            end
        end

        (any(first.(cv)), pos_scope, neg_scope)
    end

    foldup_aggregate(circuit, f_con, f_lit, f_a, f_o, Tuple{Bool, BitSet, BitSet})

    candidates, variable_scope
end

function clone_candidates(circuit::LogicCircuit)::Dict{LogicCircuit, Vector{LogicCircuit}}
    candidates = Dict{LogicCircuit, Vector{LogicCircuit}}()
    parents = Dict{LogicCircuit, Vector{LogicCircuit}}()

    f_con(n) = begin
        false
    end
    f_lit(n) = begin
        false
    end
    f_a(n, cv) = begin
        for c in children(n)
            if !is⋁gate(c)
                continue
            end

            if c in keys(parents)
                push!(parents[c], n)
            else
                parents[c] = [n]
            end
        end
        false
    end
    f_o(n, cv) = begin
        if !(n in keys(parents))
            parents[n] = []
        end
        true
    end

    foldup_aggregate(circuit, f_con, f_lit, f_a, f_o, Bool)

    # Find candidates
    candidates = filter(p->(length(last(p)) == 2), parents) # Set of AND gates shared by exactly 2 OR gates
    candidates
end

"""
Randomly picking egde and variable from candidates
"""
function random_split(circuit::LogicCircuit)
    candidates, variable_scope = split_candidates(circuit)
    if isempty(candidates) return nothing end
    or, and = rand(candidates)
    vars = Var.(collect(variable_scope[and]))
    var = rand(vars)
    (or, and), var
end

"""
Split step
"""
function split_step(circuit::LogicCircuit; loss=random_split, depth=0, sanity_check=true)
    res = loss(circuit)
    if isnothing(res) return nothing end
    edge, var = loss(circuit)
    split(circuit, edge, var; depth=depth, sanity_check=sanity_check)
end

"""
Structure learning manager
"""
function struct_learn(circuit::LogicCircuit;
    primitives=[split_step],
    kwargs=Dict(split_step=>(loss=random_split, depth=0)),
    maxiter=typemax(Int), stop::Function=x->false, verbose = true)

    valid_prims = copy(primitives)
    m = length(valid_prims)

    for iter in 1 : maxiter
        which_primitive = rand(1:m)
        primitive_step = valid_prims[which_primitive]
        kwarg = kwargs[primitive_step]
        r = primitive_step(circuit; kwarg...)
        if isnothing(r)
            verbose && println("No more candidates to ", primitive_step, ". Skipping...")
            deleteat!(valid_prims, which_primitive)
            m -= 1
            if m == 0
                verbose && println("No more transformations available.")
                return circuit
            end
        end
        c2, _ = r
        if stop(c2)
            return c2
        end
        circuit = c2
    end
    circuit
end

"""
Clone the `or` node and redirect one of its parents to the new copy
"""
function clone(root::LogicCircuit, and1::LogicCircuit, and2::LogicCircuit, or::LogicCircuit; depth=1)
    # sanity check
    @assert depth >= 1
    @assert is⋀gate(and1) && is⋀gate(and2) && is⋁gate(or)
    @assert or in children(and1) && or in children(and2)
    @assert and1 in and_nodes(root) && and2 in and_nodes(root) && or in or_nodes(root)

    # clone
    new_or = deepcopy(or, depth; cache=false)
    new_and2 = conjoin([[new_or]; filter(c -> c != or, children(and2))])

    replace_node(root, and2, new_and2)
end

import Base.merge
"""
    merge(root::LogicCircuit, or1::LogicCircuit, or2::LogicCircuit)

Merge two circuits.
"""
function merge(root::LogicCircuit, or1::LogicCircuit, or2::LogicCircuit)
    # sanity check
    @assert is⋁gate(or1) && is⋁gate(or2)
    @assert or1 in or_nodes(root) && or2 in or_nodes(root)

    # Simply replace the occurrencce of the larger OR with the smaller OR
    n_or1 = num_nodes(or1)
    n_or2 = num_nodes(or2)

    if n_or1 >= n_or2
        replace_node(root, or1, or2)
    else
        replace_node(root, or2, or1)
    end
end

"""
Replace node `old` with node `new` in circuit `root`
"""
function replace_node(root::LogicCircuit, old::LogicCircuit, new::LogicCircuit; callback::Function=(x, y, z) -> nothing)
    # @assert GateType(old) == GateType(new) # this assertion was removed because it seems unnecessary?
    f_con(n) = old == n ? new : n
    f_lit = f_con
    f_a(n, cns) = old == n ? new : conjoin([cns...]; reuse=n)
    f_o(n, cns) = old == n ? new : begin
            new_n = disjoin([cns...]; reuse=n)
            callback(new_n, n, trues(length(cns)))
            new_n
        end
    foldup_aggregate(root, f_con, f_lit, f_a, f_o, LogicCircuit)
end

"""
Standardize the circuit:

1. Children of or nodes are and nodes.
2. Children of and nodes are or nodes.
3. Each and node has exactly two children.

Note: for circuits with parameters this function will not keep the
      parameters equavalent.
"""
function standardize_circuit(circuit::LogicCircuit)
    compile_and_node(left_child::LogicCircuit, right_child::LogicCircuit) = begin
        if is⋀gate(left_child)
            left_child = disjoin(left_child)
        end
        if is⋀gate(right_child)
            right_child = disjoin(right_child)
        end
        disjoin(conjoin(left_child, right_child))
    end

    f_con(node) = compile(node, constant(node))
    f_lit(node) = compile(node, literal(node))
    f_a(node, cs) = begin
        node_children::Array{LogicCircuit} = Array{LogicCircuit}(undef, 0)

        @assert length(cs) > 1 "⋀gate $(node) should contain at least two children, but only has one."

        # Flatten the children by ensuring all children of `node` are ⋁gate or leaf gate.
        for child_node in cs
            if is⋀gate(child_node)
                for n in children(child_node)
                    push!(node_children, n)
                end
            else
                push!(node_children, child_node)
            end
        end

        if length(node_children) > 2
            n = compile_and_node(node_children[end - 1], node_children[end])
            for idx = length(node_children) - 2 : -1 : 2
                n = compile_and_node(node_children[idx], n)
            end
            node_children = [first(node_children), n]
        end

        conjoin(node_children...)
    end
    f_o(node, cs) = begin
        node_children::Array{LogicCircuit} = Array{LogicCircuit}(undef, 0)

        for child_node in cs
            if is⋁gate(child_node)
                for n in children(child_node)
                    push!(node_children, n)
                end
            else # is⋀gate(child_node) || isleafgate(child_node)
                push!(node_children, child_node)
            end
        end

        disjoin(node_children...)
    end

    foldup_aggregate(circuit, f_con, f_lit, f_a, f_o, LogicCircuit)
end

"""
Make all variables in this circuit contiguously numbered. Return new circuit and the variable mapping.
"""
function make_vars_contiguous(root::LogicCircuit)
    var_bijection = [(v, i) for (i,v) in enumerate(variables(root))]
    reIndex_vars(root, var_bijection), var_bijection
end

function reIndex_vars(root::LogicCircuit, var_bijection::AbstractVector)
    var_dict = Dict(var_bijection)
    var2lits = Dict(map(var_bijection) do (v,i)
        pos_lit = compile(typeof(root),  Lit(i))
        neg_lit = compile(typeof(root), -Lit(i))
        (v => (pos_lit, neg_lit))
    end)
    f_con(n) = n
    f_lit(n) = begin
        if var_dict[variable(n)] == variable(n)
            n
        else
            lits = var2lits[variable(n)]
            ispositive(n) ? lits[1] : lits[2]
        end
    end
    f_a(n, cn) = conjoin([cn...]; reuse=n)
    f_o(n, cn) = disjoin([cn...]; reuse=n)
    root2 = foldup_aggregate(root, f_con, f_lit, f_a, f_o, LogicCircuit)
    root2
end


text:
@testset "model.jl" begin
    @testset "lastVar" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 6, "x", trailer)
        y = SeaPearl.IntVar(2, 6, "y", trailer)

        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)

        valueselection = SeaPearl.BasicHeuristic()
        @test isnothing(model.statistics.lastVar)
        v = valueselection(SeaPearl.DecisionPhase, model, x)
        @test model.statistics.lastVar.id=="x"
        SeaPearl.assign!(x, v)
        v = valueselection(SeaPearl.DecisionPhase, model, y)
        @test model.statistics.lastVar.id=="y"
        SeaPearl.assign!(y, v)

    end

    @testset "addVariable!()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(2, 6, "x", trailer)
        y = SeaPearl.IntVar(2, 6, "y", trailer)

        model = SeaPearl.CPModel(trailer)

        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)

        @test length(model.variables) == 2
        @test SeaPearl.branchable_variables(model) == Dict{String, SeaPearl.AbstractVar}(["x" => x, "y" => y])

        z = SeaPearl.IntVar(2, 6, "y", trailer)

        @test_throws AssertionError SeaPearl.addVariable!(model, z)

        # Not branching
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(2, 6, "x", trailer)
        y = SeaPearl.IntVar(2, 6, "y", trailer)

        model = SeaPearl.CPModel(trailer)

        SeaPearl.addVariable!(model, x; branchable=false)
        SeaPearl.addVariable!(model, y)

        @test length(model.variables) == 2
        @test SeaPearl.branchable_variables(model) == Dict{String, SeaPearl.AbstractVar}(["y" => y])

        # Trying to branch on Set variable
        trailer = SeaPearl.Trailer()
        y = SeaPearl.IntSetVar(2, 6, "y", trailer)

        model = SeaPearl.CPModel(trailer)

        @test_throws AssertionError SeaPearl.addVariable!(model, y)
    end

    @testset "merge!()" begin
        test1 = SeaPearl.CPModification("x" => [2, 3, 4],"z" => [11, 12, 13, 14, 15],"y" => [7, 8])
        test2 = SeaPearl.CPModification("x" => [5],"y" => [7, 8])

        SeaPearl.merge!(test1, test2)

        @test test1 == SeaPearl.CPModification("x" => [2, 3, 4, 5],"z" => [11, 12, 13, 14, 15],"y" => [7, 8, 7, 8])
    end

    @testset "addToPrunedDomains!()" begin
        test1 = SeaPearl.CPModification("x" => [2, 3, 4],"z" => [11, 12, 13, 14, 15],"y" => [7, 8])
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(2, 6, "x", trailer)
        t = SeaPearl.IntVar(2, 6, "t", trailer)
        b = SeaPearl.BoolVar("b", trailer)

        SeaPearl.addToPrunedDomains!(test1, x, [5, 6])

        @test test1 == SeaPearl.CPModification("x" => [2, 3, 4, 5, 6],"z" => [11, 12, 13, 14, 15],"y" => [7, 8])

        SeaPearl.addToPrunedDomains!(test1, t, [5, 6])

        @test test1 == SeaPearl.CPModification("x" => [2, 3, 4, 5, 6],"z" => [11, 12, 13, 14, 15],"y" => [7, 8], "t" => [5, 6])

        SeaPearl.addToPrunedDomains!(test1, b, [true])

        @test test1 == SeaPearl.CPModification("x" => [2, 3, 4, 5, 6],"z" => [11, 12, 13, 14, 15],"y" => [7, 8], "t" => [5, 6], "b" => [true])

    end

    @testset "solutionFound()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(2, 6, "x", trailer)
        y = SeaPearl.IntVar(2, 6, "y", trailer)

        model = SeaPearl.CPModel(trailer)

        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)

        @test !SeaPearl.solutionFound(model)

        constraint = SeaPearl.EqualConstant(x, 3, trailer)
        constraint2 = SeaPearl.Equal(x, y, trailer)
        SeaPearl.addConstraint!(model, constraint)
        SeaPearl.addConstraint!(model, constraint2)

        SeaPearl.fixPoint!(model)

        @test SeaPearl.solutionFound(model)
    end

    @testset "triggerFoundSolution!()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(2, 2, "x", trailer)
        y = SeaPearl.IntVar(3, 3, "y", trailer)

        model = SeaPearl.CPModel(trailer)

        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addObjective!(model, y)

        SeaPearl.triggerFoundSolution!(model)

        @test length(model.statistics.solutions) == 1
        @test model.statistics.solutions[1] == Dict("x" => 2,"y" => 3)
        @test model.statistics.nodevisitedpersolution[1] == 0
        @test model.objectiveBound == 2
        @test model.statistics.numberOfSolutions == 1
        @test model.statistics.objectives[1] == 3
        @test model.statistics.numberOfInfeasibleSolutions == 0
        @test model.statistics.numberOfSolutionsBeforeRestart == 1
        @test model.statistics.numberOfInfeasibleSolutionsBeforeRestart == 0
        @test model.statistics.numberOfNodesBeforeRestart == 0
    end

    @testset "tightenObjective!()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(2, 2, "x", trailer)
        y = SeaPearl.IntVar(3, 3, "y", trailer)

        model = SeaPearl.CPModel(trailer)

        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addObjective!(model, y)

        @test isnothing(model.objectiveBound)

        SeaPearl.tightenObjective!(model)

        @test model.objectiveBound == 2
    end

    @testset "belowLimits()" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        model.statistics.numberOfNodes = 1500
        model.statistics.numberOfSolutions = 15
        SeaPearl.tic()
        @test SeaPearl.belowLimits(model)

        model.limit.numberOfNodes = 1501
        model.limit.numberOfSolutions = 16
        model.limit.searchingTime = 1
        SeaPearl.tic()
        sleep(0.1)
        @test SeaPearl.belowLimits(model)

        model.statistics.numberOfNodes = 1501
        @test !SeaPearl.belowLimits(model)

        model.statistics.numberOfNodes = 1500
        model.statistics.numberOfSolutions = 16
        @test !SeaPearl.belowLimits(model)

        SeaPearl.tic()
        sleep(1)
        @test !SeaPearl.belowLimits(model)
    end

    @testset "belowNodeLimit()" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        model.statistics.numberOfNodes = 1500

        @test SeaPearl.belowNodeLimit(model)

        model.limit.numberOfNodes = 1501
        @test SeaPearl.belowNodeLimit(model)

        model.statistics.numberOfNodes = 1501
        @test !SeaPearl.belowNodeLimit(model)
    end

    @testset "belowSolutionLimits()" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        model.statistics.numberOfSolutions = 15

        @test SeaPearl.belowSolutionLimit(model)

        model.limit.numberOfSolutions = 16
        @test SeaPearl.belowSolutionLimit(model)

        model.statistics.numberOfSolutions = 16
        @test !SeaPearl.belowSolutionLimit(model)
    end

    @testset "belowTimeLimits()" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        SeaPearl.tic()
        sleep(0.1)
        @test SeaPearl.belowTimeLimit(model)

        model.limit.searchingTime = 1
        SeaPearl.tic()
        sleep(0.1)

        @test SeaPearl.belowTimeLimit(model)

        model.limit.searchingTime = 0
        SeaPearl.tic()
        sleep(0.1)
        @test !SeaPearl.belowTimeLimit(model)
    end

    @testset "Base.isempty()" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        @test isempty(model)

        x = SeaPearl.IntVar(2, 2, "x", trailer)
        SeaPearl.addVariable!(model, x)

        @test !isempty(model)
    end

    @testset "Base.empty!()" begin

        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 2, "x", trailer)
        SeaPearl.addVariable!(model, x)

        empty!(model)

        @test isempty(model)
    end

    @testset "reset_model!()" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(2, 5, "x", trailer)
        SeaPearl.addVariable!(model, x)
        SeaPearl.addObjective!(model, x)

        SeaPearl.assign!(x, 3)
        SeaPearl.fixPoint!(model)
        SeaPearl.triggerFoundSolution!(model)

        @test SeaPearl.length(x.domain) == 1
        @test model.objectiveBound == 2
        SeaPearl.reset_model!(model)
        @test SeaPearl.length(x.domain) == 4
        @test isnothing(model.objectiveBound)
    end
    @testset "restart_search" begin

        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(1, 2, "x", trailer)
        y = SeaPearl.IntVar(1, 2, "y", trailer)
        z = SeaPearl.IntVar(1, 2, "z", trailer)

        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addVariable!(model, z)
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x, y, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(y, z, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x, z, trailer))

        SeaPearl.search!(model, SeaPearl.DFSearch(), SeaPearl.MinDomainVariableSelection(), SeaPearl.BasicHeuristic())
        @test model.statistics.numberOfInfeasibleSolutionsBeforeRestart == 2
        @test model.statistics.numberOfNodesBeforeRestart == 3
        @test model.statistics.numberOfSolutionsBeforeRestart == 0

        SeaPearl.restart_search!(model)

        @test model.statistics.numberOfInfeasibleSolutionsBeforeRestart == 0
        @test model.statistics.numberOfNodesBeforeRestart == 0
        @test model.statistics.numberOfSolutionsBeforeRestart == 0
    end

    @testset "addKnownObjective!" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        SeaPearl.addKnownObjective!(model, 1)
        @test model.knownObjective == 1

    end
    @testset "triggerInfeasible!" begin

        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        x = SeaPearl.IntVar(1, 2, "x", trailer)
        y = SeaPearl.IntVar(1, 2, "y", trailer)
        z = SeaPearl.IntVar(1, 2, "z", trailer)

        SeaPearl.addVariable!(model, x)
        SeaPearl.addVariable!(model, y)
        SeaPearl.addVariable!(model, z)

        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(x, y, trailer))
        SeaPearl.addConstraint!(model, SeaPearl.NotEqual(y, z, trailer))

        @test model.statistics.infeasibleStatusPerVariable["x"] == 1
        @test model.statistics.infeasibleStatusPerVariable["y"] == 2
        @test model.statistics.infeasibleStatusPerVariable["z"] == 1

        constraint = model.constraints[1]
        SeaPearl.triggerInfeasible!(constraint, model)
        @test model.statistics.infeasibleStatusPerVariable["x"] == 2
        @test model.statistics.infeasibleStatusPerVariable["y"] == 3

    end
end

ext: jl


text:
"""
    generate_gradient(E, P, ::Type{T}, access_input)

Generate the statements for the evaluation of the polynomial with exponents `E`
and paramter exponents `P`.
This assumes that `E` and `P` are in reverse lexicographic order.
"""
function generate_gradient(E, P, ::Type{T}, access_input,coeffs=map(i -> :(c[$i]), 1:size(E, 2))) where T
    nvars = size(E,1)
    if P !== nothing
        nvars += size(P, 1)
    end

    exprs = []
    # dvals is a list of tuples (s, iszero) where s is a symbol or expression
    # representing the partial derivative and iszero a boolean flag
    # indicating whether s is zero.
    val, dvals = partial_derivatives!(exprs, E, P, T, nvars, 1, coeffs, access_input)
    out = :($(val), SVector($(Expr(:tuple, map(first, dvals)...))))

    Expr(:block, exprs..., out)
end

function partial_derivatives!(exprs, E, P, ::Type{T}, nvar, nterm, coeffs, access_input) where T
    m, n = size(E)

    # We only have one Term remaining. So we just evaluate the term and compute all
    # partial derivatives
    if n == 1
        return term(E, P, T, nterm, coeffs, access_input)
    elseif m == 1
        return univariate!(exprs, E, P, T, nvar, nterm, coeffs, access_input)
    end

    degrees, submatrices = degrees_submatrices(E)
    values = []
    dvalues = []

    for (k, E_d) in enumerate(submatrices)
        val, dvals = partial_derivatives!(exprs, E_d, P, T, nvar - 1, nterm, coeffs, access_input)
        push!(values, val)
        push!(dvalues, dvals)

        nterm += size(E_d, 2)
    end
    # Now we have to evaluate polynomials
    # for our current variable we need our new partial derivative
    val, derivative_val = evalpoly_derivative!(exprs, T, degrees, values, access_input(nvar))
    dvals = []
    reverse!(dvalues)
    for k=1:(m-1)
        coeffs_k = []
        for dv in dvalues
            dval, iszero = dv[k]
            if !isempty(coeffs_k) || !iszero
                pushfirst!(coeffs_k, dval)
            end
        end

        degrees_k = degrees[1:length(coeffs_k)]
        # we do not need assign a new variable and call evalpoly
        # if we just have a degree 0 polynomial with already computed coefficient
        if length(degrees_k) == 1 && degrees_k[1] == 0 && coeffs_k[1] isa Symbol
            push!(dvals, (coeffs_k[1], false))
        elseif length(degrees_k) > 0
            @gensym c
            push!(exprs, :($c = $(evalpoly(T, degrees_k, coeffs_k, access_input(nvar)))))
            push!(dvals, (c, false))
        else
            push!(dvals, (:(zero($T)), true))
        end
    end
    derivative_val_iszero = isempty(degrees) || (length(degrees) == 1 && degrees[1] == 0)
    push!(dvals, (derivative_val, derivative_val_iszero))

    val, dvals
end

function term(E::AbstractMatrix, P, ::Type{T}, nterm, coeffs, access_input) where T
    @assert size(E, 2) == 1
    c = term_coefficient_with_parameters(T, P, nterm, coeffs, access_input)
    offset = P === nothing ? 0 : size(P, 1)
    monomial_product_with_derivatives(T, E, c, access_input=i -> access_input(i+offset))
end

"""
    term_coefficient_with_parameters(T, P, nterm, access_input)

Compute the coefficient of the `nterm`-th term with possible parameters.
"""
function term_coefficient_with_parameters(T, P, nterm, coeffs, access_input)
    monomial_product(T, P[:,nterm], coeffs[nterm], access_input=access_input)[1]
end
term_coefficient_with_parameters(T, ::Nothing, nterm, coeffs, access_input) = coeffs[nterm]

"""
    univariate!(exprs, E::AbstractMatrix, P, ::Type{T}, nvar, nterm, access_input)

`E` represents an univariate polynomial.
"""
function univariate!(exprs, E::AbstractMatrix, P, ::Type{T}, nvar, nterm, coeffs, access_input) where T
    @assert size(E, 1) == 1
    n = size(E, 2)
    # If we have parameters it can happen that there are duplicates E
    # For each group of duplicates we have to evaluate the coefficient polynomial
    rec_coeffs, E_filtered = coefficients_with_parameters!(exprs, E, P, T, nterm, coeffs, access_input)

    val, dval = evalpoly_derivative!(exprs, T, E_filtered, rec_coeffs, access_input(nvar))
    val, [(dval, n == 1)]
end

"""
    coefficients_with_parameters!(exprs, E, P, ::Type{T}, nterm, coeffperm, access_input)

If we have parameters it can happen that there are duplicates in E.
For each group of duplicates we have to evaluate the coefficient polynomial.
If there are no parameters this is simply the coefficients vector.
"""
function coefficients_with_parameters!(exprs, E, ::Nothing, ::Type{T}, nterm, coeffs, access_input) where T
    n = size(E, 2)
    coeffs[nterm:nterm+n-1], vec(E)
end
function coefficients_with_parameters!(exprs, E, P, ::Type{T}, nterm, coeffs, access_input) where T
    n = size(E, 2)
    pcoeffs = []
    last = 1
    j = 1
    E_filtered = [E[1,1]]
    while j ≤ n
        if E[1,j] !== E_filtered[end]
            push!(E_filtered, E[1,j])
        end
        if j < n && E[1, j + 1] == E[1, j]
            j += 1
            continue
        end
        j += 1

        nterm_sub = nterm+last-1

        nduplicates = j - last - 1
        if nduplicates == 0 # nothing happened
            c, _ = monomial_product(T, P[:,nterm_sub], coeffs[nterm_sub], access_input=access_input)
            push!(pcoeffs, c)
        else
            # we have duplicates, so we consider the rest (the parameters) as an polynomial
            # and pass it to evaluate
            P_sub = @view P[:,nterm_sub:nterm+j-2]
            c = generate_evaluate!(exprs, P_sub, T, size(P, 1), nterm_sub, coeffs, access_input)
            push!(pcoeffs, c)
        end
        last = j
    end
    pcoeffs, E_filtered
end

function generate_differentiate_parameters(E, P, ::Type{T}, access_input) where T
    @assert P !== nothing
    nvars = size(E,1) + size(P, 1)

    # E and P are sorted with respect to [P; E] but since we want to derivate wrt
    # the parameters we habe to reorder them
    p = revlexicographic_cols_perm([E; P])
    E = E[:, p]
    P = P[:, p]
    coeffs = map(i -> :(c[$i]), p)
    # Since we are only interested in the derivative we can by hand eliminate all
    # terms where no parameter occurs (these are constant terms in our new setting)
    # Since we have a revlexicographic order of the columns, the all 0 column
    # is the leading columns

    # We first handle the case that P is all 0. Then the derivative is simply
    # the zero vector
    if all(iszero, P)
        zero_tuple = Expr(:tuple, (:(zero($T)) for _=1:size(P,1))...)
        return quote
            SVector($zero_tuple)
        end
    end

    if all(iszero, @view P[:, 1]) # we have at least one constant term
        nconstants = 1
        while nconstants < size(P,2)
            if all(iszero, @view P[:, nconstants + 1])
                nconstants += 1
            else
                break
            end
        end
        E = E[:,nconstants+1:end]
        P = P[:,nconstants+1:end]
        coeffs = coeffs[nconstants+1:end]
    end
    # we can derivate wrt the parameters by changing the roles
    E, P = P, E

    exprs = []
    # dvals is a list of tuples (s, iszero) where s is a symbol or expression
    # representing the partial derivative and iszero a boolean flag
    # indicating whether s is zero.
    _, dvals = partial_derivatives!(exprs, E, P, T, nvars, 1, coeffs, access_input)
    quote
        $(exprs...)
        SVector($(Expr(:tuple, map(first, dvals)...)))
    end
end


text:
using InteractiveDynamics, DynamicalBilliards, GLMakie

# For 3b1b style
# BLUE = "#7BC3DC"
# BROWN = "#8D6238"

PURPLE = "#9134F1"
TEAL =  "#25D5DB"
GREEN = "#49C544"
RED = "#C32D2D"

TEALGREEN = "#20D2A0"
GREY = "#9a9a9a"

COLORS = to_color.([PURPLE, TEALGREEN, RED, GREY])
CGRAD = cgrad(COLORS)

# Overwrite default color of obstacles to white (to fit with black background)
InteractiveDynamics.obcolor(::Obstacle) = RGBf0(1, 1, 1)
InteractiveDynamics.obfill(::Obstacle) = RGBAf0(1, 1, 1, 0.4)
InteractiveDynamics.oblw(::Obstacle) = 4

videodir(args...) = projectdir("video", args...)

MAXRES = (1500, 900) # depends on monitor
MAXPS = 1.8

include("set_black_theme.jl")

max_stars_repo_name: JuliaDynamics/ChaosThroughBilliards

max_issues_repo_name: JuliaDynamics/ChaosThroughBilliards

max_forks_repo_name: JuliaDynamics/ChaosThroughBilliards


text:
mutable struct Node
    name::String
    parents::Array{Node,1}
    children::Dict{Node,Int64}
end
Base.show(io::IO, z::Node) = print(io, z.name)

rule_pattern = r"(\S+ \S+) bags contain ((?:(?:no other|\d+ \S+ \S+) bags?(?:, |\.))+)"

rules = open("baggage_rules.txt") do f
    map(rule -> match(rule_pattern, rule).captures, readlines(f))
end

function find_or_create(name::AbstractString)::Node
    i = findfirst(b -> b.name == name, bags)
    if i === nothing
        n = Node(name, [], Dict())
        push!(bags, n)
        return n
    else
        return bags[i]
    end
end

function climb(node::Node, parents = Node[])::Array{Node,1}
    if isempty(node.parents)
        return parents
    else
        visited = Node[]
        for parent in node.parents
            path = climb(parent, push!(copy(parents), parent))
            visited = union(visited, path)
        end
        return visited
    end
end

function count(node::Node)::Int64
    if isempty(node.children)
        return 0
    else
        child_sum = 0
        for (child, cnt) in node.children
            child_sum += cnt + cnt * count(child)
        end
        return child_sum
    end
end

bags = Node[]
for r in rules
    parent = find_or_create(r[1])
    for i in eachmatch(r"(\d+) (\S+ \S+)", r[2])
        child = find_or_create(i.captures[2])
        child_count = parse(Int, i.captures[1])
        parent.children[child] = child_count
        push!(child.parents, parent)
    end
end

shiny_gold = first(filter(bag -> bag.name == "shiny gold", bags))

println(length(climb(shiny_gold)))

println(count(shiny_gold))


text:
mutable struct CCMaximinModel <: AbstractModel
    settings::Settings
    constraints::Vector{Constraint}
    obj::CCMaximinObjective
    output::Output
    CCMaximinModel(settings, constraints, obj, output) =
        new(settings, constraints, obj, output)
    CCMaximinModel() = new(Settings(), Constraint[], CCMaximinObjective(), Output())
end


text:
using CartesianGeneticProgramming
using LightGraphs
using MetaGraphs
using TikzGraphs
using TikzPictures
using LaTeXStrings
using Printf
#using Base.Test

function to_graph(c::Chromosome; active_outputs=trues(c.nout))
    actives = [n.active for n in c.nodes]
    actives[1:c.nin] = true
    vids = find(actives)
    pos = get_positions(c)
    mg = MetaDiGraph(SimpleDiGraph())
    add_vertices!(mg, length(vids)+sum(active_outputs))#c.nout)
    for i in 1:c.nin
        set_prop!(mg, i, :name, LaTeXString(string("\$in_{", i, "}\$")))
        set_prop!(mg, i, :type, 0)
    end
    for vi in (c.nin+1):length(vids)
        v = vids[vi]
        n = c.nodes[v]
        f_name = split(split(repr(n.f), ".")[end], "_")[end]
        if f_name == "const"
            set_prop!(mg, vi, :name, LaTeXString(@sprintf("%0.2f", n.p)))
        else
            set_prop!(mg, vi, :name, LaTeXString(f_name))
        end
        set_prop!(mg, vi, :function, n.f)
        set_prop!(mg, vi, :type, 2)
        set_prop!(mg, vi, :param, n.p)
        cx = findfirst(vids .== n.connections[1])
        cy = findfirst(vids .== n.connections[2])
        if cx == cy
            add_edge!(mg, Edge(cx, vi))
            set_prop!(mg, cx, vi, :ci, 3)
        else
            add_edge!(mg, Edge(cx, vi))
            set_prop!(mg, cx, vi, :ci, 1)
            add_edge!(mg, Edge(cy, vi))
            set_prop!(mg, cy, vi, :ci, 2)
        end
    end
    nid_count = 1
    for o in 1:c.nout
        if active_outputs[o]
            nid = length(vids)+nid_count
            set_prop!(mg, nid, :name, LaTeXString(string("\$out_{", o, "}\$")))
            set_prop!(mg, nid, :type, 1)
            oid = findfirst(vids .== c.outputs[o])
            add_edge!(mg, Edge(oid, nid))
            set_prop!(mg, nid, oid, :ci, 0)
            nid_count += 1
        end
    end
    mg
end

function chromo_draw(c::Chromosome, file::String="graph.pdf"; active_outputs=trues(c.nout))
    mg = to_graph(c, active_outputs=active_outputs)
    names = map(x->get_prop(mg, x, :name), 1:nv(mg))
    t = TikzGraphs.plot(mg.graph, names)
    TikzPictures.save(TikzPictures.PDF(file), t)
    nothing
end

"""
Convert a connection position x from node with position
pos to the gene equivalent (in [0.0, 1.0])
"""
function con_to_gene(x::Float64, pos::Float64)::Float64
    gene = (x - CGP.Config.input_start) / (
        (((1.0 - pos) * CGP.Config.recurrency) + pos) -
        CGP.Config.input_start)
    gene
end

function get_graph_genes(mg::MetaDiGraph, nin::Int, nout::Int)
    genes = rand(nin)
    sort!(genes, rev=true)
    node_positions = rand(nv(mg) - nout - nin)
    sort!(node_positions)
    positions = [CGP.Config.input_start .* genes; node_positions]
    outputs = Array{Float64}(nout)
    for oi in 1:nout
        vi = oi + (nv(mg) - nout)
        @assert get_prop(mg, vi, :type) == 1
        for ed in edges(mg)
            if ed.dst == vi
                outputs[oi] = positions[ed.src]
            end
        end
    end
    genes = [genes; outputs]
    for vi in (nin+1):(nv(mg)-nout)
        vprop = props(mg, vi)
        @assert vprop[:type] == 2
        ngenes = -ones(5)#Array{Float64}(5)
        ngenes[1] = positions[vi]
        for ed in edges(mg)
            if ed.dst == vi
                if get_prop(mg, ed, :ci) == 1
                    ngenes[2] = con_to_gene(positions[ed.src], positions[vi])
                end
                if get_prop(mg, ed, :ci) == 2
                    ngenes[3] = con_to_gene(positions[ed.src], positions[vi])
                end
                if get_prop(mg, ed, :ci) == 3
                    ngenes[2] = con_to_gene(positions[ed.src], positions[vi])
                    ngenes[3] = con_to_gene(positions[ed.src], positions[vi])
                end
            end
        end
        fid = findfirst(CGP.Config.functions .== vprop[:function])
        ngenes[4] = fid / length(CGP.Config.functions)
        if haskey(vprop, :param)
            ngenes[5] = (vprop[:param] + 1.0) / 2.0
        else
            ngenes[5] = rand()
        end
        if ~(all(ngenes .>= 0.0) && all(ngenes .<= 1.0))
            println(ngenes)
            println(vi)
            println(vprop)
        end
        genes = [genes; ngenes]
    end
    example = CGP.PCGPChromo(nin, nout)
    genes = [genes; rand(length(example.genes) - length(genes))]
    return genes
end

function to_chromo(mg::MetaDiGraph)::PCGPChromo
    types = [get_prop(mg, i, :type) for i in vertices(mg)]
    nin = sum(types .== 0)
    nout = sum(types .== 1)
    genes = get_graph_genes(mg, nin, nout)
    CGP.PCGPChromo(genes, nin, nout)
end

function get_graph(nin::Int, nout::Int, nnodes::Int)
    mg = MetaDiGraph(SimpleDiGraph())
    add_vertices!(mg, nin + nout + nnodes)
    types = [0*ones(nin); 2*ones(nnodes); 1*ones(nout)]
    for i in eachindex(types)
        set_prop!(mg, i, :type, types[i])
    end
    mg
end

function add_node!(mg::MetaDiGraph, i::Int, f::Function;
                   x::Int=rand(1:(i-1)), y::Int=rand(1:(i-1)), p::Float64=rand())
    set_prop!(mg, i, :function, f)
    set_prop!(mg, i, :param, p)
    if x != y
        add_edge!(mg, Edge(x, i))
        set_prop!(mg, x, i, :ci, 1)
        add_edge!(mg, Edge(y, i))
        set_prop!(mg, y, i, :ci, 2)
    else
        add_edge!(mg, Edge(x, i))
        set_prop!(mg, x, i, :ci, 3)
    end
    mg
end

function set_outputs!(mg::MetaDiGraph, nin::Int, nout::Int, nnodes::Int,
                      outputs::Array{Int})
    for i in eachindex(outputs)
        inp = outputs[i]
        out = nin + nnodes + i
        add_edge!(mg, Edge(inp, out))
        set_prop!(mg, inp, out, :ci, 0)
    end
    mg
end


text:
using LinearAlgebra

mutable struct Bird <: AbstractAgent
    id::Int
    pos::NTuple{2,Float64}
    vel::NTuple{2,Float64}
    speed::Float64
    cohere_factor::Float64
    separation::Float64
    separate_factor::Float64
    match_factor::Float64
    visual_distance::Float64
end

"""
```julia
flocking(;
    n_birds = 100,
    speed = 1.0,
    cohere_factor = 0.25,
    separation = 4.0,
    separate_factor = 0.25,
    match_factor = 0.01,
    visual_distance = 5.0,
    extent = (100, 100),
    spacing = visual_distance / 1.5
)
```
Same as in [Flock model](@ref).
"""
function flocking(;
    n_birds = 100,
    speed = 1.0,
    cohere_factor = 0.25,
    separation = 4.0,
    separate_factor = 0.25,
    match_factor = 0.01,
    visual_distance = 5.0,
    extent = (100, 100),
    spacing = visual_distance / 1.5,
)
    space2d = ContinuousSpace(extent, spacing)
    model = ABM(Bird, space2d, scheduler = random_activation)
    for _ in 1:n_birds
        vel = Tuple(rand(2) * 2 .- 1)
        add_agent!(
            model,
            vel,
            speed,
            cohere_factor,
            separation,
            separate_factor,
            match_factor,
            visual_distance,
        )
    end
    return model, flocking_agent_step!, dummystep
end

function flocking_agent_step!(bird, model)
    neighbor_ids = nearby_ids(bird, model, bird.visual_distance)
    N = 0
    match = separate = cohere = (0.0, 0.0)
    for id in neighbor_ids
        N += 1
        neighbor = model[id].pos
        heading = neighbor .- bird.pos

        cohere = cohere .+ heading
        if edistance(bird.pos, neighbor, model) < bird.separation
            separate = separate .- heading
        end
        match = match .+ model[id].vel
    end
    N = max(N, 1)
    cohere = cohere ./ N .* bird.cohere_factor
    separate = separate ./ N .* bird.separate_factor
    match = match ./ N .* bird.match_factor
    bird.vel = (bird.vel .+ cohere .+ separate .+ match) ./ 2
    bird.vel = bird.vel ./ norm(bird.vel)
    move_agent!(bird, model, bird.speed)
end


text:
abstract type DenseNode{T, N, p} <: AbstractNode{T, N, p} end

Base.IndexStyle(::Type{<: DenseNode}) = IndexLinear()

@inline function Base.getindex(x::DenseNode, i::Int)
    @boundscheck checkbounds(x, i)
    # Sometimes `i` is inactive, but corresponding data is not null
    # because mask is only set to `false` when the entry is deactivated (as long as call `cleanup!`)
    # This should be very careful because returned child node can be unexpectedly changed.
    @inbounds isnull(x) ? null(childtype(x)) : unsafe_getindex(x, i)
end

@inline function Base.setindex!(x::DenseNode, v, i::Int)
    @boundscheck checkbounds(x, i)
    # The mask is activated in MaskedDenseArray.
    # Allowing mask to be active only when setting corresponding data.
    # This ensures that when mask is active, its value is always valid
    # Similarly if null at `i`, then corresponding mask is always false (see also `delete!(node, i)`)
    isnull(v) && throw(ArgumentError("trying to set null node, call `delete!(node, i)` instead"))
    @inbounds x.data[i] = v
    x
end

@inline function Base.setindex!(x::DenseNode, ::Nothing, i::Int)
    @boundscheck checkbounds(x, i)
    @inbounds begin
        # Make mask false at i (don't delete corresponding data)
        # The children are also deactivated
        if isactive(x, i)
            x.data[i] = nothing # make `mask` `false` at `i`
            deactivate!(unsafe_getindex(x, i))
        end
    end
    x
end

@inline function deactivate!(x::DenseNode)
    isnull(x) && return x
    fill!(getmask(x), false)
    for i in eachindex(x)
        @inbounds deactivate!(unsafe_getindex(x, i))
    end
    x
end

function Base.delete!(x::DenseNode, i)
    @boundscheck checkbounds(x, i)
    @inbounds begin
        # only this case allows to set null node in `x`
        x.data[i] = nothing # make mask false
        unsafe_setindex!(x, null(childtype(x)), i)
    end
    x
end

function allocate!(x::DenseNode{T}, i::Int) where {T}
    @boundscheck checkbounds(x, i)
    @inbounds begin
        childnode = unsafe_getindex(x, i)
        if isnull(childnode) # && !isactive(x, i)
            childnode = T()
        end
        x[i] = childnode # activated in setindex!
    end
    Activated(x, i)
end

isallocated(x::DenseNode, i::Int) = unsafe_getindex(x, i) !== null(childtype(x))

struct Node{T <: AbstractNode, N, p} <: DenseNode{T, N, p}
    data::MaskedDenseArray{T, N}
    Node{T, N, p}(data) where {T, N, p} = new(data)
    Node{T, N, p}(::UndefInitializer) where {T, N, p} = new()
end

function Node{T, N, p}() where {T, N, p}
    dims = size(Node{T, N, p})
    data = MaskedDenseArray([null(T) for I in CartesianIndices(dims)])
    Node{T, N, p}(data)
end

Base.size(x::Node) = size(typeof(x))

struct DynamicNode{T <: AbstractNode, N} <: DenseNode{T, N, Dynamic()}
    data::MaskedDenseArray{T, N}
    dims::NTuple{N, Int}
    DynamicNode{T, N}(data, dims) where {T, N} = new(data, dims)
    DynamicNode{T, N}(::UndefInitializer) where {T, N} = new()
end

function DynamicNode{T, N}(dims::Int...) where {T, N}
    data = MaskedDenseArray([null(T) for I in CartesianIndices(dims)])
    DynamicNode{T, N}(data, dims)
end

Base.size(x::DynamicNode) = x.dims


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_mono_opt_rho1_9(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_mono_opt_rho1_9!(A_copy)
end

@inline function exp_mono_opt_rho1_9!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=8
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 T2k8
    # Computing Bb2 = x*I+x*A
    coeff1=0.6595448673434312
    coeff2=0.5318876318059366
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.6595448673434312
    coeff2=0.5318876318059366
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=-0.06493442247909989
    coeff2=1.1048500063419517
    coeff3=0.44991499306855187
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=0.8983444723631948
    coeff2=0.7575062720368287
    coeff3=0.1964211963869549
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=0.0035047541206341434
    coeff2=-0.07204369508927326
    coeff3=0.9994455608498845
    coeff4=0.9445089631406487
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=0.5339967206451267
    coeff2=0.9385072995625688
    coeff3=0.4235591394074441
    coeff4=0.07906900652506986
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-1.3902904224681105e-5
    coeff2=-0.03077595210974021
    coeff3=0.19354756127290104
    coeff4=0.2854511786078302
    coeff5=1.1399849812558005
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.32830291408653384
    coeff2=1.1581294214848838
    coeff3=-0.09605114401110983
    coeff4=0.03710206901207318
    coeff5=0.00037129332631147163
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-0.0002657509549628772
    coeff2=0.001194898640982042
    coeff3=-0.0030005620929433826
    coeff4=0.014125915324654426
    coeff5=0.2922686929991238
    coeff6=1.0802537414931435
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-2.124769904586236e-6
    coeff2=0.9651956304399404
    coeff3=0.2777182079782203
    coeff4=0.022365502448840103
    coeff5=0.0005202065060444869
    coeff6=1.8010644736165e-6
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing T2k8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.8942411291354833
    coeff2=0.4949663171297191
    coeff3=0.14588105730612425
    coeff4=0.20752006665209352
    coeff5=0.037061842849061506
    coeff6=0.002778168784120125
    coeff7=0.0004022829910706948
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    return memslots1 # Returning T2k8
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
# This file runs different solvers and strategies for fitting data
# subject to noise and outliers. The codes were used from known
# sources and from different languages using Julia's support to Python
# and C.
#
# The following algorithms/strategies are tested
#
# - From Scipy: Least Squares, L1, Huber and Cauchy loss functions
# - From Theia: RANSAC implementations

# In order to run Scipy tests, one should install
# --------------------------
#
# Scipy 1.3.0 and Numpy 1.17.0 in the usual Python3 environment
#
#     pip3 install scipy numpy
#
# PyCall in Julia
#
#     ENV["PYTHON"] = "/usr/bin/python3"
#     ] add PyCall

# In order to run Theia test, we have developed a C++ program that
# runs RANSAC together with CERES solver from Google. The code and
# more details are available in `../theia` directory.

using DelimitedFiles
using PyCall
using PyPlot
using RAFF
using Printf
using Random
using Base.GC

# Load libraries for drawing solutions
include("draw.jl")

"""

    run_comparative_fitting()

Run all the regular tests with one-dimensional problems and all the
solvers.

"""
function run_comparative_fitting()

    large_number = 179424673

    # First call. Just compilation.
    compare_fitting(1, [1.0 1.0 0.0; 2.0 2.0 0.0])

    for (model_str, sol) in [ ("linear", [-200.0, 1000.0]), ("cubic", [0.5, -20.0, 300.0, 1000.0]),
                              ("expon", [5000.0, 4000.0, 0.2]),
                              ("logistic", [6000.0, -5000, -0.2, -3.7]) ]

        for (np, p) in [(10, 9), (10, 8), (100, 99), (100, 90)]

            n, model, mstr = RAFF.model_list[model_str]

            # Define seed for this run. The same seed for all instances.
            Random.seed!(large_number + 300)

            generate_test_problems("/tmp/output.txt", "/tmp/sol.txt", model, mstr, n, np, p;
                                   θSol=sol, x_interval=(1.0, 30.0), std=100.0)

            cla()

            compare_fitting(model_str=model_str, MAXMS=100, theia_ms=10,
                            outimage="/tmp/$(model_str)_$(np)_$(p).png")

        end

    end

    # Run cluster tests

    for (model_str, sol) in [ ("linear", [-200.0, 1000.0]), ("cubic", [0.5, -20.0, 300.0, 1000.0]),
                              ("expon", [5000.0, 4000.0, 0.2]),
                              ("logistic", [6000.0, -5000, -0.2, -3.7]) ]

        for (np, p) in [(10, 8), (100, 90)]

            n, model, mstr = RAFF.model_list[model_str]

            # Define seed for this run. The same seed for all instances.
            Random.seed!(large_number + 300)

            generate_test_problems("/tmp/output.txt", "/tmp/sol.txt", model, mstr, n, np, p,
                                   (1.0, 30.0), (5.0, 10.0); θSol=sol, std=100.0)

            cla()

            compare_fitting(model_str=model_str, MAXMS=100, theia_ms=10,
                            outimage="/tmp/$(model_str)_$(np)_$(p)_cluster.png")

        end

    end

end

"""

    ls_measure(modl, n, data)

Compute the "Least Square" measure, used to compare different
algorithms. By definition is the square root of the sum of the squared
difference between the value predicted by the model and the actual
value only for *non-outliers*.

"""

function ls_measure(modl, n, data)

    sr = 0.0

    np, = size(data)

    for i = 1:np

        (data[i, n + 2] == 0.0) && (sr += (modl(@view(data[i, 1:n])) - data[i, n + 1])^2)

    end

    return sqrt(sr)

end

"""

    compare_fitting(;model_str="linear", kwargs...)

Run several tests from Scipy.optimize library for fitting data.

"""
function compare_fitting(N, data; kwargs...)

    open("/tmp/output.txt", "w") do fp

        @printf(fp, "%d\n", N)

        l, c = size(data)

        for i = 1:l

            for j = 1:c

                @printf(fp, "%f ", data[i, j])

            end

            @printf(fp, "\n")

        end

    end

    compare_fitting(kwargs...)

end

function compare_fitting(filename="/tmp/output.txt"; model_str="linear", MAXMS=1, outimage="/tmp/scipy.png",
                         theia_str_args="2000 -ft 0.1", theia_ms=10, kwargs...)

    n, modl, = RAFF.model_list[model_str]

    optimize = pyimport("scipy.optimize")

    # SEEDMS is an argument for RAFF, the seed for multi-start strategy
    SEEDMS = 123456789

    # PyPlot configuration
    PyPlot.rc("font", family="Helvetica")

    PyPlot.rc("font", size=10)

    # Open file
    data = Array{Float64, 2}(undef, 0, 0)

    N = 1

    open(filename) do fp

        N = parse(Int, readline(fp))

        data = readdlm(fp)

    end

    # Create objective function for Python calls
    f = let

        x = copy(data[:, 1])

        y = copy(data[:, 2])

        c = zip(x, y)

        # z = zeros(length(x))

        θ -> map(t->modl(t[1], θ) .- t[2], c)
    end

    draw_problem(data)

    t = minimum(data[:, 1]):0.01:maximum(data[:, 1])

    tmpv = @view(data[:, 3])

    open("table.txt", "a") do fp

        @printf(fp, "\\multirow{5}{*}{\$(%10s, %4d, %4d)\$} \n",
                model_str, length(tmpv), count(tmpv .== 0.0))

    end

    # Run all fitting tests from Scipy. Not using 'arctan', due to its
    # terrible results.

    for (i, (loss, line)) in enumerate(zip(
        ["linear", "soft_l1", "huber", "cauchy"],
        ["--", "-.", ":", "-"]
    ))

        seedMS = MersenneTwister(SEEDMS)

        fbest = Inf
        sbest = Nothing
        nfev  = 0

        # Eliminate effects of garbage collector
        GC.gc()

        tm = @elapsed for i = 1:MAXMS

            initguess = randn(seedMS, Float64, n)

            sol = optimize.least_squares(f, initguess, loss=loss)

            nfev += sol["nfev"]

            if sol["cost"] < fbest

                fbest = sol["cost"]
                sbest = sol

            end

        end

        # This is a tentative of comparing different algorithms.  We
        # compute the sum of the squared differences only for
        # NON-outliers.

        modl2 = (x) -> modl(x, sbest["x"])

        fmeas = ls_measure(modl2, N, data)

        open("table.txt", "a") do fp

            @printf(fp, "  & %10s & %10.3e & %8.4f & %8d & ", loss, fmeas, tm, nfev)

            for j = 1:n

                @printf(fp, "\$%15.5f\$, ", sbest["x"][j])

            end

            @printf(fp, " \\\\\n")

        end

        PyPlot.plot(t, modl2.(t), color=PyPlot.cm."Set1"(i/9.0),
                    linestyle=line, label=loss)

    end

    # Run Theia

    fname = "/tmp/" * Random.randstring(10) * ".txt"

    run(pipeline(`../theia/ransac $(model_str) $(theia_str_args) -ms $(theia_ms) -f $filename`, fname))

    theia = DelimitedFiles.readdlm(fname)

    tt = minimum(data[:, 1]):1:maximum(data[:, 1])

    tm, = size(theia)

    for i = 1:Int(tm / theia_ms)

        b, bpos = findmin(@view(theia[(i - 1) * theia_ms + 1:i * theia_ms, 4]))

        total_time = sum(@view(theia[(i - 1) * theia_ms + 1:i * theia_ms, 2]))

        tvec = @view theia[(i - 1) * theia_ms + bpos, :]

        modl2 = (x) -> modl(x, tvec[end - n + 1:end])

        fmeas = ls_measure(modl2, N, data)

        open("table.txt", "a") do fp

            @printf(fp, "  & %10s & %10.3e & %8.4f & %8d & ", tvec[1], fmeas, total_time, 0)

            for j = 1:n

                @printf(fp, "\$%15.5f\$, ", tvec[end - n + j])

            end

            @printf(fp, " \\\\\n")

        end

        PyPlot.plot(tt, modl2.(tt), color=PyPlot.cm."Set1"((5 + i)/9.0),
                    linestyle="-", marker=3 + i, label=tvec[1])

    end

    # Run RAFF

    initguess = zeros(Float64, n)

    # Eliminate effects of garbage collector
    GC.gc()

    rsol, tm = @timed praff(modl, data[:, 1:end - 1], n; kwargs..., initguess=initguess,
                            SEEDMS=SEEDMS, MAXMS=MAXMS)

    modl2 = (x) -> modl(x, rsol.solution)

    fmeas = ls_measure(modl2, N, data)

    open("table.txt", "a") do fp

        @printf(fp, "  & %10s & %10.3e & %8.4f & %8d & ", "RAFF.jl", fmeas, tm, rsol.nf)

        for j = 1:n

            @printf(fp, "\$%15.5f\$, ", rsol.solution[j])

        end

        @printf(fp, " \\\\\n")

    end

    PyPlot.plot(t, modl2.(t), color=PyPlot.cm."Set1"(2.0/9.0),
                linestyle="-", label="RAFF")

    PyPlot.legend(loc="best")

    PyPlot.show()

    PyPlot.savefig(outimage, DPI=600, bbox_inches="tight")

    PyPlot.savefig("/tmp/scipy.png", DPI=150)

end

function compare_circle_fitting(filename="/tmp/output.txt";MAXMS=1, outimage="/tmp/scipy.png",
                                theia_str_args="10 -ft 0.1", theia_ms=10, kwargs...)

    n, modl, = RAFF.model_list["circle"]

    optimize = pyimport("scipy.optimize")

    # SEEDMS is an argument for RAFF, the seed for multi-start strategy
    SEEDMS = 123456789

    # PyPlot configuration
    PyPlot.rc("font", family="Helvetica")

    PyPlot.rc("font", size=10)

    # Open file
    data = Array{Float64, 2}(undef, 0, 0)

    N = 1

    open(filename) do fp

        N = parse(Int, readline(fp))

        data = readdlm(fp)

    end

    f = let

        l, = size(data)

        x = copy(data[:, 1:2])

        x = transpose(x)

        θ -> [modl(@view(x[:, i]), θ) for i = 1:l]
    end

    draw_circle_sol(data)

    t = 0:0.1:2.1 * π

    ptx = (α, ρ, d) -> ρ * cos(α) + d[1]
    pty = (α, ρ, d) -> ρ * sin(α) + d[2]

    tmpv = @view(data[:, 4])

    @printf("\\multirow{5}{*}{\$(%10s, %4d, %4d)\$} \n", "circle", length(tmpv), count(tmpv .== 0.0))

    # Run all fitting tests from Scipy. Not using 'arctan', due to its
    # terrible results.

    for (i, (loss, line)) in enumerate(zip(
        ["linear", "soft_l1", "huber", "cauchy"],
        ["--", "-.", ":", "-"]
    ))

        initguess = ones(Float64, n)

        seedMS = MersenneTwister(SEEDMS)

        fbest = Inf
        sbest = Nothing
        nfev  = 0

        tm = @elapsed for i = 1:MAXMS

            initguess = 1.0 .+ randn(seedMS, Float64, n)

            sol = optimize.least_squares(f, initguess, loss=loss)

            nfev += sol["nfev"]

            if sol["cost"] < fbest

                fbest = sol["cost"]
                sbest = sol

            end

        end

        # This is a tentative of comparing different algorithms.  We
        # compute the sum of the squared differences only for
        # NON-outliers.

        modl2 = (x) -> modl(x, sbest["x"])

        fmeas = ls_measure(modl2, N, data)

        @printf("  & %10s & %10.3e & %8.4f & %8d & ", loss, fmeas, tm, nfev)

        for j = 1:n

            @printf("\$%15.5f\$, ", sbest["x"][j])

        end

        @printf(" \\\\\n")

        modl1x = (α) -> ptx(α, sbest["x"][3], sbest["x"][1:2])
        modl1y = (α) -> pty(α, sbest["x"][3], sbest["x"][1:2])

        PyPlot.plot(modl1x.(t), modl1y.(t), color=PyPlot.cm."Set1"(i/9.0),
                    linestyle=line, label=loss)

    end

    # Run Theia

    fname = "/tmp/" * Random.randstring(10) * ".txt"

    run(pipeline(`../theia/ransac circle $(theia_str_args) -ms $(theia_ms)`, fname))

    theia = DelimitedFiles.readdlm(fname)

    tt = 0:0.2:2.1 * π

    tm, = size(theia)

    for i = 1:Int(tm / theia_ms)

        b, bpos = findmin(@view(theia[(i - 1) * theia_ms + 1:i * theia_ms, 4]))

        total_time = sum(@view(theia[(i - 1) * theia_ms + 1:i * theia_ms, 2]))

        tvec = @view theia[(i - 1) * theia_ms + bpos, :]

        modl2 = (x) -> modl(x, tvec[end - n + 1:end])

        fmeas = ls_measure(modl2, N, data)

        @printf("  & %10s & %10.3e & %8.4f & %8d & ", tvec[1], fmeas, total_time, 0)

        for j = 1:n

            @printf("\$%15.5f\$, ", tvec[end - n + j])

        end

        @printf(" \\\\\n")

        modl1x = (α) -> ptx(α, tvec[end], tvec[end - 2:end - 1])
        modl1y = (α) -> pty(α, tvec[end], tvec[end - 2:end - 1])

        PyPlot.plot(modl1x.(tt), modl1y.(tt), color=PyPlot.cm."Set1"((5 + i)/9.0),
                    linestyle="-", marker=i, label=tvec[1])

    end

    # Run RAFF

    initguess = ones(Float64, n)

    rsol, tm = @timed praff(modl, data[:, 1:end - 1], n; kwargs..., MAXMS=MAXMS,
                            initguess=initguess, SEEDMS=SEEDMS)

    fSol = rsol.solution

    modl2 = (x) -> modl(x, rsol.solution)

    fmeas = ls_measure(modl2, N, data)

    @printf("  & %10s & %10.3e & %8.4f & %8d & ", "RAFF.jl", fmeas, tm, rsol.nf)

    for j = 1:n

        @printf("\$%15.5f\$, ", fSol[j])

    end

    @printf(" \\\\\n")

    modl1x = (α) -> ptx(α, fSol[3], fSol[1:2])
    modl1y = (α) -> pty(α, fSol[3], fSol[1:2])

    PyPlot.plot(modl1x.(t), modl1y.(t), color=PyPlot.cm."Set1"(2.0/9.0),
                linestyle="-", label="RAFF")

    PyPlot.legend(bbox_to_anchor=(0., 1.01, 1., .111), loc="lower left",
           ncol=3, mode="expand", borderaxespad=0., numpoints=1)

    PyPlot.show()

    PyPlot.savefig(outimage, DPI=600, bbox_inches="tight")

    PyPlot.savefig("/tmp/scipy.png", DPI=150)

end

function run_only_praff_circle(filename="/tmp/output.txt"; MAXMS=1, kwargs...)

    n, modl, = RAFF.model_list["circle"]

    # SEEDMS is an argument for RAFF, the seed for multi-start strategy
    SEEDMS = 123456789

    # Open file
    data = Array{Float64, 2}(undef, 0, 0)

    N = 1

    open(filename) do fp

        N = parse(Int, readline(fp))

        data = readdlm(fp)

    end

    initguess = ones(Float64, n)

    rsol, tm = @timed praff(modl, data[:, 1:end - 1], n; kwargs..., MAXMS=MAXMS,
                            initguess=initguess, SEEDMS=SEEDMS)

    fSol = rsol.solution

    modl2 = (x) -> modl(x, rsol.solution)

    fmeas = ls_measure(modl2, N, data)

    @printf("  & %10s & %10.3e & %8.4f & %8d & ", "RAFF.jl", fmeas, tm, rsol.nf)

    for j = 1:n

        @printf("\$%15.5f\$, ", fSol[j])

    end

    @printf(" \\\\\n")

end

"""

    run_circle_comparative_fitting()

Run all the regular tests with one-dimensional problems and all the
solvers.

"""
function run_circle_comparative_fitting()

    large_number = 179424673

    # First call. Just compilation.
    data, v = gen_circle(10, 10, θSol=[1.0, 1.0, 1.0])

    compare_circle_fitting(MAXMS=1, theia_ms=1)

    np = 100

    n, model, mstr = RAFF.model_list["circle"]

    for f = 0.1:0.1:0.9

        # Define seed for this run. The same seed for all instances.
        Random.seed!(large_number + 300)

        p = Int(round(f * np))

        data, v = gen_circle(np, p, θSol=[- 10.0, 30.0, 2.0],
                             std=0.1, outTimes=20.0)

        cla()

        compare_circle_fitting(MAXMS=100, theia_ms=10, ftrusted=min(f, 0.5),
                               outimage="/tmp/circle_$(np)_$(p).png")

        run_only_praff_circle(MAXMS=100,
                              ftrusted=(max(0.0, f - 0.3), min(1.0, f + 0.3)))

    end

    np = 300

    for f = 0.1:0.1:0.9

        # Define seed for this run. The same seed for all instances.
        Random.seed!(large_number + 300)

        p = Int(round(f * np))

        data, v = gen_ncircle(np, p, θSol=[- 10.0, 30.0, 2.0],
                              std=0.1)

        cla()

        compare_circle_fitting(MAXMS=100, theia_ms=10, ftrusted=min(f, 0.5),
                               outimage="/tmp/ncircle_$(np)_$(p).png")

        run_only_praff_circle(MAXMS=100,
                              ftrusted=(max(0.0, f - 0.3), min(1.0, f + 0.3)))

    end

end

function praff_versus_gauss_newton(P; filename="/tmp/output.txt", model_str="circle",
                                   MAXMS=10, MAXMSGN=100, outimage="/tmp/gn.png", kwargs...)

    n, modl, = RAFF.model_list[model_str]

    # Open file
    data = Array{Float64, 2}(undef, 0, 0)

    # t = 0:0.1:2.1 * π

    # ptx = (α, ρ, d) -> ρ * cos(α) + d[1]
    # pty = (α, ρ, d) -> ρ * sin(α) + d[2]

    N = 1

    open(filename) do fp

        N = parse(Int, readline(fp))

        data = readdlm(fp)

    end

    # draw_circle_sol(data)

    @printf("%12s ", "")

    for p in P

        @printf("& %10.3e ", p)

    end

    @printf("\\\\\n%12s ", "\\raff")

    fSol = RAFFOutput()

    for p in P

        # SEEDMS is an argument for RAFF, the seed for multi-start strategy
        SEEDMS = 123456789

        initguess = ones(Float64, n)

        rsol, tm = @timed raff(modl, data[:, 1:end - 1], n; kwargs..., MAXMS=MAXMS,
                               ftrusted=(max(0.0, p - 0.3), min(1.0, p + 0.3)),
                               initguess=initguess, SEEDMS=SEEDMS, inner_solver=lmlovo)

        fSol = rsol.solution

        modl2 = (x) -> modl(x, rsol.solution)

        fmeas = ls_measure(modl2, N, data)

        @printf("& %10.1e ", fmeas)

    end

    @printf("\\\\\n")

    # modl1x = (α) -> ptx(α, fSol[3], fSol[1:2])
    # modl1y = (α) -> pty(α, fSol[3], fSol[1:2])

    # PyPlot.plot(modl1x.(t), modl1y.(t), color=PyPlot.cm."Set1"(2.0/9.0),
    #             linestyle="-", label="RAFF")

    @printf("%12s ", "\\texttt{GN}")

    for p in P

        # SEEDMS is an argument for RAFF, the seed for multi-start strategy
        SEEDMS = 123456789

        initguess = ones(Float64, n)

        rsol, tm = @timed raff(modl, data[:, 1:end - 1], n; kwargs...,
                               ftrusted=(p, p), MAXMS=MAXMSGN,
                               initguess=initguess, SEEDMS=SEEDMS,
                               inner_solver=gnlslovo)

        fSol = rsol.solution

        modl2 = (x) -> modl(x, rsol.solution)

        fmeas = ls_measure(modl2, N, data)

        @printf("& %10.1e ", fmeas)

    end

    @printf("\\\\\n")

    # modl1x = (α) -> ptx(α, fSol[3], fSol[1:2])
    # modl1y = (α) -> pty(α, fSol[3], fSol[1:2])

    # PyPlot.plot(modl1x.(t), modl1y.(t), color=PyPlot.cm."Set1"(3.0/9.0),
    #             linestyle="-", label="GN")

    # PyPlot.legend(bbox_to_anchor=(0., 1.01, 1., .111), loc="lower left",
    #               ncol=3, mode="expand", borderaxespad=0., numpoints=1)

    # PyPlot.show()

    # PyPlot.savefig(outimage, DPI=600, bbox_inches="tight")

    # PyPlot.savefig("/tmp/gn.png", DPI=150)

end

function run_praff_vs_gn()

    large_number = 179424673

    # Run cluster tests

    for (model_str, sol) in [ ("linear", [-200.0, 1000.0]), ("cubic", [0.5, -20.0, 300.0, 1000.0]),
                              ("expon", [5000.0, 4000.0, 0.2]),
                              ("logistic", [6000.0, -5000, -0.2, -3.7]) ]

        for (np, p) in [(100, 90)]

            n, model, mstr = RAFF.model_list[model_str]

            # Define seed for this run. The same seed for all instances.
            Random.seed!(large_number + 300)

            generate_test_problems("/tmp/output.txt", "/tmp/sol.txt", model, mstr, n, np, p,
                                   (1.0, 30.0), (5.0, 10.0); θSol=sol, std=100.0)

            praff_versus_gauss_newton([0.1, 0.2, 0.5, 0.6, 0.7, 0.95, 0.98, 0.99]; model_str=model_str, MAXMS=30)

        end

    end

    # Circle

    generate_ncircle("/tmp/output.txt", 300, 100; θSol=[- 10.0, 30.0, 2.0],
                     std=0.1, leftd=5, lngth=7)

    praff_versus_gauss_newton([0.01, 0.02, 0.05, 0.1, 0.4, 0.5, 0.6, 0.7]; MAXMS=30)

end


text:
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("colorscheme.jl"))

using InteractiveDynamics
using DynamicalSystems, OrdinaryDiffEq
import GLMakie

function rossler_driving(u, p, t)
    k = p[1]
    x_1, x_2, x_3, y_1, y_2, y_3, z_1, z_2, z_3 = u

    dx_1  =  - x_2 -x_3
    dx_2  =   x_1 + 0.2x_2
    dx_3  =  0.2 + x_3*(x_1-5.7)
    (dx_1, dx_2, dx_3) = 6 .* (dx_1, dx_2, dx_3)

    dy_1  =  10*(- y_1 + y_2)
    dy_2  =  28y_1 - y_2 - y_1*y_3 + k*x_2
    dy_3  =  y_1*y_2 - 2.666y_3

    dz_1  =  10*(- z_1 + z_2)
    dz_2  =  28z_1 - z_2 - z_1*z_3 + k*x_2
    dz_3  =  z_1*z_2 - 2.666z_3

    return SVector(dx_1, dx_2, dx_3, dy_1, dy_2, dy_3, dz_1, dz_2, dz_3)
end

x0 = [-1, -2, 0.1]
y0 = [5.0, 5.0, 20.0]
z0 = [-5.0, -5.0, 20.0]

u0 = vcat(x0, y0, z0)
k = 5.0

ds = ContinuousDynamicalSystem(rossler_driving, u0, [k])

ps = Dict(
    1 => 0.0:0.01:12,
)
pnames = Dict(
    1 => "k",
)

u0s =  [ds.u0]

fig, obs = interactive_evolution_timeseries(
    ds, u0s, ps; tail = 500, idxs = [4,7], pnames,
    lims = ((-20,20), (-20,20))
)

GLMakie.Label(fig[:, :, GLMakie.Top()], "Roessler subsystem driving two Lorenz subsystems: y, z. Plotted are the first variables of each subsystem. See Eq.(9.9)", valign = :bottom, padding = (0, 0, 10, 0))
GLMakie.content(fig[1,1]).xlabel = "y₁"
GLMakie.content(fig[1,1]).ylabel = "z₁"
GLMakie.content(fig[:, 2][1,1]).ylabel = "y₁"
GLMakie.content(fig[:, 2][2,1]).ylabel = "z₁"
GLMakie.content(fig[:, 2][2,1]).xlabel = "t"

# %%
record_interaction(fig, string(@__FILE__)[1:end-2]*"mp4"; total_time = 20, sleep_time = 2)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
module dev

using RMB
import RMB: inventorymodel, pidc1, pidc2, detensor, explicitinv

include("RMBplots")

# inventorymodel((dd,dc,ds),pricelevel,spotbuy,wholesaleprices)
# dims=(5,3,10)

# m1=inventorymodel(dims,5.0,20.0,[0.1,1.0,1.1])
# pi!(m1,pidc1(m1))
# # pid,pic=detensor(m1)   #for debugging
# u!(m1,[0,0.0])    #free stocking, free order
# solvedp!(m1)
# ka=rand(m1,200)
# tsplot(m1,ka)

# m2=inventorymodel(dims,5.0,20.0,[0.8,1.0,1.2])
# pi!(m2,pidc1(m2))
# u!(m2,[0.0,1.0])    #free stocking, expensive order
# solvedp!(m2)
# ka=rand(m2,200)
# tsplot(m2,ka)

# m3=inventorymodel(dims,8.0,20.0,[0.90,0.9,0.9])
# pi!(m3,pidc2(m3))
# u!(m3,[0.05,0.1])    #free stocking, free order
# solvedp!(m3)
# ka=rand(m3,300)
# tsplot(m3,ka)

end #module


text:
let
    trajdata = get_test_trajdata()

    Δt = 0.1
    rec = SceneRecord(5, Δt)
    @test record_length(rec) == 5
    @test length(rec) == 0
    @test !pastframe_inbounds(rec, 0)
    @test !pastframe_inbounds(rec, -1)
    @test !pastframe_inbounds(rec, 1)

    scene = get!(Scene(), trajdata, 1)
    update!(rec, scene)
    @test length(rec) == 1
    @test pastframe_inbounds(rec, 0)
    @test !pastframe_inbounds(rec, -1)
    @test !pastframe_inbounds(rec, 1)
    @test isapprox(get_elapsed_time(rec, 0), Δt)
    @test rec[1,0].state == get_vehiclestate(trajdata, 1, 1)
    @test rec[1,0].def == get_vehicledef(trajdata, 1)
    @test rec[2,0].state == get_vehiclestate(trajdata, 2, 1)
    @test rec[2,0].def == get_vehicledef(trajdata, 2)
    show(IOBuffer(), rec)

    get!(scene, trajdata, 2)
    update!(rec, scene)
    @test length(rec) == 2
    @test pastframe_inbounds(rec, 0)
    @test pastframe_inbounds(rec, -1)
    @test !pastframe_inbounds(rec, 1)
    @test isapprox(get_elapsed_time(rec,  0),  Δt)
    @test isapprox(get_elapsed_time(rec, -1), 2Δt)
    @test isapprox(get_elapsed_time(rec, -1, 0), Δt)
    @test rec[1,0].state == get_vehiclestate(trajdata, 1, 2)
    @test rec[1,0].def == get_vehicledef(trajdata, 1)
    @test rec[2,0].state == get_vehiclestate(trajdata, 2, 2)
    @test rec[2,0].def == get_vehicledef(trajdata, 2)
    @test rec[1,-1].state == get_vehiclestate(trajdata, 1, 1)
    @test rec[1,-1].def == get_vehicledef(trajdata, 1)
    @test rec[2,-1].state == get_vehiclestate(trajdata, 2, 1)
    @test rec[2,-1].def == get_vehicledef(trajdata, 2)

    @test  iscarinframe(rec, 1, 0)
    @test !iscarinframe(rec, 10, 0)

    @test get_index_of_first_vehicle_with_id(rec, 1, 0) == 1
    @test get_index_of_first_vehicle_with_id(rec, 2, 0) == 2
    @test get_index_of_first_vehicle_with_id(rec, 5, 0) == 0

    scene2 = get!(Scene(), rec)
    @test scene2[1].state == get_vehiclestate(trajdata, 1, 2)
    @test scene2[1].def == get_vehicledef(trajdata, 1)
    @test scene2[2].state == get_vehiclestate(trajdata, 2, 2)
    @test scene2[2].def == get_vehicledef(trajdata, 2)

    get!(scene2, rec, -1)
    @test scene2[1].state == get_vehiclestate(trajdata, 1, 1)
    @test scene2[1].def == get_vehicledef(trajdata, 1)
    @test scene2[2].state == get_vehiclestate(trajdata, 2, 1)
    @test scene2[2].def == get_vehicledef(trajdata, 2)

    empty!(rec)
    @test length(rec) == 0

    test_veh_state = VehicleState(VecSE2(7.0,7.0,2.0), trajdata.roadway, 10.0)
    test_veh_def = VehicleDef(999, AgentClass.CAR, 5.0, 3.0)
    test_veh = Vehicle(test_veh_state, test_veh_def)
    rec[1,-1] = test_veh
    @test get_vehiclestate(rec,999,-1) == test_veh_state
end


text:
# + is not defined between SingleVariable
function MP.polynomial(p::GramMatrix{MOI.SingleVariable})
    Q = convert(Vector{MOI.ScalarAffineFunction{Float64}}, p.Q.Q)
    MP.polynomial(GramMatrix(SymMatrix(Q, p.Q.n), p.x))
end
function MP.polynomial(p::GramMatrix{F}) where {F <: MOI.AbstractFunction}
    MP.polynomial(p, MOIU.promote_operation(+, Float64, F, F))
end

function primal_value(model, p::GramMatrix{MOI.SingleVariable})
    # TODO [perf] use MOI typed mapped array
    Q = MOI.get(model, MOI.VariablePrimal(),
                MOI.VariableIndex[sv.variable for sv in p.Q.Q])
    return GramMatrix(SymMatrix(Q, p.Q.n), p.x)
end

function add_matrix_variable_bridge(
    model::MOI.ModelLike, MCT,
    side_dimension::Integer, T::Type)
    mat_cone = matrix_cone(MCT, side_dimension)
    VB = variable_bridge_type(typeof(mat_cone), T)
    return add_variable_bridge(VB, model, mat_cone)
end
function union_vector_bridge_types(
    MCT, T::Type)
    return Union{variable_bridge_type(typeof(matrix_cone(MCT, 0)), T),
                 variable_bridge_type(typeof(matrix_cone(MCT, 1)), T),
                 variable_bridge_type(typeof(matrix_cone(MCT, 2)), T),
                 variable_bridge_type(typeof(matrix_cone(MCT, 3)), T)}
end
function append_added_constraint_types(
    added, MCT, T::Type)
    append!(added, MOIB.added_constraint_types(variable_bridge_type(typeof(matrix_cone(MCT, 0)), T)))
    append!(added, MOIB.added_constraint_types(variable_bridge_type(typeof(matrix_cone(MCT, 1)), T)))
    append!(added, MOIB.added_constraint_types(variable_bridge_type(typeof(matrix_cone(MCT, 2)), T)))
    append!(added, MOIB.added_constraint_types(variable_bridge_type(typeof(matrix_cone(MCT, 3)), T)))
    return added
end


text:
# Defining all the First order derivatives for the immersed boundaries

@inline conditional_x_derivative_f(LY, LZ, i, j, k, ibg::IBG{FT}, deriv, args...) where FT = ifelse(inactive_node(c, LY, LZ, i, j, k, ibg) | inactive_node(c, LY, LZ, i-1, j, k, ibg), zero(FT), deriv(i, j, k, ibg.underlying_grid, args...))
@inline conditional_x_derivative_c(LY, LZ, i, j, k, ibg::IBG{FT}, deriv, args...) where FT = ifelse(inactive_node(f, LY, LZ, i, j, k, ibg) | inactive_node(f, LY, LZ, i+1, j, k, ibg), zero(FT), deriv(i, j, k, ibg.underlying_grid, args...))
@inline conditional_y_derivative_f(LX, LZ, i, j, k, ibg::IBG{FT}, deriv, args...) where FT = ifelse(inactive_node(LX, c, LZ, i, j, k, ibg) | inactive_node(LX, c, LZ, i, j-1, k, ibg), zero(FT), deriv(i, j, k, ibg.underlying_grid, args...))
@inline conditional_y_derivative_c(LX, LZ, i, j, k, ibg::IBG{FT}, deriv, args...) where FT = ifelse(inactive_node(LX, f, LZ, i, j, k, ibg) | inactive_node(LX, f, LZ, i, j+1, k, ibg), zero(FT), deriv(i, j, k, ibg.underlying_grid, args...))
@inline conditional_z_derivative_f(LX, LY, i, j, k, ibg::IBG{FT}, deriv, args...) where FT = ifelse(inactive_node(LX, LY, c, i, j, k, ibg) | inactive_node(LX, LY, c, i, j, k-1, ibg), zero(FT), deriv(i, j, k, ibg.underlying_grid, args...))
@inline conditional_z_derivative_c(LX, LY, i, j, k, ibg::IBG{FT}, deriv, args...) where FT = ifelse(inactive_node(LX, LY, f, i, j, k, ibg) | inactive_node(LX, LY, f, i, j, k+1, ibg), zero(FT), deriv(i, j, k, ibg.underlying_grid, args...))

∂xᶠᶜᶜ(i, j, k, ibg::IBG, args...) = conditional_x_derivative_f(c, c, i, j, k, ibg, ∂xᶠᶜᶜ, args...)
∂xᶠᶜᶠ(i, j, k, ibg::IBG, args...) = conditional_x_derivative_f(c, f, i, j, k, ibg, ∂xᶠᶜᶠ, args...)
∂xᶠᶠᶜ(i, j, k, ibg::IBG, args...) = conditional_x_derivative_f(f, c, i, j, k, ibg, ∂xᶠᶠᶜ, args...)
∂xᶠᶠᶠ(i, j, k, ibg::IBG, args...) = conditional_x_derivative_f(f, f, i, j, k, ibg, ∂xᶠᶠᶠ, args...)

∂xᶜᶜᶜ(i, j, k, ibg::IBG, args...) = conditional_x_derivative_c(c, c, i, j, k, ibg, ∂xᶜᶜᶜ, args...)
∂xᶜᶜᶠ(i, j, k, ibg::IBG, args...) = conditional_x_derivative_c(c, f, i, j, k, ibg, ∂xᶜᶜᶠ, args...)
∂xᶜᶠᶜ(i, j, k, ibg::IBG, args...) = conditional_x_derivative_c(f, c, i, j, k, ibg, ∂xᶜᶠᶜ, args...)
∂xᶜᶠᶠ(i, j, k, ibg::IBG, args...) = conditional_x_derivative_c(f, f, i, j, k, ibg, ∂xᶜᶠᶠ, args...)

∂yᶜᶠᶜ(i, j, k, ibg::IBG, args...) = conditional_y_derivative_f(c, c, i, j, k, ibg, ∂yᶜᶠᶜ, args...)
∂yᶜᶠᶠ(i, j, k, ibg::IBG, args...) = conditional_y_derivative_f(c, f, i, j, k, ibg, ∂yᶜᶠᶠ, args...)
∂yᶠᶠᶜ(i, j, k, ibg::IBG, args...) = conditional_y_derivative_f(f, c, i, j, k, ibg, ∂yᶠᶠᶜ, args...)
∂yᶠᶠᶠ(i, j, k, ibg::IBG, args...) = conditional_y_derivative_f(f, f, i, j, k, ibg, ∂yᶠᶠᶠ, args...)

∂yᶜᶜᶜ(i, j, k, ibg::IBG, args...) = conditional_y_derivative_c(c, c, i, j, k, ibg, ∂yᶜᶜᶜ, args...)
∂yᶜᶜᶠ(i, j, k, ibg::IBG, args...) = conditional_y_derivative_c(c, f, i, j, k, ibg, ∂yᶜᶜᶠ, args...)
∂yᶠᶜᶜ(i, j, k, ibg::IBG, args...) = conditional_y_derivative_c(f, c, i, j, k, ibg, ∂yᶠᶜᶜ, args...)
∂yᶠᶜᶠ(i, j, k, ibg::IBG, args...) = conditional_y_derivative_c(f, f, i, j, k, ibg, ∂yᶠᶜᶠ, args...)

∂zᶜᶜᶠ(i, j, k, ibg::IBG, args...) = conditional_z_derivative_f(c, c, i, j, k, ibg, ∂zᶜᶜᶠ, args...)
∂zᶜᶠᶠ(i, j, k, ibg::IBG, args...) = conditional_z_derivative_f(c, f, i, j, k, ibg, ∂zᶜᶠᶠ, args...)
∂zᶠᶜᶠ(i, j, k, ibg::IBG, args...) = conditional_z_derivative_f(f, c, i, j, k, ibg, ∂zᶠᶜᶠ, args...)
∂zᶠᶠᶠ(i, j, k, ibg::IBG, args...) = conditional_z_derivative_f(f, f, i, j, k, ibg, ∂zᶠᶠᶠ, args...)

∂zᶜᶜᶜ(i, j, k, ibg::IBG, args...) = conditional_z_derivative_c(c, c, i, j, k, ibg, ∂zᶜᶜᶜ, args...)
∂zᶜᶠᶜ(i, j, k, ibg::IBG, args...) = conditional_z_derivative_c(c, f, i, j, k, ibg, ∂zᶜᶠᶜ, args...)
∂zᶠᶜᶜ(i, j, k, ibg::IBG, args...) = conditional_z_derivative_c(f, c, i, j, k, ibg, ∂zᶠᶜᶜ, args...)
∂zᶠᶠᶜ(i, j, k, ibg::IBG, args...) = conditional_z_derivative_c(f, f, i, j, k, ibg, ∂zᶠᶠᶜ, args...)

using  Oceananigans.Operators
import Oceananigans.Operators: Γᶠᶠᶜ

# Circulation equal to zero on a solid nodes
@inline Γᶠᶠᶜ(i, j, k, ibg::IBG, u, v) =
    conditional_x_derivative_f(f, c, i, j, k, ibg, δxᶠᵃᵃ, Δy_qᶜᶠᶜ, v) - conditional_y_derivative_f(f, c, i, j, k, ibg, δyᵃᶠᵃ, Δx_qᶠᶜᶜ, u)