text:
#=
Time steps for the MixedSolver module.
Since the FE and FV portions may use different steppers, this must work for both,
but they can be done separately by only passing one set of variables.
For all of them:
input
- fe_var
- fv_var
- vol_lhs
- vol_rhs
- surf_lhs
- surf_rhs
- step_info = (t, dt, dofs_per_node_fe, dofs_per_loop_fe, assemble_func_fe, dofs_per_node_fv, dofs_per_loop_fv, assemble_func_fv)
- fe_storage = (A, b, fe_sol, allocated_vecs, b_order, b_sizes)
- fv_storage = (fv_sol, allocated_vecs_fv)
=#

function mixed_euler_explicit_step(fe_var, fv_var, vol_lhs, vol_rhs, surf_lhs, surf_rhs, t, step_info, fe_storage, fv_storage)
    fe_stepper = step_info[1];
    fv_stepper = step_info[2];
    dofs_per_node_fe = step_info[3];
    dofs_per_loop_fe = step_info[4];
    fe_assemble_func = step_info[5];
    dofs_per_node_fv = step_info[6];
    dofs_per_loop_fv = step_info[7];
    fv_assemble_func = step_info[8];
    dt = fe_stepper.dt;
    assemble_t = 0;
    linsolve_t = 0;

    if !(fv_var === nothing) && config.num_partitions > 1 exchange_ghosts(fv_var, fv_grid, i); end
    pre_step_function();

    ### First step the FE part #######################################################################################
    if !(fe_var === nothing)
        # Unpack things
        A = fe_storage[1];
        b = fe_storage[2];
        fe_sol = fe_storage[3];
        allocated_vecs = fe_storage[4];
        b_order = fe_storage[5];
        b_sizes = fe_storage[6];

        N1 = size(grid_data.allnodes,2);

        assemble_t += @elapsed begin
                                    b = assemble(fe_var, nothing, vol_rhs[1], allocated_vecs, dofs_per_node_fe, dofs_per_loop_fe, t, dt; rhs_only = true, assemble_loops=fe_assemble_func);
                                    b = gather_system(nothing, b, N1, dofs_per_node_fe, b_order, b_sizes);
                                end

        linsolve_t += @elapsed(tmpvec = distribute_solution( linear_system_solve(A,b) , N1, dofs_per_node_fe, b_order, b_sizes));

        # At this point tmpvec holds the boundary values
        # directly write them to the variable values and zero sol.
        copy_bdry_vals_to_variables(fe_var, tmpvec, grid_data, dofs_per_node_fe, zero_vals=true);

        for i=1:length(fe_sol)
            fe_sol[i] = fe_sol[i] + dt * tmpvec[i];
        end

        copy_bdry_vals_to_vector(fe_var, fe_sol, grid_data, dofs_per_node_fe);
        place_sol_in_vars(fe_var, fe_sol);

    end

    ### Then step the FV part ######################################################################################
    if !(fv_var === nothing)
        fv_sol = fv_storage[1];
        fv_allocated_vecs = fv_storage[2];

        assemble_t += @elapsed(tmpvec = fv_assemble(fv_var, vol_lhs[2], vol_rhs[2], surf_lhs[2], surf_rhs[2], fv_allocated_vecs, dofs_per_node_fv, dofs_per_loop_fv, t, dt, assemble_loops=fv_assemble_func));

        for i=1:length(fv_sol)
            fv_sol[i] = fv_sol[i] + dt * tmpvec[i];
        end

        FV_copy_bdry_vals_to_vector(fv_var, fv_sol, fv_grid, dofs_per_node_fv);
        place_sol_in_vars(fv_var, fv_sol);
    end

    post_step_function();

    return (assemble_t, linsolve_t);
end

function mixed_euler_implicit_step(fe_var, fv_var, vol_lhs, vol_rhs, surf_lhs, surf_rhs, t, step_info, fe_storage, fv_storage)
    fe_stepper = step_info[1];
    fv_stepper = step_info[2];
    dofs_per_node_fe = step_info[3];
    dofs_per_loop_fe = step_info[4];
    fe_assemble_func = step_info[5];
    dofs_per_node_fv = step_info[6];
    dofs_per_loop_fv = step_info[7];
    fv_assemble_func = step_info[8];
    dt = fe_stepper.dt;
    assemble_t = 0;
    linsolve_t = 0;

    if !(fv_var === nothing) && config.num_partitions > 1 exchange_ghosts(fv_var, fv_grid, i); end
    pre_step_function();

    ### First step the FE part #######################################################################################
    if !(fe_var === nothing)
        # Unpack things
        A = fe_storage[1];
        b = fe_storage[2];
        fe_sol = fe_storage[3];
        allocated_vecs = fe_storage[4];
        b_order = fe_storage[5];
        b_sizes = fe_storage[6];

        N1 = size(grid_data.allnodes,2);

        assemble_t += @elapsed begin
                                    b = assemble(fe_var, nothing, vol_rhs[1], allocated_vecs, dofs_per_node_fe, dofs_per_loop_fe, t, dt; rhs_only = true, assemble_loops=fe_assemble_func);
                                    b = gather_system(nothing, b, N1, dofs_per_node_fe, b_order, b_sizes);
                                end

        linsolve_t += @elapsed(tmpvec = distribute_solution( linear_system_solve(A,b) , N1, dofs_per_node_fe, b_order, b_sizes));

        # At this point tmpvec holds the boundary values
        # directly write them to the variable values and zero sol.
        copy_bdry_vals_to_variables(fe_var, tmpvec, grid_data, dofs_per_node_fe, zero_vals=true);

        for i=1:length(fe_sol)
            fe_sol[i] = tmpvec[i];
        end

        copy_bdry_vals_to_vector(fe_var, fe_sol, grid_data, dofs_per_node_fe);
        place_sol_in_vars(fe_var, fe_sol);
    end

    ### Then step the FV part ######################################################################################
    if !(fv_var === nothing)
        fv_sol = fv_storage[1];
        fv_allocated_vecs = fv_storage[2];

        ## TODO

        FV_copy_bdry_vals_to_vector(fv_var, fv_sol, fv_grid, dofs_per_node_fv);
        place_sol_in_vars(fv_var, fv_sol);
    end

    post_step_function();

    return (assemble_t, linsolve_t);
end

function mixed_crank_nicholson_step(fe_var, fv_var, vol_lhs, vol_rhs, surf_lhs, surf_rhs, t, step_info, fe_storage, fv_storage)
    fe_stepper = step_info[1];
    fv_stepper = step_info[2];
    dofs_per_node_fe = step_info[3];
    dofs_per_loop_fe = step_info[4];
    fe_assemble_func = step_info[5];
    dofs_per_node_fv = step_info[6];
    dofs_per_loop_fv = step_info[7];
    fv_assemble_func = step_info[8];
    dt = fe_stepper.dt;
    assemble_t = 0;
    linsolve_t = 0;

    if !(fv_var === nothing) && config.num_partitions > 1 exchange_ghosts(fv_var, fv_grid, i); end
    pre_step_function();

    ### First step the FE part #######################################################################################
    if !(fe_var === nothing)
        # Unpack things
        A = fe_storage[1];
        b = fe_storage[2];
        fe_sol = fe_storage[3];
        allocated_vecs = fe_storage[4];
        b_order = fe_storage[5];
        b_sizes = fe_storage[6];

        N1 = size(grid_data.allnodes,2);

        assemble_t += @elapsed begin
            b = assemble(fe_var, nothing, vol_rhs[1], allocated_vecs, dofs_per_node_fe, dofs_per_loop_fe, t, dt; rhs_only = true, assemble_loops=fe_assemble_func);
            b = gather_system(nothing, b, N1, dofs_per_node_fe, b_order, b_sizes);
        end

        linsolve_t += @elapsed(tmpvec = distribute_solution( linear_system_solve(A,b) , N1, dofs_per_node_fe, b_order, b_sizes));

        # At this point tmpvec holds the boundary values
        # directly write them to the variable values and zero sol.
        copy_bdry_vals_to_variables(fe_var, tmpvec, grid_data, dofs_per_node_fe, zero_vals=true);

        for i=1:length(fe_sol)
        fe_sol[i] = tmpvec[i];
        end

        copy_bdry_vals_to_vector(fe_var, fe_sol, grid_data, dofs_per_node_fe);
        place_sol_in_vars(fe_var, fe_sol);
    end

    ### Then step the FV part ######################################################################################
    if !(fv_var === nothing)
        fv_sol = fv_storage[1];
        fv_allocated_vecs = fv_storage[2];

        ## TODO

        FV_copy_bdry_vals_to_vector(fv_var, fv_sol, fv_grid, dofs_per_node_fv);
        place_sol_in_vars(fv_var, fv_sol);
    end

    post_step_function();

    return (assemble_t, linsolve_t);
end

function mixed_lsrk4_step(fe_var, fv_var, vol_lhs, vol_rhs, surf_lhs, surf_rhs, t, step_info, fe_storage, fv_storage)
    fe_stepper = step_info[1];
    fv_stepper = step_info[2];
    dofs_per_node_fe = step_info[3];
    dofs_per_loop_fe = step_info[4];
    fe_assemble_func = step_info[5];
    dofs_per_node_fv = step_info[6];
    dofs_per_loop_fv = step_info[7];
    fv_assemble_func = step_info[8];
    dt = fe_stepper.dt;
    assemble_t = 0;
    linsolve_t = 0;

    ### First step the FE part #######################################################################################
    if !(fe_var === nothing)
        # Unpack things
        A = fe_storage[1];
        b = fe_storage[2];
        fe_sol = fe_storage[3];
        allocated_vecs = fe_storage[4];
        b_order = fe_storage[5];
        b_sizes = fe_storage[6];
        tmppi = fe_storage[7];
        tmpki = fe_storage[8];

        N1 = size(grid_data.allnodes,2);

        # Low storage RK4:
        # p0 = u
        #   ki = ai*k(i-1) + dt*f(p(i-1), t+ci*dt)
        #   pi = p(i-1) + bi*ki
        # u = p5

        tmppi = get_var_vals(fe_var, tmppi);
        for rki=1:fe_stepper.stages
            rktime = t + fe_stepper.c[rki]*fe_stepper.dt;
            # p(i-1) is currently in u
            pre_step_function();

            ### First step the FE part #######################################################################################
            assemble_t += @elapsed begin
                                b = assemble(fe_var, nothing, vol_rhs[1], allocated_vecs, dofs_per_node_fe, dofs_per_loop_fe, rktime, fe_stepper.dt; rhs_only = true, assemble_loops=fe_assemble_func);
                                b = gather_system(nothing, b, N1, dofs_per_node_fe, b_order, b_sizes);
                            end

            linsolve_t += @elapsed(fe_sol = distribute_solution( linear_system_solve(A,b) , N1, dofs_per_node_fe, b_order, b_sizes));

            # At this point sol holds the boundary values
            # directly write them to the variable values and zero sol.
            copy_bdry_vals_to_variables(fe_var, fe_sol, grid_data, dofs_per_node_fe, zero_vals=true);

            if rki == 1 # because a1 == 0
                tmpki = fe_stepper.dt .* fe_sol;
            else
                tmpki = fe_stepper.a[rki].*tmpki + fe_stepper.dt.*fe_sol;
            end
            tmppi = tmppi + fe_stepper.b[rki].*tmpki

            copy_bdry_vals_to_vector(fe_var, tmppi, grid_data, dofs_per_node_fe);
            place_sol_in_vars(fe_var, tmppi);

            post_step_function();
        end

    end

    ### Then step the FV part ######################################################################################
    if !(fv_var === nothing)
        fv_sol = fv_storage[1];
        fv_allocated_vecs = fv_storage[2];
        fv_tmppi = fv_storage[3];
        fv_tmpki = fv_storage[4];

        fv_tmppi = get_var_vals(fv_var, fv_tmppi);
        for rki=1:fv_stepper.stages
            rktime = t + fv_stepper.c[rki]*fv_stepper.dt;
            # p(i-1) is currently in u
            if config.num_partitions > 1 exchange_ghosts(fv_var, fv_grid, i); end
            pre_step_function();

            assemble_t += @elapsed(fv_sol = fv_assemble(fv_var, vol_lhs[2], vol_rhs[2], surf_lhs[2], surf_rhs[2], fv_allocated_vecs, dofs_per_node_fv, dofs_per_loop_fv, rktime, dt, assemble_loops=fv_assemble_func));

            if rki == 1 # because a1 == 0
                fv_tmpki = fv_stepper.dt .* fv_sol;
            else
                fv_tmpki = fv_stepper.a[rki].*fv_tmpki + fv_stepper.dt.*fv_sol;
            end
            fv_tmppi = fv_tmppi + fv_stepper.b[rki].*fv_tmpki

            FV_copy_bdry_vals_to_vector(fv_var, fv_tmppi, fv_grid, dofs_per_node_fv);
            place_sol_in_vars(fv_var, fv_tmppi);

            post_step_function();
        end
    end

    return (assemble_t, linsolve_t);
end

# This may be removed.
function mixed_rk4_step(fe_var, fv_var, vol_lhs, vol_rhs, surf_lhs, surf_rhs, t, step_info, fe_storage, fv_storage)
    return mixed_multistage_step(fe_var, fv_var, vol_lhs, vol_rhs, surf_lhs, surf_rhs, t, step_info, fe_storage, fv_storage);
end

function mixed_multistage_step(fe_var, fv_var, vol_lhs, vol_rhs, surf_lhs, surf_rhs, t, step_info, fe_storage, fv_storage)
    fe_stepper = step_info[1];
    fv_stepper = step_info[2];
    dofs_per_node_fe = step_info[3];
    dofs_per_loop_fe = step_info[4];
    fe_assemble_func = step_info[5];
    dofs_per_node_fv = step_info[6];
    dofs_per_loop_fv = step_info[7];
    fv_assemble_func = step_info[8];
    dt = fe_stepper.dt;
    assemble_t = 0;
    linsolve_t = 0;

    ### First step the FE part #######################################################################################
    if !(fe_var === nothing)
        # Unpack things
        A = fe_storage[1];
        b = fe_storage[2];
        fe_sol = fe_storage[3];
        allocated_vecs = fe_storage[4];
        b_order = fe_storage[5];
        b_sizes = fe_storage[6];
        tmpresult = fe_storage[7];
        tmpki = fe_storage[8];

        N1 = size(grid_data.allnodes,2);

        # solution will be placed in var.values for each stage
        for stage=1:fe_stepper.stages
            stime = t + fe_stepper.c[stage]*fe_stepper.dt;

            # Update the values in vars to be used in this stage
            if stage > 1
                initialized_tmpresult = false;
                for j=1:stage
                    if fe_stepper.a[stage, j] > 0
                        if !initialized_tmpresult
                            initialized_tmpresult = true;
                            for k=1:length(fe_sol)
                                tmpresult[k] = fe_sol[k] + fe_stepper.dt * fe_stepper.a[stage, j] * tmpki[k,j];
                            end
                        else
                            for k=1:length(fe_sol)
                                tmpresult[k] += fe_stepper.dt * fe_stepper.a[stage, j] * tmpki[k,j];
                            end
                        end
                    end
                end

                if initialized_tmpresult
                    copy_bdry_vals_to_vector(fe_var, tmpresult, grid_data, dofs_per_node_fe);
                    place_sol_in_vars(fe_var, tmpresult);
                end
                post_step_function(); # seems weird, but imagine this is happening after stage-1
            end

            pre_step_function();

            ### First step the FE part #######################################################################################
            assemble_t += @elapsed begin
                                b = assemble(fe_var, nothing, vol_rhs[1], allocated_vecs, dofs_per_node_fe, dofs_per_loop_fe, stime, fe_stepper.dt; rhs_only = true, assemble_loops=fe_assemble_func);
                                b = gather_system(nothing, b, N1, dofs_per_node_fe, b_order, b_sizes);
                            end

            linsolve_t += @elapsed(tmpki[:,stage] = distribute_solution( linear_system_solve(A,b) , N1, dofs_per_node_fe, b_order, b_sizes));

            # At this point tmpki[:,stage] holds the boundary values
            # directly write them to the variable values and zero sol.
            copy_bdry_vals_to_variables(fe_var, tmpki[:,stage], grid_data, dofs_per_node_fe, zero_vals=true);

        end
        for i=1:length(fe_sol)
            for stage=1:fe_stepper.stages
                fe_sol[i] += fe_stepper.dt * fe_stepper.b[stage] * tmpki[i, stage];
            end
        end

        copy_bdry_vals_to_vector(fe_var, fe_sol, grid_data, dofs_per_node_fe);
        place_sol_in_vars(fe_var, fe_sol);

        post_step_function();
    end

    ### Then step the FV part ######################################################################################
    if !(fv_var === nothing)
        fv_sol = fv_storage[1];
        fv_allocated_vecs = fv_storage[2];
        fv_tmpresult = fv_storage[3];
        fv_tmpki = fv_storage[4];

        # solution will be placed in var.values for each stage
        for stage=1:fv_stepper.stages
            stime = t + fv_stepper.c[stage]*fv_stepper.dt;

            # Update the values in vars to be used in this stage
            if stage > 1
                initialized_tmpresult = false;
                for j=1:stage
                    if fv_stepper.a[stage, j] > 0
                        if !initialized_tmpresult
                            initialized_tmpresult = true;
                            for k=1:length(fv_sol)
                                fv_tmpresult[k] = fv_sol[k] + fv_stepper.dt * fv_stepper.a[stage, j] * fv_tmpki[k,j];
                            end
                        else
                            for k=1:length(fv_sol)
                                fv_tmpresult[k] += fv_stepper.dt * fv_stepper.a[stage, j] * fv_tmpki[k,j];
                            end
                        end
                    end
                end

                if initialized_tmpresult
                    FV_copy_bdry_vals_to_vector(fv_var, fv_tmpresult, fv_grid, dofs_per_node_fv);
                    place_sol_in_vars(fv_var, fv_tmpresult);
                end
                post_step_function(); # seems weird, but imagine this is happening after stage-1
            end

            if config.num_partitions > 1 exchange_ghosts(fv_var, fv_grid, i); end
            pre_step_function();

            tmpki[:,stage] = fv_assemble(fv_var, vol_lhs[2], vol_rhs[2], surf_lhs[2], surf_rhs[2], fv_allocated_vecs, dofs_per_node_fv, dofs_per_loop_fv, stime, fv_stepper.dt, assemble_loops=fv_assemble_func);
        end
        for i=1:length(fv_sol)
            for stage=1:fv_stepper.stages
                fv_sol[i] += fv_stepper.dt * fv_stepper.b[stage] * fv_tmpki[i, stage];
            end
        end
        FV_copy_bdry_vals_to_vector(fv_var, fv_sol, fv_grid, dofs_per_node_fv);
        place_sol_in_vars(fv_var, fv_sol);

        post_step_function();
    end

    return (assemble_t, linsolve_t);
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
module TestingQuicjPriorityQueues

using Test

using Random

using Base: @propagate_inbounds, lt,
    ReverseOrdering, Reverse,
    IteratorEltype, HasEltype,
    IteratorSize, HasLength

using QuickHeaps
using QuickHeaps:
    AbstractPriorityQueue, PriorityQueue, FastPriorityQueue,
    AbstractNode, getkey, getval,
    isheap, index, nodes, ordering, in_range, heap_index

function test_queue!(A::AbstractPriorityQueue{K,V},
                     key_list::AbstractVector{K},
                     val_list::AbstractVector{V}) where {K,V}
    # Check arguments.
    axes(val_list) == axes(key_list) || error("incompatible indices")
    n = length(key_list)
    axes(key_list) == (1:n,) || error("non-standard indices")
    o = ordering(A)

    # Dummy private function used to "label" the tests.
    check(flag::Bool, comment) = flag

    # Checks keytype, etc.
    @test keytype(A) == K
    @test keytype(typeof(A)) == K
    @test IteratorSize(keys(A)) == HasLength()
    @test IteratorSize(typeof(keys(A))) == HasLength()
    @test IteratorEltype(keys(A)) == HasEltype()
    @test IteratorEltype(typeof(keys(A))) == HasEltype()
    @test eltype(keys(A)) == K
    @test eltype(typeof(keys(A))) == K

    # Checks valtype, etc.
    @test valtype(A) == V
    @test valtype(typeof(A)) == V
    @test IteratorSize(values(A)) == HasLength()
    @test IteratorSize(typeof(values(A))) == HasLength()
    @test IteratorEltype(values(A)) == HasEltype()
    @test IteratorEltype(typeof(values(A))) == HasEltype()
    @test eltype(values(A)) == V
    @test eltype(typeof(values(A))) == V

    # Check that `enqueue!` maintains the binary-heap structure.
    R = Dict{K,V}() # reference dictionary
    test_1 = true
    test_2 = true
    for (k, v) in zip(key_list, val_list)
        enqueue!(A, k, v)
        test_1 &= isheap(o, nodes(A))
        test_2 &= !haskey(R, k) # unique key?
        R[k] = v
    end
    @test check(test_1, "binary-heap structure after `enqueue!`")
    @test check(test_2, "unique keys after `enqueue!`")
    @test length(A) == n

    # Check `first`, `peek`, `keys`, and `values`.
    k, v = first(A)
    @test v == (isa(o, ReverseOrdering) ? maximum : minimum)(values(A))
    @test v == first(values(A))
    @test k == first(keys(A))
    x = peek(AbstractNode, A)
    @test v == getval(x)
    @test k == getkey(x)

    # Check `getindex` in random order.
    test_1 = true
    for i in randperm(n)
        k = key_list[i]
        test_1 &= (haskey(A,k) && A[k] == R[k])
    end
    @test check(test_1, "`getindex` in random order")

    # Test that `keys` and `values` yield all elements in the same order as
    # `iterate`.
    test_1 = true
    test_2 = true
    for (k,v,kv) in zip(keys(A), values(A), A)
        test_1 &= (k == kv.first)
        test_2 &= (v == kv.second)
    end
    @test check(test_1, "keys in heap order")
    @test check(test_2, "values in heap order")

    # Check `copy`.
    B = copy(A)
    @test typeof(B) === typeof(A)
    @test length(B) == length(A)
    @test nodes(A) == nodes(B)
    @test index(A) == index(B)
    S = Set(keys(B))
    @test length(S) == length(B) # keys are unique?
    test_1 = true
    for k in S
        test_1 &= (haskey(A, k) && haskey(B, k) && A[k] == B[k])
    end
    @test check(test_1, "`copy` yields same nodes")

    # Check `delete!` and result of addressing a non-existing key.
    k = key_list[rand(1:n)]
    delete!(B, k)
    @test !haskey(B, k)
    @test get(B, k, nothing) === nothing
    @test_throws ArgumentError B[k]

    # Check `isempty`, `empty!`, etc.
    empty!(B)
    @test isempty(B)
    @test length(B) == 0
    @test length(A) == n # no side effects on A

    # Check `setindex!` in heap order.
    length(B) > 1 && empty!(B)
    test_1 = true
    test_2 = true
    for i in randperm(n)
        k, v = key_list[i], val_list[i]
        B[k] = v
        test_1 &= isheap(o, nodes(B))
        test_2 &= (haskey(B, k) && B[k] == R[k])
    end
    @test check(test_1, "heap structure preserved by `setindex!`")
    @test check(test_2, "same value for key after `setindex!`")
    @test length(B) == length(A)
    test_1 = true
    test_2 = true
    for i in randperm(n)
        k = key_list[i]
        test_1 &= haskey(B, k)
        test_2 &= B[k] == R[k]
    end
    @test check(test_1, "no missing keys in random order")
    @test check(test_2, "same values in random order")

    # Check that `delete!` maintains the binary-heap structure.
    test_1 = true
    test_2 = true
    test_3 = true
    for i in randperm(n)
        k = key_list[i]
        test_1 &= haskey(B, k)
        delete!(B, k)
        test_2 &= isheap(o, nodes(B))
        test_3 &= !haskey(B, k)
    end
    @test check(test_1, "keys exist before `delete!`")
    @test check(test_2, "heap structure preserved by `delete!`")
    @test check(test_3, "keys do not exist after `delete!`")
    @test isempty(B)

    # Check that `pop!` extracts nodes in order and maintains the binary-heap
    # structure.
    B = copy(A)
    test_1 = true
    test_2 = true
    prev = (isa(o, ReverseOrdering) ? typemax : typemin)(V)
    for i in 1:n
        k, v = pop!(B)
        test_1 &= isheap(o, nodes(B))
        test_2 &= !lt(o, v, prev)
        prev = v
    end
    @test check(test_1, "heap structure preserved by `pop!`")
    @test check(test_2, "`pop!` yields keys in order")
    @test isempty(B)

    # Tests for fast priority queues.
    if isa(A, FastPriorityQueue)
        cartesian_index = CartesianIndices(index(A))
        linear_index = LinearIndices(index(A))
        test_1 = true
        test_2 = true
        test_3 = true
        test_4 = true
        test_5 = true
        test_6 = true
        test_7 = true
        test_8 = true
        for i in randperm(n)
            k = key_list[i]
            c = cartesian_index[k]
            inds = c.I
            test_1 &= haskey(A, c)
            test_2 &= haskey(A, inds)
            test_3 &= (A[k] == A[c])
            test_4 &= (A[k] == A[inds...])
            v = A[k]
            test_5 &= !haskey(delete!(A, c), k)
            A[c] = v
            test_6 &= (haskey(A, k) && A[k] == v)
            test_7 &= !haskey(delete!(A, inds), k)
            A[inds...] = v
            test_8 &= (haskey(A, k) && A[k] == v)
        end
        @test check(test_1, "`haskey` with Cartesian index")
        @test check(test_2, "`haskey` with multi-dimensional index")
        @test check(test_3, "`getindex` with Cartesian index")
        @test check(test_4, "`getindex` with multi-dimensional index")
        @test check(test_5, "`delete!` with Cartesian index")
        @test check(test_6, "`setindex!` with Cartesian index")
        @test check(test_7, "`delete!` with multi-dimensional index")
        @test check(test_8, "`setindex!` with multi-dimensional index")
    end
end

@testset "Priority queues       " begin
    K, V, n = Int, Float64, 20
    key_list = map(K, 1:n)
    val_list = rand(V, n)
    test_queue!(PriorityQueue{K,V}(), key_list, val_list)
    test_queue!(PriorityQueue{K,V}(Reverse), key_list, val_list)
end

@testset "Fast priority queues  " begin
    V, dims = Float32, (2,3,4)
    n = prod(dims)
    m = div(9n + 5, 10) # keep ~90% of indices
    key_list = randperm(n)[1:m]
    val_list = rand(V, m)
    test_queue!(FastPriorityQueue{V}(dims), key_list, val_list)
    test_queue!(FastPriorityQueue{V}(Reverse, dims), key_list, val_list)
end

end # module

ext: jl


text:
# Implements the mutual information objective for the target tracking problem

using LinearAlgebra
using Base.Iterators
using Statistics
using Random

export finite_horizon_information, entropy, TimedObservation

import Base.length

nan_to_zero(x::Float64) = ifelse(isnan(x), 0.0, x)

entropy_sum(prior::AbstractFilter) = sum(x -> nan_to_zero(-x * log(2, x)),
                                         HistogramFilters.get_values(prior))
function entropy(prior::Filter)
  sum = 0.0
  @inbounds @simd for x = HistogramFilters.get_values(prior)
    sum += nan_to_zero(-x * log(2, x))
  end
  sum
end
function entropy(prior::SparseFilter)
  sum = 0.0
  @inbounds @simd for x = nonzeros(HistogramFilters.get_values(prior))
    sum += nan_to_zero(-x * log(2, x))
  end
  sum
end

function entropy(prior::Array{<:AbstractFilter})
  sum(entropy, prior)
end

#
# Measurement representation and updates
#
# Most of the sequential planning utilizes the trajectory objects. I prefer to
# avoid further processing of the trajectories as this code runs a lot.
#
# I also include a more general representation TimedObservation to support
# mutual information for collections of observations at specific times
#

#
# Timed observations allow flexibility in scheduling observations
#
struct TimedObservation
  step::Int64
  state::State
end
# Object that manages and sorts observations
mutable struct TimedObservations
  observations::Vector{TimedObservation}
  index::Int64

  function TimedObservations(observations::Vector{TimedObservation})
    new(sort(observations, by=x->x.step), 1)
  end
end
TimedObservations() = TimedObservations(TimedObservation[])
length(x::TimedObservations) = length(x.observations)
access(x::TimedObservations) = x.observations[x.index]
advance!(x::TimedObservations) = (x.index += 1; nothing)
all_processed(x::TimedObservations) = x.index > length(x)

#
# Combine sets of observations (so that we can obtain the set of all
# observations from the union of the prior and posterior)
#
function combine_observations(x::Vararg{X}) where
  X <: Union{Vector{Trajectory}, Vector{TimedObservation}}

  vcat(x...)
end

#
# Simulate observations and update the filter for a given time-step
#

# Assumption: all trajectories are at least as long as "step"
function simulate_update_filter!(grid::Grid, filter::AbstractFilter,
                                 sensor::RangingSensor, target_state::State,
                                 trajectories::Vector{Trajectory}, step;
                                 rng=Random.GLOBAL_RNG,
                                 likelihood_buffer =
                                  likelihoods_buffer(get_states(grid, filter))
                                )

    for trajectory in trajectories
      robot_state = trajectory[step]

      range_observation = generate_observation(grid, sensor, robot_state,
                                               target_state; rng=rng)

      measurement_update!(filter, robot_state, get_states(grid, filter), sensor,
                          grid, range_observation, buffer=likelihood_buffer)
    end
end
# Precondition: all observations before "step" have already been processed.
# (this is reasonable because the filter update design would prevent processing
# those observations afterward. As such, observations that precede the scope of
# updates are malformed)
function simulate_update_filter!(grid::Grid, filter::AbstractFilter,
                                 sensor::RangingSensor, target_state::State,
                                 observations::TimedObservations, step;
                                 rng=Random.GLOBAL_RNG,
                                 likelihood_buffer =
                                  likelihoods_buffer(get_states(grid, filter))
                                )

  # Process all observations at the current step
  #
  # Note that observations are sorted and that observations before the current
  # steps *have already been processed*
  while !all_processed(observations) && access(observations).step == step
    robot_state = access(observations).state

    range_observation = generate_observation(grid, sensor, robot_state,
                                             target_state; rng = rng)

    measurement_update!(filter, robot_state, get_states(grid, filter), sensor,
                        grid, range_observation, buffer=likelihood_buffer)

    advance!(observations)
  end

  nothing
end

# Sum mutual information objective as by Ryan and Hedrick.
# Computes:
#
# sum_i=1:t I(X_i;Y_1,...,Y_t)
#
# with observations for multiple robots. This is not the exact mutual
# information for the horizon but an upper bound.
#
# prior: prior at planning time (the first observation in the horizon,
# occurs after the first process update)
#
# sensor: sensor model (common to robots)
#
# trajectories: array of trajectories over the current horizon. Each trajectory
# is an array of states (trajectories then an array of arrays). The first state
# in each trajectory corresponds to the robot's position at the next time-step
# (after one process update on the prior)
#
const default_num_information_samples = 1000

# Package a single trajectory in an array
function finite_horizon_information(grid::Grid, prior::AbstractFilter,
                                    sensor::RangingSensor,
                                    observation::Trajectory;
                                    kwargs...
                                   )

  finite_horizon_information(grid, prior, sensor, [observation]; kwargs...)
end
# Pull the number of time-steps from a vector of trajectories as appropriate
# and drop that into the keyword arguments
function finite_horizon_information(grid::Grid, prior::AbstractFilter,
                                    sensor::RangingSensor,
                                    observations::Vector{Trajectory};
                                    prior_observations::Vector{Trajectory}=
                                      Trajectory[],
                                    kwargs...
                                   )

  if length(observations) == 0
    error("Please supply at least one trajectory for posterior")
  end

  horizon = length(observations[1])

  if horizon == 0
    error("Information horizon is zero")
  end

  for trajectories in (observations, prior_observations)
    for trajectory in trajectories
      if length(trajectory) != horizon
        @show observations
        @show prior_observations
        error("Information trajectory lengths do not match")
      end
    end
  end

  finite_horizon_information(grid, prior, sensor, observations, horizon;
                             prior_observations=prior_observations,
                             kwargs...)
end

# Base implementation of the information objective for finite horizons that is
# called by all others
#
# (Note from above that this is sum of conditional mutual informations at each
# step)
information_counter = Threads.Atomic{Int64}(0)
const samples_per_gc = 100_000
function finite_horizon_information(grid::Grid, prior::AbstractFilter,
                                    sensor::RangingSensor,
                                    posterior_observations::O,
                                    horizon::Integer;
                                    prior_observations::O = O(),
                                    num_samples::Integer =
                                      default_num_information_samples,
                                    kwargs...
                                   ) where O <: Union{Vector{TimedObservation},
                                                      Vector{Trajectory}}

  # add one to the ccounter, and add one to the old value
  old = Threads.atomic_add!(information_counter, num_samples)
  new = old + num_samples
  # Run the garbage collector if we pass a sample boundary
  if div(new, samples_per_gc) > div(old, samples_per_gc)
    GC.gc()
  end

  prior_buffer = duplicate(prior)
  reset_prior() = copy_filter!(prior, out=prior_buffer)

  # Compute entropies over the horizon, conditional on the prior trajectories
  # to compute the entropies:
  #
  # sum_i=1^t H(Xi|Y_{prior_observations}, belief)
  prior_entropies = mean(1:num_samples) do _
    reset_prior()
    sample_finite_horizon_entropy!(grid, prior_buffer, sensor,
                                   prior_observations, horizon; kwargs...)
  end

  # Compute entropies over the horizon conditional on the input trajectory
  # (produces an array with one entry per time-step) to compute entropies:
  #
  # sum_i=1^t H(Xi|Y_{all_trajectories}, belief)
  all_observations = combine_observations(prior_observations,
                                          posterior_observations)
  conditional_entropies = mean(1:num_samples) do _
    reset_prior()
    sample_finite_horizon_entropy!(grid, prior_buffer, sensor,
                                   all_observations, horizon; kwargs...)
  end

  (
   reward = sum(prior_entropies - conditional_entropies),
   mutual_information = prior_entropies - conditional_entropies,
   prior_entropies = prior_entropies,
   conditional_entropies = conditional_entropies
  )
end

# We compute the mutual information by sampling the conditional entropy given
# ranging observations.
#
# This method samples the observations and computes the entropy
#
# Returns an array of entropies with entries for each step
# *on which entropy is evaluated*
#
# The "entropy_only_at_end" keyword argument allows the designer
# to evaluate entropy at only one time-step (the last in the horizon)
#
#
# WARNING: We operate on the prior/filter in place! It is up to the user to
# maintain the original.

# Package timed observations before passing them into the main implementation
function sample_finite_horizon_entropy!(grid::Grid, prior::AbstractFilter,
                                       sensor::RangingSensor,
                                       observations::Vector{TimedObservation},
                                       b...;
                                       kwargs...)

  sample_finite_horizon_entropy!(grid, prior, sensor,
                                TimedObservations(observations), b...;
                                kwargs...)
end
function sample_finite_horizon_entropy!(grid::Grid, prior::AbstractFilter,
                                       sensor::RangingSensor,
                                       observations::Union{Vector{Trajectory},
                                                           TimedObservations},
                                       horizon::Integer;
                                       rng = Random.GLOBAL_RNG,
                                       entropy_only_at_end=false
                                      )

  target_state = sample_state(grid, prior; rng=rng)

  # Simulate trajectories and compute updates
  #
  # Sample ranging observations
  # (by also sampling sampling target trajectories)
  # Note: we will reuse samples accross the horizon

  conditional_entropies = Float64[]

  neighbor_buffer = neighbors_buffer()
  likelihood_buffer = likelihoods_buffer(get_states(grid, prior))

  for step = 1:horizon
    # Update from prior state to the first time-step in the horizon
    target_state = target_dynamics(grid, target_state, buffer=neighbor_buffer)

    process_update!(prior, transition_matrix(grid))

    # Sample observations based on the target state and each robot's state at
    # the given time-step, and then perform measurement updates
    simulate_update_filter!(grid, prior, sensor, target_state, observations,
                            step; rng=rng, likelihood_buffer=likelihood_buffer)

    # Evaluate entropy if evaluating entropy at all time-steps or if at the last
    # time-step
    if !entropy_only_at_end || step == horizon
      push!(conditional_entropies, entropy(prior))
    end
  end

  conditional_entropies
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
using Test
using LogicCircuits
using ProbabilisticCircuits
using DataFrames: DataFrame
using CUDA

include("../helper/gpu.jl")

@testset "MAP" begin
    prob_circuit = zoo_psdd("little_4var.psdd");

    data_full = generate_data_all(num_variables(prob_circuit))

    map, mappr = MAP(prob_circuit, data_full)

    @test map == data_full

    evipr = EVI(prob_circuit, data_full)
    @test mappr ≈ evipr atol=1e-6

    data_marg = DataFrame([false false false false;
                      false true true false;
                      false false true true;
                      false false false missing;
                      missing true false missing;
                      missing missing missing missing;
                      false missing missing missing])

    map, mappr = MAP(prob_circuit, data_marg)

    @test all(zip(eachcol(map), eachcol(data_marg))) do (cf,cm)
        all(zip(cf, cm)) do (f,m)
            ismissing(m) || f == m
        end
    end

    mar = MAR(prob_circuit, data_marg)

    @test all(mar .> mappr .- 1e-6)

    map, mappr = MAP(prob_circuit, false, false, false, missing)
    @test map == [false, false, false, true]
    @test mappr ≈ -1.2729657

    # same MAP states on CPU and GPU
    cpu_gpu_agree(data_full) do d
        MAP(prob_circuit, d)[1]
    end

    # same MAP probabilities on CPU and GPU
    cpu_gpu_agree_approx(data_full) do d
        MAP(prob_circuit, d)[2]
    end

    # same MAP states on CPU and GPU
    cpu_gpu_agree(data_marg) do d
        MAP(prob_circuit, d)[1]
    end

    # same MAP probabilities on CPU and GPU
    cpu_gpu_agree_approx(data_marg) do d
        MAP(prob_circuit, d)[2]
    end

end

@testset "MAP upward pass" begin
    a,b = pos_literals(ProbCircuit, 2)
    circuit = 0.6 * (a * (.5 * b + .5 * -b)) + .4 * (-a * (0.9 * b + .1 * -b))
    data_marg = DataFrame([missing missing])
    map, mappr = MAP(circuit, data_marg)
    data_marg = DataFrame([true true; true false; false true; false false])
    mar = MAR(circuit, data_marg)
    @test all(mar .> mappr .- 1e-6)
end


text:
"""
    DiscretisedFluidQueue{T <: Mesh}

A container to hold a BoundedFluidQueue and a discretisation scheme (a mesh).
Parametric type depending on the type of discretisation scheme (<:Mesh)
# Arguments:
- `model::Model`: a BoundedFluidQueue
- `mesh::T`: either a DGMesh, FVMesh or FRAPMesh specifying the discretisation scheme for the
    continuous level of the fluid queue, X(t).
"""
struct DiscretisedFluidQueue{T<:Mesh}
    model::Model
    mesh::T
    function DiscretisedFluidQueue{T}(model::Model,mesh::T) where T
        (model.b!=mesh.nodes[end])&&throw(DomainError("mesh.nodes[end] must equal model.b"))
        return new{T}(model::Model,mesh::T)
    end
end
DiscretisedFluidQueue(model::Model,mesh::T) where T =
    DiscretisedFluidQueue{T}(model,mesh)

n_phases(dq::DiscretisedFluidQueue) = n_phases(dq.model)
rates(dq::DiscretisedFluidQueue,i::Int) = rates(dq.model,i)
rates(dq::DiscretisedFluidQueue) = rates(dq.model)
phases(dq::DiscretisedFluidQueue) = phases(dq.model)
N₋(dq::DiscretisedFluidQueue) = N₋(dq.model)
N₊(dq::DiscretisedFluidQueue) = N₊(dq.model)

n_intervals(dq::DiscretisedFluidQueue) = n_intervals(dq.mesh)
Δ(dq::DiscretisedFluidQueue) = Δ(dq.mesh)
Δ(dq::DiscretisedFluidQueue,k::Int) = Δ(dq.mesh,k)
"""
    n_bases_per_cell(m) =

Return the number of basis functions on each cell of a given DiscretisedFluidQueue or Mesh.
"""
n_bases_per_cell(dq::DiscretisedFluidQueue) = n_bases_per_cell(dq.mesh)
n_bases_per_phase(dq::DiscretisedFluidQueue) = n_bases_per_phase(dq.mesh)
"""
    total_n_bases(m)

Return the number of basis functions of a given DiscretisedFluidQueue or Mesh.
"""
total_n_bases(dq::DiscretisedFluidQueue) = n_bases_per_phase(dq)*n_phases(dq)
cell_nodes(dq::DiscretisedFluidQueue) = cell_nodes(dq.mesh)

n_bases_per_level(dq::DiscretisedFluidQueue) = n_phases(dq)*n_bases_per_cell(dq)

function qbd_idx(dq::DiscretisedFluidQueue)
    ## Make QBD index
    model = dq.model
    mesh = dq.mesh

    c = N₋(dq)
    n₊ = N₊(dq)
    n₋ = N₋(dq)
    QBDidx = zeros(Int, total_n_bases(dq) + n₊ + n₋)
    for k = 1:n_intervals(dq), i = 1:n_phases(dq), n = 1:n_bases_per_cell(dq)
        c += 1
        QBDidx[c] = (i - 1) * n_bases_per_phase(dq) + (k - 1) * n_bases_per_cell(dq) + n + n₋
    end
    QBDidx[1:n₋] = 1:n₋
    QBDidx[(end-n₊+1):end] = (total_n_bases(dq) + n₋) .+ (1:n₊)

    return QBDidx
end


text:
using MAT

W = [    0    0.3314    0.3314    0.1657    0.1482    0.1482;
    0.3314         0    0.1657    0.1482    0.1657    0.1172;
    0.3314    0.1657         0    0.1482    0.1172    0.1657;
    0.1657    0.1482    0.1482         0    0.3314    0.3314;
    0.1482    0.1657    0.1172    0.3314         0    0.1657;
    0.1482    0.1172    0.1657    0.3314    0.1657         0]

function load_sample(mat_file_name::String, var_name::String = "spat_data")
    file = matopen(mat_file_name)
    var = read(file, var_name) # note that this does NOT introduce a variable ``varname`` into scope
    close(file)
    # If No var, then broken
    keywd = split(mat_file_name, "/")[end]
    start_idx = 3
    if keywd == "gaussian_diff"
        start_idx = 4
    end

    SPtrend = 3.6814e-04
    spat_data = var[:, start_idx:end]' .- SPtrend
end


text:
using PackageCompiler
create_sysimage([:TestCall, :SpecialFunctions], replace_default=true, precompile_execution_file="precompile_testcall.jl")

max_stars_repo_name: el-oso/TestCall

max_issues_repo_name: el-oso/TestCall

max_forks_repo_name: el-oso/TestCall


text:
@testset "701.insert-into-a-binary-search-tree.jl" begin
    @test insert_into_bst(TreeNode{Int}([4, 2, 7, 1, 3]), 5) ==
        TreeNode{Int}([4, 2, 7, 1, 3, 5])
    @test insert_into_bst(TreeNode{Int}([40, 20, 60, 10, 30, 50, 70]), 25) ==
        TreeNode{Int}([40, 20, 60, 10, 30, 50, 70, nothing, nothing, 25])
    @test insert_into_bst(
        TreeNode{Int}([
            4, 2, 7, 1, 3, nothing, nothing, nothing, nothing, nothing, nothing
        ]),
        5,
    ) == TreeNode{Int}([4, 2, 7, 1, 3, 5])
end


text:
"""
FinEtools (C) 2017-2019, Petr Krysl

Finite Element tools.  Julia implementation  of the finite element method
for continuum mechanics. Package for heat diffusion problems.
"""
module FinEtoolsHeatDiff

__precompile__(true)

include("allmodules.jl")

# Exports follow:

###########################################################################
# Heat diffusion functionality
###########################################################################
using FinEtoolsHeatDiff.MatHeatDiffModule: MatHeatDiff
# Exported: type of heat-diffusion  material
export MatHeatDiff

using FinEtoolsHeatDiff.FEMMHeatDiffModule: FEMMHeatDiff, conductivity, nzebcloadsconductivity, energy, inspectintegpoints, capacity
# Exported: type  for linear heat diffusion and discretization methods
export FEMMHeatDiff, conductivity, nzebcloadsconductivity, energy, inspectintegpoints, capacity

using FinEtoolsHeatDiff.FEMMHeatDiffSurfModule: FEMMHeatDiffSurf, surfacetransfer, surfacetransferloads, nzebcsurfacetransferloads
# Exported: type  for linear heat diffusion boundary conditions and discretization methods
export FEMMHeatDiffSurf, surfacetransfer, surfacetransferloads, nzebcsurfacetransferloads

end # module

max_stars_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff

max_issues_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff

max_forks_repo_name: PetrKryslUCSD/FinEtoolsHeatDiff


text:
module DiracNotation

using Requires
include("prettyprint.jl")

function __init__()
    @require QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c" include("qo.jl")
end

end # module


text:
module Matmul

using ..NeuralAttentionlib: @imexport

@imexport import ..NeuralAttentionlib: matmul,
    collapsed_size, noncollapsed_size,
    CollapsedDimArray, collapseddim, unwrap_collapse

"""
    collapseddim(x::AbstractArray, xi, xj)

Reshape `x` into 3 dim array, equivalent to `reshape(x, collapsed_size(x, xi, xj))`

See also: [`collapsed_size`](@ref)
"""
collapseddim(x::AbstractArray, xi, xj)

"""
    collapseddim(ca::CollapsedDimArray)

remove the wrapper and really reshape it.

See also: [`CollapsedDimArray`](@ref), [`unwrap_collapse`](@ref)
"""
collapseddim(ca::CollapsedDimArray)

"""
    unwrap_collapse(ca::CollapsedDimArray)

Return the underlying array of `CollapsedDimArray`, otherwise just return the input.
"""
unwrap_collapse

"""
    matmul(a::AbstractArray, b::AbstractArray, s::Number = 1)

Equivalent to `s .* (a * b)` if `a` and `b` are `Vector` or `Matrix`. For array with higher dimension,
 it will convert `a` and `b` to [`CollapsedDimArray`](@ref) and perform batched matrix multiplication, and then
 return the result as `CollapsedDimArray`. This is useful for preserving the dimensionality. If the batch dimension
 of `a` and `b` have different shape, it pick the shape of `b` for batch dimension. Work with `NNlib.batch_transpose`
 and `NNlib.batch_adjoint`.

# Example

```julia
# b-dim shape: (6,)
julia> a = CollapsedDimArray(randn(3,4,2,3,6), 3, 5); size(a)
(12, 6, 6)

# b-dim shape: (3,1,2)
julia> b = CollapsedDimArray(randn(6,2,3,1,2), 2, 3); size(b)
(6, 2, 6)

julia> c = matmul(a, b); size(c), typeof(c)
((12, 2, 6), CollapsedDimArray{Float64, Array{Float64, 6}, Static.StaticInt{3}, Static.StaticInt{4}, Static.False})

# b-dim shape: (3,1,2)
julia> d = unwrap_collapse(c); size(d), typeof(d)
((3, 4, 2, 3, 1, 2), Array{Float64, 6})

# equivanlent to `batched_mul` but preserve shape
julia> NNlib.batched_mul(collapseddim(a), collapseddim(b)) == collapseddim(matmul(a, b))
true

```

See also: [`CollapsedDimArray`](@ref), [`unwrap_collapse`](@ref), [`collapseddim`](@ref)
"""
matmul

end


text:
using Documenter, IsoRank

makedocs(
           format = :html,
           sitename = "Isorank",
           pages = [
                    "index.md", "Functions" => "funs.md"
           ]
       )

deploydocs(
           repo = "github.com/vvjn/IsoRank.jl.git",
           target = "build",
           deps   = nothing,
           make   = nothing,
           julia = "0.6"
)


text:
module aoc_10
    using Test
    export read_input, part1, part2
    export test_part1, test_part2, test
    export input, demo

    const input_r = r""

    function read_input(filename = "input.txt")
        [parse_line(line) for line in eachline(filename)]
    end

    function parse_line(line)
        parse(Int,line)
    end

    const input = read_input("input.txt")
    const demo = read_input("demo.txt")
    const demo2 = read_input("demo2.txt")

    function part1(input = input)
        diffs = diff(sort([input; 0]))
        accum = zeros(maximum(diffs))
        for d in diffs
            accum[d] += 1
        end
        accum[3] += 1
        accum
        # accum[1] * accum[3]
    end

    function test_part1()
        @testset "Part 1" begin
            @test true
        end
        part1( demo )
        #part1( demo2 )
    end

    "Calculate length of differences by only one"
    function calculate_run_lengths(input = input)
        # We need to add zero to the beginning
        # but not the device adapter since that has to remain
        sorted = sort([input; 0])
        diffs = diff(sorted)
        total = 0
        inrun = false
        runlength = 0
        runs = Vector{Int}()
        for i=1:length(diffs)
            if inrun
                if diffs[i] == 1
                    # println(sorted[i+1])
                    total += 1
                    runlength += 1
                else
                    push!(runs, runlength)
                    inrun = false
                end
            else
                if diffs[i] == 1
                    inrun = true
                    # println(sorted[i+1])
                    runlength = 1
                else
                    inrun = false
                end
            end
        end
        if inrun
            push!(runs, runlength)
        end
        runs
    end

    function run_lengths_to_combinations(runs)
        # After having calculated the runs of differences of 1,
        # we have to map them to combinations

        # Generally, the the combinations is 2^(length-1)
        runs2mul = zeros(Int,4)
        runs2mul[4] = 7 # We remove the possibility of three numbers being removed
        runs2mul[3] = 4 # The rest are multiples
        runs2mul[2] = 2
        runs2mul[1] = 1
        runs2mul[runs]
    end

    function part2(input = input)
        runs = calculate_run_lengths(input)
        combinations = run_lengths_to_combinations(runs)
        prod(combinations)
    end

    function test_part2()
        @testset "Part 2" begin
            @test part2(demo) == 8
            @test part2(demo2) == 19208
            @test part2() == 1727094849536
        end
        part2( demo2 )
    end

    function test()
        test_part1()
        test_part2()
    end
end


text:
# ATTENTION: This file was generated, do not edit by hand.
module GAP_pkg_factint

using Pkg.Artifacts
using GAP

function __init__()
    # ensure GAP will load the correct instance of this package
    path = artifact"factint"
    @debug "GAP package 'factint' path = " * path
    GAP.Globals.SetPackagePath(GapObj("factint"), GapObj(path))
end

function load()
    res = GAP.Globals.LoadPackage(GapObj("factint"))
    if res != true
        error("failed to load GAP package factint")
    end
end

end # module


text:
import Base: diff

## diff(ex, x)         f'
## what is the rest of the interface. This does:
## diff(ex, x, n)      f^(n)
## diff(ex, x, y, ...) f_{xy...}  # also diff(ex, (x,y))
## no support for diff(ex, x,n1, y,n2, ...), but can do diff(ex, (x,y), (n1, n2))

function diff(b1::SymbolicType, b2::BasicType{Val{:Symbol}})
    a = Basic()
    ret = ccall((:basic_diff, libsymengine), Int, (Ref{Basic}, Ref{Basic}, Ref{Basic}), a, b1, b2)
    return a
end

diff(b1::SymbolicType, b2::BasicType) =
    throw(ArgumentError("Second argument must be of Symbol type"))

function diff(b1::SymbolicType, b2::SymbolicType, n::Integer=1)
    n < 0 && throw(DomainError("n must be non-negative integer"))
    n==0 && return b1
    n==1 && return diff(b1, BasicType(b2))
    n > 1 && return diff(diff(b1, BasicType(b2)), BasicType(b2), n-1)
end

function diff(b1::SymbolicType, b2::SymbolicType, b3::SymbolicType)
    isa(BasicType(b3), BasicType{Val{:Integer}}) ? diff(b1, b2, N(b3)) : diff(b1, (b2, b3))
end

diff(b1::SymbolicType, b2::SymbolicType, b3::SymbolicType, b4::SymbolicType, b5...) =
    diff(b1, (b2,b3,b4,b5...))

## mixed partials
diff(ex::SymbolicType, bs::Tuple) = reduce((ex, x) -> diff(ex, x), bs, init=ex)
diff(ex::SymbolicType, bs::Tuple, ns::Tuple) =
    reduce((ex, x) -> diff(ex, x[1],x[2]), zip(bs, ns), init=ex)

diff(b1::SymbolicType, x::Union{String,Symbol}) = diff(b1, Basic(x))

"""
Series expansion to order `n` about point `x0`
"""
function series(ex::SymbolicType, x::SymbolicType, x0=0, n::Union{Integer, Basic}=6)
    (!isa(N(n), Integer) || n < 0) && throw(DomainError("n must be non-negative integer"))

    fc = subs(ex, x, x0)
    n==0 && return fc

    fp = ex
    nfact = Basic(1)
    for k in 1:n
        fp = diff(fp, x)
        nfact = k * nfact
        fc = fc + subs(fp, x, x0) * (x-x0)^k / nfact
    end

    fc
end


text:
#######################
## Log PDF functions ##
#######################

function _logpdf(::D, x::V, θ...) where {D<:Normal,V<:AbstractVector{<:Real}}
    (μ, σ) = θ
    return normlogpdf.(μ, σ, x)
end

function _logpdf(::D, x::V, w::T) where {D<:Bernoulli,V<:AbstractVector{<:Real},T<:Real}
    return bernlogpdf.(w, x)
end

function bernlogpdf(w::T, x::V) where {T<:Real,V<:Real}
    return isnan(x) ? 0.0 : (x == one(x) ? log(w) : log(1-w))
end

function catlogpdf(lw::T, x::V) where {T<:AbstractVector{<:Real},V<:Real}
    return 0 < x <= length(lw) ? (isinteger(x) ? lw[Int(x)] : (isnan(x) ? 0.0 : -Inf)) : -Inf
end

function _logpostpdf(d::Normal{T}, x::AbstractVector{<:Real}, θ::Tuple) where {T<:Real}

    (μ, ν, α, β) = θ

    dof = 2*α
    loc = μ
    scale = sqrt( (β*(1+1/ν)) / α )
    # Compute posterior predictive shifted and scaled, see: https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.t.html
    return @inbounds _logpostpdf.(Ref(d), x, dof, loc, scale)
end

function _logpostpdf(::Normal{T}, x::V, dof, loc, scale) where {T<:Real, V<:Real}
    return isnan(x) ? 0.0 : tdistlogpdf(dof, (x-loc)/scale) - log(scale)
end

function _logpostpdf(d::Categorical{T}, x::AbstractVector{<:Real}, θ::AbstractArray) where {T<:Real}
    α = [i for i in θ]
    lw = log.(α) .- log(sum(α))
    return _logpostpdf.(Ref(d), x, Ref(lw))
end

function _logpostpdf(::Categorical{T}, x::V, lw::AbstractVector) where {T<:Real,V<:Real}
    if isnan(x)
        return zero(V)
    else
        return isinteger(x) ? (length(lw) >= x > 0) ? lw[map(Int,x)] : map(V,-Inf) : map(V,-Inf)
    end
end

## Bayesian Poisson Gamma
function _logpostpdf(d::Poisson{V}, x::AbstractVector{<:Real}, θ::Tuple) where {V<:Real}
    α, β = θ
    dist = NegativeBinomial(α, β / (1+β))
    return _logpostpdf.(Ref(d), x, Ref(dist))
end

function _logpostpdf(::Poisson{V}, x::T, dist::NegativeBinomial) where {V<:Real,T<:Real}
    return isnan(x) ? zero(T) : Distributions.logpdf(dist, x)
end

## Bayesian Exponential Gamma
function _logpostpdf(d::Exponential{V}, x::AbstractVector{<:Real}, θ::Tuple) where {V<:Real}
    (λ, α) = θ
    return _logpostpdf.(Ref(d), x, α, λ)
end

function _logpostpdf(::Exponential{V}, x::T, α, β) where {V<:Real,T<:Real}
    return isnan(x) ? zero(T) : lomax(x, α, β)
end

## Bayesian Bernoulli Beta
function _logpostpdf(d::Bernoulli{V}, x::AbstractVector{<:Real}, θ::Tuple) where {V<:Real}
    (α, β) = θ
    p = α / (α+β)
    return _logpostpdf.(Ref(d), x, log(p), log(1-p))
end

function _logpostpdf(::D, x::T, lw0::V, lw1::V) where {T<:Real,V<:Real,D<:Bernoulli}
    return isnan(x) ? 0.0 : x == one(x) ? lw0 : lw1
end

### ################################# ###
### Marginal log likelihood functions ###
### ################################# ###
function _logmllh(d::Bernoulli{V}, prior::Distribution, x::AbstractVector{<:Real}, θ::Tuple) where {V<:Real}
    (α, β) = θ
    if isempty(x)
        θ_ = Distributions.rand(prior)
        return Distributions.logpdf(prior, θ_)
    end

    p = α / (α+β)
    return sum(_logpostpdf.(Ref(d), x, log(p), log(1-p)))
end

function _logmllh(d::Categorical{T}, prior::Distribution, x::AbstractVector{<:Real}, θ::AbstractArray) where {V<:Real, T<:Real}
    α = [i for i in θ]
    if isempty(x)
        θ_ = Distributions.rand(prior)
        return Distributions.logpdf(prior, θ_)
    end

    lw = log.(α) .- log(sum(α))
    return sum(_logpostpdf.(Ref(d), x, Ref(lw)))
end

function _logmllh(d::Normal{T}, prior::Distribution, x::AbstractVector{<:Real}, θ::Tuple) where {T<:Real}

    (μ, ν, α, β) = θ

    dof = 2*α
    loc = μ
    scale = sqrt( (β*(1+1/ν)) / α )

    if isempty(x)
        θ_ = rand(prior)
        return logpdf(prior, θ_...)
    end

    return @inbounds sum(_logpostpdf.(Ref(d), x, dof, loc, scale))
end

function _logmllh(d::Exponential{V}, prior::Distribution, x::AbstractVector{<:Real}, θ::Tuple) where {V<:Real}
    λ, α = θ
    if isempty(x)
        θ_ = Distributions.rand(prior)
        return Distributions.logpdf(prior, θ_...)
    end

    return @inbounds sum(_logpostpdf.(Ref(d), x, α, λ))
end

function _logmllh(d::Poisson{V}, prior::Distribution, x::AbstractVector{<:Real}, θ::Tuple) where {V<:Real}
    α, β = θ

    if isempty(x)
        θ_ = Distributions.rand(prior)
        return Distributions.logpdf(prior, θ_...)
    end

    dist = NegativeBinomial(α, 1 / (1+β))

    return @inbounds sum(_logpostpdf.(Ref(d), x, Ref(dist)))
end

function invgmllh(μ0, ν0, a0, b0, μn, νn, an, bn, n)
    r = ( lgamma(an) - lgamma(a0) ) + 1/2 * (log(ν0) - log(νn) ) + ( a0*log(2*b0) - an*log(2*bn) ) - ( n/2 * log(π))
    return r
end

function categoricalMLL(γ::AbstractVector{U}, T::AbstractVector) where {U<:Real}

    if length(γ) == 1
        lT = length(T)
        γsum = γ[1] * lT
        return (lgamma(γsum) - sum(lgamma.(ones(U, lT) * γ[1]))) + (sum(lgamma.(T .+ γ[1])) - lgamma(sum(T) + γsum))
    else
        return (lgamma(sum(γ)) - sum(lgamma.(γ))) + (sum(lgamma.(T + γ)) - lgamma(sum(T) + sum(γ)))
    end
end

function categoricalLPP(γ, T, x)
    return categoricalLPP(γ + T, x)
end

function categoricalLPP(γ::Vector, x::Int)
    return log( γ[x] ) - log( sum(γ) )
end

function dirichlet_postparam(α::AbstractVector, T::AbstractVector)
    return α + T
end

function dirichlet_postparam(α::AbstractFloat, T::AbstractVector)
    return T .+ (α/length(T))
end

max_stars_repo_name: crinders/BayesianSumProductNetworks

max_issues_repo_name: crinders/BayesianSumProductNetworks

max_forks_repo_name: crinders/BayesianSumProductNetworks


text:
module Plots

export Plot, Histogram, Histogram2, BarChart, Linear, Linear3, Image, Patch2D, Contour, Scatter, Quiver, Node, Circle, Ellipse, Command

using ..ColorMaps
using Discretizers
using StatsBase
using Distributed

const RealRange = Tuple{Real,Real}

include("ndgrid.jl")

abstract type Plot end

mutable struct Linear <: Plot
    data::AbstractMatrix{Real}
    mark
    markSize
    style
    legendentry
    onlyMarks
    errorBars
    Linear(data::AbstractMatrix{T}; mark=nothing, markSize=nothing, style=nothing, legendentry=nothing, onlyMarks=nothing, errorBars=nothing) where {T <: Real} = new(data, mark, markSize, style, legendentry, onlyMarks, errorBars)
end
Linear(x::AbstractVector{A}, y::AbstractVector{B}; kwargs...) where {A<:Real, B<:Real} = Linear(hcat(x, y)'; kwargs...)
Linear(data::AbstractVector{A}; kwargs...) where {A<:Real} = Linear(collect(1:length(data)), data; kwargs...)

mutable struct Linear3 <: Plot
    data::AbstractMatrix{Real}
    mark
    markSize
    style
    legendentry
    onlyMarks
    Linear3(data::AbstractMatrix{T}; mark=nothing, markSize=nothing, style=nothing, legendentry=nothing, onlyMarks=nothing) where {T<:Real} = new(data, mark, markSize, style, legendentry, onlyMarks)
end
Linear3(x::AbstractVector{A}, y::AbstractVector{B}, z::AbstractVector{C}; kwargs...) where {A<:Real, B<:Real, C<:Real} = Linear3(hcat(x, y, z)'; kwargs...)

const THRESHOLD_NSAMPLES_DISC_OURSELVES = 1000 # if we have more samples than this we discretize ourselves
function _construct_histogram_linear_data(
    data::Vector{Q},
    binedges::Vector{R},
    density::Bool, # If true, the bar height will be based on the probability density - otherwise directly on counts
    cumulative::Bool, # A cumulative histogram uses the sum of all previous bins and the current one as final value.
    ) where {Q<:Real,R<:Real}

    n = length(binedges)
    disc = LinearDiscretizer(binedges)
    counts = get_discretization_counts(disc, data)
    if cumulative
        cumsum!(counts, counts)
    end

    arr_x = convert(Vector{Float64}, binedges)
    arr_y = convert(Vector{Float64}, counts)
    if density
        arr_y ./= sum(counts)
        arr_y ./= binwidths(disc)
    end
    push!(arr_y, arr_y[end])

    Linear(hcat(arr_x, arr_y)', style="ybar interval,fill=blue!10, draw=blue", mark="none")
end
mutable struct Histogram <: Plot
    data::AbstractVector{Real}
    bins::Integer
    density::Bool
    cumulative::Bool
    style::AbstractString
    discretization::Symbol
    Histogram(data; bins=10, discretization=:default, density=false, cumulative=false, style="fill=blue!10") = new(data,bins,density,cumulative,style,discretization)
end

mutable struct BarChart <: Plot
    keys # symbolic x coords
    values
    style
    legendentry
    errorBars

    BarChart(keys::AbstractVector, values::AbstractVector{R}; style=nothing, legendentry=nothing, errorBars=nothing) where {R<:Real} = new(keys, values, style, legendentry, errorBars)
end
function BarChart(
    values::AbstractVector{R};
    kwargs...) where {R<:Real}

    keys = [string(i) for i in 1 : length(values)]
    return BarChart(keys, values; kwargs...)
end
function BarChart(
    values::AbstractVector{S};
    kwargs...) where {S<:AbstractString}

    dict = Dict{S,Int}()
    for v in values
        dict[v] = get(dict, v, 0) + 1
    end

    keys = S[]
    values = Int[]
    for (k,v) in dict
        push!(keys, k)
        push!(values, v)
    end

    return BarChart(keys, values; kwargs...)
end
function BarChart(
    values::AbstractVector{V},
    disc::CategoricalDiscretizer{N,D}; kwargs...) where {N,D,V}

    subkeys = collect(keys(disc.n2d))
    subvalues = encode(disc, values)
    return BarChart(keys, values; kwargs...)
end
function symbolic_x_coords(keys::AbstractVector)
    retval = ""
    for (i,k) in enumerate(keys)
        retval *= string(k)
        if i != length(keys)
            retval *= ", "
        end
    end
    return retval
end
symbolic_x_coords(p::BarChart) = symbolic_x_coords(p.keys)

mutable struct Contour <: Plot
    data::AbstractMatrix
    xbins
    ybins
    style
    contour_style
    number
    levels
    labels
    Contour(data, xbins, ybins; style=nothing, contour_style=nothing, number=nothing, levels=nothing, labels=nothing) = new(data, xbins, ybins, style, contour_style, number, levels, labels)
    function Contour(f::Function, xrange::RealRange, yrange::RealRange; xbins=40, ybins=40, style=nothing, contour_style=nothing, number=nothing, levels=nothing, labels=nothing)
        x = range(xrange[1], stop=xrange[2], length=xbins)
        y = range(yrange[1], stop=yrange[2], length=ybins)
        A = zeros(xbins, ybins)
        try
            A = Float64[f(xi, yi) for xi in x, yi in y]
        catch
            A = Float64[f([xi,yi]) for xi in x, yi in y]
        end
        new(A, x, y, style, contour_style, number, levels, labels)
    end
end

mutable struct Scatter <: Plot
    data::AbstractMatrix{Any}
    mark
    markSize
    style
    legendentry
    onlyMarks
    scatterClasses

    function Scatter(data::AbstractMatrix; mark=nothing, markSize=nothing, style=nothing, onlyMarks=true, legendentry=nothing, scatterClasses=nothing)
        new(data, mark, markSize, style, legendentry, onlyMarks, scatterClasses)
    end
end
Scatter(x::AbstractVector{A}, y::AbstractVector{B}; kwargs...) where {A<:Real, B<:Real} = Scatter(hcat(x, y)'; kwargs...)
Scatter(x::AbstractVector{A}, y::AbstractVector{B}, f::AbstractVector{C}; kwargs...) where {A<:Real, B<:Real, C<:Any} = Scatter(permutedims(hcat(x, y, f), [2,1]); kwargs...)
Scatter(x::A, y::B; kwargs...) where {A<:Real, B<:Real} = Scatter(hcat(x, y)'; kwargs...)
Scatter(x::A, y::B, f; kwargs...) where {A<:Real, B<:Real} = Scatter(hcat(x, y, f)'; kwargs...)

mutable struct Quiver <: Plot
    data::AbstractMatrix{Real}
    style
    legendentry
    Quiver(data::AbstractMatrix{T}; style=nothing, legendentry=nothing) where {T<:Real} = new(data, style, legendentry)
end
function Quiver(f::Function, xrange::RealRange, yrange::RealRange; style=nothing, legendentry=nothing, samples=15, normalize=true)
    x = range(xrange[1], stop=xrange[2], length=samples)
    y = range(yrange[1], stop=yrange[2], length=samples)
    (X, Y) = meshgrid(x, y)
    n = length(X)
    U = zeros(n)
    V = zeros(n)
    for i = 1:n
        (U[i], V[i]) = f(X[i], Y[i])
    end
    if normalize
        r = max(maximum(U),maximum(V))
        r /= min(minimum(diff(x)),minimum(diff(y)))
        U /= r
        V /= r
    end
    Quiver(X[:], Y[:], U, V, style=style, legendentry=legendentry)
end
function Quiver(
    x::Vector{A},
    y::Vector{B},
    u::Vector{C},
    v::Vector{D};
    kwargs...
    ) where {A<:Real,B<:Real,C<:Real,D<:Real}
    return Quiver(hcat(x, y, u, v)'; kwargs...)
end

mutable struct Node <: Plot
    data
    style
    x
    y
    axis # `nothing` will default to "axis cs", other options include "axis description cs", "xticklabel cs", etc.
    Node(data, x, y; style=nothing, axis=nothing) = new(data, style, x, y, axis)
end

mutable struct Circle <: Plot
    xc
    yc
    radius
    style
    Circle(xc=0,yc=0,radius=1;style=nothing) = new(xc,yc,radius,style)
end

mutable struct Ellipse <: Plot
    xc
    yc
    xradius
    yradius
    style
    Ellipse(xc=0,yc=0,xradius=1,yradius=1;style=nothing) = new(xc,yc,xradius,yradius,style)
end

mutable struct Command <: Plot
    cmd::AbstractString
    Command(cmd::AbstractString) = new(cmd)
end

global _imgid = 1

mutable struct Image <: Plot
    filename::AbstractString
    xmin::Real
    xmax::Real
    ymin::Real
    ymax::Real
    zmin::Real
    zmax::Real
    colorbar::Bool
    colormap::ColorMaps.ColorMap
    style
    function Image(
        A::AbstractMatrix{T},
        xrange::RealRange,
        yrange::RealRange;
        filename=nothing,
        colorbar=true,
        colormap=ColorMaps.GrayMap(),
        zmin=nothing,
        zmax=nothing,
        style=nothing,
        ) where {T <: Real}

        global _imgid
        if filename == nothing
            id=myid()*10000000000000+_imgid
            filename = "tmp_$(id).png"
            _imgid += 1
        end
        if zmin == nothing
            zmin = minimum(A)
        end
        if zmax == nothing
            zmax = maximum(A)
        end
        if zmin == zmax
            zmin -= 1.
            zmax += 1.
        end
        A = clamp.(A, zmin, zmax)
        A = A .- zmin
        A = A ./ (zmax - zmin)
        if isa(colormap, ColorMaps.ColorMap)
            write(colormap, A, filename)
        else
            write(ColorMaps.RGBArrayMap(colormap), A, filename)
        end
        new(filename, xrange[1], xrange[2], yrange[1], yrange[2], zmin, zmax, colorbar, colormap, style)
    end
    function Image(f::Function, xrange::RealRange, yrange::RealRange; filename=nothing, colorbar=true, colormap=ColorMaps.GrayMap(), zmin=nothing, zmax=nothing, xbins=100, ybins=100, style=nothing)
        x = range(xrange[1], stop=xrange[2], length=xbins)
        y = range(yrange[1], stop=yrange[2], length=ybins)
        (X, Y) = meshgrid(x, y)
        A = map(f, X, Y)
        A = reverse(A, dims=1)
        Image(A, xrange, yrange, filename=filename, colorbar=colorbar, colormap=colormap, zmin=zmin, zmax=zmax, style=style)
    end
end

mutable struct Patch2D <: Plots.Plot
    data::AbstractMatrix{Real} # d × m⋅n
                    # where m = 4 for rect and m = 3 for triangle (defaults to triangle)
                    #   and n = number of patches
                    #   and d = {x,y,color} or {x,y}
    style
    patch_type
    shader
    legendentry
end
Patch2D(data::AbstractMatrix; style="patch", patch_type=nothing, shader=nothing, legendentry=nothing) = Patch2D(data, style, patch_type, shader, legendentry)

function Histogram2(
    x::Vector{A},
    y::Vector{B};
    xmin=minimum(x),
    xmax=maximum(x),
    ymin=minimum(y),
    ymax=maximum(y),
    xbins=50,
    ybins=50,
    density=false,
    filename=nothing,
    colorbar=true,
    colormap=ColorMaps.GrayMap(),
    zmin=nothing,
    zmax=nothing,
    style=nothing,
    ) where {A<:Real, B<:Real}

    ex = range(xmin, stop=xmax, length=xbins+1)
    ey = range(ymin, stop=ymax, length=ybins+1)
    h = fit(StatsBase.Histogram, (y, x), (ey, ex), closed=:left)
    ex, ey, M = h.edges[1], h.edges[2], h.weights
    M = reverse(M, dims=1)
    if density
        scale =  xbins * ybins / ((xmax-xmin) * (ymax-ymin) * sum(M))
        M = M * scale
    end
    Image(M, (xmin, xmax), (ymin, ymax), filename=filename, colorbar=colorbar, colormap=colormap, zmin=zmin, zmax=zmax, style=style)
end
function Histogram2(
    x::Vector{A},
    y::Vector{B},
    edges_x::AbstractVector{C},
    edges_y::AbstractVector{C};
    density=false,
    style=nothing,
    ) where {A<:Real, B<:Real, C<:Real}

    h = fit(StatsBase.Histogram, (x, y), (edges_x, edges_y), closed=:left)
    ex, ey, M = h.edges[1], h.edges[2], h.weights
    m = length(ex)-1
    n = length(ey)-1
    scale =  m*n / ((ex[end]-ex[1]) * (ey[end]-ey[1]) * sum(M))

    patchdata = Array{Float64}(undef, 3, 4*n*m)
    patchidx = 0
    for i in 1 : m
        x₁, x₂ = ex[i], ex[i+1]
        for j in 1 : n
            y₁, y₂ = ey[j], ey[j+1]
            c = M[i,j]
            if density
                c /= (scale*(x₂-x₁)*(y₂-y₁))
            end

            patchidx += 1
            patchdata[1, patchidx] = x₁
            patchdata[2, patchidx] = y₁
            patchdata[3, patchidx] = c

            patchidx += 1
            patchdata[1, patchidx] = x₂
            patchdata[2, patchidx] = y₁
            patchdata[3, patchidx] = c

            patchidx += 1
            patchdata[1, patchidx] = x₂
            patchdata[2, patchidx] = y₂
            patchdata[3, patchidx] = c

            patchidx += 1
            patchdata[1, patchidx] = x₁
            patchdata[2, patchidx] = y₂
            patchdata[3, patchidx] = c
        end
    end

    if isa(style, Nothing)
        style = "patch"
    elseif isa(style, String)
        style = "patch" * (isempty(style) ? "" : (", "*style))
    end

    Patch2D(patchdata, style=style, patch_type="rectangle")
end

end # end plot module


text:
using PyPlot

function bc04plotfct(PDATA::PAINTER.PAINTER_Data,OIDATA::PAINTER.PAINTER_Input)
    nx=OIDATA.nx
    nw=OIDATA.nw
    wvl=OIDATA.wvl
    FOV=OIDATA.FOV
    x=PDATA.x
    w=PDATA.w
    if sum(w)==0
      w=1.
    end

    indpix=linspace(-(FOV / 2), (FOV / 2), nx)
    pos=round(Int, [1, round(Int,nx / 4), round(Int,nx / 2), round(Int,nx * 3 / 4), nx])
    X3D= squeeze(sqrt( x .* max(0, w) .* max(x,0) ),3)

    imshow( rotr90(X3D)  , origin = "lower")

    xticks(collect(pos - 1), round(Int,indpix[pos] * 1000) )
    xlabel("FOV (mas)")

    yticks(collect(pos - 1), round(Int,indpix[pos] * 1000) )
    ylabel("FOV (mas)")
end

# To change size of the simulation
# nx pixels

println("")
println(" -----------------------------------------")
println("| 1 Wavelength DEMO Beauty Constest 2004 |")
println(" -----------------------------------------")
println("")

Folder= joinpath(dirname(@__FILE__), "OIFITS_bc04")

savepath="bc04.jld"
CountPlot=25
nbitermax= 2000
PlotFct=bc04plotfct
aff=true
FOV = 0.05 * 1e-3 * 242

# pretty and working, maybe try to adjust a bit, but slowly :)
nx=256

rho_y    = .01
rho_spat = .001
rho_ps=rho_spat

lambda_spat=  .5 * 1e-4
lambda_L1 = 1e-3 * 2

rho_spec=0 # 1 wavelength -> No spectral regularization
epsilon=1e-6 # enough for the pb

## For fast result
xinit3D = PAINTER.mask( nx,round(Int, nx/3 ), choice="disk")
mask3D  = PAINTER.mask( nx,round(Int, nx/2 ))

# initialize algorithm and run admm
OIDATA, PDATA=PAINTER.painter(nbitermax=nbitermax,nx=nx,FOV=FOV,
  rho_y=rho_y, rho_spat=rho_spat,rho_spec=rho_spec,rho_ps=rho_ps,
  lambda_spat=lambda_spat, lambda_L1=lambda_L1,
  xinit3D=xinit3D,mask3D=mask3D,
  PlotFct=PlotFct,aff=aff,CountPlot=CountPlot,Folder=Folder)

# save data struture in .jld files
println("save results of BC2004 data")
PAINTER.paintersave(savepath,PDATA,OIDATA)

# load data struture in .jld files
println("load results of BC2004 data")
PDATA, OIDATA=PAINTER.painterload(savepath)

nothing


text:
using LightXML
using JuliaPBF.StructPBF

function Parsing_input_file(inputdir::String)
    # if xml
    #   return Parsing_xml_file(inpudir)
    # else if json
    #   return Parsing_json_file(inputdir)
end

# PRIVATE

function Parsing_xml_file(xmldir::String)

    # OPEN XML OBJECT
    xdoc = parse_file(xmldir) # xmldir = "./input.xml"
    xroot = root(xdoc)
    if (name(xroot) != "AnalysisSetting")
        println("CHECK THE XML")
    end

    in_timeStep         = TimeStep()
    in_EndTime          = 0.0
    in_gravity          = Vec2()
    in_analysisBox      = Box2()
    in_tensileInstability = TensileInstability(0.0,0.0,0.0)
    in_nMaxParticles    = 0
    in_maxParticleSize  = 0.0
    in_a_initialBoxes   = []
    in_a_phases         = []

    # READ XML NODES
    for depth1 in child_nodes(xroot)
    # AnalysisData

        if name(depth1) == "TimeStep"
            in_type = parse(Int64, attribute(XMLElement(depth1),"type"))
            if in_type == 1 # Fixed Time Step
                in_dt    = parse(Float64, attribute(XMLElement(depth1),"dt"))
                in_cfl   = -1
                in_dtMin = -1
                in_dtMax = -1
            elseif type == 2 # Varies Time Step
                in_dt     = -1
                in_cfl    = parse(Float64, attribute(XMLElement(depth1),"cfl"))
                in_dtMin  = parse(Float64, attribute(XMLElement(depth1),"minDt"))
                in_dtMax  = parse(Float64, attribute(XMLElement(depth1),"maxDt"))
            else
                ErrorException("Check Tims Step")
            end
            in_timeStep = TimeStep(in_type,in_dt,in_cfl,in_dtMin,in_dtMax)
        end

        if name(depth1) == "EndTime"
            in_EndTime = parse(Float64, attribute(XMLElement(depth1),"sec"))
        end

        if name(depth1) == "Gravity"
            in_gx = parse(Float64, attribute(XMLElement(depth1),"x"))
            in_gy = parse(Float64, attribute(XMLElement(depth1),"y"))
            in_gravity = Vec2(in_gx,in_gy)
        end

        if name(depth1) == "Particles"
            in_nMaxParticles   = parse(Int64, attribute(XMLElement(depth1),"nMaxParticles"))
            in_maxParticleSize = parse(Float64, attribute(XMLElement(depth1),"maxParticleSize"))
        end

        if name(depth1) == "TensileInstability"
            in_tensileInstability.k   = parse(Float64, attribute(XMLElement(depth1),"k"))
            in_tensileInstability.dq  = parse(Float64, attribute(XMLElement(depth1),"dq"))
            in_tensileInstability.n   = parse(Float64, attribute(XMLElement(depth1),"n"))
        end

        if name(depth1) == "AnalysisBox"
            in_startPoint = Vec2(parse(Float64, attribute(XMLElement(depth1),"startX")), parse(Float64, attribute(XMLElement(depth1),"startY")))
            in_endPoint   = Vec2(parse(Float64, attribute(XMLElement(depth1),"endX")), parse(Float64, attribute(XMLElement(depth1),"endY")))
            in_analysisBox = Box2(in_startPoint, in_endPoint)
        end

        if name(depth1) == "InitialBox"
            for depth2 in child_nodes(depth1)
                if name(depth2) == "initBox"
                    in_phaseID     = parse(Int64, attribute(XMLElement(depth2),"phaseID"))
                    in_startPoint  = Vec2(parse(Float64, attribute(XMLElement(depth2),"startX")), parse(Float64, attribute(XMLElement(depth2),"startY")))
                    in_endPoint    = Vec2(parse(Float64, attribute(XMLElement(depth2),"endX")),   parse(Float64, attribute(XMLElement(depth2),"endY")))
                    in_initVelocity= Vec2(parse(Float64, attribute(XMLElement(depth2),"velX")),   parse(Float64, attribute(XMLElement(depth2),"velY")))
                    in_initialBox  = Box2(in_startPoint, in_endPoint)
                    in_a_initialBoxes = push!(in_a_initialBoxes, InitialBox(in_phaseID,in_initialBox,in_initVelocity))
                end
            end
        end

        if name(depth1) == "PhaseProperties"
            for depth2 in child_nodes(depth1)
                if name(depth2) == "phase"
                    if attribute(XMLElement(depth2),"type") == "Fluid"
                        in_type = attribute(XMLElement(depth2),"type")
                        in_name = attribute(XMLElement(depth2),"name")
                        in_size = parse(Float64, attribute(XMLElement(depth2),"size"))
                        in_density = parse(Float64, attribute(XMLElement(depth2),"density"))
                        in_viscosity = parse(Float64, attribute(XMLElement(depth2),"viscosity"))
                        in_a_phases = push!(in_a_phases, PhaseData(in_type, FluidPropertyData(in_name,in_size,in_density,in_viscosity), SolidPropertyData()))
                    elseif attribute(XMLElement(depth2),"type") == "Solid"
                        in_type = attribute(XMLElement(depth2),"type")
                        in_name = attribute(XMLElement(depth2),"name")
                        in_size = parse(Float64, attribute(XMLElement(depth2),"size"))
                        in_density = parse(Float64, attribute(XMLElement(depth2),"density"))
                        in_a_phases = push!(in_a_phases, PhaseData(in_type, FluidPropertyData(), SolidPropertyData(in_name,in_size,in_density)))
                    end
                end
            end
        end

        # 작업 중
        if name(depth1) == "Boundaries"
            # for depth2 in child_nodes(depth1)
            #     if name(depth2) == "Line"
            #         in_startPoint = Vec2(parse(Float64, attribute(XMLElement(depth1),"startX")), parse(Float64, attribute(XMLElement(depth1),"startY")))
            #         in_endPoint   = Vec2(parse(Float64, attribute(XMLElement(depth1),"endX")), parse(Float64, attribute(XMLElement(depth1),"endY")))
            #         # in_BoundaryLines = Box2(in_startPoint, in_endPoint)
            #     end
            #     if name(depth2) == "Box"
            #         in_startPoint = Vec2(parse(Float64, attribute(XMLElement(depth1),"startX")), parse(Float64, attribute(XMLElement(depth1),"startY")))
            #         in_endPoint   = Vec2(parse(Float64, attribute(XMLElement(depth1),"endX")), parse(Float64, attribute(XMLElement(depth1),"endY")))
            #         # in_boundaryBox = Box2(in_startPoint, in_endPoint)
            #     end
            # end

        end

    end

    # OUTPUT
    return AnalysisDataStruct(
                        in_timeStep,               # Struct TimeStep(type, dt, cfl, dtMin, dtMax)
                        in_EndTime,                # Float64
                        in_gravity,                # Struct Gravity(x,y)
                        in_analysisBox,            # Struct Box2(staPoint, endPoint)
                        in_nMaxParticles,          # Int64
                        in_maxParticleSize,        # Float64
                        in_maxParticleSize*1.04,   # Float64 (Kernel Radius)
                        in_tensileInstability,     # TensileInstability
                        in_a_initialBoxes,         # Vector{Box2}       // Struct Box2(staPoint, endPoint)
                        in_a_phases                # Vector{PhaseData}} // Struct PhaseData(type, Fluid, Solid)
                        )

end # function parsing_xml_file

function Parsing_json_file(jsondir::String)
    # Need to develop it
end # function Parsing_json_file

max_stars_repo_name: Hipgineer/JuliaPBF

max_issues_repo_name: Hipgineer/JuliaPBF

max_forks_repo_name: Hipgineer/JuliaPBF


text:
using Janitor
using Documenter

makedocs(;
    modules=[Janitor],
    authors="Brad <bradley.lindblad@gmail.com> and contributors",
    repo="https://github.com/bradlindblad/Janitor.jl/blob/{commit}{path}#L{line}",
    sitename="Janitor.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://bradlindblad.github.io/Janitor.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/bradlindblad/Janitor.jl",
)
# using DocumenterTools
# DocumenterTools.genkeys(Janitor)


text:
# This file includes interpolant for interplation of sampled inputs.

"""
    Interpolant(tinit, tfinal, coefinit, coeffinal)

Constructs a linnear interpolant that interpolates between the poinsts `(tinit, coefinit)` and `(tfinal, coeffinal)`.
"""
mutable struct Interpolant{TMB, INB, ITP}
    timebuf::TMB
    databuf::INB
    itp::ITP
    function Interpolant(nt::Int, nd::Int)
        timebuf = Buffer(nt)
        databuf = nd == 1 ? Buffer(nt) : Buffer(nd, nt)
        itp = [interpolation(zeros(1), zeros(1)) for i in 1 : nd]
        new{typeof(timebuf), typeof(databuf), typeof(itp)}(timebuf, databuf, itp)
    end
end

show(io::IO, interpolant::Interpolant) = print(io, "Interpolant(timebuf:$(interpolant.timebuf), ",
    "databuf:$(interpolant.databuf), itp:$(interpolant.itp))")

# Callling interpolant.
getindex(interpolant::Interpolant, idx::Int) = interpolant.itp[idx]

# Update of interpolant. That is, reinterpolation.
"""
    update!(intepolant)

Updates `interpolant` using the data in `timebuf` and `databuf` of `interpolant`.
"""
update!(interpolant::Interpolant{T1, <:AbstractVector, T2}) where {T1,T2} = interpolant.itp[1] = _update!(interpolant)
update!(interpolant::Interpolant{T1, <:AbstractMatrix, T2}) where {T1,T2} = interpolant.itp = _update!(interpolant)
_update!(interpolant) = interpolation(content(interpolant.timebuf, flip=true), content(interpolant.databuf, flip=true))

interpolation(tt, uu::AbstractMatrix) = map(row -> interpolation(tt, row), eachrow(uu))
interpolation(tt, uu::AbstractVector) = CubicSplineInterpolation(getranges(tt, uu)...; extrapolation_bc=Line())

function getranges(tt, uu)
    length(tt) < 2 && return (range(tt[1], length=2, step=eps()), range(uu[1], length=2, step=eps()))
    return range(tt[1], tt[end], length=length(tt)), uu
end


text:
module GeneralizedSDistributions

import Distributions: UnivariateDistribution,
                      ContinuousUnivariateDistribution,
                      DiscreteUnivariateDistribution
import Distributions: rand, sampler, logpdf, cdf, quantile, minimum, maximum, insupport,
                      mean, var, std, modes, mode, skewness, kurtosis, entropy, mgf, cf,
                      pdf, params

import LsqFit: curve_fit, coef
import HypergeometricFunctions: _₂F₁
import QuadGK: quadgk
import OrdinaryDiffEq: ODEProblem, solve, AutoTsit5, Rosenbrock23

export GSDist

include("GSDist.jl")
include("show.jl")
include("utils.jl")

end


text:
using JSON

abstract TimberTruck

log(t::TimberTruck, a::Dict) = error("please implement `log(truck::$(typeof(t)), args::Dict)`")

function configure(t::TimberTruck; mode = nothing)
    !in(:_mode, fieldnames(t)) && return
    t._mode = mode
end

# -------

type CommonLogTruck <: TimberTruck
    out::IO

    # for use by the framework, will be
    # ignored if absent or set to nothing
    _mode

    CommonLogTruck(out::IO, mode = nothing) = new(out, mode)

    function CommonLogTruck(filename::AbstractString, mode = nothing)
        file = open(filename, "a")
        truck = new(file, mode)
        finalizer(truck, (t)->close(t.out))
        truck
    end
end

function log(truck::CommonLogTruck, l::Dict)
    println(truck.out, "$(l[:remotehost]) $(l[:rfc931]) $(l[:authuser]) $(l[:date]) \"$(l[:request])\" $(l[:status]) $(l[:bytes])")
    flush(truck.out)
end

# -------

type LumberjackTruck <: TimberTruck
    out::IO
    _mode
    opts::Dict

    function LumberjackTruck(out::IO, mode = nothing, opts = Dict())
        setup_opts(opts)
        new(out, mode, opts)
    end

    function LumberjackTruck(filename::AbstractString, mode = nothing, opts = Dict())
        file = open(filename, "a")
        setup_opts(opts)
        truck = new(file, mode, opts)
        finalizer(truck, (t)->close(t.out))
        truck
    end

    function setup_opts(opts)
        if haskey(opts, :colors)
            opts[:is_colorized] = true
        elseif (!haskey(opts, :colors) && haskey(opts, :is_colorized) && opts[:is_colorized])
            # set default colors
            opts[:colors] = Dict{ASCIIString,Symbol}("debug" => :cyan, "info" => :blue, "warn" => :yellow, "error" => :red)
        else
            opts[:is_colorized] = false
        end

        if (!haskey(opts, :uppercase))
            opts[:uppercase] = false
        end

        opts
    end
end

function log(truck::LumberjackTruck, l::Dict)
    l = copy(l)

    date_stamp = get(l, :date, nothing)
    record = date_stamp == nothing ? "" : "$date_stamp - "

    lookup = get(l, :lookup, nothing)
    if !is(lookup, nothing)
        # lookup is a StackFrame
        name, file, line = l[:lookup].func, l[:lookup].file, l[:lookup].line
        lookup_str = "$(name)@$(basename(string(file))):$(line) - "
        record = record * lookup_str
    end

    mode = l[:mode]
    if (truck.opts[:uppercase])
        l[:mode] = uppercase(l[:mode])
    end

    record = string(record, "$(l[:mode]): $(l[:msg])")

    stacktrace = get(l, :stacktrace, nothing)
    if !is(stacktrace, nothing)
        # stacktrace is a vector of StackFrames
        record = record * string(" stack:[",
            join(
                map(f->"$(f.func)@$(basename(string(f.file))):$(f.line)", stacktrace), ", "
            ), "]"
        )
    end

    delete!(l, :date)
    delete!(l, :lookup)
    delete!(l, :stacktrace)
    delete!(l, :mode)
    delete!(l, :msg)

    for (k, v) in l
        record = string(record, " $k: $(repr(v))")
    end

    if isa(truck.out, Syslog)
        # Syslog needs to be explicitly told what the error level is.
        println(truck.out, mode, record)
    elseif (truck.opts[:is_colorized])
        # check if color has been defined for key
        if (haskey(truck.opts[:colors], mode))
            print_with_color(truck.opts[:colors][mode], truck.out, string(record,"\n"))
        # if not, don't apply colors
        else
            println(truck.out, record)
        end
    else
        println(truck.out, record)
    end
    flush(truck.out)
end

# -------

type JsonTruck <: TimberTruck
    out::IO
    _mode
end

JsonTruck(out::IO) = JsonTruck(out, nothing)

function log(truck::JsonTruck, l::Dict)
    l = copy(l)

    if haskey(l, :date)
        l[:date] = string(l[:date])
    end

    if haskey(l, :lookup)
        # lookup is a StackFrame
        l[:lookup] = Dict(
            :name => l[:lookup].func, :file => basename(string(l[:lookup].file)),
            :line => l[:lookup].line
        )
    end

    if haskey(l, :stacktrace)
        # stacktrace is a vector of StackFrames
        l[:stacktrace] = map(
            f -> Dict(:name => f.func, :file => basename(string(f.file)), :line => f.line),
            l[:stacktrace]
        )
    end

    if isa(truck.out, Syslog)
        # Syslog needs to be explicitly told what the error level is.
        println(truck.out, l[:mode], json(l))
    else
        println(truck.out, json(l))
        flush(truck.out)
    end
end


text:
using ViewerGL
GL = ViewerGL

GL.VIEW([
	GL.GLFrame2,
 	GL.GLText("PLaSM", GL.COLORS[1])
]);


text:
module ExtendableSparseTest
using ExtendableSparse
using SparseArrays
using Printf

function test_assembly(;m=1000,n=1000,xnnz=5000,nsplice=1)
    println("test_assembly( m=$(m), n=$(n), xnnz=$(xnnz), nsplice=$(nsplice))")
    A=ExtendableSparseMatrix{Float64,Int64}(m,n)
    m,n=size(A)
    S=spzeros(m,n)
    for isplice=1:nsplice
        for inz=1:xnnz
            i=rand((1:m))
            j=rand((1:n))
            a=1.0+rand(Float64)
            S[i,j]+=a
            A[i,j]+=a
        end
        flush!(A)
        for j=1:n
            @assert(issorted(A.cscmatrix.rowval[A.cscmatrix.colptr[j]:A.cscmatrix.colptr[j+1]-1]))
        end
        @assert(nnz(S)==nnz(A))

        (I,J,V)=findnz(S)
        for inz=1:nnz(S)
            @assert(A[I[inz],J[inz]]==V[inz])
        end

        (I,J,V)=findnz(A)
        for inz=1:nnz(A)
            @assert(S[I[inz],J[inz]]==V[inz])
        end

    end
    return true
end

function test_transient_construction(;m=10,n=10,d=0.1) where {Tv,Ti}
    println("test_transient_construction( m=$(m), n=$(n), d=$(d))")
    csc=sprand(m,n,d)
    lnk=SparseMatrixLNK(csc)
    csc2=SparseMatrixCSC(lnk)
    return csc2==csc
end

function test_addition(;m=10,n=10,d=0.1) where {Tv,Ti}
    println("test_addition( m=$(m), n=$(n), d=$(d))")
    csc=sprand(m,n,d)
    lnk=SparseMatrixLNK(csc)
    csc2=csc+lnk
    csc2==2*csc
end

function test_timing(;n=10000,m=10000,nnz=50000)
    println("test_timing( m=$(m), n=$(n), nnz=$(nnz))")

    println("SparseMatrixCSC:")
    mat=spzeros(Float64,Int64,m,n)
    @time sprand!(mat,nnz)

    println("SparseMatrixLNK:")
    extmat=SparseMatrixLNK{Float64,Int64}(m,n)
    @time sprand!(extmat,nnz)

    println("ExtendableSparseMatrix:")
    xextmat=ExtendableSparseMatrix{Float64,Int64}(m,n)
    @time begin
        sprand!(xextmat,nnz)
        @inbounds flush!(xextmat)
    end
    true
end

function test_constructors()
    println("test_constructors()")
    m1=ExtendableSparseMatrix(10,10)
    m2=ExtendableSparseMatrix(Float16,10,10)
    m3=ExtendableSparseMatrix(Float32,Int16,10,10)

    csc=sprand(10,10,0.1)
    m4=ExtendableSparseMatrix(csc)
    sprand!(m4,6)
    flush!(m4)

    true
end

function test_operations(n)
    println("test_operations(n=$(n))")
    A=ExtendableSparseMatrix(n,n)
    sprand_sdd!(A)
    b=rand(n)
    @time x=A\b
    Ax=A*x
    b≈ Ax
end
end


text:
#!/usr/bin/env julia --project
include("../shared_source/boilerplate.jl")
include("../three_dimensional/ThreeDimensionalCompressibleNavierStokesEquations.jl")
include("sphere_helper_functions.jl")

ClimateMachine.init()

#! format: off
refVals = (
[
 [ "state",     "ρ",   1.99999999999600046319e-02,  1.00000000000004396483e+00,  5.10000000000001008083e-01,  2.91613390275248574035e-01 ],
 [ "state", "ρu[1]",  -3.13311246519756673238e-11,  3.34851852426747048856e-11,  1.50582245544617117858e-13,  4.32064975692569532007e-12 ],
 [ "state", "ρu[2]",  -3.98988708940419829624e-11,  4.12795751448340701883e-11, -1.66842288197274268178e-14,  4.34989190229689036913e-12 ],
 [ "state", "ρu[3]",  -4.16490644835465690341e-11,  4.18355107496431387673e-11,  5.18313726923951858909e-14,  4.62898170296138141882e-12 ],
 [ "state",    "ρθ",  -4.90000000000050732751e+01, -9.79999999998041548821e-01, -2.49900000000000446221e+01,  1.42890561234872102148e+01 ],
],
[
 [ "state",     "ρ",    12,    12,    12,    12 ],
 [ "state", "ρu[1]",     0,     0,     0,     0 ],
 [ "state", "ρu[2]",     0,     0,     0,     0 ],
 [ "state", "ρu[3]",     0,     0,     0,     0 ],
 [ "state",    "ρθ",    12,    12,    12,    12 ],
],
)
#! format: on

#################
# RUN THE TESTS #
#################
@testset "$(@__FILE__)" begin

    ########
    # Setup physical and numerical domains
    ########
    Ω = AtmosDomain(radius = 6e6, height = 1e5)
    grid = DiscretizedDomain(
        Ω;
        elements = (vertical = 4, horizontal = 4),
        polynomial_order = (vertical = 1, horizontal = 3),
        overintegration_order = (vertical = 1, horizontal = 1),
    )

    ########
    # Define physical parameters and parameterizations
    ########
    parameters = (
        ρₒ = 1, # reference density
        cₛ = 100.0, # sound speed
        ν = 1e-5,
        ∂θ = 0.98 / 1e5,
        α = 2e-4,
        g = 10.0,
        power = 1,
    )

    ########
    # Define timestepping parameters
    ########
    Δtᴬ = min_node_distance(grid.numerical) / parameters.cₛ * 0.25
    Δtᴰ = min_node_distance(grid.numerical)^2 / parameters.ν * 0.25
    Δt = minimum([Δtᴬ, Δtᴰ])
    start_time = 0
    end_time = 86400 * 0.5

    method = SSPRK22Heuns
    timestepper = TimeStepper(method = method, timestep = Δt)
    callbacks = (Info(), StateCheck(10))# , VTKState(iteration = 2880, filepath = "."))

    physics = FluidPhysics(;
        orientation = SphericalOrientation(),
        advection = NonLinearAdvectionTerm(),
        dissipation = ConstantViscosity{Float64}(
            μ = 0,
            ν = parameters.ν,
            κ = 0.0,
        ),
        coriolis = SphereCoriolis{Float64}(Ω = 2π / 86400),
        buoyancy = Buoyancy{Float64}(α = parameters.α, g = parameters.g),
    )

    ########
    # Define boundary conditions (west east are the ones that are enforced for a sphere)
    ########
    ρu_bcs = (bottom = Impenetrable(NoSlip()), top = Impenetrable(NoSlip()))
    ρθ_bcs = (bottom = Insulating(), top = Insulating())
    BC = (ρθ = ρθ_bcs, ρu = ρu_bcs)

    ########
    # Define initial conditions
    ########
    # Earth Spherical Representation
    # longitude: λ ∈ [-π, π), λ = 0 is the Greenwich meridian
    # latitude:  ϕ ∈ [-π/2, π/2], ϕ = 0 is the equator
    # radius:    r ∈ [Rₑ - hᵐⁱⁿ, Rₑ + hᵐᵃˣ], Rₑ = Radius of sphere; hᵐⁱⁿ, hᵐᵃˣ ≥ 0

    ρ₀(p, λ, ϕ, r) =
        (1 - p.∂θ * (r - 6e6)^p.power / p.power * 1e5^(1 - p.power)) * p.ρₒ
    ρuʳᵃᵈ(p, λ, ϕ, r) = 0.0
    ρuˡᵃᵗ(p, λ, ϕ, r) = 0.0
    ρuˡᵒⁿ(p, λ, ϕ, r) = 0.0
    ρθ₀(p, λ, ϕ, r) =
        -ρ₀(p, λ, ϕ, r) *
        p.∂θ *
        (r - 6e6)^(p.power - 1) *
        1e5^(1 - p.power) *
        (p.cₛ)^2 / (p.α * p.g)

    # Cartesian Representation (boiler plate really)
    ρ₀ᶜᵃʳᵗ(p, x...) = ρ₀(p, lon(x...), lat(x...), rad(x...))
    ρu⃗₀ᶜᵃʳᵗ(p, x...) = (
        ρuʳᵃᵈ(p, lon(x...), lat(x...), rad(x...)) * r̂(x...) +
        ρuˡᵃᵗ(p, lon(x...), lat(x...), rad(x...)) * ϕ̂(x...) +
        ρuˡᵒⁿ(p, lon(x...), lat(x...), rad(x...)) * λ̂(x...)
    )
    ρθ₀ᶜᵃʳᵗ(p, x...) = ρθ₀(p, lon(x...), lat(x...), rad(x...))

    initial_conditions = (ρ = ρ₀ᶜᵃʳᵗ, ρu = ρu⃗₀ᶜᵃʳᵗ, ρθ = ρθ₀ᶜᵃʳᵗ)

    ########
    # Create the things
    ########
    model = SpatialModel(
        balance_law = Fluid3D(),
        physics = physics,
        numerics = (flux = RoeNumericalFlux(),),
        grid = grid,
        boundary_conditions = BC,
        parameters = parameters,
    )

    simulation = Simulation(
        model = model,
        initial_conditions = initial_conditions,
        timestepper = timestepper,
        callbacks = callbacks,
        time = (; start = start_time, finish = end_time),
    )

    ########
    # Run the model
    ########
    tic = Base.time()
    evolve!(simulation, model, refDat = refVals)
    toc = Base.time()
    time = toc - tic
    println(time)

end


text:
# Pricing Engines module
# module PricingEngines
struct NullPricingEngine <: PricingEngine end

_calculate!(pe::NullPricingEngine, inst::Instrument) = error("You must set a valid pricing engine")


text:
struct SortedDict{K,V} <: AbstractDict{K,V}
    keys :: Vector{K}
    values :: Vector{V}
end

SortedDict{K,V}() where{K,V} = SortedDict(Vector{K}(), Vector{V}())
function SortedDict{K,V}(ps::Vector{Pair{K,V}}) where {K,V}
    if !issorted(ps, by=first)
        ps = sort(ps, by=first)
    end
    SortedDict{K,V}(map(first, ps), map(last, ps))
end
SortedDict(kv::Vector{Pair{K,V}}) where{K,V} = SortedDict{K,V}(kv)

SortedDict(ps::Pair...)= SortedDict(ps)
SortedDict(ps::Base.Generator) = SortedDict([ps...])
SortedDict(ps::Base.Iterators.Zip) = SortedDict([k=>v for (k,v) in ps])
SortedDict{K,V}(ps::Pair{K,V}...) where {K,V} = SortedDict{K,V}(ps)
function SortedDict{K,V}(kv) where {K,V}
    d = SortedDict{K,V}()
    sizehint!(d, length(kv))
    for (k,v) in kv
        push!(d, k=>v)
    end
    return d
end

@inline Base.keytype(d::SortedDict{K,V}) where{K,V} = K
#@inline Base.valuetype(d::SortedDict{K,V}) where{K,V} = V

Base.length(d::SortedDict) = length(d.keys)
function Base.sizehint!(d::SortedDict, sz)
    sizehint!(d.keys, sz)
    sizehint!(d.values, sz)
    d
end

function _searchsortedfirst(v::Vector, k)
    i = 1
    @inbounds while i <= length(v) && isless(v[i], k)
        i += 1
    end
    return i
end
@inline function findkey(d::SortedDict, k)
    key = convert(keytype(d), k)
    if !isequal(key, k)
        throw(ArgumentError("$k is not a valid key for type $(keytype(d))"))
    end
    i = _searchsortedfirst(d.keys, key)
    @inbounds i<=length(d.keys) && d.keys[i] == k && return i, true
    i, false
end

@inline function haskey(d::SortedDict, k)
    i, exactfound = findkey(d, k)
    exactfound
end

@inline function Base.get(d::SortedDict, k, _default)
    i, exactfound = findkey(d, k)
    if !exactfound
        return _default
    end
    d.values[i]
end

@inline function Base.getindex(d::SortedDict, k)
    i, exactfound = findkey(d, k)
    if !exactfound
        throw(KeyError(k))
    end
    d.values[i]
end

function Base.setindex!(d::SortedDict, v, k)
    i, exactfound = findkey(d, k)
    if exactfound
        d.values[i] = v
    else
        insert!(d.keys, i, convert(keytype(d), k))
        insert!(d.values, i, v)
    end
    d
end

function Base.iterate(d::SortedDict, i = 1)
    @inbounds if i > length(d)
        return nothing
    else
        return (d.keys[i] => d.values[i]), i+1
    end
end

function onlysemitokens(d::SortedDict)
    return d.keys
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps