text:
utc = FixedTimeZone("UTC")
warsaw = resolve("Europe/Warsaw", tzdata["europe"]...)

# Converting a ZonedDateTime into a DateTime
dt = DateTime(2015, 1, 1, 0)
zdt = ZonedDateTime(dt, warsaw)
@test DateTime(zdt) == dt

# Converting from ZonedDateTime to DateTime isn't possible as it is always inexact.
@test_throws MethodError convert(DateTime, zdt)

# Vectorized accessors
arr = repmat([zdt], 10)
@test Dates.DateTime(arr) == repmat([dt], 10)

# now function
dt = Dates.unix2datetime(time())  # Base.now in UTC
zdt = now(warsaw)
@test zdt.timezone == warsaw
@test isapprox(map(Dates.datetime2unix, [dt, TimeZones.utc(zdt)])...)

# Changing time zones
dt = DateTime(2015, 1, 1, 0)
zdt_utc = ZonedDateTime(dt, utc; from_utc=true)
zdt_warsaw = ZonedDateTime(dt, warsaw; from_utc=true)

# Identical since ZonedDateTime is immutable
@test astimezone(zdt_utc, warsaw) === zdt_warsaw
@test astimezone(zdt_warsaw, utc) === zdt_utc


text:
using Test
using DandelionWebSockets
using DandelionWebSockets.Proxy: WebSocketsHandlerProxy, stopproxy
import DandelionWebSockets: on_text, on_binary
import DandelionWebSockets: state_connecting, state_open, state_closing, state_closed

struct WebSocketsHandlerNotification
    method_name::String
    payload::Any
end

struct FakeWebSocketConnection <: WebSocketConnection end

"""
This is a mock WebSocketHandler which is designed to verify that calls are made to the on_*
and state_* methods in a task separate from the client logic (which is the test in this
scenario). To verify this, the WebSocketHandler uses an unbuffered channel to inform the test
of callbacks. Therefore, if the test runs the `on_text()` method in the same task, then it will
block, because it cannot simultaneously take that notification from the channel.
"""
struct MockWebSocketsHandler <: WebSocketHandler
    notification::Channel{WebSocketsHandlerNotification}

    MockWebSocketsHandler() = new(Channel{WebSocketsHandlerNotification}(0))
end

takenotification!(h::MockWebSocketsHandler) = take!(h.notification)

on_text(h::MockWebSocketsHandler, text::String) = put!(h.notification, WebSocketsHandlerNotification("on_text", text))
on_binary(h::MockWebSocketsHandler, data::AbstractVector{UInt8}) = put!(h.notification, WebSocketsHandlerNotification("on_binary", data))
state_connecting(h::MockWebSocketsHandler, conn::WebSocketConnection) = put!(h.notification, WebSocketsHandlerNotification("state_connecting", ""))
state_open(h::MockWebSocketsHandler) = put!(h.notification, WebSocketsHandlerNotification("state_open", ""))
state_closing(h::MockWebSocketsHandler) = put!(h.notification, WebSocketsHandlerNotification("state_closing", ""))
state_closed(h::MockWebSocketsHandler) = put!(h.notification, WebSocketsHandlerNotification("state_closed", ""))

@testset "User handler proxy     " begin
    @testset "Callbacks are done in a separate task" begin
        @testset "on_text callback" begin
            # Arrange
            handler = MockWebSocketsHandler()
            handlerproxy = WebSocketsHandlerProxy(handler)

            # Act
            on_text(handlerproxy, "Some text")

            # Assert
            notification = takenotification!(handler)
            @test notification.method_name == "on_text"
            @test notification.payload == "Some text"
        end

        @testset "on_binary callback" begin
            # Arrange
            handler = MockWebSocketsHandler()
            handlerproxy = WebSocketsHandlerProxy(handler)

            # Act
            on_binary(handlerproxy, b"Some binary")

            # Assert
            notification = takenotification!(handler)
            @test notification.method_name == "on_binary"
            @test notification.payload == b"Some binary"
        end

        @testset "state_connecting callback" begin
            # Arrange
            handler = MockWebSocketsHandler()
            handlerproxy = WebSocketsHandlerProxy(handler)

            # Act
            state_connecting(handlerproxy, FakeWebSocketConnection())

            # Assert
            notification = takenotification!(handler)
            @test notification.method_name == "state_connecting"
        end

        @testset "state_open callback" begin
            # Arrange
            handler = MockWebSocketsHandler()
            handlerproxy = WebSocketsHandlerProxy(handler)

            # Act
            state_open(handlerproxy)

            # Assert
            notification = takenotification!(handler)
            @test notification.method_name == "state_open"
        end

        @testset "state_closing callback" begin
            # Arrange
            handler = MockWebSocketsHandler()
            handlerproxy = WebSocketsHandlerProxy(handler)

            # Act
            state_closing(handlerproxy)

            # Assert
            notification = takenotification!(handler)
            @test notification.method_name == "state_closing"
        end

        @testset "state_closed callback" begin
            # Arrange
            handler = MockWebSocketsHandler()
            handlerproxy = WebSocketsHandlerProxy(handler)

            # Act
            state_closed(handlerproxy)

            # Assert
            notification = takenotification!(handler)
            @test notification.method_name == "state_closed"
        end

        @testset "stop the proxy" begin
            # Arrange
            handler = MockWebSocketsHandler()
            handlerproxy = WebSocketsHandlerProxy(handler)

            # Act
            stopproxy(handlerproxy)

            # Assert
            @test isopen(handlerproxy.callbacks) == false
        end
    end
end

ext: jl


text:
import MPI

# min x1^2 + x2^2 + x3^2 + x4^2 + x5^2 + x6^2
# st.
#     x1 + x2 = 100
#     0< x2 + x3 + x4  < 500
#     0< x1 + x5 + x6  < 500
# x free variables

include("../../JuliaInterface/ParPipsNlp.jl")

using ParPipsNlp

function str_init_x0(nodeid, x0)
    assert(length(x0) ==2)
    if(nodeid == 0)
        x0[1] = 1.0
        x0[2] = 1.0
    elseif nodeid == 1
        x0[1] = 1.0
        x0[2] = 1.0
    elseif nodeid == 2
        x0[1] = 1.0
        x0[2] = 1.0
    end
end

function str_prob_info(nodeid,mode,col_lb,col_ub,row_lb,row_ub)

    if(mode == :Values)
        #setting col lb ub, row lb ub
        if(nodeid==0)
            assert(length(row_lb) == length(row_ub))
            assert(length(col_lb) == length(col_ub))
            fill!(col_lb,-Inf)
            fill!(col_ub, Inf)
            fill!(row_lb, 100.0)
            fill!(row_ub, 100.0)
        elseif (nodeid == 1 )
            assert(length(row_lb) == length(row_ub))
            assert(length(col_lb) == length(col_ub))
            fill!(col_lb, -Inf)
            fill!(col_ub, Inf)
            fill!(row_lb, 0)
            fill!(row_ub, 500)
        elseif (nodeid == 2)
            assert(length(row_lb) == length(row_ub))
            assert(length(col_lb) == length(col_ub))
            fill!(col_lb, -Inf)
            fill!(col_ub, Inf)
            fill!(row_lb, 0)
            fill!(row_ub, 500)
        else
            assert(false)
        end

        # @show col_lb
        # @show col_ub
        # @show row_lb
        # @show row_ub
    end

    if(nodeid==0)
        return (2,1)
    elseif(nodeid == 1 || nodeid == 2)
        return (2,1)
    else
        assert(false)
    end
end

function str_eval_f(nodeid,x0,x1)
    @show nodeid
    # @show x0
    # @show x1
    fval = 0.0
    if nodeid ==  0
        x11 = x0[1]
        x22 = x0[2]
        fval = x11*x11 + x22*x22
    elseif nodeid == 1
        x33 = x1[1]
        x44 = x1[2]
        fval = x33*x33 + x44*x44
    elseif nodeid == 2
        x55 = x1[1]
        x66 = x1[1]
        fval = x55*x55 + x66*x66
    else
        assert(false)
    end
    return fval
end

function str_eval_g(nodeid,x0,x1,new_eq_g, new_inq_g)
    # @show nodeid
    x11 = x0[1]
    x22 = x0[2]
    if nodeid == 0
        assert(x0 == x1)
        assert(length(new_inq_g) == 0 )
        assert(length(new_eq_g) == 1)
        new_eq_g[1] = x11 + x22
    elseif nodeid == 1
        x33 = x1[1]
        x44 = x1[2]
        assert(length(new_inq_g) == 1 )
        assert(length(new_eq_g) == 0)
        new_inq_g[1] = x22 + x33 + x44
    elseif nodeid == 2
        x55 = x1[1]
        x66 = x1[2]
        assert(length(new_inq_g) == 1 )
        assert(length(new_eq_g) == 0)
        new_inq_g[1] = x11 + x55 + x66
    else
        assert(false)
    end
end

function str_eval_grad_f(rowid,colid,x0,x1,new_grad_f)
    # @show rowid,colid
    # @show x0
    # @show x1
    x11 = x0[1]
    x22 = x0[2]
    if (0,0) == (rowid, colid)
        assert(length(new_grad_f) == 2)
        assert(x1 == x0)
        new_grad_f[1] = 2.0*x11
        new_grad_f[2] = 2.0*x22
     elseif (1,1) == (rowid, colid)
        assert(length(new_grad_f) == 2)
        x33 = x1[1]
        x44 = x1[2]
        new_grad_f[1] = 2.0*x33
        new_grad_f[2] = 2.0*x44
    elseif (2,2) == (rowid, colid)
        assert(length(new_grad_f) == 2)
        x55 = x1[1]
        x66 = x1[2]
        new_grad_f[1] = 2.0*x55
        new_grad_f[2] = 2.0*x66
    elseif (1,0) == (rowid, colid)
        assert(length(new_grad_f) == 2)
        new_grad_f[1] = 0.0
        new_grad_f[2] = 0.0
    elseif (2,0) == (rowid, colid)
        assert(length(new_grad_f) == 2)
        new_grad_f[1] = 0.0
        new_grad_f[2] = 0.0
    else
        assert(false)
    end
end

function str_eval_jac_g(rowid,colid,x0,x1,mode,e_rowidx,e_colptr,e_values,i_rowidx,i_colptr,i_values)
    if(mode == :Structure)
        if (rowid,colid) == (0,0)
            return (2,0)
        elseif (rowid, colid) == (1,1)
            return (0,2)
        elseif (rowid, colid) == (2,2)
            return (0,2)
        elseif (rowid, colid) == (1,0)
            return (0,1)
        elseif (rowid, colid) == (2,0)
            return (0,1)
        else
            assert(false)
        end
    else
        x11 = x0[1]
        x22 = x0[2]
        if (rowid,colid) == (0,0)
            assert(length(e_rowidx) == length(e_values) == 2)
            assert(length(i_rowidx) == length(i_values) == 0)
            e_rowidx[1] = 1
            e_rowidx[2] = 1
            e_colptr[1] = 1
            e_colptr[2] = 2
            e_colptr[3] = 3
            e_values[1] = 1.0
            e_values[2] = 1.0
        elseif (rowid, colid) == (1,1)
            assert(length(e_rowidx) == length(e_values) == 0)
            assert(length(i_rowidx) == length(i_values) == 2)
            i_rowidx[1] = 1
            i_rowidx[2] = 1
            i_colptr[1] = 1
            i_colptr[2] = 2
            i_colptr[3] = 3
            i_values[1] = 1.0
            i_values[2] = 1.0
        elseif (rowid, colid) == (2,2)
            assert(length(e_rowidx) == length(e_values) == 0)
            assert(length(i_rowidx) == length(i_values) == 2)
            i_rowidx[1] = 1
            i_rowidx[2] = 1
            i_colptr[1] = 1
            i_colptr[2] = 2
            i_colptr[3] = 3
            i_values[1] = 1.0
            i_values[2] = 1.0
        elseif (rowid, colid) == (1,0)
            assert(length(e_rowidx) == length(e_values) == 0)
            assert(length(i_rowidx) == length(i_values) == 1)
            i_rowidx[1] = 1
            i_colptr[1] = 1
            i_colptr[2] = 1
            i_colptr[3] = 2
            i_values[1] = 1.0
        elseif (rowid, colid) == (2,0)
            assert(length(e_rowidx) == length(e_values) == 0)
            assert(length(i_rowidx) == length(i_values) == 1)
            i_rowidx[1] = 1
            i_colptr[1] = 1
            i_colptr[2] = 2
            i_colptr[3] = 2
            i_values[1] = 1.0
        else
            assert(false)
        end
        convert_to_c_idx(i_rowidx)
        convert_to_c_idx(i_colptr)
        convert_to_c_idx(e_rowidx)
        convert_to_c_idx(e_colptr)
    end
end

function str_eval_h(rowid,colid,x0,x1,obj_factor,lambda,mode,rowidx,colptr,values)
    # @show obj_factor
    if(mode == :Structure)
        if(rowid,colid) == (0,0)
            return 2
        elseif (rowid,colid) == (1,1)
            return 2
        elseif (rowid,colid) == (2,2)
            return 2
        elseif (rowid,colid) == (1,0)
            return 2
        elseif (rowid,colid) == (2,0)
            return 2
        elseif (rowid,colid) == (0,1)
            return 0
        elseif (rowid,colid) == (0,2)
            return 0
        else
            assert(false)
        end
    else
        assert(length(lambda) == 1)
        x11 = x0[1]
        x22 = x0[2]
        # @show (rowid,colid)
        if(rowid,colid) == (0,0)
            # @show "diag 0,0 "
            colptr[1] = 1
            colptr[2] = 2
            colptr[3] = 3
            rowidx[1] = 1
            rowidx[2] = 2
            values[1] = 2.0 * obj_factor
            values[2] = 2.0 * obj_factor
        elseif (rowid,colid) == (1,1)
            # @show "1,1"
            colptr[1] = 1
            colptr[2] = 2
            colptr[3] = 3
            rowidx[1] = 1
            rowidx[2] = 2
            values[1] = 2.0 * obj_factor
            values[2] = 2.0 * obj_factor
        elseif (rowid,colid) == (2,2)
            # @show "2,2"
            colptr[1] = 1
            colptr[2] = 2
            colptr[3] = 3
            rowidx[1] = 1
            rowidx[2] = 2
            values[1] = 2.0 * obj_factor
            values[2] = 2.0 * obj_factor
        elseif (rowid,colid) == (1,0) #the diagnal contribution to 1st stage
            # @show "diag  contr - 1, 0 "
            colptr[1] = 1
            colptr[2] = 2
            colptr[3] = 3
            rowidx[1] = 1
            rowidx[2] = 2
            values[1] = 0.0
            values[2] = 0.0
        elseif (rowid,colid) == (2,0) #the diagnal contribution  to 1st stage
            # @show "diag contr - 2, 0 "
            colptr[1] = 1
            colptr[2] = 2
            colptr[3] = 3
            rowidx[1] = 1
            rowidx[2] = 2
            values[1] = 0.0
            values[2] = 0.0
        elseif (rowid,colid) == (0,1) #liking border
            # @show "linking border - 0 ,1 "
            # colptr[1] = 1
            # colptr[2] = 2
            # colptr[3] = 3
            # rowidx[1] = 1
            # rowidx[2] = 1
            # values[1] = 1.0 * obj_factor + 1.0 * lambda[1]
            # values[2] = 1.0 * obj_factor
        elseif (rowid,colid) == (0,2) #liking border
            # @show "linking border - ", (rowid,colid)
            # colptr[1] = 1
            # colptr[2] = 2
            # colptr[3] = 3
            # rowidx[1] = 1
            # rowidx[2] = 1
            # values[1] = 1.0 * obj_factor + 1.0 * lambda[1]
            # values[2] = 1.0 * obj_factor
        else
            assert(false)
        end
        convert_to_c_idx(colptr)
        convert_to_c_idx(rowidx)
        # @show colptr
        # @show rowidx
        # @show values
    end
end

function convert_to_c_idx(indicies)
    for i in 1:length(indicies)
        indicies[i] = indicies[i] - 1
    end
end

function create()
    comm = MPI.COMM_WORLD
    println("[$(MPI.Comm_rank(comm))/$(MPI.Comm_size(comm))] create problem ")
    # @show comm
    model = FakeModel(:Min,0, 2,
        str_init_x0, str_prob_info, str_eval_f, str_eval_g, str_eval_grad_f, str_eval_jac_g, str_eval_h)
    prob = createProblemStruct(comm, model)
    # println("end create problem ")
    return prob
end

function solve(prob)
    println("Solve problem")
    # @show prob
    ret = solveProblemStruct(prob)
    # println("end solve problem")
end

function main()
    MPI.Init()

    prob = create()

    solve(prob)

    MPI.Finalize()
end

main()

max_stars_repo_name: jalving/PIPS

max_issues_repo_name: jalving/PIPS

max_forks_repo_name: jalving/PIPS


text:
using Test, MutableNamedTuples

@testset "tests" begin
    mnt = MutableNamedTuple(a=1, b="hi")
    @test mnt isa MutableNamedTuple
    mnt.a = 2

    @test mnt.a == 2
    @test NamedTuple(mnt) == (;a=2, b="hi")
    @test collect(mnt) == [2;"hi"]
    @test length(mnt) == 2
    @test mnt[1] == 2
    @test mnt[2] == "hi"
    @test mnt[:a] == 2
    @test mnt[:b] == "hi"

    mnt2 = MutableNamedTuple{(:a,:b)}((2,"hi"))
    @test NamedTuple(mnt2) == NamedTuple(mnt)
end


text:
using Test
using Distances: Euclidean, CosineDist
using Random
using SparseArrays
using LinearAlgebra
using UMAP
using UMAP: initialize_embedding, fuzzy_simplicial_set, compute_membership_strengths, smooth_knn_dists, smooth_knn_dist, spectral_layout, optimize_embedding, knn_search, combine_fuzzy_sets, fit_ab, SMOOTH_K_TOLERANCE

include("utils_tests.jl")
include("umap_tests.jl")


text:
include("graphs.jl")
using DataFrames
using Plots
using Statistics
using HypothesisTests
using CSV
using Weave
using IJulia
using Random

function benchmark(exp::Experiment)
    workspace = exp.workspace

    if !ispath(workspace * "/notebooks")
        mkpath(workspace * "/notebooks")
    end

    convert_doc(
        "./src/benchmark/notebooks/overview.jmd",
        workspace * "/notebooks/overview.ipynb",
    )

    IJulia.jupyterlab(dir = workspace * "/notebooks")
end

function benchmark(path::String = "_results")
    df = DataFrames.DataFrame()
    for item in readdir(path)
        println(item)
        maindf = DataFrame(CSV.File(path * '/' * item))
        rename!(maindf, :iteration => :runtime, :size => :problem_size)
        df = vcat(df, maindf[:, [:algorithm, :runtime, :problem_size]])
    end
    resultpvalue(df)
    return benchmark(df)
end

function benchmark(traces::Array{Trace,1})
    df = tracestodataframe(traces)
    resultpvalue(df)
    return benchmark(df)
end

function tracestodataframe(traces::Array{Trace,1})
    result = DataFrames.DataFrame()
    for trace in traces
        push!(
            result,
            (
                algorithm = string(typeof(trace.algorithm)),
                runtime = trace.optimum[2],
                problem_size = length(trace.individual.problem),
            ),
        )
    end
    result
end

function jobstodataframe(jobs::Array{Job,1}) end

function benchmark(data::DataFrames.DataFrame)
    algorithms = keys(groupby(data, "algorithm"))
    plts = []

    append!(plts, [plot()])
    for algo in algorithms
        classical_0 = filter(row -> row["algorithm"] == algo[1], data)
        classical_0 = combine(groupby(classical_0, :problem_size), :runtime => mean)
        y = classical_0[!, :runtime_mean]
        x = classical_0[!, :problem_size]
        # v = classical_0[!,4]
        plot!(
            plts[end],
            x,
            y,
            fillalpha = 0.3,
            formatter = :plain,
            dpi = 300,
            fontfamily = "Times New Roman Bold",
            legend = :right,
            linewidth = 2,
            markersize = 3,
            markershape = :circle,
            linealpha = 0.5,
            size = (600, 300),
            label = algo,
        )
    end
    return plts[1]
end

function resultpvalue(df::DataFrames.DataFrame)
    algo = keys(groupby(df, :algorithm))
    problem_size = keys(groupby(df, :ProblemParam_n))
    allowmissing!(df)

    for item in problem_size
        println("Problem Size: ", item["ProblemParam_n"])
        data2 = groupby(
            filter(row -> row["ProblemParam_n"] == item["ProblemParam_n"] && !ismissing(row["runtime"]), df),
            :algorithm,
        )
        for item in data2
            for item2 in data2
                p_value = pvalue(MannWhitneyUTest(item[!, "runtime"], item2[!, "runtime"]))
                println(item[1][1], " vs ", item2[1][1], ", Pvalue = ", p_value)
            end
        end
    end
end

function runtimes(jobs::Array{Job,1})
    result = DataFrames.DataFrame()
    allowmissing!(result)
    for job in jobs
        trace = job.trace
        if !isnothing(trace) && job.status == finished
            thing = (
                seed = trace.seed,
                algorithm = name(trace.algorithm),
                problem = string(trace.individual.problem),
                runtime = trace.optimum[2],
            )
            push!(result, thing)
        end
    end

    for job in shuffle(jobs)
        trace = job.trace
        if !isnothing(trace) && job.status == finished && size(trace.rows)[1]>1
            thing = (seed = trace.seed, algorithm = name(trace.algorithm))
            paramgroups = (
                ProblemParam = deepparameters(trace.individual.problem),
                AlgorithmParams = deepparameters(trace.algorithm),
            )
            for paramgroupkey in keys(paramgroups)
                for p in paramgroups[paramgroupkey]
                    label = Symbol(string(paramgroupkey) * "_" * string(p[1]))
                    if string(label) in names(result)
                        result[result[!,:seed].==trace.seed, label] .= p[2]
                    else
                        df = DataFrame([(seed = trace.seed, l = p[2])])
                        rename!(df, :l => label)
                        result = outerjoin(result, df, on = [:seed], makeunique = true)
                    end
                end
            end
        end
    end

    result
end

runtimes(exp::Experiment) = runtimes(exp.jobs)

max_stars_repo_name: DTUComputeTONIA/SDRLS

max_issues_repo_name: DTUComputeTONIA/SDRLS

max_forks_repo_name: DTUComputeTONIA/SDRLS


text:
using JLD2
using Statistics

mutable struct TTMazeErrorMap <: LoggerKeyData
    error_map::Array{Any}
    eval_set::Dict
    log_interval::Int
    est::Array{Any}

    function TTMazeErrorMap(logger_init_info)
        # eval_set = @load string(pwd(),"/src/data/TTMazeEvalSet.jld2") TTMazeEvalSet
        eval_set = @load string(pwd(),"/src/data/TTMazeUniformEvalSet.jld2") TTMazeUniformEvalSet
        num_logged_steps = fld(logger_init_info[LoggerInitKey.TOTAL_STEPS], logger_init_info[LoggerInitKey.INTERVAL])

        new([], TTMazeUniformEvalSet, logger_init_info[LoggerInitKey.INTERVAL], [])
    end
end

function lg_step!(self::TTMazeErrorMap, env, agent, s, a, s_next, r, is_terminal, cur_step_in_episode, cur_step_total)
    if rem(cur_step_total, self.log_interval) == 0
        ind = fld(cur_step_total, self.log_interval)
        Q_est = hcat([get_demon_prediction(agent, state, Int(action)) for (state, action) in zip(self.eval_set["states"], self.eval_set["actions"])]...)
        true_values = TabularTMazeUtils.get_true_values(env, self.eval_set["ests"])
        err = (Q_est - true_values) .^ 2
        push!(self.error_map, err)
        push!(self.est, Q_est)
    end
end

function lg_episode_end!(self::TTMazeErrorMap, cur_step_in_episode, cur_step_total)
end

function save_log(self::TTMazeErrorMap, save_dict::Dict)
    save_dict[:ttmaze_error_map] = self.error_map
    save_dict[:ttmaze_error_map_eval_set] = self.eval_set
    save_dict[:ttmaze_error_map_est] = self.est
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
module WinApi

export IMAGE_BITMAP, IMAGE_ICON, IMAGE_CURSOR
export LR_CREATEDIBSECTION, LR_DEFAULTCOLOR, LR_DEFAULTSIZE, LR_LOADFROMFILE,
    LR_LOADMAP3DCOLORS, LR_LOADTRANSPARENT, LR_MONOCHROME, LR_SHARED, LR_VGACOLOR
export WM_CLOSE, WM_DESTROY, WM_LBUTTONDOWN, WM_CREATE, WM_COMMAND
export WS_EX_CLIENTEDGE, WS_EX_WINDOWEDGE, WS_EX_TOOLWINDOW, WS_EX_TOPMOST,
    WS_EX_OVERLAPPEDWINDOW, WS_EX_PALETTEWINDOW
export WS_OVERLAPPED, WS_CAPTION, WS_SYSMENU, WS_THICKFRAME, WS_MINIMIZEBOX,
    WS_MAXIMIZEBOX, WS_OVERLAPPEDWINDOW
export SW_SHOWDEFAULT
export CW_USEDEFAULT
export MF_BITMAP, MF_CHECKED, MF_DISABLED, MF_ENABLED, MF_GRAYED, MF_MENUBARBREAK,
    MF_MENUBREAK, MF_OWNERDRAW, MF_POPUP, MF_SEPARATOR, MF_STRING, MF_UNCHECKED
export POINT, MSG, @winproc
export getlasterror, defwindowproc, destroywindow, post_quitmessage, registerclass,
    unregisterclass, createwindow, showwindow, updatewindow, getmessage,
    translatemessage, dispatchmessage, iswindow, createmenu, createpopupmenu,
    appendmenu, setmenu, loadimage,
    showwindow_wait, openwindow

###### Begin Constants #####

##### Image types #####
const IMAGE_BITMAP = convert(Cint, 0)
const IMAGE_ICON = convert(Cint, 1)
const IMAGE_CURSOR = convert(Cint, 2)

##### Resource load flags #####
const LR_CREATEDIBSECTION = convert(Cuint, 0x00002000)
const LR_DEFAULTCOLOR = convert(Cuint, 0x00000000)
const LR_DEFAULTSIZE = convert(Cuint, 0x00000040)
const LR_LOADFROMFILE = convert(Cuint, 0x00000010)
const LR_LOADMAP3DCOLORS = convert(Cuint, 0x00001000)
const LR_LOADTRANSPARENT = convert(Cuint, 0x00000020)
const LR_MONOCHROME = convert(Cuint, 0x00000001)
const LR_SHARED = convert(Cuint, 0x00008000)
const LR_VGACOLOR = convert(Cuint, 0x00000080)

###### Window Messages #####
const WM_CLOSE   = convert(Cuint,0x0010)
const WM_DESTROY = convert(Cuint,0x0002)
const WM_LBUTTONDOWN = convert(Cuint, 0x0201)
const WM_CREATE = convert(Cuint, 0x0001)
const WM_COMMAND = convert(Cuint, 0x0111)

##### Extended Window Styles #####
const WS_EX_CLIENTEDGE = convert(Culong,0x00000200)
const WS_EX_WINDOWEDGE = convert(Culong,0x00000100)
const WS_EX_TOOLWINDOW = convert(Culong,0x00000080)
const WS_EX_TOPMOST = convert(Culong,0x00000008)
const WS_EX_OVERLAPPEDWINDOW = WS_EX_WINDOWEDGE | WS_EX_CLIENTEDGE
const WS_EX_PALETTEWINDOW = WS_EX_WINDOWEDGE | WS_EX_TOOLWINDOW | WS_EX_TOPMOST

##### Window Styles #####
const WS_OVERLAPPED = convert(Culong,0x00000000)
const WS_CAPTION = convert(Culong,0x00C00000)
const WS_SYSMENU = convert(Culong,0x00080000)
const WS_THICKFRAME = convert(Culong,0x00040000)
const WS_MINIMIZEBOX = convert(Culong,0x00020000)
const WS_MAXIMIZEBOX = convert(Culong,0x00010000)
const WS_OVERLAPPEDWINDOW = WS_OVERLAPPED | WS_CAPTION | WS_SYSMENU | WS_THICKFRAME | WS_MINIMIZEBOX | WS_MAXIMIZEBOX

##### ShowWindow Parameters #####
const SW_SHOWDEFAULT = convert(Cint,10)

##### CreateWindow Parameters #####
const CW_USEDEFAULT = convert(Cint,-2147483648)

##### WinApi Error Codes #####
const ERROR_INVALID_WINDOW_HANDLE = convert(Culong,0x00000578)
const ERROR_CLASS_ALREADY_EXISTS = convert(Culong,0x00000582)

##### Menu creation flags #####
const MF_BITMAP = convert(Cuint, 0x00000004)
const MF_CHECKED = convert(Cuint, 0x00000008)
const MF_DISABLED = convert(Cuint, 0x00000002)
const MF_ENABLED = convert(Cuint, 0x00000000)
const MF_GRAYED = convert(Cuint, 0x00000001)
const MF_MENUBARBREAK = convert(Cuint, 0x00000020)
const MF_MENUBREAK = convert(Cuint, 0x00000040)
const MF_OWNERDRAW = convert(Cuint, 0x00000100)
const MF_POPUP = convert(Cuint, 0x00000010)
const MF_SEPARATOR = convert(Cuint, 0x00000800)
const MF_STRING = convert(Cuint, 0x00000000)
const MF_UNCHECKED = convert(Cuint, 0x00000000)

##### End Constants #####

##### Begin Structs #####

mutable struct WNDCLASSEXW
    cbSize::Cuint
    style::Cuint
    lpfnWndProc::Ptr{Cvoid}
    cbClsExtra::Cint
    cbWndExtra::Cint
    hInstance::Ptr{Cvoid}
    hIcon::Ptr{Cvoid}
    hCursor::Ptr{Cvoid}
    hbrBackground::Ptr{Cvoid}
    lpszMenuName::Cwstring
    lpszClassName::Cwstring
    hIconSm::Ptr{Cvoid}
end

mutable struct POINT
    x::Clong
    y::Clong
end

mutable struct MSG
    hwnd::Ptr{Cvoid}
    message::Cuint
    wParam::Culonglong
    lParam::Clonglong
    time::Culong
    pt::POINT
end

MSG(hwnd, msg, wParam, lParam) = MSG(hwnd, msg, wParam, lParam, 0, POINT(0, 0))
MSG() = MSG(C_NULL,0,0,0)

##### End Structs #####

macro winproc(proc)
    return :( @cfunction($proc, Clonglong, (Ptr{Cvoid},Cuint,Culonglong,Clonglong,)) )
end

function cwstring(val)
    val == C_NULL || isa(val, Cwstring) ? val :
        Base.unsafe_convert(Cwstring,Base.cconvert(Cwstring,string(val)))
end

function getlasterror()
    ccall((:GetLastError, "kernel32"), Culong, ())
end

function defwindowproc(hwnd, msg, wParam, lParam)
    ccall((:DefWindowProcW,"user32"),Clonglong,(Ptr{Cvoid},Cuint,Culonglong,Clonglong),hwnd,msg,wParam,lParam)
end

function destroywindow(hwnd)
    ret = ccall((:DestroyWindow,"user32"),Cint,(Ptr{Cvoid},),hwnd)
    error = ret == C_NULL ? getlasterror() : 0
    (ret, error)
end

function post_quitmessage(exitcode)
    ccall((:PostQuitMessage,"user32"),Cvoid,(Cint,),exitcode)
end

function windowproc(hwnd, msg, wParam, lParam)
    if msg == WM_CLOSE
        destroywindow(hwnd)
    elseif msg == WM_DESTROY
        post_quitmessage(0)
    else
        return defwindowproc(hwnd, msg, wParam, lParam)
    end
    return 0
end

windowproc_c = @winproc(windowproc)

function registerclass(classname;
    proc = C_NULL,
    instance = C_NULL,
    icon = C_NULL,
    cursor = C_NULL,
    bg = C_NULL,
    menu = C_NULL,
    iconSm = C_NULL)

    wcex = WNDCLASSEXW(
        sizeof(WNDCLASSEXW),
        0,
        proc == C_NULL ? windowproc_c : proc,
        0,
        0,
        instance,
        icon,
        cursor,
        bg,
        cwstring(menu),
        cwstring(classname),
        iconSm
    )

    atom = ccall((:RegisterClassExW, "user32"), Cushort, (Ptr{WNDCLASSEXW},), pointer_from_objref(wcex))
    error = atom == 0 ? getlasterror() : 0

    (atom, error)
end

function unregisterclass(classname)
    ret = ccall((:UnregisterClassW, "user32"), Cint, (Cwstring, Ptr{Cvoid}), classname, C_NULL)
    error = ret == 0 ? getlasterror() : 0
    (ret, error)
end

function createwindow(;
    classname = :WC_DEFAULT,
    title = "Untitled",
    x = CW_USEDEFAULT,
    y = CW_USEDEFAULT,
    width = CW_USEDEFAULT,
    height = CW_USEDEFAULT,
    exStyle = 0,
    style = WS_OVERLAPPEDWINDOW)

    # If the window is created with the default class, that class is auto-registered
    if classname == :WC_DEFAULT
        (_, err) = registerclass(:WC_DEFAULT)
        if err != 0 && err != ERROR_CLASS_ALREADY_EXISTS
            return (0, err)
        end
    end

    hwnd = ccall(
        (:CreateWindowExW, "user32"),
        Ptr{Cvoid},
        (Culong, Cwstring, Cwstring, Culong, Cint, Cint, Cint, Cint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
        convert(Culong,exStyle),
        string(classname),
        string(title),
        convert(Culong,style),
        convert(Cint,x),
        convert(Cint,y),
        convert(Cint,width),
        convert(Cint,height),
        C_NULL,
        C_NULL,
        C_NULL,
        C_NULL)

    error = hwnd == C_NULL ? getlasterror() : 0

    (hwnd, error)
end

function showwindow(hwnd, nCmdShow = SW_SHOWDEFAULT)
    ccall((:ShowWindow, "user32"), Cint, (Ptr{Cvoid}, Cint), hwnd, nCmdShow)
end

function updatewindow(hwnd)
    ccall((:UpdateWindow, "user32"), Cint, (Ptr{Cvoid},), hwnd)
end

function getmessage(ptr_msg, hwnd = C_NULL)
    ccall((:GetMessageW,"user32"),Cint,(Ptr{MSG},Ptr{Cvoid},Cuint,Cuint),ptr_msg,C_NULL,0,0)
end

function translatemessage(ptr_msg)
    ccall((:TranslateMessage,"user32"),Cint,(Ptr{MSG},),ptr_msg)
end

function dispatchmessage(ptr_msg)
    ccall((:DispatchMessageW,"user32"),Clonglong,(Ptr{MSG},),ptr_msg)
end

function iswindow(hwnd)
    ret = ccall((:IsWindow, "user32"), Cint, (Ptr{Cvoid},), hwnd)
    ret != 0
end

function createmenu()
    ret = ccall((:CreateMenu, "user32"), Ptr{Cvoid}, ())
    err = ret == C_NULL ? getlasterror() : 0
    (ret, err)
end

function createpopupmenu()
    ret = ccall((:CreatePopupMenu, "user32"), Ptr{Cvoid}, ())
    err = ret == C_NULL ? getlasterror() : 0
    (ret, err)
end

function appendmenu(hmenu, flags, idNewItem, newItem)
    ret = ccall((:AppendMenuW, "user32"), Cint,
        (Ptr{Cvoid}, Cuint, Culonglong, Cwstring),
        hmenu, flags, idNewItem, newItem)
    err = ret == 0 ? getlasterror() : 0
    (ret, err)
end

function setmenu(hwnd, hmenu)
    ret = ccall((:SetMenu, "user32"), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), hwnd, hmenu)
    err = ret == 0 ? getlasterror() : 0
    (ret, err)
end

function loadimage(name, imgType, cx, cy, loadFlags)
    handle = ccall((:LoadImageW, "user32"), Ptr{Cvoid},
        (Ptr{Cvoid}, Cwstring, Cuint, Cint, Cint, Cuint),
        C_NULL, name, imgType, cx, cy, loadFlags)

    error = handle == C_NULL ? getlasterror() : 0

    (handle, error)
end

function showwindow_wait(hwnd)
    if !iswindow(hwnd)
        return (0, ERROR_INVALID_WINDOW_HANDLE)
    end

    showwindow(hwnd)
    updatewindow(hwnd)

    msg = MSG()
    ptr_msg = pointer_from_objref(msg)
    while getmessage(ptr_msg) > 0
        translatemessage(ptr_msg)
        dispatchmessage(ptr_msg)
    end
    return (msg.wParam, 0)
end

function openwindow(;kwargs...)
    (hwnd, err) = createwindow(;kwargs...)

    if err != 0
        return (hwnd, err)
    end

    showwindow_wait(hwnd)
end

end


text:
module POMDPDistributions

using Compat

import Base: size, eltype, length, full, convert, show, getindex, scale, rand, rand!
import Base: +, -, .+, .-

export
     # generic types
    VariateForm,
    ValueSupport,
    Univariate,
    Multivariate,
    Matrixvariate,
    Discrete,
    Continuous,
    Sampleable,
    Distribution,
    UnivariateDistribution,
    MultivariateDistribution,
    MatrixDistribution,
    NoncentralHypergeometric,
    NonMatrixDistribution,
    DiscreteDistribution,
    ContinuousDistribution,
    DiscreteUnivariateDistribution,
    DiscreteMultivariateDistribution,
    DiscreteMatrixDistribution,
    ContinuousUnivariateDistribution,
    ContinuousMultivariateDistribution,
    ContinuousMatrixDistribution,
    SufficientStats,
    AbstractMvNormal,
    AbstractMixtureModel,
    UnivariateMixture,
    MultivariateMixture,

    # distribution types
    Categorical,

    # methods
    rand,
    set_prob!

### source files

# type system
include("common.jl")

# implementation helpers
include("utils.jl")

# specific samplers and distributions
include("univariates.jl")

end # module


text:
#=
Code generation functions for Julia.
=#

# Extract the input args. This must match the arguments passed by the solver.
function handle_input_args_dg_julia(lorr, vors)
    code =
"var =       args[1];
eid =       args[2];
fid =       args[3];
grid =      args[4];
geo_facs =  args[5];
refel =     args[6]
time =      args[7];
dt =        args[8]
"
    # A trick for uniform grids to avoid repeated work
    if vors == "volume" && config.mesh_type == UNIFORM_GRID && config.geometry == SQUARE
        code *= "stiffness = args[9]; # set of stiffness matrices for each dimension\n";
        code *= "mass =      args[10]; # mass matrix\n";
    end

    return code;
end

# If needed, build derivative matrices
function build_derivative_matrices_dg_julia(lorr, vors)
    return "";
    code =
"
# Note on derivative matrices:
# RQn are vandermond matrices for the derivatives of the basis functions
# with Jacobian factors. They are made like this.
# |RQ1|   | rx sx tx || Qx |
# |RQ2| = | ry sy ty || Qy |
# |RQ3|   | rz sz tz || Qz |

"
    use_full_deriv_mat = lorr==RHS;
    if config.dimension == 1
        if vors == "volume"
            code *= "(RQ1, RD1) = build_deriv_matrix(refel, J);\n";
            code *= "TRQ1 = RQ1';\n"
        else
            code *= "RQ1_1 = refel.surf_Qr[frefelind[1]] .* vol_J1.rx[1]";
            code *= "RQ1_2 = refel.surf_Qr[frefelind[2]] .* vol_J2.rx[1]";
            code *= "TRQ1_1 = RQ1_1'";
            code *= "TRQ1_2 = RQ1_2'";
        end

    elseif config.dimension == 2
        if vors == "volume"
            code *= "(RQ1, RQ2, RD1, RD2) = build_deriv_matrix(refel, J);\n";
            code *= "(TRQ1, TRQ2) = (RQ1', RQ2');\n"
        else
            code *= "(RQ1_1,RQ2_1,RD1_1,RD2_1) = build_face_deriv_matrix(refel, frefelind[1], vol_J1, $use_full_deriv_mat)";
            code *= "(RQ1_2,RQ2_2,RD1_2,RD2_2) = build_face_deriv_matrix(refel, frefelind[2], vol_J2, $use_full_deriv_mat)";
            code *= "(TRQ1_1,TRQ1_2) = (RQ1_1',RQ1_2')";
            code *= "(TRQ2_1,TRQ2_2) = (RQ2_1',RQ2_2')";
        end

    elseif config.dimension == 3
        if vors == "volume"
            code *= "(RQ1, RQ2, RQ3, RD1, RD2, RD3) = build_deriv_matrix(refel, J);\n";
            code *= "(TRQ1, TRQ2, TRQ3) = (RQ1', RQ2', RQ3');\n"
        else
            code *= "(RQ1_1,RQ2_1,RQ3_1,RD1_1,RD2_1,RD3_1) = build_face_deriv_matrix(refel, frefelind[1], vol_J1, $use_full_deriv_mat)";
            code *= "(RQ1_2,RQ2_2,RQ3_2,RD1_2,RD2_2,RD3_2) = build_face_deriv_matrix(refel, frefelind[2], vol_J2, $use_full_deriv_mat)";
            code *= "(TRQ1_1,TRQ2_1,TRQ3_1) = (RQ1_1',RQ2_1',RQ3_1')";
            code *= "(TRQ1_2,TRQ2_2,TRQ3_2) = (RQ1_2',RQ2_2',RQ3_2')";
        end
    end

    return code;
end

# Allocate, compute, or fetch all needed values
function prepare_needed_values_dg_julia(entities, var, lorr, vors)
    # Only gather the pieces that are used. See the end of this function.
    if vors == "volume"
        # el, loc2glb, nodex, detj, J
        piece_needed = zeros(Bool, 5);
    else # surface
        # (els, loc2glb, nodex, frefelind, face2glb, facex, normal, fdetJ, vol_J)
        piece_needed = zeros(Bool, 9);
    end
    need_deriv_matrix = false;

    code = "";
    for i=1:length(entities)
        cname = make_entity_name(entities[i]);
        side1_entity = copy(entities[i]);
        side2_entity = copy(entities[i]);
        push!(side1_entity.flags, "DGSIDE1");
        push!(side2_entity.flags, "DGSIDE2");
        cnameside1 = make_entity_name(side1_entity);
        cnameside2 = make_entity_name(side2_entity);
        if is_test_function(entities[i])
            # Assign it a transpose quadrature matrix
            if vors == "volume"
                if length(entities[i].derivs) > 0
                    xyzchar = ["x","y","z"];
                    for di=1:length(entities[i].derivs)
                        code *= cname * " = TRQ"*string(entities[i].derivs[di])*"; # d/d"*xyzchar[entities[i].derivs[di]]*" of test function\n";
                    end
                    need_deriv_matrix = true;
                else
                    code *= cname * " = refel.Q'; # test function.\n";
                end
            else
                dgside = 0;
                if "DGSIDE1" in entities[i].flags
                    dgside = 1;
                elseif "DGSIDE2" in entities[i].flags
                    dgside = 2;
                end
                if length(entities[i].derivs) > 0
                    xyzchar = ["x","y","z"];
                    for di=1:length(entities[i].derivs)
                        if dgside == 0
                            code *= cnameside1 * " = TRQ"*string(entities[i].derivs[di])*"_1; # d/d"*xyzchar[entities[i].derivs[di]]*" of test function on side 1\n";
                            code *= cnameside2 * " = TRQ"*string(entities[i].derivs[di])*"_2; # d/d"*xyzchar[entities[i].derivs[di]]*" of test function on side 2\n";
                        elseif dgside == 1
                            code *= cname * " = TRQ"*string(entities[i].derivs[di])*"_1; # d/d"*xyzchar[entities[i].derivs[di]]*" of test function on side 1\n";
                        elseif dgside == 2
                            code *= cname * " = TRQ"*string(entities[i].derivs[di])*"_2; # d/d"*xyzchar[entities[i].derivs[di]]*" of test function on side 2\n";
                        end
                        need_deriv_matrix = true;
                    end
                else
                    if dgside == 0
                        code *= cnameside1 * " = (refel.surf_Q[frefelind[1]])'; # test function side 1.\n";
                        code *= cnameside2 * " = (refel.surf_Q[frefelind[2]])'; # test function side 2.\n";
                    elseif dgside == 1
                        code *= cname * " = (refel.surf_Q[frefelind[1]])'; # test function side 1.\n";
                    elseif dgside == 2
                        code *= cname * " = (refel.surf_Q[frefelind[2]])'; # test function side 2.\n";
                    end
                end
            end

        elseif is_unknown_var(entities[i], var) && lorr == LHS
            if vors == "volume"
                if length(entities[i].derivs) > 0
                    xyzchar = ["x","y","z"];
                    for di=1:length(entities[i].derivs)
                        code *= cname * " = RQ"*string(entities[i].derivs[di])*"; # d/d"*xyzchar[entities[i].derivs[di]]*" of trial function\n";
                    end
                    need_deriv_matrix = true;
                else
                    code *= cname * " = refel.Q; # trial function.\n";
                end
            else
                dgside = 0;
                if "DGSIDE1" in entities[i].flags
                    dgside = 1;
                elseif "DGSIDE2" in entities[i].flags
                    dgside = 2;
                end
                if length(entities[i].derivs) > 0
                    xyzchar = ["x","y","z"];
                    for di=1:length(entities[i].derivs)
                        if dgside == 0
                            code *= cnameside1 * " = RQ"*string(entities[i].derivs[di])*"_1; # d/d"*xyzchar[entities[i].derivs[di]]*" of trial function on side 1\n";
                            code *= cnameside2 * " = RQ"*string(entities[i].derivs[di])*"_2; # d/d"*xyzchar[entities[i].derivs[di]]*" of trial function on side 2\n";
                        elseif dgside == 1
                            code *= cname * " = RQ"*string(entities[i].derivs[di])*"_1; # d/d"*xyzchar[entities[i].derivs[di]]*" of trial function on side 1\n";
                        elseif dgside == 2
                            code *= cname * " = RQ"*string(entities[i].derivs[di])*"_2; # d/d"*xyzchar[entities[i].derivs[di]]*" of trial function on side 2\n";
                        end
                        need_deriv_matrix = true;
                    end
                else
                    if dgside == 0
                        code *= cnameside1 * " = refel.surf_Q[frefelind[1]]; # trial function side 1.\n";
                        code *= cnameside2 * " = refel.surf_Q[frefelind[2]]; # trial function side 2.\n";
                    elseif dgside == 1
                        code *= cname * " = refel.surf_Q[frefelind[1]]; # trial function side 1.\n";
                    elseif dgside == 2
                        code *= cname * " = refel.surf_Q[frefelind[2]]; # trial function side 2.\n";
                    end
                end
            end

        else
            # Is coefficient(number or function) or variable(array)?
            (ctype, cval) = get_coef_val(entities[i]);
            if ctype == -1
                # It was a special symbol like dt
                if entities[i].name == "FACENORMAL1"
                    code *= cname * " = normal["*string(entities[i].index)*"]; # normal vector component\n"
                    piece_needed[7] = true;
                else entities[i].name == "FACENORMAL2"
                    code *= cname * " = -normal["*string(entities[i].index)*"]; # reverse normal vector component\n"
                    piece_needed[7] = true;
                end
            elseif ctype == 0
                # It was a number, do nothing?
            elseif ctype == 1 # a constant wrapped in a coefficient
                # This generates something like: coef_k_i = 4;
                if length(entities[i].derivs) > 0
                    code *= cname * " = 0; # NOTE: derivative applied to constant coefficient = 0\n";
                else
                    code *= cname * " = " * string(cval) * ";\n";
                end

            elseif ctype == 2 # a coefficient function
                # This generates something like:
                ######################################
                # coef_n_i = zeros(refel.Np);
                # for coefi = 1:refel.Np
                #     coef_k_i[coefi] = (Finch.genfunctions[cval]).func(x[1,coefi], x[2,coefi],x[3,coefi],time);
                # end
                ######################################
                # Note: for LHS, this is only Nfp values. For RHS, this is Np values, but only evaluated at face
                # cargs = "(nodex[coefi], 0, 0, time)";
                # if config.dimension == 2
                #     cargs = "(nodex[1, coefi], nodex[2, coefi], 0, time)";
                # elseif config.dimension == 3
                #     cargs = "(nodex[1, coefi], nodex[2, coefi], nodex[3, coefi], time)";
                # end
                # cargs1 = "(nodex[1][coefi], 0, 0, time)";
                # if config.dimension == 2
                #     cargs1 = "(nodex[1][1, coefi], nodex[1][2, coefi], 0, time)";
                # elseif config.dimension == 3
                #     cargs1 = "(nodex[1][1, coefi], nodex[1][2, coefi], nodex[1][3, coefi], time)";
                # end
                # cargs2 = "(nodex[2][coefi], 0, 0, time)";
                # if config.dimension == 2
                #     cargs2 = "(nodex[2][1, coefi], nodex[2][2, coefi], 0, time)";
                # elseif config.dimension == 3
                #     cargs2 = "(nodex[2][1, coefi], nodex[2][2, coefi], nodex[2][3, coefi], time)";
                # end
                # cargsf = "(facex[coefi], 0, 0, time)";
                # if config.dimension == 2
                #     cargsf = "(facex[1, coefi], facex[2, coefi], 0, time)";
                # elseif config.dimension == 3
                #     cargsf = "(facex[1, coefi], facex[2, coefi], facex[3, coefi], time)";
                # end

                nodesymbol = "nodex";
                nodesymbol1 = "nodex[1]";
                nodesymbol2 = "nodex[2]";
                nodesymbolf = "facex";
                coef_index = get_coef_index(entities[i]);

                if vors == "volume"
                    code *= cname * " = zeros(refel.Np);\n";
                    # code *= "for coefi = 1:refel.Np " * cname * "[coefi] = (Finch.genfunctions["*string(cval)*"]).func" * cargs * " end\n";
                    code *= "for coefi = 1:refel.Np " * cname * "[coefi] = evaluate_coefficient(Finch.coefficients["*string(coef_index)*"], "*string(entities[i].index)*", "*nodesymbol*"[:,coefi], time, eid, fid) end\n";
                    # Apply any needed derivative operators. Interpolate at quadrature points.
                    if length(entities[i].derivs) > 0
                        xyzchar = ["x","y","z"];
                        for di=1:length(entities[i].derivs)
                            code *= cname * " = RQ"*string(entities[i].derivs[di])*" * " * cname *
                                    "; # Apply d/d"*xyzchar[entities[i].derivs[di]]*" and interpolate at quadrature points.\n";
                        end
                        need_deriv_matrix = true;
                    else
                        code *= cname * " = refel.Q * " * cname * "; # Interpolate at quadrature points.\n";
                    end
                    piece_needed[3] = true;

                else # surface
                    dgside = 0;
                    if "DGSIDE1" in entities[i].flags
                        dgside = 1;
                    elseif "DGSIDE2" in entities[i].flags
                        dgside = 2;
                    end
                    # If derivatives are needed, must evaluate at all volume nodes.
                    if length(entities[i].derivs) > 0
                        code *= cname * " = zeros(refel.Np);\n";
                        if dgside == 2
                            # code *= "for coefi = 1:refel.Np " * cname * "[coefi] = (Finch.genfunctions["*string(cval)*"]).func" * cargs2 * " end\n";
                            code *= "for coefi = 1:refel.Np " * cname * "[coefi] = evaluate_coefficient(Finch.coefficients["*string(coef_index)*"], "*string(entities[i].index)*", "*nodesymbol2*"[:,coefi], time, eid, fid) end\n";
                        else
                            # code *= "for coefi = 1:refel.Np " * cname * "[coefi] = (Finch.genfunctions["*string(cval)*"]).func" * cargs1 * " end\n";
                            code *= "for coefi = 1:refel.Np " * cname * "[coefi] = evaluate_coefficient(Finch.coefficients["*string(coef_index)*"], "*string(entities[i].index)*", "*nodesymbol1*"[:,coefi], time, eid, fid) end\n";
                        end
                        xyzchar = ["x","y","z"];
                        for di=1:length(entities[i].derivs)
                            if dgside == 2
                                code *= cname * " = RD"*string(entities[i].derivs[di])*"_2 * " * cname *
                                    "; # Apply d/d"*xyzchar[entities[i].derivs[di]]*".\n";
                            else
                                code *= cname * " = RD"*string(entities[i].derivs[di])*"_1 * " * cname *
                                    "; # Apply d/d"*xyzchar[entities[i].derivs[di]]*".\n";
                            end
                        end
                        if dgside == 2
                            code *= cname * " = " * cname * "[refel.face2local[frefelind[2]]; # extract face values only.";
                        else
                            code *= cname * " = " * cname * "[refel.face2local[frefelind[1]]; # extract face values only.";
                        end
                        need_deriv_matrix = true;
                        piece_needed[3] = true;

                    else # no derivatives, only need surface nodes
                        code *= cname * " = zeros(Nfp);\n";
                        # code *= "for coefi = 1:Nfp " * cname * "[coefi] = (Finch.genfunctions["*string(cval)*"]).func" * cargsf * " end\n";
                        code *= "for coefi = 1:refel.Nfp " * cname * "[coefi] = evaluate_coefficient(Finch.coefficients["*string(coef_index)*"], "*string(entities[i].index)*", "*nodesymbolf*"[:,coefi], time, eid, fid) end\n";
                        piece_needed[6] = true;
                    end

                    # Interpolate at surface quadrature points
                    if dgside == 2
                        code *= cname * " = refel.surf_Q[frefelind[2]][:,refel.face2local[frefelind[2]]] * " * cname * "; # Interpolate at quadrature points.\n";
                    else
                        code *= cname * " = refel.surf_Q[frefelind[1]][:,refel.face2local[frefelind[1]]] * " * cname * "; # Interpolate at quadrature points.\n";
                    end
                end

            elseif ctype == 3 # a known variable value
                # This generates something like: coef_u_1 = copy((Finch.variables[1]).values[1, loc2glb])
                if vors == "volume"
                    code *= cname * " = copy((Finch.variables["*string(cval)*"]).values["*string(entities[i].index)*", loc2glb]);\n";
                    # Apply any needed derivative operators.
                    if length(entities[i].derivs) > 0
                        xyzchar = ["x","y","z"];
                        for di=1:length(entities[i].derivs)
                            code *= cname * " = RQ"*string(entities[i].derivs[di])*" * " * cname *
                                    "; # Apply d/d"*xyzchar[entities[i].derivs[di]]*"\n";
                        end
                        need_deriv_matrix = true;
                    else
                        code *= cname * " = refel.Q * " * cname * "; # Interpolate at quadrature points.\n";
                    end
                    piece_needed[2] = true;

                else # surface
                    dgside = 0;
                    if "DGSIDE1" in entities[i].flags
                        dgside = 1;
                    elseif "DGSIDE2" in entities[i].flags
                        dgside = 2;
                    end
                    # If derivatives are needed, must evaluate at all volume nodes.
                    if length(entities[i].derivs) > 0
                        if dgside == 0
                            code *= cnameside1 * " = copy((Finch.variables["*string(cval)*"]).values["*string(entities[i].index)*", loc2glb[1]]);\n";
                            code *= cnameside2 * " = copy((Finch.variables["*string(cval)*"]).values["*string(entities[i].index)*", loc2glb[2]]);\n";
                        elseif dgside == 1
                            code *= cname * " = copy((Finch.variables["*string(cval)*"]).values["*string(entities[i].index)*", loc2glb[1]]);\n";
                        else
                            code *= cname * " = copy((Finch.variables["*string(cval)*"]).values["*string(entities[i].index)*", loc2glb[2]]);\n";
                        end
                        xyzchar = ["x","y","z"];
                        for di=1:length(entities[i].derivs)
                            if dgside == 0
                                code *= cnameside1 * " = RD"*string(entities[i].derivs[di])*"_1 * " * cnameside1 *
                                    "; # Apply d/d"*xyzchar[entities[i].derivs[di]]*".\n";
                                code *= cnameside2 * " = RD"*string(entities[i].derivs[di])*"_2 * " * cnameside2 *
                                    "; # Apply d/d"*xyzchar[entities[i].derivs[di]]*".\n";
                            elseif dgside == 1
                                code *= cname * " = RD"*string(entities[i].derivs[di])*"_1 * " * cname *
                                    "; # Apply d/d"*xyzchar[entities[i].derivs[di]]*".\n";
                            elseif dgside == 2
                                code *= cname * " = RD"*string(entities[i].derivs[di])*"_2 * " * cname *
                                    "; # Apply d/d"*xyzchar[entities[i].derivs[di]]*".\n";
                            end
                        end
                        # if dgside == 0
                        #     code *= cnameside1 * " = " * cnameside1 * "[refel.face2local[frefelind[1]]; # extract face values only.";
                        #     code *= cnameside2 * " = " * cnameside2 * "[refel.face2local[frefelind[2]]; # extract face values only.";
                        # elseif dgside == 1
                        #     code *= cname * " = " * cname * "[refel.face2local[frefelind[1]]; # extract face values only.";
                        # elseif dgside == 2
                        #     code *= cname * " = " * cname * "[refel.face2local[frefelind[2]]; # extract face values only.";
                        # end
                        need_deriv_matrix = true;
                        piece_needed[2] = true;

                    else # no derivatives
                        if dgside == 0
                            code *= cnameside1 * " = copy((Finch.variables["*string(cval)*"]).values["*string(entities[i].index)*", loc2glb[1]]);\n";
                            code *= cnameside2 * " = copy((Finch.variables["*string(cval)*"]).values["*string(entities[i].index)*", loc2glb[2]]);\n";
                        elseif dgside == 1
                            code *= cname * " = copy((Finch.variables["*string(cval)*"]).values["*string(entities[i].index)*", loc2glb[1]]);\n";
                        else
                            code *= cname * " = copy((Finch.variables["*string(cval)*"]).values["*string(entities[i].index)*", loc2glb[2]]);\n";
                        end
                        # piece_needed[5] = true;
                        piece_needed[2] = true;
                    end

                    # # Interpolate at surface quadrature points
                    # if dgside == 0
                    #     code *= cnameside1 * " = refel.surf_Q[frefelind[1]][:,refel.face2local[frefelind[1]]] * " * cnameside1 * "; # Interpolate at quadrature points.\n";
                    #     code *= cnameside2 * " = refel.surf_Q[frefelind[2]][:,refel.face2local[frefelind[2]]] * " * cnameside2 * "; # Interpolate at quadrature points.\n";
                    # elseif dgside == 1
                    #     code *= cname * " = refel.surf_Q[frefelind[1]][:,refel.face2local[frefelind[1]]] * " * cname * "; # Interpolate at quadrature points.\n";
                    # else
                    #     code *= cname * " = refel.surf_Q[frefelind[2]][:,refel.face2local[frefelind[2]]] * " * cname * "; # Interpolate at quadrature points.\n";
                    # end
                end

            end
        end # if coefficient
        #code *= "println(\"did \"*string("*string(i)*"))\n";
    end # entity loop

    # Add the needed pieces
    if vors == "volume"
        # el, loc2glb, nodex, detj, J
        piece_needed[4] = true; #I'm basically always going to need this.
        needed_pieces = "";
        if piece_needed[1]
            needed_pieces *= "el = eid;\n"
        end
        if piece_needed[2] || piece_needed[3]
            needed_pieces *= "loc2glb = grid.loc2glb[:,eid];    # local to global map\n"
        end
        if piece_needed[3]
            needed_pieces *= "nodex = grid.allnodes[:,loc2glb]; # node coordinates\n"
        end
        if piece_needed[4]
            needed_pieces *= "detj = geo_factors.detJ[eid];     # geometric factors\n"
            needed_pieces *= "wdetj = refel.wg .* detj;         # quadrature weights\n"
        end
        if piece_needed[5] || need_deriv_matrix
            needed_pieces *= "J = geo_factors.J[eid];           # geometric factors\n"
        end

        if need_deriv_matrix
            needed_pieces *=
"
# Note on derivative matrices:
# RQn are quadrature matrices for the derivatives of the basis functions
# with Jacobian factors. They are made like this.
# |RQ1|   | rx sx tx || Qx |
# |RQ2| = | ry sy ty || Qy |
# |RQ3|   | rz sz tz || Qz |

"
            use_full_deriv_mat = lorr==RHS;
            if config.dimension == 1
                needed_pieces *= "(RQ1, RD1) = build_deriv_matrix(refel, J);\n";
                needed_pieces *= "TRQ1 = RQ1';\n"

            elseif config.dimension == 2
                needed_pieces *= "(RQ1, RQ2, RD1, RD2) = build_deriv_matrix(refel, J);\n";
                needed_pieces *= "(TRQ1, TRQ2) = (RQ1', RQ2');\n"

            elseif config.dimension == 3
                needed_pieces *= "(RQ1, RQ2, RQ3, RD1, RD2, RD3) = build_deriv_matrix(refel, J);\n";
                needed_pieces *= "(TRQ1, TRQ2, TRQ3) = (RQ1', RQ2', RQ3');\n"
            end
        end

    else # surface
        piece_needed[4] = true; #I'm basically always going to need this.
        piece_needed[8] = true; #I'm basically always going to need this.
        needed_pieces =
"if grid.face2element[1, fid] == eid # The normal points out of e"*(piece_needed[7] ? "\n\tnormal = grid.facenormals[:, fid];" : "")*"
    neighbor = grid.face2element[2, fid];
    frefelind = [grid.faceRefelInd[1,fid], grid.faceRefelInd[2,fid]]; # refel based index of face in both elements
else # The normal points into e"*(piece_needed[7] ? "\n\tnormal = -grid.facenormals[:, fid];" : "")*"
    neighbor = grid.face2element[1, fid];
    frefelind = [grid.faceRefelInd[2,fid], grid.faceRefelInd[1,fid]]; # refel based index of face in both elements
end
if neighbor == 0 # This is a boundary face. For now, just compute as if neighbor is identical. BCs handled later.
    neighbor = eid;
    frefelind[2] = frefelind[1];
end
Nfp = refel.Nfp[frefelind[1]]; # number of face nodes\n
";
        # (els, loc2glb, nodex, frefelind, face2glb, facex, normal, fdetJ, vol_J)
        if piece_needed[1]
            needed_pieces *= "els = (eid, neighbor); # indices of elements on both sides\n"
        end
        if piece_needed[2] || piece_needed[3]
            needed_pieces *= "loc2glb = (grid.loc2glb[:,eid], grid.loc2glb[:, neighbor]); # volume local to global\n"
        end
        if piece_needed[3]
            needed_pieces *= "nodex = (grid.allnodes[:,loc2glb[1][:]], grid.allnodes[:,loc2glb[2][:]]); # volume node coordinates\n"
        end
        if piece_needed[5] || piece_needed[6]
            needed_pieces *= "face2glb = grid.face2glb[:,:,fid];         # global index for face nodes for each side of each face\n"
            needed_pieces *=
"if neighbor == eid
    face2glb[:,2] = face2glb[:,1];
end\n"
        end
        if piece_needed[6]
            needed_pieces *= "facex = grid.allnodes[:, face2glb[:, 1]];  # face node coordinates\n"
        end
        if piece_needed[8]
            needed_pieces *= "face_detJ = geo_facs.face_detJ[fid];       # detJ on face\n"
            needed_pieces *= "wdetj = refel.surf_wg[1] .* face_detJ;     # quadrature weights\n"
        end
        if piece_needed[9] || need_deriv_matrix
            needed_pieces *= "vol_J = (geo_facs.J[eid], geo_facs.J[neighbor]);\n"
        end

        if need_deriv_matrix
            needed_pieces *=
"
# Note on derivative matrices:
# RQn are vandermond matrices for the derivatives of the basis functions
# with Jacobian factors. They are made like this.
# |RQ1|   | rx sx tx || Qx |
# |RQ2| = | ry sy ty || Qy |
# |RQ3|   | rz sz tz || Qz |

"
            use_full_deriv_mat = lorr==RHS;
            if config.dimension == 1
                needed_pieces *= "RQ1_1 = refel.surf_Qr[frefelind[1]] .* vol_J[1].rx[1]\n";
                needed_pieces *= "RQ1_2 = refel.surf_Qr[frefelind[2]] .* vol_J[2].rx[1]\n";
                needed_pieces *= "TRQ1_1 = RQ1_1'\n";
                needed_pieces *= "TRQ1_2 = RQ1_2'\n";

            elseif config.dimension == 2
                needed_pieces *= "(RQ1_1,RQ2_1,RD1_1,RD2_1) = build_face_deriv_matrix(refel, frefelind[1], vol_J[1], $use_full_deriv_mat)\n";
                needed_pieces *= "(RQ1_2,RQ2_2,RD1_2,RD2_2) = build_face_deriv_matrix(refel, frefelind[2], vol_J[2], $use_full_deriv_mat)\n";
                needed_pieces *= "(TRQ1_1,TRQ1_2) = (RQ1_1',RQ1_2')\n";
                needed_pieces *= "(TRQ2_1,TRQ2_2) = (RQ2_1',RQ2_2')\n";

            elseif config.dimension == 3
                needed_pieces *= "(RQ1_1,RQ2_1,RQ3_1,RD1_1,RD2_1,RD3_1) = build_face_deriv_matrix(refel, frefelind[1], vol_J[1], $use_full_deriv_mat)\n";
                needed_pieces *= "(RQ1_2,RQ2_2,RQ3_2,RD1_2,RD2_2,RD3_2) = build_face_deriv_matrix(refel, frefelind[2], vol_J[2], $use_full_deriv_mat)\n";
                needed_pieces *= "(TRQ1_1,TRQ2_1,TRQ3_1) = (RQ1_1',RQ2_1',RQ3_1')\n";
                needed_pieces *= "(TRQ1_2,TRQ2_2,TRQ3_2) = (RQ1_2',RQ2_2',RQ3_2')\n";
            end
        end
    end

    code = needed_pieces * "\n" * code;

    return code;
end

function make_elemental_computation_dg_julia(terms, var, dofsper, offset_ind, lorr, vors)
    # Here is where I make some assumption about the form of the expression.
    # Since it was expanded by the parser it should look like a series of terms: t1 + t2 + t3...
    # Where each term is multiplied by one test function component, and if LHS, involves one unknown component.
    # The submatrix modified by a term is determined by these, so go through the terms and divide them
    # into their submatrix expressions.
    # Each term will look something like
    # LHS: test_part * diagm(weight_part .* coef_part) * trial_part
    # RHS: test_part * (weight_part .* coef_part)
    code = "";

    # Allocate the vector or matrix to be returned if needed
    # There will be four pieces for LHS and two for RHS.
    if dofsper > 1
        if vors == "volume"
            if lorr == RHS
                code *= "element_vector = zeros(refel.Np * "*string(dofsper)*"); # Allocate the returned vector.\n"
            else
                code *= "element_matrix = zeros(refel.Np * "*string(dofsper)*", refel.Np * "*string(dofsper)*"); # Allocate the returned matrix.\n"
            end
        else # surface
            if lorr == RHS
                code *= "element_vectorL = zeros(Nfp * "*string(dofsper)*"); # Allocate return vector for element 1.\n"
                code *= "element_vectorR = zeros(Nfp * "*string(dofsper)*"); # Allocate return vector for element 2.\n"
            else
                code *= "element_matrixLL = zeros(Nfp * "*string(dofsper)*", Nfp * "*string(dofsper)*"); # Allocate return matrix for LL.\n"
                code *= "element_matrixLR = zeros(Nfp * "*string(dofsper)*", Nfp * "*string(dofsper)*"); # Allocate return matrix for LR.\n"
                code *= "element_matrixRL = zeros(Nfp * "*string(dofsper)*", Nfp * "*string(dofsper)*"); # Allocate return matrix for RL.\n"
                code *= "element_matrixRR = zeros(Nfp * "*string(dofsper)*", Nfp * "*string(dofsper)*"); # Allocate return matrix for RR.\n"
            end
        end

    else
        if vors == "volume"
            if lorr == RHS
                code *= "element_vector = zeros(refel.Np); # Allocate the returned vector.\n"
            else
                code *= "element_matrix = zeros(refel.Np, refel.Np); # Allocate the returned matrix.\n"
            end
        else # surface
            if lorr == RHS
                code *= "element_vectorL = zeros(Nfp); # Allocate return vector for element 1.\n"
                code *= "element_vectorR = zeros(Nfp); # Allocate return vector for element 2.\n"
            else
                code *= "element_matrixLL = zeros(Nfp, Nfp); # Allocate return matrix for LL.\n"
                code *= "element_matrixLR = zeros(Nfp, Nfp); # Allocate return matrix for LR.\n"
                code *= "element_matrixRL = zeros(Nfp, Nfp); # Allocate return matrix for RL.\n"
                code *= "element_matrixRR = zeros(Nfp, Nfp); # Allocate return matrix for RR.\n"
            end
        end

    end

    # Separate the factors of each term into test, trial, coef and form the calculation
    if vors == "volume"
        if lorr == LHS
            submatrices = Array{String, 2}(undef, dofsper, dofsper);
        else # RHS
            submatrices = Array{String, 1}(undef, dofsper);
        end
        for smi=1:length(submatrices)
            submatrices[smi] = "";
        end
    else # surface
        if lorr == LHS
            submatricesll = Array{String, 2}(undef, dofsper, dofsper);
            submatriceslr = Array{String, 2}(undef, dofsper, dofsper);
            submatricesrl = Array{String, 2}(undef, dofsper, dofsper);
            submatricesrr = Array{String, 2}(undef, dofsper, dofsper);
            submatrices = [submatricesll, submatriceslr, submatricesrl, submatricesrr];
        else # RHS
            submatricesl = Array{String, 1}(undef, dofsper);
            submatricesr = Array{String, 1}(undef, dofsper);
            submatrices = [submatricesl, submatricesr];
        end
        for smi=1:length(submatrices)
            for smj=1:length(submatrices[smi])
                submatrices[smi][smj] = "";
            end
        end
    end

    if typeof(var) <: Array # multiple variables
        for vi=1:length(var) # variables
            # Process the terms for this variable
            for ci=1:length(terms[vi]) # components
                for i=1:length(terms[vi][ci])
                    (term_result, test_ind, trial_ind) = generate_term_calculation_dg_julia(terms[vi][ci][i], var, lorr, vors);

                    # println(terms[vi][ci][i])
                    # println(term_result * " : "*string(test_ind)*", "*string(trial_ind))

                    # Find the appropriate submatrix for this term
                    if vors == "volume"
                        submati = offset_ind[vi] + test_ind;
                        submatj = trial_ind;
                        if lorr == LHS
                            submat_ind = submati + dofsper * (submatj-1);
                        else
                            submat_ind = submati;
                        end

                        if length(submatrices[submat_ind]) > 1
                            submatrices[submat_ind] *= " .+ " * term_result;
                        else
                            submatrices[submat_ind] = term_result;
                        end

                    else # surface
                        submati = offset_ind[vi] + test_ind;
                        submatj = trial_ind;
                        if lorr == LHS
                            submat_ind = submati + dofsper * (submatj-1);
                            smparts = 4;
                        else
                            submat_ind = submati;
                            smparts = 2;
                        end

                        for smi=1:smparts
                            if length(submatrices[smi][submat_ind]) > 1
                                if length(term_result[smi]) > 0
                                    submatrices[smi][submat_ind] *= " .+ " * term_result[smi];
                                end
                            else
                                submatrices[smi][submat_ind] = term_result[smi];
                            end
                        end
                    end
                end
            end
        end # vi

    else # only one variable
        # Process the terms for this variable
        for ci=1:length(terms) # components
            for i=1:length(terms[ci])
                (term_result, test_ind, trial_ind) = generate_term_calculation_dg_julia(terms[ci][i], var, lorr, vors);

                # Find the appropriate submatrix for this term
                if vors == "volume"
                    if lorr == LHS
                        submat_ind = test_ind + dofsper * (trial_ind-1);
                    else
                        submat_ind = test_ind;
                    end

                    if length(submatrices[submat_ind]) > 1
                        submatrices[submat_ind] *= " .+ " * term_result;
                    else
                        submatrices[submat_ind] = term_result;
                    end

                else # surface
                    if lorr == LHS
                        submat_ind = test_ind + dofsper * (trial_ind-1);
                        smparts = 4;
                    else
                        submat_ind = test_ind;
                        smparts = 2;
                    end

                    for smi=1:smparts
                        if length(submatrices[smi][submat_ind]) > 1
                            if length(term_result[smi]) > 0
                                submatrices[smi][submat_ind] *= " .+ " * term_result[smi];
                            end
                        else
                            submatrices[smi][submat_ind] = term_result[smi];
                        end
                    end
                end
            end
        end
    end

    # Put the submatrices together into element_matrix or element_vector
    if lorr == LHS
        if vors == "volume"
            for emi=1:dofsper
                for emj=1:dofsper
                    if length(submatrices[emi, emj]) > 1
                        rangei = "(("*string(emi-1)*")*refel.Np + 1):("*string(emi)*"*refel.Np)";
                        rangej = "(("*string(emj-1)*")*refel.Np + 1):("*string(emj)*"*refel.Np)";
                        code *= "element_matrix["*rangei*", "*rangej*"] = " * submatrices[emi,emj] * "\n";
                    end
                end
            end
            code *= "return element_matrix;\n"

        else # surface
            mat_tags = ["LL", "LR", "RL", "RR"];
            frefs = [["1","1"],["1","2"],["2","1"],["2","2"]];
            for emi=1:dofsper
                for emj=1:dofsper
                    for smpart=1:4
                        if length(submatrices[smpart][emi, emj]) > 1
                            rangei = "(("*string(emi-1)*")*Nfp + 1):("*string(emi)*"*Nfp)";
                            rangej = "(("*string(emj-1)*")*Nfp + 1):("*string(emj)*"*Nfp)";
                            code *= "element_matrix"*mat_tags[smpart]*"["*rangei*", "*rangej*"] = (" * submatrices[smpart][emi,emj] * ")[refel.face2local[frefelind["*frefs[smpart][1]*"]], refel.face2local[frefelind["*frefs[smpart][2]*"]]]\n";
                        end
                    end
                end
            end
            code *= "return [element_matrixLL, element_matrixLR, element_matrixRL, element_matrixRR];\n"
        end

    else # RHS
        if vors == "volume"
            for emi=1:dofsper
                if length(submatrices[emi]) > 1
                    rangei = "(("*string(emi-1)*")*refel.Np + 1):("*string(emi)*"*refel.Np)";
                    code *= "element_vector["*rangei*"] = " * submatrices[emi] * "\n";
                end
            end
            code *= "return element_vector;\n"

        else # surface
            mat_tags = ["L", "R"];
            frefs = ["1","2"];
            for emi=1:dofsper
                for smpart=1:2
                    if length(submatrices[smpart][emi]) > 1
                        rangei = "(("*string(emi-1)*")*Nfp + 1):("*string(emi)*"*Nfp)";
                        code *= "element_vector"*mat_tags[smpart]*"["*rangei*"] = (" * submatrices[smpart][emi] * ")[refel.face2local[frefelind["*frefs[smpart]*"]]]\n";
                    end
                end
            end
            code *= "return [element_vectorL, element_vectorR];\n"
        end
    end

    return code;
end

function generate_term_calculation_dg_julia(term, var, lorr, vors)
    result = "";

    if vors == "volume"
        if lorr == LHS
            (test_part, trial_part, coef_part, test_ind, trial_ind) = separate_factors(term, var);
            # LHS: test_part * diagm(weight_part .* coef_part) * trial_part
            if !(coef_part === nothing)
                result = string(replace_entities_with_symbols(test_part)) * " * diagm(wdetj .* " *
                        string(replace_entities_with_symbols(coef_part)) * ") * " *
                        string(replace_entities_with_symbols(trial_part));
            else # no coef_part
                result = string(replace_entities_with_symbols(test_part)) * " * diagm(wdetj) * " *
                        string(replace_entities_with_symbols(trial_part));
            end
        else
            (test_part, trial_part, coef_part, test_ind, trial_ind) = separate_factors(term);
            # RHS: test_part * (weight_part .* coef_part)
            if !(coef_part === nothing)
                result = string(replace_entities_with_symbols(test_part)) * " * (wdetj .* " *
                        string(replace_entities_with_symbols(coef_part)) * ")";
            else
                result = string(replace_entities_with_symbols(test_part)) * " * (wdetj)";
            end
        end

    else # surface
        if lorr == LHS
            result = ["","","",""];
            (test_part, trial_part, coef_part, test_ind, trial_ind) = separate_factors(term, var);

            # Determine which result, LL, LR, RL, RR to add to based on test and trial parts.
            # For example, SIDE1 test -> L and SIDE2 trial -> R will contribute to LR.
            # If there's no SIDE flag, diagonals are possible.
            test_side = 0;
            trial_side = 0;
            if "DGSIDE1" in test_part.flags
                test_side = 1;
            elseif "DGSIDE2" in test_part.flags
                test_side = 2;
            else
                test_part_side1 = copy(test_part);
                test_part_side2 = copy(test_part);
                push!(test_part_side1.flags, "DGSIDE1");
                push!(test_part_side2.flags, "DGSIDE2");
            end
            if "DGSIDE1" in trial_part.flags
                trial_side = 1;
                # also apply side1 to known variables
                coef_part_side1 = add_flag_to_var_entities(copy(coef_part), variables, "DGSIDE1", nevermind="DGSIDE", copy_ent=true);
            elseif "DGSIDE2" in trial_part.flags
                trial_side = 2;
                # also apply side2 to known variables
                coef_part_side2 = add_flag_to_var_entities(copy(coef_part), variables, "DGSIDE2", nevermind="DGSIDE", copy_ent=true);
            else
                trial_part_side1 = copy(trial_part);
                trial_part_side2 = copy(trial_part);
                push!(trial_part_side1.flags, "DGSIDE1");
                push!(trial_part_side2.flags, "DGSIDE2");
                # also apply side to known variables
                coef_part_side1 = add_flag_to_var_entities(copy(coef_part), variables, "DGSIDE1", nevermind="DGSIDE", copy_ent=true);
                coef_part_side2 = add_flag_to_var_entities(copy(coef_part), variables, "DGSIDE2", nevermind="DGSIDE", copy_ent=true);
            end

            # LHS: test_part * diagm(weight_part .* coef_part) * trial_part
            if test_side == 1
                if trial_side == 1
                    # LL only
                    if !(coef_part === nothing)
                        result[1] = string(replace_entities_with_symbols(test_part)) * " * diagm(wdetj .* " *
                                string(replace_entities_with_symbols(coef_part_side1)) * ") * " *
                                string(replace_entities_with_symbols(trial_part));
                    else # no coef_part
                        result[1] = string(replace_entities_with_symbols(test_part)) * " * diagm(wdetj) * " *
                                string(replace_entities_with_symbols(trial_part));
                    end

                elseif trial_side == 2
                    # LR only
                    if !(coef_part === nothing)
                        result[2] = string(replace_entities_with_symbols(test_part)) * " * diagm(wdetj .* " *
                                string(replace_entities_with_symbols(coef_part_side2)) * ") * " *
                                string(replace_entities_with_symbols(trial_part));
                    else # no coef_part
                        result[2] = string(replace_entities_with_symbols(test_part)) * " * diagm(wdetj) * " *
                                string(replace_entities_with_symbols(trial_part));
                    end

                else
                    # LL only and must change trial symbols to side1
                    if !(coef_part === nothing)
                        result[1] = string(replace_entities_with_symbols(test_part)) * " * diagm(wdetj .* " *
                                string(replace_entities_with_symbols(coef_part_side1)) * ") * " *
                                string(replace_entities_with_symbols(trial_part_side1));
                    else # no coef_part
                        result[1] = string(replace_entities_with_symbols(test_part)) * " * diagm(wdetj) * " *
                                string(replace_entities_with_symbols(trial_part_side1));
                    end
                end

            elseif test_side == 2
                if trial_side == 1
                    # RL only
                    if !(coef_part === nothing)
                        result[3] = string(replace_entities_with_symbols(test_part)) * " * diagm(wdetj .* " *
                                string(replace_entities_with_symbols(coef_part_side1)) * ") * " *
                                string(replace_entities_with_symbols(trial_part));
                    else # no coef_part
                        result[3] = string(replace_entities_with_symbols(test_part)) * " * diagm(wdetj) * " *
                                string(replace_entities_with_symbols(trial_part));
                    end

                elseif trial_side == 2
                    # RR only
                    if !(coef_part === nothing)
                        result[4] = string(replace_entities_with_symbols(test_part)) * " * diagm(wdetj .* " *
                                string(replace_entities_with_symbols(coef_part_side2)) * ") * " *
                                string(replace_entities_with_symbols(trial_part));
                    else # no coef_part
                        result[4] = string(replace_entities_with_symbols(test_part)) * " * diagm(wdetj) * " *
                                string(replace_entities_with_symbols(trial_part));
                    end

                else
                    # RR only and must change trial symbols to side2
                    if !(coef_part === nothing)
                        result[4] = string(replace_entities_with_symbols(test_part)) * " * diagm(wdetj .* " *
                                string(replace_entities_with_symbols(coef_part_side2)) * ") * " *
                                string(replace_entities_with_symbols(trial_part_side2));
                    else # no coef_part
                        result[4] = string(replace_entities_with_symbols(test_part)) * " * diagm(wdetj) * " *
                                string(replace_entities_with_symbols(trial_part_side2));
                    end
                end

            else # no test side specified. Diagonals are possible
                if trial_side == 1
                    # LL only
                    if !(coef_part === nothing)
                        result[1] = string(replace_entities_with_symbols(test_part_side1)) * " * diagm(wdetj .* " *
                                string(replace_entities_with_symbols(coef_part_side1)) * ") * " *
                                string(replace_entities_with_symbols(trial_part));
                    else # no coef_part
                        result[1] = string(replace_entities_with_symbols(test_part_side1)) * " * diagm(wdetj) * " *
                                string(replace_entities_with_symbols(trial_part));
                    end

                elseif trial_side == 2
                    # RR only
                    if !(coef_part === nothing)
                        result[4] = string(replace_entities_with_symbols(test_part_side2)) * " * diagm(wdetj .* " *
                                string(replace_entities_with_symbols(coef_part_side2)) * ") * " *
                                string(replace_entities_with_symbols(trial_part));
                    else # no coef_part
                        result[4] = string(replace_entities_with_symbols(test_part_side2)) * " * diagm(wdetj) * " *
                                string(replace_entities_with_symbols(trial_part));
                    end

                else
                    # LL and RR and must change both symbols
                    if !(coef_part === nothing)
                        result[1] = string(replace_entities_with_symbols(test_part_side1)) * " * diagm(wdetj .*" *
                                string(replace_entities_with_symbols(coef_part_side1)) * ") * " *
                                string(replace_entities_with_symbols(trial_part_side1));
                        result[4] = string(replace_entities_with_symbols(test_part_side2)) * " * diagm(wdetj .* " *
                                string(replace_entities_with_symbols(coef_part_side2)) * ") * " *
                                string(replace_entities_with_symbols(trial_part_side2));
                    else # no coef_part
                        result[1] = string(replace_entities_with_symbols(test_part_side1)) * " * diagm(wdetj) * " *
                                string(replace_entities_with_symbols(trial_part_side1));
                        result[4] = string(replace_entities_with_symbols(test_part_side2)) * " * diagm(wdetj) * " *
                                string(replace_entities_with_symbols(trial_part_side2));
                    end
                end
            end

        else # RHS
            result = ["",""];
            (test_part, trial_part, coef_part, test_ind, trial_ind) = separate_factors(term);
            # figure out what side, if any, the coefficient part is on
            function get_coef_side(ex)
                side = 0;
                if typeof(ex) == Expr
                    for i=1:length(ex.args)
                        side = max(side, get_coef_side(ex.args[i]));
                    end
                elseif typeof(ex) == SymEntity
                    if "DGSIDE1" in ex.flags
                        side = 1;
                    elseif "DGSIDE2" in ex.flags
                        side = 2;
                    end
                end
                return side;
            end
            coef_side = get_coef_side(coef_part);
            if coef_side == 0
                coef_part_side1 = add_flag_to_var_entities(copy(coef_part), variables, "DGSIDE1", nevermind="DGSIDE", copy_ent=true);
                coef_part_side2 = add_flag_to_var_entities(copy(coef_part), variables, "DGSIDE2", nevermind="DGSIDE", copy_ent=true);
            end
            test_side = 0;
            if "DGSIDE1" in test_part.flags
                test_side = 1;
            elseif "DGSIDE2" in test_part.flags
                test_side = 2;
            else
                test_part_side1 = copy(test_part);
                test_part_side2 = copy(test_part);
                push!(test_part_side1.flags, "DGSIDE1");
                push!(test_part_side2.flags, "DGSIDE2");
            end

            # RHS: test_part * (weight_part .* coef_part)
            if test_side == 1
                if !(coef_part === nothing)
                    if coef_side == 1
                        result[1] = string(replace_entities_with_symbols(test_part)) * " * (wdetj .* (refel.surf_Q[frefelind[1]] * " *
                            string(replace_entities_with_symbols(coef_part)) * "))";
                    elseif coef_side == 2
                        result[1] = string(replace_entities_with_symbols(test_part)) * " * (wdetj .* (refel.surf_Q[frefelind[2]] * " *
                            string(replace_entities_with_symbols(coef_part)) * "))";
                    else
                        result[1] = string(replace_entities_with_symbols(test_part)) * " * (wdetj .* (refel.surf_Q[frefelind[1]] * " *
                            string(replace_entities_with_symbols(coef_part_side1)) * " .+ refel.surf_Q[frefelind[2]] * " *
                            string(replace_entities_with_symbols(coef_part_side2)) * "))";
                    end

                else
                    result[1] = string(replace_entities_with_symbols(test_part)) * " * (wdetj .* (refel.surf_Q[frefelind[1]] .+ refel.surf_Q[frefelind[2]]))";
                end
            elseif test_side == 2
                if !(coef_part === nothing)
                    if coef_side == 1
                        result[2] = string(replace_entities_with_symbols(test_part)) * " * (wdetj .* (refel.surf_Q[frefelind[1]] * " *
                            string(replace_entities_with_symbols(coef_part)) * "))";
                    elseif coef_side == 2
                        result[2] = string(replace_entities_with_symbols(test_part)) * " * (wdetj .* (refel.surf_Q[frefelind[2]] * " *
                            string(replace_entities_with_symbols(coef_part)) * "))";
                    else
                        result[2] = string(replace_entities_with_symbols(test_part)) * " * (wdetj .* (refel.surf_Q[frefelind[1]] * " *
                            string(replace_entities_with_symbols(coef_part_side1)) * " .+ refel.surf_Q[frefelind[2]] * " *
                            string(replace_entities_with_symbols(coef_part_side2)) * "))";
                    end
                else
                    result[2] = string(replace_entities_with_symbols(test_part)) * " * (wdetj .* (refel.surf_Q[frefelind[1]] .+ refel.surf_Q[frefelind[2]]))";
                end
            else # no test side specified
                if !(coef_part === nothing)
                    if coef_side == 1
                        result[1] = string(replace_entities_with_symbols(test_part_side1)) * " * (wdetj .* (refel.surf_Q[frefelind[1]] * " *
                            string(replace_entities_with_symbols(coef_part)) * "))";
                        result[2] = string(replace_entities_with_symbols(test_part_side2)) * " * (wdetj .* (refel.surf_Q[frefelind[1]] * " *
                            string(replace_entities_with_symbols(coef_part)) * "))";
                    elseif coef_side == 2
                        result[1] = string(replace_entities_with_symbols(test_part_side1)) * " * (wdetj .* (refel.surf_Q[frefelind[2]] * " *
                            string(replace_entities_with_symbols(coef_part)) * "))";
                        result[2] = string(replace_entities_with_symbols(test_part_side2)) * " * (wdetj .* (refel.surf_Q[frefelind[2]] * " *
                            string(replace_entities_with_symbols(coef_part)) * "))";
                    else
                        result[1] = string(replace_entities_with_symbols(test_part_side1)) * " * (wdetj .* (refel.surf_Q[frefelind[1]] * " *
                            string(replace_entities_with_symbols(coef_part_side1)) * "))";
                        result[2] = string(replace_entities_with_symbols(test_part_side2)) * " * (wdetj .* (refel.surf_Q[frefelind[2]] * " *
                            string(replace_entities_with_symbols(coef_part_side2)) * "))";
                    end
                else
                    result[1] = string(replace_entities_with_symbols(test_part_side1)) * " * (wdetj .* (refel.surf_Q[frefelind[1]] .+ refel.surf_Q[frefelind[2]]))";
                    result[2] = string(replace_entities_with_symbols(test_part_side2)) * " * (wdetj .* (refel.surf_Q[frefelind[1]] .+ refel.surf_Q[frefelind[2]]))";
                end
            end

        end
    end

    return (result, test_ind, trial_ind);
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
## A small driver file to run the 3 part tests for CLT and combine
#run like this
#julia -p 4 -L testCLTHarness.jl tests_3partCLT.jl
#pass arguments for things via ARGS[1] is numRun per batch

spath = "./Results/3PartCLT_plot"
N_grid = collect(1:20)
const n = 2^17
numRuns = parse(Int, ARGS[1])
dist_type = ARGS[2]

tic()
a = @spawn test_3PartCLT(spath, numRuns, n, N_grid, 8675309, dist_type)
b = @spawn test_3PartCLT(spath, numRuns, n, N_grid, 5164174290, dist_type)
c = @spawn test_3PartCLT(spath, numRuns, n, N_grid, 5167462266, dist_type)
d = @spawn test_3PartCLT(spath, numRuns, n, N_grid, 123456, dist_type)

######
file_a = fetch(a)
file_b = fetch(b)
file_c = fetch(c)
file_d = fetch(d)

time_stamp = toc()

##read everyone in, throw away a line
data, header = readdlm(file_a, ',', header=true)

data_t = readdlm(file_b, ',', skipstart=1)
data_t[:, 1] += numRuns
data = vcat(data, data_t)

data_t = readdlm(file_c, ',', skipstart=1)
data_t[:, 1] += 2numRuns
data = vcat(data, data_t)

data_t = readdlm(file_d, ',', skipstart=1)
data_t[:, 1] += 3numRuns
data = vcat(data, data_t)

#strip the name of file_a to make the numbers better
println("Filea \t", file_a)
println("spath \t", spath)

indx = search(file_a, spath)[end] + 1

f = open("$(spath)_$(file_a[indx:end])_full_$(4numRuns).csv", "w")
writedlm(f,  header)
writedlm(f,  data)
close(f)

println("Num Paths: \t $(numRuns) \t Time:", time_stamp )


text:
import BenchmarkTools as BT
import Dates
import SimpleDraw as SD

const SHAPE_TYPES = [
                     SD.Point,
                     SD.Background,
                     SD.Line,
                     SD.ThickLine,
                     SD.Circle,
                     SD.FilledCircle,
                     SD.ThickCircle,
                     SD.Rectangle,
                     SD.FilledRectangle,
                     SD.ThickRectangle,
                     SD.Character,
                     SD.TextLine,
                    ]

const SIZES = [64, 256, 1024]

get_shape(::Type{SD.Point}, n) = SD.Point(n ÷ 2 + 1, n ÷ 2 + 1)
get_shape(::Type{SD.Background}, n) = SD.Background()
get_shape(::Type{SD.Line}, n) = SD.Line(SD.Point(n ÷ 8 + 1, 2), SD.Point((7 * n) ÷ 8, n - 1))
get_shape(::Type{SD.ThickLine}, n) = SD.ThickLine(SD.Point(n ÷ 8 + 1, n ÷ 8 + 1), SD.Point((7 * n) ÷ 8, (7 * n) ÷ 8), (n ÷ 8) - 1)
get_shape(::Type{SD.Circle}, n) = SD.Circle(SD.Point(2, 2), n - 2)
get_shape(::Type{SD.FilledCircle}, n) = SD.FilledCircle(SD.Point(2, 2), n - 2)
get_shape(::Type{SD.ThickCircle}, n) = SD.ThickCircle(SD.Point(2, 2), n - 2, n ÷ 4)
get_shape(::Type{SD.Rectangle}, n) = SD.Rectangle(SD.Point(2, 2), n - 1, n - 1)
get_shape(::Type{SD.FilledRectangle}, n) = SD.FilledRectangle(SD.Point(2, 2), n - 1, n - 1)
get_shape(::Type{SD.ThickRectangle}, n) = SD.ThickRectangle(SD.Point(2, 2), n - 1, n - 1, n ÷ 4)
get_shape(::Type{SD.Character}, n) = SD.Character(SD.Point(2, 2), 'A', SD.TERMINUS_32_16)

function get_shape(::Type{SD.TextLine}, n)
    str = ""
    for j in 1 : n ÷ 16
        str = str * Char(rand(65:90))
    end

    return SD.TextLine(SD.Point(1, 1), str, SD.TERMINUS_32_16)
end

function get_benchmarks(shape_types, sizes)
    benchmarks = Dict()

    color = 0x00ffffff
    for n in sizes
        image = zeros(typeof(color), n, n)
        for ShapeType in shape_types
            shape_name = nameof(ShapeType)
            shape = get_shape(ShapeType, n)
            benchmark = BT.@benchmark SD.draw!($(Ref(image))[], $(Ref(shape))[], $(Ref(color))[])
            memory, median_time = get_summary(benchmark)
            benchmarks[(shape_name, n)] = (shape, memory, median_time)
            println("############################")
            println("shape type = $(shape_name)")
            println("image height = $(n)")
            println("image width = $(n)")
            println("median_time = $(median_time)")
            println("memory = $(memory)")
            println("shape = $(Base.limitrepr(shape))")
        end
        fill!(image, zero(color))
    end

    return benchmarks
end

function get_summary(trial)
    median_trial = BT.median(trial)
    memory = BT.prettymemory(median_trial.memory)
    median_time = BT.prettytime(median_trial.time)
    return memory, median_time
end

function get_table(shape_types, sizes, benchmarks)
    table = "|shape type|"
    table = table * "image height|"
    table = table * "image width|"
    table = table * "median time|"
    table = table * "memory|"
    table = table * "shape|"
    table = table * "\n"

    table = table * ":---|"
    table = table * ":---|"
    table = table * ":---|"
    table = table * ":---|"
    table = table * ":---|"
    table = table * ":---|"
    table = table * "\n"

    for ShapeType in shape_types
        for n in sizes
            shape_name = nameof(ShapeType)
            shape, memory, median_time = benchmarks[(shape_name, n)]

            table = table * "|"
            table = table * "$(shape_name)|"
            table = table * "$(n)|"
            table = table * "$(n)|"
            table = table * "$(median_time)|"
            table = table * "$(memory)|"
            table = table * "$(Base.limitrepr(shape))|"
            table = table * "\n"
        end
    end

    return table
end

function generate_benchmark_file(; shape_types = SHAPE_TYPES, sizes = SIZES, file_name = nothing)
    date = Dates.format(Dates.now(), "yyyy_mm_dd_HH_MM_SS")

    if isnothing(file_name)
        file_name = date * ".md"
    end

    io = open(file_name, "w")

    println(io, "Timestamp: $(date) (yyyy_mm_dd_HH_MM_SS)")
    println(io)
    println(io, "Shapes are drawn on an image of type `Matrix{UInt32}` with a color of type `UInt32`")
    println(io)

    benchmarks = get_benchmarks(shape_types, sizes)

    table = get_table(shape_types, sizes, benchmarks)

    println(io, table)

    close(io)

    return nothing
end