text:
using Test
using SparseArrays

@testset "Model Generation" begin
    G = [0.0 0.0 1.0 1.0;  0.0 1.0 0.0 1.0]
    T1 = SparseArrays.sparse(Array{Int8,2}([1 1 0 0;  0 1 1 0;  0 0 1 1]))
    T2 = SparseArrays.sparse(Array{Int8,2}([1 1 1]))

    @test isnothing(CAGD.Model())

    m = CAGD.Model(G)
    @test size(m, 0, 1) == 2
    @test size(m, 0, 2) == 4
    @test size(m, 1, 1) == 0
    @test size(m, 1, 2) == 4
    @test size(m, 2, 1) == 0
    @test size(m, 2, 2) == 0
    @test_throws BoundsError size(m, 3, 1)

    m = CAGD.Model(G, [T1, T2])
    @test size(m, 0, 1) == 2
    @test size(m, 0, 2) == 4
    @test size(m, 1, 1) == 3
    @test size(m, 1, 2) == 4
    @test size(m, 2, 1) == 1
    @test size(m, 2, 2) == 3
    @test_throws BoundsError size(m, 3, 1)

    @test_throws ArgumentError CAGD.Model(G, [T1])
    @test_throws ArgumentError CAGD.Model(G, [T1, T1])
    @test_throws ArgumentError CAGD.Model(G, [T1, T1, T1])

end

@testset "Model Cloning" begin
    G = [0.0 0.0 1.0 1.0;  0.0 1.0 0.0 1.0]
    T1 = SparseArrays.sparse(Array{Int8,2}([1 1 0 0;  0 1 1 0;  0 0 1 1]))
    T2 = SparseArrays.sparse(Array{Int8,2}([1 1 1]))
    m1 = CAGD.Model(G, [T1, T2])
    m2 = deepcopy(m1)
    m3 = copy(m2)
    CAGD.addModelCell!(m2, 1, [1, 3])

    @test m1.T[1] != m2.T[1]
    @test m2.T[1] == m3.T[1]
end

@testset "Topology Manipulation" begin

    G = [
        0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0
        0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0
        0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0
    ]
    m = CAGD.Model(G)

    @testset "add Model Cells" begin
        CAGD.addModelCell!(m, 1, SparseArrays.sparse(Int8[-1; 1; 0; 0;  0; 0; 0; 0]))
        @test_throws MethodError CAGD.addModelCell!(m, 1, SparseArrays.sparse(Int8[-1 0 1 0  0 0 0 0]))
        @test_throws ArgumentError CAGD.addModelCell!(m, 1, SparseArrays.sparse(Int8[-1; 1; 0; 0;  0; 0; 0]))

        CAGD.addModelCells!(m, 1, SparseArrays.sparse(Int8[-1 0 1 0  0 0 0 0]))
        @test_throws ArgumentError CAGD.addModelCells!(m, 1, SparseArrays.sparse(Int8[-1 0 1 0  0 0 0]))
        @test size(m, 1) == (2, 8)

        @test_throws ArgumentError CAGD.addModelCells!(m, 0, SparseArrays.sparse(Int8[-1 0 1 0  0 0 0 0]))
        @test_throws ArgumentError CAGD.addModelCells!(m, 4, SparseArrays.sparse(Int8[-1 0 1 0  0 0 0 0]))
    end

    @testset "add Model Cells (EV)" begin
        CAGD.addModelCell!(m, 1, [2; 3], signed = true)
        @test m.T[1][3, :].nzval == [-1; 1]
        @test_throws ArgumentError CAGD.addModelCell!(m, 2, [1, 2, 4])
        CAGD.addModelCells!(m, 2, [[1, 2, 3]], signed = true)
        @test Array(m.T[2][1, :]) == [1; -1; 1]
        CAGD.addModelCells!(m, 1, [[2, 4], [3, 4]], signed = true)
        @test size(m, 2) == (1, 5)
        CAGD.addModelCell!(m, 2, [1; 2; 4; 5])
        @test m.T[2][2, :].nzval == [1; 1; 1; 1]
        CAGD.addModelCells!(m, 2, [[3; 4; 5]], signed = true)
        @test size(m, 3) == (0, 3)
    end

    @testset "get Model Lo Cell" begin
        @test CAGD.getModelLoCell(m, 2, 1) == [1; 2; 3]
        @test CAGD.getModelLoCell(m, 2, 2) == [1; 2; 4; 5]
        @test CAGD.getModelLoCell(m, 1, 2) == [1; 3]
        @test_throws ArgumentError CAGD.getModelLoCell(m, 0, 2)
        @test_throws ArgumentError CAGD.getModelLoCell(m, 4, 2)
        @test_throws ArgumentError CAGD.getModelLoCell(m, 3, 2)
        @test_throws ArgumentError CAGD.getModelLoCell(m, 2, -1)

        @test CAGD.getModelLoCell(m, 2, [1;2]) == [1; 2; 3; 4; 5]
        @test CAGD.getModelLoCell(m, 1, [1;2]) == [1; 2; 3]
        @test_throws ArgumentError CAGD.getModelLoCell(m, 0, [1; 2])
        @test_throws ArgumentError CAGD.getModelLoCell(m, 4, [1; 2])
        @test_throws ArgumentError CAGD.getModelLoCell(m, 2, [6; 2])
        @test_throws ArgumentError CAGD.getModelLoCell(m, 2, [2; 6])
        @test_throws ArgumentError CAGD.getModelLoCell(m, 2, [-1; 1])
    end

    @testset "get Model Cell Vertices" begin
        @test CAGD.getModelCellVertices(m, 2, 3) == [
            [0.0, 0.0, 1.0],
            [0.0, 1.0, 0.0],
            [0.0, 1.0, 1.0]
        ]

        @test CAGD.getModelCellVertices(m, 2, [1; 3]) == [
            [0.0, 0.0, 0.0],
            [0.0, 0.0, 1.0],
            [0.0, 1.0, 0.0],
            [0.0, 1.0, 1.0]
        ]

        @test CAGD.getModelCellVertices(m, 2, [3], true)[2] == [2; 3; 4]

        @test_throws ArgumentError CAGD.getModelCellVertices(m, 0, [1; 2])
        @test_throws ArgumentError CAGD.getModelCellVertices(m, 4, [1; 2])
        @test_throws ArgumentError CAGD.getModelCellVertices(m, 2, [6; 2])
        @test_throws ArgumentError CAGD.getModelCellVertices(m, 2, [2; 6])
        @test_throws ArgumentError CAGD.getModelCellVertices(m, 2, [-1; 1])
    end

    @testset "get Model Cell Geometry" begin
        @test CAGD.getModelCellGeometry(m, 2, 3) == [
            0.0 0.0 0.0
            0.0 1.0 1.0
            1.0 0.0 1.0
        ]

        @test CAGD.getModelCellGeometry(m, 2, [1; 3]) == [
            0.0 0.0 0.0 0.0
            0.0 0.0 1.0 1.0
            0.0 1.0 0.0 1.0
        ]

        @test CAGD.getModelCellGeometry(m, 2, [3], true)[2] == [2; 3; 4]

        @test_throws ArgumentError CAGD.getModelCellGeometry(m, 0, [1; 2])
        @test_throws ArgumentError CAGD.getModelCellGeometry(m, 4, [1; 2])
        @test_throws ArgumentError CAGD.getModelCellGeometry(m, 2, [6; 2])
        @test_throws ArgumentError CAGD.getModelCellGeometry(m, 2, [2; 6])
        @test_throws ArgumentError CAGD.getModelCellGeometry(m, 2, [-1; 1])
    end

    @testset "delete Model Cells" begin
        CAGD.deleteModelCell!(m, 2, 2)
        @test size(m, 2) == (2, 5)
        CAGD.addModelCell!(m, 2, [1; 2; 4; 5])
        CAGD.deleteModelCells!(m, 1, [4; 5])
        @test size(m, 2) == (1, 3)
        @test_throws ArgumentError CAGD.deleteModelCells!(m, 1, [4; 5])
        @test_throws ArgumentError CAGD.deleteModelCells!(m, 0, [4; 5])
        @test_throws ArgumentError CAGD.deleteModelCells!(m, 4, [4; 5])
    end

end

@testset "Geometry Manipulation" begin
    G1 = [0. 1.; 0. 1.];
    G2 = [1. 0.; 0. 1.];
    m = CAGD.Model(G1);

    @testset "add model vertices" begin
        CAGD.addModelCell!(m, 1, [1, 2]);
        CAGD.addModelVertices!(m, G2)
        @test m.G == [G1 G2]
        @test size(m, 1) == (1, 4)
        @test m.T[1] == SparseArrays.sparse(Int8[1 1 0 0])
    end

    @testset "delete model vertices" begin
        CAGD.deleteModelVertex!(m, 3)
        @test m.G == [0. 1. 0.; 0. 1. 1.];
        @test size(m, 1) == (1, 3)

        CAGD.deleteModelVertices!(m, [1, 3]);
        @test m.G == [1.; 1.][:, :]
        @test size(m, 1) == (0, 1)

        @test_throws ArgumentError CAGD.deleteModelVertices!(m, [1, 3])
        @test_throws ArgumentError CAGD.deleteModelVertices!(m, [0, 1])

    end
end

@testset "Model Manipulation" begin

    @testset "Purge Model" begin
        m = CAGD.Model([0.0 0.0 1.0 0.0 2.0; 0.0 1.0 0.0 2.0 0.0], [
            SparseArrays.sparse(Int8[-1 1 0 0 0; -1 0 1 0 0; 0 -1 1 0 0; 0 -1 0 1 0]),
            SparseArrays.sparse(Int8[1 -1 1 0 ])
        ])
        CAGD.modelPurge!(m)
        @test size(m, 1) == (4, 4)

        CAGD.modelPurge!(m, depth = 1)
        @test size(m, 1) == (3, 3)
    end

    @testset "Merge Model Vertices" begin
        m = CAGD.Model([
                0.0  0.0  1.0  0.0   0.0  -1.0  0.0  -1.0  -1.0
                0.0  1.0  0.0  0.0  -1.0   0.0  1.0   0.0   1.0
            ], [
                SparseArrays.sparse(Int8[
                    -1   1  0   0   0  0   0   0  0
                    -1   0  1   0   0  0   0   0  0
                     0  -1  1   0   0  0   0   0  0
                     0   0  0  -1   1  0   0   0  0
                     0   0  0  -1   0  1   0   0  0
                     0   0  0   0  -1  1   0   0  0
                     0   0  0   0   0  0  -1   1  0
                     0   0  0   0   0  0  -1   0  1
                     0   0  0   0   0  0   0  -1  1
                ]),
                SparseArrays.sparse(Int8[
                    1  -1  1  0   0  0  0   0  0
                    0   0  0  1  -1  1  0   0  0
                    0   0  0  0   0  0  1  -1  1
                ]),
        ])

        m1 = CAGD.mergeModelVertices(m, signed_merge = false);
        m2 = CAGD.mergeModelVertices(m, signed_merge = true);

        @test size(m1, 1) == (9, 6)
        @test size(m2, 1) == (9, 6)
        @test size(m1, 2) == (3, 9)
        @test size(m2, 2) == (3, 9)

        @test m1.G == m2.G
        @test sortslices(m1.T[1], dims = 1) == sortslices(abs.(m2.T[1]), dims = 1)

        @test m1.G == [
            0.0  0.0  1.0   0.0  -1.0  -1.0
            0.0  1.0  0.0  -1.0   0.0   1.0
        ]

        @test Matrix(m2.T[1]) == [
            -1   0  0   0   1  0
            -1   0  0   1   0  0
            -1   0  1   0   0  0
            -1   1  0   0   0  0
             0  -1  0   0   0  1
             0  -1  0   0   1  0
             0  -1  1   0   0  0
             0   0  0  -1   1  0
             0   0  0   0  -1  1
        ]

        @test Matrix(m2.T[2]) == [
            1  -1  0   0  0   0   0  -1   0
            0   0  1  -1  0   0  -1   0   0
            0   0  0   0  1  -1   0   0  -1
        ]
    end
end

@testset "Models Interaction" begin
    @testset "Unite Models" begin
        m1 = CAGD.Model([0.0 0.0 1.0; 0.0 1.0 0.0], [
            SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
            SparseArrays.sparse(Int8[1 -1 1])
        ])

        m2 = CAGD.Model([0.0 0.0 -1.0; 0.0 -1.0 0.0], [
            SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
            SparseArrays.sparse(Int8[1 -1 1])
        ])

        m3 = CAGD.Model([0.0 -1.0 -1.0; 1.0 0.0 1.0], [
            SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
            SparseArrays.sparse(Int8[1 -1 1])
        ])

        m4 = CAGD.uniteModels(m1, m2)
        CAGD.uniteModels!(m1, m2)

        @test m1 == m4
        @test size(m1, 1) == (6, 6)
        @test size(m1, 2) == (2, 6)

        models = [
            CAGD.Model([0.0 0.0 1.0; 0.0 1.0 0.0], [
                SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
                SparseArrays.sparse(Int8[1 -1 1])
            ]);

            CAGD.Model([0.0 0.0 -1.0; 0.0 -1.0 0.0], [
                SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
                SparseArrays.sparse(Int8[1 -1 1])
            ]);

            CAGD.Model([0.0 -1.0 -1.0; 1.0 0.0 1.0], [
                SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
                SparseArrays.sparse(Int8[1 -1 1])
            ])
        ];

        CAGD.uniteModels!(m1, m3)

        @test m1 == CAGD.uniteMultipleModels(models)
    end

    @testset "Merge Models" begin
        models = [
            CAGD.Model([0.0 0.0 1.0; 0.0 1.0 0.0], [
                SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
                SparseArrays.sparse(Int8[1 -1 1])
            ]);

            CAGD.Model([0.0 0.0 -1.0; 0.0 -1.0 0.0], [
                SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
                SparseArrays.sparse(Int8[1 -1 1])
            ]);

            CAGD.Model([0.0 -1.0 -1.0; 1.0 0.0 1.0], [
                SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
                SparseArrays.sparse(Int8[1 -1 1])
            ])
        ]

        m = CAGD.mergeMultipleModels(models, signed_merge = true)

        @test size(m, 1) == (9, 6)
        @test size(m, 2) == (3, 9)

        @test m.G == [
            0.0  0.0  1.0   0.0  -1.0  -1.0
            0.0  1.0  0.0  -1.0   0.0   1.0
        ]

        @test Matrix(m.T[1]) == [
            -1   0  0   0   1  0
            -1   0  0   1   0  0
            -1   0  1   0   0  0
            -1   1  0   0   0  0
             0  -1  0   0   0  1
             0  -1  0   0   1  0
             0  -1  1   0   0  0
             0   0  0  -1   1  0
             0   0  0   0  -1  1
        ]

        @test Matrix(m.T[2]) == [
            1  -1  0   0  0   0   0  -1   0
            0   0  1  -1  0   0  -1   0   0
            0   0  0   0  1  -1   0   0  -1
        ]
    end
end

@testset "Model Export" begin
    modelOut = CAGD.Model([0.0 0.0 1.0; 0.0 1.0 0.0], [
        SparseArrays.sparse(Int8[-1 1 0; -1 0 1; 0 -1 1]),
        SparseArrays.sparse(Int8[1 -1 1])
    ]);

    #CAGD.jlexportModel(modelOut, "./testModelExport.jl", "modelIn");

    #include("./testModelExport.jl");

    #@test modelIn == modelOut

    #rm("./testModelExport.jl")
end


text:
function discretizebranches!(system::CVSystem,old=Dict("a"=>0),restart="no")
    h = [];

    # branch grid spacing
    for i in 1:length(system.branches.ID)
        push!(system.branches.k,system.branches.lengthincm[i]*cmTom/
            (system.solverparams.JL-1));
        push!(h,system.solverparams.CFL*system.branches.k[i]/
            system.branches.c0[i][end]);
    end

    # time step guaranteed to satisfy CFL for all branches
    # system.solverparams.h = minimum(h);

    # fix time step size to couple with 3D liver model
    system.solverparams.h = 1e-4;

    if restart == "no"
        system.solverparams.numsteps = ceil(system.heart.activation.th[1]/
            system.solverparams.h);

        # allocate space for 1D domain solution variables
        for i in 1:length(system.branches.ID)
            push!(system.branches.A,
                zeros(system.solverparams.JL,system.solverparams.numsteps+1));
            push!(system.branches.Q,
                zeros(system.solverparams.JL,system.solverparams.numsteps+1));
            push!(system.branches.P,
                zeros(system.solverparams.JL,system.solverparams.numsteps+1));
            push!(system.branches.Fp,
                zeros(2*system.solverparams.JL));
            push!(system.branches.Fbarforward,
                zeros(2*system.solverparams.JL-4));
            push!(system.branches.Fbarbackward,
                zeros(2*system.solverparams.JL-4));
            push!(system.branches.Abackward,
                zeros(system.solverparams.JL-2));
            push!(system.branches.Aforward,
                zeros(system.solverparams.JL-2));
            push!(system.branches.Qbackward,
                zeros(system.solverparams.JL-2));
            push!(system.branches.Qforward,
                zeros(system.solverparams.JL-2));
            push!(system.branches.W1end,0.);
            push!(system.branches.W1,0.);
            push!(system.branches.W2,0.);
        end

        # discretize time for first cardiac cycle
        system.t = system.solverparams.h*[0:1:size(system.branches.A[1],2)-1;];
    elseif restart == "yes"
        # determine time shift
        temp = old["heart"];
        temp = temp["activation"];
        system.solverparams.tshift = old["t"][end] - sum(temp["th"][1:end-1]);

        # change grid spacing (e.g., for grid refinement study)
        JLnew = 11;

        if JLnew != system.solverparams.JL
            system.solverparams.JL = JLnew;
            for i = 1:length(system.branches.ID)
                system.branches.k[i] = system.branches.lengthincm[i]*cmTom/
                    (system.solverparams.JL-1);
                h[i] = system.solverparams.CFL*system.branches.k[i]/
                    system.branches.c0[i][end];
            end

            # fix other index ranges
            system.solverparams.acols = [1:system.solverparams.JL;];
            system.solverparams.qcols = [system.solverparams.JL+1:2*system.solverparams.JL;];
            system.solverparams.acolspre = [2:system.solverparams.JL-1;];
            system.solverparams.qcolspre = [system.solverparams.JL+2:2*system.solverparams.JL-1;];
            system.solverparams.acolscor = [1:system.solverparams.JL-2;];
            system.solverparams.qcolscor = [system.solverparams.JL-1:2*system.solverparams.JL-4;];
            system.solverparams.colsint = [2:system.solverparams.JL-1;];

            # time step guaranteed to satisfy CFL for all branches
            system.solverparams.h = minimum(h);
        end

        # update total number of time steps
        ntoadd = ceil((system.heart.activation.th[end]-
            system.solverparams.tshift)/system.solverparams.h);
        system.solverparams.numsteps=ntoadd;

        # update discrete times
        ttoadd = [0:1:ntoadd;]*system.solverparams.h + system.solverparams.tshift;
        system.t = ttoadd;

        # allocate space for 1D domain solution variables
        for i in 1:length(system.branches.ID)
            push!(system.branches.A,
                zeros(system.solverparams.JL,system.solverparams.numsteps+1));
            push!(system.branches.Q,
                zeros(system.solverparams.JL,system.solverparams.numsteps+1));
            push!(system.branches.P,
                zeros(system.solverparams.JL,system.solverparams.numsteps+1));
            push!(system.branches.Fp,
                zeros(2*system.solverparams.JL));
            push!(system.branches.Fbarforward,
                zeros(2*system.solverparams.JL-4));
            push!(system.branches.Fbarbackward,
                zeros(2*system.solverparams.JL-4));
            push!(system.branches.Abackward,
                zeros(system.solverparams.JL-2));
            push!(system.branches.Aforward,
                zeros(system.solverparams.JL-2));
            push!(system.branches.Qbackward,
                zeros(system.solverparams.JL-2));
            push!(system.branches.Qforward,
                zeros(system.solverparams.JL-2));
            push!(system.branches.W1end,0.);
            push!(system.branches.W1,0.);
            push!(system.branches.W2,0.);
        end
    end
    print("Number of nodes per branch: ")
    println(system.solverparams.JL)
    return system
end

ext: jl

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
"""
This submodule contains example systems
"""
module PBTExample

using ..ProBeTune
using LightGraphs

relu(x) = max.(0., x)

"""
Random fourier modes sampler
"""
function random_fourier_modes_sampler(rng, t_span, n; N_frequencies = 10, periods = 2)
    a = randn(rng, N_frequencies)
    theta = 2*pi*rand(rng, N_frequencies)
    ω = periods * 2π/(t_span[2] - t_span[1])
    return t -> sum([a[n]*cos(n*t+theta[n]) for n in 1:N_frequencies])
end

struct nl_diff_dyn{T}
    L::T
end

function (dd::nl_diff_dyn)(dx, x, i, p, t)
    flows = dd.L * x
    @. dx = x - relu(p) * x^3 - flows
    dx[1] += i - x[1]
    nothing
end

mutable struct StandardOutputMetric
    n_sys
    n_spec
end

function (som::StandardOutputMetric)(sol_sys, sol_spec)
    if sol_sys.retcode == :Success && sol_spec.retcode == :Success
        return sum((sol_sys[som.n_sys, :] .- sol_spec[som.n_spec, :]) .^ 2)
    else
        return Inf # Solvers failing is bad.
    end
end

mutable struct NoTransientOutputMetric
    n_sys
    n_spec
    transient_time
end

function (som::NoTransientOutputMetric)(sol_sys, sol_spec)
    if sol_sys.retcode == :Success && sol_spec.retcode == :Success
        return sum((sol_sys[som.n_sys, som.transient_time:end] .- sol_spec[som.n_spec, som.transient_time:end]) .^ 2)
    else
        return Inf # Solvers failing is bad.
    end
end

function create_diffusive_nl_example(dim_sys, av_deg, tsteps, N_samples)
    g_spec = SimpleGraph([Edge(1 => 2)])
    g_sys = barabasi_albert(dim_sys, av_deg)

    f_spec = nl_diff_dyn(laplacian_matrix(g_spec))
    f_sys = nl_diff_dyn(laplacian_matrix(g_sys))

    t_span = (tsteps[1], tsteps[end])
    PBTProblem(
        f_spec = f_spec,
        f_sys = f_sys,
        input_sampler = ContinousInputSampler(random_fourier_modes_sampler, t_span),
        output_metric = StandardOutputMetric(1, 1),
        N_samples = N_samples,
        size_p_spec = 2,
        size_p_sys = dim_sys,
        y0_spec = zeros(2),
        y0_sys = zeros(dim_sys),
        tsteps = tsteps,
        t_span = t_span
    )
end

mutable struct kuramoto_osc{w, N, K}
    w::w
    N::N
    K_av::K
end

function (dd::kuramoto_osc)(dx, x, i, p, t)
    # x -> Theta, p -> K
    K_total = 0.
    p_ind = 1
    for k in 1:dd.N
        dx[k] = x[k + dd.N]
    end
    for k in 1:dd.N
        for j in 1:(k-1)
            dx[k+dd.N] -= (relu(p[p_ind]) + 1.) * sin(x[k] - x[j])
            dx[j+dd.N] -= (relu(p[p_ind]) + 1.) * sin(x[j] - x[k])
            K_total += relu(p[p_ind]) + 1.
            p_ind += 1
        end
    end
    for k in 1:dd.N
      dx[k+dd.N] *= dd.K_av / K_total * (p_ind - 1) # (p_ind - 1) counts the number of edges
      dx[k+dd.N] += dd.w[k] - (relu(p[p_ind]) + 1.) * x[k+dd.N]
      p_ind += 1
    end
    # dx[1+dd.N] += dd.K_av * i
    dx[1+dd.N] += 1. * i
    nothing
end
## The specification, a kuramoto with inertia

function kuramoto_spec(dx, x, i, p, t)
    dx[1] = x[2]
    dx[2] = p[1] - (relu(p[2]) + 1.) * x[2]  + 1. * i# inertial node with a slackbus
    nothing
end

function kuramoto_out_metric(sol_sys, sol_spec; n_transient = length(0.:0.1:2pi))
    if sol_sys.retcode == :Success && sol_spec.retcode == :Success
        return sum((sol_sys[1, n_transient:end] .- sol_spec[1, n_transient:end]) .^ 2)
    else
        return Inf # Solvers failing is bad.
    end
end

function create_kuramoto_example(w, N_osc, dim_p_spec, K,  tsteps, N_samples)
    f_sys = kuramoto_osc(w[1:N_osc], N_osc, K)
    t_span = (tsteps[1], tsteps[end])
    PBTProblem(
        f_spec = kuramoto_spec,
        f_sys = f_sys,
        input_sampler = ContinousInputSampler(random_fourier_modes_sampler, t_span),
        output_metric = kuramoto_out_metric, # phase at interface node
        N_samples = N_samples,
        size_p_spec = dim_p_spec, # Parameters in the spec
        size_p_sys = N_osc * (N_osc - 1) ÷ 2 + N_osc, # Parameters in the sys
        y0_spec = zeros(2), # Initial state conditions for spec ...
        y0_sys = zeros(2 * N_osc), # ... and sys
        tsteps = tsteps,
        t_span = t_span)
end

end


text:
#######################################
#
# elp file
#
#######################################

@doc """
Read elp file

#### Input
* `fname`: Name or path for the sfp file
* `coordinate`: Coordinate system for electrode location
* `r`: Radius for converting spherical coords

#### Output
* `elecs`: Array of electrode objects
""" ->
function read_elp(fname::AbstractString; coordinate=Talairach, r::Real=90)
    # This does not work yet, need to convert to 3d coord system

    info("Reading elp file = $fname")

    # Create an empty electrode set
    elecs = Electrode[]

    # Read file
    df = readtable(fname, header = false, separator = ' ')

    phi = df[:x2]
    theta = df[:x3]

    # Save locations
    x = r .* sin(phi*(pi/180)) .* cos(theta*(pi/180))
    y = r .* sin(phi*(pi/180)) .* sin(theta*(pi/180)) - 17.5
    z = r .* cos(phi*(pi/180))

    # Convert label to ascii and remove '
    labels = df[:x1]
    for i = 1:length(labels)
        push!(elecs, Electrode(replace(labels[i], "'", "" ), coordinate(x[i], y[i], z[i]), Dict()))
    end

    debug("Imported $(length(elecs)) electrodes")

    return elecs
end


text:
import Base.full
#=
full : return A in dense String matrix form
=#
function full(A::Assoc)

    output = full(A.A)
    if !(A.val ==[1.0]) #Require Val mapping
        output = Array{AbstractString,2}([
            if A.A[x,y] != 0
                string(A.val[A.A[x,y]])
                else
                ""
                end
                for x = 1:length(A.row), y = 1:length(A.col)
        ])
    end
    output = vcat(A.col',output)
    output = hcat(vcat("",A.row),output)
    return output
end


text:
using Base.Test, DataFrames, MixedModels

contra = readtable(joinpath(dirname(@__FILE__), "data", "Contraception.csv.gz"))
for c in [:district, :use, :urban, :livch, :urbdist]
    contra[c] = pool(contra[c])
end
contra[:age2] = abs2(contra[:age])
gm1 = fit!(glmm(use01 ~ 1 + age + age2 + urban + livch + (1 | urbdist), contra, Bernoulli()));
@test lowerbd(gm1) == push!(fill(-Inf, 7), 0.)
@test isapprox(LaplaceDeviance(gm1), 2361.5457555; atol = 0.0001)
@test isapprox(loglikelihood(gm1), -1180.7729; atol = 0.001)
# There may be multiple optima here.
#@show logdet(gm1), sumabs2(gm1.u[1]), sum(gm1.devresid)
#@test isapprox(logdet(gm1), 75.717; atol = 0.001)
#@test isapprox(sumabs2(gm1.u[1]), 48.475; atol = 0.001)
#@test isapprox(sum(gm1.devresid), 2237.354; atol = 0.001)
#@test isapprox(fixef(gm1), [-1.0592,0.00328263,-0.00447496,0.770853,0.834578,0.913848,0.927232];
#    atol = 0.001)
show(IOBuffer(), gm1)

cbpp = readtable(joinpath(dirname(@__FILE__), "data", "cbpp.csv.gz"))
for c in [:herd, :period]
    cbpp[c] = pool(oftype(Int8[], cbpp[c]))
end
for c in [:incidence, :size]
    cbpp[c] = oftype(Int8[], cbpp[c])
end
cbpp[:prop] = cbpp[:incidence] ./ cbpp[:size]
gm2 = fit!(glmm(prop ~ 1 + period + (1 | herd), cbpp, Binomial(), LogitLink(); wt = cbpp[:size]));

@test isapprox(LaplaceDeviance(gm2), 100.095856; atol = 0.0001)
@test isapprox(sumabs2(gm2.u[1]), 9.72305; atol = 0.0001)
@test isapprox(logdet(gm2), 16.901054; atol = 0.0001)
@test isapprox(sum(gm2.devresid), 73.47143; atol = 0.001)
@test isapprox(loglikelihood(gm2), -92.02628; atol = 0.001)
@test sdest(gm2) == 1
@test varest(gm2) == 1
#VerbAgg = rcopy("lme4::VerbAgg")
#VerbAgg[:r201] = [Int8(x == "N" ? 0 : 1) for x in VerbAgg[:r2]]
#m3 = glmm(r201 ~ 1 + Anger + Gender + btype + situ + (1 | id) + (1 | item), VerbAgg, Bernoulli());
#f = r201 ~ 1 + Anger + Gender + btype + situ + (1 | id)


text:
s = open("input.txt") do file
    readlines(file);
end

function ParseLine(line)
    m = match(r"mem\[(\d+)\] = (\d+)", line);
    if m != nothing
        return map(x->parse(UInt64,x), m.captures);
    end

    m = match(r"mask = ([01X]+)", line);
    if m != nothing
        return (m.captures|>collect);
    end

    return ();
end

function SetAllAddresses(address, val, mask, mem)
    numX = sum(mask);
    toSet = findall(x->x==1, mask).+28;
    bAddress = bitstring(address)|>collect;

    # println(mask);
    # println(toSet)
    # println(numX);

    for i::UInt64 = 1:2^numX
        bAddress = bitstring(address)|>collect;
        bAddress[toSet] .= collect(bitstring(i))[end-numX+1:end];
        # println(String(bAddress));
        tempAddress = parse(UInt64, String(bAddress), base=2);
        mem[tempAddress] = val;
    end
    # println("Set to $(val)");
end

function SumMem(program)
    mem = Dict();
    oneMask = 0;
    xMask = -1;
    for line in program
        parsed = ParseLine(line);
        if length(parsed) == 0
            println(line);
            continue;
        end

        if length(parsed) == 1
            # new mask
            oneMask = parse(UInt64, String(map(x->ifelse(x=='1', '1', '0'), parsed[1])), base=2);
            xMask = map(x->ifelse(x=='X', 1, 0), parsed[1]|>collect);
            # println(parsed[1]);
        else
            # put in mem
            address = parsed[1];
            address |= oneMask;

            SetAllAddresses(address, parsed[2], xMask, mem);
        end
    end

    return sum(map(UInt64,values(mem)));
end

println(Int(SumMem(s[1:2])))
println(Int(SumMem(s)))


text:
export SeisEvent

@doc """
    SeisEvent

A structure for discrete seismic events, comprising three structures:
* :hdr, a SeisHdr for the event descriptor
* :source, a SeisSrc for descrition of the seismic source process
* :data, an EventTraceData structure for channel data, including phases

See also: `SeisHdr`, `SeisSrc`, `EventTraceData`
"""
mutable struct SeisEvent
  hdr::SeisHdr
  source::SeisSrc
  data::EventTraceData

  SeisEvent(hdr::SeisHdr, source::SeisSrc, data::EventTraceData) = new(hdr, source, data)
end

function SeisEvent(;
                    hdr ::SeisHdr     = SeisHdr(),
                    source::SeisSrc   = SeisSrc(),
                    data::T           = EventTraceData()
                    ) where {T<:GphysData}

  if T != EventTraceData
    return SeisEvent(hdr, source, convert(EventTraceData, data))
  else
    return SeisEvent(hdr, source, data)
  end
end

# =============================================================================
# Methods from Base
isequal(S::SeisEvent, T::SeisEvent) = min(isequal(S.hdr, T.hdr),
                                          isequal(S.source, T.source),
                                          isequal(S.data, T.data))
==(S::SeisEvent, T::SeisEvent) = isequal(S,T)
sizeof(Ev::SeisEvent) = 24 + sizeof(Ev.hdr) + sizeof(Ev.source) + sizeof(Ev.data)

# SeisEvent
write(io::IO, W::SeisEvent) = ( write(io, getfield(W, :hdr));
                                write(io, getfield(W, :source));
                                write(io, getfield(W, :data))
                                )

read(io::IO, ::Type{SeisEvent}) = SeisEvent(read(io, SeisHdr),
                                            read(io, SeisSrc),
                                            read(io, EventTraceData))

summary(V::SeisEvent) = string("Event ", V.hdr.id, ": SeisEvent with ",
  V.data.n, " channel", V.data.n == 1 ? "" : "s")

function show(io::IO, S::SeisEvent)
  println(io, summary(S))
  println(io, "\n(.hdr)")
  show(io, getfield(S, :hdr))
  println(io, "\n(.source)")
  show(io, getfield(S, :source))
  println(io, "\n(.data)")
  println(io, summary(getfield(S, :data)))
  return nothing
end
show(S::SeisEvent) = show(stdout, S)


text:
using MixedModels
using MixedModels: dataset
using MixedModelsExtras
using Statistics
using Test

function _adjust(r2, model)
    n = nobs(model)
    p = dof(model)
    return 1 - (1 - r2) * (n - 1) / (n - p)
end

@testset "LMM" begin
    model = fit(MixedModel, @formula(reaction ~ 1 + (1|subj)), dataset(:sleepstudy))

    @test r2(model;conditional=true) ≈ cor(response(model), fitted(model))^2
    @test r2(model;conditional=false) ≈ cor(response(model), model.X * model.β)^2

    @test adjr2(model;conditional=true) ≈ _adjust(r2(model; conditional=true), model)
    @test adjr2(model;conditional=false) ≈ _adjust(r2(model; conditional=false), model)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""a metric value for some object

    IoK8sApiCustomMetricsV1beta1MetricValue(;
        apiVersion=nothing,
        kind=nothing,
        describedObject=nothing,
        metricName=nothing,
        timestamp=nothing,
        windowSeconds=nothing,
        value=nothing,
    )

    - apiVersion::String : APIVersion defines the versioned schema of this representation of an object. Servers should convert recognized schemas to the latest internal value, and may reject unrecognized values. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#resources
    - kind::String : Kind is a string value representing the REST resource this object represents. Servers may infer this from the endpoint the client submits requests to. Cannot be updated. In CamelCase. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - describedObject::IoK8sApiCoreV1ObjectReference : a reference to the described object
    - metricName::String : the name of the metric
    - timestamp::IoK8sApimachineryPkgApisMetaV1Time : indicates the time at which the metrics were produced
    - windowSeconds::Int64 : indicates the window ([Timestamp-Window, Timestamp]) from which these metrics were calculated, when returning rate metrics calculated from cumulative metrics (or zero for non-calculated instantaneous metrics).
    - value::IoK8sApimachineryPkgApiResourceQuantity : the value of the metric for this
"""
mutable struct IoK8sApiCustomMetricsV1beta1MetricValue <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    describedObject::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ObjectReference } # spec name: describedObject
    metricName::Any # spec type: Union{ Nothing, String } # spec name: metricName
    timestamp::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: timestamp
    windowSeconds::Any # spec type: Union{ Nothing, Int64 } # spec name: windowSeconds
    value::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApiResourceQuantity } # spec name: value

    function IoK8sApiCustomMetricsV1beta1MetricValue(;apiVersion=nothing, kind=nothing, describedObject=nothing, metricName=nothing, timestamp=nothing, windowSeconds=nothing, value=nothing)
        o = new()
        validate_property(IoK8sApiCustomMetricsV1beta1MetricValue, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiCustomMetricsV1beta1MetricValue, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiCustomMetricsV1beta1MetricValue, Symbol("describedObject"), describedObject)
        setfield!(o, Symbol("describedObject"), describedObject)
        validate_property(IoK8sApiCustomMetricsV1beta1MetricValue, Symbol("metricName"), metricName)
        setfield!(o, Symbol("metricName"), metricName)
        validate_property(IoK8sApiCustomMetricsV1beta1MetricValue, Symbol("timestamp"), timestamp)
        setfield!(o, Symbol("timestamp"), timestamp)
        validate_property(IoK8sApiCustomMetricsV1beta1MetricValue, Symbol("windowSeconds"), windowSeconds)
        setfield!(o, Symbol("windowSeconds"), windowSeconds)
        validate_property(IoK8sApiCustomMetricsV1beta1MetricValue, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        o
    end
end # type IoK8sApiCustomMetricsV1beta1MetricValue

const _property_map_IoK8sApiCustomMetricsV1beta1MetricValue = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("describedObject")=>Symbol("describedObject"), Symbol("metricName")=>Symbol("metricName"), Symbol("timestamp")=>Symbol("timestamp"), Symbol("windowSeconds")=>Symbol("windowSeconds"), Symbol("value")=>Symbol("value"))
const _property_types_IoK8sApiCustomMetricsV1beta1MetricValue = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("describedObject")=>"IoK8sApiCoreV1ObjectReference", Symbol("metricName")=>"String", Symbol("timestamp")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("windowSeconds")=>"Int64", Symbol("value")=>"IoK8sApimachineryPkgApiResourceQuantity")
Base.propertynames(::Type{ IoK8sApiCustomMetricsV1beta1MetricValue }) = collect(keys(_property_map_IoK8sApiCustomMetricsV1beta1MetricValue))
Swagger.property_type(::Type{ IoK8sApiCustomMetricsV1beta1MetricValue }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCustomMetricsV1beta1MetricValue[name]))}
Swagger.field_name(::Type{ IoK8sApiCustomMetricsV1beta1MetricValue }, property_name::Symbol) =  _property_map_IoK8sApiCustomMetricsV1beta1MetricValue[property_name]

function check_required(o::IoK8sApiCustomMetricsV1beta1MetricValue)
    true
end

function validate_property(::Type{ IoK8sApiCustomMetricsV1beta1MetricValue }, name::Symbol, val)
end


text:
function new_variables_to_attach(
    formulater::DisjunctiveFormulater{DisjunctiveConstraints.NaiveBigM},
)
    return fill(
        ZO(),
        DisjunctiveConstraints.num_alternatives(formulater.disjunction.s),
    )
end

# NOTE: This will ignore any general constraints in `node`.
function compute_disjunction_activity(
    form::DMIPFormulation,
    formulater::DisjunctiveFormulater{DisjunctiveConstraints.NaiveBigM},
    node::Node,
    config::AlgorithmConfig,
)
    cvis = form.disjunction_formulaters[formulater]
    @assert length(cvis) ==
            DisjunctiveConstraints.num_alternatives(formulater.disjunction.s)
    lbs = Dict{CVI,Float64}()
    ubs = Dict{CVI,Float64}()
    for cvi in cvis
        @assert !haskey(lbs, cvi)
        @assert !haskey(ubs, cvi)
        l, u = get_bounds(form, cvi)
        lbs[cvi] = l
        ubs[cvi] = u
    end
    for lt_bound in node.lt_bounds
        cvi = lt_bound.cvi
        if haskey(ubs, cvi)
            ubs[cvi] = min(ubs[cvi], lt_bound.s.upper)
        end
    end
    for gt_bound in node.gt_bounds
        cvi = gt_bound.cvi
        if haskey(lbs, cvi)
            lbs[cvi] = max(lbs[cvi], gt_bound.s.lower)
        end
    end

    activity = Bool[]
    proven_active = Int[]
    ϵ_int = config.int_tol
    for (i, cvi) in enumerate(cvis)
        l = lbs[cvi]
        u = ubs[cvi]
        # Approximate version of: l == 1 == u
        if (abs(l - 1) ≤ ϵ_int) && (abs(u - 1) ≤ ϵ_int)
            push!(proven_active, i)
        end
        # Approximate version of: l <= 1 <= u
        push!(activity, (l ≤ 1 + ϵ_int) & (1 <= u + ϵ_int))
    end
    if length(proven_active) == 1
        # We have proven one alternative must be active; therefore, we can
        # infer that all others must be inactive.
        fill!(activity, false)
        activity[proven_active[1]] = true
    elseif length(proven_active) > 1
        # We have proven that two alternatives must be active. This is not
        # possible; therefore, the problem must be infeasible.
        fill!(activity, false)
    end
    return activity
end

function delete_all_constraints!(
    model::Gurobi.Optimizer,
    disjunction_state::DisjunctiveConstraints.NaiveBigMState,
)
    for lt_ci in disjunction_state.lt_cis
        if lt_ci !== nothing
            MOI.delete(model, lt_ci)
        end
    end
    for gt_ci in disjunction_state.gt_cis
        if gt_ci !== nothing
            MOI.delete(model, gt_ci)
        end
    end
    return nothing
end


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libhelicsSharedLib"], :libhelicsSharedLib),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/GMLC-TDC/HELICSBuilder/releases/download/v2.3.0"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.aarch64-linux-gnu-gcc7.tar.gz", "2479213ebdf0dec96912bb7ab67f3fc451be5c01657d5d04725e178245b813c4"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.aarch64-linux-gnu-gcc8.tar.gz", "9e71d255f654eb0a1202214268eee4effd248706287f6372daf71690e6e96ad7"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.aarch64-linux-musl-gcc7.tar.gz", "9f683ab8f01168011ed191ad34e915043048d5f1dc4cddd1d58f564523cfcd86"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.aarch64-linux-musl-gcc8.tar.gz", "3226531ddf33d7993de8a94efebf3d4b27d4cb573520b03ab7db3909dccd02e9"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.arm-linux-gnueabihf-gcc7.tar.gz", "2a24a02a238b98d9deb92d328224d1033d604ccc45207cce1425e034de203df9"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.arm-linux-gnueabihf-gcc8.tar.gz", "f00aa104df7a2519c9cf0ebc05d55c05b276e8c5a02c890c2c2696edc21dada0"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.arm-linux-musleabihf-gcc7.tar.gz", "421aaeb0267a2803ac002520ab682fe4ee9687a2a2d92e844e711d44f881c3cf"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.arm-linux-musleabihf-gcc8.tar.gz", "259d2f68689ea096c424b55b88222d96bd36aad5289ba0f9cc73706ee5d9c856"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.i686-linux-gnu-gcc7.tar.gz", "66820dd0eaacdd6e90e7c45cac5d20b6b1298bb22c6a33ca6ca635071b4f50db"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.i686-linux-gnu-gcc8.tar.gz", "9b976365906c47e47094094a373a43a54c78903e92293298d07f1e617be4ffd6"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.i686-linux-musl-gcc7.tar.gz", "49614b1a32db889bff71463e4420167b46e662749f612439cb3ced5e6d1d6ebe"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.i686-linux-musl-gcc8.tar.gz", "80c8e6c5b800bba8e6a9868db98e818895ac6a5b6738a892b3c28a3304120973"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.powerpc64le-linux-gnu-gcc7.tar.gz", "9d63e2344141fc37e007c7d63fd93c97a4c97453625ec2462ba594db1bc7f495"),
    Linux(:powerpc64le, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.powerpc64le-linux-gnu-gcc8.tar.gz", "b884bc78954ce102f8828b43c4112ced7d22fefa514a7c3f709a0036bd993103"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.x86_64-apple-darwin14-gcc7.tar.gz", "bad98d9bee1fb87b2d5df2db34ac3461e2fb555eb713e9a87f236c48fa17bd60"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.x86_64-apple-darwin14-gcc8.tar.gz", "ae0f1597e9906e8f024a4ae503c43dd24ec39ce8af25c6e449b2d65eb1d915db"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.x86_64-linux-gnu-gcc7.tar.gz", "0c2728c7a4635244b7c24ecd0cbcc6739ecdabb712307ee89a8fbb3da48276be"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.x86_64-linux-gnu-gcc8.tar.gz", "39b1ac0b0e6ff9848c8497144f2fe83f96cff0054073198e92c9d9f5237038d6"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.x86_64-linux-musl-gcc7.tar.gz", "9593958ade5e0e689fcce19cd0436e61d65f2db2cab4d40ed30d423db96ef012"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.x86_64-linux-musl-gcc8.tar.gz", "dea3206779ccb310dc4f674dba568791945083bad8ef5f9cd91ffb7bc8ccbddb"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc7)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.x86_64-w64-mingw32-gcc7.tar.gz", "f74f0fe20116b76502e50bce47885d230d1b501a4c9817653e7286b35719c8dd"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc8)) => ("$bin_prefix/libhelicsSharedLib.v2.3.0.x86_64-w64-mingw32-gcc8.tar.gz", "fb3264a7e92497f5dee0607765031672e8ae0ca5b7393729b281f83eb02893dc"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    install(dl_info...; prefix=prefix, force=true, verbose=verbose)
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products, verbose=verbose)


text:
global srvr

function start_service(s)
    global srvr
    srvr = s
    serve(srvr)
end

function stop_service()
    global srvr
    close(srvr)
end


text:
function extract_surface(_dataset::AbstractVTKStructuredData)
    if dim(_dataset) == 2
        return VTKPolyData(_dataset)
    end

    dataset = VTKStructuredData(_dataset)
    cextents = cell_extents(dataset)
    faces = Vector{Int}[]

    _cell_data = empty(_dataset.cell_data)
    for m in keys(dataset.cell_data)
        _cell_data[m] = Float64[]
    end

    for cind in surface_cell_inds(dataset)
        hexa_cc = cell_connectivity(dataset, cind)
        if cind[1] == 1
            push!(faces, [hexa_cc[1], hexa_cc[4], hexa_cc[8], hexa_cc[5]])
            for m in keys(dataset.cell_data)
                add_cell_data!(dataset, _cell_data, m, cind)
            end
        end
        if cind[1] == cextents[1]
            push!(faces, [hexa_cc[2], hexa_cc[3], hexa_cc[7], hexa_cc[6]])
            for m in keys(dataset.cell_data)
                add_cell_data!(dataset, _cell_data, m, cind)
            end
        end

        if cind[2] == 1
            push!(faces, [hexa_cc[1], hexa_cc[2], hexa_cc[6], hexa_cc[5]])
            for m in keys(dataset.cell_data)
                add_cell_data!(dataset, _cell_data, m, cind)
            end
        end
        if cind[2] == cextents[2]
            push!(faces, [hexa_cc[3], hexa_cc[4], hexa_cc[8], hexa_cc[7]])
            for m in keys(dataset.cell_data)
                add_cell_data!(dataset, _cell_data, m, cind)
            end
        end

        if cind[3] == 1
            push!(faces, [hexa_cc[1], hexa_cc[2], hexa_cc[3], hexa_cc[4]])
            for m in keys(dataset.cell_data)
                add_cell_data!(dataset, _cell_data, m, cind)
            end
        end
        if cind[3] == cextents[3]
            push!(faces, [hexa_cc[5], hexa_cc[6], hexa_cc[7], hexa_cc[8]])
            for m in keys(dataset.cell_data)
                add_cell_data!(dataset, _cell_data, m, cind)
            end
        end
    end

    _cell_types = fill(9, length(faces))

    for m in keys(dataset.cell_data)
        _var_dim = var_dim(dataset, m, "Cell")
        if _var_dim > 1
            _cell_data[m] = reshape(_cell_data[m], (_var_dim, length(faces)))
        end
    end

    _point_coords = Float64[]
    point_inds = Int[]
    _point_data = empty(_dataset.point_data)
    for m in keys(dataset.point_data)
        _point_data[m] = Float64[]
    end
    pextents = extents(dataset)

    k = 0
    for _face in faces
        for (i,p) in enumerate(_face)
            _i = findfirst(point_inds, p)
            if _i > 0
                _face[i] = _i
            else
                k += 1
                append!(_point_coords, dataset.point_coords[:,ind2sub(pextents, p)...])
                push!(point_inds, p)
                for m in keys(dataset.point_data)
                    _var_dim = var_dim(dataset, m, "Point")
                    if _var_dim == 1
                        push!(_point_data[m], dataset.point_data[m][ind2sub(pextents, p)...])
                    else
                        append!(_point_data[m], dataset.point_data[m][:, ind2sub(pextents, p)...])
                    end
                end
                _face[i] = k
            end
        end
    end

    _point_coords = reshape(_point_coords, (3, k))
    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim == 1
            _point_data[m] = reshape(_point_data[m], (k,))
        else
            _point_data[m] = reshape(_point_data[m], (_var_dim, k))
        end
    end

    return VTKPolyData(_point_coords, _cell_types, faces, _point_data, _cell_data)
    #=
    _point_coords = reshape(dataset.point_coords, (3, num_of_points(dataset)))
    _point_data = empty(_dataset.point_data)
    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim == 1
            _point_data[m] = reshape(dataset.point_data[m], (num_of_points(dataset),))
        else
            _point_data[m] = reshape(dataset.point_data[m], (_var_dim, num_of_points(dataset)))
        end
    end
    return VTKPolyData(_point_coords, _cell_types, faces, _point_data, _cell_data)
    =#
end

function add_cell_data!(dataset, _cell_data, m, cind)
    _var_dim = var_dim(dataset, m, "Cell")
    if _var_dim == 1
        _cd = dataset.cell_data[m][cind...]
        push!(_cell_data[m], _cd)
    else
        _cd = dataset.cell_data[m][:, cind...]
        append!(_cell_data[m], _cd)
    end
    return
end

function extract_surface(dataset::T) where {T<:AbstractVTKUnstructuredData}
    if T <: VTKPolyData
        return dataset
    end
    if dim(dataset) == 2
        return VTKPolyData(dataset)
    end

    #Sorted inds and cell type are keys
    #Cell connectivity, cell count, and cell data are the values

    cell_register = Dict{Tuple{Vector{Int}, Int}, Tuple{Vector{Int}, _Counter, typeof(dataset.cell_data)}}()
    for i in 1:length(dataset.cell_connectivity)
        if dataset.cell_types[i] ∉ POLY_CELLS
            _cells, _types = decompose_cell(dataset.cell_connectivity[i], dataset.cell_types[i], target="Faces")
        else
            _cells, _types = [dataset.cell_connectivity[i]], [dataset.cell_types[i]]
        end
        for j in 1:length(_cells)
            _key = (sort(_cells[j]), _types[j])
            if haskey(cell_register, _key)
                cell_register[_key][2].a += 1
            else
                cell_register[_key] = (_cells[j], _Counter(1), empty(dataset.cell_data))
                for m in keys(dataset.cell_data)
                    _var_dim = var_dim(dataset, m, "Cell")
                    if _var_dim == 1
                        cell_register[_key][3][m] = [dataset.cell_data[m][i]]
                    else
                        cell_register[_key][3][m] = dataset.cell_data[m][:,i]
                    end
                end
            end
        end
    end

    _cell_data = empty(dataset.cell_data)
    _cell_connectivity = Vector{Int}[]
    _cell_types = Int[]

    for m in keys(dataset.cell_data)
        _cell_data[m] = Float64[]
    end

    for (i, kv) in enumerate(cell_register)
        k, v = kv
        if v[2].a == 1
            push!(_cell_connectivity, v[1])
            push!(_cell_types, k[2])
            for m in keys(dataset.cell_data)
                append!(_cell_data[m], v[3][m])
            end
        end
    end

    for m in keys(dataset.cell_data)
        _var_dim = var_dim(dataset, m, "Cell")
        if _var_dim > 1
            _cell_data[m] = reshape(_cell_data[m], (_var_dim, length(_cell_connectivity)))
        end
    end

    _point_coords = Float64[]
    point_inds = Int[]
    _point_data = empty(dataset.point_data)
    for m in keys(dataset.point_data)
        _point_data[m] = Float64[]
    end

    k = 0
    for _cc in _cell_connectivity
        for (i,p) in enumerate(_cc)
            _i = findfirst(point_inds, p)
            if _i > 0
                _cc[i] = _i
            else
                k += 1
                append!(_point_coords, dataset.point_coords[:,p])
                push!(point_inds, p)
                for m in keys(dataset.point_data)
                    _var_dim = var_dim(dataset, m, "Point")
                    if _var_dim == 1
                        push!(_point_data[m], dataset.point_data[m][p])
                    else
                        append!(_point_data[m], dataset.point_data[m][:,p])
                    end
                end
                _cc[i] = k
            end
        end
    end

    _point_coords = reshape(_point_coords, (3, k))
    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim > 1
            _point_data[m] = reshape(_point_data[m], (_var_dim, k))
        end
    end

    return VTKPolyData(_point_coords, _cell_types, _cell_connectivity, _point_data, _cell_data)
end


text:
module TestFormula
    using Base.Test
    using DataFrames
    using Compat

    # TODO:
    # - grouped variables in formulas with interactions
    # - is it fast?  Can expand() handle DataFrames?
    # - deal with intercepts
    # - implement ^2 for datavector's
    # - support more transformations with I()?

    ## Formula parsing
    import DataFrames.Terms

    ## totally empty
    t = Terms(Formula(nothing, 0))
    @test t.response == false
    @test t.intercept == false
    @test t.terms == []
    @test t.eterms == []

    ## empty RHS
    t = Terms(@formula(y ~ 0))
    @test t.intercept == false
    @test t.terms == []
    @test t.eterms == [:y]
    t = Terms(@formula(y ~ -1))
    @test t.intercept == false
    @test t.terms == []

    ## intercept-only
    t = Terms(@formula(y ~ 1))
    @test t.response == true
    @test t.intercept == true
    @test t.terms == []
    @test t.eterms == [:y]

    ## terms add
    t = Terms(@formula(y ~ 1 + x1 + x2))
    @test t.intercept == true
    @test t.terms == [:x1, :x2]
    @test t.eterms == [:y, :x1, :x2]

    ## implicit intercept behavior:
    t = Terms(@formula(y ~ x1 + x2))
    @test t.intercept == true
    @test t.terms == [:x1, :x2]
    @test t.eterms == [:y, :x1, :x2]

    ## no intercept
    t = Terms(@formula(y ~ 0 + x1 + x2))
    @test t.intercept == false
    @test t.terms == [:x1, :x2]

    t = Terms(@formula(y ~ -1 + x1 + x2))
    @test t.intercept == false
    @test t.terms == [:x1, :x2]

    t = Terms(@formula(y ~ x1 & x2))
    @test t.terms == [:(x1 & x2)]
    @test t.eterms == [:y, :x1, :x2]

    ## `*` expansion
    t = Terms(@formula(y ~ x1 * x2))
    @test t.terms == [:x1, :x2, :(x1 & x2)]
    @test t.eterms == [:y, :x1, :x2]

    ## associative rule:
    ## +
    t = Terms(@formula(y ~ x1 + x2 + x3))
    @test t.terms == [:x1, :x2, :x3]

    ## &
    t = Terms(@formula(y ~ x1 & x2 & x3))
    @test t.terms == [:((&)(x1, x2, x3))]
    @test t.eterms == [:y, :x1, :x2, :x3]

    ## distributive property of + and &
    t = Terms(@formula(y ~ x1 & (x2 + x3)))
    @test t.terms == [:(x1&x2), :(x1&x3)]

    ## FAILS: ordering of expanded interaction terms is wrong
    ## (only has an observable effect when both terms are categorical and
    ## produce multiple model matrix columns that are multiplied together...)
    ##
    ## t = Terms(@formula(y ~ (x2 + x3) & x1))
    ## @test t.terms == [:(x2&x1), :(x3&x1)]

    ## three-way *
    t = Terms(@formula(y ~ x1 * x2 * x3))
    @test t.terms == [:x1, :x2, :x3,
                      :(x1&x2), :(x1&x3), :(x2&x3),
                      :((&)(x1, x2, x3))]
    @test t.eterms == [:y, :x1, :x2, :x3]

    ## Interactions with `1` reduce to main effect.  All fail at the moment.
    ## t = Terms(@formula(y ~ 1 & x1))
    ## @test t.terms == [:x1]              # == [:(1 & x1)]
    ## @test t.eterms == [:y, :x1]

    ## t = Terms(@formula(y ~ (1 + x1) & x2))
    ## @test t.terms == [:x2, :(x1&x2)]    # == [:(1 & x1)]
    ## @test t.eterms == [:y, :x1, :x2]

    ## Tests for constructing ModelFrame and ModelMatrix

    sparsetype = SparseMatrixCSC{Float64,Int}

    d = DataFrame()
    d[:y] = [1:4;]
    d[:x1] = [5:8;]
    d[:x2] = [9:12;]
    d[:x3] = [13:16;]
    d[:x4] = [17:20;]

    x1 = [5.:8;]
    x2 = [9.:12;]
    x3 = [13.:16;]
    x4 = [17.:20;]
    f = @formula(y ~ x1 + x2)
    mf = ModelFrame(f, d)
    ## @test mm.response_colnames == ["y"] # nope: no response_colnames
    @test coefnames(mf) == ["(Intercept)","x1","x2"]
    ## @test model_response(mf) == transpose([1. 2 3 4]) # fails: Int64 vs. Float64
    mm = ModelMatrix(mf)
    smm = ModelMatrix{sparsetype}(mf)
    @test mm.m[:,1] == ones(4)
    @test mm.m[:,2:3] == [x1 x2]
    @test mm.m == smm.m

    @test isa(mm.m, Matrix{Float64})
    @test isa(smm.m, sparsetype)
    @test isa(ModelMatrix{DataMatrix{Float64}}(mf).m, DataMatrix{Float64})

    #test_group("expanding a PooledVec into a design matrix of indicators for each dummy variable")

    d[:x1p] = PooledDataArray(d[:x1])
    mf = ModelFrame(@formula(y ~ x1p), d)
    mm = ModelMatrix(mf)

    @test mm.m[:,2] == [0, 1., 0, 0]
    @test mm.m[:,3] == [0, 0, 1., 0]
    @test mm.m[:,4] == [0, 0, 0, 1.]
    @test coefnames(mf)[2:end] == ["x1p: 6", "x1p: 7", "x1p: 8"]
    @test mm.m == ModelMatrix{sparsetype}(mf).m

    #test_group("create a design matrix from interactions from two DataFrames")
    ## this was removed in commit dead4562506badd7e84a2367086f5753fa49bb6a

    ## b = DataFrame()
    ## b["x2"] = DataVector(x2)
    ## df = interaction_design_matrix(a,b)
    ## @test df[:,1] == DataVector([0, 10., 0, 0])
    ## @test df[:,2] == DataVector([0, 0, 11., 0])
    ## @test df[:,3] == DataVector([0, 0, 0, 12.])

    #test_group("expanding an singleton expression/symbol into a DataFrame")
    ## generalized expand was dropped, too
    ## df = deepcopy(d)
    ## r = expand(:x2, df)
    ## @test isa(r, DataFrame)
    ## @test r[:,1] == DataVector([9,10,11,12])  # TODO: test float vs int return

    ## df = deepcopy(d)
    ## ex = :(log(x2))
    ## r = expand(ex, df)
    ## @test isa(r, DataFrame)
    ## @test r[:,1] == DataVector(log([9,10,11,12]))

    # ex = :(x1 & x2)
    # r = expand(ex, df)
    # @test isa(r, DataFrame)
    # @test ncol(r) == 1
    # @test r[:,1] == DataArray([45, 60, 77, 96])

    ## r = expand(:(x1 + x2), df)
    ## @test isa(r, DataFrame)
    ## @test ncol(r) == 2
    ## @test r[:,1] == DataVector(df["x1"])
    ## @test r[:,2] == DataVector(df["x2"])

    ## df["x1"] = PooledDataArray(x1)
    ## r = expand(:x1, df)
    ## @test isa(r, DataFrame)
    ## @test ncol(r) == 3
    ## @test r == expand(PooledDataArray(x1), "x1", DataFrame())

    ## r = expand(:(x1 + x2), df)
    ## @test isa(r, DataFrame)
    ## @test ncol(r) == 4
    ## @test r[:,1:3] == expand(PooledDataArray(x1), "x1", DataFrame())
    ## @test r[:,4] == DataVector(df["x2"])

    ## df["x2"] = PooledDataArray(x2)
    ## r = expand(:(x1 + x2), df)
    ## @test isa(r, DataFrame)
    ## @test ncol(r) == 6
    ## @test r[:,1:3] == expand(PooledDataArray(x1), "x1", DataFrame())
    ## @test r[:,4:6] == expand(PooledDataArray(x2), "x2", DataFrame())

    #test_group("Creating a model matrix using full formulas: y ~ x1 + x2, etc")

    df = deepcopy(d)
    f = @formula(y ~ x1 & x2)
    mf = ModelFrame(f, df)
    mm = ModelMatrix(mf)
    @test mm.m == [ones(4) x1.*x2]
    @test mm.m == ModelMatrix{sparsetype}(mf).m

    f = @formula(y ~ x1 * x2)
    mf = ModelFrame(f, df)
    mm = ModelMatrix(mf)
    @test mm.m == [ones(4) x1 x2 x1.*x2]
    @test mm.m == ModelMatrix{sparsetype}(mf).m

    df[:x1] = PooledDataArray(x1)
    x1e = [[0, 1, 0, 0] [0, 0, 1, 0] [0, 0, 0, 1]]
    f = @formula(y ~ x1 * x2)
    mf = ModelFrame(f, df)
    mm = ModelMatrix(mf)
    @test mm.m == [ones(4) x1e x2 [0, 10, 0, 0] [0, 0, 11, 0] [0, 0, 0, 12]]
    @test mm.m == ModelMatrix{sparsetype}(mf).m

    #test_group("Basic transformations")

    ## the log(x2) appears to be broken (again, in The Great Purge, by the
    ## removal of the x -> with(d, x) in the ModelFrame constructor)
    ## df = deepcopy(d)
    ## f = y ~ x1 + log(x2)
    ## mf = ModelFrame(f, df)
    ## mm = ModelMatrix(mf)
    ## @test mm.m == [ones(4) x1 log(x2)]

    ## df = deepcopy(d)
    ## df["x1"] = PooledDataArray([5:8])
    ## f = Formula(:(y ~ x1 * (log(x2) + x3)))
    ## mf = ModelFrame(f, df)
    ## mm = ModelMatrix(mf)
    ## @test mm.model_colnames == [
    ##  "(Intercept)"
    ##  "x1:6"
    ##  "x1:7"
    ##  "x1:8"
    ##  "log(x2)"
    ##  "x3"
    ##  "x1:6&log(x2)"
    ##  "x1:6&x3"
    ##  "x1:7&log(x2)"
    ##  "x1:7&x3"
    ##  "x1:8&log(x2)"
    ##  "x1:8&x3" ]

    #test_group("Model frame response variables")

    ## also does not work, not sure what used to happen but now it seems that
    ## the LHS is assumed to be just a single symbol that indexes a DF column

    ## f = x1 + x2 ~ y + x3
    ## mf = ModelFrame(f, d)
    ## @test mf.y_indexes == [1, 2]
    ## @test isequal(mf.formula.lhs, [:(x1 + x2)])
    ## @test isequal(mf.formula.rhs, [:(y + x3)])

    # unique_symbol tests
    #@test unique_symbols(:(x1 + x2)) ==  {"x1"=>:x1, "x2"=>:x2}
    #unique_symbols(:(y ~ x1 + x2 + x3))

    # additional tests from Tom
    y = [1., 2, 3, 4]
    mf = ModelFrame(@formula(y ~ x2), d)
    mm = ModelMatrix(mf)
    @test mm.m == [ones(4) x2]
    @test mm.m == ModelMatrix{sparsetype}(mf).m
    ## @test model_response(mf) == y''     # fails: Int64 vs. Float64

    df = deepcopy(d)
    df[:x1] = PooledDataArray(df[:x1])

    f = @formula(y ~ x2 + x3 + x3*x2)
    mm = ModelMatrix(ModelFrame(f, df))
    @test mm.m == [ones(4) x2 x3 x2.*x3]
    mm = ModelMatrix(ModelFrame(@formula(y ~ x3*x2 + x2 + x3), df))
    @test mm.m == [ones(4) x3 x2 x2.*x3]
    mm = ModelMatrix(ModelFrame(@formula(y ~ x1 + x2 + x3 + x4), df))
    @test mm.m[:,2] == [0, 1., 0, 0]
    @test mm.m[:,3] == [0, 0, 1., 0]
    @test mm.m[:,4] == [0, 0, 0, 1.]
    @test mm.m[:,5] == x2
    @test mm.m[:,6] == x3
    @test mm.m[:,7] == x4

    mm = ModelMatrix(ModelFrame(@formula(y ~ x2 + x3 + x4), df))
    @test mm.m == [ones(4) x2 x3 x4]
    mm = ModelMatrix(ModelFrame(@formula(y ~ x2 + x2), df))
    @test mm.m == [ones(4) x2]
    mm = ModelMatrix(ModelFrame(@formula(y ~ x2*x3 + x2&x3), df))
    @test mm.m == [ones(4) x2 x3 x2.*x3]
    mm = ModelMatrix(ModelFrame(@formula(y ~ x2*x3*x4), df))
    @test mm.m == [ones(4) x2 x3 x4 x2.*x3 x2.*x4 x3.*x4 x2.*x3.*x4]
    mm = ModelMatrix(ModelFrame(@formula(y ~ x2&x3 + x2*x3), df))
    @test mm.m == [ones(4) x2 x3 x2.*x3]

    f = @formula(y ~ x2 & x3 & x4)
    mf = ModelFrame(f, df)
    mm = ModelMatrix(mf)
    @test mm.m == [ones(4) x2.*x3.*x4]
    @test mm.m == ModelMatrix{sparsetype}(mf).m

    f = @formula(y ~ x1 & x2 & x3)
    mf = ModelFrame(f, df)
    mm = ModelMatrix(mf)
    @test mm.m[:, 2:end] == diagm(x2.*x3)
    @test mm.m == ModelMatrix{sparsetype}(mf).m

    #test_group("Column groups in formulas")
    ## set_group was removed in The Great Purge (55e47cd)

    ## set_group(d, "odd_predictors", ["x1","x3"])
    ## @test expand(:odd_predictors, d) == d["odd_predictors"]
    ## mf = ModelFrame(Formula(:(y ~ odd_predictors)), d)
    ## @test mf.df[:,1] == d["y"]
    ## @test mf.df[:,2] == d["x1"]
    ## @test mf.df[:,3] == d["x3"]
    ## @test ncol(mf.df) == 3
    ## mf = ModelFrame(Formula(:(y ~ odd_predictors * x2)), d)
    ## mm = ModelMatrix(mf)
    ## @test mm.model == [ones(4) x1 x3 x2 x1.*x2 x3.*x2]

    ## Interactions between three PDA columns
    ##
    ## FAILS: behavior is wrong when no lower-order terms (1+x1+x2+x1&x2...)
    ##
    ## df = DataFrame(y=1:27,
    ##                x1 = PooledDataArray(vec([x for x in 1:3, y in 4:6, z in 7:9])),
    ##                x2 = PooledDataArray(vec([y for x in 1:3, y in 4:6, z in 7:9])),
    ##                x3 = PooledDataArray(vec([z for x in 1:3, y in 4:6, z in 7:9])))
    ## f = y ~ x1 & x2 & x3
    ## mf = ModelFrame(f, df)
    ## @test coefnames(mf)[2:end] ==
    ##     vec([string("x1: ", x, " & x2: ", y, " & x3: ", z) for
    ##          x in 2:3,
    ##          y in 5:6,
    ##          z in 8:9])

    ## mm = ModelMatrix(mf)
    ## @test mm.m[:,2] == 0. + (df[:x1] .== 2) .* (df[:x2] .== 5) .* (df[:x3].==8)
    ## @test mm.m[:,3] == 0. + (df[:x1] .== 3) .* (df[:x2] .== 5) .* (df[:x3].==8)
    ## @test mm.m[:,4] == 0. + (df[:x1] .== 2) .* (df[:x2] .== 6) .* (df[:x3].==8)
    ## @test mm.m[:,5] == 0. + (df[:x1] .== 3) .* (df[:x2] .== 6) .* (df[:x3].==8)
    ## @test mm.m[:,6] == 0. + (df[:x1] .== 2) .* (df[:x2] .== 5) .* (df[:x3].==9)
    ## @test mm.m[:,7] == 0. + (df[:x1] .== 3) .* (df[:x2] .== 5) .* (df[:x3].==9)
    ## @test mm.m[:,8] == 0. + (df[:x1] .== 2) .* (df[:x2] .== 6) .* (df[:x3].==9)
    ## @test mm.m[:,9] == 0. + (df[:x1] .== 3) .* (df[:x2] .== 6) .* (df[:x3].==9)

    ## Distributive property of :& over :+
    df = deepcopy(d)
    f = @formula(y ~ (x1+x2) & (x3+x4))
    mf = ModelFrame(f, df)
    mm = ModelMatrix(mf)
    @test mm.m == hcat(ones(4), x1.*x3, x1.*x4, x2.*x3, x2.*x4)
    @test mm.m == ModelMatrix{sparsetype}(mf).m

    ## Condensing nested :+ calls
    f = @formula(y ~ x1 + (x2 + (x3 + x4)))
    @test ModelMatrix(ModelFrame(f, df)).m == hcat(ones(4), x1, x2, x3, x4)

    ## Extra levels in categorical column
    mf_full = ModelFrame(@formula(y ~ x1p), d)
    mm_full = ModelMatrix(mf_full)
    @test size(mm_full) == (4,4)

    mf_sub = ModelFrame(@formula(y ~ x1p), d[2:4, :])
    mm_sub = ModelMatrix(mf_sub)
    ## should have only three rows, and only three columns (intercept plus two
    ## levels of factor)
    @test size(mm_sub) == (3,3)

    ## Missing data
    d[:x1m] = @data [5, 6, NA, 7]
    mf = ModelFrame(@formula(y ~ x1m), d)
    mm = ModelMatrix(mf)
    @test mm.m[:, 2] == d[completecases(d), :x1m]
    @test mm.m == ModelMatrix{sparsetype}(mf).m

    ## Same variable on left and right side
    mf = ModelFrame(@formula(x1 ~ x1), df)
    mm = ModelMatrix(mf)
    mm.m == float(model_response(mf))

## Promote non-redundant categorical terms to full rank

d = DataFrame(x = Compat.repeat([:a, :b], outer = 4),
              y = Compat.repeat([:c, :d], inner = 2, outer = 2),
              z = Compat.repeat([:e, :f], inner = 4))
[pool!(d, name) for name in names(d)]
cs = Dict([Pair(name, EffectsCoding()) for name in names(d)])
d[:n] = 1.:8

## No intercept
mf = ModelFrame(@formula(n ~ 0 + x), d, contrasts=cs)
mm = ModelMatrix(mf)
@test mm.m == [1 0
               0 1
               1 0
               0 1
               1 0
               0 1
               1 0
               0 1]
@test mm.m == ModelMatrix{sparsetype}(mf).m
@test coefnames(mf) == ["x: a", "x: b"]

## No first-order term for interaction
mf = ModelFrame(@formula(n ~ 1 + x + x&y), d, contrasts=cs)
mm = ModelMatrix(mf)
@test mm.m[:, 2:end] == [-1 -1  0
                         1  0 -1
                         -1  1  0
                         1  0  1
                         -1 -1  0
                         1  0 -1
                         -1  1  0
                         1  0  1]
@test mm.m == ModelMatrix{sparsetype}(mf).m
@test coefnames(mf) == ["(Intercept)", "x: b", "x: a & y: d", "x: b & y: d"]

## When both terms of interaction are non-redundant:
mf = ModelFrame(@formula(n ~ 0 + x&y), d, contrasts=cs)
mm = ModelMatrix(mf)
@test mm.m == [1 0 0 0
               0 1 0 0
               0 0 1 0
               0 0 0 1
               1 0 0 0
               0 1 0 0
               0 0 1 0
               0 0 0 1]
@test mm.m == ModelMatrix{sparsetype}(mf).m
@test coefnames(mf) == ["x: a & y: c", "x: b & y: c",
                        "x: a & y: d", "x: b & y: d"]

# only a three-way interaction: every term is promoted.
mf = ModelFrame(@formula(n ~ 0 + x&y&z), d, contrasts=cs)
mm = ModelMatrix(mf)
@test mm.m == eye(8)
@test mm.m == ModelMatrix{sparsetype}(mf).m

# two two-way interactions, with no lower-order term. both are promoted in
# first (both x and y), but only the old term (x) in the second (because
# dropping x gives z which isn't found elsewhere, but dropping z gives x
# which is found (implicitly) in the promoted interaction x&y).
mf = ModelFrame(@formula(n ~ 0 + x&y + x&z), d, contrasts=cs)
mm = ModelMatrix(mf)
@test mm.m == [1 0 0 0 -1  0
               0 1 0 0  0 -1
               0 0 1 0 -1  0
               0 0 0 1  0 -1
               1 0 0 0  1  0
               0 1 0 0  0  1
               0 0 1 0  1  0
               0 0 0 1  0  1]
@test mm.m == ModelMatrix{sparsetype}(mf).m
@test coefnames(mf) == ["x: a & y: c", "x: b & y: c",
                        "x: a & y: d", "x: b & y: d",
                        "x: a & z: f", "x: b & z: f"]

# ...and adding a three-way interaction, only the shared term (x) is promoted.
# this is because dropping x gives y&z which isn't present, but dropping y or z
# gives x&z or x&z respectively, which are both present.
mf = ModelFrame(@formula(n ~ 0 + x&y + x&z + x&y&z), d, contrasts=cs)
mm = ModelMatrix(mf)
@test mm.m == [1 0 0 0 -1  0  1  0
               0 1 0 0  0 -1  0  1
               0 0 1 0 -1  0 -1  0
               0 0 0 1  0 -1  0 -1
               1 0 0 0  1  0 -1  0
               0 1 0 0  0  1  0 -1
               0 0 1 0  1  0  1  0
               0 0 0 1  0  1  0  1]
@test mm.m == ModelMatrix{sparsetype}(mf).m
@test coefnames(mf) == ["x: a & y: c", "x: b & y: c",
                        "x: a & y: d", "x: b & y: d",
                        "x: a & z: f", "x: b & z: f",
                        "x: a & y: d & z: f", "x: b & y: d & z: f"]

# two two-way interactions, with common lower-order term. the common term x is
# promoted in both (along with lower-order term), because in every case, when
# x is dropped, the remaining terms (1, y, and z) aren't present elsewhere.
mf = ModelFrame(@formula(n ~ 0 + x + x&y + x&z), d, contrasts=cs)
mm = ModelMatrix(mf)
@test mm.m == [1 0 -1  0 -1  0
               0 1  0 -1  0 -1
               1 0  1  0 -1  0
               0 1  0  1  0 -1
               1 0 -1  0  1  0
               0 1  0 -1  0  1
               1 0  1  0  1  0
               0 1  0  1  0  1]
@test mm.m == ModelMatrix{sparsetype}(mf).m
@test coefnames(mf) == ["x: a", "x: b",
                        "x: a & y: d", "x: b & y: d",
                        "x: a & z: f", "x: b & z: f"]

## FAILS: When both terms are non-redundant and intercept is PRESENT
## (not fully redundant). Ideally, would drop last column. Might make sense
## to warn about this, and suggest recoding x and y into a single variable.
# mf = ModelFrame(n ~ 1 + x&y, d[1:4, :], contrasts=cs)
# @test ModelMatrix(mf).m == [1 1 0 0
#                             1 0 1 0
#                             1 0 0 1
#                             1 0 0 0]
# @test coefnames(mf) == ["x: a & y: c", "x: b & y: c",
#                         "x: a & y: d", "x: b & y: d"]

## note that R also does not detect this automatically. it's left to glm et al.
## to detect numerically when the model matrix is rank deficient, which is hard
## to do correctly.
# > d = data.frame(x = factor(c(1, 2, 1, 2)), y = factor(c(3, 3, 4, 4)))
# > model.matrix(~ 1 + x:y, d)
#   (Intercept) x1:y3 x2:y3 x1:y4 x2:y4
# 1           1     1     0     0     0
# 2           1     0     1     0     0
# 3           1     0     0     1     0
# 4           1     0     0     0     1

# Ensure that random effects terms are dropped from coefnames
df = DataFrame(x = [1,2,3], y = [4,5,6])
mf = ModelFrame(@formula(y ~ 1 + (1 | x)), df)
@test coefnames(mf) == ["(Intercept)"]

mf = ModelFrame(@formula(y ~ 0 + (1 | x)), df)
@test_throws ErrorException ModelMatrix(mf)
@test coefnames(mf) == Vector{Compat.UTF8String}()

# Ensure X is not a view on df column
df = DataFrame(x = [1.0,2.0,3.0], y = [4.0,5.0,6.0])
mf = ModelFrame(@formula(y ~ 0 + x), df)
X = ModelMatrix(mf).m
X[1] = 0.0
@test mf.df[1, :x] == 1.0
end


text:
using StaticArrays

# Karlgaard's formulation of the relative equations of motion
function sphericalkarlgaard(t, x::SVector, u::SVector, v::SVector, μ::Float64)

# Pull out parameters
ρ = x(1)
θ = x(3)
ϕ = x(5)
r = x(7)
f = x(9) # what is f?
ω = x(10)

dρ = x(2)
dθ = x(4)
dϕ = x(6)
dr = x(8)

# Define mu as gravparameter type?
Ga = [1 0 0; 0 1/ρ 0; 0 0 1/ρ]

# Common denominator
Ψ = (r^2 + ρ^2 + 2*r*ρ*cos(θ)*cos(ϕ))^(3/2)

# Spherical equations of motion
fρ = (ω^2 + 2*ω*dθ + dθ^2)*ρ*cos(ϕ)^2 + ρ*dϕ^2 - μ*(r*cos(θ)*cos(ϕ) + ρ)/Ψ + μ*cos(θ)*cos(ϕ)/r^2
fθ = 2*(ω + dθ)*dϕ*tan(ϕ) - dω - 2*(ω + dθ)*dρ/ρ + μ*(r*sin(θ)*sec(ϕ))/ρ*Ψ - μ*sin(θ)*sec(ϕ)/ρ*r^2
fϕ = -0.5*sin(2*ϕ)*(ω + dθ)^2 - 2*dρ*dϕ/ρ - μ*cos(θ)*sin(ϕ)/ρ*r^2 + μ*r*cos(θ)*sin(ϕ)/ρ*Ψ
fa = [fρ;fθ;fϕ]

# Need control noise term here

# Process noise term - Why did Jimmy multiply by gaussian random variables?
v = diagm(randn(3))*sqrt.(v)

# Full relative motion acceleration with controls, η = {ρ;θ;ϕ}
ddη = fa + Ga*(u*v)
ddρ = ddη(1)
ddθ = ddη(2)
ddϕ = ddη(3)
ddr = r*ω^2 - μ/r^2
dω = -2*dr*ω/r

# State deriviative (velocity and acceleration)
dx = [dρ;ddρ;dθ;ddθ;dϕ;ddϕ;dr;ddr;ω;dω]
return dx

end

max_stars_repo_name: glemieux/RelativeEstimationProject

max_issues_repo_name: glemieux/RelativeEstimationProject

max_forks_repo_name: glemieux/RelativeEstimationProject


text:
module DICOMClient

import Base64
import DICOM
import DocStringExtensions
import FHIRClient
import FileIO
import HTTP
import ImageMagick
import Images
import JSON3
import Pkg
import StructTypes

include("types.jl")
include("version.jl")

include("authentication.jl")
include("authentication-headers.jl")
include("base-url.jl")
include("client.jl")
include("headers.jl")
include("image-utils.jl")
include("multipart-images.jl")
include("requests.jl")
include("uri.jl")
include("version-utils.jl")

end # end module DICOMClient


text:
export SrcDistRowMatrix, DistRowMatrix, RowMatrix
export isFillActive, isLocallyIndexed
export getGraph, getGlobalRowCopy, getLocalRowCopy, getGlobalRowView, getLocalRowView, getLocalDiagCopy, leftScale!, rightScale!
export localApply, mul!

"""
RowMatrix is the base type for all row oriented Petra matrices.
RowMatrix fufils both the Operator and DistObject interfaces.

    getGraph(mat::RowMatrix)::RowGraph
Returns the graph that represents the structure of the row matrix

    getLocalRowCopy!(copy::Tuple{<:AbstractVector{<:Integer}, <:AbstractVector{Data}}, matrix::RowMatrix{Data, GID, PID, LID}, localRow::LID)::Integer
Copies the given row into the provided arrays and returns the number of elements in that row using local indices

    getGlobalRowCopy!(copy::Tuple{<:AbstractVector{<:Integer}, <:AbstractVector{Data}}, matrix::RowMatrix{Data, GID, PID, LID}, globalRow::GID)::Integer
Copies the given row into the provided arrays and returns the number of elements in that row using global indices

    getGlobalRowView(matrix::RowMatrix{Data, GID, PID, LID}, globalRow::Integer)::Tuple{AbstractArray{GID, 1}, AbstractArray{Data, 1}}
Returns a view to the given row using global indices

    getLocalRowView(matrix::RowMatrix{Data, GID, PID, LID},localRow::Integer)::Tuple{AbstractArray{GID, 1}, AbstractArray{Data, 1}}
Returns a view to the given row using local indices

    getLocalDiagCopy!(copy::MultiVector{Data, GID, PID, LID}, matrix::RowMatrix{Data, GID, PID, LID})::MultiVector{Data, GID, PID, LID}
Copies the local diagonal into the given `MultiVector` then returns the `MultiVector`

    leftScale!(matrix::RowMatrix{Data, GID, PID, LID}, X::AbstractArray{Data, 1})
Scales matrix on the left with X

    rightScale!(matrix::RowMatrix{Data, GID, PID, LID}, X::AbstractArray{Data, 1})
Scales matrix on the right with X

`getMap(...)`, as required by SrcDistObject, is implemented by calling `getRowMap(...)`

`apply!(...)`, as required by Operator, is implemented, but can be optimized by overrideing the following method
    localApply(Y::MultiVector, A::RowMatrix, X::MultiVector, ::TransposeMode, α::Data, β::Data)
Does the computations for `Y = β⋅Y + α⋅A⋅X`, `X` and `Y` match the row map and column map, depending on the transpose mode

The following methods are currently implemented as no-ops, but can be overridden to improve performance.

    setColumnMapMultiVector(::RowMatrix{Data, GID, PID, LID}, ::Union{MultiVector{Data, GID, PID, LID}, Nothing})
Caches a `MultiVector` that uses the matrix's column map.

    getColumnMapMultiVector(::RowMatrix{Data, GID, PID, LID})::Union{MultiVector{Data, GID, PID, LID}, Nothing}
Fetches any cached `MultiVector` that uses the matrix's column map.

    setRowMapMultiVector(::RowMatrix{Data, GID, PID, LID}, ::Union{MultiVector{Data, GID, PID, LID}, Nothing})
Caches a `MultiVector` that uses the matrix's row map.

    getRowMapMultiVector(::RowMatrix{Data, GID, PID, LID})::Union{MultiVector{Data, GID, PID, LID}, Nothing}
Fetches any cached `MultiVector` that uses the matrix's row map.

Some pre-implemented methods can be optimized by providing specialized implementations
`apply!`, as mentioned above
All `RowMatrix` methods that are also implemented by `RowGraph` are implemented using `getGraph`.
`pack` is implemented using `getLocalRowCopy`
`getGlobalRowCopy!` is implemented by calling `getLocalRowCopy!` and remapping the values using `gid(::BlockMap, ::Integer)`

Additionally, Julia's `mul!` and `*` functions are implemented for `RowMatrix`-`MultiVector` products
"""
abstract type RowMatrix{Data <: Number, GID <: Integer, PID <: Integer, LID <: Integer} <: AbstractArray{Data, 2}
end

isFillActive(matrix::RowMatrix) = !isFillComplete(matrix)
isLocallyIndexed(matrix::RowMatrix) = !isGloballyIndexed(matrix)

function leftScale!(matrix::RowMatrix{Data, GID, PID, LID}, X::MultiVector{Data, GID, PID, LID}) where {
        Data <: Number, GID <: Integer, PID <: Integer, LID <: Integer}
    if numVectors(X) != 1
        throw(InvalidArgumentError("Can only scale row matrix with column vector, not multi vector"))
    end
    leftScale!(matrix, getLocalArray(X))
end

function rightScale!(matrix::RowMatrix{Data, GID, PID, LID}, X::MultiVector{Data, GID, PID, LID}) where {
        Data <: Number, GID <: Integer, PID <: Integer, LID <: Integer}
    if numVectors(X) != 1
        throw(InvalidArgumentError("Can only scale row matrix with column vector, not multi vector"))
    end
    rightScale!(matrix, getLocalArray(X))
end

#for SrcDistObject
getMap(matrix::RowMatrix) = getRowMap(matrix)

#TODO document
function getLocalDiagCopyWithoutOffsetsNotFillComplete(A::RowMatrix{Data, GID, PID, LID})::MultiVector{Data, GID, PID, LID} where {Data, GID, PID, LID}

    localRowMap = getLocalMap(getRowMap(A))
    localColMap = getLocalMap(getColMap(A))
    sorted = isSorted(A.myGraph)

    localNumRows = getLocalNumRows(A)
    diag = MultiVector{Data, GID, PID, LID}(getRowMap(A), 1)
    diagLocal1D = getVectorView(diag, 1)

    range = LID(1):localNumRows
    for localRowIndex in range
        diagLocal1D[localRowIndex] = 0
        globalIndex = gid(localRowMap, localRowIndex)
        localColIndex = lid(localColMap, globalIndex)
        if localColIndex != 0
            indices, values = getLocalRowView(A, localRowIndex)

            if !sorted
                offset = findfirst(indices, localColumnIndex)
            else
                offset = searchsorted(indices, localColumnIndex)
            end

            if offset <= length(indices)
                diagLocal1D[localRowIndex] = values[offset]
            end
        end
    end
    diag
end

"""
    getGlobalRowCopy(matrix::RowMatrix{Data, GID, PID, LID}, globalRow::Integer)::Tuple{AbstractArray{GID, 1}, AbstractArray{Data, 1}}

Returns a copy of the given row using global indices
"""
Base.@propagate_inbounds function getGlobalRowCopy(matrix::RowMatrix{Data, GID, PID, LID}, globalRow::Integer) where {Data, GID, PID, LID}
    numEntries = getNumEntriesInGlobalRow(matrix, globalRow)
    copy = (Vector{GID}(undef, numEntries), Vector{Data}(undef, numEntries))
    getGlobalRowCopy!(copy, matrix, globalRow)
    copy
end

"""
    getLocalRowCopy(matrix::RowMatrix{Data, GID, PID, LID}, localRow::Integer)::Tuple{AbstractArray{LID, 1}, AbstractArray{Data, 1}}

Returns a copy of the given row using local indices
"""
Base.@propagate_inbounds function getLocalRowCopy(matrix::RowMatrix{Data, GID, PID, LID}, localRow::Integer
        )::Tuple{AbstractArray{LID, 1}, AbstractArray{Data, 1}} where {Data, GID, PID, LID}
    numEntries = getNumEntriesInGlobalRow(matrix, localRow)
    copy = (Vector{LID}(undef, numEntries), Vector{Data}(undef, numEntries))
    getLocalRowCopy!(copy, matrix, localRow)
    copy
end

"""
    getGlobalRowCopy!(copy::Tuple{<:AbstractVector{<:Integer}, <:AbstractVector{Data}}, matrix::RowMatrix{Data, GID, PID, LID}, globalRow::Integer)::Integer

Copies the given row into the provided arrays and returns the number of elements in that row using global indices
"""
Base.@propagate_inbounds function getGlobalRowCopy!(copy::Tuple{<:AbstractVector{<:Integer}, <:AbstractVector{Data}},
            matrix::RowMatrix{Data, GID}, globalRow::Integer) where {Data, GID}
    getGlobalRowCopy!(copy, matrix, GID(globalRow))
end

Base.@propagate_inbounds function getGlobalRowCopy!(copy::Tuple{<:AbstractVector{<:Integer}, <:AbstractVector{Data}},
            matrix::RowMatrix{Data, GID}, globalRow::GID) where {Data, GID}
    rowMap = getRowMap(matrix)
    numElts = getLocalRowCopy!(copy, matrix, lid(rowMap, globalRow))
    inds = copy[1]
    @. inds = gid(rowMap, inds)
    numElts
end

"""
    getLocalRowCopy!(copy::Tuple{<:AbstractVector{LID}, <:AbstractVector{Data}}, matrix::RowMatrix{Data, GID, PID, LID}, localRow::Integer)::Integer

Copies the given row into the provided arrays and returns the number of elements in that row using local indices
"""
Base.@propagate_inbounds function getLocalRowCopy!(copy::Tuple{<:AbstractVector{LID}, <:AbstractVector{Data}},
            matrix::RowMatrix{Data, GID, PID, LID}, localRow::Integer) where {Data, GID, PID, LID}
    getLocalRowCopy!(copy, matrix, LID(localRow))
end

"""
    getLocalDiagCopy(matrix::RowMatrix{Data, GID, PID, LID})::MultiVector{Data, GID, PID, LID}

Returns a copy of the diagonal elements on the calling processor
"""
Base.@propagate_inbounds function getLocalDiagCopy(matrix::RowMatrix{Data, GID, PID, LID}) where {Data, GID, PID, LID}
    copy = DenseMultiVector{Data}(getRowMap(matrix), 1, false)
    getLocalDiagCopy!(copy, matrix)
    copy
end

function packAndPrepare(source::RowMatrix{Data, GID, PID, LID},
        target::RowMatrix{Data, GID, PID, LID}, exportLIDs::AbstractArray{LID, 1},
        distor::Distributor{GID, PID, LID})::AbstractArray where {Data, GID, PID, LID}
    pack(source, exportLIDs, distor)
end

"""
    pack(::RowMatrix{Data, GID, PID, LID}, exportLIDs::AbstractVector{LID}, distor::Distributor{GID, PID, LID})::AbstractArray{Tuple{AbstractVector{GID}, AbstractVector{Data}}}
Packs this object's data for import or export
"""
function pack(source::RowMatrix{Data, GID, PID, LID}, exportLIDs::AbstractVector{LID},
        distor::Distributor{GID, PID, LID})::AbstractArray{Tuple{AbstractVector{GID}, AbstractVector{Data}}} where{Data, GID, PID, LID}
    srcMap = getMap(source)
    map(lid->getGlobalRowCopy(source, gid(srcMap, lid)), exportLIDs)
end

"""
    createColumnMapMultiVector(mat::RowMatrix, x::MultiVector; force=false)::Union{MultiVector, Nothing}

Returns a `MultiVector` that uses the matrix's column map or `nothing` if `x` can be used
If `getImporter(mat)` is null (ie a trivial import), then the multivector will only be created if `force` is true.
"""
function createColumnMapMultiVector(mat::RowMatrix{Data, GID, PID, LID}, X::MultiVector{Data, GID, PID, LID}; force = false) where {Data, GID, PID, LID}
    if !isFillComplete(mat)
        throw(InvalidStateError("Can only call createColumnMapMultiVector if the matrix is fill active"))
    end
    if !hasColMap(mat)
        throw(InvalidStateError("Can only call createColumnMapMultiVector with a matrix that has a column map"))
    end

    numVecs = numVectors(X)
    importer = getImporter(mat)
    colMap = getColMap(mat)

    #if import object is trivial, don't need a seperate column map multivector
    if importer != nothing || force
        importMV = getColumnMapMultiVector(mat)
        if importMV == nothing || numVectors(importMV) != numVecs
            importMV = DenseMultiVector{Data}(colMap, numVecs)
            setColumnMapMultiVector(mat, importMV)
        end
        importMV
    else
        nothing
    end
end

"""
    createRowMapMultiVector(mat::RowMatrix, x::MultiVector; force=false)::Union{MultiVector, Nothing}

Returns a `MultiVector` that uses the matrix's row map or `nothing` if `x` can be used
If `getExporter(mat)` is null (ie a trivial export), then the multivector will only be created if `force` is true.
"""
function createRowMapMultiVector(mat::RowMatrix{Data, GID, PID, LID}, Y::MultiVector{Data, GID, PID, LID}; force = false) where {Data, GID, PID, LID}
    if !isFillComplete(mat)
        throw(InvalidStateError("Cannot call createRowMapMultiVector if the matrix is fill active"))
    end

    numVecs = numVectors(Y)
    exporter = getExporter(operator)
    rowMap = getRowMap(mat)

    if exporter != nothing || force
        exportMV = getRowMapMultiVector(mat)
        if exportMV == nothing || getNumVectors(exportMV) != numVecs
            exportMV = DenseMultiVector{Data}(rowMap, numVecs)
        end
        exportMV
    else
        nothing
    end
end

function apply!(Y::MultiVector{Data, GID, PID, LID},
        operator::RowMatrix{Data, GID, PID, LID}, X::MultiVector{Data, GID, PID, LID},
        mode::TransposeMode, alpha::Data, beta::Data) where {Data, GID, PID, LID}

    ZERO = Data(0)

    if isFillActive(operator)
        throw(InvalidStateError("Cannot call apply(...) until fillComplete(...)"))
    end

	if alpha == ZERO
        if beta == ZERO
            fill!(Y, ZERO)
        elseif beta != Data(1)
            scale!(Y, beta)
        end
        return Y
    end

    importer = getImporter(operator)
    exporter = getExporter(operator)

    YIsReplicated = !distributedGlobal(Y)
    YIsOverwritted = (beta == ZERO)
    if YIsReplicated && myPid(getComm(operator)) != 1
        beta = ZERO
    end

    if mode == NO_TRANS
        if importer == nothing
            XColMap = X
        else
            #need to import source multivector
            XColMap = createColumnMapMultiVector(operator, X)
            doImport(X, XColMap, importer, INSERT)
        end

        if exporter != nothing
            YRowMap = createRowMapMultiVector(operator, Y)
            localApply(YRowMap, operator, XColMap, NO_TRANS, alpha, ZERO)

            if YIsOverwritten
                fill!(Y, ZERO)
            else
                scale!(Y, beta)
            end

            doExport(YRowMap, Y, exporter, ADD)
        else
            #don't do export row Map and range map are the same
            if XColMap === Y
                YRowMap = createRowMapMultiVector(operator, Y; force=true)

                if beta != 0
                    copy!(YRowMap, Y)
                end

                localApply(YRowMap, operator, XColmap, NO_TRANS, alpha, ZERO)
                copy!(Y, YRowMap)
            else
                localApply(Y, operator, XColMap, NO_TRANS, alpha, beta)
            end
        end
    else
        if exporter == nothing
            XRowMap = X
        else
            rowMapMV = createRowMapMultiVector(mat, X)
            doImport(X, rowMapMV, exporter, INSERT)
            XRowMap = rowMapMV
        end

        if importer != nothing
            YColMap = createColumnMapMultiVector(mat, X)
            localApply(YColMap, operator, XRowMap, mode, alpha, ZERO)

            if YIsOverwritten
                fill!(Y, ZERO)
            else
                scale!(Y, beta)
            end
            doExport(YColMap, Y, importer, ADD)
        else
            if XRowMap === Y
                YCopy = copy(Y)
                localApply(YCopy, operator, XRowMap, mode, alpha, beta)
                copy!(Y, YCopy)
            else
                localApply(Y, operator, XRowMap, mode, alpha, beta)
            end
        end
    end
    if YIsReplicated
        commReduce(Y)
    end
    Y
end

function localApply(Y::MultiVector{Data, GID, PID, LID},
        A::RowMatrix{Data, GID, PID, LID}, X::MultiVector{Data, GID, PID, LID},
        mode::TransposeMode, alpha::Data, beta::Data) where {Data, GID, PID, LID}

    rawY = getLocalArray(Y)
    rawX = getLocalArray(Y)

    maxElts = getLocalMaxNumRowEntries(A)
    indices = Vector{LID}(undef, maxElts)
    values = Vector{LID}(undef, maxElts)

    if !isTransposed(mode)
        numRows = getLocalNumRows(A)
        for vect = LID(1):numVectors(Y)
            @inbounds for row = LID(1):numRows
                sum::Data = Data(0)
                numElts = getLocalRowCopy!((indices, values), A, row)
                for i in LID(1):numElts
                    ind::LID = indices[i]
                    val::Data = values[i]
                    sum += val*rawX[ind, vect]
                end
                sum = applyConjugation(mode, sum*alpha)
                @inbounds rawY[row, vect] *= beta
                @inbounds rawY[row, vect] += sum
            end
        end
    else
        rawY[:, :] *= beta
        numRows = getLocalNumRows(A)
        for vect = LID(1):numVectors(Y)
            @inbounds for mRow in LID(1):numRows
                numElts = getLocalRowCopy!((indices, values), A, row)
                for i in LID(1):numElts
                    ind::LID = indices[i]
                    val::Data = values[i]
                    rawY[ind, vect] += applyConjugation(mode, alpha*rawX[mRow, vect]*val)
                end
            end
        end
    end

    Y
end

#### Julia LA API wrappers ####

function LinearAlgebra.mul!(Y::MultiVector{Data, GID, PID, LID},
        	operator::RowMatrix{Data, GID, PID, LID},
			X::MultiVector{Data, GID, PID, LID}) where {Data, GID, PID, LID}
	apply!(Y, operator, X)
end

function Base.:*(operator::RowMatrix{Data, GID, PID, LID},
			X::MultiVector{Data, GID, PID, LID}) where {Data, GID, PID, LID}
	Y = DenseMultiVector{Data}(getRowMap(operator), numVectors(X), true)
	mul!(Y, operator, X)
end

#### default implementations ####

"""
    setColumnMapMultiVector(::RowMatrix{Data, GID, PID, LID}, ::Union{MultiVector{Data, GID, PID, LID}, Nothing})

Caches a `MultiVector` that uses the matrix's column map.
"""
function setColumnMapMultiVector(::RowMatrix{Data, GID, PID, LID}, ::Union{MultiVector{Data, GID, PID, LID}, Nothing}) where{Data, GID, PID, LID}
    nothing
end

"""
    getColumnMapMultiVector(::RowMatrix{Data, GID, PID, LID})::Union{MultiVector{Data, GID, PID, LID}, Nothing}

Fetches any cached `MultiVector` that uses the matrix's column map.
"""
function getColumnMapMultiVector(::RowMatrix{Data, GID, PID, LID}) where{Data, GID, PID, LID}
    nothing
end

"""
    setRowMapMultiVector(::RowMatrix{Data, GID, PID, LID}, ::Union{MultiVector{Data, GID, PID, LID}, Nothing})

Caches a `MultiVector` that uses the matrix's row map.
"""
function setRowMapMultiVector(::RowMatrix{Data, GID, PID, LID}, ::Union{MultiVector{Data, GID, PID, LID}, Nothing}) where{Data, GID, PID, LID}
    nothing
end

"""
    getRowMapMultiVector(::RowMatrix{Data, GID, PID, LID})::Union{MultiVector{Data, GID, PID, LID}, Nothing}

Fetches any cached `MultiVector` that uses the matrix's row map.
"""
function getRowMapMultiVector(::RowMatrix{Data, GID, PID, LID}) where{Data, GID, PID, LID}
    nothing
end

#### default implementations using getGraph(...) ####
"""
    isFillComplete(mat::RowMatrix)

Whether `fillComplete(...)` has been called
"""
isFillComplete(mat::RowMatrix) = isFillComplete(getGraph(mat))

"""
    getRowMap(::RowMatrix{Data, GID, PID, LID})::BlockMap{GID, PID, LID}

Gets the row map for the container
"""
getRowMap(mat::RowMatrix) = getRowMap(getGraph(mat))

"""
    getColMap(::RowMatrix{Data, GID, PID, LID})::BlockMap{GID, PID, LID}

Gets the column map for the container
"""
getColMap(mat::RowMatrix) = getColMap(getGraph(mat))

"""
    hasColMap(::RowMatrix)::Bool

Whether the container has a well-defined column map
"""
hasColMap(mat::RowMatrix) = hasColMap(getGraph(mat))

"""
    getImporter(::RowMatrix{Data, GID, PID, LID})::Union{Import{GID, PID, LID}, Nothing}

Gets the `Import` object for the matrix
"""
getImporter(mat::RowMatrix) = getImporter(getGraph(mat))

"""
    getExporter(::RowMatrix{Data, GID, PID, LID})::Union{Export{GID, PID, LID}, Nothing}

Gets the `Export` object for the matrix
"""
getExporter(mat::RowMatrix) = getExporter(getGraph(mat))

"""
    isGloballyIndexed(mat::RowMatrix)

Whether the matrix stores indices with global indexes
"""
isGloballyIndexed(mat::RowMatrix) = isGloballyIndexed(getGraph(mat))

"""
    getGlobalNumRows(mat::RowMatrix)

Returns the number of rows across all processors
"""
getGlobalNumRows(mat::RowMatrix) = getGlobalNumRows(getGraph(mat))

"""
    getGlobalNumCols(mat::RowMatrix)

Returns the number of columns across all processors
"""
getGlobalNumCols(mat::RowMatrix) = getGlobalNumCols(getGraph(mat))

"""
    getLocalNumRows(mat::RowMatrix)

Returns the number of rows on the calling processor
"""
getLocalNumRows(mat::RowMatrix) = getLocalNumRows(getGraph(mat))

"""
    getLocalNumCols(mat::RowMatrix)

Returns the number of columns on the calling processor
"""
getLocalNumCols(mat::RowMatrix) = getLocalNumCols(getGraph(mat))

"""
    getGlobalNumEntries(mat::RowMatrix)

Returns the number of entries across all processors
"""
getGlobalNumEntries(mat::RowMatrix) =  getGlobalNumEntries(getGraph(mat))

"""
    getLocalNumEntries(mat::RowMatrix)

Returns the number of entries on the calling processor
"""
getLocalNumEntries(mat::RowMatrix) = getLocalNumEntries(getGraph(mat))

"""
    getNumEntriesInGlobalRow(mat::RowMatrix, globalRow)

Returns the number of entries on the local processor in the given row
"""
Base.@propagate_inbounds getNumEntriesInGlobalRow(mat::RowMatrix, globalRow) = getNumEntriesInGlobalRow(getGraph(mat), globalRow)

"""
    getNumEntriesInLocalRow(mat::RowMatrix, localRow)

Returns the number of entries on the local processor in the given row
"""
Base.@propagate_inbounds getNumEntriesInLocalRow(mat::RowMatrix, localRow) = getNumEntriesInLocalRow(getGraph(mat), localRow)

"""
    getGlobalNumDiags(mat::RowMatrix)

Returns the number of diagonal elements across all processors
"""
getGlobalNumDiags(mat::RowMatrix) = getGlobalNumDiags(getGraph(mat))

"""
    getLocalNumDiags(mat::RowMatrix)

Returns the number of diagonal element on the calling processor
"""
getLocalNumDiags(mat::RowMatrix) = getLocalNumDiags(getGraph(mat))

"""
    getGlobalMaxNumRowEntries(mat::RowMatrix)

Returns the maximum number of row entries across all processors
"""
getGlobalMaxNumRowEntries(mat::RowMatrix) = getGlobalMaxNumRowEntries(getGraph(mat))

"""
    getLocalMaxNumRowEntries(mat::RowMatrix)

Returns the maximum number of row entries on the calling processor
"""
getLocalMaxNumRowEntries(mat::RowMatrix) = getLocalMaxNumRowEntries(getGraph(mat))

"""
    isLowerTriangular(mat::RowMatrix)

Whether the matrix is lower triangular
"""
isLowerTriangular(mat::RowMatrix) = isLowerTriangular(getGraph(mat))

"""
    isUpperTriangular(mat::RowMatrix)

Whether the matrix is upper triangular
"""
isUpperTriangular(mat::RowMatrix) = isUpperTriangular(getGraph(mat))

getDomainMap(mat::RowMatrix) = getDomainMap(getGraph(mat))
getRangeMap(mat::RowMatrix) = getRangeMap(getGraph(mat))

#### required method documentation stubs ####

"""
    getGraph(mat::RowMatrix)

Returns the graph that represents the structure of the row matrix
"""
function getGraph end

"""
    getGlobalRowView(matrix::RowMatrix{Data, GID, PID, LID},globalRow::Integer)::Tuple{AbstractArray{GID, 1}, AbstractArray{Data, 1}}

Returns a view to the given row using global indices
"""
function getGlobalRowView end

"""
    getLocalRowView(matrix::RowMatrix{Data, GID, PID, LID},localRow::Integer)::Tuple{AbstractArray{GID, 1}, AbstractArray{Data, 1}}

Returns a view to the given row using local indices
"""
function getLocalRowView end

"""
    getLocalDiagCopy!(copy::MultiVector{Data, GID, PID, LID}, matrix::RowMatrix{Data, GID, PID, LID})::MultiVector{Data, GID, PID, LID}

Copies the local diagonal into the given `MultiVector` then returns the `MultiVector`
"""
function getLocalDiagCopy! end

"""
    leftScale!(matrix::RowMatrix{Data, GID, PID, LID}, X::AbstractArray{Data})

Scales matrix on the left with X
"""
function leftScale! end

"""
    rightScale!(matrix::RowMatrix{Data, GID, PID, LID}, X::AbstractArray{Data})

Scales matrix on the right with X
"""
function rightScale! end

### Julia Array functions ###
Base.size(mat::RowMatrix) = (getGlobalNumRows(mat), getGlobalNumCols(mat))

#TODO this might break for funky maps, however indices needs to return a unit range
Base.axes(A::RowMatrix{GID}) where GID = if hasColMap(A)
        (minMyGID(getRowMap(A)):maxMyGID(getRowMap(A)), minMyGID(getColMap(A)):maxMyGID(getColMap(A)))
    else
        (minMyGID(getRowMap(A)):maxMyGID(getRowMap(A)), GID(1):getGlobalNumCols(A))
    end

function Base.getindex(A::RowMatrix, I::Vararg{Int, 2})
    if isGloballyIndexed(A)
        @boundscheck begin
            (n, m) = size(A)
            if I[1] > n || I[1] < 1 || I[2] > m || I[2] < 1
                throw(BoundsError(A, I))
            end
        end
        (rowInds, rowVals) = getGlobalRowView(A, I[0])
        for i in 1:length(rowInds)
            if rowInds[i] == I[1]
                return rowVals[i]
            end
        end
    else
        lRow = lid(getMap(A), I[1])
        lCol = lid(getMap(A), I[2])
        (rowInds, rowVals) = getLocalRowView(A, lRow)
        for i in 1:length(rowInds)
            if rowInds[i] == lCol
                return rowVals[i]
            end
        end
    end
    return 0
end

#TODO look into setindex!

function Base.show(io::IO, mat::RowMatrix)
    print(io, "$(typeof(mat)) on process $(myPid(getComm(mat)))")
end