text:
function cumul_grasp_prob(arm::Arm, ref_time::Float64, ie::InteractionEvent)

    true_ref_time = (ref_time < ie.timestamps[START]) ? ie.timestamps[START] : ref_time
    grasp_time = convert(Int64, floor(ie.timestamps[FINISH] - arm.downtime/2.0 - true_ref_time))

    return cumprob(arm.grasp_model, grasp_time)
end

# Utility of a successful pickup is just 1
pickup_util(ie::InteractionEvent) = 1.0

"""
For conveyor belt, the CG is a directed chain, so we do not need any higher level coordination.
A later arm only needs to consider the previous arm's policy tree to identify objects
that are either inevitable or potentially likely to enter its own workspace.
"""
function identify_inevitable_adversarial_entries(tree::SearchTree, workspace_objs::Set{String}, decision_idx::Int64)

    if decision_idx == -1
        return workspace_objs, ""
    end

    # Now only explore a subtree if there IS a decision node that will be taken
    # N.B - This will be the chosen option in most cases
    @assert typeof(tree.nodes[decision_idx]) == SCoBASolver.DecisionNode

    # Attempted box may fail - adversarial entry
    adversarial_entry = tree.nodes[decision_idx].task_name

    possibly_considered = Set{String}()

    # Traverse the subtree and note which objects are possibly considered
    # Then add the remainder to entries
    subtree_fringe = Deque{Int64}()
    push!(subtree_fringe, decision_idx)

    while !isempty(subtree_fringe)

        stf_top = popfirst!(subtree_fringe)

        if typeof(tree.nodes[stf_top]) == SCoBASolver.DecisionNode &&
            tree.nodes[stf_top].attempt == true
            push!(possibly_considered, tree.nodes[stf_top].task_name)
        end

        if haskey(tree.child_ids, stf_top)
            for stfc in tree.child_ids[stf_top]
                push!(subtree_fringe, stfc)
            end
        end
    end

    # Get all the objs not considered during subtree search
    not_considered = setdiff(workspace_objs, possibly_considered)

    return not_considered, adversarial_entry
end

"""
Run SCoBA on the conveyor belt domain.
"""
function scoba_conveyor!(server::ArmAllocationServer, cbelt::ConveyorBelt)

    first_arm = server.agent_ordering[1]

    if server.agent_prop_set[first_arm].obj_grabbed == ""
        # Create wrapper for success probability

        @debug first_arm
        @debug server.agent_prop_set[first_arm].workspace_objs

        success_prob_fn_first(rt, ie) = cumul_grasp_prob(server.agent_set[first_arm], rt, ie)
        generate_search_tree!(server, first_arm, server.agent_prop_set[first_arm].workspace_objs,
                             success_prob_fn_first, pickup_util,server.agent_set[first_arm].downtime)

        tree = server.agent_prop_set[first_arm].tree
        entries = Set{String}()
        adv_entry = ""
        if ~(isempty(tree))
            dec_idx = get_next_attempt_idx(tree)

            if dec_idx != -1
                dec_node = tree.nodes[dec_idx]
                @debug dec_node
                grasp_coords = get_box_point(cbelt, dec_node.task_name, dec_node.timeval)
                @assert grasp_coords[1] != Inf

                server.agent_task_allocation[first_arm] = (name=dec_node.task_name, time=dec_node.timeval)
            end

            entries, adv_entry = identify_inevitable_adversarial_entries(tree,
                                    server.agent_prop_set[first_arm].workspace_objs, dec_idx)
        end
    else
        # Need to populate entries and adv_entry
        entries = server.agent_prop_set[first_arm].workspace_objs
        adv_entry = ""
    end

    for arm_name in server.agent_ordering[2:end]

        if server.agent_prop_set[arm_name].obj_grabbed != ""
            continue
        end

        @debug arm_name
        # Objs to consider is the union of inevitable entries and workspace objs
        objs_to_consider = union(entries, server.agent_prop_set[arm_name].workspace_objs)
        @debug objs_to_consider

        # If the above yields nothing, THEN just consider adversarial_entry
        if isempty(objs_to_consider)
            if adv_entry != ""
                objs_to_consider = Set{String}()
                push!(objs_to_consider, adv_entry)
            end
        end

        if !isempty(objs_to_consider)

            success_prob_fn_rest(rt, ie) = cumul_grasp_prob(server.agent_set[arm_name], rt, ie)
            generate_search_tree!(server, arm_name, objs_to_consider,
                                 success_prob_fn_rest, pickup_util,
                                 server.agent_set[arm_name].downtime)
            tree = server.agent_prop_set[arm_name].tree
            dec_idx = get_next_attempt_idx(tree)

            if dec_idx != -1
                dec_node = tree.nodes[dec_idx]
                @debug dec_node
                grasp_coords = get_box_point(cbelt, dec_node.task_name, dec_node.timeval)
                @assert grasp_coords[1] != Inf

                server.agent_task_allocation[arm_name] = (name=dec_node.task_name, time=dec_node.timeval)
            end

            entries, adv_entry = identify_inevitable_adversarial_entries(tree,
                            server.agent_prop_set[arm_name].workspace_objs, dec_idx)
        end

        # If empty, then all relevant things would anyway be empty
    end # end for arm_name

    @debug server.agent_task_allocation
end


text:
export select_device

import MPI
using CUDA

"""
    select_device()

Select the device (GPU) corresponding to the node-local MPI rank and return its ID.

See also: [`init_global_grid`](@ref)
"""
function select_device()
	if cuda_enabled()
	    check_initialized();
	    @assert CUDA.functional(true)
	    comm_l = MPI.Comm_split_type(comm(), MPI.MPI_COMM_TYPE_SHARED, me())
		if (MPI.Comm_size(comm_l) > length(CUDA.devices())) error("More processes have been launched per node than there are GPUs available."); end
		me_l = MPI.Comm_rank(comm_l)
	    CUDA.device!(me_l)
	    return me_l
	else
		error("Cannot select a device because CUDA was detected not functional when the ImplicitGlobalGrid module was loaded.")
	end
end

_select_device() = select_device()

ext: jl


text:
using Test, LinearAlgebra, LazyArrays, StaticArrays, FillArrays
import LazyArrays: MulAdd, MemoryLayout, DenseColumnMajor, DiagonalLayout, SymTridiagonalLayout, Add, AddArray,
                    MulAddStyle, Applied, ApplyStyle, LmulStyle, Lmul, QLayout, ApplyArrayBroadcastStyle, DefaultArrayApplyStyle,
                    FlattenMulStyle, RmulStyle, Rmul, ApplyLayout, arguments
import Base.Broadcast: materialize, materialize!, broadcasted

@testset "Matrix * Vector" begin
    @testset "eltype" begin
        @test @inferred(eltype(Mul(zeros(Int,2,2), zeros(Float64,2)))) == Float64
        @test @inferred(eltype(Mul(zeros(ComplexF16,2,2),zeros(Int,2,2),zeros(Float64,2)))) == ComplexF64

        v = Mul(zeros(Int,2,2), zeros(Float64,2))
        A = Mul(zeros(Int,2,2), zeros(Float64,2,2))
        @test @inferred(axes(v)) == (@inferred(axes(v,1)),) == (Base.OneTo(2),)
        @test @inferred(size(v)) == (@inferred(size(v,1)),) == (2,)
        @test @inferred(axes(A)) == (@inferred(axes(A,1)),@inferred(axes(A,2))) == (Base.OneTo(2),Base.OneTo(2))
        @test @inferred(size(A)) == (@inferred(size(A,1)),size(A,2)) == (2,2)

        Ã = materialize(A)
        @test Ã isa Matrix{Float64}
        @test Ã == zeros(2,2)

        c = similar(v)
        fill!(c,NaN)
        @test copyto!(c,v) == zeros(2)
        fill!(c,NaN)
        c .= v
        @test c == zeros(2)

        A = randn(6,5); B = randn(5,5); C = randn(5,6)
        M = Mul(A,B,C)
        @test @inferred(eltype(M)) == Float64
        @test materialize(M) == A*B*C
    end

    @testset "gemv Float64" begin
        for A in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:),
                  view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)),
            b in (randn(5), view(randn(5),:), view(randn(5),1:5), view(randn(9),1:2:9))
            c = similar(b);

            c .= applied(*,A,b)
            @test applied(*,A,b) isa Mul{MulAddStyle}
            @test all(c .=== A*b .=== BLAS.gemv!('N', 1.0, A, b, 0.0, similar(c)))
            copyto!(c, applied(*,A,b))
            @test all(c .=== A*b .=== BLAS.gemv!('N', 1.0, A, b, 0.0, similar(c)))
            c .= @~ A*b
            @test all(c .=== A*b .=== BLAS.gemv!('N', 1.0, A, b, 0.0, similar(c)))

            b̃ = copy(b)
            copyto!(b̃, applied(*,A,b̃))
            @test all(c .=== b̃)

            M = applied(*, 2.0, A,b)
            @test M isa Applied{MulAddStyle}
            c .= M
            @test all(c .=== BLAS.gemv!('N', 2.0, A, b, 0.0, similar(c)))
            copyto!(c, M)
            @test all(c .=== BLAS.gemv!('N', 2.0, A, b, 0.0, similar(c)))
            c .= @~ 2.0*A*b
            @test all(c .=== BLAS.gemv!('N', 2.0, A, b, 0.0, similar(c)))

            c = copy(b)
            M = applied(+, applied(*,A,b), c)
            @test M isa Applied{MulAddStyle}
            copyto!(c, M)
            @test all(c .=== BLAS.gemv!('N', 1.0, A, b, 1.0, copy(b)))
            c = copy(b)
            c .= applied(+, applied(*,A,b), c)
            @test all(c .=== BLAS.gemv!('N', 1.0, A, b, 1.0, copy(b)))
            c = copy(b)
            c .= @~ A*b + c
            @test all(c .=== BLAS.gemv!('N', 1.0, A, b, 1.0, copy(b)))

            c = copy(b)
            @test (@~ A*b + 2.0*c) isa Applied{MulAddStyle}
            c .= @~ A*b + 2.0*c
            @test all(c .=== BLAS.gemv!('N', 1.0, A, b, 2.0, copy(b)))

            c = copy(b)
            @test (@~ 2.0*A*b + c) isa Applied{MulAddStyle}
            c .= @~ 2.0*A*b + c
            @test all(c .=== BLAS.gemv!('N', 2.0, A, b, 1.0, copy(b)))

            c = copy(b)
            @test ( @~ 3.0*A*b + 2.0*c) isa Applied{MulAddStyle}
            c .= @~ 3.0*A*b + 2.0*c
            @test all(c .=== BLAS.gemv!('N', 3.0, A, b, 2.0, copy(b)))

            d = similar(c)
            c = copy(b)
            d .= @~ 3.0*A*b + 2.0*c
            @test all(d .=== BLAS.gemv!('N', 3.0, A, b, 2.0, copy(b)))
        end
    end

    @testset "gemv mixed array types" begin
        (A, b, c) = (randn(5,5), randn(5), 1.0:5.0)
        d = similar(b)
        d .= @~ A*b + c
        @test all(d .=== BLAS.gemv!('N', 1.0, A, b, 1.0, Vector{Float64}(c)))

        d .= @~ A*b + 2.0*c
        @test all(d .=== BLAS.gemv!('N', 1.0, A, b, 2.0, Vector{Float64}(c)))

        d .= @~ 2.0*A*b + c
        @test all(d .=== BLAS.gemv!('N', 2.0, A, b, 1.0, Vector{Float64}(c)))

        d .= @~ 3.0*A*b + 2.0 * c
        @test all(d .=== BLAS.gemv!('N', 3.0, A, b, 2.0, Vector{Float64}(c)))

        @test (similar(b) .= @~ 1 * A * b) ≈
            (similar(b) .= @~ 1 * A *b + 0 * b) ≈
            A*b
    end

    @testset "gemv Complex" begin
        for T in (ComplexF64,),
                A in (randn(T,5,5), view(randn(T,5,5),:,:)),
                b in (randn(T,5), view(randn(T,5),:))
            c = similar(b);

            c .= @~ A*b
            @test all(c .=== BLAS.gemv!('N', one(T), A, b, zero(T), similar(c)))

            b .= @~ A*b
            @test all(c .=== b)

            c .= applied(*, 2one(T), A, b)
            @test all(c .=== BLAS.gemv!('N', 2one(T), A, b, zero(T), similar(c)))

            c = copy(b)
            c .= @~ A*b + c
            @test all(c .=== BLAS.gemv!('N', one(T), A, b, one(T), copy(b)))

            c = copy(b)
            c .= applied(+, @~(A*b), applied(*, 2one(T), c))
            @test all(c .=== BLAS.gemv!('N', one(T), A, b, 2one(T), copy(b)))

            c = copy(b)
            c .= applied(+, applied(*, 2one(T), A, b), c)
            @test all(c .=== BLAS.gemv!('N', 2one(T), A, b, one(T), copy(b)))

            c = copy(b)
            c .= applied(+, applied(*, 3one(T), A, b), applied(*, 2one(T), c))
            @test all(c .=== BLAS.gemv!('N', 3one(T), A, b, 2one(T), copy(b)))

            d = similar(c)
            c = copy(b)
            d .= applied(+, applied(*, 3one(T), A, b), applied(*, 2one(T), c))
            @test all(d .=== BLAS.gemv!('N', 3one(T), A, b, 2one(T), copy(b)))
        end
    end
end
@testset "Matrix * Matrix" begin
    @testset "gemm" begin
        for A in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:),
                  view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)),
            B in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:),
                      view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5))
            C = similar(B);

            C .= @~ A*B
            @test all(C .=== BLAS.gemm!('N', 'N', 1.0, A, B, 0.0, similar(C)))

            B .= @~ A*B
            @test all(C .=== B)

            C .= @~ 2.0 * A*B
            @test all(C .=== BLAS.gemm!('N', 'N', 2.0, A, B, 0.0, similar(C)))

            C = copy(B)
            C .= @~ A*B + C
            @test all(C .=== BLAS.gemm!('N', 'N', 1.0, A, B, 1.0, copy(B)))

            C = copy(B)
            C .= @~ A*B + 2.0 * C
            @test all(C .=== BLAS.gemm!('N', 'N', 1.0, A, B, 2.0, copy(B)))

            C = copy(B)
            C .= @~ 2.0 * A*B + C
            @test all(C .=== BLAS.gemm!('N', 'N', 2.0, A, B, 1.0, copy(B)))

            C = copy(B)
            C .= @~ 3.0 * A*B + 2.0 * C
            @test all(C .=== BLAS.gemm!('N', 'N', 3.0, A, B, 2.0, copy(B)))

            d = similar(C)
            C = copy(B)
            d .= @~ 3.0 * A*B + 2.0 * C
            @test all(d .=== BLAS.gemm!('N', 'N', 3.0, A, B, 2.0, copy(B)))
        end
    end

    @testset "gemm Complex" begin
        for T in (ComplexF64,),
                A in (randn(T,5,5), view(randn(T,5,5),:,:)),
                B in (randn(T,5,5), view(randn(T,5,5),:,:))
            C = similar(B);

            C .= @~ A*B
            @test all(C .=== BLAS.gemm!('N', 'N', one(T), A, B, zero(T), similar(C)))

            B .= @~ A*B
            @test all(C .=== B)

            C .= applied(*, 2one(T), A, B)
            @test all(C .=== BLAS.gemm!('N', 'N', 2one(T), A, B, zero(T), similar(C)))

            C = copy(B)
            C .= @~ A*B + C
            @test all(C .=== BLAS.gemm!('N', 'N', one(T), A, B, one(T), copy(B)))

            C = copy(B)
            C .= applied(+, @~(A*B), applied(*, 2one(T), C))
            @test all(C .=== BLAS.gemm!('N', 'N', one(T), A, B, 2one(T), copy(B)))

            C = copy(B)
            C .= applied(+, applied(*, 2one(T), A,B), C)
            @test all(C .=== BLAS.gemm!('N', 'N', 2one(T), A, B, one(T), copy(B)))

            C = copy(B)
            C .= applied(+, applied(*, 3one(T),A,B), applied(*, 2one(T), C))
            @test all(C .=== BLAS.gemm!('N', 'N', 3one(T), A, B, 2one(T), copy(B)))

            d = similar(C)
            C = copy(B)
            d .= applied(+, applied(*,3one(T),A,B), applied(*,2one(T), C))
            @test all(d .=== BLAS.gemm!('N', 'N', 3one(T), A, B, 2one(T), copy(B)))
        end
    end

    @testset "gemm mixed array types" begin
        (A, B, C) = (randn(5,5), randn(5,5), reshape(1.0:25.0,5,5))
        D = similar(B)
        D .= @~ A*B + C
        @test all(D .=== BLAS.gemm!('N', 'N', 1.0, A, B, 1.0, Matrix{Float64}(C)))

        D .= @~ A*B + 2.0 * C
        @test all(D .=== BLAS.gemm!('N', 'N', 1.0, A, B, 2.0, Matrix{Float64}(C)))

        D .= @~ 2.0 * A*B + C
        @test all(D .=== BLAS.gemm!('N', 'N', 2.0, A, B, 1.0, Matrix{Float64}(C)))

        D .= @~ 3.0 * A*B + 2.0 * C
        @test all(D .=== BLAS.gemm!('N', 'N', 3.0, A, B, 2.0, Matrix{Float64}(C)))

        @test (similar(B) .= 1 * A*B) ≈
            (similar(B) .= 1 * A*B + 0 * B) ≈
            A*B
    end
end

@testset "adjtrans" begin
    @testset "gemv adjtrans" begin
        for A in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:),
                  view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5)),
                    b in (randn(5), view(randn(5),:), view(randn(5),1:5), view(randn(9),1:2:9)),
                    Ac in (transpose(A), A')
            c = similar(b);

            c .= @~ Ac*b
            @test all(c .=== BLAS.gemv!('T', 1.0, A, b, 0.0, similar(c)))

            b .= @~ Ac*b
            @test all(c .=== b)

            c .= @~ 2.0 * Ac*b
            @test all(c .=== BLAS.gemv!('T', 2.0, A, b, 0.0, similar(c)))

            c = copy(b)
            c .= @~ Ac*b + c
            @test all(c .=== BLAS.gemv!('T', 1.0, A, b, 1.0, copy(b)))

            c = copy(b)
            c .= @~ Ac*b + 2.0 * c
            @test all(c .=== BLAS.gemv!('T', 1.0, A, b, 2.0, copy(b)))

            c = copy(b)
            c .= @~ 2.0 * Ac*b + c
            @test all(c .=== BLAS.gemv!('T', 2.0, A, b, 1.0, copy(b)))

            c = copy(b)
            c .= @~ 3.0 * Ac*b + 2.0 * c
            @test all(c .=== BLAS.gemv!('T', 3.0, A, b, 2.0, copy(b)))

            d = similar(c)
            c = copy(b)
            d .= @~ 3.0 * Ac*b + 2.0 * c
            @test all(d .=== BLAS.gemv!('T', 3.0, A, b, 2.0, copy(b)))
        end
    end

    @testset "gemv adjtrans mixed types" begin
        (A, b, c) = (randn(5,5), randn(5), 1.0:5.0)
        for Ac in (transpose(A), A')
            d = similar(b)
            d .= @~ Ac*b + c
            @test all(d .=== BLAS.gemv!('T', 1.0, A, b, 1.0, Vector{Float64}(c)))

            d .= @~ Ac*b + 2.0 * c
            @test all(d .=== BLAS.gemv!('T', 1.0, A, b, 2.0, Vector{Float64}(c)))

            d .= @~ 2.0 * Ac*b + c
            @test all(d .=== BLAS.gemv!('T', 2.0, A, b, 1.0, Vector{Float64}(c)))

            d .= @~ 3.0 * Ac*b + 2.0 * c
            @test all(d .=== BLAS.gemv!('T', 3.0, A, b, 2.0, Vector{Float64}(c)))
        end
    end

    @testset "gemv adjtrans Complex" begin
        for T in (ComplexF64,),
                A in (randn(T,5,5), view(randn(T,5,5),:,:)),
                b in (randn(T,5), view(randn(T,5),:)),
                (Ac,trans) in ((transpose(A),'T'), (A','C'))
            c = similar(b);

            c .= @~ Ac*b
            @test all(c .=== BLAS.gemv!(trans, one(T), A, b, zero(T), similar(c)))

            b .= @~ Ac*b
            @test all(c .=== b)

            c .= applied(*,2one(T),Ac,b)
            @test all(c .=== BLAS.gemv!(trans, 2one(T), A, b, zero(T), similar(c)))

            c = copy(b)
            c .= @~ Ac*b + c
            @test all(c .=== BLAS.gemv!(trans, one(T), A, b, one(T), copy(b)))

            c = copy(b)
            c .= applied(+, applied(*,Ac,b), applied(*,2one(T), c))
            @test all(c .=== BLAS.gemv!(trans, one(T), A, b, 2one(T), copy(b)))

            c = copy(b)
            c .= applied(+,applied(*,2one(T),Ac,b), c)
            @test all(c .=== BLAS.gemv!(trans, 2one(T), A, b, one(T), copy(b)))

            c = copy(b)
            c .= applied(+,applied(*,3one(T),Ac,b), applied(*,2one(T), c))
            @test all(c .=== BLAS.gemv!(trans, 3one(T), A, b, 2one(T), copy(b)))

            d = similar(c)
            c = copy(b)
            d .= applied(+,applied(*,3one(T),Ac,b), applied(*,2one(T), c))
            @test all(d .=== BLAS.gemv!(trans, 3one(T), A, b, 2one(T), copy(b)))
        end
    end

    @testset "gemm adjtrans" begin
        for A in (randn(5,5), view(randn(5,5),:,:)),
            B in (randn(5,5), view(randn(5,5),1:5,:))
            for Ac in (transpose(A), A')
                C = similar(B)
                C .= @~ Ac*B
                @test all(C .=== BLAS.gemm!('T', 'N', 1.0, A, B, 0.0, similar(C)))

                B .= @~ Ac*B
                @test all(C .=== B)

                C .= @~ 2.0 * Ac*B
                @test all(C .=== BLAS.gemm!('T', 'N', 2.0, A, B, 0.0, similar(C)))

                C = copy(B)
                C .= @~ Ac*B + C
                @test all(C .=== BLAS.gemm!('T', 'N', 1.0, A, B, 1.0, copy(B)))

                C = copy(B)
                C .= @~ Ac*B + 2.0 * C
                @test all(C .=== BLAS.gemm!('T', 'N', 1.0, A, B, 2.0, copy(B)))

                C = copy(B)
                C .= @~ 2.0 * Ac*B + C
                @test all(C .=== BLAS.gemm!('T', 'N', 2.0, A, B, 1.0, copy(B)))

                C = copy(B)
                C .= @~ 3.0 * Ac*B + 2.0 * C
                @test all(C .=== BLAS.gemm!('T', 'N', 3.0, A, B, 2.0, copy(B)))

                d = similar(C)
                C = copy(B)
                d .= @~ 3.0 * Ac*B + 2.0 * C
                @test all(d .=== BLAS.gemm!('T', 'N', 3.0, A, B, 2.0, copy(B)))
            end
            for Bc in (transpose(B), B')
                C = similar(B)
                C .= @~ A*Bc
                @test all(C .=== BLAS.gemm!('N', 'T', 1.0, A, B, 0.0, similar(C)))

                B .= @~ A*Bc
                @test all(C .=== B)

                C .= @~ 2.0 * A*Bc
                @test all(C .=== BLAS.gemm!('N', 'T', 2.0, A, B, 0.0, similar(C)))

                C = copy(B)
                C .= @~ A*Bc + C
                @test all(C .=== BLAS.gemm!('N', 'T', 1.0, A, B, 1.0, copy(B)))

                C = copy(B)
                C .= @~ A*Bc + 2.0 * C
                @test all(C .=== BLAS.gemm!('N', 'T', 1.0, A, B, 2.0, copy(B)))

                C = copy(B)
                C .= @~ 2.0 * A*Bc + C
                @test all(C .=== BLAS.gemm!('N', 'T', 2.0, A, B, 1.0, copy(B)))

                C = copy(B)
                C .= @~ 3.0 * A*Bc + 2.0 * C
                @test all(C .=== BLAS.gemm!('N', 'T', 3.0, A, B, 2.0, copy(B)))

                d = similar(C)
                C = copy(B)
                d .= @~ 3.0 * A*Bc + 2.0 * C
                @test all(d .=== BLAS.gemm!('N', 'T', 3.0, A, B, 2.0, copy(B)))
            end
            for Ac in (transpose(A), A'), Bc in (transpose(B), B')
                C = similar(B)
                C .= @~ Ac*Bc
                @test all(C .=== BLAS.gemm!('T', 'T', 1.0, A, B, 0.0, similar(C)))

                B .= @~ Ac*Bc
                @test all(C .=== B)

                C .= @~ 2.0 * Ac*Bc
                @test all(C .=== BLAS.gemm!('T', 'T', 2.0, A, B, 0.0, similar(C)))

                C = copy(B)
                C .= @~ Ac*Bc + C
                @test all(C .=== BLAS.gemm!('T', 'T', 1.0, A, B, 1.0, copy(B)))

                C = copy(B)
                C .= @~ Ac*Bc + 2.0 * C
                @test all(C .=== BLAS.gemm!('T', 'T', 1.0, A, B, 2.0, copy(B)))

                C = copy(B)
                C .= @~ 2.0 * Ac*Bc + C
                @test all(C .=== BLAS.gemm!('T', 'T', 2.0, A, B, 1.0, copy(B)))

                C = copy(B)
                C .= @~ 3.0 * Ac*Bc + 2.0 * C
                @test all(C .=== BLAS.gemm!('T', 'T', 3.0, A, B, 2.0, copy(B)))

                d = similar(C)
                C = copy(B)
                d .= @~ 3.0 * Ac*Bc + 2.0 * C
                @test all(d .=== BLAS.gemm!('T', 'T', 3.0, A, B, 2.0, copy(B)))
            end
        end
    end

    @testset "symv adjtrans" begin
        A = randn(100,100)
        x = randn(100)

        @test all( (similar(x) .= applied(*,Symmetric(A),x)) .===
                    (similar(x) .= applied(*,Hermitian(A),x)) .===
                    (similar(x) .= applied(*,Symmetric(A)',x)) .===
                    (similar(x) .= applied(*,Symmetric(view(A,:,:)',:L),x)) .===
                    (similar(x) .= applied(+,applied(*,1.0,Symmetric(A),x), applied(*,0.0,similar(x)))) .===
                    (similar(x) .= applied(*,view(Symmetric(A),:,:),x)) .===
                    BLAS.symv!('U', 1.0, A, x, 0.0, similar(x)) )

        y = copy(x)
        y .= applied(*,Symmetric(A), y)
        @test all( (similar(x) .= Mul(Symmetric(A),x)) .=== y)

        @test all( (similar(x) .= Mul(Symmetric(A,:L),x)) .===
                    (similar(x) .= Mul(Hermitian(A,:L),x)) .===
                    (similar(x) .= Mul(Symmetric(A,:L)',x)) .===
                    (similar(x) .= Mul(Symmetric(view(A,:,:)',:U),x)) .===
                    (fill!(similar(x),NaN) .= applied(+,applied(*,1.0,Symmetric(A,:L),x), applied(*,0.0,similar(x)))) .===
                    (fill!(similar(x),NaN) .= Mul(view(Hermitian(A,:L),:,:),x)) .===
                    BLAS.symv!('L', 1.0, A, x, 0.0, similar(x)) )
        T = ComplexF64
        A = randn(T,100,100)
        x = randn(T,100)
        @test all( (similar(x) .= Mul(Symmetric(A),x)) .===
                    (similar(x) .= Mul(transpose(Symmetric(A)),x)) .===
                    (similar(x) .= Mul(Symmetric(transpose(view(A,:,:)),:L),x)) .===
                    (similar(x) .= applied(+,applied(*,one(T),Symmetric(A),x), applied(*,zero(T),similar(x)))) .===
                    (similar(x) .= Mul(view(Symmetric(A),:,:),x)) .===
                    BLAS.symv!('U', one(T), A, x, zero(T), similar(x)) )

        @test all( (similar(x) .= Mul(Symmetric(A,:L),x)) .===
                    (similar(x) .= Mul(transpose(Symmetric(A,:L)),x)) .===
                    (similar(x) .= Mul(Symmetric(transpose(view(A,:,:)),:U),x)) .===
                    (similar(x) .= applied(+,applied(*,one(T),Symmetric(A,:L),x), applied(*,zero(T),similar(x)))) .===
                    (similar(x) .= Mul(view(Symmetric(A,:L),:,:),x)) .===
                    BLAS.symv!('L', one(T), A, x, zero(T), similar(x)) )

        @test all( (similar(x) .= Mul(Hermitian(A),x)) .===
                    (similar(x) .= Mul(Hermitian(A)',x)) .===
                    (similar(x) .= applied(+,applied(*,one(T),Hermitian(A),x), applied(*,zero(T),similar(x)))) .===
                    (similar(x) .= Mul(view(Hermitian(A),:,:),x)) .===
                    BLAS.hemv!('U', one(T), A, x, zero(T), similar(x)) )

        @test all( (similar(x) .= Mul(Hermitian(A,:L),x)) .===
                    (similar(x) .= Mul(Hermitian(A,:L)',x)) .===
                    (similar(x) .= applied(+,applied(*,one(T),Hermitian(A,:L),x), applied(*,zero(T),similar(x)))) .===
                    (similar(x) .= Mul(view(Hermitian(A,:L),:,:),x)) .===
                    BLAS.hemv!('L', one(T), A, x, zero(T), similar(x)) )

        y = copy(x)
        y .= Mul(Hermitian(A), y)
        @test all( (similar(x) .= Mul(Hermitian(A),x)) .=== y)
    end
end

@testset "Mul" begin
    @testset "Mixed types" begin
        A = randn(5,6)
        b = rand(Int,6)
        c = Array{Float64}(undef, 5)
        c .= applied(*,A,b)

        @test_throws DimensionMismatch (similar(c,3) .= Mul(A,b))
        @test_throws DimensionMismatch (c .= Mul(A,similar(b,2)))

        d = similar(c)
        mul!(d, A, b)
        @test all(c .=== d)

        copyto!(d, MulAdd(1, A, b, 0.0, d))
        @test d == copyto!(similar(d), MulAdd(1, A, b, 0.0, d)) ≈ A*b
        @test copyto!(similar(d), MulAdd(1, A, b, 1.0, d)) ≈ A*b + d

        @test all((similar(d) .= MulAdd(1, A, b, 1.0, d)) .=== copyto!(similar(d), MulAdd(1, A, b, 1.0, d)))

        B = rand(Int,6,4)
        C = Array{Float64}(undef, 5, 4)
        C .= Mul(A,B)

        @test_throws DimensionMismatch materialize!(MulAdd(1,A,B,0,similar(C,4,4)))
        @test_throws DimensionMismatch (similar(C,4,4) .= Mul(A,B))
        @test_throws DimensionMismatch (C .= Mul(A,similar(B,2,2)))

        D = similar(C)
        mul!(D, A, B)
        @test all(C .=== D)

        A = randn(Float64,20,22)
        B = randn(ComplexF64,22,24)
        C = similar(B,20,24)
        @test all((C .= Mul(A,B)  ) .=== copyto!(similar(C), MulAdd(1.0, A, B, 0.0, C)))
        @test copyto!(similar(C), MulAdd(1.0, A, B, 0.0, C)) ≈ A*B
    end

    @testset "no allocation" begin
        function blasnoalloc(c, α, A, x, β, y)
            c .= @~ A*x
            c .= @~ α*A*x
            c .= @~ A*x + y
            c .= @~ α * A*x + y
            c .= @~ A*x + β * y
            c .= @~ α * A*x + β * y
        end

        A = randn(5,5); x = randn(5); y = randn(5); c = similar(y);

        if VERSION ≥ v"1.1"
            @inferred(MulAdd(@~ A*x + y))
            @test blasnoalloc(c, 2.0, A, x, 3.0, y) === c
            @test @allocated(blasnoalloc(c, 2.0, A, x, 3.0, y)) == 0
            Ac = A'
            blasnoalloc(c, 2.0, Ac, x, 3.0, y)
            @test @allocated(blasnoalloc(c, 2.0, Ac, x, 3.0, y)) == 0
            Aa = ApplyArray(+, A, Ac)
            blasnoalloc(c, 2.0, Aa, x, 3.0, y)
			@test_broken @allocated(blasnoalloc(c, 2.0, Aa, x, 3.0, y)) == 0
		end
    end

    @testset "multi-argument mul" begin
        A = randn(5,5)
        B = apply(*,A,A,A)
        @test B isa Matrix{Float64}
        @test all(B .=== (A*A)*A)
    end

    @testset "#14" begin
        A = ones(1,1) * 1e200
        B = ones(1,1) * 1e150
        C = ones(1,1) * 1e-300

        @test apply(*, A, applied(*,B,C)) == A*(B*C)
        @test apply(*, A , applied(*,B,C), C) == A * (B*C) * C
    end

    @testset "#15" begin
        N = 10
        A = randn(N,N); B = randn(N,N); C = randn(N,N); R1 = similar(A); R2 = similar(A)
        M = Mul(A, Mul(B, C))
        @test axes(M) == (Base.OneTo(N),Base.OneTo(N))
        @test ndims(M) == ndims(typeof(M)) == 2
        @test eltype(M) == Float64
        @test all(copyto!(R1, M) .=== A*(B*C) .=== (R2 .= M))
    end

    @testset "broadcasting" begin
        A = randn(5,5); B = randn(5,5); C = randn(5,5)
        C .= NaN
        C .= @~ 1.0 * A*B + 0.0 * C
        @test C == A*B
    end

    @testset "BigFloat" begin
        A = BigFloat.(randn(5,5))
        x = BigFloat.(randn(5))
        @test A*x == apply(*,A,x) == copyto!(similar(x), applied(*,A,x))
        @test_throws UndefRefError materialize!(MulAdd(1.0,A,x,0.0,similar(x)))
    end

    @testset "Scalar * Vector" begin
        A, x =  [1 2; 3 4] , [[1,2],[3,4]]
        @test apply(*,A,x) == A*x
    end

    @testset "Complex broadcast" begin
        A = randn(5,5) .+ im*randn(5,5)
        x = randn(5) .+ im*randn(5)
        y = randn(5) .+ im*randn(5)
        z = similar(x)
        @test all((z .= @~ (2.0+0.0im)*A*x + (3.0+0.0im)*y) .=== BLAS.gemv!('N',2.0+0.0im,A,x,3.0+0.0im,y))
    end
end

@testset "Lmul/Rmul" begin
    @testset "tri Lmul" begin
        @testset "Float * Float vector" begin
            A = randn(Float64, 100, 100)
            x = randn(Float64, 100)

            L = Lmul(UpperTriangular(A),x)
            @test size(L) == (size(L,1),) == (100,)
            @test axes(L) == (axes(L,1),) == (Base.OneTo(100),)
            @test eltype(L) == Float64
            @test length(L) == 100

            @test similar(L) isa Vector{Float64}
            @test similar(L,Int) isa Vector{Int}

            @test applied(*, UpperTriangular(A), x) isa Applied{LmulStyle}
            @test similar(applied(*, UpperTriangular(A), x), Float64) isa Vector{Float64}

            @test ApplyStyle(*, typeof(UpperTriangular(A)), typeof(x)) isa LmulStyle

            @test all((y = copy(x); y .= Mul(UpperTriangular(A),y) ) .===
                        (similar(x) .= Mul(UpperTriangular(A),x)) .===
                        copy(Lmul(UpperTriangular(A),x)) .===
                        materialize!(Lmul(UpperTriangular(A),copy(x))) .===
                        copyto!(similar(x),Lmul(UpperTriangular(A),x)) .===
                        UpperTriangular(A)*x .===
                        BLAS.trmv!('U', 'N', 'N', A, copy(x)))
            @test all((y = copy(x); y .= Mul(UnitUpperTriangular(A),y) ) .===
                        (similar(x) .= Mul(UnitUpperTriangular(A),x)) .===
                        UnitUpperTriangular(A)*x .===
                        BLAS.trmv!('U', 'N', 'U', A, copy(x)))
            @test all((y = copy(x); y .= Mul(LowerTriangular(A),y) ) .===
                        (similar(x) .= Mul(LowerTriangular(A),x)) .===
                        LowerTriangular(A)*x .===
                        BLAS.trmv!('L', 'N', 'N', A, copy(x)))
            @test all((y = copy(x); y .= Mul(UnitLowerTriangular(A),y) ) .===
                        (similar(x) .= Mul(UnitLowerTriangular(A),x)) .===
                        UnitLowerTriangular(A)x .===
                        BLAS.trmv!('L', 'N', 'U', A, copy(x)))

            @test all((y = copy(x); y .= Mul(UpperTriangular(A)',y) ) .===
                        (similar(x) .= Mul(UpperTriangular(A)',x)) .===
                        (similar(x) .= Mul(LowerTriangular(A'),x)) .===
                        UpperTriangular(A)'x .===
                        BLAS.trmv!('U', 'T', 'N', A, copy(x)))
            @test all((y = copy(x); y .= Mul(UnitUpperTriangular(A)',y) ) .===
                        (similar(x) .= Mul(UnitUpperTriangular(A)',x)) .===
                        (similar(x) .= Mul(UnitLowerTriangular(A'),x)) .===
                        UnitUpperTriangular(A)'*x .===
                        BLAS.trmv!('U', 'T', 'U', A, copy(x)))
            @test all((y = copy(x); y .= Mul(LowerTriangular(A)',y) ) .===
                        (similar(x) .= Mul(LowerTriangular(A)',x)) .===
                        (similar(x) .= Mul(UpperTriangular(A'),x)) .===
                        LowerTriangular(A)'*x .===
                        BLAS.trmv!('L', 'T', 'N', A, copy(x)))
            @test all((y = copy(x); y .= Mul(UnitLowerTriangular(A)',y) ) .===
                        (similar(x) .= Mul(UnitLowerTriangular(A)',x)) .===
                        (similar(x) .= Mul(UnitUpperTriangular(A'),x)) .===
                        UnitLowerTriangular(A)'*x .===
                        BLAS.trmv!('L', 'T', 'U', A, copy(x)))
        end

        @testset "Float * Complex vector"  begin
            T = ComplexF64
            A = randn(T, 100, 100)
            x = randn(T, 100)

            @test all((y = copy(x); y .= Mul(UpperTriangular(A),y) ) .===
                        (similar(x) .= Mul(UpperTriangular(A),x)) .===
                        UpperTriangular(A)x .===
                        BLAS.trmv!('U', 'N', 'N', A, copy(x)))
            @test all((y = copy(x); y .= Mul(UnitUpperTriangular(A),y) ) .===
                        (similar(x) .= Mul(UnitUpperTriangular(A),x)) .===
                        UnitUpperTriangular(A)x .===
                        BLAS.trmv!('U', 'N', 'U', A, copy(x)))
            @test all((y = copy(x); y .= Mul(LowerTriangular(A),y) ) .===
                        (similar(x) .= Mul(LowerTriangular(A),x)) .===
                        LowerTriangular(A)x .===
                        BLAS.trmv!('L', 'N', 'N', A, copy(x)))
            @test all((y = copy(x); y .= Mul(UnitLowerTriangular(A),y) ) .===
                        (similar(x) .= Mul(UnitLowerTriangular(A),x)) .===
                        UnitLowerTriangular(A)x .===
                        BLAS.trmv!('L', 'N', 'U', A, copy(x)))
            LowerTriangular(A')  == UpperTriangular(A)'

            @test all((y = copy(x); y .= Mul(transpose(UpperTriangular(A)),y) ) .===
                        (similar(x) .= Mul(transpose(UpperTriangular(A)),x)) .===
                        (similar(x) .= Mul(LowerTriangular(transpose(A)),x)) .===
                        transpose(UpperTriangular(A))x .===
                        BLAS.trmv!('U', 'T', 'N', A, copy(x)))
            @test all((y = copy(x); y .= Mul(transpose(UnitUpperTriangular(A)),y) ) .===
                        (similar(x) .= Mul(transpose(UnitUpperTriangular(A)),x)) .===
                        (similar(x) .= Mul(UnitLowerTriangular(transpose(A)),x)) .===
                        transpose(UnitUpperTriangular(A))x .===
                        BLAS.trmv!('U', 'T', 'U', A, copy(x)))
            @test all((y = copy(x); y .= Mul(transpose(LowerTriangular(A)),y) ) .===
                        (similar(x) .= Mul(transpose(LowerTriangular(A)),x)) .===
                        (similar(x) .= Mul(UpperTriangular(transpose(A)),x)) .===
                        transpose(LowerTriangular(A))x .===
                        BLAS.trmv!('L', 'T', 'N', A, copy(x)))
            @test all((y = copy(x); y .= Mul(transpose(UnitLowerTriangular(A)),y) ) .===
                        (similar(x) .= Mul(transpose(UnitLowerTriangular(A)),x)) .===
                        (similar(x) .= Mul(UnitUpperTriangular(transpose(A)),x)) .===
                        transpose(UnitLowerTriangular(A))x .===
                        BLAS.trmv!('L', 'T', 'U', A, copy(x)))

            @test all((y = copy(x); y .= Mul(UpperTriangular(A)',y) ) .===
                        (similar(x) .= Mul(UpperTriangular(A)',x)) .===
                        (similar(x) .= Mul(LowerTriangular(A'),x)) .===
                        UpperTriangular(A)'x .===
                        BLAS.trmv!('U', 'C', 'N', A, copy(x)))
            @test all((y = copy(x); y .= Mul(UnitUpperTriangular(A)',y) ) .===
                        (similar(x) .= Mul(UnitUpperTriangular(A)',x)) .===
                        (similar(x) .= Mul(UnitLowerTriangular(A'),x)) .===
                        UnitUpperTriangular(A)'x .===
                        BLAS.trmv!('U', 'C', 'U', A, copy(x)))
            @test all((y = copy(x); y .= Mul(LowerTriangular(A)',y) ) .===
                        (similar(x) .= Mul(LowerTriangular(A)',x)) .===
                        (similar(x) .= Mul(UpperTriangular(A'),x)) .===
                        LowerTriangular(A)'x .===
                        BLAS.trmv!('L', 'C', 'N', A, copy(x)))
            @test all((y = copy(x); y .= Mul(UnitLowerTriangular(A)',y) ) .===
                        (similar(x) .= Mul(UnitLowerTriangular(A)',x)) .===
                        (similar(x) .= Mul(UnitUpperTriangular(A'),x)) .===
                        UnitLowerTriangular(A)'x .===
                        BLAS.trmv!('L', 'C', 'U', A, copy(x)))
        end

        @testset "Float * Float Matrix" begin
            A = randn(Float64, 100, 100)
            x = randn(Float64, 100, 100)

            @test UpperTriangular(A)*x ≈ (similar(x) .= Mul(UpperTriangular(A), x))
        end

        @testset "adjtrans" begin
            for T in (Float64, ComplexF64)
                A = randn(T,100,100)
                b = randn(T,100)

                @test all(apply(*, UpperTriangular(A)', b) .=== UpperTriangular(A)'b)
                @test all(apply(*, UnitUpperTriangular(A)', b) .=== UnitUpperTriangular(A)'b)
                @test all(apply(*, LowerTriangular(A)', b) .=== LowerTriangular(A)'b)
                @test all(apply(*, UnitLowerTriangular(A)', b) .=== UnitLowerTriangular(A)'b)

                @test all(apply(*, transpose(UpperTriangular(A)), b) .=== transpose(UpperTriangular(A))b)
                @test all(apply(*, transpose(UnitUpperTriangular(A)), b) .=== transpose(UnitUpperTriangular(A))b)
                @test all(apply(*, transpose(LowerTriangular(A)), b) .=== transpose(LowerTriangular(A))b)
                @test all(apply(*, transpose(UnitLowerTriangular(A)), b) .=== transpose(UnitLowerTriangular(A))b)

                B = randn(T,100,100)

                @test all(apply(*, UpperTriangular(A)', B) .=== UpperTriangular(A)'B)
                @test all(apply(*, UnitUpperTriangular(A)', B) .=== UnitUpperTriangular(A)'B)
                @test all(apply(*, LowerTriangular(A)', B) .=== LowerTriangular(A)'B)
                @test all(apply(*, UnitLowerTriangular(A)', B) .=== UnitLowerTriangular(A)'B)

                @test all(apply(*, transpose(UpperTriangular(A)), B) .=== transpose(UpperTriangular(A))B)
                @test all(apply(*, transpose(UnitUpperTriangular(A)), B) .=== transpose(UnitUpperTriangular(A))B)
                @test all(apply(*, transpose(LowerTriangular(A)), B) .=== transpose(LowerTriangular(A))B)
                @test all(apply(*, transpose(UnitLowerTriangular(A)), B) .=== transpose(UnitLowerTriangular(A))B)
            end

            for T in (Float64, ComplexF64)
                A = big.(randn(T,100,100))
                b = big.(randn(T,100))

                @test all(apply(*, UpperTriangular(A)', b) ≈ UpperTriangular(A)'b)
                @test all(apply(*, UnitUpperTriangular(A)', b) ≈ UnitUpperTriangular(A)'b)
                @test all(apply(*, LowerTriangular(A)', b) ≈ LowerTriangular(A)'b)
                @test all(apply(*, UnitLowerTriangular(A)', b) ≈ UnitLowerTriangular(A)'b)

                @test all(apply(*, transpose(UpperTriangular(A)), b) ≈ transpose(UpperTriangular(A))b)
                @test all(apply(*, transpose(UnitUpperTriangular(A)), b) ≈ transpose(UnitUpperTriangular(A))b)
                @test all(apply(*, transpose(LowerTriangular(A)), b) ≈ transpose(LowerTriangular(A))b)
                @test all(apply(*, transpose(UnitLowerTriangular(A)), b) ≈ transpose(UnitLowerTriangular(A))b)

                B = big.(randn(T,100,100))

                @test all(apply(*, UpperTriangular(A)', B) ≈ UpperTriangular(A)'B)
                @test all(apply(*, UnitUpperTriangular(A)', B) ≈ UnitUpperTriangular(A)'B)
                @test all(apply(*, LowerTriangular(A)', B) ≈ LowerTriangular(A)'B)
                @test all(apply(*, UnitLowerTriangular(A)', B) ≈ UnitLowerTriangular(A)'B)

                @test all(apply(*, transpose(UpperTriangular(A)), B) ≈ transpose(UpperTriangular(A))B)
                @test all(apply(*, transpose(UnitUpperTriangular(A)), B) ≈ transpose(UnitUpperTriangular(A))B)
                @test all(apply(*, transpose(LowerTriangular(A)), B) ≈ transpose(LowerTriangular(A))B)
                @test all(apply(*, transpose(UnitLowerTriangular(A)), B) ≈ transpose(UnitLowerTriangular(A))B)
            end
        end
    end

    @testset "tri Rmul" begin
        for T in (Float64, ComplexF64)
            A = randn(T, 100,100)
            B = randn(T, 100,100)
            R = Rmul(copy(A), UpperTriangular(B))
            @test size(R) == (size(R,1),size(R,2)) == (100,100)
            @test axes(R) == (axes(R,1),axes(R,2)) == (Base.OneTo(100),Base.OneTo(100))
            @test eltype(R) == T
            @test length(R) == 100^2

            @test similar(R) isa Matrix{T}
            @test similar(R,Int) isa Matrix{Int}

            @test applied(*, A, UpperTriangular(B)) isa Applied{RmulStyle}
            @test similar(applied(*, A, UpperTriangular(B)), Float64) isa Matrix{Float64}

            R2 = deepcopy(R)
            Ap = applied(*, copy(A), UpperTriangular(B))
            Ap2 = applied(*, copy(A), UpperTriangular(B))
            @test all(BLAS.trmm('R', 'U', 'N', 'N', one(T), B, A) .=== apply(*, A, UpperTriangular(B)) .===
                    copyto!(similar(Ap),Ap) .=== materialize!(Ap2) .=== copyto!(similar(R2), R2) .=== materialize!(R))
            @test R.A ≠ A
            @test all(BLAS.trmm('R', 'U', 'T', 'N', one(T), B, A) .=== apply(*, A, transpose(UpperTriangular(B))) .=== A*transpose(UpperTriangular(B)))
            @test all(BLAS.trmm('R', 'U', 'N', 'U', one(T), B, A) .=== apply(*, A, UnitUpperTriangular(B)) .=== A*UnitUpperTriangular(B))
            @test all(BLAS.trmm('R', 'U', 'T', 'U', one(T), B, A) .=== apply(*, A, transpose(UnitUpperTriangular(B))) .=== A*transpose(UnitUpperTriangular(B)))
            @test all(BLAS.trmm('R', 'L', 'N', 'N', one(T), B, A) .=== apply(*, A, LowerTriangular(B)) .=== A*LowerTriangular(B))
            @test all(BLAS.trmm('R', 'L', 'T', 'N', one(T), B, A) .=== apply(*, A, transpose(LowerTriangular(B))) .=== A*transpose(LowerTriangular(B)))
            @test all(BLAS.trmm('R', 'L', 'N', 'U', one(T), B, A) .=== apply(*, A, UnitLowerTriangular(B)) .=== A*UnitLowerTriangular(B))
            @test all(BLAS.trmm('R', 'L', 'T', 'U', one(T), B, A) .=== apply(*, A, transpose(UnitLowerTriangular(B))) .=== A*transpose(UnitLowerTriangular(B)))

            if T <: Complex
                @test all(BLAS.trmm('R', 'U', 'C', 'N', one(T), B, A) .=== apply(*, A, UpperTriangular(B)') .=== A*UpperTriangular(B)')
                @test all(BLAS.trmm('R', 'U', 'C', 'U', one(T), B, A) .=== apply(*, A, UnitUpperTriangular(B)') .=== A*UnitUpperTriangular(B)')
                @test all(BLAS.trmm('R', 'L', 'C', 'N', one(T), B, A) .=== apply(*, A, LowerTriangular(B)') .=== A*LowerTriangular(B)')
                @test all(BLAS.trmm('R', 'L', 'C', 'U', one(T), B, A) .=== apply(*, A, UnitLowerTriangular(B)') .=== A*UnitLowerTriangular(B)')
            end
        end

        T = Float64
        A = big.(randn(100,100))
        B = big.(randn(100,100))
        materialize!(Rmul(A,UpperTriangular(B)))
    end

    @testset "Diagonal and SymTridiagonal" begin
        A = randn(5,5)
        B = Diagonal(randn(5))
        @test MemoryLayout(typeof(B)) == DiagonalLayout{DenseColumnMajor}()
        @test ApplyStyle(*, typeof(A), typeof(B)) == RmulStyle()
        @test apply(*,A,B) == A*B == materialize!(Rmul(copy(A),B))

        @test ApplyStyle(*, typeof(B), typeof(A)) == LmulStyle()
        @test apply(*,B,A) == B*A

        @test ApplyStyle(*, typeof(B), typeof(B)) == LmulStyle()
        @test apply(*,B,B) == B*B
        @test apply(*,B,B) isa Diagonal

        A = randn(5,5)
        B = SymTridiagonal(randn(5),randn(4))
        @test MemoryLayout(typeof(B)) == SymTridiagonalLayout{DenseColumnMajor}()
        @test apply(*,A,B) == A*B
    end
end

@testset "Factorizations" begin
    @testset "QR" begin
        A = randn(5,3)
        b = randn(3)
        B = randn(3,3)
        Q,R = qr(A)
        @test MemoryLayout(typeof(Q)) == QLayout()
        @test all(Q*b .=== apply(*,Q,b))
        @test all(Q*B .=== copyto!(similar(B,5,3),Lmul(Q,B)) .=== copyto!(similar(B,5,3),applied(*,Q,B)) .=== apply(*,Q,B))
        @test all(Q*B*b .=== apply(*,Q,B,b))
        @test_throws DimensionMismatch apply(*, Q, randn(4))
        @test_throws DimensionMismatch apply(*, Q, randn(4,3))
        dest = fill(NaN,5)
        @test copyto!(dest, applied(*,Q,b)) == Q*b

        b = randn(5)
        B = randn(5,5)
        @test all(Q*b .=== apply(*,Q,b))
        @test all(Q*b .=== apply(*,Q,b))
        @test all(Q*B .=== apply(*,Q,B))
        @test all(Q*B*b .=== apply(*,Q,B,b))

        @test MemoryLayout(typeof(Q')) == QLayout()
        @test all(Q'b .=== apply(*,Q',b))
        @test all(Q'B .=== apply(*,Q',B))
        @test all(Q'B*b .=== apply(*,Q',B,b))
        @test_throws DimensionMismatch apply(*,Q',randn(3))
        @test_throws DimensionMismatch apply(*,Q',randn(3,3))
        dest = fill(NaN,5)
        @test all(copyto!(dest, applied(*,Q',b)) .=== Q'*b)
    end

    @testset "applied axes" begin
        T = Float64
        A = randn(5,5)
        x = randn(5)
        M =  applied(+, applied(*,T(1.0),A,x), applied(*,T(0.0),x))
        @test axes(M) == axes(M.args[1]) == (Base.OneTo(5),)
    end

    @testset "Broadcast" begin
        A = randn(5,5)
        b = randn(5)
        M = Mul(A,b)
        @test Base.BroadcastStyle(typeof(M)) isa ApplyArrayBroadcastStyle{1}
        @test M .+ 1 ≈ A*b .+ 1

        @test Base.BroadcastStyle(ApplyArrayBroadcastStyle{1}(), Broadcast.DefaultArrayStyle{1}()) == Broadcast.DefaultArrayStyle{1}()
        @test Base.BroadcastStyle(ApplyArrayBroadcastStyle{1}(), Broadcast.DefaultArrayStyle{2}()) == Broadcast.DefaultArrayStyle{2}()
    end

    @testset "Diagonal" begin
       @test Diagonal(Fill(2,10))  * Fill(3,10) ≡ Fill(6,10)
       @test apply(*, Diagonal(Fill(2,10)), Fill(3,10)) ≡ Fill(6,10)
       @test_broken Diagonal(Fill(2,10))  * Fill(3,10,3) ≡ Fill(6,10)
       @test apply(*, Diagonal(Fill(2,10)), Fill(3,10,3)) ≡ Fill(6,10,3)
       @test apply(*,Fill(3,10,10),Fill(3,10)) ≡ Fill(9,10)
       @test apply(*, Eye(10), Ones(10)) == Ones(10)
       @test apply(*, Eye(10), Eye(10)) == Eye(10)
    end

    @testset "ApplyArray MulTest" begin
        A = ApplyArray(*,randn(2,2), randn(2,2))
        @test ApplyStyle(*,typeof(A),typeof(randn(2,2))) == FlattenMulStyle()
        @test ApplyArray(*,Diagonal(Fill(2,10)), Fill(3,10,10))*Fill(3,10) ≡ Fill(18,10)
        @test ApplyArray(*,Diagonal(Fill(2,10)), Fill(3,10,10))*ApplyArray(*,Diagonal(Fill(2,10)), Fill(3,10,10)) == Fill(36,10,10)
    end
end

@testset "MulAdd" begin
    A = randn(5,5)
    B = randn(5,4)
    C = randn(5,4)
    b = randn(5)
    c = randn(5)

    M = MulAdd(2.0,A,B,3.0,C)
    @test size(M) == size(C)
    @test size(M,1) == size(C,1)
    @test size(M,2) == size(C,2)
    @test_broken size(M,3) == size(C,3)
    @test length(M) == length(C)
    @test axes(M) == axes(C)
    @test eltype(M) == Float64
    @test materialize(M) ≈ 2.0A*B + 3.0C

    @test_throws DimensionMismatch materialize(MulAdd(2.0,A,randn(3),1.0,B))
    @test_throws DimensionMismatch materialize(MulAdd([1,2],A,B,[1,2],C))
    @test_throws DimensionMismatch materialize(MulAdd(2.0,A,B,3.0,randn(3,4)))
    @test_throws DimensionMismatch materialize(MulAdd(2.0,A,B,3.0,randn(5,5)))

    B = randn(5,5)
    C = randn(5,5)
    @test materialize(MulAdd(2.0,Diagonal(A),Diagonal(B),3.0,Diagonal(C))) ==
          materialize!(MulAdd(2.0,Diagonal(A),Diagonal(B),3.0,Diagonal(copy(C)))) == 2.0Diagonal(A)*Diagonal(B) + 3.0*Diagonal(C)
    @test_broken materialize(MulAdd(2.0,Diagonal(A),Diagonal(B),3.0,Diagonal(C))) isa Diagonal

    @test materialize(MulAdd(1.0, Eye(5), A, 3.0, C)) == materialize!(MulAdd(1.0, Eye(5), A, 3.0, copy(C))) == A + 3.0C
    @test materialize(MulAdd(1.0, A, Eye(5), 3.0, C)) == materialize!(MulAdd(1.0, A, Eye(5), 3.0, copy(C))) == A + 3.0C
end

@testset "MulArray" begin
  @testset "Basic" begin
        A = randn(5,5)
        M = ApplyArray(*,A,A)
        @test M[5] == M[5,1]
        @test M[6] == M[1,2]
        @test Matrix(M) ≈ A^2
        x = randn(5)
        @test x'M ≈ transpose(x)*M ≈ x'Matrix(M)
        @test_throws DimensionMismatch materialize(applied(*, randn(5,5), randn(4)))
        @test_throws DimensionMismatch ApplyArray(*, randn(5,5), randn(4))
    end

    @testset "Bug in getindex" begin
        M = ApplyArray(*,[1,2,3],Ones(1,20))
        @test M[1,1] == 1
        @test M[2,1] == 2
        M = Applied(*,[1 2; 3 4], [1 2; 3 4])
        @test M[1] == 7
    end

    @testset "Views" begin
        A = randn(500,500)
        b = randn(500)
        M = ApplyArray(*,A,b)

        V = view(M,2:300)
        @test MemoryLayout(typeof(V)) isa ApplyLayout{typeof(*)}
        @test arguments(V) == (view(A,2:300,Base.OneTo(500)),view(b, Base.OneTo(500)))
        @test Applied(V) isa Applied{MulAddStyle}
        @test ApplyArray(V) ≈ (A*b)[2:300]
        c = similar(V)
        copyto!(c,Applied(V))
        VERSION ≥ v"1.2" && @test @allocated(copyto!(c,Applied(V))) ≤ 200
        copyto!(c, V)
        VERSION ≥ v"1.2" && @test @allocated(copyto!(c, V)) ≤ 200
        @test all(c .=== apply(*, arguments(V)...))

        B = randn(500,500)
        M = ApplyArray(*,A,B)
        V = view(M,2:300,3:400)
        @test MemoryLayout(typeof(V)) isa ApplyLayout{typeof(*)}
        @test arguments(V) == (view(A,2:300,Base.OneTo(500)),view(B, Base.OneTo(500),3:400))
        @test Applied(V) isa Applied{MulAddStyle}
        c = similar(V)
        copyto!(c,Applied(V))
        VERSION ≥ v"1.2" && @test @allocated(copyto!(c,Applied(V))) ≤ 1000
        copyto!(c, V)
        VERSION ≥ v"1.2" && @test @allocated(copyto!(c, V)) ≤ 1000
        @test all(c .=== apply(*, arguments(V)...))
    end

    @testset "* algebra" begin
        A = ApplyArray(*,[1 2; 3 4], Vcat(Fill(1,1,3),Fill(2,1,3)))
        @test 2.0A isa ApplyArray
        @test 2.0\A isa ApplyArray
        @test A/2 isa ApplyArray
        @test (2.0A) == 2.0Array(A)
        @test (2.0\A) == 2.0\Array(A)
        @test A/2.0 == Array(A)/2.0
    end

    @testset "* trans" begin
        A = ApplyArray(*,[1 2; 3 4], Vcat(Fill(1,1,3),Fill(2,1,3)))
        V = view(A,1:2,1:2)
        @test MemoryLayout(typeof(V)) isa ApplyLayout{typeof(*)}
        @test V == ApplyArray(V)
        @test MemoryLayout(typeof(V')) isa ApplyLayout{typeof(*)}
        @test V' == ApplyArray(V')
        @test MemoryLayout(typeof(transpose(V))) isa ApplyLayout{typeof(*)}
        @test transpose(V) == ApplyArray(transpose(V))
    end
end


text:
using BinaryBuilder

name = "unpaper"
version = v"6.1.100" # <--- This version number is a lie, (it is v6.1) we just need to bump it to build for experimental platforms

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://www.flameeyes.com/files/unpaper-6.1.tar.xz",
                  "237c84f5da544b3f7709827f9f12c37c346cdf029b1128fb4633f9bafa5cb930"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/unpaper-*/
if [[ "${target}" == *-mingw* ]]; then
    # FFMPEG_jll installs the pkgconfig files in the wrong directory for Windows
    export PKG_CONFIG_PATH="${PKG_CONFIG_PATH}:${libdir}/pkgconfig"
    # Give some hints to the linker
    export LDFLAGS="-L${libdir}"
    export LIBAV_LIBS="-lavformat -lavutil -lavcodec"
fi
update_configure_scripts
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> arch(p) != "armv6l", supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    ExecutableProduct("unpaper", :unpaper),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # Offer a native xsltproc
    HostBuildDependency("XSLT_jll"),
    Dependency("FFMPEG_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
# FFMPEG uses `preferred_gcc_version=v"8"`.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"8", julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using DynamicHMCModels
using Documenter

append!(page_list, [Pair("Functions", "index.md")])

makedocs(root = DOC_ROOT,
    modules = Module[],
    sitename = "DynamicHMCModels.jl",
    authors = "Rob Goedman, Richard Torkar, and contributors.",
    pages = page_list
)

deploydocs(root = DOC_ROOT,
    repo = "github.com/StatisticalRethinkingJulia/DynamicHMCModels.jl.git",
 )


text:
module LinearModel

export jacobian_motion, equation_motion,
jacobian_measurement, equation_measurement, inverse_measurement

H = 0.74
L = 2.83
a = 0.95 + L
b = 0.50

function jacobian_motion(θ,v,α,dt)
    [[1 0 -dt*(v*sin(θ)+v/L*tan(α)*(a*cos(θ)-b*sin(θ)))];
    [0 1 dt*(v*cos(θ)-v/L*tan(α)*(a*sin(θ)+b*cos(θ)))];
    [0 0 1]]
end

function equation_motion(θ,v,α,dt)
    v /= 1-tan(α)*H/L
    [dt*(v*cos(θ)-v/L*tan(α)*(a*sin(θ)+b*cos(θ))),
    dt*(v*sin(θ)+v/L*tan(α)*(a*cos(θ)-b*sin(θ))),
    dt*v/L*tan(α)]
end

function jacobian_measurement(μ, ind)
    n = length(ind)

    x = μ[2*ind .+ 2]
    y = μ[2*ind .+ 3]

    δ = [x y] .- μ[1:2]'
    q = sum(δ.^2,dims=2)

    J = [1 0;0 -1;0 1;1 0]*δ'
    J = [-J;zeros(1,n);-ones(1,n);J]
    J[1:2:end,:] .*= 1 ./ sqrt.(q')
    J[[2,4,8,10],:] .*= 1 ./ q'

    return reshape(J,(2,5,n))
end

function equation_measurement(μ,ind)
    n = length(ind)

    x = μ[2*ind .+ 2]
    y = μ[2*ind .+ 3]

    δ = [x y] .- μ[1:2]'
    q = sqrt.(sum(δ.^2,dims=2))

    [q atan.(δ[:,2],δ[:,1]).-μ[3].+π/2]'
end

function inverse_measurement(μ,z)
    c = cos.(μ[3] .+ z[2,:])
    s = sin.(μ[3] .+ z[2,:])

    return μ[1:2] .+ [(z[1,:].*s)'; (-z[1,:].*c)']
end

end


text:
function test_mat_interface()

  if MPI.Comm_size(MPI.COMM_WORLD) != 1
    return nothing
  end

  @testset "----- Testing Mat interface -----" begin
    vals = PetscScalar[1 2 3; 4 5 6; 7 8 9]
    vals2 = zeros(vals)
    idx = PetscInt[1, 2, 3]
    idy = copy(idx)
    A_j = zeros(3,3)
    A_s = sparse(A_j)
    A_p = PetscMat(3, 3, PETSc2.MATMPIAIJ, MPI.COMM_WORLD)
    SetUp(A_p)

    test_setvalues1(A_j, idx, idy, vals)
    test_setvalues1(A_s, idx, idy, vals)
    test_setvalues1(A_p, idx, idy, vals)

    PetscView(A_j)
    PetscView(A_s)

    test_sizes(A_j)
    test_sizes(A_s)
    test_sizes(A_p)

    @test ( local_indices(A_j) )== 1:3
    @test ( local_indices(A_s) )== 1:3
    @test ( local_indices(A_p) )== 1:3

    fill_zero!(A_j)
    @test isapprox( norm(A_j, Inf), 0.0)
    fill_zero!(A_s)
    @test isapprox( norm(A_s, Inf), 0.0)
    fill_zero!(A_p)
    @test isapprox( norm(A_p, Inf), 0.0)

    set_values1!(A_p, idx, idy, vals)
    assembly_begin(A_p, MAT_FINAL_ASSEMBLY)
    assembly_end(A_p, MAT_FINAL_ASSEMBLY)
    scale!(A_p, 2.0)
    get_values1!(A_p, idx, idy, vals2)
    @test isapprox( norm(vals2 - 2*vals), 0.0) atol=1e-14

    set_values1!(A_p, idx, idy, vals)
    assembly_begin(A_p, MAT_FINAL_ASSEMBLY)
    assembly_end(A_p, MAT_FINAL_ASSEMBLY)

    @test isapprox( vecnorm(A_p), vecnorm(vals)) atol=1e-13

    # test multiplication
    set_values1!(A_p, idx, idy, vals.')
    assembly_begin(A_p, MAT_FINAL_ASSEMBLY)
    assembly_end(A_p, MAT_FINAL_ASSEMBLY)

    x = PetscVec(3, VECMPI, MPI.COMM_WORLD)
    b = PetscVec(3, VECMPI, MPI.COMM_WORLD)
    vals_vec = PetscScalar[1, 2, 3]
    vals_vec2 = zeros(vals_vec)
    set_values1!(x, idx, vals_vec)

    A_mul_B!(b, A_p, x)
    get_values1!(b, idx, vals_vec2)
    b2 = vals*vals_vec
    @test isapprox( norm(vals_vec2 - b2), 0.0) atol=1e-14

    At_mul_B!(b, A_p, x)
    get_values1!(b, idx, vals_vec2)
    b2 = vals.'*vals_vec

    @test isapprox( norm(vals_vec2 - b2), 0.0) atol=1e-14

    # test diagonal_shift
    set_values1!(A_p, idx, idy, vals.')
    assembly_begin(A_p, MAT_FINAL_ASSEMBLY)
    assembly_end(A_p, MAT_FINAL_ASSEMBLY)
    A_j[:, :] = vals
    A_s[:, :] = vals

    set_values1!(x, idx, vals_vec)
    alpha = 2
    b2 = vals*vals_vec + alpha*vals_vec  # (alpha * I + A)*vals_vec
    diagonal_shift!(A_p, alpha)
    A_mul_B!(b, A_p, x)
    get_values1!(b, idx, vals_vec2)
    @test isapprox( norm(vals_vec2 - b2), 0.0) atol=1e-14

    diagonal_shift!(A_j, alpha)
    b3 = A_j*vals_vec
    @test isapprox( norm(vals_vec2 - b3), 0.0) atol=1e-14

    diagonal_shift!(A_s, alpha)
    b4 = A_s*vals_vec
    @test isapprox( norm(vals_vec2 - b4), 0.0) atol=1e-14

    PetscDestroy(A_p)
    PetscDestroy(x)
    PetscDestroy(b)
  end

  return nothing
end

function test_setvalues1(A, idx, idy, vals)

  vals2 = zeros(vals)

  # check set_values1!, both modes
  set_values1!(A, idx, idy, vals)
  assembly_begin(A, MAT_FINAL_ASSEMBLY)
  assembly_end(A, MAT_FINAL_ASSEMBLY)
  get_values1!(A, idx, idy, vals2)
  @test isapprox( norm(vals - vals2), 0.0)

  set_values1!(A, idx, idy, vals)
  assembly_begin(A, MAT_FINAL_ASSEMBLY)
  assembly_end(A, MAT_FINAL_ASSEMBLY)
  get_values1!(A, idx, idy, vals2)
  @test isapprox( norm(vals - vals2), 0.0)

  set_values1!(A, idx, idy, vals, ADD_VALUES)
  assembly_begin(A, MAT_FINAL_ASSEMBLY)
  assembly_end(A, MAT_FINAL_ASSEMBLY)
  get_values1!(A, idx, idy, vals2)
  @test isapprox( norm(2*vals - vals2), 0.0)

  return nothing
end

function test_sizes(A)

  @test ( size_local(A) )== (3, 3)
  @test ( size_global(A) )== (3, 3)

  return nothing
end

test_mat_interface()


text:
"""
    create_fibonacci_sphere(npoint)

Create `npoint` evenly distributed points on the sphere.
"""
function create_fibonacci_sphere(npoint)
    φ = (1 + √5) / 2
    inds = collect(0:npoint-1)'
    θ = 2π / φ * inds
    ϕ = acos.(1 .- 2 .* (inds .+ 0.5) ./ npoint)

    # Create Fibonacci directions (uniformly distributed on unit sphere)
    [
        cos.(θ) .* sin.(ϕ)
        sin.(θ) .* sin.(ϕ)
        cos.(ϕ)
    ]
end


text:
# function sor_cyl_gpu!(
@kernel function sor_cyl_gpu!(
        potential::AbstractArray{T, 4},
        point_types::AbstractArray{PointType, 4},
        volume_weights::AbstractArray{T, 4},
        q_eff_imp::AbstractArray{T, 4},
        q_eff_fix::AbstractArray{T, 4},
        ϵ_r::AbstractArray{T, 3},
        geom_weights::NTuple{3, <:AbstractArray{T, 2}},
        sor_const::AbstractArray{T, 1},
        update_even_points::Bool,
        depletion_handling_enabled::Bool,
        is_weighting_potential::Bool,
        only2d::Bool
    ) where {T}
    # # GPU Indexing
    # linear_idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x # CUDA.jl
    linear_idx = @index(Global) # KernelAbstractions.jl

    eff_size = broadcast(idim -> 2:size(potential, idim)-1, (1, 2, 3))
    if linear_idx <= prod(length.(eff_size))
        i1, i2, i3 = Tuple(CartesianIndices(eff_size)[linear_idx])
        # Comparison to CPU indecies: (Cyl / Car)
        # i3 <-> idx3 / ir / iz
        # i2 <-> idx2 / iφ / iy
        # i1 <-> idx1 / iz / ix
        in3 = i3 - 1
        in2 = i2 - 1
        in1 = nidx(i1, update_even_points, iseven(i2 + i3))
        ir = i3
        inr = in3
        iφ = i2
        inφ = in2
        iz = i1
        inz = in1
        izr = get_rbidx_right_neighbour(iz, update_even_points, iseven(i2 + i3))

        rb_tar_idx, rb_src_idx = update_even_points ? (rb_even::Int, rb_odd::Int) : (rb_odd::Int, rb_even::Int)

        pwwrr               = geom_weights[1][1, inr]
        pwwrl               = geom_weights[1][2, inr]
        r_inv_pwΔmpr        = geom_weights[1][3, inr]
        Δr_ext_inv_r_pwmprr = geom_weights[1][4, inr]
        Δr_ext_inv_l_pwmprl = geom_weights[1][5, inr]
        Δmpr_squared        = geom_weights[1][6, inr]

        pwwφr        = geom_weights[2][1, inφ]
        pwwφl        = geom_weights[2][2, inφ]
        pwΔmpφ       = geom_weights[2][3, inφ]
        Δφ_ext_inv_r = geom_weights[2][4,  iφ]
        Δφ_ext_inv_l = geom_weights[2][4, inφ]

        if inr == 1
            pwwφr = T(0.5)
            pwwφl = T(0.5)
            pwΔmpφ = T(2π)
            Δφ_ext_inv_r = T(1/pwΔmpφ)
            Δφ_ext_inv_l = Δφ_ext_inv_r
        end
        pwwrr_pwwφr = pwwrr * pwwφr
        pwwrr_pwwφl = pwwrr * pwwφl
        pwwrl_pwwφr = pwwrl * pwwφr
        pwwrl_pwwφl = pwwrl * pwwφl

        Δr_ext_inv_r_pwmprr_pwΔmpφ = Δr_ext_inv_r_pwmprr * pwΔmpφ
        Δr_ext_inv_l_pwmprl_pwΔmpφ = Δr_ext_inv_l_pwmprl * pwΔmpφ
        pwΔmpφ_Δmpr_squared = pwΔmpφ * Δmpr_squared
        r_inv_pwΔmpr_Δφ_ext_inv_r = r_inv_pwΔmpr * Δφ_ext_inv_r
        r_inv_pwΔmpr_Δφ_ext_inv_l = r_inv_pwΔmpr * Δφ_ext_inv_l

        pwwzr::T        = geom_weights[3][1, inz]
        pwwzl::T        = geom_weights[3][2, inz]
        pwΔmpz::T       = geom_weights[3][3, inz]
        Δz_ext_inv_r::T = geom_weights[3][4, inz + 1]
        Δz_ext_inv_l::T = geom_weights[3][4, inz]

        ϵ_rrr::T = ϵ_r[  ir,  iφ, inz + 1]
        ϵ_rlr::T = ϵ_r[  ir, inφ, inz + 1]
        ϵ_rrl::T = ϵ_r[  ir,  iφ, inz ]
        ϵ_rll::T = ϵ_r[  ir, inφ, inz ]
        ϵ_lrr::T = ϵ_r[ inr,  iφ, inz + 1]
        ϵ_llr::T = ϵ_r[ inr, inφ, inz + 1]
        ϵ_lrl::T = ϵ_r[ inr,  iφ, inz ]
        ϵ_lll::T = ϵ_r[ inr, inφ, inz ]

        pwwφr_pwwzr::T = pwwφr * pwwzr
        pwwφl_pwwzr::T = pwwφl * pwwzr
        pwwφr_pwwzl::T = pwwφr * pwwzl
        pwwφl_pwwzl::T = pwwφl * pwwzl
        pwwrl_pwwzr::T = pwwrl * pwwzr
        pwwrr_pwwzr::T = pwwrr * pwwzr
        pwwrl_pwwzl::T = pwwrl * pwwzl
        pwwrr_pwwzl::T = pwwrr * pwwzl

        # right weight in r: wrr
        wrr::T = ϵ_rrr * pwwφr_pwwzr
        wrr    = muladd(ϵ_rlr, pwwφl_pwwzr, wrr)
        wrr    = muladd(ϵ_rrl, pwwφr_pwwzl, wrr)
        wrr    = muladd(ϵ_rll, pwwφl_pwwzl, wrr)
        # left weight in r: wrr
        wrl::T = ϵ_lrr * pwwφr_pwwzr
        wrl    = muladd(ϵ_llr, pwwφl_pwwzr, wrl)
        wrl    = muladd(ϵ_lrl, pwwφr_pwwzl, wrl)
        wrl    = muladd(ϵ_lll, pwwφl_pwwzl, wrl)
        # right weight in φ: wφr
        wφr::T = ϵ_lrr * pwwrl_pwwzr
        wφr    = muladd(ϵ_rrr, pwwrr_pwwzr, wφr)
        wφr    = muladd(ϵ_lrl, pwwrl_pwwzl, wφr)
        wφr    = muladd(ϵ_rrl, pwwrr_pwwzl, wφr)
        # left weight in φ: wφl
        wφl::T = ϵ_llr * pwwrl_pwwzr
        wφl    = muladd(ϵ_rlr, pwwrr_pwwzr, wφl)
        wφl    = muladd(ϵ_lll, pwwrl_pwwzl, wφl)
        wφl    = muladd(ϵ_rll, pwwrr_pwwzl, wφl)
        # right weight in z: wzr
        wzr::T = ϵ_rrr * pwwrr_pwwφr
        wzr    = muladd(ϵ_rlr, pwwrr_pwwφl, wzr)
        wzr    = muladd(ϵ_lrr, pwwrl_pwwφr, wzr)
        wzr    = muladd(ϵ_llr, pwwrl_pwwφl, wzr)
        # left weight in z: wzr
        wzl::T = ϵ_rrl * pwwrr_pwwφr
        wzl    = muladd(ϵ_rll, pwwrr_pwwφl, wzl)
        wzl    = muladd(ϵ_lrl, pwwrl_pwwφr, wzl)
        wzl    = muladd(ϵ_lll, pwwrl_pwwφl, wzl)

        wrr *= Δr_ext_inv_r_pwmprr_pwΔmpφ * pwΔmpz
        wrl *= Δr_ext_inv_l_pwmprl_pwΔmpφ * pwΔmpz
        wφr *= r_inv_pwΔmpr_Δφ_ext_inv_r * pwΔmpz
        wφl *= r_inv_pwΔmpr_Δφ_ext_inv_l * pwΔmpz
        wzr *= Δz_ext_inv_r * pwΔmpφ_Δmpr_squared
        wzl *= Δz_ext_inv_l * pwΔmpφ_Δmpr_squared

        old_potential::T = potential[iz, iφ, ir, rb_tar_idx]

        vrr::T = potential[     iz,     iφ, ir + 1, rb_src_idx]
        vrl::T = potential[     iz,     iφ,    inr, rb_src_idx]
        vφr::T = only2d ? old_potential : potential[ iz, iφ + 1, ir, rb_src_idx]
        vφl::T = only2d ? old_potential : potential[ iz,    inφ, ir, rb_src_idx]
        vzr::T = potential[    izr,     iφ,     ir, rb_src_idx]
        vzl::T = potential[izr - 1,     iφ,     ir, rb_src_idx]

        new_potential::T = is_weighting_potential ? zero(T) : (q_eff_imp[iz, iφ, ir, rb_tar_idx] + q_eff_fix[iz, iφ, ir, rb_tar_idx])

        new_potential = muladd( wrr, vrr, new_potential)
        new_potential = muladd( wrl, vrl, new_potential)
        new_potential = muladd( wφr, vφr, new_potential)
        new_potential = muladd( wφl, vφl, new_potential)
        new_potential = muladd( wzr, vzr, new_potential)
        new_potential = muladd( wzl, vzl, new_potential)

        new_potential *= volume_weights[iz, iφ, ir, rb_tar_idx]

        new_potential -= old_potential
        new_potential = muladd(new_potential, sor_const[inr], old_potential)

        if depletion_handling_enabled
            if inr == 1 vrl = vrr end
            vmin::T = min(vrr, vrl, vφr, vφl, vzr, vzl)
            vmax::T = max(vrr, vrl, vφr, vφl, vzr, vzl)
            new_point_type = point_types[iz, iφ, ir, rb_tar_idx]
            if new_potential < vmin || new_potential > vmax
                new_potential -= q_eff_imp[iz, iφ, ir, rb_tar_idx] * volume_weights[iz, iφ, ir, rb_tar_idx] * sor_const[inr]
                if (point_types[iz, iφ, ir, rb_tar_idx] & undepleted_bit == 0 &&
                    point_types[iz, iφ, ir, rb_tar_idx] & pn_junction_bit > 0) # undepleted_bit
                    # point_types[iz, iφ, ir, rb_tar_idx] += undepleted_bit
                    new_point_type += undepleted_bit
                end # mark this point as undepleted
            elseif point_types[iz, iφ, ir, rb_tar_idx] & undepleted_bit > 0
                new_point_type -= undepleted_bit
            end
            point_types[iz, iφ, ir, rb_tar_idx] = new_point_type
        end

        potential[iz, iφ, ir, rb_tar_idx]::T = ifelse(point_types[iz, iφ, ir, rb_tar_idx] & update_bit > 0, new_potential, old_potential)
    end
end


text:
# Import seeding generator for reproducibility
using StableRNGs

"""
    Sigmoid activation function
"""
function sigmoid(Z)
    A = 1 ./ (1 .+ exp.(.-Z))
    return (A = A, Z = Z)
end

"""
    ReLU activation function
"""
function relu(Z)
    A = max.(0, Z)
    return (A = A, Z = Z)
end

"""
    Funtion to initialise the parameters or weights of the desired network.
"""
function initialise_model_weights(layer_dims, seed)
    params = Dict()

    # Build a dictionary of initialised weights and bias units
    for l=2:length(layer_dims)
        params[string("W_", (l-1))] = rand(StableRNG(seed), layer_dims[l], layer_dims[l-1]) * sqrt(2 / layer_dims[l-1])
        params[string("b_", (l-1))] = zeros(layer_dims[l], 1)
    end

    return params
end

"""
    Make a linear forward calculation
"""
function linear_forward(A, W, b)
    # Make a linear forward and return inputs as cache
    Z = (W * A) .+ b
    cache = (A, W, b)

    @assert size(Z) == (size(W, 1), size(A, 2))

    return (Z = Z, cache = cache)
end

"""
    Make a forward activation from a linear forward.
"""
function linear_forward_activation(A_prev, W, b, activation_function="relu")
    @assert activation_function ∈ ("sigmoid", "relu")
    Z, linear_cache = linear_forward(A_prev, W, b)

    if activation_function == "sigmoid"
        A, activation_cache = sigmoid(Z)
    end

    if activation_function == "relu"
        A, activation_cache = relu(Z)
    end

    cache = (linear_step_cache=linear_cache, activation_step_cache=activation_cache)

    @assert size(A) == (size(W, 1), size(A_prev, 2))

    return A, cache
end

"""
    Forward the design matrix through the network layers using the parameters.
"""
function forward_propagate_model_weights(DMatrix, parameters)
    master_cache = []
    A = DMatrix
    L = Int(length(parameters) / 2)

    # Forward propagate until the last (output) layer
    for l = 1 : (L-1)
        A_prev = A
        A, cache = linear_forward_activation(A_prev,
                                             parameters[string("W_", (l))],
                                             parameters[string("b_", (l))],
                                             "relu")
        push!(master_cache , cache)
    end

    # Make predictions in the output layer
    Ŷ, cache = linear_forward_activation(A,
                                         parameters[string("W_", (L))],
                                         parameters[string("b_", (L))],
                                         "sigmoid")
    push!(master_cache, cache)

    return Ŷ, master_cache
end

"""
    Computes the log loss (binary cross entropy) of the current predictions.
"""
function calculate_cost(Ŷ, Y)
    m = size(Y, 2)
    epsilon = eps(1.0)

    # Deal with log(0) scenarios
    Ŷ_new = [max(i, epsilon) for i in Ŷ]
    Ŷ_new = [min(i, 1-epsilon) for i in Ŷ_new]

    cost = -sum(Y .* log.(Ŷ_new) + (1 .- Y) .* log.(1 .- Ŷ_new)) / m
    return cost
end

"""
    Derivative of the Sigmoid function.
"""
function sigmoid_backwards(∂A, activated_cache)
    s = sigmoid(activated_cache).A
    ∂Z = ∂A .* s .* (1 .- s)

    @assert (size(∂Z) == size(activated_cache))
    return ∂Z
end

"""
    Derivative of the ReLU function.
"""
function relu_backwards(∂A, activated_cache)
    return ∂A .* (activated_cache .> 0)
end

"""
    Partial derivatives of the components of linear forward function
    using the linear output (∂Z) and caches of these components (cache).
"""
function linear_backward(∂Z, cache)
    # Unpack cache
    A_prev , W , b = cache
    m = size(A_prev, 2)

    # Partial derivates of each of the components
    ∂W = ∂Z * (A_prev') / m
    ∂b = sum(∂Z, dims = 2) / m
    ∂A_prev = (W') * ∂Z

    @assert (size(∂A_prev) == size(A_prev))
    @assert (size(∂W) == size(W))
    @assert (size(∂b) == size(b))

    return ∂W , ∂b , ∂A_prev
end

"""
    Unpack the linear activated caches (cache) and compute their derivatives
    from the applied activation function.
"""
function linear_activation_backward(∂A, cache, activation_function="relu")
    @assert activation_function ∈ ("sigmoid", "relu")

    linear_cache , cache_activation = cache

    if (activation_function == "relu")

        ∂Z = relu_backwards(∂A , cache_activation)
        ∂W , ∂b , ∂A_prev = linear_backward(∂Z , linear_cache)

    elseif (activation_function == "sigmoid")

        ∂Z = sigmoid_backwards(∂A , cache_activation)
        ∂W , ∂b , ∂A_prev = linear_backward(∂Z , linear_cache)

    end

    return ∂W , ∂b , ∂A_prev
end

"""
    Compute the gradients (∇) of the parameters (master_cache) of the constructed model
    with respect to the cost of predictions (Ŷ) in comparison with actual output (Y).
"""
function back_propagate_model_weights(Ŷ, Y, master_cache)
    # Initiate the dictionary to store the gradients for all the components in each layer
    ∇ = Dict()

    L = length(master_cache)
    Y = reshape(Y , size(Ŷ))

    # Partial derivative of the output layer
    ∂Ŷ = (-(Y ./ Ŷ) .+ ((1 .- Y) ./ ( 1 .- Ŷ)))
    current_cache = master_cache[L]

    # Backpropagate on the layer preceeding the output layer
    ∇[string("∂W_", (L))], ∇[string("∂b_", (L))], ∇[string("∂A_", (L-1))] = linear_activation_backward(∂Ŷ,
                                                                                current_cache,
                                                                                "sigmoid")
    # Go backwards in the layers and compute the partial derivates of each component.
    for l=reverse(0:L-2)
        current_cache = master_cache[l+1]
        ∇[string("∂W_", (l+1))], ∇[string("∂b_", (l+1))], ∇[string("∂A_", (l))] = linear_activation_backward(∇[string("∂A_", (l+1))],
                                                                                current_cache,
                                                                                "relu")
    end

    # Return the gradients of the network
    return ∇
end

"""
    Update the paramaters of the model using the gradients (∇)
    and the learning rate (η).
"""
function update_model_weights(parameters, ∇, η)

    L = Int(length(parameters) / 2)

    # Update the parameters (weights and biases) for all the layers
    for l = 0: (L-1)
        parameters[string("W_", (l + 1))] -= η .* ∇[string("∂W_", (l + 1))]
        parameters[string("b_", (l + 1))] -= η .* ∇[string("∂b_", (l + 1))]
    end

    return parameters
end

"""
    Check the accuracy between predicted values (Ŷ) and the true values(Y).
"""
function assess_accuracy(Ŷ , Y)
    @assert size(Ŷ) == size(Y)
    return sum((Ŷ .> 0.5) .== Y) / length(Y)
end

"""
    Train the network using the desired architecture that best possible
    matches the training inputs (DMatrix) and their corresponding ouptuts(Y)
    over some number of iterations (epochs) and a learning rate (η).
"""
function train_network(layer_dims , DMatrix, Y;  η=0.001, epochs=1000, seed=2020, verbose=true)
    # Initiate an empty container for cost, iterations, and accuracy at each iteration
    costs = []
    iters = []
    accuracy = []

    # Initialise random weights for the network
    params = initialise_model_weights(layer_dims, seed)

    # Train the network
    for i = 1:epochs

        Ŷ , caches  = forward_propagate_model_weights(DMatrix, params)
        cost = calculate_cost(Ŷ, Y)
        acc = assess_accuracy(Ŷ, Y)
        ∇  = back_propagate_model_weights(Ŷ, Y, caches)
        params = update_model_weights(params, ∇, η)

        if verbose
            println("Iteration -> $i, Cost -> $cost, Accuracy -> $acc")
        end

        # Update containers for cost, iterations, and accuracy at the current iteration (epoch)
        push!(iters , i)
        push!(costs , cost)
        push!(accuracy , acc)
    end
        return (cost = costs, iterations = iters, accuracy = accuracy, parameters = params)
end

max_stars_repo_name: PyDataBlog/ML-From-Scratch-With-Julia

max_issues_repo_name: PyDataBlog/ML-From-Scratch-With-Julia

max_forks_repo_name: PyDataBlog/ML-From-Scratch-With-Julia


text:
#--- check multi-species
maxdeg = 5
ord = 3
Pr = transformed_jacobi(maxdeg, trans, rcut; pcut = 2)
species = [:C, :O, :H]
P1 = ACE.RnYlm1pBasis(Pr; species = [:C, :O, :H], D = D)
dagbasis = ACE.PIBasis(P1, ord, D, maxdeg)
basis = standardevaluator(dagbasis)
Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, species)
AA = evaluate(basis, Rs, Zs, z0)
println(@test(length(basis) == length(AA)))
dAA = evaluate_d(basis, Rs, Zs, z0)
println(@test(length(basis) == length(AA)))

# construct a basis with dag-evaluator, and check they are identical!!!
AAdag = evaluate(dagbasis, Rs, Zs, z0)
println(@test(AA ≈ AAdag))

dAA = evaluate_d(basis, Rs, Zs, z0)
dAAdag = evaluate_d(dagbasis, Rs, Zs, z0)
println(@test dAA ≈ dAAdag)

println(@test all(JuLIP.Testing.test_fio(dagbasis)))
#---

@info("Check several properties of PIBasis")
for species in (:X, :Si, [:C, :O, :H]), N = 1:5
   local AA, AAdag, dAA, dAAdag, dagbasis, basis, Rs, Zs, z0
   maxdeg = 7
   Nat = 15
   P1 = ACE.RnYlm1pBasis(Pr; species = species)
   dagbasis = ACE.PIBasis(P1, N, D, maxdeg)
   basis = standardevaluator(dagbasis)
   @info("species = $species; N = $N; length = $(length(basis))")
   @info("test (de-)serialisation")
   # only require dagbasis to deserialize correctly
   println(@test all(JuLIP.Testing.test_fio(dagbasis)))
   @info("Check Permutation invariance")
   for ntest = 1:20
      Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, species)
      p = randperm(length(Rs))
      print_tf(@test(evaluate(basis, Rs, Zs, z0) ≈
                     evaluate(basis, Rs[p], Zs[p], z0)))
   end
   println()
   @info("Check gradients")
   for ntest = 1:20
      Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, species)
      AA = evaluate(basis, Rs, Zs, z0)
      dAA = evaluate_d(basis, Rs, Zs, z0)
      Us = [ rand(eltype(Rs)) .- 0.5 for _=1:length(Rs) ]
      dAA_dUs = transpose.(dAA) * Us
      errs = []
      for p = 2:12
         h = 0.1^p
         AA_h = evaluate(basis, Rs + h * Us, Zs, z0)
         dAA_h = (AA_h - AA) / h
         # @show norm(dAA_h - dAA_dUs, Inf)
         push!(errs, norm(dAA_h - dAA_dUs, Inf))
      end
      success = (/(extrema(errs)...) < 1e-3) || (minimum(errs) < 1e-10)
      print_tf(@test success)
   end

   println()
   @info("Check Standard=DAG Evaluator")
   for ntest = 1:20
      Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, species)
      AA = evaluate(basis, Rs, Zs, z0)
      AAdag = evaluate(dagbasis, Rs, Zs, z0)
      print_tf(@test AA ≈ AAdag)

      dAA = evaluate_d(basis, Rs, Zs, z0)
      dAAdag = evaluate_d(dagbasis, Rs, Zs, z0)
      print_tf(@test dAA ≈ dAAdag)
   end
   println()
end
println()


text:
using CSyntax
using STBImage.LibSTBImage
using Dates

@static if Sys.isapple()
    const VERSION_MAJOR = 4
    const VERSION_MINOR = 1
end

include(joinpath(@__DIR__, "glutils.jl"))
include(joinpath(@__DIR__, "camera.jl"))

# init window
width, height = 640, 480
window = startgl(width, height)

glEnable(GL_DEPTH_TEST)
glDepthFunc(GL_LESS)

# camera
camera = PerspectiveCamera()
setposition!(camera, [0.0, 0.0, 2.0])

# screen capture
function screen_capture()
    buffer = zeros(Cuchar, 3 * width * height)
    glReadPixels(0, 0, width, height, GL_RGB, GL_UNSIGNED_BYTE, buffer)
    filename = joinpath(@__DIR__, "screenshot_"*"$(Dates.format(Dates.now(), "yyyy_mm_ddTH_M_S"))"*".png")
    stbi_flip_vertically_on_write(true)
    stbi_write_png(filename, width, height, 3, buffer, 3 * width) != 0 && return true
    return false
end

# load texture
function load_texture(path::AbstractString)
    x, y, n = Cint(0), Cint(0), Cint(0)
    force_channels = 4
    stbi_set_flip_vertically_on_load(true)
    tex_data = @c stbi_load(path, &x, &y, &n, force_channels)
    if tex_data == C_NULL
        @error "could not load $path."
        return nothing
    end
    ( ( x & ( x - 1 ) ) != 0 || ( y & ( y - 1 ) ) != 0 ) && @warn "texture $path is not power-of-2 dimensions."
    id = GLuint(0)
    @c glGenTextures(1, &id)
    glActiveTexture(GL_TEXTURE0)
    glBindTexture(GL_TEXTURE_2D, id)
    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, x, y, 0, GL_RGBA, GL_UNSIGNED_BYTE, tex_data)
    glGenerateMipmap(GL_TEXTURE_2D)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_LINEAR)
    return id
end
load_texture(joinpath(@__DIR__, "skulluvmap.png"))

# vertex data
points = GLfloat[-0.5, -0.5, 0.0,
                  0.5, -0.5, 0.0,
                  0.5,  0.5, 0.0,
                  0.5,  0.5, 0.0,
                 -0.5,  0.5, 0.0,
                 -0.5, -0.5, 0.0]

texcoords = GLfloat[0.0, 0.0,
                    1.0, 0.0,
                    1.0, 1.0,
                    1.0, 1.0,
                    0.0, 1.0,
                    0.0, 0.0]

# create VBO
points_vbo = GLuint(0)
@c glGenBuffers(1, &points_vbo)
glBindBuffer(GL_ARRAY_BUFFER, points_vbo)
glBufferData(GL_ARRAY_BUFFER, sizeof(points), points, GL_STATIC_DRAW)

texcoords_vbo = GLuint(0)
@c glGenBuffers(1, &texcoords_vbo)
glBindBuffer(GL_ARRAY_BUFFER, texcoords_vbo)
glBufferData(GL_ARRAY_BUFFER, sizeof(texcoords), texcoords, GL_STATIC_DRAW)

# create VAO
vao = GLuint(0)
@c glGenVertexArrays(1, &vao)
glBindVertexArray(vao)
glBindBuffer(GL_ARRAY_BUFFER, points_vbo)
glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 0, C_NULL)
glBindBuffer(GL_ARRAY_BUFFER, texcoords_vbo)
glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 0, C_NULL)
glEnableVertexAttribArray(0)
glEnableVertexAttribArray(1)

# load and compile shaders from file
vert_shader = createshader(joinpath(@__DIR__, "screencap.vert"), GL_VERTEX_SHADER)
frag_shader = createshader(joinpath(@__DIR__, "screencap.frag"), GL_FRAGMENT_SHADER)

# link program
shader_prog = createprogram(vert_shader, frag_shader)

view_loc = glGetUniformLocation(shader_prog, "view")
proj_loc = glGetUniformLocation(shader_prog, "proj")
glUseProgram(shader_prog)
glUniformMatrix4fv(view_loc, 1, GL_FALSE, get_view_matrix(camera))
glUniformMatrix4fv(proj_loc, 1, GL_FALSE, get_projective_matrix(window, camera))

# enable cull face
glEnable(GL_CULL_FACE)
glCullFace(GL_BACK)
glFrontFace(GL_CCW)
# set background color to gray
glClearColor(0.2, 0.2, 0.2, 1.0)

let
updatefps = FPSCounter()
# render
while !GLFW.WindowShouldClose(window)
    updatefps(window)
    # clear drawing surface
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glViewport(0, 0, GLFW.GetFramebufferSize(window)...)
    # drawing
    glUseProgram(shader_prog)
    glBindVertexArray(vao)
    glDrawArrays(GL_TRIANGLES, 0, 6)
    # check and call events
    GLFW.PollEvents()
    # screen capture
    if GLFW.GetKey(window, GLFW.KEY_SPACE)
        println("screen captured")
        screen_capture()
    end
    # move camera
    updatecamera!(window, camera)
    glUniformMatrix4fv(view_loc, 1, GL_FALSE, get_view_matrix(camera))
    # swap the buffers
    GLFW.SwapBuffers(window)
end
end # let

GLFW.DestroyWindow(window)

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
#=

 Test cases for reverse interval absolute value function.

 Copyright 2015-2016 Oliver Heimlich (oheim@posteo.de)

 Copying and distribution of this file, with or without modification,
 are permitted in any medium without royalty provided the copyright
 notice and this notice are preserved.  This file is offered as-is,
 without any warranty.

=#
#Language imports

#Test library imports
using Test

#Arithmetic library imports
using IntervalArithmetic

#Preamble
setprecision(53)
setprecision(Interval, Float64)
setrounding(Interval, :tight)
# Set full format, and show decorations
@format full

@testset "minimal.absRevBin_test" begin

end
# FactCheck.exitstatus()


text:
const PROJ_IDER = "iJR904"
const BIOMASS_IDER = "BiomassEcoli"
const EX_GLC_IDER = "EX_glc_LPAREN_e_RPAREN_"
const ATPM_IDER = "ATPM"
const COST_IDER = "tot_cost"
const ABS_MAX_BOUND = 100.0
const MAX_CONC = 9999.0


text:
using StateSpaceReconstruction
import DelayEmbeddings.Dataset

import StateSpaceReconstruction:
    MutablePointTriangulation, PointTriangulation,
    MutableDatasetTriangulation, DatasetTriangulation,
    MutableMatrixTriangulation, MatrixTriangulation,
    DelaunayTriangulation,
    triangulate, triangulate_m

import StaticArrays: SVector, MVector

using DelayEmbeddings

pts = [rand(3) for i = 1:30];

@testset "Triangulation structs" begin

    @test triangulate(Dataset(pts)) isa DatasetTriangulation{3,Float64}
    @test triangulate([SVector{3, Float64}(x) for x in pts]) isa PointTriangulation{3, Float64}
    @test triangulate([MVector{3, Float64}(x) for x in pts]) isa PointTriangulation{3, Float64}
    @test triangulate(hcat(pts...,)) isa MatrixTriangulation{3,Float64}
    @test triangulate(transpose(hcat(pts...,))) isa MatrixTriangulation{3,Float64}

    @test triangulate_m(Dataset(pts)) isa MutableDatasetTriangulation{3,Float64}
    @test triangulate_m([SVector{3, Float64}(x) for x in pts]) isa MutablePointTriangulation{3, Float64}
    @test triangulate_m([MVector{3, Float64}(x) for x in pts]) isa MutablePointTriangulation{3, Float64}
    @test triangulate_m(hcat(pts...,)) isa MutableMatrixTriangulation{3,Float64}
    @test triangulate_m(transpose(hcat(pts...,))) isa MutableMatrixTriangulation{3,Float64}
end


text:
#datalakedeletefile
type DataLakeDeleteFileRequest
    DataStoreName::String
    TableName::String
    FilePath::String
end
function delete_datalake_file!(datastorename::String,tablename::String,filename_including_partition::AbstractString)
    table_def=get_table_definition(storename=datastorename,tablename=tablename)
    schematype=filter(tuple->last(tuple)==table_def.SchemaType,collect(schema_type_map))[1][1]
    folderpath=lowercase("/managed/$schematype/table/$tablename/")
    request=DataLakeDeleteFileRequest(datastorename,tablename,folderpath*strip(filename_including_partition,'/'))
    response=NeuroJulia.neurocall("80","DataMovementService","DataLakeDeleteFile",request)
    check_request=Dict("JobId"=>response["JobId"])
    status=0
    errormsg=""
    while(status==0)
        sleep(1)
        response_c=NeuroJulia.neurocall("80","DataMovementService","CheckJob",check_request)
        status=response_c["Status"]
        if status>1
            errormsg=response_c["Message"]
        end
    end

    NeuroJulia.neurocall("80","DataMovementService","FinaliseJob",check_request)

    if status!=1
        error("Neuroverse error: " * errormsg)
    end
    return nothing
end
#listdatalaketablefiles
type ListDataLakeTableFilesRequest
    DataStoreName::String
    TableName::String
end

function list_datalake_table_files_with_partitions(datastorename::String,tablename::String)
    request=ListDataLakeTableFilesRequest(datastorename,tablename)
    files=NeuroJulia.neurocall("80","DataMovementService","ListDataLakeTableFiles",request)["Files"]
    return [split(files[i],lowercase(tablename))[2] for i=1:length(files)]
end

type GetLinesInDataLakeCsvFileRequest
    DataStoreName::String
    TableName::String
    FilePath::String
end

function get_lines_in_datalake_csv(datastorename::String,tablename::String,filename_including_partition::AbstractString)
    table_def=get_table_definition(storename=datastorename,tablename=tablename)
    schematype=filter(tuple->last(tuple)==table_def.SchemaType,collect(schema_type_map))[1][1]
    folderpath=lowercase("/managed/$schematype/table/$tablename/")
    request=GetLinesInDataLakeCsvFileRequest(datastorename,tablename,folderpath*strip(filename_including_partition,'/'))
    response=NeuroJulia.neurocall("80","DataMovementService","GetLinesInDataLakeCsvFile",request)
    check_request=Dict("JobId"=>response["JobId"])
    status=0
    errormsg=""
    while(status==0)
        sleep(1)
        response_c=NeuroJulia.neurocall("80","DataMovementService","CheckJob",check_request)
        status=response_c["Status"]
        if status>=1
            errormsg=response_c["Message"]
        end
    end

    NeuroJulia.neurocall("80","DataMovementService","FinaliseJob",check_request)

    if status!=1
        error("Neuroverse error: " * errormsg)
    end
    return parse(errormsg)
end

type DataLakeReChunkCsvFileRequest
    FromDataStoreName::String
    FromTableName::String
    FilePath::String
    ToTableName::String
end

function rechunk_datalake_csv(datastorename::String,fromtablename::String,filename_including_partition::AbstractString,totablename::String)
    table_def=get_table_definition(storename=datastorename,tablename=fromtablename)
    schematype=filter(tuple->last(tuple)==table_def.SchemaType,collect(schema_type_map))[1][1]
    folderpath=lowercase("/managed/$schematype/table/$fromtablename/")
    request=DataLakeReChunkCsvFileRequest(datastorename,fromtablename,folderpath*strip(filename_including_partition,'/'),totablename)
    response=NeuroJulia.neurocall("80","DataMovementService","DataLakeReChunkCsvFile",request)
    check_request=Dict("JobId"=>response["JobId"])
    status=0
    errormsg=""
    while(status==0)
        sleep(1)
        response_c=NeuroJulia.neurocall("80","DataMovementService","CheckJob",check_request)
        status=response_c["Status"]
        if status>1
            errormsg=response_c["Message"]
        end
    end

    NeuroJulia.neurocall("80","DataMovementService","FinaliseJob",check_request)

    if status!=1
        error("Neuroverse error: " * errormsg)
    end

    files=list_datalake_table_files_with_partitions(datastorename,totablename)

    outfiles=String[]

    for file in files
        if contains(file,response["JobId"])
            push!(outfiles,file)
        end
    end

    return outfiles
end

max_stars_repo_name: SnowdenNeuroverse/NeuroJulia

max_issues_repo_name: SnowdenNeuroverse/NeuroJulia

max_forks_repo_name: SnowdenNeuroverse/NeuroJulia


text:
export objective, error_metric, impute, impute_missing

### OBJECTIVE FUNCTION EVALUATION FOR MPCA
function objective(glrm::GLRM, X::Array{Float64,2}, Y::Array{Float64,2},
                   XY::Array{Float64,2};
                   yidxs = get_yidxs(glrm.losses), # mapping from columns of A to columns of Y; by default, the identity
                   include_regularization=true)
    m,n = size(glrm.A)
    @assert(size(XY)==(m,yidxs[end][end]))
    @assert(size(Y)==(glrm.k,yidxs[end][end]))
    @assert(size(X)==(glrm.k,m))
    err = 0.0
    for j=1:n
        for i in glrm.observed_examples[j]
            err += evaluate(glrm.losses[j], XY[i,yidxs[j]], glrm.A[i,j])
        end
    end
    # add regularization penalty
    if include_regularization
        err += calc_penalty(glrm,X,Y; yidxs = yidxs)
    end
    return err
end
function row_objective(glrm::AbstractGLRM, i::Int, x::AbstractArray, Y::Array{Float64,2} = glrm.Y;
                   yidxs = get_yidxs(glrm.losses), # mapping from columns of A to columns of Y; by default, the identity
                   include_regularization=true)
    m,n = size(glrm.A)
    err = 0.0
    XY = x'*Y
    for j in glrm.observed_features[i]
        err += evaluate(glrm.losses[j], XY[1,yidxs[j]], glrm.A[i,j])
    end
    # add regularization penalty
    if include_regularization
        err += evaluate(glrm.rx[i], x)
    end
    return err
end
function col_objective(glrm::AbstractGLRM, j::Int, y::AbstractArray, X::Array{Float64,2} = glrm.X;
                   include_regularization=true)
    m,n = size(glrm.A)
    sz = size(y)
    if length(sz) == 1 colind = 1 else colind = 1:sz[2] end
    err = 0.0
    XY = X'*y
    obsex = glrm.observed_examples[j]
    @inbounds XYj = XY[obsex,colind]
    @inbounds Aj = convert(Array, glrm.A[obsex,j])
    err += evaluate(glrm.losses[j], XYj, Aj)
    # add regularization penalty
    if include_regularization
        err += evaluate(glrm.ry[j], y)
    end
    return err
end
# The user can also pass in X and Y and `objective` will compute XY for them
function objective(glrm::GLRM, X::Array{Float64,2}, Y::Array{Float64,2};
                   sparse=false, include_regularization=true,
                   yidxs = get_yidxs(glrm.losses), kwargs...)
    @assert(size(Y)==(glrm.k,yidxs[end][end]))
    @assert(size(X)==(glrm.k,size(glrm.A,1)))
    XY = Array{Float64}(undef, (size(X,2), size(Y,2)))
    if sparse
        # Calculate X'*Y only at observed entries of A
        m,n = size(glrm.A)
        err = 0.0
        for j=1:n
            for i in glrm.observed_examples[j]
                err += evaluate(glrm.losses[j], dot(X[:,i],Y[:,yidxs[j]]), glrm.A[i,j])
            end
        end
        if include_regularization
            err += calc_penalty(glrm,X,Y; yidxs = yidxs)
        end
        return err
    else
        # dense calculation variant (calculate XY up front)
        gemm!('T','N',1.0,X,Y,0.0,XY)
        return objective(glrm, X, Y, XY; include_regularization=include_regularization, yidxs = yidxs, kwargs...)
    end
end
# Or just the GLRM and `objective` will use glrm.X and .Y
objective(glrm::GLRM; kwargs...) = objective(glrm, glrm.X, glrm.Y; kwargs...)

# For shared arrays
# TODO: compute objective in parallel
objective(glrm::ShareGLRM, X::SharedArray{Float64,2}, Y::SharedArray{Float64,2}) =
    objective(glrm, X.s, Y.s)

# Helper function to calculate the regularization penalty for X and Y
function calc_penalty(glrm::AbstractGLRM, X::Array{Float64,2}, Y::Array{Float64,2};
    yidxs = get_yidxs(glrm.losses))
    m,n = size(glrm.A)
    @assert(size(Y)==(glrm.k,yidxs[end][end]))
    @assert(size(X)==(glrm.k,m))
    penalty = 0.0
    for i=1:m
        penalty += evaluate(glrm.rx[i], view(X,:,i))
    end
    for f=1:n
        penalty += evaluate(glrm.ry[f], view(Y,:,yidxs[f]))
    end
    return penalty
end

## ERROR METRIC EVALUATION (BASED ON DOMAINS OF THE DATA)
function raw_error_metric(glrm::AbstractGLRM, XY::Array{Float64,2}, domains::Array{Domain,1};
    yidxs = get_yidxs(glrm.losses))
    m,n = size(glrm.A)
    err = 0.0
    for j=1:n
        for i in glrm.observed_examples[j]
            err += error_metric(domains[j], glrm.losses[j], XY[i,yidxs[j]], glrm.A[i,j])
        end
    end
    return err
end
function std_error_metric(glrm::AbstractGLRM, XY::Array{Float64,2}, domains::Array{Domain,1};
    yidxs = get_yidxs(glrm.losses))
    m,n = size(glrm.A)
    err = 0.0
    for j=1:n
        column_mean = 0.0
        column_err = 0.0
        for i in glrm.observed_examples[j]
            column_mean += glrm.A[i,j]^2
            column_err += error_metric(domains[j], glrm.losses[j], XY[i,yidxs[j]], glrm.A[i,j])
        end
        column_mean = column_mean/length(glrm.observed_examples[j])
        if column_mean != 0
            column_err = column_err/column_mean
        end
        err += column_err
    end
    return err
end
function error_metric(glrm::AbstractGLRM, XY::Array{Float64,2}, domains::Array{Domain,1};
    standardize=false,
    yidxs = get_yidxs(glrm.losses))
    m,n = size(glrm.A)
    @assert(size(XY)==(m,yidxs[end][end]))
    if standardize
        return std_error_metric(glrm, XY, domains; yidxs = yidxs)
    else
        return raw_error_metric(glrm, XY, domains; yidxs = yidxs)
    end
end
# The user can also pass in X and Y and `error_metric` will compute XY for them
function error_metric(glrm::AbstractGLRM, X::Array{Float64,2}, Y::Array{Float64,2}, domains::Array{Domain,1}=Domain[l.domain for l in glrm.losses]; kwargs...)
    XY = Array{Float64}(undef,(size(X,2), size(Y,2)))
    gemm!('T','N',1.0,X,Y,0.0,XY)
    error_metric(glrm, XY, domains; kwargs...)
end
# Or just the GLRM and `error_metric` will use glrm.X and .Y
error_metric(glrm::AbstractGLRM, domains::Array{Domain,1}; kwargs...) = error_metric(glrm, glrm.X, glrm.Y, domains; kwargs...)
error_metric(glrm::AbstractGLRM; kwargs...) = error_metric(glrm, Domain[l.domain for l in glrm.losses]; kwargs...)

# Use impute and errors over GLRMS
impute(glrm::AbstractGLRM) = impute(glrm.losses, glrm.X'*glrm.Y)
function impute_missing(glrm::AbstractGLRM)
  Ahat = impute(glrm)
  for j in 1:size(glrm.A,2)
    for i in glrm.observed_examples[j]
      Ahat[i,j] = glrm.A[i,j]
    end
  end
  return Ahat
end


text:
# playMadoka.jl

#global flog=open("log.playMadoka","a+")

##文法
include("grammer.jl")

qvar=[誰 何 何処]

#魔法世界
const QB=:QB
const ほむら=:ほむら
const 暁美=:暁美
const まどか=:まどか
const さやか=:さやか
const 上条=:上条
const 杏子=:杏子
const まみ=:まみ
const ひとみ=:ひとみ
const 魔法少女=:魔法少女
const 魔女=:魔女
const ソウルジェム=:ソウルジェム
const ワルプルギスの夜=:ワルプルギスの夜
const 濁っている=:濁っている
const きれい=:きれい

魔法会議=[
 [暁美 は まどか が 好き],
 [杏子 は さやか が 好き],
 [さやか は 上条 が 好き],
 [ひとみ は 上条 が 好き],
 [上条 は ひとみ が 好き],
 [まどか は 魔法少女 です],
 [暁美 は 魔法少女 です],
 [杏子 は 魔法少女 です],
 [さやか は 魔法少女 です],
 [まみ は 魔法少女 です],
 [さやか は 魔女 です],
 [暁美 は ソウルジェム を 持っている],
 [さやか は ソウルジェム を 持っている],
 [まみ は ソウルジェム を 持っている],
 [QB は 魔法少女 を 作る],
 [魔法少女 が 魔女 に なる],
 [ソウルジェム は 濁っている]
]

魔法会議_0=[
 [暁美 の ソウルジェム は きれい],
 [さやか は ソウルジェム を 持っている],
 [まみ は ソウルジェム を 持っている],
 [ソウルジェム は きれい]
]
魔法会議_1=[
 [暁美 は ソウルジェム を 持っている],
 [さやか は ソウルジェム を 持っている],
 [まみ は ソウルジェム を 持っている],
 [QB は 魔法少女 を 作る],
 [魔法少女 が 魔女 に なる],
 [ソウルジェム は 濁っている]
]
魔法会議_2=[
 [魔法少女 が 魔女 に なる],
 [ソウルジェム は 濁っている]
]

qm1=[[誰 は ソウルジェム を 持っている]]
qm2=[[誰 は 誰 が 好き]]
qm3=[[誰 は 魔法少女 を 作る]]

m1 = solveq([誰], [誰 は 魔法少女 です],魔法会議)
st1 = statements([誰], [[誰 は 魔法少女 です]],m1)

print(stringT(st1))

solveq([:x,:y],[:x は :y です],魔法会議)
solveq([:x,:y],[:x は 魔法少女 です],魔法会議)
solveQ([:x,:y,:z],[[:x は :y を 作る],[:z は :y です]],魔法会議)

solveQ([:x,:y,:z],[[:x は :y を 作る],[:z は :y です]],魔法会議)

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews