text:
include("../roblib.jl")

function mark_centers(r)
    num_steps = through_rectangles_into_angle(r, (Down, Left))

    num_steps_to_right = sum(num_steps[1:2:end])
    num_steps_to_up = sum(num_steps[2:2:end])

    movements!(r, Up, num_steps_to_up)
    putmarker!(r)
    num_steps_to_down = get_num_movements!(r, Up)

    movements!(r, Right, num_steps_to_right)
    putmarker!(r)
    num_steps_to_left = get_num_movements!(r, Right)

    movements!(r, Down, num_steps_to_down)
    putmarker!(r)
    movements!(r, Down)

    movements!(r, Left, num_steps_to_left)
    putmarker!(r)
    movements!(r, Left)

    movements!(r, (Up, Right), reverse(num_steps))
end


text:
# Model detector efficiency (Window +  Detector Crystal)

struct DetectorEfficiency
    name::String
    window::AbstractWindow
    surface::Vector{Film}
    active::Film
end

Base.show(io::IO, de::DetectorEfficiency) = print(io, "$(de.name)[$(de.window)]")

SDDEfficiency(
    window::AbstractWindow;
    thickness = 0.0370,
    deadlayer = 30.0e-7,
    entrance = Film(pure(n"Al"), 10.0e-7),
) = DetectorEfficiency(
    "SDD",
    window,
    [entrance, Film(pure(n"Si"), deadlayer)],
    Film(pure(n"Si"), thickness),
)

SiLiEfficiency(
    window::AbstractWindow;
    thickness = 0.250,
    deadlayer = 30.0e-7,
    entrance = Film(pure(n"Al"), 10.0e-7),
) = DetectorEfficiency(
    "Si(Li)",
    window,
    [entrance, Film(pure(n"Si"), deadlayer)],
    Film(pure(n"Si"), thickness),
)

efficiency(aa::DetectorEfficiency, energy::Float64, angle::Float64 = π / 2) =
    transmission(aa.window, energy, angle) *
    (1.0 - transmission(aa.active, energy, angle)) *
    mapreduce(lyr -> transmission(lyr, energy, angle), *, aa.surface)

"""
    detectorresponse(det::EDSDetector, eff::DetectorEfficiency, incidence::Float64=π/2)::AbstractMatrix

Build a matrix which models the detector response including aspects like the detector efficiency, the resolution, the
escape peaks.  All the warts that can be modeled within a linear model but not things like coincidence peaks that
are non-linear.  This function can (!should!) be specialized for more sophisticated detector models that include more
warts.

Example:

    genint = computegeneratedintensity(....) # Either characteristic or Bremsstrahlung...
    det = simpleEDS(4096, 5.0, 0.0, 132.0, 10)
    eff = SDDEfficiency(AP33Tabulation(); thickness=0.0370, deadlayer=30.0e-7, entrance=Film(pure(n"Al"), 10.0e-7))
    resp = detectorresponse(det, eff)
    # finally compute the measured signal
    measured = genint*resp
"""
function detectorresponse(
    det::EDSDetector,
    eff::DetectorEfficiency,
    incidence::Float64 = π / 2,
)::AbstractMatrix
    res = zeros(Float64, (channelcount(det), channelcount(det)))
    # An x-ray with energy in ch will be dispersed among a range of channels about ch
    for ch = lld(det):channelcount(det) # X-ray energy by channel
        el, eh = energy(ch, det), energy(ch + 1, det)  # X-ray energy
        effic, fwhm =
            0.5 * (efficiency(eff, eh, incidence) + efficiency(eff, el, incidence)),
            resolution(0.5 * (eh + el), det)
        roc =
            max(
                lld(det),
                channel(el - 3.0 * fwhm, det),
            ):min(channelcount(det), channel(el + 3.0 * fwhm, det))
        res[ch, roc] = map(ch2 -> effic * profile(ch2, 0.5 * (eh + el), det), roc)
    end
    return res
end

ext: jl


text:
"""
    struct SBMLModel

Thin wrapper around the model from SBML.jl library. Allows easy conversion from
SBML to any other model format.
"""
struct SBMLModel <: MetabolicModel
    sbml::SBML.Model
end

"""
    reactions(model::SBMLModel)::Vector{String}

Get reactions from a [`SBMLModel`](@ref).
"""
reactions(model::SBMLModel)::Vector{String} = [k for k in keys(model.sbml.reactions)]

"""
    metabolites(model::SBMLModel)::Vector{String}

Get metabolites from a [`SBMLModel`](@ref).
"""
metabolites(model::SBMLModel)::Vector{String} = [k for k in keys(model.sbml.species)]

"""
    n_reactions(model::SBMLModel)::Int

Efficient counting of reactions in [`SBMLModel`](@ref).
"""
n_reactions(model::SBMLModel)::Int = length(model.sbml.reactions)

"""
    n_metabolites(model::SBMLModel)::Int

Efficient counting of metabolites in [`SBMLModel`](@ref).
"""
n_metabolites(model::SBMLModel)::Int = length(model.sbml.species)

"""
    stoichiometry(model::SBMLModel)::SparseMat

Recreate the stoichiometry matrix from the [`SBMLModel`](@ref).
"""
function stoichiometry(model::SBMLModel)::SparseMat
    _, _, S = SBML.getS(model.sbml)
    return S
end

"""
    bounds(model::SBMLModel)::Tuple{SparseVec,SparseVec}

Get the lower and upper flux bounds of model [`SBMLModel`](@ref). Throws `DomainError` in
case if the SBML contains mismatching units.
"""
function bounds(model::SBMLModel)::Tuple{SparseVec,SparseVec}
    lbu = SBML.getLBs(model.sbml)
    ubu = SBML.getUBs(model.sbml)

    unit = lbu[1][2]
    getvalue = (val, _)::Tuple -> val
    getunit = (_, unit)::Tuple -> unit

    allunits = unique([getunit.(lbu) getunit.(ubu)])
    length(allunits) == 1 || throw(
        DomainError(
            allunits,
            "The SBML file uses multiple units; loading needs conversion",
        ),
    )

    return sparse.((getvalue.(lbu), getvalue.(ubu)))
end

"""
    balance(model::SBMLModel)::SparseVec

Balance vector of a [`SBMLModel`](@ref). This is always zero.
"""
balance(model::SBMLModel)::SparseVec = spzeros(n_metabolites(model))

"""
    objective(model::SBMLModel)::SparseVec

Objective of the [`SBMLModel`](@ref).
"""
objective(model::SBMLModel)::SparseVec = SBML.getOCs(model.sbml)

"""
    genes(model::SBMLModel)::Vector{String}

Get genes of a [`SBMLModel`](@ref).
"""
genes(model::SBMLModel)::Vector{String} = [k for k in keys(model.sbml.gene_products)]

"""
    n_genes(model::SBMLModel)::Int

Get number of genes in [`SBMLModel`](@ref).
"""
n_genes(model::SBMLModel)::Int = length(model.sbml.gene_products)

"""
    reaction_gene_association(model::SBMLModel, rid::String)::Maybe{GeneAssociation}

Retrieve the [`GeneAssociation`](@ref) from [`SBMLModel`](@ref).
"""
reaction_gene_association(model::SBMLModel, rid::String)::Maybe{GeneAssociation} =
    _maybemap(_parse_grr, model.sbml.reactions[rid].gene_product_association)

"""
    metabolite_formula(model::SBMLModel, mid::String)::Maybe{MetaboliteFormula}

Get [`MetaboliteFormula`](@ref) from a chosen metabolite from [`SBMLModel`](@ref).
"""
metabolite_formula(model::SBMLModel, mid::String)::Maybe{MetaboliteFormula} =
    _maybemap(_parse_formula, model.sbml.species[mid].formula)

"""
    metabolite_charge(model::SBMLModel, mid::String)::Maybe{Int}

Get charge of a chosen metabolite from [`SBMLModel`](@ref).
"""
metabolite_charge(model::SBMLModel, mid::String)::Maybe{Int} =
    model.sbml.species[mid].charge

function _sbml_export_annotation(annotation)::Maybe{String}
    if isnothing(annotation) || isempty(annotation)
        nothing
    elseif length(annotation) != 1 || first(annotation).first != ""
        @_io_log @warn "Possible data loss: annotation cannot be exported to SBML" annotation
        join(["$k: $v" for (k, v) in annotation], "\n")
    else
        first(annotation).second
    end
end

_sbml_export_notes = _sbml_export_annotation

"""
    reaction_stoichiometry(model::SBMLModel, rxn_id::String)::Dict{String, Float64}

Return the reaction equation of reaction with id `rxn_id` in model. The reaction
equation maps metabolite ids to their stoichiometric coefficients.
"""
function reaction_stoichiometry(m::SBMLModel, rxn_id::String)::Dict{String, Float64}
    m.sbml.reactions[rxn_id].stoichiometry
end

"""
    Base.convert(::Type{SBMLModel}, mm::MetabolicModel)

Convert any metabolic model to [`SBMLModel`](@ref).
"""
function Base.convert(::Type{SBMLModel}, mm::MetabolicModel)
    if typeof(mm) == SBMLModel
        return mm
    end

    mets = metabolites(mm)
    rxns = reactions(mm)
    stoi = stoichiometry(mm)
    (lbs, ubs) = bounds(mm)
    ocs = objective(mm)
    comps = _default.("", metabolite_compartment.(Ref(mm), mets))
    compss = Set(comps)

    return SBMLModel(
        SBML.Model(
            Dict(), # parameters
            Dict("" => []), # units
            Dict([
                comp =>
                    SBML.Compartment(nothing, nothing, nothing, nothing, nothing, nothing)
                for comp in compss
            ],),
            Dict([
                mid => SBML.Species(
                    nothing, # name
                    _default("", comps[mi]), # compartment
                    nothing, # no information about boundary conditions
                    metabolite_formula(mm, mid),
                    metabolite_charge(mm, mid),
                    nothing, # initial amount
                    nothing, # only substance unit flags
                    _sbml_export_notes(metabolite_notes(mm, mid)),
                    _sbml_export_annotation(metabolite_annotations(mm, mid)),
                ) for (mi, mid) in enumerate(mets)
            ],),
            Dict([
                rid => SBML.Reaction(
                    Dict([
                        mets[i] => stoi[i, ri] for
                        i in SparseArrays.nonzeroinds(stoi[:, ri])
                    ],),
                    (lbs[ri], ""),
                    (ubs[ri], ""),
                    ocs[ri],
                    _maybemap(
                        x -> _unparse_grr(SBML.GeneProductAssociation, x),
                        reaction_gene_association(mm, rid),
                    ),
                    nothing,
                    _sbml_export_notes(reaction_notes(mm, rid)),
                    _sbml_export_annotation(reaction_annotations(mm, rid)),
                ) for (ri, rid) in enumerate(rxns)
            ],),
            Dict([
                gid => SBML.GeneProduct(
                    nothing,
                    nothing,
                    _sbml_export_notes(gene_notes(mm, gid)),
                    _sbml_export_annotation(gene_annotations(mm, gid)),
                ) for gid in genes(mm)
            ],),
            Dict(), # function definitions
        ),
    )
end


text:
function characterize_zeta2(particle_type::String,
							R::Array{Float64, 2},
							Lx::Float64,
							Ly::Float64,
							Lz::Float64,
							X::Array{Float64, 1},
							Y::Array{Float64, 1},
							Z::Array{Float64, 1},
							Q0::Array{Float64, 1},
							Q1::Array{Float64, 1},
							Q2::Array{Float64, 1},
							Q3::Array{Float64, 1},
							A11::Array{Float64, 1},
							A12::Array{Float64, 1},
							A13::Array{Float64, 1},
							A21::Array{Float64, 1},
							A22::Array{Float64, 1},
							A23::Array{Float64, 1},
							A31::Array{Float64, 1},
							A32::Array{Float64, 1},
							A33::Array{Float64, 1},
							number_of_samples::Int64,
							cell_lists::Array{Array{Int64, 1}, 3})

	# Number of cells.
	(number_of_cells_x::Int64, number_of_cells_y::Int64, number_of_cells_z::Int64) = size(cell_lists)

	current_particle_in_cell::Int64 = 0
	is_void::Bool = true

	#x0::Float64 = 0.0
	#y0::Float64 = 0.0
	#z0::Float64 = 0.0

	number_of_points::Int64 = 7
	x::Array{Float64, 1} = zeros(number_of_points)
	y::Array{Float64, 1} = zeros(number_of_points)
	z::Array{Float64, 1} = zeros(number_of_points)

	r::Float64 = 0.0
	s::Float64 = 0.0
	theta::Float64 = 0.0
	mu::Float64 = 1.0

	S2r::Float64 = 0.0
	S2s::Float64 = 0.0
	S3::Float64 = 0.0
	indicator::Array{Bool, 1} = falses(number_of_points)

	zeta2_integral::Array{Float64, 1} = zeros(2)
	integrand::Float64 = 0.0

	vxr::Float64 = 0.0
	vyr::Float64 = 0.0
	vzr::Float64 = 0.0
	vxs::Float64 = 0.0
	vys::Float64 = 0.0
	vzs::Float64 = 0.0
	magnitude::Float64 = 0.0

	current_cell_x::Int64 = 0
	current_cell_y::Int64 = 0
	current_cell_z::Int64 = 0
	number_of_particles_current_cell::Int64 = 0

	for current_sample = 1:number_of_samples
		#if mod(current_sample, 1000) == 0
		#	println(current_sample)
		#end

		# Pick random intial sample position for S2r.
		x[1] = Lx * rand()
		y[1] = Ly * rand()
		z[1] = Lz * rand()

		# Pick random intial sample position for S2s.
		x[3] = Lx * rand()
		y[3] = Ly * rand()
		z[3] = Lz * rand()

		# Pick random intial sample position for S3.
		x[5] = Lx * rand()
		y[5] = Ly * rand()
		z[5] = Lz * rand()

		# Pick random radii.
		r = mu * randexp()
		s = mu * randexp()

		# Pick random direction and position for r in S2r.
		vxr = randn()
		vyr = randn()
		vzr = randn()
		magnitude = sqrt(vxr * vxr + vyr * vyr + vzr * vzr)
		vxr /= magnitude
		vyr /= magnitude
		vzr /= magnitude
		x[2] = x[1] + r * vxr
		y[2] = y[1] + r * vyr
		z[2] = z[1] + r * vzr

		# Pick random direction and position for s in S2s.
		vxs = randn()
		vys = randn()
		vzs = randn()
		magnitude = sqrt(vxs * vxs + vys * vys + vzs * vzs)
		vxs /= magnitude
		vys /= magnitude
		vzs /= magnitude
		x[4] = x[3] + s * vxs
		y[4] = y[3] + s * vys
		z[4] = z[3] + s * vzs

		# Pick random 'r' sample position in random direction for S3.
		vxr = randn()
		vyr = randn()
		vzr = randn()
		magnitude = sqrt(vxr * vxr + vyr * vyr + vzr * vzr)
		vxr /= magnitude
		vyr /= magnitude
		vzr /= magnitude
		x[6] = x[5] + r * vxr
		y[6] = y[5] + r * vyr
		z[6] = z[5] + r * vzr

		# Pick random 's' sample position in random direction for S3.
		vxs = randn()
		vys = randn()
		vzs = randn()
		magnitude = sqrt(vxs * vxs + vys * vys + vzs * vzs)
		vxs /= magnitude
		vys /= magnitude
		vzs /= magnitude
		x[7] = x[5] + s * vxs
		y[7] = y[5] + s * vys
		z[7] = z[5] + s * vzs

		# Compute angle between the two vectors.
		theta = acos(vxs * vxr + vys * vyr + vzs * vzr)

		# Convert to absolute coordinates (w.r.t. to Lx, Ly, Lz)
		for current_point = 1:number_of_points
			x[current_point] = position_mod(x[current_point], Lx)
			y[current_point] = position_mod(y[current_point], Ly)
			z[current_point] = position_mod(z[current_point], Lz)
		end

		# Check whether void or solid for each sample point.
		for current_point = 1:number_of_points

			current_cell_x = convert(Int64, ceil(x[current_point] / Lx * convert(Float64, number_of_cells_x)))
			current_cell_y = convert(Int64, ceil(y[current_point] / Ly * convert(Float64, number_of_cells_y)))
			current_cell_z = convert(Int64, ceil(z[current_point] / Lz * convert(Float64, number_of_cells_z)))
			number_of_particles_current_cell = length(cell_lists[current_cell_x, current_cell_y, current_cell_z])
			current_cell_list = cell_lists[current_cell_x, current_cell_y, current_cell_z]

			is_void = true
			current_particle_in_cell = 0
			while current_particle_in_cell < number_of_particles_current_cell && is_void
				current_particle_in_cell += 1

				vx = signed_distance_mod(x[current_point], X[current_cell_list[current_particle_in_cell]], Lx)
				vy = signed_distance_mod(y[current_point], Y[current_cell_list[current_particle_in_cell]], Ly)
				vz = signed_distance_mod(z[current_point], Z[current_cell_list[current_particle_in_cell]], Lz)

				if particle_type == "sphere"
					if vx^2 + vy^2 + vz^2 <= R[current_cell_list[current_particle_in_cell], 1]^2
						is_void = false
					end
				elseif particle_type == "ellipse"
					# Not supported for S2.
				elseif particle_type == "ellipsoid"
					if vx * (A11[current_cell_list[current_particle_in_cell]] * vx + A12[current_cell_list[current_particle_in_cell]] * vy + A13[current_cell_list[current_particle_in_cell]] * vz) + vy * (A21[current_cell_list[current_particle_in_cell]] * vx + A22[current_cell_list[current_particle_in_cell]] * vy + A23[current_cell_list[current_particle_in_cell]] * vz) + vz * (A31[current_cell_list[current_particle_in_cell]] * vx + A32[current_cell_list[current_particle_in_cell]] * vy + A33[current_cell_list[current_particle_in_cell]] * vz) <= 1.0
						is_void = false
					end
				elseif particle_type == "cuboid"
					(vx, vy, vz) = (A11[current_cell_list[current_particle_in_cell]] * vx + A12[current_cell_list[current_particle_in_cell]] * vy + A13[current_cell_list[current_particle_in_cell]] * vz,
									A21[current_cell_list[current_particle_in_cell]] * vx + A22[current_cell_list[current_particle_in_cell]] * vy + A23[current_cell_list[current_particle_in_cell]] * vz,
									A31[current_cell_list[current_particle_in_cell]] * vx + A32[current_cell_list[current_particle_in_cell]] * vy + A33[current_cell_list[current_particle_in_cell]] * vz)

					if abs(vx) <= R[current_cell_list[current_particle_in_cell], 1] && abs(vy) <= R[current_cell_list[current_particle_in_cell], 2] && abs(vz) <= R[current_cell_list[current_particle_in_cell], 3]
						is_void = false
					end
				elseif particle_type == "superellipsoid"
					(vx, vy, vz) = (A11[current_cell_list[current_particle_in_cell]] * vx + A12[current_cell_list[current_particle_in_cell]] * vy + A13[current_cell_list[current_particle_in_cell]] * vz,
									A21[current_cell_list[current_particle_in_cell]] * vx + A22[current_cell_list[current_particle_in_cell]] * vy + A23[current_cell_list[current_particle_in_cell]] * vz,
									A31[current_cell_list[current_particle_in_cell]] * vx + A32[current_cell_list[current_particle_in_cell]] * vy + A33[current_cell_list[current_particle_in_cell]] * vz)

					if (abs(vx)/R[current_cell_list[current_particle_in_cell], 1])^R[current_cell_list[current_particle_in_cell], 4] +
						(abs(vy)/R[current_cell_list[current_particle_in_cell], 2])^R[current_cell_list[current_particle_in_cell], 4] +
						(abs(vz)/R[current_cell_list[current_particle_in_cell], 3])^R[current_cell_list[current_particle_in_cell], 4] <= 1.0
						is_void = false
					end
				end
			end

			if is_void
				indicator[current_point] = true
			else
				indicator[current_point] = false
			end
		end

		# Compute S2 & S3.
		S2r = 0.0
		S2s = 0.0
		S3 = 0.0

		#integrand = 0.0

		if indicator[1] & indicator[2]
			S2r = 1.0
		end
		if indicator[3] & indicator[4]
			S2s = 1.0
		end
		if indicator[5] & indicator[6] & indicator[7]
			S3 = 1.0
		end

		integrand1 = 0.5 * (3.0 * cos(theta)^2 - 1.0) * S3
		integrand2 = 0.5 * (3.0 * cos(theta)^2 - 1.0) * (- S2r * S2s)

		integrand1 /= 1.0 / mu^2 * exp(-r/mu) * exp(-s/mu) * 0.5 # This is the 'weighting' compensation in importance sampling.
		integrand2 /= 1.0 / mu^2 * exp(-r/mu) * exp(-s/mu) * 0.5 # This is the 'weighting' compensation in importance sampling.

		zeta2_integral[1] += integrand1
		zeta2_integral[2] += integrand2

	end

	return zeta2_integral
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module LibGit2OnlineTests

using Test
import LibGit2
using Random

function transfer_progress(progress::Ptr{LibGit2.TransferProgress}, payload::Dict)
    status = payload[:transfer_progress]
    progress = unsafe_load(progress)

    status[] = (current=progress.received_objects, total=progress.total_objects)

    return Cint(0)
end

#########
# TESTS #
#########
# init & clone
mktempdir() do dir
    repo_url = "https://github.com/JuliaLang/Example.jl"

    @testset "Cloning repository" begin
        @testset "HTTPS protocol" begin
            repo_path = joinpath(dir, "Example.HTTPS")
            c = LibGit2.CredentialPayload(allow_prompt=false, allow_git_helpers=false)
            repo = LibGit2.clone(repo_url, repo_path, credentials=c)
            try
                @test isdir(repo_path)
                @test isdir(joinpath(repo_path, ".git"))
            finally
                close(repo)
            end
        end

        @testset "Transfer progress callbacks" begin
            status = Ref((current=0, total=-1))
            callbacks = LibGit2.Callbacks(
                :transfer_progress => (
                    @cfunction(transfer_progress, Cint, (Ptr{LibGit2.TransferProgress}, Any)),
                    status,
                )
            )

            repo_path = joinpath(dir, "Example.TransferProgress")
            c = LibGit2.CredentialPayload(allow_prompt=false, allow_git_helpers=false)
            repo = LibGit2.clone(repo_url, repo_path, credentials=c, callbacks=callbacks)
            try
                @test isdir(repo_path)
                @test isdir(joinpath(repo_path, ".git"))

                @test status[].total >= 0
                @test status[].current == status[].total
            finally
                close(repo)
            end
        end

        @testset "Incorrect URL" begin
            repo_path = joinpath(dir, "Example.IncorrectURL")
            # credentials are required because github tries to authenticate on unknown repo
            cred = LibGit2.UserPasswordCredential("JeffBezanson", "hunter2") # make sure Jeff is using a good password :)
            c = LibGit2.CredentialPayload(cred, allow_prompt=false, allow_git_helpers=false)
            try
                LibGit2.clone(repo_url*randstring(10), repo_path, credentials=c)
                error("unexpected")
            catch ex
                @test isa(ex, LibGit2.Error.GitError)
                @test ex.code == LibGit2.Error.EAUTH
            end
        end

        @testset "Empty Credentials" begin
            repo_path = joinpath(dir, "Example.EmptyCredentials")
            # credentials are required because github tries to authenticate on unknown repo
            cred = LibGit2.UserPasswordCredential("","") # empty credentials cause authentication error
            c = LibGit2.CredentialPayload(cred, allow_prompt=false, allow_git_helpers=false)
            try
                LibGit2.clone(repo_url*randstring(10), repo_path, credentials=c)
                error("unexpected")
            catch ex
                @test isa(ex, LibGit2.Error.GitError)
                @test ex.code == LibGit2.Error.EAUTH
            end
        end
    end
end

end # module

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
# This is a simple `@debug` macro that we can use in the code
# without it slowing the code down, unlike `Base.@debug`.

const DEBUG_LEVEL = Ref(0)

function setdebug!(level::Int)
    DEBUG_LEVEL[] = level
    return nothing
end

"""
    withdebug(level::Int) do
        func()
    end
"""
function withdebug(f, level)
    lvl = DEBUG_LEVEL[]
    try
        setdebug!(level)
        f()
    finally
        setdebug!(lvl)
    end
end

"""
    @debug 1 "msg"
"""
macro debug(level, msg)
    esc(quote
        if DEBUG_LEVEL[] >= $level
            println(string("DEBUG: ", $(QuoteNode(__source__.file)), ":", $(QuoteNode(__source__.line)), " ", $msg))
        end
    end)
end


text:
using irtfun
using RCall, DataFrames

# Single group
Science = R"""
library(mirt)
Science
"""
Science = convert(Matrix, Science)
R"""
coef(mirt(Science, 1))
"""
@profview estimate(Science, 1, graded(), E = Estep(quadpts = 21))
@time estimate(Science, 1, graded(), E = Estep(quadpts = 61))
est_gh = estimate(Science, 1, graded(), E = Estep(quadpts = 5, dentype = :GH))
# mirt
est_nd = estimate(Science, 1, graded(), E = Estep(quadpts = 61))
est_nd = estimate(Science, 2, graded(), E = Estep(quadpts = 5, dentype = :GH))

# Multple group
using irtfun
using RCall, DataFrames
mgdata = R"""
library(mirt)
set.seed(12345)
a <- matrix(abs(rnorm(15,1,.3)), ncol=1)
d <- matrix(rnorm(15,0,.7),ncol=1)
itemtype <- rep('2PL', nrow(a))
N <- 1000
dataset1 <- simdata(a, d, N, itemtype)
dataset2 <- simdata(a, d, N, itemtype, mu = .1, sigma = matrix(1.5))
dat <- rbind(dataset1, dataset2)
group <- c(rep('D1', N), rep('D2', N))
list(dat, group)
"""
R"""
models <- 'F1 = 1-15'
mod_fullconstrain <- multipleGroup(dat, models, group = group,
invariance=c('slopes', 'intercepts', 'free_mean', 'free_var'))
coef(mod_fullconstrain)
"""
data = DataFrame(convert(Matrix, mgdata[1]))
group = convert(Vector, mgdata[2])

@time estimate(data, 1, graded(); E = Estep(quadpts = 21), M = Mstep(N = 5))
@time estimate(data, 1, graded(); group = group, E = Estep(quadpts = 61, dentype = :DN))
@time estimate(data, 1, LPE(); group = group, E = Estep(quadpts = 61, dentype = :DN))
@time estimate(data, 1, graded(); group = group, E = Estep(quadpts = 21, dentype = :GH))

# Guessing model
using RCall, DataFrames
gdata = R"""
library(mirt)
dat <- expand.table(LSAT7)
"""
R"""
coef(mirt(dat, 1, '3PL'))
"""
data = convert(DataFrame, gdata)
@time estimate(data, 1, guessing(); E = Estep(quadpts = 61), M = Mstep(N = 5))


text:
module MyPhysicalConstants

# exporting constants
export c0,
       sigma,
       h,
       ħ,
       kb,
       vf_au,
       vf_ag,
       electron,
       epsilon0,
       m0

# Physical constants

" Speed of light in vacuum in m/s "
const c0    = 299792458.0

" Stefan-Boltzmann's constant in W m^{-2} K^{-4} "
const sigma = 5.670367e-8

" Planck's constant in J s "
const h = 6.6260704e-34
const ħ = 1.0545718e-34

" Boltzmann constant in m^2 kg s^{-2} K^{-1} "
const kb = 1.38064852e-23

" Fermi velocity for gold in m s^{-1} "
const vf_au = 1.4e6

" Fermi velocity for silver in m s^{-1} "
const vf_ag = 1.36e6

" Elementary charge in C "
const electron = 1.60217662e-19

" Permittivity of vacumm in m^{-3} kg^{-1} s^4 A^2"
const epsilon0  = 8.85418782e-12

" Free electron mass in kg"
const m0  = 9.10938356e-31

end # module


text:
"""
    STRF_response(STRF, stimulus; kwargs...)

Computes the "raw" response rate of a neuron with spatio-temporal field given by
`STRF` in response to `stimulus`. Returns a pair `(r,tau)` where `r` is a vector
of rate values and `tau` is the temporal resolution (i.e. the width in seconds
of each time bin represented by `r`). For now, I assume the STRF and the
stimulus use the same frame rate, meaning `tau = frame_time(STRF)`. I'm leaving
the door open to needing to fiddle with this in the future, though (i.e. if
we're comparing an STRF computed at 30Hz to a stimulus recorded at 60Hz, that
sort of thing.)

The "raw" response rate is given by convolving `STRF` and `stimulus` as
functions of time.

"""
function STRF_response(STRF::AbstractStimulus, stimulus::AbstractStimulus;
    flip_STRF_time=false, flip_stimulus_time=false,
    kwargs...)

    M_strf = matrix_form(STRF); M_stim = matrix_form(stimulus)
    if flip_STRF_time
        M_strf = flipdim(M_strf, 2)
    end
    if flip_stimulus_time
        M_stim = flipdim(M_stim, 2)
    end
    dx = prod(STRF.d .* STRF.mm_per_px)
    # dx = prod(patch_size(STRF))
    N_bits_per_frame, N_frames = size(M_stim)
    r = zeros(N_frames)
    for bit = 1:N_bits_per_frame
        r[:] += conv(M_strf[bit,:], M_stim[bit,:])[1:length(r)]
    end
    return r * dx * frame_time(STRF), frame_time(STRF)
end

"""
    scale_response(r, n, phi, Q; objective=:min, kwargs...)

Finds parameters for the function `phi` to minimize penalty function `Q`.
Returns the scaled output, the optimal parameters, and the value of `Q`.

Mathematically, I'm thinking of ``ϕ(x; θ)`` as the transfer function which I'm
applying to `r`, and `Q` as the penalty function for how close `ϕ(r(t); θ)` is
to `n(t)`.

Programmatically, this means:
 * `r::Vector{Float64}` is probably a rate computed from `STRF_response`
 * `n::Vector{Float64}` is probably a single neuron's spike histogram
 * `phi(x, theta)` should accept vectors `x` and `theta` and return a vector the same length as `x`
 * `Q(L,n)` should accept two vectors and return a single number

Some keyword arguments:
 * `objective=:min` whether to maximize or minimize `Q`
 * `d` is the dimensionality of the search space (if using brute force search)
 * `theta0` is an initial guess for the value of the parameters (if using some other algorithm)
 * `algorithm=:brute_force_search` how to search for the optimum

Currently, only supports brute force search, but I'm leaving the door open to
expand this if it seems necessary. (If we're being honest this function is
mostly unnecessary anyway since all it's doing is acting as a wrapper to...
another function. But it'll make my scripts readable so whatever.)

See the `brute_force_optimizer` documentation for some necessary kwargs.

"""
function scale_response(r::Vector{Float64}, n::Vector{Float64}, phi, Q;
    algorithm=:brute_force_search, kwargs...)

    # dkwargs = Dict(kwargs)
    if algorithm != :brute_force_search
        #MAYBEDO expand to other options
        error("scale_response: No, seriously, the only algorithm option is :brute_force_search (got $algorithm)")
    end
    # d = pop!(dkwargs, :d, 3)
    obj_fun(x::Vector) = Q(phi(r, x), n)
    Q_opt, theta_opt, criteria = brute_force_optimizer3(obj_fun, 3; kwargs...)
    return phi(r, theta_opt), theta_opt, Q_opt
end


text:
charlie_mapping = Dict{String,Int}(
    "Nooit" => 1,
    "Bijna nooit" => 2,
    "Soms" => 3,
    "Vaak" => 4,
    "Altijd" => 5
)

charlie2int(s::AbstractString)::Int = charlie_mapping[s]

"""
    charlie2int(q::Int, s::AbstractString)

Convert `s` to Int for question number `q`.
"""
function charlie2int(q::Int, s::AbstractString)::Int
    value = charlie2int(s)
    reversed_questions = [2, 4, 6]
    if q in reversed_questions
        value = 6 - value
    end
    value
end

"""
    resilience2scores(df::DataFrame)

Return the score for questionnaire charlie.
"""
function resilience2scores(df::DataFrame)
    df = dropmissing(df)
    nq = 6
    T = ["v$i" => x -> charlie2int.(i, x) for i in 1:nq]
    transform!(df, :, T...; renamecols=false)
    cols = ["v$i" for i in 1:nq]
    transform!(df, :, cols => ByRow(+) => :score)
    select!(df, :id, :completed_at, :score)
end


text:
# #############################################################################
# #############################################################################
#Interactive PSS plot
function plotpss_interactive(resultsdir::AbstractString;kwargs...)

    # #########################################################################
    #Initiate plot

    #Prepare plot
    fig = plt[:figure](figsize=(9,9))
    axes = fig[:gca](aspect="equal")

    #Plot the non-interactive version
    #and save PSS plot state
    pssarray::Array{Float64,2},residarray::Array{Int64,2} = plotpss(resultsdir;kwargs...)

    # #########################################################################
    #Interactions tasks

    #Highlight bouncepoints on PSS
    #This makes the PSS orbit from a single set of initial conditions stand out
    #from other points, and any periodicity in the orbit to be more obvious
    function pss_highlight(xdata::Float64,ydata::Float64)
        #Get appropriate entry in PSS data arrays based on position on PSS
        range::Array{Float64,1} = [plt[:axis]()...]
        thetaindex::Int64 = getbinindex(xdata,range[1],range[2],size(pssarray,2))
        chiindex::Int64 = getbinindex(ydata,range[4],range[3],size(pssarray,1))

        #Check that entry index is within bounds
        if thetaindex != 0 && chiindex != 0
            #Read info from data arrays and recorded results
            score::Float64 = pssarray[chiindex,thetaindex]
            resultid::Int64 = residarray[chiindex,thetaindex]
            if resultid > 0
                bouncepoints::Array{Float64,2} = readresults(resultsdir,resultid,"bouncepoints")

                #Plot bouncepoints on PSS as red squares with thick black borders for visibility
                plt[:plot](bouncepoints[:,1],bouncepoints[:,2],
                           ".",mec="k",mfc="w",mew=2,ms=10,axes=axes)
                plt[:axis](range)
                plt[:draw]()
            end
        end
    end

    #Dehighlight selected bouncepoints
    function pss_dehighlight()
        #Get lines in axes and remove them all
        lines = axes[:get_lines]()
        for line in lines
            line[:remove]()
        end
        plt[:draw]()
    end

    #Clear PSS plot and replot
    function pss_replot(newrange::Array{Float64,1})
        axes[:cla]()
        pssarray,residarray = plotpss(resultsdir::AbstractString;axes=axes,range=newrange)
        plt[:axis](newrange)
        plt[:draw]()
    end

    #Pan
    function pss_pan(dxdata::Float64,dydata::Float64)
        #Get transformed range
        range::Array{Float64,1} = [plt[:axis]()...]
        newrange::Array{Float64,1} = [range[1]-dxdata, range[2]-dxdata,
                                      range[3]-dydata, range[4]-dydata]
        pss_replot(newrange)
    end

    #Zoom
    function pss_zoom(zoomsteps::Int64)
        #Calculate scaling (8 steps to scale by factor of 2)
        scale::Float64 = 2.0^(zoomsteps/8)

        #Get transformed range
        range::Array{Float64,1} = [plt[:axis]()...]
        centerx::Float64,centery::Float64 = (range[2]+range[1])/2, (range[4]+range[3])/2
        newdx::Float64,newdy::Float64 = scale*(range[2]-range[1])/2, scale*(range[4]-range[3])/2
        newrange::Array{Float64,1} = [centerx-newdx, centerx+newdx,
                                      centery-newdy, centery+newdy]
        pss_replot(newrange)
    end

    # #########################################################################
    #Event handlers

    #Left clicking selects the mode to be emphasized
    #Right clicking deselects any emphasized mode
    #Middle clicking and dragging moves PSS
    #Scrolling zooms in and out

    #Mouse press-drag-release state
    pss_dragstart::Array{Float64,1} = [0.0,0.0]
    pss_button::Int64 = 1
    pss_pressed::Bool = false

    #Press event
    function onpress(mev::PyObject)
        #Ignore all mouse presses outside axes, or when a button is already pressed
        if !pss_pressed && mev[:inaxes] != nothing
            #Store new press-drag-release state
            pss_pressed = true
            pss_button = mev[:button]
            pss_dragstart = Float64[mev[:xdata],mev[:ydata]]
        end
    end

    #Release event
    function onrelease(mev::PyObject)
        #All button releases other than from the one first pressed are ignored
        if pss_pressed && mev[:button] == pss_button
            #Toggle pressed state
            pss_pressed = false

            #Ensure release is in axes
            axes = mev[:inaxes]
            if axes != nothing
                #Choose task depending on button pressed
                if pss_button == 1
                    #Left release
                    #Only clicks without drags are handled
                    if mev[:xdata] == pss_dragstart[1] && mev[:ydata] == pss_dragstart[2]
                        #Select mode with point at click location
                        pss_dehighlight()
                        pss_highlight(mev[:xdata],mev[:ydata])
                    end
                elseif pss_button == 3
                    #Right release
                    #Only clicks without drags are handled
                    if mev[:xdata] == pss_dragstart[1] && mev[:ydata] == pss_dragstart[2]
                        #Deselect mode
                        pss_dehighlight()
                    end
                elseif pss_button == 2
                    #Middle release
                    pss_pan(mev[:xdata]-pss_dragstart[1],mev[:ydata]-pss_dragstart[2])
                end
            end
        end
    end

    #Mouse scroll state
    pss_zoomsteps::Int64 = 0

    #Scroll event
    function onscroll(mev::PyObject)
        axes = mev[:inaxes]
        if axes != nothing
            #Zoom on scroll
            pss_zoomsteps += mev[:step]
        end
    end

    #PSS closed state
    pss_closed::Bool = false

    #Close event
    function onclose(cev::PyObject)
        pss_closed = true
    end

    # #########################################################################
    #Connect handlers and show plot

    #Connect event handlers
    fig[:canvas][:mpl_connect]("button_press_event",onpress)
    fig[:canvas][:mpl_connect]("button_release_event",onrelease)
    fig[:canvas][:mpl_connect]("scroll_event",onscroll)
    fig[:canvas][:mpl_connect]("close_event",onclose)

    #Show plot in the background
    plt[:ion]()
    plt[:show]()
    plt[:ioff]()

    #Update zoom level at regular intervals
    while !pss_closed
        if pss_zoomsteps != 0
            zoomsteps::Int64 = pss_zoomsteps
            pss_zoomsteps = 0
            pss_zoom(zoomsteps)
        end
        sleep(1)
    end

end


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

const SDL_MIX_MAXVOLUME = 128
const SDL_MAJOR_VERSION = 2
const SDL_MINOR_VERSION = 0
const SDL_PATCHLEVEL = 7

# Skipping MacroDefinition: SDL_VERSION ( x ) \
#{ ( x ) -> major = SDL_MAJOR_VERSION ; ( x ) -> minor = SDL_MINOR_VERSION ; ( x ) -> patch = SDL_PATCHLEVEL ; \
#}
# Skipping MacroDefinition: SDL_VERSIONNUM ( X , Y , Z ) ( ( X ) * 1000 + ( Y ) * 100 + ( Z ) )
# Skipping MacroDefinition: SDL_COMPILEDVERSION SDL_VERSIONNUM ( SDL_MAJOR_VERSION , SDL_MINOR_VERSION , SDL_PATCHLEVEL )
# Skipping MacroDefinition: SDL_VERSION_ATLEAST ( X , Y , Z ) ( SDL_COMPILEDVERSION >= SDL_VERSIONNUM ( X , Y , Z ) )

const SDL_MIXER_MAJOR_VERSION = 2
const SDL_MIXER_MINOR_VERSION = 0
const SDL_MIXER_PATCHLEVEL = 2

# Skipping MacroDefinition: SDL_MIXER_VERSION ( X ) \
#{ ( X ) -> major = SDL_MIXER_MAJOR_VERSION ; ( X ) -> minor = SDL_MIXER_MINOR_VERSION ; ( X ) -> patch = SDL_MIXER_PATCHLEVEL ; \
#}

const MIX_MAJOR_VERSION = SDL_MIXER_MAJOR_VERSION
const MIX_MINOR_VERSION = SDL_MIXER_MINOR_VERSION
const MIX_PATCHLEVEL = SDL_MIXER_PATCHLEVEL

# Skipping MacroDefinition: MIX_VERSION ( X ) SDL_MIXER_VERSION ( X )
# Skipping MacroDefinition: SDL_MIXER_COMPILEDVERSION SDL_VERSIONNUM ( SDL_MIXER_MAJOR_VERSION , SDL_MIXER_MINOR_VERSION , SDL_MIXER_PATCHLEVEL )
# Skipping MacroDefinition: SDL_MIXER_VERSION_ATLEAST ( X , Y , Z ) ( SDL_MIXER_COMPILEDVERSION >= SDL_VERSIONNUM ( X , Y , Z ) )

const MIX_CHANNELS = 8
const MIX_DEFAULT_FREQUENCY = 22050
const MIX_DEFAULT_FORMAT = AUDIO_S16LSB
const MIX_DEFAULT_CHANNELS = 2
const MIX_MAX_VOLUME = SDL_MIX_MAXVOLUME

# Skipping MacroDefinition: Mix_LoadWAV ( file ) Mix_LoadWAV_RW ( SDL_RWFromFile ( file , "rb" ) , 1 )
Mix_LoadWAV(file::String) = Mix_LoadWAV_RW(RWFromFile( file , "rb" ),Int32(1))

const MIX_CHANNEL_POST = -2
const MIX_EFFECTSMAXSPEED = "MIX_EFFECTSMAXSPEED"

# Skipping MacroDefinition: Mix_PlayChannel ( channel , chunk , loops ) Mix_PlayChannelTimed ( channel , chunk , loops , - 1 )
# Skipping MacroDefinition: Mix_FadeInChannel ( channel , chunk , loops , ms ) Mix_FadeInChannelTimed ( channel , chunk , loops , ms , - 1 )
Mix_PlayChannel(channel , chunk , loops) = Mix_PlayChannelTimed( channel , chunk , loops , Int32(-1) )
Mix_FadeInChannel(channel , chunk , loops , ms) = Mix_FadeInChannelTimed( channel , chunk , loops , ms , Int32(-1) )

# const Mix_SetError = SDL_SetError
# const Mix_GetError = SDL_GetError
# const Mix_ClearError = SDL_ClearError

mutable struct Mix_Chunk
    allocated::Cint
    abuf::Ptr{Uint8}
    alen::Uint32
    volume::Uint8
end

# begin enum ANONYMOUS_9
const ANONYMOUS_9 = UInt32
const MIX_NO_FADING = (UInt32)(0)
const MIX_FADING_OUT = (UInt32)(1)
const MIX_FADING_IN = (UInt32)(2)
# end enum ANONYMOUS_9

Mix_Fading = UInt32

# begin enum ANONYMOUS_10
const ANONYMOUS_10 = UInt32
const MUS_NONE = (UInt32)(0)
const MUS_CMD = (UInt32)(1)
const MUS_WAV = (UInt32)(2)
const MUS_MOD = (UInt32)(3)
const MUS_MID = (UInt32)(4)
const MUS_OGG = (UInt32)(5)
const MUS_MP3 = (UInt32)(6)
const MUS_MP3_MAD_UNUSED = (UInt32)(7)
const MUS_FLAC = (UInt32)(8)
const MUS_MODPLUG_UNUSED = (UInt32)(9)
# end enum ANONYMOUS_10

Mix_MusicType = UInt32

mutable struct _Mix_Music
end

Mix_Music = _Mix_Music
Mix_EffectFunc_t = Ptr{Cvoid}
Mix_EffectDone_t = Ptr{Cvoid}

# begin enum ANONYMOUS_8
const ANONYMOUS_8 = UInt32
const MIX_INIT_FLAC = (UInt32)(1)
const MIX_INIT_MOD = (UInt32)(2)
const MIX_INIT_MP3 = (UInt32)(8)
const MIX_INIT_OGG = (UInt32)(16)
const MIX_INIT_MID = (UInt32)(32)
# end enum ANONYMOUS_8

MIX_InitFlags = UInt32

max_stars_repo_name: Kyjor/julia-game

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: Kyjor/julia-game

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: Kyjor/julia-game

max_forks_repo_licenses: MIT Unlicense


text:
using InteractMithril
using Test

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


text:
import Base.show
import Base.isapprox
using Folds

export
    Factor,

    show,
    product,
    sum_over,
    nextkey,
    normalize

"""
    struct Factor{N}

Representation of a factor over `N` instances.

# arguments
    dims a tuple of dimensions of the instances
    keys ids of the instances
    entries a vector of factor values
"""
struct Factor{N}
    dims :: Union{NTuple{N, Int}, Nothing}
    keys :: Union{NTuple{N, Int}, Nothing}
    entries :: Array{<:Real, 1}
end

"""
    function normalize(f::Factor)

Return a new factor equal to the given factor except that entries sum to 1
"""
function normalize(f::Factor)
    z = sum(f.entries)
    new_entries = [e / z for e in f.entries]
    return Factor(f.dims, f.keys, new_entries)
end

"""
    nextkey()::Int

Produce a fresh instance id that does not conflict with an existing id.
"""
nextkey = begin
    count = 0
    f() = begin
        global count
        count += 1
        return count
    end
    f
end

function Base.isapprox(fact1 :: Factor, fact2 :: Factor, epsilon = 0.0001)
    if fact1.dims != fact2.dims || fact1.keys != fact2.keys
        return false
    end
    len = length(fact1.entries)
    return all(map(i -> isapprox(fact2.entries[i], fact1.entries[i]; atol = epsilon), 1:len))
end

###########################
#                         #
# Pretty printing factors #
#                         #
###########################
function show_boundary(N, space)
    total = (space + 1) * (N + 1) + 1
    for i in 1:total
        print("#")
    end
    println()
end

get_content(content :: String) = content

get_content(content) = sprint(show, content)

function show_spaced(space, content)
    c = get_content(content)
    init_gap = div(space - length(c), 2)
    for i in 1:init_gap
        print(' ')
    end
    print(c)
    for i in 1:space - length(c) - init_gap
        print(' ')
    end
end

function show_content(N, space, content)
    for i in 1:N
        print('#')
        show_spaced(space, content[i])
    end
    print('#')
    for i in 1:space
        print(' ')
    end
    print('#')
    println()
end

"""
    function show(f::Factor)

Print the factor in an easy to read format.
"""
function show(f :: Factor{0})
    # This factor should have been produced by summing out all the variables
    # and should have one row
    space = length(string(f.entries[1])) + 2
    show_boundary(0, space)
    print('#')
    show_spaced(space, f.entries[1])
    println('#')
    show_boundary(0, space)
end

function show(f :: Factor{N}) where N
    max1 = max(map(k -> length(string(k)), f.keys)...)
    max2 = max(map(k -> length(string(k)), f.entries)...)
    space = max(max1, max2) + 2
    show_boundary(N, space)
    show_content(N, space, f.keys)
    show_boundary(N, space)
    mults = Array{Int64}(undef, N)
    mults[N] = 1
    for k in N-1:-1:1
        mults[k] = mults[k+1] * f.dims[k+1]
    end
    for i in eachindex(f.entries)
        vs = []
        j = i-1
        for k in 1:N
            (r, j) = divrem(j, mults[k])
            print("#")
            show_spaced(space, r + 1)
        end
        print("#")
        show_spaced(space, f.entries[i])
        println("#")
    end
    show_boundary(N, space)
end

##################
#                #
# Factor product #
#                #
##################

# The product computation is written so that it can be compiled and optimized
# once the sizes of the input factors (N1 and N2),
# dimensions of the variables (D1 and D2),
# and join indices (J) are known.
struct ProdOp{N1, N2, D1, D2, J} end

function mults(dims)
    n = length(dims)
    result = Array{Int64}(undef, n)
    result[n] = 1
    for k in n-1:-1:1
        result[k] = result[k+1] * dims[k+1]
    end
    return result
end

function do_prod(op :: ProdOp{N1, N2, D1, D2, J}, f1, f2) where {N1, N2, D1, D2, J}
    # The result factor is constructed so that the variables in the second
    # input come first, in order, and then any variables in the first factor
    # that are not in the join come in order.
    # We construct DR, the dimensions of variables in the result factor,
    # as well as the keys of the result factor.
    # We work out the instructions,
    # which will be used to indicate which rows of the two input factors
    # any row of the result factor will be composed from.
    instructions = []
    NR = N2
    rda = [d for d in D2]
    rk = [k for k in f2.keys]
    for j = 1:N2
        push!(instructions, [(2,j)])
    end

    for (i, j) in J
        if j == 0
            NR += 1
            push!(instructions, [(1, i)])
            push!(rda, D1[i])
            push!(rk, f1.keys[i])
        else
            push!(instructions[j], (1, i))
        end
    end

    DR = NTuple{NR, Int}(rda)
    mults1 = mults(f1.dims)
    mults2 = mults(f2.dims)
    multsr = mults(DR)
    rkeys = NTuple{NR, Int}(rk)
    result = Array{Float64}(undef, Folds.reduce(*, DR))
    for i in eachindex(result)
        # Cartesian indices don't work here because you can't get inside them
        # So we have to build up the index ourselves
        idx1 = Array{Int16}(undef, N1)
        idx2 = Array{Int16}(undef, N2)
        j = i-1
        for k = 1:NR
            instr = instructions[k]
            (r, j) = divrem(j, multsr[k])
            e = r + 1
            if length(instr) == 2
                idx2[instr[1][2]] = e
                idx1[instr[2][2]] = e
            elseif instr[1][1] == 1
                idx1[instr[1][2]] = e
            else
                idx2[instr[1][2]] = e
            end
        end
        x1 = sum(map(q -> (idx1[q] - 1) * mults1[q], 1:N1)) + 1
        x2 = sum(map(q -> (idx2[q] - 1) * mults2[q], 1:N2)) + 1
        result[i] = f1.entries[x1] * f2.entries[x2]
    end

    return Factor(DR, rkeys, result)
end

function product(f1 :: Factor{N1}, f2 :: Factor{N2}) where {N1, N2}
    js = map(k -> findfirst(x -> x == k, f2.keys), f1.keys)
    ijs :: Array{Tuple{Int16, Int16}, 1}= []
    for i = length(f1.keys):-1:1
        y = isnothing(js[i]) ? 0 : js[i]
        push!(ijs, (i, y))
    end
    matches = NTuple{N1, Tuple{Int16, Int16}}(ijs)
    op = ProdOp{N1, N2, f1.dims, f2.dims, matches}()
    return do_prod(op, f1, f2)
end

###########################
#                         #
# Summing over a variable #
#                         #
###########################

struct SumOp{N, D, I} end

# Special case: Summing out the only variable in a factor
function do_sum(op :: SumOp{1, D, 1}, f :: Factor) where D
    x = sum(f.entries)
    dims = nothing
    keys = nothing
    return Factor{0}(dims, keys, [x])
end

function do_sum(op :: SumOp{N, D, I}, f :: Factor) where {N, D, I}
    rdims = []
    rkeys = []
    for i in 1:I-1
        push!(rdims, f.dims[i])
        push!(rkeys, f.keys[i])
    end
    for i in I+1:N
        push!(rdims, f.dims[i])
        push!(rkeys, f.keys[i])
    end
    rdims = Tuple(rdims)
    rkeys = Tuple(rkeys)

    # The summation is efficeintly implemented by computing a pattern of
    # moves through the input factor's entries as we accumulate the result
    # factor's entries.
    # For any given entry of the result, we add a number of entries of the
    # input. inner_skip is the gap between those entries.
    # We also do inner_skip parallel sums per block, before skipping to the next
    # block. outer_skip is the amount to skip to the next block.
    inner_skip = 1
    for k = I+1:N
        inner_skip *= D[k]
    end
    outer_skip = inner_skip * D[I]
    result_size = inner_skip
    for k = 1:I-1
        result_size *= D[k]
    end
    result = Array{Float64}(undef, result_size)
    if result_size == 0 # This can legitimately happen for unexpanded Expander
        return Factor(rdims, rkeys, result)
    end

    section_start = 1
    result_index = 1
    for j = 1 : div(result_size, inner_skip)
        start = section_start
        for k = 1 : inner_skip
            total = 0.0
            orig_index = start
            for l = 1 : D[I]
                total += f.entries[orig_index]
                orig_index += inner_skip
            end
            result[result_index] = total
            result_index += 1
            start += 1
        end
        section_start += outer_skip
    end

    return Factor(rdims, rkeys, result)
end

function sum_over(f :: Factor{N}, key :: Int) where N
    for i = 1:N
        if f.keys[i] == key
            return do_sum(SumOp{N, f.dims, i}(), f)
        end
    end
end


text:
# Prototype - copied from pluto_RayTracingWeekend.jl
# Adapted from [Ray Tracing In One Weekend by Peter Shirley](https://raytracing.github.io/books/RayTracingInOneWeekend.html) and
# [cshenton's Julia implementation](https://github.com/cshenton/RayTracing.jl)"

using RayTracingWeekend

using BenchmarkTools
using Images
using InteractiveUtils
using LinearAlgebra
using StaticArrays
using LoopVectorization: @turbo

# Float32/Float64
ELEM_TYPE = Float64

t_default_cam = default_camera(SA{ELEM_TYPE}[0,0,0])

t_cam1 = default_camera([13,2,3], [0,0,0], [0,1,0], 20, 16/9, 0.1, 10.0; elem_type=ELEM_TYPE)

t_cam2 = default_camera([3,3,2], [0,0,-1], [0,1,0], 20, 16/9, 2.0, norm([3,3,2]-[0,0,-1]);
						elem_type=ELEM_TYPE)

# After some optimization:
#  46.506 ms (917106 allocations: 16.33 MiB)
# Using convert(Float32, ...) instead of MyFloat(...):
#  14.137 ms (118583 allocations: 4.14 MiB)
# Don't specify return value of Option{HitRecord} in hit()
#  24.842 ms (527738 allocations: 16.63 MiB)
# Don't initialize unnecessary elements in HitRecord():
#  14.862 ms (118745 allocations: 4.15 MiB)  (but computer was busy...)
# Replace MyFloat by Float32:
#  11.792 ms ( 61551 allocations: 2.88 MiB)
# Remove ::HitRecord return value in remaining hit() method:
#  11.545 ms ( 61654 allocations: 2.88 MiB)
# with mutable HitRecord
#  11.183 ms ( 61678 allocations: 2.88 MiB) (insignificant?)
# @inline tons of stuff. Note: render() uses `for i in 1:image_height, j in 1:image_width`,
#    i.e. iterating 1 row at time!
#   8.129 ms ( 61660 allocations: 2.88 MiB)
# Using in render(): `for j in 1:image_width, i in 1:image_height # iterate over each column`
#  10.489 ms ( 61722 allocations: 2.88 MiB) (consistently slower!)
# ... sticking with `for i in 1:image_height, j in 1:image_width # iterate over each row` for now...
# Using `get_ray(cam, u+δu, v+δv)` (fixes minor bug, extract constants outside inner loop):
# ... performance appears equivalent, maybe a tiny bit faster on avg (1%?)
# Re-measured:
#   8.077 ms (61610 allocations: 2.88 MiB)
# After parameterized Vec3{T}:
# Float64: 8.344 ms (61584 allocations: 4.82 MiB)
# Float32: 8.750 ms (123425 allocations: 3.83 MiB)
# rand() using MersenneTwister _rng w/ Float64:
#   6.967 ms (61600 allocations: 4.82 MiB)
# rand() using Xoroshiro128Plus _rng w/ Float64:
#   6.536 ms (61441 allocations: 4.81 MiB)
# Above was all using 1 single thread. With 16 threads:
#   4.414 ms (61673 allocations: 4.82 MiB)
# Above was all using max bounces=4, since this looked fine to me (except the negatively scaled sphere).
# Switching to max bounces=16 to match C++ version decreased performance by 7.2%:
#   4.465 ms (65680 allocations: 5.13 MiB)
# Lots of optimizations... ending with make HitRecord non-mutable:
#   2.225 ms (445188 allocations: 34.08 MiB)
# Using non-mutable HitRecord, Union{HitRecord,Missing}, ismissing():
#   976.365 μs (65574 allocations: 5.12 MiB)
# Using @paulmelis' style of hit(): @inbounds for i in eachindex(hittables) and Union{HitRecord, Nothing}
#   951.447 μs (65574 allocations: 5.12 MiB)
render(scene_2_spheres(; elem_type=ELEM_TYPE), t_default_cam, 96, 16) # 16 samples

# Iterate over each column: 614.820 μs
# Iterate over each row: 500.334 μs
# With Rand(Float32) everywhere:
#   489.745 μs (3758 allocations: 237.02 KiB)
# With parameterized Vec3{T}:
# Float64: 530.957 μs (3760 allocations: 414.88 KiB)
# rand() using MersenneTwister _rng w/ Float64:
#   473.672 μs (3748 allocations: 413.94 KiB)
# rand() using Xoroshiro128Plus _rng w/ Float64:
#   444.399 μs (3737 allocations: 413.08 KiB)
# Above was all using 1 single thread. With 16 threads:
#   300.438 μs (3829 allocations: 420.86 KiB)
# Above was all using max bounces=4, since this looked fine to me (except the negatively scaled sphere).
# Switching to max bounces=16 to match C++ version decreased performance by 7.2%:
#   314.094 μs (4009 allocations: 434.97 KiB)
# Lots of optimizations... ending with make HitRecord non-mutable:
#   136.388 μs (28306 allocations: 2.28 MiB)
# Using non-mutable HitRecord, Union{HitRecordMissing}, ismissing():
#   102.764 μs (4314 allocations: 459.41 KiB)
# Using @paulmelis' style of hit(): @inbounds for i in eachindex(hittables) and Union{HitRecord, Nothing}
#   101.161 μs (4314 allocations: 459.41 KiB)
render(scene_2_spheres(; elem_type=ELEM_TYPE), t_default_cam, 96, 1) # 1 sample

#render(scene_4_spheres(; elem_type=ELEM_TYPE), t_default_cam, 96, 16)

#render(scene_diel_spheres(; elem_type=ELEM_TYPE), t_default_cam, 96, 16)
#render(scene_diel_spheres(), default_camera(), 320, 32)

# Hollow Glass sphere using a negative radius
#render(scene_diel_spheres(-0.5; elem_type=ELEM_TYPE), t_default_cam, 96, 16)

#render(scene_diel_spheres(; elem_type=ELEM_TYPE), default_camera((SA{ELEM_TYPE}[-2,2,1]), (SA{ELEM_TYPE}[0,0,-1]),
#																 (SA{ELEM_TYPE}[0,1,0]), ELEM_TYPE(20)), 96, 16)

#render(scene_blue_red_spheres(; elem_type=ELEM_TYPE), t_default_cam, 96, 16)

# took ~20s (really rough timing) in REPL, before optimization
# after optimization:
#   880.997 ms (48801164 allocations: 847.10 MiB)
# after switching to Float32 + reducing allocations using rand_vec3!():
#    79.574 ms (  542467 allocations:  14.93 MiB)
# after optimizing skycolor, rand*, probably more stuff I forgot...
#    38.242 ms (  235752 allocations:   6.37 MiB)
# after removing all remaining Color, Vec3, replacing them with @SVector[]...
#    26.790 ms (   91486 allocations: 2.28 MiB)
# Using convert(Float32, ...) instead of MyFloat(...):
#    25.856 ms (   70650 allocations: 1.96 MiB)
# Don't specify return value of Option{HitRecord} in hit()
# Don't initialize unnecessary elements in HitRecord():
#    38.961 ms (   70681 allocations: 1.96 MiB) # WORSE, probably because we're writing a lot more to stack...?
# Replace MyFloat by Float32:
#    36.726 ms (13889 allocations: 724.09 KiB)
# @inline lots of stuff:
#    14.690 ms (13652 allocations: 712.98 KiB)
# rand(Float32) to avoid Float64s:
#    14.659 ms (13670 allocations: 713.84 KiB)
# Re-measured:
#    14.069 ms (13677 allocations: 714.22 KiB)
# After parameterized Vec3{T}:
# Float32: 14.422 ms (26376 allocations: 896.52 KiB) # claforte: why are there 2X the allocations as previously?
# Float64: 12.772 ms (12868 allocations: 1.08 MiB) (10% speed-up! Thanks @woclass!)
# rand() using MersenneTwister _rng w/ Float64:
#    14.092 ms (12769 allocations: 1.07 MiB) (SLOWER?!)
# rand() using Xoroshiro128Plus _rng w/ Float64:
#    12.581 ms (12943 allocations: 1.08 MiB)
# Above was all using 1 single thread. With 16 threads:
#     1.789 ms (12926 allocations: 1.08 MiB) (WOW!)
# Above was all using max bounces=4, since this looked fine to me (except the negatively scaled sphere).
# Switching to max bounces=16 to match C++ version decreased performance by 7.2%:
#     2.168 ms (13791 allocations: 1.15 MiB)
# Using  bunch of @inbounds, @simd in low-level functions
#     2.076 ms (13861 allocations: 1.15 MiB)
# Lots of optimizations, up to `Using non-mutable HitRecord, Union{HitRecordMissing}, ismissing():`
#     2.042 ms (14825 allocations: 1.23 MiB)
#@btime render(scene_random_spheres(; elem_type=ELEM_TYPE), t_cam1, 96, 1)

# took 5020s in Pluto.jl, before optimizations!
# after lots of optimizations, up to switching to Float32 + reducing allocations using rand_vec3!():
#   10.862032 seconds (70.09 M allocations: 1.858 GiB, 1.53% gc time)
# after optimizing skycolor, rand*, probably more stuff I forgot...
#    4.926 s (25,694,163 allocations: 660.06 MiB)
# after removing all remaining Color, Vec3, replacing them with @SVector[]...
#    3.541 s (9074843 allocations: 195.22 MiB)
# Using convert(Float32, ...) instead of MyFloat(...):
#    3.222 s (2055106 allocations: 88.53 MiB)
# Don't specify return value of Option{HitRecord} in hit()
# Don't initialize unnecessary elements in HitRecord():
#    5.016 s (2056817 allocations: 88.61 MiB) #  WORSE, probably because we're writing a lot more to stack...?
# Replace MyFloat by Float32:
#    5.185 s (1819234 allocations: 83.55 MiB) # Increase of 1% is probably noise
# Remove normalize() in reflect() (that function assumes the inputs are normalized)
#    5.040 s (1832052 allocations: 84.13 MiB)
# @inline lots of stuff:
#    2.110 s (1823044 allocations: 83.72 MiB)
# rand(Float32) to avoid Float64s:
#    2.063 s (1777985 allocations: 81.66 MiB)
# @woclass's "Use alias instead of new struct", i.e. `const HittableList = Vector{Hittable}`:
#    1.934 s (1796954 allocations: 82.53 MiB)
# @woclass's Vec3{T} with T=Float64: (7.8% speed-up!)
#    1.800 s (1711061 allocations: 131.03 MiB)
# rand() using Xoroshiro128Plus _rng w/ Float64:
#    1.808 s (1690104 allocations: 129.43 MiB) (i.e. rand is not a bottleneck?)
# Above was all using 1 single thread. With 16 threads:
#  265.331 ms (1645402 allocations: 126.02 MiB) (WOW!)
# Above was all using max bounces=4, since this looked fine to me (except the negatively scaled sphere).
# Switching to max bounces=16 to match C++ version decreased performance by 7.2%:
#  308.217 ms (1830162 allocations: 140.12 MiB)
# Using @inbounds, @simd in low-level functions:
#  302.952 ms (1892513 allocations: 144.88 MiB)
# Convert Camera and every Material structs to non-mutable:
#  301.042 ms (1849711 allocations: 141.61 MiB)  (i.e. - unchanged)
# Adapt @Christ_Foster's Base.getproperty w/ @inline @inbounds:
#  292.603 ms (1856398 allocations: 142.12 MiB) (ran multiple times, seems like real, 3-5% speed-up)
# Eliminate the off-by-half-a-pixel offset:
#  286.873 ms (1811412 allocations: 138.69 MiB) (ran multiple times, seems like ~2.5% speed-up)
# Fixed, per-thread RNGs with fixed seeds
#  286.575 ms (1884433 allocations: 144.26 MiB) (i.e. maybe a tiny bit faster considering this fixed seed has more allocations?)
# Make HitRecord non-mutable:
#   29.733 s (937962909 allocations: 69.88 GiB) (WTF!)
# Lots of optimizations, up to `Using non-mutable HitRecord, Union{HitRecordMissing}, ismissing():`
#  306.011 ms (1884433 allocations: 144.26 MiB) (Still slower... Hum)
# Using @paulmelis' style of hit(): @inbounds for i in eachindex(hittables) and Union{HitRecord, Nothing}
#  304.877 ms (1884433 allocations: 144.26 MiB)
# Extract the scene creation from the render() call:
#  300.344 ms (1883484 allocations: 144.21 MiB)
# Separate the module's reusable functions/classes from this proto.jl
#  296.824 ms (min) (1909996 allocs: 146.23 MiB)
print("render(scene_random_spheres(; elem_type=ELEM_TYPE), t_cam1, 200, 32):")
reseed!()
_scene_random_spheres = scene_random_spheres(; elem_type=ELEM_TYPE)
@benchmark render($_scene_random_spheres, $t_cam1, 200, 32)

# After some optimization, took ~5.6 hours:
#   20171.646846 seconds (94.73 G allocations: 2.496 TiB, 1.06% gc time)
# ... however the image looked weird... too blurry
# After removing all remaining Color, Vec3, replacing them with @SVector[]... took ~3.7 hours:
#   13326.770907 seconds (29.82 G allocations: 714.941 GiB, 0.36% gc time)
# Using convert(Float32, ...) instead of MyFloat(...):
# Don't specify return value of Option{HitRecord} in hit()
# Don't initialize unnecessary elements in HitRecord():
# Took ~4.1 hours:
#   14723.339976 seconds (5.45 G allocations: 243.044 GiB, 0.11% gc time) # WORSE, probably because we're writing a lot more to stack...?
# Replace MyFloat by Float32:
# Lots of other optimizations including @inline lots of stuff:
#    6018.101653 seconds (5.39 G allocations: 241.139 GiB, 0.21% gc time) (1.67 hour)
# @woclass's rand(Float32) to avoid Float64s: (expected to provide 2.2% speed up)
# @woclass's "Use alias instead of new struct", i.e. `const HittableList = Vector{Hittable}`
# @woclass's Vec3{T} with T=Float64: (7.8% speed-up!):
#    5268.175362 seconds (4.79 G allocations: 357.005 GiB, 0.47% gc time) (1.46 hours)
# Above was all using 1 single thread. With 16 threads: (~20 minutes)
#    1210.363539 seconds (4.94 G allocations: 368.435 GiB, 10.08% gc time)
# Above was all using max bounces=4, since this looked fine to me (except the negatively scaled sphere).
# Switching to max bounces=16 to match C++ version decreased performance by 7.2%:
#    1298.522674 seconds (5.43 G allocations: 404.519 GiB, 10.18% gc time)
# Using @inbounds, @simd in low-level functions:
#    1314.510565 seconds (5.53 G allocations: 411.753 GiB, 10.21% gc time) # NOTE: difference due to randomness?
# Adapt @Christ_Foster's Base.getproperty w/ @inline @inbounds: (expect 3-5% speed-up)
# Eliminate the off-by-half-a-pixel offset: (expect ~2.5% speed-up)
# Fixed, per-thread RNGs with fixed seeds (expecting no noticeable change in speed)
#  Using 16 threads: (21m22s)
#    1282.437499 seconds (5.53 G allocations: 411.742 GiB, 10.08% gc time) (i.e. 2.5% speed-up... currently GC- and memory-bound?)
#  Using 14 threads: (21m45s)
#    1305.767627 seconds (5.53 G allocations: 411.741 GiB, 9.97% gc time)
#print("@time render(scene_random_spheres(; elem_type=ELEM_TYPE), t_cam1, 1920, 1000):")
#@time render(scene_random_spheres(; elem_type=ELEM_TYPE), t_cam1, 1920, 1000)

# Before optimization:
#  5.993 s  (193097930 allocations: 11.92 GiB)
# after disabling: `Base.getproperty(vec::SVector{3}, sym::Symbol)`
#  1.001 s  ( 17406437 allocations: 425.87 MiB)
# after forcing Ray and point() to use Float64 instead of AbstractFloat:
#  397.905 ms (6269207 allocations: 201.30 MiB)
# after forcing use of Float32 instead of Float64:
#  487.680 ms (7128113 allocations: 196.89 MiB) # More allocations... something is causing them...
# after optimizing rand_vec3!/rand_vec2! to minimize allocations:
#  423.468 ms (6075725 allocations: 158.92 MiB)
# after optimizing skycolor, rand*, probably more stuff I forgot...
#  217.853 ms (2942272 allocations: 74.82 MiB)
# after removing all remaining Color, Vec3, replacing them with @SVector[]...
#   56.778 ms (1009344 allocations: 20.56 MiB)
# Using convert(Float32, ...) instead of MyFloat(...):
#   23.870 ms (210890 allocations: 8.37 MiB)
# Replace MyFloat by Float32:
#   22.390 ms (153918 allocations: 7.11 MiB)
# Various other changes, e.g. remove unnecessary normalize
#   20.241 ms (153792 allocations: 7.10 MiB)
# @inline lots of stuff:
#   18.065 ms (153849 allocations: 7.10 MiB)
# rand(Float32) to avoid Float64s:
#   16.035 ms (153777 allocations: 7.10 MiB)
# After @woclass's Vec3{T} with T=Float64:
#   16.822 ms (153591 allocations: 11.84 MiB)
# rand() using Xoroshiro128Plus _rng w/ Float64:
#   13.469 ms (153487 allocations: 11.83 MiB)
# Above was all using 1 single thread. With 16 threads:
#    6.537 ms (153599 allocations: 11.84 MiB)
# Above was all using max bounces=4, since this looked fine to me (except the negatively scaled sphere).
# Switching to max bounces=16 to match C++ version decreased performance by 7.2%:
#    6.766 ms (161000 allocations: 12.40 MiB)
# @inbounds and @simd in low-level functions
#    6.519 ms (160609 allocations: 12.37 MiB)
render(scene_diel_spheres(; elem_type=ELEM_TYPE), t_cam2, 96, 16)

using Profile
render(scene_random_spheres(; elem_type=ELEM_TYPE), t_cam1, 16, 1)
Profile.clear_malloc_data()
render(scene_random_spheres(; elem_type=ELEM_TYPE), t_cam1, 320, 64)


text:
dt = 0.01
s = [10.0, 28.0, 8/3]
d = 3
function flow(u::Array{Float64,1})
    sigma, rho, beta = s
    x, y, z = u[1], u[2], u[3]
    return [sigma*(y - x), x*(rho - z) - y, x*y - beta*z]
end
function next(u, sn)
    k1 = dt*flow(u)
    k2 = dt*flow(u .+ k1)
    u_next = @. u + (k1 + k2)/2
	noise = sn*randn(d)
	u_next .+= noise
    return u_next
end


text:
module TestTapirBenchmarks
using Test

@testset "$file" for file in sort([
    file for file in readdir(@__DIR__) if match(r"^test_.*\.jl$", file) !== nothing
])
    include(file)
end

end


text:
import Distances
import GeometricalPredicates
import LightGraphs

using DataFrames
using Statistics
using VoronoiDelaunay

import GeometricalPredicates.getx, GeometricalPredicates.gety

using StatsBase: countmap

struct IndexedPoint2D <: GeometricalPredicates.AbstractPoint2D
    _x::Float64
    _y::Float64
    _index::Int
    IndexedPoint2D(x::Float64,y::Float64, index::Int) = new(x, y, index)
end

IndexedPoint2D() = IndexedPoint2D(0., 0., 0)
IndexedPoint2D(x::Float64, y::Float64) = IndexedPoint2D(x, y, 0)

getx(p::IndexedPoint2D) = p._x
gety(p::IndexedPoint2D) = p._y
geti(p::IndexedPoint2D) = p._index

function adjacency_list(points::AbstractArray{T, 2} where T <: Real; filter::Bool=true, n_mads::T2 where T2 <: Real =2, k_adj::Int=5,
        adjacency_type::Symbol=:triangulation, distance::TD where TD <: Distances.SemiMetric = Euclidean())
    @assert size(points, 1) == 2

    if !in(adjacency_type, (:knn, :triangulation, :both))
        error("Unknown adjacency type: $adjacency_type")
    end

    points = deepcopy(points)

    # Rescale to match the requirement for VoronoiDelaunay.jl: min_coord=1.0+eps(Float64) and max_coord=2.0-2eps(Float64)
    points .-= minimum(points, dims=2)
    points ./= maximum(points) * 1.1
    points .+= 1.01

    # Tesselation don't work with duplicated points. And in knn 0-distance points result in loop edges
    hashes = vec(mapslices(row -> "$(row[1]) $(row[2])", round.(points, digits=10), dims=1));
    is_duplicated = get.(Ref(countmap(hashes)), hashes, 0) .> 1;
    points[:, is_duplicated] .+= (rand(Float64, (2, sum(is_duplicated))) .- 0.5) .* 2e-5;

    edge_list = Matrix{Int}(undef, 2,0)
    if (adjacency_type == :triangulation) || (adjacency_type == :both)
        points_g = [IndexedPoint2D(points[:,i]..., i) for i in 1:size(points, 2)];

        tess = DelaunayTessellation2D(length(points_g), IndexedPoint2D());
        push!(tess, points_g);

        edge_list = hcat([geti.([geta(v), getb(v)]) for v in delaunayedges(tess)]...);
    end

    if (adjacency_type == :knn) || (adjacency_type == :both)
        edge_list = hcat(hcat([hcat([[i, v] for v in x[2:end]]...)
            for (i, x) in enumerate(knn(KDTree(points), points, k_adj + 1, true)[1])]...), edge_list)
        edge_list = hcat(collect.(unique(mapslices(x -> Tuple(sort(x)), edge_list, dims=1)))...)
    end

    adj_dists = Distances.colwise(distance, points[:, edge_list[1,:]], points[:, edge_list[2,:]])

    if filter
        adj_dists_log = log10.(adj_dists)
        d_threshold = median(adj_dists_log) + n_mads * mad(adj_dists_log, normalize=true)

        filt_mask = (adj_dists_log .< d_threshold)
        edge_list = edge_list[:, filt_mask]
        adj_dists = adj_dists[filt_mask]
    end

    return edge_list, adj_dists
end

adjacency_list(spatial_df::DataFrame; kwargs...) = adjacency_list(position_data(spatial_df); kwargs...)

function connected_components(adjacent_points::Array{Array{Int, 1}, 1})
    g = LightGraphs.SimpleGraph(length(adjacent_points));
    for (v1, vs) in enumerate(adjacent_points)
        for v2 in vs
            LightGraphs.add_edge!(g, v1, v2)
        end
    end

    return LightGraphs.connected_components(g);
end

max_stars_repo_name: YinuoJin/Baysor

max_issues_repo_name: YinuoJin/Baysor

max_forks_repo_name: YinuoJin/Baysor


text:
using AdFem
using LinearAlgebra

for d in [1, 2]
    mesh = Mesh(8, 8, 1/8, degree=d)
    c0 = rand(get_ngauss(mesh))
    c = constant(c0)

    out = compute_fem_source_term1(c, mesh)
    out0 = compute_fem_source_term1(c0, mesh)
    C = compute_fem_mass_matrix1(c, mesh)
    o = C * ones(mesh.ndof)
    sess = Session(); init(sess)
    @info maximum(abs.(run(sess, out)-run(sess, o)))
    @info maximum(abs.(run(sess, c)-c0))
end


text:
# Naive implementations of multiply_dimensionwise used to demonstrate the functionality
# without performance optimizations and for testing correctness of the optimized versions
# implemented below.
function multiply_dimensionwise_naive(matrix::AbstractMatrix, data_in::AbstractArray{<:Any, 2})
  size_out = size(matrix, 1)
  size_in  = size(matrix, 2)
  n_vars   = size(data_in, 1)
  data_out = zeros(promote_type(eltype(data_in), eltype(matrix)), n_vars, size_out)

  for i in 1:size_out
    for ii in 1:size_in
      for v in 1:n_vars
        data_out[v, i] += matrix[i, ii] * data_in[v, ii]
      end
    end
  end

  return data_out
end

function multiply_dimensionwise_naive(matrix::AbstractMatrix, data_in::AbstractArray{<:Any, 3})
  size_out = size(matrix, 1)
  size_in  = size(matrix, 2)
  n_vars   = size(data_in, 1)
  data_out = zeros(promote_type(eltype(data_in), eltype(matrix)), n_vars, size_out, size_out)

  for j in 1:size_out, i in 1:size_out
    for jj in 1:size_in, ii in 1:size_in
      for v in 1:n_vars
        data_out[v, i, j] += matrix[i, ii] * matrix[j, jj] * data_in[v, ii, jj]
      end
    end
  end

  return data_out
end

function multiply_dimensionwise_naive(matrix::AbstractMatrix, data_in::AbstractArray{<:Any, 4})
  size_out = size(matrix, 1)
  size_in  = size(matrix, 2)
  n_vars   = size(data_in, 1)
  data_out = zeros(promote_type(eltype(data_in), eltype(matrix)), n_vars, size_out, size_out, size_out)

  for k in 1:size_out, j in 1:size_out, i in 1:size_out
    for kk in 1:size_in, jj in 1:size_in, ii in 1:size_in
      for v in 1:n_vars
        data_out[v, i, j, k] += matrix[i, ii] * matrix[j, jj] * matrix[k, kk] * data_in[v, ii, jj, kk]
      end
    end
  end

  return data_out
end

"""
    multiply_dimensionwise(matrix::AbstractMatrix, data_in::AbstractArray{<:Any, NDIMS+1})

Multiply the array `data_in` by `matrix` in each coordinate direction, where `data_in`
is assumed to have the first coordinate for the number of variables and the remaining coordinates
are multiplied by `matrix`.
"""
function multiply_dimensionwise(matrix::AbstractMatrix, data_in::AbstractArray{<:Any, 2})
  # 1D
  # optimized version of multiply_dimensionwise_naive
  size_out = size(matrix, 1)
  n_vars   = size(data_in, 1)
  data_out = zeros(promote_type(eltype(data_in), eltype(matrix)), n_vars, size_out)

  multiply_dimensionwise!(data_out, matrix, data_in)

  return data_out
end

function multiply_dimensionwise(matrix::AbstractMatrix, data_in::AbstractArray{<:Any, 3})
  # 2D
  # optimized version of multiply_dimensionwise_naive
  size_out = size(matrix, 1)
  n_vars   = size(data_in, 1)
  data_out = zeros(promote_type(eltype(data_in), eltype(matrix)), n_vars, size_out, size_out)

  multiply_dimensionwise!(data_out, matrix, data_in)

  return data_out
end

function multiply_dimensionwise(matrix::AbstractMatrix, data_in::AbstractArray{<:Any, 4})
  # 3D
  # optimized version of multiply_dimensionwise_naive
  size_out = size(matrix, 1)
  n_vars   = size(data_in, 1)
  data_out = zeros(promote_type(eltype(data_in), eltype(matrix)), n_vars, size_out, size_out, size_out)

  multiply_dimensionwise!(data_out, matrix, data_in)

  return data_out
end

# In the following, there are several optimized in-place versions of multiply_dimensionwise.
# These make use of the macro `@tullio` from Tullio.jl, which basically uses an Einstein
# summation convention syntax.

# 1D version
function multiply_dimensionwise!(data_out::AbstractArray{<:Any, 2}, matrix::AbstractMatrix,
                                 data_in ::AbstractArray{<:Any, 2})
  @tullio threads=false data_out[v, i] = matrix[i, ii] * data_in[v, ii]

  return nothing
end

# 1D version, apply matrixJ to data_inJ
function multiply_dimensionwise!(data_out::AbstractArray{<:Any, 2}, matrix1::AbstractMatrix,
                                 data_in1::AbstractArray{<:Any, 2}, matrix2::AbstractMatrix,
                                 data_in2::AbstractArray{<:Any, 2})
  @tullio threads=false data_out[v, i] = matrix1[i, ii] * data_in1[v, ii] + matrix2[i, ii] * data_in2[v, ii]

  return nothing
end

# 2D version
function multiply_dimensionwise!(data_out::AbstractArray{<:Any, 3}, matrix::AbstractMatrix,
                                 data_in:: AbstractArray{<:Any, 3},
                                 tmp1=zeros(eltype(data_out), size(data_out, 1), size(matrix, 1), size(matrix, 2)))

  # Interpolate in x-direction
  @tullio threads=false tmp1[v, i, j]     = matrix[i, ii] * data_in[v, ii, j]

  # Interpolate in y-direction
  @tullio threads=false data_out[v, i, j] = matrix[j, jj] * tmp1[v, i, jj]

  return nothing
end

# 2D version, apply matrixJ to dimension J of data_in
function multiply_dimensionwise!(data_out::AbstractArray{<:Any, 3},
                                 matrix1::AbstractMatrix, matrix2::AbstractMatrix,
                                 data_in:: AbstractArray{<:Any, 3},
                                 tmp1=zeros(eltype(data_out), size(data_out, 1), size(matrix1, 1), size(matrix1, 2)))

  # Interpolate in x-direction
  @tullio threads=false tmp1[v, i, j]     = matrix1[i, ii] * data_in[v, ii, j]

  # Interpolate in y-direction
  @tullio threads=false data_out[v, i, j] = matrix2[j, jj] * tmp1[v, i, jj]

  return nothing
end

# 2D version, apply matrixJ to dimension J of data_in and add the result to data_out
function add_multiply_dimensionwise!(data_out::AbstractArray{<:Any, 3},
                                     matrix1::AbstractMatrix, matrix2::AbstractMatrix,
                                     data_in:: AbstractArray{<:Any, 3},
                                     tmp1=zeros(eltype(data_out), size(data_out, 1), size(matrix1, 1), size(matrix1, 2)))

  # Interpolate in x-direction
  @tullio threads=false tmp1[v, i, j]     = matrix1[i, ii] * data_in[v, ii, j]

  # Interpolate in y-direction
  @tullio threads=false data_out[v, i, j] += matrix2[j, jj] * tmp1[v, i, jj]

  return nothing
end

# 3D version
function multiply_dimensionwise!(data_out::AbstractArray{<:Any, 4}, matrix::AbstractMatrix,
                                 data_in:: AbstractArray{<:Any, 4},
                                 tmp1=zeros(eltype(data_out), size(data_out, 1), size(matrix, 1), size(matrix, 2), size(matrix, 2)),
                                 tmp2=zeros(eltype(data_out), size(data_out, 1), size(matrix, 1), size(matrix, 1), size(matrix, 2)))

  # Interpolate in x-direction
  @tullio threads=false tmp1[v, i, j, k]     = matrix[i, ii] * data_in[v, ii, j, k]

  # Interpolate in y-direction
  @tullio threads=false tmp2[v, i, j, k]     = matrix[j, jj] * tmp1[v, i, jj, k]

  # Interpolate in z-direction
  @tullio threads=false data_out[v, i, j, k] = matrix[k, kk] * tmp2[v, i, j, kk]

  return nothing
end

# 3D version, apply matrixJ to dimension J of data_in
function multiply_dimensionwise!(data_out::AbstractArray{<:Any, 4},
                                 matrix1::AbstractMatrix, matrix2::AbstractMatrix, matrix3::AbstractMatrix,
                                 data_in:: AbstractArray{<:Any, 4},
                                 tmp1=zeros(eltype(data_out), size(data_out, 1), size(matrix1, 1), size(matrix1, 2), size(matrix1, 2)),
                                 tmp2=zeros(eltype(data_out), size(data_out, 1), size(matrix1, 1), size(matrix1, 1), size(matrix1, 2)))

  # Interpolate in x-direction
  @tullio threads=false tmp1[v, i, j, k]     = matrix1[i, ii] * data_in[v, ii, j, k]

  # Interpolate in y-direction
  @tullio threads=false tmp2[v, i, j, k]     = matrix2[j, jj] * tmp1[v, i, jj, k]

  # Interpolate in z-direction
  @tullio threads=false data_out[v, i, j, k] = matrix3[k, kk] * tmp2[v, i, j, kk]

  return nothing
end

# 3D version, apply matrixJ to dimension J of data_in and add the result to data_out
function add_multiply_dimensionwise!(data_out::AbstractArray{<:Any, 4},
                                     matrix1::AbstractMatrix, matrix2::AbstractMatrix, matrix3::AbstractMatrix,
                                     data_in:: AbstractArray{<:Any, 4},
                                     tmp1=zeros(eltype(data_out), size(data_out, 1), size(matrix1, 1), size(matrix1, 2), size(matrix1, 2)),
                                     tmp2=zeros(eltype(data_out), size(data_out, 1), size(matrix1, 1), size(matrix1, 1), size(matrix1, 2)))

  # Interpolate in x-direction
  @tullio threads=false tmp1[v, i, j, k]     = matrix1[i, ii] * data_in[v, ii, j, k]

  # Interpolate in y-direction
  @tullio threads=false tmp2[v, i, j, k]     = matrix2[j, jj] * tmp1[v, i, jj, k]

  # Interpolate in z-direction
  @tullio threads=false data_out[v, i, j, k] += matrix3[k, kk] * tmp2[v, i, j, kk]

  return nothing
end

# Calculate the Dhat matrix
function calc_dhat(nodes, weights)
  n_nodes = length(nodes)
  dhat = polynomial_derivative_matrix(nodes)
  dhat = transpose(dhat)

  for n in 1:n_nodes, j in 1:n_nodes
    dhat[j, n] *= -weights[n] / weights[j]
  end

  return dhat
end

# Calculate the Dsplit matrix for split-form differentiation: dplit = 2D - M⁻¹B
function calc_dsplit(nodes, weights)
  # Start with 2 x the normal D matrix
  dsplit = polynomial_derivative_matrix(nodes)
  dsplit = 2 .* dsplit

  # Modify to account for
  dsplit[1, 1] += 1/weights[1]
  dsplit[end, end] -= 1/weights[end]

  return dsplit
end

# Calculate the polynomial derivative matrix D
function polynomial_derivative_matrix(nodes)
  n_nodes = length(nodes)
  d = zeros(n_nodes, n_nodes)
  wbary = barycentric_weights(nodes)

  for i in 1:n_nodes, j in 1:n_nodes
    if j != i
      d[i, j] = wbary[j] / wbary[i] * 1 / (nodes[i] - nodes[j])
      d[i, i] -= d[i, j]
    end
  end

  return d
end

# Calculate and interpolation matrix (Vandermonde matrix) between two given sets of nodes
function polynomial_interpolation_matrix(nodes_in, nodes_out)
  n_nodes_in = length(nodes_in)
  n_nodes_out = length(nodes_out)
  wbary_in = barycentric_weights(nodes_in)
  vdm = zeros(n_nodes_out, n_nodes_in)

  for k in 1:n_nodes_out
    match = false
    for j in 1:n_nodes_in
      if isapprox(nodes_out[k], nodes_in[j], rtol=eps())
        match = true
        vdm[k, j] = 1
      end
    end

    if match == false
      s = 0.0
      for j in 1:n_nodes_in
        t = wbary_in[j] / (nodes_out[k] - nodes_in[j])
        vdm[k, j] = t
        s += t
      end
      for j in 1:n_nodes_in
        vdm[k, j] = vdm[k, j] / s
      end
    end
  end

  return vdm
end

# Calculate the barycentric weights for a given node distribution.
function barycentric_weights(nodes)
  n_nodes = length(nodes)
  weights = ones(n_nodes)

  for j = 2:n_nodes, k = 1:(j-1)
    weights[k] *= nodes[k] - nodes[j]
    weights[j] *= nodes[j] - nodes[k]
  end

  for j in 1:n_nodes
    weights[j] = 1 / weights[j]
  end

  return weights
end

# Calculate Lhat.
function calc_lhat(x::Float64, nodes, weights)
  n_nodes = length(nodes)
  wbary = barycentric_weights(nodes)

  lhat = lagrange_interpolating_polynomials(x, nodes, wbary)

  for i in 1:n_nodes
    lhat[i] /= weights[i]
  end

  return lhat
end

# Calculate Lagrange polynomials for a given node distribution.
function lagrange_interpolating_polynomials(x::Float64, nodes, wbary)
  n_nodes = length(nodes)
  polynomials = zeros(n_nodes)

  for i in 1:n_nodes
    if isapprox(x, nodes[i], rtol=eps(x))
      polynomials[i] = 1
      return polynomials
    end
  end

  for i in 1:n_nodes
    polynomials[i] = wbary[i] / (x - nodes[i])
  end
  total = sum(polynomials)

  for i in 1:n_nodes
    polynomials[i] /= total
  end

  return polynomials
end

# From FLUXO (but really from blue book by Kopriva)
function gauss_lobatto_nodes_weights(n_nodes::Integer)
  # From Kopriva's book
  n_iterations = 10
  tolerance = 1e-15

  # Initialize output
  nodes = zeros(n_nodes)
  weights = zeros(n_nodes)

  # Get polynomial degree for convenience
  N = n_nodes - 1

  # Calculate values at boundary
  nodes[1] = -1.0
  nodes[end] = 1.0
  weights[1] = 2 / (N * (N + 1))
  weights[end] = weights[1]

  # Calculate interior values
  if N > 1
    cont1 = pi/N
    cont2 = 3/(8 * N * pi)

    # Use symmetry -> only left side is computed
    for i in 1:(div(N + 1, 2) - 1)
      # Calculate node
      # Initial guess for Newton method
      nodes[i+1] = -cos(cont1*(i+0.25) - cont2/(i+0.25))

      # Newton iteration to find root of Legendre polynomial (= integration node)
      for k in 0:n_iterations
        q, qder, _ = calc_q_and_l(N, nodes[i+1])
        dx = -q/qder
        nodes[i+1] += dx
        if abs(dx) < tolerance * abs(nodes[i+1])
          break
        end
      end

      # Calculate weight
      _, _, L = calc_q_and_l(N, nodes[i+1])
      weights[i+1] = weights[1] / L^2

      # Set nodes and weights according to symmetry properties
      nodes[N+1-i] = -nodes[i+1]
      weights[N+1-i] = weights[i+1]
    end
  end

  # If odd number of nodes, set center node to origin (= 0.0) and calculate weight
  if n_nodes % 2 == 1
    _, _, L = calc_q_and_l(N, 0)
    nodes[div(N, 2) + 1] = 0.0
    weights[div(N, 2) + 1] = weights[1] / L^2
  end

  return nodes, weights
end

# From FLUXO (but really from blue book by Kopriva)
function calc_q_and_l(N::Integer, x::Float64)
  L_Nm2 = 1.0
  L_Nm1 = x
  Lder_Nm2 = 0.0
  Lder_Nm1 = 1.0

  local L
  for i in 2:N
    L = ((2 * i - 1) * x * L_Nm1 - (i - 1) * L_Nm2) / i
    Lder = Lder_Nm2 + (2 * i - 1) * L_Nm1
    L_Nm2 = L_Nm1
    L_Nm1 = L
    Lder_Nm2 = Lder_Nm1
    Lder_Nm1 = Lder
  end

  q = (2 * N + 1)/(N + 1) * (x * L - L_Nm2)
  qder = (2 * N + 1) * L

  return q, qder, L
end
calc_q_and_l(N::Integer, x::Real) = calc_q_and_l(N, convert(Float64, x))

# From FLUXO (but really from blue book by Kopriva)
function gauss_nodes_weights(n_nodes::Integer)
  # From Kopriva's book
  n_iterations = 10
  tolerance = 1e-15

  # Initialize output
  nodes = ones(n_nodes) * 1000
  weights = zeros(n_nodes)

  # Get polynomial degree for convenience
  N = n_nodes - 1
  if N == 0
    nodes .= 0.0
    weights .= 2.0
    return nodes, weights
  elseif N == 1
    nodes[1] = -sqrt(1/3)
    nodes[end] = -nodes[1]
    weights .= 1.0
    return nodes, weights
  else # N > 1
    # Use symmetry property of the roots of the Legendre polynomials
    for i in 0:(div(N + 1, 2) - 1)
      # Starting guess for Newton method
      nodes[i+1] = -cos(pi / (2 * N + 2) * (2 * i + 1))

      # Newton iteration to find root of Legendre polynomial (= integration node)
      for k in 0:n_iterations
        poly, deriv = legendre_polynomial_and_derivative(N + 1, nodes[i+1])
        dx = -poly / deriv
        nodes[i+1] += dx
        if abs(dx) < tolerance * abs(nodes[i+1])
          break
        end
      end

      # Calculate weight
      poly, deriv = legendre_polynomial_and_derivative(N + 1, nodes[i+1])
      weights[i+1] = (2 * N + 3) / ((1 - nodes[i+1]^2) * deriv^2)

      # Set nodes and weights according to symmetry properties
      nodes[N+1-i] = -nodes[i+1]
      weights[N+1-i] = weights[i+1]
    end

    # If odd number of nodes, set center node to origin (= 0.0) and calculate weight
    if n_nodes % 2 == 1
      poly, deriv = legendre_polynomial_and_derivative(N + 1, 0.0)
      nodes[div(N, 2) + 1] = 0.0
      weights[div(N, 2) + 1] = (2 * N + 3) / deriv^2
    end

    return nodes, weights
  end
end

# From FLUXO (but really from blue book by Kopriva)
function legendre_polynomial_and_derivative(N::Int, x::Real)
  if N == 0
    poly = 1.0
    deriv = 0.0
  elseif N == 1
    poly = convert(Float64, x)
    deriv = 1.0
  else
    poly_Nm2 = 1.0
    poly_Nm1 = convert(Float64, x)
    deriv_Nm2 = 0.0
    deriv_Nm1 = 1.0

    poly = 0.0
    deriv = 0.0
    for i in 2:N
      poly = ((2*i-1) * x * poly_Nm1 - (i-1) * poly_Nm2) / i
      deriv=deriv_Nm2 + (2*i-1)*poly_Nm1
      poly_Nm2=poly_Nm1
      poly_Nm1=poly
      deriv_Nm2=deriv_Nm1
      deriv_Nm1=deriv
    end
  end

  # Normalize
  poly = poly * sqrt(N+0.5)
  deriv = deriv * sqrt(N+0.5)

  return poly, deriv
end

# Calculate Legendre vandermonde matrix and its inverse
function vandermonde_legendre(nodes, N)
  n_nodes = length(nodes)
  n_modes = N + 1
  vandermonde = zeros(n_nodes, n_modes)

  for i in 1:n_nodes
    for m in 1:n_modes
      vandermonde[i, m], _ = legendre_polynomial_and_derivative(m-1, nodes[i])
    end
  end
  # for very high polynomial degree, this is not well conditioned
  inverse_vandermonde = inv(vandermonde)
  return vandermonde, inverse_vandermonde
end
vandermonde_legendre(nodes) = vandermonde_legendre(nodes, length(nodes) - 1)


text:
# p19.jl - 2nd-order wave eq. on Chebyshev grid (compare p6.jl)

# Time-stepping by leap frog formula:
N = 80; x = cheb(N)[2]; dt = 8/N^2;
v = @. exp(-200*x^2); vold = @. exp(-200*(x-dt)^2);
tmax = 4; tplot = .075;
plotgap = round(Int,tplot/dt); dt = tplot/plotgap;
nplots = round(Int,tmax/tplot);
plotdata = [v  zeros(N+1,nplots)]; tdata = 0;
for i = 1:nplots
    for n = 1:plotgap
        w = chebfft(chebfft(v)); w[1] = 0; w[N+1] = 0;
        vnew = 2*v - vold + dt^2*w; vold = v; v = vnew;
    end
    plotdata[:,i+1] = v; tdata=[tdata;dt*i*plotgap];
end

# Plot results:
clf(); mesh(x,tdata,plotdata',ccount=0,rcount=N+1);
axis([-1,1,0,tmax]); zlim(-2,2); view(10,70); grid(false);
xlabel("x"); ylabel("t"); zlabel("u");


text:
"""
	cubicbezier2D(curvePts::Array{Array{Float,1},1})

Produce the two ``coordinate functions`` for a ``cubic`` *Bézier curve* in ``2D``.
The input is the array `curvePts` of 4 control points with 2 coordinates.

# Example

```
julia> curvePts = [[1232.24, 1340.84],[1259.53, 1119.86],
	[1417.91, 1015.16],[1133.47, 1010.63]]

julia> bx,by = cubicbezier2D(curvePts);

julia> (bx(0), by(0))
(1232.24, 1340.84)

julia> (bx(1), by(1))
(1133.47, 1010.63)
```
"""
function cubicbezier2D( curvePts::Array{Array{Float64,1},1} )
	b1(u) = (1 - u)^3
	b2(u) = 3*u*(1 - u)^2
	b3(u) = 3*u^2*(1 - u)
	b4(u) = u^3
	cntrlverts = hcat(curvePts...)
	x = cntrlverts[1,:]
	y = cntrlverts[2,:]
	Bx = u -> x[1]*b1(u) + x[2]*b2(u) + x[3]*b3(u) + x[4]*b4(u)
	By = u -> y[1]*b1(u) + y[2]*b2(u) + y[3]*b3(u) + y[4]*b4(u)
	return Bx,By
end

"""
	const cmdsplit

Regex for splitting the `d` element of a `<path` element into
a sequence of graphics commnds.
"""
const cmdsplit = r"\s*([mMzZlLhHvVcCsSqQtTaA])\s*"

"""
	const digitRegEx

Regex for extracting numeric params from from a sequence of graphics
commnds of a `<path` element.
"""
const digitRegEx = r"\s*(-?[0-9]*\.?\d+)\s*"

"""
	pathparse(elements)

Parse graphics commands in SVG `<path` tagged element.

# Example

```julia

```
"""
function pathparse(data, cnrtlpolygon=false)
	tokens = split(data,'\"')
	pathstring = tokens[2]

	iterator = eachmatch(cmdsplit, pathstring)
	offsets = [m.offset for m in iterator]
	commands = [m.match for m in iterator]

	substrings = [pathstring[offsets[k]+1:offsets[k+1]-1]
		for k=1:length(offsets)-1]
	push!(substrings, pathstring[offsets[end]+1:end])

	params = Array{Float64,1}[]
	for substring in substrings
		numbers = []
		for match in eachmatch(digitRegEx, substring)
			push!(numbers,String(match.match))
		end
		numbers = map(x->parse(Float64,x),numbers)
		push!(params,numbers)
	end

	lines = Array{Float64,1}[]
	global startpoint, endpoint = 0.0, 0.0
	for (command, args) in zip(commands, params)
		if command == "M"
			global startpoint = args
		elseif command == "L"
			endpoint = args
			line = vcat([startpoint,endpoint]...)
			startpoint = endpoint
			push!(lines, line)
		elseif command == "C"
			contrlpt1 = args[1:2]
			contrlpt2 = args[3:4]
			endpoint = args[5:6]
			curvePts = [startpoint,contrlpt1,contrlpt2,endpoint]
			#println(curvePts)

			Bx,By = cubicbezier2D(curvePts)
			pts = [[Bx(u),By(u)] for u=0:.1:1]
			curvelines = [ vcat([pts[k],pts[k+1]]...) for k=1:length(pts)-1 ]
			for line in curvelines
				push!(lines, line)
			end

 			if cnrtlpolygon
				push!(lines,vcat(curvePts[1:2]...),vcat(curvePts[2:3]...),
					vcat(curvePts[3:4]...))
			end
			startpoint = endpoint
		end

	end
	return lines
end

"""
	lines2lar(lines)

LAR model construction from array of float quadruples.
Each `line` in input array stands for `x1,y1,x2,y2`.
"""
function lines2lar(lines)
	vertdict = OrderedDict{Array{Float64,1}, Int64}()
	EV = Array{Int64,1}[]
	idx = 0
	for h=1:size(lines,2)
		x1,y1,x2,y2 = lines[:,h]

		if ! haskey(vertdict, [x1,y1])
			idx += 1
			vertdict[[x1,y1]] = idx
		end
		if ! haskey(vertdict, [x2,y2])
			idx += 1
			vertdict[[x2,y2]] = idx
		end
		v1,v2 = vertdict[[x1,y1]],vertdict[[x2,y2]]
		push!(EV, [v1,v2])
	end
	V = hcat(collect(keys(vertdict))...)
	return V,EV
end

"""
	normalize(V::Lar.Points; flag=true::Bool)::Lar.Points

2D normalization transformation (isomorphic by defaults) of model
vertices to normalized coordinates ``[0,1]^2``. Used with SVG importing.
"""
function normalize(V::Lar.Points; flag=true)
	m,n = size(V)
	if m > n # V by rows
		V = convert(Lar.Points, V')
	end

	xmin = minimum(V[1,:]); ymin = minimum(V[2,:]);
	xmax = maximum(V[1,:]); ymax = maximum(V[2,:]);
	box = [[xmin; ymin] [xmax; ymax]]	# containment box
	aspectratio = (xmax-xmin)/(ymax-ymin)
	if flag
		if aspectratio > 1
			umin = 0; umax = 1
			vmin = 0; vmax = 1/aspectratio; ty = vmax
		elseif aspectratio < 1
			umin = 0; umax = aspectratio
			vmin = 0; vmax = 1; ty = vmax
		end
		T = Lar.t(0,ty) * Lar.s(1,-1) * Lar.s((umax-umin), (vmax-vmin)) *
			Lar.s(1/(xmax-xmin),1/(ymax-ymin)) * Lar.t(-xmin,-ymin)
	else
		# T = Lar.t(0, ymax-ymin) * Lar.s(1,-1)
		T = Lar.s(1,-1)
	end
	dim = size(V,1)
	W = T[1:dim,:] * [V;ones(1,size(V,2))]
	#V = map( x->round(x,digits=8), W )
	V = map(Lar.approxVal(8), W)

	if m > n # V by rows
		V = convert(Lar.Points, V')
	end
	return V
end

"""
	svg2lar(filename::String; flag=true)::Lar.LAR

Parse a SVG file to a `LAR` model `(V,EV)`.
Only  `<line >` and `<rect >` and `<path >` SVG primitives are currently translated.
TODO:  interpretation of transformations.
"""
function svg2lar(filename::String; flag=true)::Lar.LAR
	outlines = Array{Float64,1}[]
	matchall(r::Regex, s::AbstractString; overlap::Bool=false) =
		collect(( m.match for m=eachmatch(r,s,overlap=overlap) ));

	for line in eachline(filename)
		parts = split(lstrip(line), ' ')
		elements = [part for part in parts if part≠""]
		tag = elements[1]
		# SVG <line > primitives
		if tag == "<line"
			r = r"(.)(.)="
			regex = r"([0-9]*?\.[0-9]*)"
			prefixes = matchall(r,line)
			values = matchall(regex,line)
			outline = map(Meta.parse, values)
			#exprs = [ eval(Meta.parse(*(pre,val))) for (pre,val) in zip(prefixes, values) ]
			push!(outlines, outline)
		# SVG <rect > primitives
		elseif tag == "<rect"
			r = r"(.)(.)="
			regex = r"([0-9]*?\.[0-9]*)"
			prefixes = matchall(r,line)
			values = matchall(regex,line)
			x, y, width, height = map(Meta.parse, values)
			# regex = r"""(<rect x=")(.+?)(" y=")(.+?)(" )(.*?)( width=")(.+?)(" height=")(.+?)("/>)"""
			# coords = collect(match( regex , line)[k] for k in (4,6,8,10))
			# x, y, width, height = [ parse(Float64, string) for string in coords ]
			line1 = [ x, y, x+width, y ]
			line2 = [ x, y, x, y+height ]
			line3 = [ x+width, y, x+width, y+height ]
			line4 = [ x, y+height, x+width, y+height ]
			push!(outlines, line1, line2, line3, line4)
		# SVG <path  > primitives
		# see https://github.com/chebfun/chebfun/issues/1617
		elseif tag == "<path"
			dataregex = r"""d=(".*?")(.*?)"""
			data = string(match( dataregex , line)[1])
			polyline = pathparse(data)
			append!(outlines, polyline)
		end
	end
	lines = hcat(outlines...)
	lines = map( x->round(x,sigdigits=8), lines )
	# LAR model construction
	V,EV = lines2lar(lines)
	# normalization
	V = normalize(V,flag=flag)
	return V,EV
end