text:
@testset "upbranch and downbranch" begin
    fm = _build_dmip_formulation()
    node = Cerberus.Node()
    let db = @inferred Cerberus.down_branch(node, _CVI(1), 0.5)
        @test db.lt_bounds == [Cerberus.BoundUpdate(_CVI(1), _LT(0.0))]
        @test isempty(db.gt_bounds)
        @test db.depth == 1
    end

    let ub = @inferred Cerberus.up_branch(node, _CVI(2), 0.5)
        @test isempty(ub.lt_bounds)
        @test ub.gt_bounds == [Cerberus.BoundUpdate(_CVI(2), _GT(1.0))]
        @test ub.depth == 1
    end

    # Now do general integer branchings
    let db = @inferred Cerberus.down_branch(node, _CVI(1), 3.7)
        @test db.lt_bounds == [Cerberus.BoundUpdate(_CVI(1), _LT(3.0))]
        @test isempty(db.gt_bounds)
        @test db.depth == 1
    end

    let ub = @inferred Cerberus.up_branch(node, _CVI(2), 3.7)
        @test isempty(ub.lt_bounds)
        @test ub.gt_bounds == [Cerberus.BoundUpdate(_CVI(2), _GT(4.0))]
        @test ub.depth == 1
    end
end

@testset "apply_branching!" begin
    fm = _build_dmip_formulation()
    node = Cerberus.Node()
    let f = _CVI(2), s = _LT(1.0)
        bd = Cerberus.BoundUpdate(f, s)
        @inferred Cerberus.apply_branching!(node, bd)
        @test node.lt_bounds == [Cerberus.BoundUpdate(_CVI(2), _LT(1.0))]
        @test isempty(node.gt_bounds)
        @test isempty(node.lt_general_constrs)
        @test isempty(node.gt_general_constrs)
        @test node.depth == 1
        @test node.dual_bound == -Inf
    end
    let f = _CVI(4), s = _GT(3.0)
        bd = Cerberus.BoundUpdate(f, s)
        @inferred Cerberus.apply_branching!(node, bd)
        @test node.lt_bounds == [Cerberus.BoundUpdate(_CVI(2), _LT(1.0))]
        @test node.gt_bounds == [Cerberus.BoundUpdate(_CVI(4), _GT(3.0))]
        @test isempty(node.lt_general_constrs)
        @test isempty(node.gt_general_constrs)
        @test node.depth == 2
        @test node.dual_bound == -Inf
    end
    let f = _CSAF([1.2, 3.4], [_CVI(1), _CVI(3)], 5.6), s = _LT(7.8)
        bd = Cerberus.AffineConstraint(f, s)
        @inferred Cerberus.apply_branching!(node, bd)
        @test node.lt_bounds == [Cerberus.BoundUpdate(_CVI(2), _LT(1.0))]
        @test node.gt_bounds == [Cerberus.BoundUpdate(_CVI(4), _GT(3.0))]
        @test length(node.lt_general_constrs) == 1
        lt_constr = node.lt_general_constrs[1]
        @test _is_equal(lt_constr.f, _CSAF([1.2, 3.4], [_CVI(1), _CVI(3)], 0.0))
        @test lt_constr.s == _LT(7.8 - 5.6)
        @test isempty(node.gt_general_constrs)
        @test node.depth == 3
        @test node.dual_bound == -Inf
    end
    let f = _CSAF([2.4, 4.6], [_CVI(2), _CVI(1)], 6.8), s = _GT(8.0)
        bd = Cerberus.AffineConstraint(f, s)
        @inferred Cerberus.apply_branching!(node, bd)
        @test node.lt_bounds == [Cerberus.BoundUpdate(_CVI(2), _LT(1.0))]
        @test node.gt_bounds == [Cerberus.BoundUpdate(_CVI(4), _GT(3.0))]
        @test length(node.lt_general_constrs) == 1
        @test _is_equal(
            node.lt_general_constrs[1].f,
            _CSAF([1.2, 3.4], [_CVI(1), _CVI(3)], 0.0),
        )
        @test node.lt_general_constrs[1].s == _LT(7.8 - 5.6)
        @test length(node.gt_general_constrs) == 1
        gt_constr = node.gt_general_constrs[1]
        @test _is_equal(gt_constr.f, _CSAF([2.4, 4.6], [_CVI(2), _CVI(1)], 0.0))
        @test gt_constr.s == _GT(8.0 - 6.8)
        @test node.depth == 4
        @test node.dual_bound == -Inf
    end
end

struct DummyVariableBranchingRule <: Cerberus.AbstractVariableBranchingRule end

@testset "variable branching helpers" begin
    form = _build_dmip_formulation()
    config =
        Cerberus.AlgorithmConfig(branching_rule = DummyVariableBranchingRule())
    @testset "branching_candidates" begin
        let nr = Cerberus.NodeResult(
                Cerberus.OPTIMAL_LP,
                12.3,
                [0.6, 0.7, 0.1],
                12,
                13,
                14,
            )
            bc = @inferred Cerberus.branching_candidates(form, nr, config)
            @test bc == [
                Cerberus.VariableBranchingCandidate(_CVI(1), 0.6),
                Cerberus.VariableBranchingCandidate(_CVI(3), 0.1),
            ]
        end
        let nr = Cerberus.NodeResult(
                Cerberus.OPTIMAL_LP,
                12.3,
                [0.8, 0.7, 0.0],
                12,
                13,
                14,
            )
            bc = @inferred Cerberus.branching_candidates(form, nr, config)
            @test bc == [Cerberus.VariableBranchingCandidate(_CVI(1), 0.8)]
        end
        let nr = Cerberus.NodeResult(
                Cerberus.OPTIMAL_LP,
                12.3,
                [1.0, 0.7, 0.0],
                12,
                13,
                14,
            )
            bc = @inferred Cerberus.branching_candidates(form, nr, config)
            @test isempty(bc)
        end
    end

    @testset "branch_on" begin
        candidate = @inferred Cerberus.VariableBranchingCandidate(_CVI(1), 0.6)
        let node = Cerberus.Node()
            let score = @inferred Cerberus.VariableBranchingScore(0.7, 0.8, 0.9)
                nodes = @inferred Cerberus.branch_on(node, candidate, score)
                @test length(nodes) == 2
                let node = nodes[1]
                    @test node.lt_bounds ==
                          [Cerberus.BoundUpdate(_CVI(1), _LT(0.0))]
                    @test isempty(node.gt_bounds)
                    @test isempty(node.lt_general_constrs)
                    @test isempty(node.gt_general_constrs)
                    @test node.depth == 1
                end
                let node = nodes[2]
                    @test isempty(node.lt_bounds)
                    @test node.gt_bounds ==
                          [Cerberus.BoundUpdate(_CVI(1), _GT(1.0))]
                    @test isempty(node.lt_general_constrs)
                    @test isempty(node.gt_general_constrs)
                    @test node.depth == 1
                end
            end
            let score = @inferred Cerberus.VariableBranchingScore(0.8, 0.7, 0.9)
                nodes = @inferred Cerberus.branch_on(node, candidate, score)
                @test length(nodes) == 2

                let node = nodes[1]
                    @test isempty(node.lt_bounds)
                    @test node.gt_bounds ==
                          [Cerberus.BoundUpdate(_CVI(1), _GT(1.0))]
                    @test isempty(node.lt_general_constrs)
                    @test isempty(node.gt_general_constrs)
                    @test node.depth == 1
                end
                let node = nodes[2]
                    @test node.lt_bounds ==
                          [Cerberus.BoundUpdate(_CVI(1), _LT(0.0))]
                    @test isempty(node.gt_bounds)
                    @test isempty(node.lt_general_constrs)
                    @test isempty(node.gt_general_constrs)
                    @test node.depth == 1
                end
            end
        end
        let node = Cerberus.Node(
                Cerberus.BoundUpdate{_LT}[],
                [Cerberus.BoundUpdate(_CVI(3), _GT(1.0))],
                1,
            )
            let score = @inferred Cerberus.VariableBranchingScore(0.7, 0.8, 0.9)
                nodes = @inferred Cerberus.branch_on(node, candidate, score)
                @test length(nodes) == 2
                let node = nodes[1]
                    @test node.lt_bounds ==
                          [Cerberus.BoundUpdate(_CVI(1), _LT(0.0))]
                    @test node.gt_bounds ==
                          [Cerberus.BoundUpdate(_CVI(3), _GT(1.0))]
                    @test isempty(node.lt_general_constrs)
                    @test isempty(node.gt_general_constrs)
                    @test node.depth == 2
                end
                let node = nodes[2]
                    @test isempty(node.lt_bounds)
                    @test node.gt_bounds == [
                        Cerberus.BoundUpdate(_CVI(3), _GT(1.0)),
                        Cerberus.BoundUpdate(_CVI(1), _GT(1.0)),
                    ]
                    @test isempty(node.lt_general_constrs)
                    @test isempty(node.gt_general_constrs)
                    @test node.depth == 2
                end
            end
            let score = @inferred Cerberus.VariableBranchingScore(0.8, 0.7, 0.9)
                nodes = @inferred Cerberus.branch_on(node, candidate, score)
                @test length(nodes) == 2
                let node = nodes[1]
                    @test isempty(node.lt_bounds)
                    @test node.gt_bounds == [
                        Cerberus.BoundUpdate(_CVI(3), _GT(1.0)),
                        Cerberus.BoundUpdate(_CVI(1), _GT(1.0)),
                    ]
                    @test isempty(node.lt_general_constrs)
                    @test isempty(node.gt_general_constrs)
                    @test node.depth == 2
                end
                let node = nodes[2]
                    @test node.lt_bounds ==
                          [Cerberus.BoundUpdate(_CVI(1), _LT(0.0))]
                    @test node.gt_bounds ==
                          [Cerberus.BoundUpdate(_CVI(3), _GT(1.0))]
                    @test isempty(node.lt_general_constrs)
                    @test isempty(node.gt_general_constrs)
                    @test node.depth == 2
                end
            end
        end
    end
end

@testset "MostInfeasible" begin
    mi_config = Cerberus.AlgorithmConfig(
        branching_rule = Cerberus.MostInfeasible(),
        silent = true,
    )
    @testset "branching_score" begin
        state = Cerberus.CurrentState()
        let x = [0.6, 0.5, 0.35]
            nr = Cerberus.NodeResult(Cerberus.OPTIMAL_LP, 12.3, x, 12, 13, 14)
            let bc = Cerberus.VariableBranchingCandidate(_CVI(1), 0.6)
                vbs =
                    @inferred Cerberus.branching_score(state, bc, nr, mi_config)
                @test vbs == Cerberus.VariableBranchingScore(0.6, 0.4, 0.4)
            end
            let bc = Cerberus.VariableBranchingCandidate(_CVI(3), 0.35)
                vbs =
                    @inferred Cerberus.branching_score(state, bc, nr, mi_config)
                @test vbs == Cerberus.VariableBranchingScore(0.35, 0.65, 0.35)
            end
        end
    end

    @testset "end to end" begin
        let fm = _build_dmip_formulation()
            state = Cerberus.CurrentState()
            node = Cerberus.Node()
            x = [0.6, 0.7, 0.1]
            cost = 1.2
            result =
                Cerberus.NodeResult(Cerberus.OPTIMAL_LP, cost, x, 12, 13, 14)
            n1, n2 =
                @inferred Cerberus.branch(state, fm, node, result, mi_config)
            @test isempty(n1.lt_bounds)
            @test n1.gt_bounds == [Cerberus.BoundUpdate(_CVI(1), _GT(1.0))]
            @test n1.depth == 1
            @test n1.dual_bound == -Inf

            @test n2.lt_bounds == [Cerberus.BoundUpdate(_CVI(1), _LT(0.0))]
            @test isempty(n2.gt_bounds)
            @test n2.depth == 1
            @test n2.dual_bound == -Inf

            x2 = [1.0, 0.7, 0.1]
            result.x = x2
            n3, n4 = @inferred Cerberus.branch(state, fm, n2, result, mi_config)
            @test n3.lt_bounds == [
                Cerberus.BoundUpdate(_CVI(1), _LT(0.0)),
                Cerberus.BoundUpdate(_CVI(3), _LT(0.0)),
            ]
            @test isempty(n3.gt_bounds)
            @test n3.depth == 2
            @test n3.dual_bound == -Inf

            @test n4.lt_bounds == [Cerberus.BoundUpdate(_CVI(1), _LT(0.0))]
            @test n4.gt_bounds == [Cerberus.BoundUpdate(_CVI(3), _GT(1.0))]
            @test n4.depth == 2
            @test n4.dual_bound == -Inf

            # Nothing to branch on, should throw. Really, should have pruned by integrality before.
            x3 = [1.0, 0.7, 0.0]
            result.x = x3
            @test_throws ErrorException Cerberus.branch(
                state,
                fm,
                n4,
                result,
                mi_config,
            )
        end

        # General integer branching
        let fm = _build_gi_dmip_formulation()
            state = Cerberus.CurrentState()
            node = Cerberus.Node()
            x = [0.6, 0.4, 0.7]
            cost = 1.2
            result =
                Cerberus.NodeResult(Cerberus.OPTIMAL_LP, cost, x, 12, 13, 14)
            fc, oc =
                @inferred Cerberus.branch(state, fm, node, result, mi_config)
            @test fc.lt_bounds == [Cerberus.BoundUpdate(_CVI(2), _LT(0.0))]
            @test isempty(fc.gt_bounds)
            @test fc.depth == 1
            @test isempty(oc.lt_bounds)
            @test oc.gt_bounds == [Cerberus.BoundUpdate(_CVI(2), _GT(1.0))]
            @test oc.depth == 1

            x2 = [0.6, 0.4, 2.55]
            result.x = x2
            fc_2, oc_2 =
                @inferred Cerberus.branch(state, fm, oc, result, mi_config)
            @test isempty(fc_2.lt_bounds)
            @test fc_2.gt_bounds == [
                Cerberus.BoundUpdate(_CVI(2), _GT(1.0)),
                Cerberus.BoundUpdate(_CVI(3), _GT(3.0)),
            ]
            @test fc_2.depth == 2
            @test oc_2.lt_bounds == [Cerberus.BoundUpdate(_CVI(3), _LT(2.0))]
            @test oc_2.gt_bounds == [Cerberus.BoundUpdate(_CVI(2), _GT(1.0))]
            @test oc_2.depth == 2
        end
    end
end

@testset "StrongBranching" begin
    @testset "branching_score" begin
        let sb_config = Cerberus.AlgorithmConfig(
                branching_rule = Cerberus.StrongBranching(),
                silent = true,
            )
            v = [_SV(_VI(i)) for i in 1:3]
            ac1 = Cerberus.AffineConstraint(
                4.0 * v[1] + 2.0 * v[2] + 2.0 * v[3],
                _GT(5.0),
            )
            ac2 = Cerberus.AffineConstraint(
                4.0 * v[1] + 2.0 * v[2] + 5.0 * v[3],
                _LT(5.0),
            )
            aff_constrs = [ac1, ac2]
            bounds = [_IN(0.0, Inf) for i in 1:3]
            p = Cerberus.Polyhedron(aff_constrs, bounds)
            variable_kind = [_GI() for i in 1:3]
            obj = _CSAF([-1.0, 1.0, 0.0], [_CVI(1), _CVI(2), _CVI(3)], 0.0)
            form = Cerberus.DMIPFormulation(p, variable_kind, obj)

            state = Cerberus.CurrentState()
            node = Cerberus.pop_node!(state.tree)
            nr = Cerberus.process_node!(state, form, node, sb_config)

            bc = Cerberus.VariableBranchingCandidate(_CVI(1), nr.x[1])
            vbs = @inferred Cerberus.branching_score(state, bc, nr, sb_config)
            @test vbs == Cerberus.VariableBranchingScore(0.75, Inf, Inf)
        end

        let sb_config = Cerberus.AlgorithmConfig(
                branching_rule = Cerberus.StrongBranching(),
                silent = true,
            )
            v = [_SV(_VI(i)) for i in 1:3]
            ac1 = Cerberus.AffineConstraint(
                1.0 * v[1] + 1 / 9 * v[2] + 5 / 3 * v[3],
                _LT(2.0),
            )
            ac2 = Cerberus.AffineConstraint(
                1 / 9 * v[1] + 1.0 * v[2] + 5 / 3 * v[3],
                _LT(2.0),
            )
            ac3 = Cerberus.AffineConstraint(
                1.0 * v[1] + 1.0 * v[2] + 1.0 * v[3],
                _LT(2.0),
            )
            aff_constrs = [ac1, ac2]
            bounds = [_IN(0.0, Inf) for i in 1:3]
            p = Cerberus.Polyhedron(aff_constrs, bounds)
            variable_kind = [_ZO(), _ZO(), _GI()]
            obj = _CSAF([-1.0, -1.0, -1.5], [_CVI(1), _CVI(2), _CVI(3)], 0.0)
            form = Cerberus.DMIPFormulation(p, variable_kind, obj)

            state = Cerberus.CurrentState()
            node = Cerberus.pop_node!(state.tree)
            nr = Cerberus.process_node!(state, form, node, sb_config)

            bc = Cerberus.VariableBranchingCandidate(_CVI(3), nr.x[3])
            vbs = @inferred Cerberus.branching_score(state, bc, nr, sb_config)
            @test isapprox(vbs.down_branch_score, 0.8)
            @test isapprox(vbs.up_branch_score, 0.7)
            @test isapprox(vbs.aggregate_score, 43 / 60)
        end
    end

    @testset "end-to-end" begin
        let sb_config = Cerberus.AlgorithmConfig(
                branching_rule = Cerberus.StrongBranching(),
                silent = true,
            )
            v = [_SV(_VI(i)) for i in 1:3]
            ac1 = Cerberus.AffineConstraint(
                1.0 * v[1] + 1.0 * v[2] + 2.0 * v[3],
                _GT(4.5),
            )
            ac2 = Cerberus.AffineConstraint(
                1.0 * v[1] + 2.0 * v[2] + 1.0 * v[3],
                _GT(4.5),
            )
            ac3 = Cerberus.AffineConstraint(
                2.0 * v[1] + 1.0 * v[2] + 1.0 * v[3],
                _GT(4.5),
            )
            aff_constrs = [ac1, ac2, ac3]
            bounds = [_IN(0.0, 4.0) for i in 1:3]
            p = Cerberus.Polyhedron(aff_constrs, bounds)
            variable_kind = [_GI() for i in 1:3]
            obj = _CSAF([1.0, 2.0, 3.0], [_CVI(1), _CVI(2), _CVI(3)], 0.0)
            fm = Cerberus.DMIPFormulation(p, variable_kind, obj)

            state = Cerberus.CurrentState()
            node = Cerberus.pop_node!(state.tree)
            nr = Cerberus.process_node!(state, fm, node, sb_config)

            n1, n2 = @inferred Cerberus.branch(state, fm, node, nr, sb_config)

            @test isempty(n1.gt_bounds)
            @test n1.lt_bounds == [Cerberus.BoundUpdate(_CVI(2), _LT(0.0))]
            @test n1.depth == 1
            @test n1.dual_bound == -Inf

            @test n2.gt_bounds == [Cerberus.BoundUpdate(_CVI(2), _GT(1.0))]
            @test isempty(n2.lt_bounds)
            @test n2.depth == 1
            @test n2.dual_bound == -Inf

            nr1 = Cerberus.process_node!(state, fm, n1, sb_config)
            n3, n4 = @inferred Cerberus.branch(state, fm, n1, nr1, sb_config)
            @test n3.lt_bounds == [Cerberus.BoundUpdate(_CVI(2), _LT(0.0))]
            @test n3.gt_bounds == [Cerberus.BoundUpdate(_CVI(3), _GT(1.0))]
            @test n3.depth == 2
            @test n3.dual_bound == -Inf

            @test n4.lt_bounds == [
                Cerberus.BoundUpdate(_CVI(2), _LT(0.0)),
                Cerberus.BoundUpdate(_CVI(3), _LT(0.0)),
            ]
            @test isempty(n4.gt_bounds)
            @test n4.depth == 2
            @test n4.dual_bound == -Inf
        end
    end
end

struct DummyBranchingRule <: Cerberus.AbstractBranchingRule end
struct DummyBranchingCandidate <: Cerberus.AbstractBranchingCandidate
    cvi::_CVI
end
struct DummyBranchingScore <: Cerberus.AbstractBranchingScore
    val::Any
end
function Cerberus.branching_candidates(
    ::Cerberus.DMIPFormulation,
    ::Cerberus.NodeResult,
    config::Cerberus.AlgorithmConfig{DummyBranchingRule},
)
    @assert config.branching_rule === DummyBranchingRule()
    return [DummyBranchingCandidate(_CVI(1)), DummyBranchingCandidate(_CVI(3))]
end
function Cerberus.branching_score(
    ::Cerberus.CurrentState,
    dbc::DummyBranchingCandidate,
    ::Cerberus.NodeResult,
    config::Cerberus.AlgorithmConfig{DummyBranchingRule},
)
    @assert config.branching_rule === DummyBranchingRule()
    return DummyBranchingScore(Cerberus.index(dbc.cvi))
end
Cerberus.aggregate_score(::DummyBranchingScore) = 42.0
function Cerberus.branch_on(
    node::Cerberus.Node,
    ::DummyBranchingCandidate,
    ::DummyBranchingScore,
)
    return (
        Cerberus.Node(
            vcat(node.lt_bounds, Cerberus.BoundUpdate(_CVI(3), _LT(1.0))),
            vcat(node.gt_bounds, Cerberus.BoundUpdate(_CVI(3), _GT(1.0))),
            node.depth + 1,
        ),
        Cerberus.Node(
            vcat(node.lt_bounds, Cerberus.BoundUpdate(_CVI(3), _LT(2.0))),
            vcat(node.gt_bounds, Cerberus.BoundUpdate(_CVI(3), _GT(2.0))),
            node.depth + 1,
        ),
        Cerberus.Node(
            vcat(node.lt_bounds, Cerberus.BoundUpdate(_CVI(3), _LT(3.0))),
            vcat(node.gt_bounds, Cerberus.BoundUpdate(_CVI(3), _GT(3.0))),
            node.depth + 1,
        ),
    )
end

@testset "branch" begin
    form = _build_dmip_formulation()
    config = Cerberus.AlgorithmConfig(
        branching_rule = DummyBranchingRule(),
        silent = true,
    )
    state = Cerberus.CurrentState()
    node = Cerberus.Node(
        Cerberus.BoundUpdate{_LT}[],
        [Cerberus.BoundUpdate(_CVI(1), _GT(1.0))],
        1,
    )
    nr = Cerberus.NodeResult(
        Cerberus.OPTIMAL_LP,
        12.3,
        [0.6, 0.5, 1.0],
        12,
        13,
        14,
    )
    nodes = @inferred Cerberus.branch(state, form, node, nr, config)
    @test length(nodes) == 3
    for i in 1:3
        node = nodes[i]
        @test node.lt_bounds == [Cerberus.BoundUpdate(_CVI(3), _LT(i))]
        @test node.gt_bounds == [
            Cerberus.BoundUpdate(_CVI(1), _GT(1.0)),
            Cerberus.BoundUpdate(_CVI(3), _GT(i)),
        ]
        @test isempty(node.lt_general_constrs)
        @test isempty(node.gt_general_constrs)
        @test node.depth == 2
    end
end


text:
abstract type AbstractObjective end

# Given callables to calculate objectives and partial first derivatives
# create a function that calculates both.
function make_fdf(x, F, f!, j!)
    function fj!(fx, jx, x)
        j!(jx, x)
        return f!(fx, x)
    end
end
function make_fdf(x, F::Number, f, g!)
    function fg!(gx, x)
        g!(gx, x)
        return f(x)
    end
end

# Initialize an n-by-n Jacobian
alloc_DF(x, F) = fill(eltype(x)(NaN), length(F), length(x))
# Initialize a gradient shaped like x
alloc_DF(x, F::Number) = x_of_nans(x)
# Initialize an n-by-n Hessian
alloc_H(x) = alloc_DF(x, x)

ext: jl


text:
export Partition, Partitions, YoungTableau, SkewDiagram
export dim, has_left_neighbor, has_bottom_neighbor, inskewdiag, isrimhook,
       hooklength, leglength, matrix_repr, partitionseq

##############################################################################
#
#   Partition type, AbstractVector interface
#
##############################################################################

length(p::Partition) = length(p.part)

size(p::Partition) = size(p.part)

Base.IndexStyle(::Type{T}) where {T <: Partition} = Base.IndexLinear()

getindex(p::Partition, i::Int) = p.part[i]

function setindex!(p::Partition, v::Int, i::Int)
   prev = Inf
   nex = 1
   if i == length(p)
      prev = p[i-1]
   elseif i == 1
      nex = p[2]
   elseif 1 < i < length(p)
      prev = p[i-1]
      nex = p[i+1]
   end
   nex <= v <= prev || throw("Partition must be positive and non-increasing")
   p.n += v - p.part[i]
   p.part[i] = v
   return p
end

==(p::Partition, m::Partition) = p.n ==m.n && p.part == m.part
hash(p::Partition, h::UInt) = hash(p.part, hash(Partition, h))

convert(::Type{Partition}, p::Vector{Int}) = Partition(p)

##############################################################################
#
#   IO for Partition
#
##############################################################################

function subscriptify(n::Int)
   subscript_0 = Int(0x2080) # Char(0x2080) -> subscript 0
   return join([Char(subscript_0 + i) for i in reverse(digits(n))])
end

function show(io::IO, p::Partition)
   uniq = unique(p.part)
   mults = [count(i -> i == u, p.part) for u in uniq]
   str = join((string(u)*subscriptify(m) for (u,m) in zip(uniq, mults)))
   print(io, str)
end

##############################################################################
#
#   Iterator interface for Integer Partitions
#
##############################################################################

# Implemented following RuleAsc (Algorithm 3.1) from
#    "Generating All Partitions: A Comparison Of Two Encodings"
# by Jerome Kelleher and Barry O’Sullivan, ArXiv:0909.2331

function Base.start(parts::Partitions)
    if parts.n < 1
        return (Int[], 0)
    elseif parts.n == 1
        return ([1], 0)
    else
        p = zeros(Int, parts.n)
        p[2] = parts.n
        return (p, 2)
    end
end

Base.next(parts::Partitions, state) = nextpart_asc(state...)
Base.done(parts::Partitions, state) = state[2] == 1
Base.eltype(::Type{Partitions}) = Partition
length(parts::Partitions) = BigInt(numpart(parts.n))

function nextpart_asc(part, k)
    if k == 0
        return Partition(part, false), (part, 1)
    end
    y = part[k] - 1
    k -= 1
    x = part[k] + 1
    while x <= y
        part[k] = x
        y -= x
        k += 1
    end
    part[k] = x + y
    return Partition(reverse(part[1:k]), false), (part, k)
end

doc"""
    conj(part::Partition)
> Returns the conjugated partition of `part`, i.e. the partition corresponding
> to the Young tableau of `part` reflected through the main diagonal.
"""
function conj(part::Partition)
    p = Int[]
    for i in 1:part.n
        n = sum(part .>= i)
        n == 0 && break
        push!(p, n)
    end
    return Partition(p, false)
end

##############################################################################
#
#   Partition sequences and Murnaghan-Nakayama formula
#
##############################################################################

doc"""
    partitionseq(lambda::Partition)
> Returns a sequence (as `BitVector`) of `false`s and `true`s constructed from
> `lambda`: tracing the lower contour of the Young Diagram associated to
> `lambda` from left to right a `true` is inserted for every horizontal and
> `false` for every vertical step. The sequence always starts with `true` and
> ends with `false`.
"""
function partitionseq(lambda::Partition)
   seq = trues(maximum(lambda) + length(lambda))
   j = lambda[end]
   for i in (length(lambda)-1):-1:1
      seq[j+1] = false
      j += lambda[i] - lambda[i+1] + 1
   end
   seq[j+1] = false
   return seq
end

partitionseq(v::Vector{Int}) = partitionseq(Partition(v))

doc"""
    partitionseq(seq::BitVector)
> Returns the essential part of the sequence `seq`, i.e. a subsequence starting
> at first `true` and ending at last `false`.
"""
partitionseq(seq::BitVector) = seq[findfirst(seq, true):findlast(seq, false)]

doc"""
    isrimhook(R::BitVector, idx::Int, len::Int)
> `R[idx:idx+len]` forms a rim hook in the Young Diagram of parition
> corresponding to `R` iff `R[idx] == true` and `R[idx+len] == false`.
"""
function isrimhook(R::BitVector, idx::Int, len::Int)
   return (R[idx+len] == false) && (R[idx] == true)
end

const _charvalsTable = Dict{Tuple{BitVector,Vector{Int}}, BigInt}()

doc"""
    MN1inner(R::BitVector, mu::Partition, t::Int, [charvals])
> Returns the value of $\lambda$-th irreducible character on conjugacy class of
> permutations represented by partition `mu`, where `R` is the (binary)
> partition sequence representing $\lambda$. Values already computed are stored
> in `charvals::Dict{Tuple{BitVector,Vector{Int}}, Int}`.
> This is an implementation (with slight modifications) of the
> Murnaghan-Nakayama formula as described in
>
>     Dan Bernstein,
>     "The computational complexity of rules for the character table of Sn"
>     _Journal of Symbolic Computation_, 37(6), 2004, p. 727-748.
"""
function MN1inner(R::BitVector, mu::Partition, t::Int,
        charvals=Dict{Tuple{BitVector,Vector{Int}}, BigInt}())
    if t > length(mu)
        chi = 1
    elseif mu[t] > length(R)
        chi = 0
    else
        chi = 0
        sgn = false

        for j in 1:mu[t]-1
            if R[j] == false
                sgn = !sgn
            end
        end
        for i in 1:length(R)-mu[t]
            if R[i] != R[i+mu[t]-1]
                sgn = !sgn
            end
            if isrimhook(R, i, mu[t])
                R[i], R[i+mu[t]] = R[i+mu[t]], R[i]
                essR = (partitionseq(R), mu[t+1:end])
                if !haskey(charvals, essR)
                    charvals[essR] = MN1inner(R, mu, t+1, charvals)
                end
                chi += (-1)^Int(sgn)*charvals[essR]
                R[i], R[i+mu[t]] = R[i+mu[t]], R[i]
            end
        end
    end
    return chi
end

##############################################################################
#
#   YoungTableau type, AbstractVector interface
#
##############################################################################

function YoungTableau(part::Partition, fill::Vector{Int}=collect(1:part.n))
   part.n == length(fill) || throw("Can't fill Young digaram of $part with $fill: different number of elemnets.")
   tab = zeros(Int, length(part), maximum(part))
   k=1
   for (idx, p) in enumerate(part)
      tab[idx, 1:p] = fill[k:k+p-1]
      k += p
   end
   return YoungTableau(part, tab)
end

YoungTableau(p::Vector{Int}) = YoungTableau(Partition(p))

size(Y::YoungTableau) = size(Y.tab)

Base.IndexStyle(::Type{T}) where {T <: YoungTableau} = Base.IndexLinear()

getindex(Y::YoungTableau, i::Int) = Y.tab[i]

function ==(Y1::YoungTableau,Y2::YoungTableau)
   Y1.part == Y2.part || return false
   Y1.tab == Y2.tab || return false
   return true
end

hash(Y::YoungTableau, h::UInt) = hash(Y.part, hash(Y.tab, hash(YoungTableau, h)))

doc"""
    conj(Y::YoungTableau)
> Returns the conjugated tableau, i.e. the tableau reflected through the main
> diagonal.
"""
conj(Y::YoungTableau) = YoungTableau(conj(Y.part), transpose(Y.tab))

##############################################################################
#
#   Misc functions for YoungTableaux
#
##############################################################################

rowlen(Y::YoungTableau, i, j) = sum(Y[i, j:end] .> 0)
collen(Y::YoungTableau, i, j) = sum(Y[i:end, j] .> 0)

doc"""
    hooklength(Y::YoungTableau, i, j)
> Returns the hooklength of an element in `Y` at position `(i,j)`. `hooklength`
> will return `0` for `(i,j)` not in the tableau `Y`.
"""
function hooklength(Y::YoungTableau, i, j)
   if Y[i,j] == 0
      return 0
   else
      return rowlen(Y, i, j) + collen(Y, i, j) - 1
   end
end

doc"""
    dim(Y::YoungTableau)
> Returns the dimension of the irreducible representation of
> `PermutationGroup(n)` associated to YoungTableau `Y` of a partition of `n`.
"""
function dim(Y::YoungTableau)
   n, m = size(Y)
   num = fac(maximum(Y))
   den = reduce(*, 1, BigInt(hooklength(Y,i,j)) for i in 1:n, j in 1:m if j <= Y.part[i])
   return divexact(num, den)
end

##############################################################################
#
#   SkewDiagrams
#
##############################################################################

SkewDiagram(lambda::Vector{Int}, mu::Vector{Int}) = SkewDiagram(Partition(lambda), Partition(mu))

/(lambda::Partition, mu::Partition) = SkewDiagram(lambda, mu)

==(xi::SkewDiagram, psi::SkewDiagram) = xi.lam == psi.lam && xi.mu == psi.mu
hash(xi::SkewDiagram, h::UInt) = hash(xi.lam, hash(xi.mu, hash(SkewDiagram, h)))

###############################################################################
#
#   String I/O
#
###############################################################################

doc"""
    matrix_repr(xi::SkewDiagram)
> Returns a binary representation of the diagram `xi`, i.e. a binary array `A`
> where `A[i,j] == 1` if and only if `(i,j)` is in `xi.lam` but not in `xi.mu`.
"""
function matrix_repr(xi::SkewDiagram)
   ydiag = zeros(Int,length(xi.lam), maximum(xi.lam))
   for i in 1:length(xi.mu)
      ydiag[i, xi.mu[i]+1:xi.lam[i]] .= 1
   end
   for i in length(xi.mu)+1:length(xi.lam)
      ydiag[i,1:xi.lam[i]] .= 1
   end
   return ydiag
end

show(io::IO, xi::SkewDiagram) = show(io, MIME("text/plain"), matrix_repr(xi))

##############################################################################
#
#   Misc functions for SkewDiagrams
#
##############################################################################

doc"""
    inskewdiag(xi::SkewDiagram, i::Int, j::Int)
> Checks if box at position `(i,j)` belongs to the skew diagram `xi`.
"""
function inskewdiag(xi::SkewDiagram, i::Int, j::Int)
   if i <= 0 || j <= 0
      return false
   elseif i > length(xi.lam) || j > maximum(xi.lam)
      return false
   elseif length(xi.mu) >= i
      return xi.mu[i] < j ≤ xi.lam[i]
   else
      return j ≤ xi.lam[i]
   end
end

doc"""
    has_left_neighbor(xi::SkewDiagram, i::Int, j::Int)
> Checks if box at position `(i,j)` has neighbour in `xi` to the left.
"""
function has_left_neighbor(xi::SkewDiagram, i::Int, j::Int)
   if j == 1
      return false
   else
      return inskewdiag(xi, i, j) && inskewdiag(xi, i, j-1)
   end
end

doc"""
    has_bottom_neighbor(xi::SkewDiagram, i::Int, j::Int)
> Checks if box at position `(i,j)` has neighbour in `xi` below.
"""
function has_bottom_neighbor(xi::SkewDiagram, i::Int, j::Int)
   if i == length(xi.lam)
      return false
   else
      return inskewdiag(xi, i, j) && inskewdiag(xi, i+1, j)
   end
end

doc"""
    isrimhook(xi::SkewDiagram)
> Checks if `xi` represents a rim-hook diagram, i.e. its diagram is
> edge-connected and contains no $2\times 2$ squares.
"""
function isrimhook(xi::SkewDiagram)
   i = 1
   j = xi.lam[1]
   while i != length(xi.lam) && j != 1
      left = has_left_neighbor(xi, i,j)
      bottom = has_bottom_neighbor(xi, i,j)
      if left && bottom # there is 2×2 square in xi
         return false
      elseif left
         j -= 1
      elseif bottom
         i += 1
      else
         lam_tail = xi.lam[i+1:end]
         mu_tail = zeros(Int, length(lam_tail))
         mu_tail[1:length(xi.mu)-i] = xi.mu[i+1:end]

         if any(lam_tail .- mu_tail .> 0)
            return false # xi is disconnected
         else
            return true # we arrived at the end of xi
         end
      end
   end
   return true
end

doc"""
    leglength(xi::SkewDiagram, check::Bool=true)
> Computes the leglength of a rim-hook `xi`, i.e. the number of rows with
> non-zero entries minus one. If `check` is `false` function will not check
> whether `xi` is actually a rim-hook.
"""
function leglength(xi::SkewDiagram, check::Bool=true)
   if check
      isrimhook(xi) || throw("$xi is not a rimhook. leglength is defined only for rim hooks")
   end
   m = zeros(length(xi.lam))
   m[1:length(xi.mu)] = xi.mu
   return sum((xi.lam .- m) .> 0) - 1
end


text:
const packages = [
    "AverageShiftedHistograms",
    "BDF",
    "BrainWave",
    "Diversity",
    "IterativeSolvers",
    "RCall",
    "SDT",
    "Sims",
    "TargetedLearning"
    ]

const failures = Set()

for package in packages
    print(" - ", package)
    try
        require(package)
        println(" ✓")
    catch err
        push!(failures, (package, err))
        println(" ×")
    end
end

println("-"^50)

for (package, err) in failures
    println("""

    # $(package)

    $(err)

    ---
    """)
end

failures


text:
# function Modelling example

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
# In this example, a simple homogeneous elastic modelling is shown with proper#
# settings. The example is a multisource example.                             #
# To run the example, copy the following code into Julia interface. You may   #
# need to check if the path has been added
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#

    using seisimu

    pvel = 3500
    svel = nothing
    rho = 2.8
    DZ = 1000
    HX = 1000
    pkf = 20
    T = 0.5
    ext = 10
    iflag = 2
    dz = nothing
    dx = nothing
    dt= nothing
    nT= nothing
    mor = false

    #=== elastic ===#
    # model = initmodel(pvel,svel,rho,DZ,HX,pkf,T);

    #=== acoustic ===#
    model = initmodel(pvel,rho,DZ,HX,pkf,T)

    sn = 2
    loc = [500 500; 800 300]
    ot = [0.05, 0.08]
    tp = ["expl", "expl"]
    wavelet = Ricker(pkf,model.medium.dt)

    waveform = [wavelet wavelet]

    sou = initsource(sn,loc,ot,tp,waveform,model)

    @time extrap2d(model,sou)


text:
module MixedAnova
using DataFrames, Statistics, StatsBase, LinearAlgebra, Distributions, Reexport, Printf
@reexport using StatsModels
import StatsBase: fit!, fit, dof, dof_residual, deviance, nobs, coefnames
import StatsModels: TableStatisticalModel, TableRegressionModel, vectorize, kron_insideout,
                    ModelFrame, ModelMatrix, response, columntable, asgn, collect_matrix_terms
import LinearAlgebra.BlasReal
import Tables.istable
import Base: show

export
    # wrappers
    AnovaResult,

    # anova functions
    anova, anova_lm, anova_lme, anova_glm, anova_lfe,

    # GoodnessOfFit
    GoodnessOfFit, FTest, LikelihoodRatioTest, LRT, canonicalgoodnessoffit,

    # Others
    lme, glme, lfe, nestedmodels, calcdof, formula, #getterms,
    teststat, pval, anova_test, anova_type, to_trm,

    # init
    glm_init, mm_init, fem_init

# Test
"""
    abstract type GoodnessOfFit end
"""
abstract type GoodnessOfFit end
"""
    struct FTest <: GoodnessOfFit end

ANOVA by F-test. It can be the first argument or keyword argument `test`.
"""
struct FTest <: GoodnessOfFit end
"""
    struct LikelihoodRatioTest <: GoodnessOfFit end
    const LRT = LikelihoodRatioTest

ANOVA by likelihood-ratio test. It can be the first argument or keyword argument `test`.
"""
struct LikelihoodRatioTest <: GoodnessOfFit end
"""
    const LRT = LikelihoodRatioTest

See [`LikelihoodRatioTest`](@ref).
"""
const LRT = LikelihoodRatioTest

# Wrapper for ANOVA
"""
    AnovaResult{M, T, N}

Returned object of `anova`.
* `M` is a subtype of `Tuple` if multiple models are provided; otherwise, a typeof model.
* `T` is a subtype of `GoodnessOfFit`; either `FTest` or `LRT`.
* `N` is the length of parameters.

## Fields

* `model`: full model or tuple of tested models.
* `type`: type of `anova`.
* `dof`: degree of freedoms of models or factors.
* `deviance`: deviance(s) for calculating test statistics. See `deviance` for more details.
* `teststat`: value(s) of test statiscics.
* `pval`: p-value(s) of test statiscics.
* `tests`: `NamedTuple` contained extra statistics.
"""
struct AnovaResult{M, T, N}
    model::M
    type::Int
    dof::NTuple{N, Int}
    deviance::NTuple{N, Float64}
    teststat::NTuple{N, Float64}
    pval::NTuple{N, Float64}
    tests::NamedTuple
end

AnovaResult{T}(model::M,
                type::Int,
                dof::NTuple{N, Int},
                deviance::NTuple{N, Float64},
                teststat::NTuple{N, Float64},
                pval::NTuple{N, Float64},
                tests::NamedTuple) where {M, N, T <: GoodnessOfFit} =
    AnovaResult{M, T, N}(model, type, dof, deviance, teststat, pval, tests)

include("api.jl")
include("termIO.jl")
const path = @__DIR__()

function glm_init()
    include(joinpath(path, "GLM", "anova.jl"))
    include(joinpath(path, "GLM", "io.jl"))
    include(joinpath(path, "GLM", "fit.jl"))
    return
end

function mm_init()
    include(joinpath(path, "MixedModels", "anova.jl"))
    include(joinpath(path, "MixedModels", "io.jl"))
    include(joinpath(path, "MixedModels", "fit.jl"))
    return
end

function fem_init()
    include(joinpath(path, "FixedEffectModels", "anova.jl"))
    include(joinpath(path, "FixedEffectModels", "io.jl"))
    include(joinpath(path, "FixedEffectModels", "fit.jl"))
    return
end

end


text:
# Tests of intstats.jl

using MLBase
using Test

# classify

ss = rand(8, 50)
for i = 1:size(ss,2)
    ss_i = ss[:,i]
    kmax = argmax(ss_i)
    kmin = argmin(ss_i)
    vmax = ss_i[kmax]
    vmin = ss_i[kmin]

    @test classify(ss_i) == kmax
    @test classify(ss_i, Forward) == kmax
    @test classify(ss_i, Reverse) == kmin

    @test classify(ss_i, 0.8) == (vmax >= 0.8 ? kmax : 0)
    @test classify(ss_i, 0.8, Forward) == (vmax >= 0.8 ? kmax : 0)
    @test classify(ss_i, 0.2, Reverse) == (vmin <= 0.2 ? kmin : 0)

    @test classify_withscore(ss_i) == (kmax, ss_i[kmax])
    @test classify_withscore(ss_i, Forward) == (kmax, ss_i[kmax])
    @test classify_withscore(ss_i, Reverse) == (kmin, ss_i[kmin])
end

rmax = Int[argmax(ss[:,i]) for i = 1:size(ss,2)]
rmin = Int[argmin(ss[:,i]) for i = 1:size(ss,2)]
vmax = ss[LinearIndices(size(ss))[CartesianIndex.(rmax, 1:size(ss,2))]]
vmin = ss[LinearIndices(size(ss))[CartesianIndex.(rmin, 1:size(ss,2))]]

trmax = copy(rmax); trmax[vmax .< 0.8] .= 0
trmin = copy(rmin); trmin[vmin .> 0.2] .= 0

@test classify(ss) == rmax
@test classify(ss, Forward) == rmax
@test classify(ss, Reverse) == rmin

@test classify_withscores(ss) == (rmax, vmax)
@test classify_withscores(ss, Forward) == (rmax, vmax)
@test classify_withscores(ss, Reverse) == (rmin, vmin)

@test classify(ss, 0.8) == trmax
@test classify(ss, 0.8, Forward) == trmax
@test classify(ss, 0.2, Reverse) == trmin

# labelmap & labelencode

xs = ["a", "a", "b", "b", "a", "b", "c", "a"]
labels = [1, 1, 2, 2, 1, 2, 3, 1]
lmap = labelmap(xs)

@test keys(lmap) == ["a", "b", "c"]
@test labelencode(lmap, xs) == labels
@test labeldecode(lmap, labels) == xs

# groupindices

gs = Any[[1,2,5,8],[3,4,6],[7]]

@test groupindices(3, labels) == gs
@test groupindices(lmap, xs) == gs


text:
using OctreeBH
using StaticArrays

const BOXSIZE = 1.0

const N = 3 #spatial dimension
const T = Float64

Npart = 100000 #number of particles

#desired number of neighbor particles (ngbs) within a search radius (hsml)
Nngb0 = 32

#searching radius that should contain roughly Nngb0 neighboring particles
hsml0 = BOXSIZE * (Nngb0/(4*pi/3*Npart))^(1/N)

@testset "ngbs" begin

#randomly distributing particles
X = [@SVector rand(N) for _ in 1:Npart]

topnode_length = @SVector(ones(N)) * BOXSIZE  #actual length of tree

center = topnode_length .* 0.5

#individual smoothing length (only relevant for scatter ngbs)
δ = 0.0 #add a small perturbation
hsml = ones(Npart) .* hsml0 .* ( 1 .+ δ .* (rand(Npart) .- 0.5) )

mass = ones(Npart) #particle mass

part = [Data{N,T}(SVector(X[i]), i, hsml[i], mass[i]) for i in eachindex(X)]

#build the tree
tree = buildtree(part, center, topnode_length);

#search center
x = @SVector rand(N)

boxsizes = @SVector(ones(N)) * BOXSIZE  #for periodic B.C.

#get indices of the gather neighbors
idx_ngbs_g = get_gather_ngb_tree(x, hsml0, tree, boxsizes)

#get indices of the scatter neighbors
idx_ngbs_s = get_scatter_ngb_tree(x, tree, boxsizes)

#for constant hsml, the two ngbs are identical
@test all(sort(idx_ngbs_g) .== sort(idx_ngbs_s))

end


text:
export setindex

struct SetIndex{T,N} <: ArrayOp{T,N}
    input::ArrayOp{T,N}
    idx::Tuple
    val
end

function setindex(x::ArrayOp, val, idxs...)
    SetIndex(x, idxs, val)
end

function stage(ctx, sidx::SetIndex)
    inp = cached_stage(ctx, sidx.input)

    dmn = domain(inp)
    idxs = [if isa(sidx.idx[i], Colon)
        indexes(dmn)[i]
    else
        sidx.idx[i]
    end for i in 1:length(sidx.idx)]

    ps = Array{Any}(size(chunks(inp)))
    ps[:] = chunks(inp)
    subdmns = domainchunks(inp)
    d = ArrayDomain(idxs)

    groups = map(group_indices, subdmns.cumlength, indexes(d))
    sz = map(length, groups)
    pieces = Array{Union{Chunk, Thunk}}(sz)
    for i = CartesianRange(sz)
        idx_and_dmn = map(getindex, groups, i.I)
        idx = map(x->x[1], idx_and_dmn)
        local_dmn = ArrayDomain(map(x->x[2], idx_and_dmn))
        s = subdmns[idx...]
        part_to_set = sidx.val
        ps[idx...] = Thunk(ps[idx...]) do p
            q = copy(p)
            q[indexes(project(s, local_dmn))...] = part_to_set
            q
        end
    end
    inp.chunks = ps
    inp
end


text:
module SketchyCGAL

using LinearAlgebra
using SparseArrays
# using Arpack
# using ArnoldiMethod
using Random
using Printf

include("utils.jl")
include("types.jl")
include("sketch.jl")
include("eig.jl")
include("cgal.jl")
include("scgal.jl")
include("scgal-maxcut.jl")

export scgal_full
export construct_Xhat

end


text:
"""
    JETBenchmarkUtils

This module exports utility macros to execute and benchmark arbitrary code in new Julia process:
- [`@freshexec ex`](@ref): executes `ex` in fresh process and restore its result back into the caller runtime
- [`@benchmark_freshexec ex`](@ref): returns statistics of execution of `ex` in a fresh process

!!! note
    `@benchmark_freshexec` is not integrated with [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl)'s CI infrastructure.
    `JETBenchmarkTools` provides utilties for that.
"""
module JETBenchmarkUtils

import Base.Meta: isexpr

"""
    ret = @freshexec [setup_ex] ex

Runs `ex` in an external process and gets back the final result (, which is supposed to be
  such a simple Julia object that we can restore it from its string representation).
Running in external process can be useful for testing JET analysis, because:
- the first time analysis is not affected by native code cache and JET's global report cache
- we can do something that might break a Julia process without wondering it breaks the
  original test process later on
The optional positional argument `setup_ex` runs before each execution of `ex` and defaults
  to `JET_LOAD_EX`, which just loads JET into the process.
"""
macro freshexec(args...)
    args = map(a->Expr(:quote,a), args)
    return Expr(:escape, Expr(:call, GlobalRef(@__MODULE__, :freshexec), __module__, args...))
end
freshexec(mod, ex) = freshexec(mod, JET_LOAD_EX, ex)
freshexec(args...) = _freshexec(args..., collect_last_result)

function collect_last_result(exs)
    lines = string.(exs)
    lines[end] = "ret = $(lines[end])"
    return join(lines, '\n')
end

"""
    stats = @benchmark_freshexec [ntimes = 5] [setup_ex] ex

Runs `ex` in an external process multiple times (which can be configured by the optional
  keyword argument `ntimes`), and collects execution statistics from [`@timed`](@ref).
  The final result is the statistics from a trial with minimum execution time.
  The statistics are generated by taking the `mean` of all the trials.
Each "evaluation" of `ex` is done in new Julia process in order to benchmark the performance
  of the "first-time analysis", where the native code cache and JET's global report cache
  have no effect for execution (of course, this is very time-consuming though ...).
The optional positional argument `setup_ex` runs before each execution of `ex` and its
  execution statistis are not included in the benchmark result; it defaults to
  `JET_WARMUP_EX`, which loads JET and runs a warm up analysis `@report_call identity(nothing)`.
"""
macro benchmark_freshexec(args...)
    isn(x) = isexpr(x, :(=)) && first(x.args) === :ntimes
    i = findfirst(isn, args)
    ntimes = i === nothing ? 5 : Expr(:quote, last(args[i].args))
    args = map(a->Expr(:quote,a), filter(!isn, args))
    return Expr(:escape, Expr(:call, GlobalRef(@__MODULE__, :benchmark_freshexec), ntimes, __module__, args...))
end
benchmark_freshexec(ntimes, mod, ex) = benchmark_freshexec(ntimes, mod, JET_WARMUP_EX, ex)
function benchmark_freshexec(ntimes, args...)
    stats = [_freshexec(args..., collect_statistics) for _ in 1:ntimes]
    return stats[argmin(getproperty.(stats, :time))]
end

function collect_statistics(exs)
    return """
    ret = @timed begin
        $(join(exs, '\n'))
        nothing # ensure `stats` can be parsed
    end
    """
end

function _freshexec(mod, setup_ex, ex, exs2script)
    # we need to flatten block expression into a toplevel expression to correctly handle
    # e.g. macro expansions
    setup_exs = isexpr(setup_ex, :block) ? setup_ex.args : [setup_ex]
    setup_script = join(string.(setup_exs), '\n')
    exs = isexpr(ex, :block) ? ex.args : [ex]
    script = exs2script(exs)

    prog = """
    old = stdout
    rw, wr = redirect_stdout()
    $(setup_script)
    $(script)
    redirect_stdout(old)
    close(rw); close(wr)
    println(stdout, repr(ret))
    """

    cmd = Cmd([JULIA_BIN, "-e", prog])
    io = IOBuffer()
    run(pipeline(cmd; stdout = io))

    return Core.eval(mod, Meta.parse(String(take!(io))))
end

const JET_LOAD_EX = :(using JET)
const JET_WARMUP_EX = quote
    using JET
    @report_call identity(nothing) # warm-up for JET analysis
end
const JULIA_BIN = normpath(Sys.BINDIR, "julia")

export
    @freshexec,
    @benchmark_freshexec

end


text:
using ThreeBodyDecay

let
    ev = LinRange(2.0, 10.0, 100)
    cal = iRhoQTB.(ev.^2, 3.0, 1.0, 0.2, 1.0)
    # plot(ev, [real.(cal) imag(cal)], lab="")
end

# iRhoQTB(5.0^2, 3, 1, 0.2, 2)

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
module BVH

using CthulhuVision.Light
using CthulhuVision.Random
using CthulhuVision.Materials
using CthulhuVision.Worlds
using CthulhuVision.Math
using CthulhuVision.Triangles

using StaticArrays
using CUDAnative

export Hitable, BVHNode, BVHAcceleration, AABB
export buildacceleration, hitacceleration

struct AABB
    min::Vector3
    max::Vector3

    function AABB(v1::Vector3, v2::Vector3, v3::Vector3)
        mini = Vector3(
            min(v1.x, v2.x, v3.x) - 0.01f0,
            min(v1.y, v2.y, v3.y) - 0.01f0,
            min(v1.z, v2.z, v3.z) - 0.01f0,
        )
        maxi = Vector3(
            max(v1.x, v2.x, v3.x) + 0.01f0,
            max(v1.y, v2.y, v3.y) + 0.01f0,
            max(v1.z, v2.z, v3.z) + 0.01f0,
        )

        new(mini, maxi)
    end

    function AABB(box1::AABB, box2::AABB)
        mini = Vector3(
            min(box1.min.x, box2.min.x),
            min(box1.min.y, box2.min.y),
            min(box1.min.z, box2.min.z),
        )
        maxi = Vector3(
            max(box1.max.x, box2.max.x),
            max(box1.max.y, box2.max.y),
            max(box1.max.z, box2.max.z),
        )
        new(mini, maxi)
    end

    function AABB()
        new(Vector3(0f0, 0f0, 0f0), Vector3(0f0, 0f0, 0f0))
    end
end

@inline function hitbox(aabb::AABB, ray::Ray, tmin::Float32, tmax::Float32) :: Bool
    inversedistancex = 1.0f0 / direction(ray).x
    tx0 = (aabb.min.x - origin(ray).x) * inversedistancex
    tx1 = (aabb.max.x - origin(ray).x) * inversedistancex
    if inversedistancex < 0.0f0
        tx0, tx1 = tx1, tx0
    end

    inversedistancey = 1.0f0 / direction(ray).y
    ty0 = (aabb.min.y - origin(ray).y) * inversedistancey
    ty1 = (aabb.max.y - origin(ray).y) * inversedistancey
    if inversedistancey < 0.0f0
        ty0, ty1 = ty1, ty0
    end

    inversedistancez = 1.0f0 / direction(ray).z
    tz0 = (aabb.min.z - origin(ray).z) * inversedistancez
    tz1 = (aabb.max.z - origin(ray).z) * inversedistancez
    if inversedistancez < 0.0f0
        tz0, tz1 = tz1, tz0
    end

    t0 = max(tx0, ty0, tz0, tmin)
    t1 = min(tx1, ty1, tz1, tmax)

    t0 < t1
end

# Data structures:
# Hitable:
#   Contains enough data an object to use in BVH building. Contains:
#   - A bounding box
#   - An index to the owning triangle mesh
#   - An index to the triangle inside the mesh
struct Hitable
    box::AABB
    triangleindex::UInt32
    instanceindex::UInt32
end

# BVHNode:
#   Is a node in the hierarchical BVH tree. Is either a leaf node or a parent node.
#   - Bounding box
#   - Mesh index/left index
#   - Triangle index/right index

struct BVHNode
    box::AABB
    left::UInt32
    right::UInt32
    isleafnode::Bool

    function BVHNode(hitable::Hitable)
        left = hitable.instanceindex
        right = hitable.triangleindex
        new(hitable.box, left, right, true)
    end

    function BVHNode(nodes::Vector{BVHNode}, leftindex::UInt32, rightindex::UInt32)
        leftbox = nodes[leftindex].box
        rightbox = nodes[rightindex].box
        box = AABB(leftbox, rightbox)
        new(box, leftindex, rightindex, false)
    end

    function BVHNode()
        new(AABB(), UInt32(0), UInt32(0), false)
    end
end

@inline isleaf(node::BVHNode) :: Bool = node.isleafnode
@inline getinstanceindex(node::BVHNode) :: UInt32 = node.left
@inline gettriangleindex(node::BVHNode) :: UInt32 = node.right
@inline getleft(node::BVHNode) :: UInt32 = node.left
@inline getright(node::BVHNode) :: UInt32 = node.right

# BVHAcceleration:
#   Contains all accelration data and can be ask if it was hit.
#   Contains only enough information to find the leaf triangle index, but does
#   not contain the triangle data itself. It must be used in conjuction with a World
#   object.
#   Contains:
#   - A list of BVHNodes

struct BVHAcceleration
    nodes::CuDeviceArray{BVHNode, 1, CUDAnative.AS.Global}
end

# VisitationStack is a stack wrapper around a thread local array
#   - array: A thread local array of fixed size
#   - len: length
#   Methods:
#       + push!: a node index
#       + pop!: a node index
mutable struct VisitationStack
    array::MVector{30, UInt32}
    len::UInt32

    VisitationStack() = new(zeros(MVector{30, UInt32}), 0)
end

@inline function pushnode!(s::VisitationStack, index::UInt32)
    s.len += 1
    @inbounds s.array[s.len] = index
end

@inline function popnode!(s::VisitationStack) :: UInt32
    s.len -= 1
    @inbounds s.array[s.len + 1]
end

@inline isempty(s::VisitationStack) = s.len == 0

# Building root node algorithm:
# Input:
#   - List of all Hitables in the scene
#   - PRNG
# Output: A BVHAcceleration object
# Algorithm:
#   - If only one object, create a leaf node:
#       + Create leaf node with bounding box, triangle mesh index, triangle index
#       + Return leaf node
#   - Randomly decide on an axis to sort on
#   - Sort all input Hitables by its axis value (x, y, or z)
#   - Split in half
#   - Call itself recursively on left and right
#   - Create BVHNode with left and right
#   - Return BVHNode

function allocatebvhnode(nodes::Vector{BVHNode}) :: UInt32
    push!(nodes, BVHNode())
    length(nodes)
end

function allocatebvhnode(nodes::Vector{BVHNode}, node::BVHNode) :: UInt32
    push!(nodes, node)
    length(nodes)
end

function buildbvhnode(hitables::AbstractVector{Hitable}, nodes::Vector{BVHNode}, start::UInt32, last::UInt32, rng::UniformRNG) :: UInt32
    isthisaleaf = start == last
    if isthisaleaf
        hitable = hitables[start]
        node = BVHNode(hitable)
        leafindex = allocatebvhnode(nodes, node)
        leafindex
    else
        # Allocate an empty node for this parent node, so it
        # ends up in the list before the children. That way, the first
        # node is the root node.
        thisindex = allocatebvhnode(nodes)

        # Choose axis randomly to sort by.
        axischoice = next(rng)
        byaxis = if axischoice < 0.33f0
            h -> h.box.min.x
        elseif axischoice < 0.66f0
            h -> h.box.min.y
        else
            h -> h.box.min.z
        end

        # Sort the hitables vector from start to last, based on the randomly
        # chosen axis above.
        hitablesview = view(hitables, start:last)
        sort!(hitablesview, by=byaxis)

        # Split the remaining nodes into two parts
        n = last - start + 1
        leftlast = UInt32(start + ceil(UInt32, n / 2) - 1)
        rightstart = UInt32(leftlast + 1)

        # Recursively create the child nodes.
        leftindex = buildbvhnode(hitables, nodes, start, leftlast, rng)
        rightindex = buildbvhnode(hitables, nodes, rightstart, last, rng)

        # Finally create the parent node and write it to the pre-allocated
        # node index.
        thisnode = BVHNode(nodes, leftindex, rightindex)
        nodes[thisindex] = thisnode

        thisindex
    end
end

function buildacceleration(hitables::AbstractVector{Hitable}, rng::UniformRNG) :: AbstractVector{BVHNode}
    nodes = Vector{BVHNode}()

    buildbvhnode(hitables, nodes, UInt32(1), UInt32(length(hitables)), rng)

    nodes
end

# Hit algorithm:
# Input:
#   - Ray
#   - BVHAcceleration
#   - World
#   - VisitationStack
# Output: HitRecord
#
# Algorithm:
#   - Store nearest hit, initialized to none
#   - While local list is not empty:
#       - Get next node
#       - If leaf node:
#           + If is a hit, based on bounding box:
#               * Has mesh index and triangle index
#               * Get transformed Triangle from World
#               * If Triangle is hit by Ray:
#                   - If nearer than so far, store as nearest hit
#       - If parent node:
#           + If is a hit:
#               * Push index to left node
#               * Push index to right node
#   - Return nearest hit record

@inline function hitacceleration(acceleration::BVHAcceleration, tmin::Float32, tmax::Float32, ray::Ray, world::World) :: HitRecord
    visit = VisitationStack()
    pushnode!(visit, UInt32(1))

    closesthit = HitRecord()

    while !isempty(visit)
        # Get the next BVHNode to visit.
        thisindex = popnode!(visit)
        @inbounds node = acceleration.nodes[thisindex]

        if hitbox(node.box, ray, tmin, tmax)
            if isleaf(node)
                instanceindex = getinstanceindex(node)
                triangleindex = gettriangleindex(node)
                triangle = gettriangle(world, instanceindex, triangleindex)
                thisrec = hittriangle(triangle, tmin, tmax, ray)

                # Store the hit if it was closer than the current closest hit.
                if thisrec.ishit && thisrec.t < closesthit.t
                    closesthit = thisrec
                end
            else
                pushnode!(visit, getleft(node))
                pushnode!(visit, getright(node))
            end
        end
    end

    closesthit
end

end


text:
@testset "utils" begin
    @testset "single_dim_size" begin

        @test (1,1,1,3) ==IndexFunArrays.single_dim_size(4, 3)

        @test (1,1,1,5) == IndexFunArrays.single_dim_size(4, 5)

        @test (1,5) == IndexFunArrays.single_dim_size(2, 5)
    end

    @testset "Test own size method" begin
        x = ones((2,4,6,8, 10));
        @test selectsizes(x, (2, 3, 4)) == (1, 4, 6, 8, 1)
        @test selectsizes(x, (4, 3, 2), keep_dims = false) == (8, 6, 4)

        x = ones((10));
        @test selectsizes(x, (1,), keep_dims=false) == (10,)
        @test selectsizes(x, (1,), keep_dims=true) == (10,)
        @test selectsizes(x, 1, keep_dims=false) == (10,)
        @test selectsizes(x, 1, keep_dims=true) == (10,)

        x = ones((10, 10));
        @test selectsizes(x, (1, 2), keep_dims=false) == (10, 10)
        @test selectsizes(x, (1, 2), keep_dims=true) == (10, 10)

        x = ones((1,1,3));
        @test selectsizes(x, (1,2), keep_dims=false) == (1,1)
        @test selectsizes(x, (3,3), keep_dims=false) == (3,3)
        @test selectsizes(x, (3,3), keep_dims=true) == (1, 1, 3)

    end

end


text:
println("==============flag! MOI_wrapper>nlp.jl======================================")
"""
MOI NLPBlock
"""

MOI.supports(::Optimizer, ::MOI.NLPBlock) = true

function MOI.set(model::Optimizer, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData)
    model.nlp_data = nlp_data
    return
end


text:
using MCpack
using Test

@testset "MCpack.jl" begin
    # Write your tests here.
end


text:
### A Pluto.jl notebook ###
# v0.12.18

using Markdown
using InteractiveUtils

# ╔═╡ 373d0f40-265b-11eb-24a5-55e733605914
begin
	using Pkg
    Pkg.activate(pwd())
	Pkg.instantiate()
end

# ╔═╡ 864115a0-20d4-11eb-16ff-71d198f9342e
begin
	using Plots
	using LinearAlgebra
	using Distributions

	gr();
end

# ╔═╡ a33c14a0-20c2-11eb-0bf8-837173d9d91c
md"""
# Assignment 2: Implement Gradient Descent
### Aadam (CS1945)
### CSE 662 - Deep Learning by Dr. Masroor Hussain

---

## Assignment Statement:
Write down the code to optimize $f(x) = \frac{1}{2} \|Ax-b\|^2_2$ using gradient descent method in C/C++/Python/Matlab. See equation 4.21.

---
"""

# ╔═╡ cf5c9eb0-20c2-11eb-05ba-2f5dd23fdaf4
md"""
# Gradient Descent
Gradient Descent is a well-known iterative optimization algorithm, widely used in Machine Learning and Deep Learning. It is a way to minimize an objective function $f(x)$ parameterized by a model's parameters $x$ by updating the parameters in the opposite direction of the gradient of the objective function $\nabla_xf(x)$ w.r.t the parameters.

For this assignment, we'll use the following algorithm, from the Deep Learning book:
![algo](https://i.imgur.com/XSwselr.png)
"""

# ╔═╡ b6ed7820-20ce-11eb-0584-9130387e39ef
md"""
# Solution:
"""

# ╔═╡ 2782f610-262d-11eb-0017-6314af9f27a0
md"""
First, we'll randomly initialize $A$ (a matrix) and $b$ (a vector) from a normal distribution. $A$ will be of size $m \times n$ and $b$ is a vector of length $m$.
"""

# ╔═╡ 667a0680-2595-11eb-2859-258f6407b790
rnorm = Normal();

# ╔═╡ 39217590-262d-11eb-0858-5106bc2498d8
m, n = 2, 5;

# ╔═╡ e009fab0-2594-11eb-1711-475e4d2bd945
𝐀 = rand(rnorm, m, n)

# ╔═╡ 7e8c2280-2595-11eb-1eed-cd25325faa84
𝒃 = rand(rnorm, m)

# ╔═╡ b9f447b0-262d-11eb-0444-4582e4bb3e37
md"""
Now, let's define our objective function, i.e. $f(x) = \frac{1}{2} \|Ax-b\|^2_2$.
"""

# ╔═╡ 0a73bcbe-2584-11eb-398f-57c8bc991795
𝑓(𝒙) = 1/2 * norm(𝐀 * 𝒙 .- 𝒃)^2

# ╔═╡ 0a59cc20-2584-11eb-3ac8-d3f4c1a325ec
𝑓(0)

# ╔═╡ e3246890-262d-11eb-1ab3-eb263ec2ae51
md"""
The gradient of this function would be:
$
\nabla_xf(x) = A^T(Ax - b) = A^TAx - A^Tb
$
"""

# ╔═╡ 0a3db8a0-2584-11eb-0508-d935b4b32358
∇ₓ𝑓(𝒙) = 𝐀'𝐀 * 𝒙 .- 𝐀'𝒃

# ╔═╡ 3f354f30-25c7-11eb-2d7b-29d7a44aadd7
∇ₓ𝑓(0)

# ╔═╡ 81edab80-262e-11eb-129c-55f41ac9c662
md"""
Let's define our step size $(\epsilon)$ and tolerance $(\delta)$.
"""

# ╔═╡ 09bdda40-2584-11eb-0a1d-5535b8cbfffb
ϵ = 0.01;

# ╔═╡ 10c38d80-2598-11eb-1c62-3bc6c09ef195
δ = 1e-10;

# ╔═╡ 3e2aa142-265c-11eb-2731-5bb42f1843f5
md"""
Here, $x$ is a vector of length $n$, starting from an arbitrary value, randomly taken from a normal distribution.
"""

# ╔═╡ 5d45dce0-25b5-11eb-3ea7-4964249ce9a5
𝒙 = rand(rnorm, n)

# ╔═╡ 10df79f0-2598-11eb-39a5-476ac6e96ff5
function gradient_descent(𝒙)
	i = 0
	fxvals = []
	while norm(∇ₓ𝑓(𝒙)) > δ && i < 10000
		push!(fxvals, 𝑓(𝒙))
		𝒙 = 𝒙 - ϵ * (∇ₓ𝑓(𝒙))
		i += 1
		println("$i : $𝒙")
	end
	𝒙, fxvals
end

# ╔═╡ 58c08130-25c8-11eb-0efd-81830bcfaecd
minx, fxvals = gradient_descent(𝒙);

# ╔═╡ 5a03b4c0-25b6-11eb-3ad0-753d1aec5cfb
𝑓(minx)

# ╔═╡ cf29ec50-265c-11eb-1796-53065b4536e4
plot(
	fxvals,
	title = "Value of f(x)",
	label = "f(x)",
	xaxis="Steps",
	yaxis="f(x)",
	lw = 2
)

# ╔═╡ 645934a0-20fb-11eb-1750-cf7ca28ab9e6
md"""
## Conclusion
In this notebook, we showed how we can use **Gradient Descent** algorithm to find the optimum value of $x$ that minimizes the function $f(x) = \frac{1}{2} \|Ax-b\|^2_2$.
"""

# ╔═╡ Cell order:
# ╟─a33c14a0-20c2-11eb-0bf8-837173d9d91c
# ╟─cf5c9eb0-20c2-11eb-05ba-2f5dd23fdaf4
# ╟─b6ed7820-20ce-11eb-0584-9130387e39ef
# ╠═373d0f40-265b-11eb-24a5-55e733605914
# ╠═864115a0-20d4-11eb-16ff-71d198f9342e
# ╟─2782f610-262d-11eb-0017-6314af9f27a0
# ╠═667a0680-2595-11eb-2859-258f6407b790
# ╠═39217590-262d-11eb-0858-5106bc2498d8
# ╠═e009fab0-2594-11eb-1711-475e4d2bd945
# ╠═7e8c2280-2595-11eb-1eed-cd25325faa84
# ╟─b9f447b0-262d-11eb-0444-4582e4bb3e37
# ╠═0a73bcbe-2584-11eb-398f-57c8bc991795
# ╠═0a59cc20-2584-11eb-3ac8-d3f4c1a325ec
# ╟─e3246890-262d-11eb-1ab3-eb263ec2ae51
# ╠═0a3db8a0-2584-11eb-0508-d935b4b32358
# ╠═3f354f30-25c7-11eb-2d7b-29d7a44aadd7
# ╟─81edab80-262e-11eb-129c-55f41ac9c662
# ╠═09bdda40-2584-11eb-0a1d-5535b8cbfffb
# ╠═10c38d80-2598-11eb-1c62-3bc6c09ef195
# ╟─3e2aa142-265c-11eb-2731-5bb42f1843f5
# ╠═5d45dce0-25b5-11eb-3ea7-4964249ce9a5
# ╠═10df79f0-2598-11eb-39a5-476ac6e96ff5
# ╠═58c08130-25c8-11eb-0efd-81830bcfaecd
# ╠═5a03b4c0-25b6-11eb-3ad0-753d1aec5cfb
# ╠═cf29ec50-265c-11eb-1796-53065b4536e4
# ╟─645934a0-20fb-11eb-1750-cf7ca28ab9e6


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#
using DiscreteEvents, Random, Distributions

println("... basic tests: processes ...")
simex = DiscreteEvents.PrcException(DiscreteEvents.Stop(), nothing)
@test simex.event == DiscreteEvents.Stop()
@test simex.value === nothing

# ===== test process registration
ch1 = Channel(32)
ch2 = Channel(32)

resetClock!(𝐶)
incr(𝐶, c1::Channel, c2::Channel, a) = (a+1, yield())
a = [1,1,3.0,3.0,"A","A","A","A"]
b = [1,2,3.0,nextfloat(3.0),"A","A#1","A#2","A#3"]
for i in 1:8
    @test process!(Prc(a[i], incr, ch1, ch2, 1)) == b[i]
end
for i in 1:8
    @test 𝐶.processes[b[i]].id == b[i]
end
@test process!(Prc((1,2), incr, ch1, ch2, 1)) == (1,2)
@test process!(Prc((1,2), incr, ch1, ch2, 1)).state === :failed # Julia 1.0

println("... test channel 4 example ...")
A = []

function simple(c::Clock, input::Channel, output::Channel, name, id, op)
    token = take!(input)         # take something from the input
    push!(A, (tau(c), name, id, token))
    delay!(c, rand())           # after a delay
    put!(output, op(token, id))  # put it out with some op applied
end

resetClock!(𝐶)
Random.seed!(123)

for i in 1:2:8    # create, register and start 8 SimProcesses Prc
    process!(Prc(i, simple, ch1, ch2, "foo", i, +))
    process!(Prc(i+1, simple, ch2, ch1, "bar", i+1, *))
end

@test length(𝐶.processes) == 8
for p in values(𝐶.processes)
    @test istaskstarted(p.task)
end

put!(ch1, 1)
sleep(0.01)
run!(𝐶, 10)

@test length(A) > 20
p = [i[3] for i in A]
for i in 1:8
    @test i ∈ p  # all processes did something
end

println("... test errors and stop ...")
schedule(𝐶.processes[8].task, ErrorException, error=true)
sleep(0.1)
@test 𝐶.processes[8].task.state == :failed
delete!(𝐶.processes, 8)

for p in values(𝐶.processes)
    @test istaskstarted(p.task)
    DiscreteEvents.stop!(p, DiscreteEvents.Stop())
    @test istaskdone(p.task)
end

println("... test wait! ...")

ch1 = Channel(1)
ch2 = Channel(1)

a = 1
b = 1
res = []
incb() = global b +=1
checktime(clk, x) = tau(clk) ≥ x
checka(x) = a == x
checkb(x) = b ≥ x

function testwait(clk::Clock, c1::Channel, c2::Channel)
    wait!(clk, (fun(checktime, clk, 2), fun(checka, 1)))
    push!(res, (tau(clk), 1, a, b))
    wait!(clk, fun(isa, a, Int)) # must return immediately
    push!(res, (tau(clk), 2, a, b))
    periodic!(clk, fun(incb))
    wait!(clk, fun(checkb, 201))
    push!(res, (tau(clk), 3, a, b))
    take!(c1)
end

resetClock!(𝐶)
process!(Prc(1, testwait, ch1, ch2))

run!(𝐶, 10)
r = [i[1] for i in res]
@test r[1] ≈ 2
@test r[2] ≈ 2
@test r[3] ≈ 4
@test res[3][4] == 201
@test b == 801

a = [1]
function testdelay(clk::Clock, x)
    delay!(clk, at, 2)     # assertion error bad timing
    x[1] += 10
end

function testdelay2(clk::Clock, x)
    delay!(clk, until, 5)
    x[1] += 1
end

function testdelay3(clk::Clock, x)
    delay!(clk, until, Normal(5,1))
    x[1] += 1
end

resetClock!(𝐶)
process!(Prc(1, testdelay, a), 3)
process!(Prc(2, testdelay2, a), 3)
process!(Prc(3, testdelay3, a), 3)
run!(𝐶, 10)

@test 𝐶.processes[1].task.state == :failed
@test a[1] == 7

function testnow(c, x)
    delay!(c, 1);
    x[1] += 1;
    println(c, "$(tau(c)): ", "x is $(x[1])")
end

resetClock!(𝐶)
process!(Prc(1, testnow, a), 3)
run!(𝐶, 5)
@test a[1] == 10


text:
# Activate environment

using ClusterManagers,Distributed

# Add in the cores allocated by the scheduler as workers
addprocs(SlurmManager(parse(Int,ENV["SLURM_NTASKS"])-1))
print("Added workers: ")
println(nworkers())

# Load required packages and helpers on all processes
@everywhere begin
    using DistributedArrays,Statistics
    include("../parallelhelpers.jl")
end

# Load the file names
dataLoc = "../../../../data/word_count/";
fnames = dataLoc.*readdir(dataLoc)

# Create a distributed array with the file names
dfnames = distribute(fnames)

# Calculate the word counts
# Broadcast over a distributed array automatically does the operation on the correct worker
allcounts =  getcounts.(dfnames)

# Collect the distributed array onto the leader process
allcountsL = convert(Array,allcounts)

# Calculate the overall counts
overallcounts = getoverallcount(allcountsL)

# Share the overall counts with each of the workers
@everywhere overallcounts=$overallcounts

# Calculate the normalized counts on each worker
normcount = getnormcount.(allcounts)

# Print out the top 5 words for each file and their counts on the leader process
for i=1:length(normcount)
    println(sort(collect(normcount[i]), by=x->x[2], rev=true)[1:min(5,length(normcount[i])),:])
end

max_stars_repo_name: llsc-supercloud/teaching-examples

max_issues_repo_name: llsc-supercloud/teaching-examples

max_forks_repo_name: llsc-supercloud/teaching-examples


text:
using Base.Threads, Distributed
addprocs(Sys.CPU_THREADS - 1)
@everywhere using SharedArrays
@everywhere using GR, LinearAlgebra, Random, Distributions
@everywhere function pointsDelaunay(dimBase::Float64, dim::Int64)
    points = rand(dim, 2)
    if dimBase != 1.0
        points[:,2] /= dimBase
    end
    n, tri = GR.delaunay(points[:,1], points[:,2])
    n, tri, points
end
@everywhere function triDelaunay(points)
    n, tri = GR.delaunay(points[:,1], points[:,2])
    n, tri, points
end
@everywhere function rcentersGeqOneLune(beta::Float64, x1::Float64,
        y1::Float64, x2::Float64, y2::Float64)
    c1 = [beta.*x1./2.0 .+ (1.0 .- beta./2).*x2 beta.*y1./2.0 .+ (1.0 .- beta./2.0).*y2]
    c2 = [(1.0 .- beta./2.0).*x1 .+ beta.*x2./2.0 (1.0 .- beta./2).*y1 .+ beta.*y2./2]
    r = beta.*sqrt((x2 .- x1).^2 + (y2 .- y1).^2)./2.0
    r, c1, c2
end
@everywhere function rcentersLeqOne(beta::Float64, x1::Float64, y1::Float64,
        x2::Float64, y2::Float64)
    dot_v1 = [x1, y1]
    dot_v2 = [x2, y2]
    vec_v1v2 = dot_v2 .- dot_v1
    rotvec = reverse(vec_v1v2).*[-1,1]
    radOut = sqrt.(sum(vec_v1v2.^2))./(2*beta)
    b_out = sqrt.(1 .- beta^2)./(2*beta)
    p_c1 = (dot_v1 .+ dot_v2)./2 .+ b_out.*rotvec
    p_c2 = (dot_v1 .+ dot_v2)./2 .- b_out.*rotvec
    radOut, p_c1', p_c2'
end
@everywhere function rcentersGeqOneCircle(beta::Float64, x1::Float64,
        y1::Float64, x2::Float64, y2::Float64)
    dot_v1 = [x1, y1]
    dot_v2 = [x2, y2]
    vec_v1v2 = dot_v2 .- dot_v1
    rotvec = reverse(vec_v1v2).*[-1,1]
    radOut = beta.*sqrt.(sum(vec_v1v2.^2))./2
    b_out = sqrt.(beta^2 .- 1)./2
    p_c1 = (dot_v1 .+ dot_v2)./2 .+ b_out.*rotvec
    p_c2 = (dot_v1 .+ dot_v2)./2 .- b_out.*rotvec
    radOut, p_c1', p_c2'
end
@everywhere function betaRandSkeletonGeqOneLune(beta::Float64, dim::Int64,
        plotSkeleton = false, dimBase::Float64 = 1.0)
    n, tri, points = pointsDelaunay(dimBase, dim)
    mAdjacency = diagm(0 => randn(dim))
    indices = [1, 1, 2]
    j = [2, 3, 3]
    for trindx in 1:n
        for (indx, i) in enumerate(indices)
            indx1 = tri[trindx,:][i]
            indx2 = tri[trindx,:][j[indx]]
            x1, y1 = points[indx1, :]
            x2, y2 = points[indx2, :]
            r, c1, c2 = rcentersGeqOneLune(beta, x1, y1, x2, y2)
            full_1 = sqrt.(sum((points .- c1).^2, dims = 2)) .< r
            full_2 = sqrt.(sum((points .- c2).^2, dims = 2)) .< r
            full_1[[indx1, indx2]] .= false
            if any(full_1 .+ full_2 .== 2) == false
                rnd = randn()./sqrt(2)
                mAdjacency[indx1, indx2] = rnd
                mAdjacency[indx2, indx1] = rnd
            end
        end
    end
    plotSkeleton == true ? (mAdjacency, points) : mAdjacency
end
@everywhere function betaRandSkeletonGeqOneCircle(beta::Float64, dim::Int64,
        plotSkeleton = false, dimBase::Float64 = 1.0)
    n, tri, points = pointsDelaunay(dimBase, dim)
    mAdjacency = diagm(0 => randn(dim))
    indices = [1, 1, 2]
    j = [2, 3, 3]
    for trindx in 1:n
        for (indx, i) in enumerate(indices)
            indx1 = tri[trindx,:][i]
            indx2 = tri[trindx,:][j[indx]]
            x1, y1 = points[indx1, :]
            x2, y2 = points[indx2, :]
            r, c1, c2 = rcentersGeqOneCircle(beta, x1, y1, x2, y2)
            full_1 = sqrt.(sum((points .- c1).^2, dims = 2)) .< r
            full_2 = sqrt.(sum((points .- c2).^2, dims = 2)) .< r
            full_1[[indx1, indx2]] .= false
            full_2[[indx1, indx2]] .= false
            if any(full_1 .| full_2) == false
                rnd = randn()./sqrt(2)
                mAdjacency[indx1, indx2] = rnd
                mAdjacency[indx2, indx1] = rnd
            end
        end
    end
    plotSkeleton == true ? (mAdjacency, points, tri) : mAdjacency
end
@everywhere function betaRandSkeletonLeqOne(beta::Float64, dim::Int64, plotSkeleton = false, dimBase::Float64 = 1.0)
    mAdjacency = diagm(0 => randn(dim))
    points = rand(dim, 2)
    if dimBase != 1.0
        points[:,2] /= dimBase
    end
    for indx1 in 1:dim-1
        for indx2 in indx1 + 1:dim
            x1, y1 = points[indx1, :]
            x2, y2 = points[indx2, :]
            r, c1, c2 = rcentersLeqOne(beta, x1, y1, x2, y2)
            full_1 = sqrt.(sum((points .- c1).^2, dims = 2)) .< r
            full_2 = sqrt.(sum((points .- c2).^2, dims = 2)) .< r
            full_1[[indx1, indx2]] .= false
            if any(full_1 .+ full_2 .== 2) == false
                rnd = randn()./sqrt(2)
                mAdjacency[indx1, indx2] = rnd
                mAdjacency[indx2, indx1] = rnd
            end
        end
    end
    plotSkeleton == true ? (mAdjacency, points) : mAdjacency
end
@everywhere function betaSkeletonGeqOneLune(beta::Float64, points, plotSkeleton = false)
    n, tri, points = triDelaunay(points)
    dim = length(points[:,1])
    mAdjacency = diagm(0 => randn(dim))
    indices = [1, 1, 2]
    j = [2, 3, 3]
    for trindx in 1:n
        for (indx, i) in enumerate(indices)
            indx1 = tri[trindx,:][i]
            indx2 = tri[trindx,:][j[indx]]
            x1, y1 = points[indx1, :]
            x2, y2 = points[indx2, :]
            r, c1, c2 = rcentersGeqOneLune(beta, x1, y1, x2, y2)
            full_1 = sqrt.(sum((points .- c1).^2, dims = 2)) .< r
            full_2 = sqrt.(sum((points .- c2).^2, dims = 2)) .< r
            full_1[[indx1, indx2]] .= false
            if any(full_1 .+ full_2 .== 2) == false
                rnd = randn()./sqrt(2)
                mAdjacency[indx1, indx2] = rnd
                mAdjacency[indx2, indx1] = rnd
            end
        end
    end
    plotSkeleton == true ? (mAdjacency, points) : mAdjacency
end

@everywhere function betaSkeletonGeqOneCircle(beta::Float64, points, plotSkeleton = false)
    n, tri, points = triDelaunay(points)
    dim = length(points[:,1])
    mAdjacency = diagm(0 => randn(dim))
    indices = [1, 1, 2]
    j = [2, 3, 3]
    for trindx in 1:n
        for (indx, i) in enumerate(indices)
            indx1 = tri[trindx,:][i]
            indx2 = tri[trindx,:][j[indx]]
            x1, y1 = points[indx1, :]
            x2, y2 = points[indx2, :]
            r, c1, c2 = rcentersGeqOneCircle(beta, x1, y1, x2, y2)
            full_1 = sqrt.(sum((points .- c1).^2, dims = 2)) .< r
            full_2 = sqrt.(sum((points .- c2).^2, dims = 2)) .< r
            full_1[[indx1, indx2]] .= false
            full_2[[indx1, indx2]] .= false
            if any(full_1 .| full_2) == false
                rnd = randn()./sqrt(2)
                mAdjacency[indx1, indx2] = rnd
                mAdjacency[indx2, indx1] = rnd
            end
        end
    end
    plotSkeleton == true ? (mAdjacency, points, tri) : mAdjacency
end
@everywhere function betaSkeletonLeqOne(beta::Float64, points)
    dim = length(points[:,1])
    mAdjacency = diagm(0 => rand(dim))
    for indx1 in 1:dim-1
        for indx2 in indx1 + 1:dim
            x1, y1 = points[indx1, :]
            x2, y2 = points[indx2, :]
            r, c1, c2 = rcentersLeqOne(beta, x1, y1, x2, y2)
            full_1 = sqrt.(sum((points .- c1).^2, dims = 2)) .< r
            full_2 = sqrt.(sum((points .- c2).^2, dims = 2)) .< r
            full_1[[indx1, indx2]] .= false
            if any(full_1 .+ full_2 .== 2) == false
                rnd = randn()./sqrt(2)
                mAdjacency[indx1, indx2] = rnd
                mAdjacency[indx2, indx1] = rnd
            end
        end
    end
    mAdjacency
end
@everywhere function betaSkeleton(beta::Float64, dim::Int64, dimBase::Float64 = 1.0)
    if beta < 1.0
        madea, points = betaRandSkeletonLeqOne(beta, dim, true, dimBase)
        return points, madea
    elseif beta >= 1.0
        madeaLune, points =  betaRandSkeletonGeqOneLune(beta, dim, true, dimBase)
        madeaCircle = betaSkeletonGeqOneCircle(beta, points)
        return points, madeaLune, madeaCircle
    end
end

@everywhere function betaSkeletonV2(beta::Float64, dim::Int64, dimBase::Float64 = 1.0)
    if beta < 1.0
        madea, points = betaRandSkeletonLeqOne(beta, dim, true, dimBase)
        return points, madea
    elseif beta >= 1.0
        madeaLune, points =  betaRandSkeletonGeqOneLune(beta, dim, true, dimBase)
        madeaCircle = betaSkeletonGeqOneCircle(beta, points)
        return points, madeaLune, madeaCircle
    end
end

function pointsHexagonX(center, size, i)
    angledeg = 60*i #- 30
    anglerad = pi/180*angledeg
    (center[1] .+ size*cos(anglerad), center[2] .+ size*sin(anglerad))
end

max_stars_repo_name: lazarusA/betaSkeletons

max_issues_repo_name: lazarusA/betaSkeletons

max_forks_repo_name: lazarusA/betaSkeletons


text:
function getoperator(op, t)
    if op isa Tuple
        op = Semiring(op...)
    end
    #Default Semiring should be LOR_LAND for boolean
    if op == Semirings.PLUS_TIMES
        if t == Bool
            op = Semirings.LOR_LAND
        end
    end
    #Default BinaryOp should be LAND or LOR for boolean
    if op == BinaryOps.TIMES
        if t == Bool
            op = BinaryOps.LAND
        end
    elseif op == BinaryOps.PLUS
        if t == Bool
            op = BinaryOps.LOR
        end
    end
    #Default monoid should be LAND/LOR
    if op == Monoids.TIMES_MONOID
        if t == Bool
            op = Monoids.LAND_MONOID
        end
    elseif op == Monoids.PLUS_MONOID
        if t == Bool
            op = Monoids.LOR_MONOID
        end
    end

    if op isa AbstractOp
        return op[t]
    else
        return op
    end
end

getaccum(::Nothing, t) = C_NULL
getaccum(op::Function, t) = getoperator(BinaryOp(op), t)
getaccum(op, t) = getoperator(op, t)

_isloaded(op::AbstractOp) = !isempty(getfield(op, :typedops))
_isloaded(op::Union{AbstractSelectOp, AbstractDescriptor}) = getfield(op, :p) != C_NULL
"""
    validtypes(operator::AbstractOp)::Vector{DataType}
    validtypes(operator::SelectOp)::Nothing

Determine the types available as a domain for a particular operator.
Each operator is defined on a specific set of types, for instance the [`LAND_LOR`](@ref)
    semiring is only defined for `Boolean` arguments.

When applied to an `AbstractSelectOp` this will return `nothing`.
When applied to certain operators like positional semirings it will return `[Any]`.
"""
function validtypes(o::AbstractOp)
    if !_isloaded(o)
        _load(o)
    end
    return collect(keys(o.typedops))
end

function Base.getindex(o::AbstractOp, t::DataType)
    _isloaded(o) || _load(o)
    if Any ∈ keys(o.typedops)
        getindex(o.typedops, Any)
    else
        getindex(o.typedops, t)
    end
end

function Base.show(io::IO, ::MIME"text/plain", o::AbstractOp)
    print(io, o.name, ": ", validtypes(o))
end

juliaop(op...) = nothing


text:
@testset "GJS_BWG_HNE_DL" begin
    θ₀ = [0.7, 0.9, 3.3, 5.0, 4.0, 0.3, 4, 2.7, 0.55, 0.26]
    ω = 1.1
    @testset "Preallocation" begin
        model = GJS_BWG_HNE_DL{1,40}(θ₀)
        @test model.r_over_s == model.r/model.s
        @test model.ωₚ² == model.ωₚ^2
        @test model.ωₚ³ == model.ωₚ^3
        @test model.ωₚˢ == model.ωₚ^model.s
        @test model.ωₚˢ⁻¹ == model.ωₚ^(model.s-1)
        @test model.logγ == log(model.γ)
        @test model.s⁻¹ == 1/(model.s)
        @test model.logωₚminuss⁻¹ == log(model.ωₚ) - 1/model.s
        @test model.σᵣ_over3 == model.σᵣ/3
        @test model.cosϕₘ == cos(model.ϕₘ)
        @test model.sinϕₘ == sin(model.ϕₘ)
        @test model.cos2ϕₘ == cos(2model.ϕₘ)
        @test model.sin2ϕₘ == sin(2model.ϕₘ)
    end
    @testset "Gradient" begin
        @test approx_gradient(θ -> sdf(GJS_BWG_HNE_DL{1,40}(θ), ω), θ₀)      ≈ grad_sdf(GJS_BWG_HNE_DL{1,40}(θ₀), ω)
        @test approx_gradient(θ -> sdf(GJS_BWG_HNE_DL{1,40}(θ), 0), θ₀)      ≈ grad_sdf(GJS_BWG_HNE_DL{1,40}(θ₀), 0)
    end
    @testset "Error handling" begin
        @test_throws ArgumentError GJS_BWG_HNE_DL{1,40}(-0.7,0.9, 3.3, 5.0, 4.0, 0, 4, 2.7, 0.55, 0.26)
        @test_throws ArgumentError GJS_BWG_HNE_DL{1,40}(0.7,-0.9, 3.3, 5.0, 4.0, 0, 4, 2.7, 0.55, 0.26)
        @test_throws ArgumentError GJS_BWG_HNE_DL{1,40}(0.7, 0.9, 0.8, 5.0, 4.0, 0, 4, 2.7, 0.55, 0.26) # 0.8 <1
        @test_throws ArgumentError GJS_BWG_HNE_DL{1,40}(0.7, 0.9, 3.3, 0.9, 4.0, 0, 4, 2.7, 0.55, 0.26) # 0.9 < 1
        @test_throws ArgumentError GJS_BWG_HNE_DL{1,40}(0.7, 0.9, 3.3, 5.0,-4.0, 0, 4, 2.7, 0.55, 0.26)
        @test_throws ArgumentError GJS_BWG_HNE_DL{1,40}(0.7, 0.9, 3.3, 5.0, 4.0, 0,-4, 2.7, 0.55, 0.26)
        @test_throws ArgumentError GJS_BWG_HNE_DL{1,40}(0.7, 0.9, 3.3, 5.0, 4.0, 0, 4,-2.7, 0.55, 0.26)
        @test_throws ArgumentError GJS_BWG_HNE_DL{1,40}(0.7, 0.9, 3.3, 5.0, 4.0, 0, 4, 2.7,-0.55, 0.26)
        @test_throws ArgumentError GJS_BWG_HNE_DL{1,40}(0.7, 0.9, 3.3, 5.0, 4.0, 0, 4, 2.7, 0.55,-0.26)
        @test_throws ArgumentError GJS_BWG_HNE_DL{1,40}(ones(9))
        @test_throws ArgumentError GJS_BWG_HNE_DL{1,40}(ones(11))
        @test_throws ErrorException GJS_BWG_HNE_DL(ones(10))
        @test_throws ErrorException GJS_BWG_HNE_DL(0.7, 0.9, 3.3, 5.0, 4.0, 0, 4, 2.7, 0.55, 0.26)
    end
end


text:
using PyPlot, Revise, HiQGA.transD_GP, Random, Statistics
## waveform and times
ramp =  [  # -0.0400066666667    0.5
			-0.0400000000000    0.0
			-0.0399933333333   -0.5
			-0.0200066666667   -0.5
			-0.0200000000000    0.0
			-0.0199933333333    0.5
			-0.0000066666667    0.5
			 0.0000000000000    0.0
			 0.0000066666667   -0.5]

times = vec(10 .^mean(log10.([
            0.0000066667	0.0000200000
			0.0000333333	0.0000466667
			0.0000600000	0.0000733333
			0.0000866667	0.0001266667
			0.0001400000	0.0002066667
			0.0002200000	0.0003400000
			0.0003533333	0.0005533333
			0.0005666667	0.0008733333
			0.0008866667	0.0013533333
			0.0013666667	0.0021000000
			0.0021133333	0.0032733333
			0.0032866667	0.0051133333
			0.0051266667	0.0079933333
			0.0080066667	0.0123933333
			0.0124066667	0.0199933333]), dims=2))
## Set up operator
ntimesperdecade = 10
nfreqsperdecade = 5
nkᵣeval = 60
# x,y,z in z downward!!
# +x is flight direction
zTx = -120
zRx = -80
x_rx = -115.
y_rx = 0.
# these are in z upward, ga-aem convention
rx_roll = 0.
rx_pitch = 0.
rx_yaw = 0.
tx_roll = 0. # GA AEM downward pointing dipole is roll 180
tx_pitch = 90.
tx_yaw = 0.
# model
zfixed   = [-1e6,   0, 40, 60]
rho      = [1e12,   20, 1, 100]
## do it
xmin, xmax, dx = -4.25, 4.25, 0.25
ymin, ymax, dy = -3.25, 3.25, 0.25
xrx = xmin:dx:xmax
yrx = ymin:dy:ymax
u = zeros(length(yrx), length(xrx), length(times))
v = zeros(length(yrx), length(xrx), length(times))
w = zeros(length(yrx), length(xrx), length(times))
for (jx, x) in enumerate(xrx), (iy, y) in enumerate(yrx)
	tempest = transD_GP.TEMPEST1DInversion.Bfield(
  					  ntimesperdecade = ntimesperdecade,
                      nfreqsperdecade = nfreqsperdecade,
					  nkᵣeval = nkᵣeval,
                      times  = times,
                      ramp   = ramp,
					  zTx    = zTx,
					  zRx    = zRx,
					  x_rx   = x,
					  y_rx   = y,
					  rx_roll = rx_roll,
					  rx_pitch = rx_pitch,
					  rx_yaw = rx_yaw,
					  tx_roll = tx_roll,
					  tx_pitch = tx_pitch,
					  tx_yaw = tx_yaw,
					  strictgeometry = false)
	transD_GP.TEMPEST1DInversion.getfieldTD!(tempest, zfixed, rho)
	# Roll180 underneath to get back to Z down co-ordinate system, easier to plot
	# and make sense
	u[iy,jx,:] = tempest.Hx
	v[iy,jx,:] = -tempest.Hy
	w[iy,jx,:] = -tempest.Hz
end
## plot it
f, ax = plt.subplots(1,3, sharex=true, sharey=true, figsize=(16,4))

sleeptime = 1
ntimes = 15
for it = 1:ntimes
	im1 = ax[1].pcolormesh(xrx,yrx, u[:,:,it])
	cb1 = colorbar(im1,ax=ax[1])
	ax[1].quiver(xrx, yrx, u[:,:,it], zeros(size(u,1), size(u,2)), width=0.005)
	ax[1].set_xlabel("x m")
	ax[1].set_ylabel("y m")
	ax[1].set_title("Hx")
	im2 = ax[2].pcolormesh(xrx,yrx, v[:,:,it])
	cb2 = colorbar(im2,ax=ax[2])
	ax[2].quiver(xrx, yrx, zeros(size(v,1), size(v,2)), v[:,:,it], width=0.005)
	ax[2].set_xlabel("x m")
	ax[2].set_ylabel("y m")
	ax[2].set_title("Hy")
	im3 = ax[3].pcolormesh(xrx,yrx, w[:,:,it])
	cb3 = colorbar(im3,ax=ax[3])
	ax[3].set_xlabel("x m")
	ax[3].set_ylabel("y m")
	ax[3].set_title("Hz")
	plt.draw()
	sleep(sleeptime)
	if it != ntimes
		ax[1].clear()
		ax[2].clear()
		ax[3].clear()
		cb1.remove()
		cb2.remove()
		cb3.remove()
	end
end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
using Documenter, Bio, Pkg

makedocs(
    format = Documenter.HTML(),
    sitename = "Bio.jl",
    authors = replace(join(Pkg.TOML.parsefile("Project.toml")["authors"], ", "), r" <.*?>" => "" ),
    pages = [
        "Home"    => "index.md"
    ],

)

deploydocs(
    repo = "github.com/BioJulia/Bio.jl.git",
    push_preview = true,
    deps = nothing,
    make = nothing
)


text:
function fVel(x,Param)
    str = lower(Param.ProfVel)
    if str == "solidbody"
      uS=0;
      vS=0;
    elseif str == "hyperdiffcart"
      if abs(x[1]-3000. * 1.e3)<=1000 * 1.e3 && abs(x[2])<=1000 * 1.e3
        uS=sqrt(5000);
        vS=sqrt(5000);
      else
        uS=0;
        vS=0;
      end
    elseif str == "hyperdiff"
      (Lon,Lat,R)=cart2sphere(x[1],x[2],x[3]);
      GreatCircleR=acos(sin(Param.PertLat)*sin(Lat)+
        cos(Param.PertLat)*cos(Lat)*cos(Lon-Param.PertLon));
      GreatCircleR=GreatCircleR/Param.PertExpR;
      # Apply perturbation in zonal velocity
      if GreatCircleR<6.0
        uS=sqrt(Param.ValDiff);
        vS=sqrt(Param.ValDiff);
      else
        uS=0.0;
        vS=0.0;
      end
    elseif str == "barowavesphere"
      (Lon,Lat,R)=cart2sphere(x[1],x[2],x[3]);
      Z=max(R-Param.RadEarth,0);
  #     if Z>27500 && abs(Lat-pi/4)<=.1
  #       aa=3;
  #     end
      T0=0.5*(Param.T0E+Param.T0P);
      ConstA=1.0/Param.LapseRate;
      ConstB=(T0-Param.T0P)/(T0*Param.T0P);
      ConstC=0.5*(Param.K+2.0)*(Param.T0E-Param.T0P)/(Param.T0E*Param.T0P);
      ConstH=Param.Rd*T0/Param.Grav;
      ScaledZ=Z/(Param.B*ConstH);
      Tau1=ConstA*Param.LapseRate/T0*exp(Param.LapseRate/T0*Z) +
        ConstB*(1.0-2.0*ScaledZ*ScaledZ)*exp(-ScaledZ*ScaledZ);
      Tau2=ConstC*(1.0-2.0*ScaledZ*ScaledZ)*exp(-ScaledZ*ScaledZ);
      IntTau2=ConstC*Z*exp(-ScaledZ*ScaledZ);
      if Param.Deep
        RRatio= R/Param.EarthRadius;
      else
        RRatio = 1.0;
      end
      InteriorTerm=(RRatio*cos(Lat))^Param.K -
        Param.K/(Param.K+2.0)*(RRatio*cos(Lat))^(Param.K+2.0);
      Temperature=1.0/(RRatio*RRatio)/(Tau1-Tau2*InteriorTerm);

      InteriorTermU=(RRatio*cos(Lat))^(Param.K-1.0) -
         (RRatio*cos(Lat))^(Param.K+1.0);
      BigU=Param.Grav/Param.RadEarth*Param.K *
        IntTau2*InteriorTermU*Temperature;
      if Param.Deep
        RCosLat=R*cos(Lat);
      else
        RCosLat=Param.RadEarth*cos(Lat);
      end
      OmegaRCosLat =Param.Omega*RCosLat;

      # 		if (dOmegaRCosLat * dOmegaRCosLat + dRCosLat * dBigU < 0.0) {
      # 			_EXCEPTIONT("Negative discriminant detected.");
      # 		}

      uS=-OmegaRCosLat+sqrt(OmegaRCosLat*OmegaRCosLat+RCosLat*BigU);
      vS=0;
      # Exponential perturbation
      GreatCircleR=acos(sin(Param.PertLat)*sin(Lat) +
        cos(Param.PertLat)*cos(Lat)*cos(Lon-Param.PertLon));

      GreatCircleR=GreatCircleR/Param.PertExpR;

      # Tapered perturbation with height
      if Z < Param.PertZ
        PertTaper=1.0-3.0*Z*Z/(Param.PertZ * Param.PertZ)+
           2.0*Z*Z*Z/(Param.PertZ*Param.PertZ*Param.PertZ);
      else
        PertTaper=0.0;
      end

      # Apply perturbation in zonal velocity
      if GreatCircleR<1.0
        uSPert=Param.Up*PertTaper*exp(-GreatCircleR*GreatCircleR);
      else
        uSPert=0.0;
      end
      uS=uS+uSPert;

    elseif str == "barowavecart"
      eta=EtaFromZ(x[1],x[2],x[3],Param);
      uS=-Param.u0*sin(pi*x[2]/Param.Ly)^2*log(eta) *
        exp(-(log(eta)/Param.b)^2);
      if x[2]<0
        uS=-uS;
      else
        uS=uS+Param.uP*exp(-((x[1]-Param.xC)^2  +
          (x[2]-Param.yC)^2)/Param.Lp^2);
      end
      vS=0;
    elseif str == "baldauf"
      uS=0;
      vS=0;
    elseif str == "barowavesphere"
  #     # parameters

    elseif str == "gravityhill"
      uS=Param.uMax;
      vS=0;
    elseif str == "sphericalharmonics"
      (lon,lat,r)=cart2sphere(x[1],x[2],x[3]);
      (uS,vS)=harmonicVec(Param.lHar,Param.mHar,lat,lon,r);
    elseif str == "galewsky"
      (lon,lat,r)=cart2sphere(x[1],x[2],x[3]);
      if (lat<=Param.lat0G) || (lat>=Param.lat1G)
        uS=0;
      else
        uS=Param.uM/Param.eN*exp(1.0/((lat-Param.lat0G)*(lat-Param.lat1G)));
      end
      vS=0;
    elseif str == "rossbyhaurwitz"
      omega6=7.8480e-6;
      K6=7.8480e-6;
      R6=4.0;
      (lon,lat,r)=cart2sphere(x[1],x[2],x[3]);
      uS=r*omega6*cos(lat) + r*K6*(cos(lat))^(R6-1.0) *
        (R6*sin(lat)*sin(lat)-cos(lat)*cos(lat))*cos(R6*lon);
       vS=-r*K6*R6*(cos(lat))^(R6-1.0)*sin(lat)*sin(R6*lon);
    elseif str == "hill"
      uMax=20;
      rotation_angle=0;
      (lon,lat,r)=cart2sphere(x[1],x[2],x[3]);
      (rot_lon,rot_lat)=Rotate(lon,lat,rotation_angle);
      u_lat=uMax*cos(rot_lat);
      v_lat=0.0e0;
      if abs(rotation_angle)<1.0E-8
        uS = u_lat;
        vS = v_lat;
      else
        # rotate wind components
        (uS,vS) = turnwi(u_lat,v_lat,lon,lat,rot_lon,rot_lat,0.0e0,-0.5e0*pi+rotation_angle);
        if abs(uS)<1.0E-10
          uS=0.0e0;
        end
      end
    elseif str == "bickley"
      U1 = cosh(x[2])^(-2);

      # ?? = exp(-(x2 + Param.l / 10)^2 / 2p.l^2) * cos(p.k * x1) * cos(p.k * x2)
      # Vortical velocity fields
      gaussian = exp(-(x[2] + Param.l / 10)^2 / (2*Param.l^2));
      u1 = gaussian * (x[2] + Param.l / 10) / Param.l^2 * cos(Param.k * x[1]) * cos(Param.k * x[2]);
      u1 = u1+Param.k * gaussian * cos(Param.k * x[1]) * sin(Param.k * x[2]);
      u2 = -Param.k * gaussian * sin(Param.k * x[1]) * cos(Param.k * x[2]);
      uS=U1+Param.Eps*u1;
      vS=Param.Eps*u2;
    elseif str == "cosinebell"
      (lon,lat,r)=cart2sphere(x[1],x[2],x[3]);
      alpha0=pi/2;
      u0=2*pi*Param.RadEarth/(86400*12);
      uS=u0*(cos(alpha0)*cos(lat)+sin(alpha0)*cos(lon)*sin(lat));
      vS =-u0*sin(alpha0)*sin(lon);
    elseif str == "const"
      uS=Param.uMax;
      vS=Param.vMax;
    elseif str == "linear"
      uS=x[3];
      vS=0;
    end
    return (uS,vS)
end

# function [rot_lon,rot_lat]=Rotate(lon,lat,rotation_angle)
# if abs(rotation_angle)<1.0e-8
#   rot_lon = lon;
#   rot_lat = lat;
# else
#   [rot_lon,rot_lat]=regrot(lon,lat,0.0e0,-0.5e0*pi+rotation_angle);
# end
# end

# function [pxrot,pyrot]=regrot(pxreg,pyreg,pxcen,pycen)

# #----------------------------------------------------------------------
# #
# #*    conversion between regular and rotated spherical coordinates.
# #*
# #*    pxreg     longitudes of the regular coordinates
# #*    pyreg     latitudes of the regular coordinates
# #*    pxrot     longitudes of the rotated coordinates
# #*    pyrot     latitudes of the rotated coordinates
# #*              all coordinates given in degrees n (negative for s)
# #*              and degrees e (negative values for w)
# #*    pxcen     regular longitude of the south pole of the rotated grid
# #*    pycen     regular latitude of the south pole of the rotated grid
# #*
# #*    kcall=-1: find regular as functions of rotated coordinates.
# #*    kcall= 1: find rotated as functions of regular coordinates.
# #
# #-----------------------------------------------------------------------
# #

# zpih = pi*0.5d0;
# zsycen = sin((pycen+zpih));
# zcycen = cos((pycen+zpih));

# zxmxc  = pxreg - pxcen;
# zsxmxc = sin(zxmxc);
# zcxmxc = cos(zxmxc);
# zsyreg = sin(pyreg);
# zcyreg = cos(pyreg);
# zsyrot = zcycen*zsyreg - zsycen*zcyreg*zcxmxc;
# zsyrot = max(zsyrot,-1.e0);
# zsyrot = min(zsyrot,+1.e0);

# pyrot = asin(zsyrot);

# zcyrot = cos(pyrot);
# zcxrot = (zcycen*zcyreg*zcxmxc +zsycen*zsyreg)/zcyrot;
# zcxrot = max(zcxrot,-1.e0);
# zcxrot = min(zcxrot,+1.e0);
# zsxrot = zcyreg*zsxmxc/zcyrot;

# pxrot = acos(zcxrot);

# if zsxrot<0.0
#   pxrot = -pxrot;
# end
# end

# function [pures,pvres]=turnwi(puarg,pvarg,
#                   pxreg,pyreg,pxrot,pyrot,
#                   pxcen,pycen)
# #
# #-----------------------------------------------------------------------
# #
# #*    turn horizontal velocity components between regular and
# #*    rotated spherical coordinates.
# #
# #*    puarg : input u components
# #*    pvarg : input v components
# #*    pures : output u components
# #*    pvres : output v components
# #*    pa    : transformation coefficients
# #*    pb    :    -"-
# #*    pc    :    -"-
# #*    pd    :    -"-
# #*    pxreg : regular longitudes
# #*    pyreg : regular latitudes
# #*    pxrot : rotated longitudes
# #*    pyrot : rotated latitudes
# #*    kxdim              : dimension in the x (longitude) direction
# #*    kydim              : dimension in the y (latitude) direction
# #*    kx                 : number of gridpoints in the x direction
# #*    ky                 : number of gridpoints in the y direction
# #*    pxcen              : regular longitude of the south pole of the
# #*                         transformed grid
# #*    pycen              : regular latitude of the south pole of the
# #*                         transformed grid
# #*
# #*    kcall < 0          : find wind components in regular coordinates
# #*                         from wind components in rotated coordinates
# #*    kcall > 0          : find wind components in rotated coordinates
# #*                         from wind components in regular coordinates
# #*    note that all coordinates are given in degrees n and degrees e.
# #*       (negative values for s and w)
# #
# #-----------------------------------------------------------------------

#      zpih = pi*0.5e0;
#      zsyc = sin(pycen+zpih);
#      zcyc = cos(pycen+zpih);

#      zsxreg = sin(pxreg);
#      zcxreg = cos(pxreg);
#      zsyreg = sin(pyreg);
#      zcyreg = cos(pyreg);

#      zxmxc  = pxreg - pxcen;
#      zsxmxc = sin(zxmxc);
#      zcxmxc = cos(zxmxc);

#      zsxrot = sin(pxrot);
#      zcxrot = cos(pxrot);
#      zsyrot = sin(pyrot);
#      zcyrot = cos(pyrot);

#      pa = zcxmxc*zcxrot + zcyc*zsxmxc*zsxrot;
#      pb = zcyc*zsxmxc*zcxrot*zsyrot + zsyc*zsxmxc*zcyrot -
#           zcxmxc*zsxrot*zsyrot;
#      pc =-zsyc*zsxrot/zcyreg;
#      pd = (zcyc*zcyrot - zsyc*zcxrot*zsyrot)/zcyreg;
#      !
#      pures = pa*puarg + pb*pvarg;
#      pvres = pc*puarg + pd*pvarg;
# end


text:
using RNGPool

include("simpleModel.jl")
include("demo.jl")

setRNGs(12345)

println("Using a Simple Model with n = 1000, boxSize = 1.0")
simpleModel, simplelM = setupSimpleModel(1000, 1.0)

## CCBPF couples after a number of steps, less with larger values of N
visualizeCCSMC(simpleModel, simplelM, 2, :BS, true)
visualizeCCSMC(simpleModel, simplelM, 2, :BS, true, true)
visualizeCCSMC(simpleModel, simplelM, 4, :BS, true)
visualizeCCSMC(simpleModel, simplelM, 4, :BS, true, true)
visualizeCCSMC(simpleModel, simplelM, 8, :BS, true)
visualizeCCSMC(simpleModel, simplelM, 8, :BS, true, true)
visualizeCCSMC(simpleModel, simplelM, 16, :BS, true)
visualizeCCSMC(simpleModel, simplelM, 16, :BS, true, true)
visualizeCCSMC(simpleModel, simplelM, 32, :BS, true)
visualizeCCSMC(simpleModel, simplelM, 32, :BS, true, true)
visualizeCCSMC(simpleModel, simplelM, 64, :BS, true)
visualizeCCSMC(simpleModel, simplelM, 64, :BS, true, true)
visualizeCCSMC(simpleModel, simplelM, 128, :BS, true)
visualizeCCSMC(simpleModel, simplelM, 128, :BS, true, true)

## CCPF essentially will struggle to couple with 128 particles
visualizeCCSMC(simpleModel, 128, true, false, 100)
visualizeCCSMC(simpleModel, 128, true, true, 100)
## CCPF essentially will couple eventually in one shot with 2^14 particles
visualizeCCSMC(simpleModel, 2^14, true, false, 100)
visualizeCCSMC(simpleModel, 2^14, true, true, 100)

## CCAPF performs well for N ∈ {2,4} but struggles for moderate values
visualizeCCSMC(simpleModel, simplelM, 2, :AS, true)
visualizeCCSMC(simpleModel, simplelM, 2, :AS, true, true)
visualizeCCSMC(simpleModel, simplelM, 4, :AS, true)
visualizeCCSMC(simpleModel, simplelM, 4, :AS, true, true)
visualizeCCSMC(simpleModel, simplelM, 8, :AS, true)
visualizeCCSMC(simpleModel, simplelM, 8, :AS, true, true)
visualizeCCSMC(simpleModel, simplelM, 16, :AS, true, true, 1000)
visualizeCCSMC(simpleModel, simplelM, 64, :AS, true, true, 1000)
visualizeCCSMC(simpleModel, simplelM, 256, :AS, true, true, 100)
visualizeCCSMC(simpleModel, simplelM, 1024, :AS, true, true, 100)
visualizeCCSMC(simpleModel, simplelM, 4096, :AS, true, true, 100)
visualizeCCSMC(simpleModel, simplelM, 2^14, :AS, true, true, 100)


text:
"""
`RandomODeProblem` - structure to solve the random ODE

V(x)u - Δu = f

on (a, b), where A = Δ and V(x) is a random potential. Mutability allows for ease of updating the potential V(x).
"""
mutable struct RandomODEProblem{TM, TV, TLBC, TRBC, TF}
    a::Real # Left endpoint
    b::Real # Right endpoint
    nₓ::Int # Total number of interior mesh points
    Δx::Real # Mesh spacing
    x::Vector # Mesh points including boundary
    m::TM # Mean of (log of) random coefficient
    σ::Real # Standard deviation of (log of) random coefficient
    s::Real # Parameter associated with covariance operator
    V::TV # Potential
    left_bc::TLBC # Left boundary condition
    right_bc::TRBC # Right boundary condition
    A::SparseMatrixCSC # Matrix for LHS (Laplacian w/ random perturbation on diagonal)
    f::TF # RHS function
    rhs::Vector # RHS vector
    u::Vector # Solution of the problem
end

"""
`DirichletBC` - structure to house a DirichletBC boundary condition, i.e., a boundary condition for which

u(a) = uₗ, u(b) = uᵣ
"""
struct DirichletBC
    boundary_condition::Real
end


text:
module CoreMLProtobuf
  const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
  module Specification
    const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
    include("ArrayFeatureExtractor_pb.jl")
    include("BayesianProbitRegressor_pb.jl")
    include("FeatureTypes_pb.jl")
    include("DataStructures_pb.jl")
    include("CategoricalMapping_pb.jl")
    include("CustomModel_pb.jl")
    include("DictVectorizer_pb.jl")
    include("FeatureVectorizer_pb.jl")
    include("GLMClassifier_pb.jl")
    include("GLMRegressor_pb.jl")
    include("Identity_pb.jl")
    include("Imputer_pb.jl")
    include("ItemSimilarityRecommender_pb.jl")
    include("Parameters_pb.jl")
    include("LinkedModel_pb.jl")
    include("NearestNeighbors_pb.jl")
    include("NeuralNetwork_pb.jl")
    include("Normalizer_pb.jl")
    include("OneHotEncoder_pb.jl")
    include("Scaler_pb.jl")
    include("NonMaximumSuppression_pb.jl")
    include("SVM_pb.jl")
    include("TreeEnsemble_pb.jl")
    module CoreMLModels
      const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
      include("Gazetteer_pb.jl")
      include("VisionFeaturePrint_pb.jl")
      include("TextClassifier_pb.jl")
      include("WordTagger_pb.jl")
      include("WordEmbedding_pb.jl")
      include("SoundAnalysisPreprocessing_pb.jl")
    end
    include("Model_pb.jl")
  end
end