text:
#base functions for useful things

#creates evenly spaced points for qqplots
ppoints(n) = (collect(1:n) .- 0.5) ./ n

#finds closest even number above or at current number
function nearestEven(num::Union{Float64, Int64})
    return (num % 2 == 0) ? num : (num + 1);
end

"""

    qq(pvalues::AbstractArray)

# Position arguments

- `pvalues::AbstractArray`: pvalues. A one dimensional array containing pvalues to be used in the qqplot.

    qq(df::DataFrame)

# Position arguments

- `df::DataFrame`: DataFrame containing pvalues to be used in the qqplot. Note: The column of the dataframe
that indicates pvalues must be named pval (df[!, :pval] must exist)

# Keyword arguments

- `titles::AbstractString`: Title for the plot. Default is "QQ Plot of GWAS p-values".
 To have blank enter "".

- `outfile::AbstractString`: output name to save the QQplot. Name should end in format.
Default is "qqplot.png". Supports .png, .pdf, and .svg files.

- `dpi::Union{Float64, Int64}`: dots per inch to save the png file. Higher DPI results in larger file with
higher resolution. Default dpi is 350.

- `xlabel::AbstractString`: option to replace x-label text

- `ylabel::AbstractString`: option to replace y-label text

- `xmin::Union{Float64, Int64}`: Specified minimum x value to represent on the plot

- `xmax::Union{Float64, Int64}`: Specified maximum x value to represent on the plot

- `ymin::Union{Float64, Int64}`: Specified minimum y value to represent on the plot

- `ymax::Union{Float64, Int64}`: Specified maximum y value to represent on the plot

- `xstep`: Step-size for x-axis label ticks. Default value is 1.

- `ystep`: Step-size for y-axis label ticks. Default value is 2.

- `gifdist`: Distance from right edge to print the genomic inflation factor `λ`. Default is 1.0.

- `linecolor::AbstractString`: Color of "normal" line. Default color is "red".

- `dotcolor::AbstractString`: Color of the dots. Default color is "black".

- `fontsize` size of the axis labels. Default is "15pt".

- `pvalvar` variable indicating pvalue column name (for dataframes only). Default is "pval".

"""
function qq(pvalues::AbstractArray;
    titles::AbstractString = "QQ Plot of GWAS p-values",
    outfile::AbstractString = "qqplot.png",
    dpi::Int64 = 350,
    xlabel::AbstractString = "Expected -log10(p)",
    ylabel::AbstractString = "Observed -log10(p)",
    xmin::Union{Float64, Int64} = 0.0,
    ymin::Union{Float64, Int64} = 0.0,
    xmax::Union{Float64, Int64} = 0.0,
    ymax::Union{Float64, Int64} = 0.0,
    xstep::Union{Float64, Int64} = 1,
    ystep::Union{Float64, Int64} = 2,
    gifdist::Union{Float64, Int64} = 1.0,
    linecolor::AbstractString = "red",
    dotcolor::AbstractString = "black",
    fontsize = 15pt, kwargs...)

    format = lowercase(split(outfile, ".")[end])
    if !(format in ["png", "pdf", "svg"])
        throw(ArgumentError("Incorrect file format specified in outfile option.
         outfile should end in .png, .pdf, or .svg (very slow)."))
    end

    all(pvalues .>= 0.0) && all(pvalues .<= 1.0) ||
    @error("P values that are not between 0 and 1 inputted.")

    N = length(pvalues)
    up = Array{Float64}(undef, N)
    low = Array{Float64}(undef, N)
    for i in 1:N
        low[i] = -log10(quantile(Beta(i, N - i + 1), 0.975))
        up[i] = -log10(quantile(Beta(i, N - i + 1), 0.025))
    end
    up = sort!(up)
    obs = -log10.(sort(pvalues))
    obs = clamp.(obs, 0.0, -log10(floatmin(0.0)))
    expect = -log10.(ppoints(length(pvalues)))

    civals = [low expect; up sort(expect)]

    if xmax == 0.0
        xmax = ceil(maximum(expect))
    end
    if ymax == 0.0
        ymax = ceil(maximum(obs))
    end
    xticks = collect(xmin:xstep:xmax)
    yticks = collect(ymin:ystep:nearestEven(ymax))
    xλ = xmax - gifdist
    #calculate genomic inflation factor
    λ = cquantile(Chisq(1), median(pvalues)) / quantile(Chisq(1), 0.5)
    pCI = layer(x = civals[:, 2], y = civals[:, 1], Geom.polygon(fill = true,
    preserve_order = true), Theme(panel_fill = nothing, grid_line_width = 0mm,
     lowlight_color = c -> RGB{Float32}(0.0, 0.0, 0.0), alphas = fill(0.2, N),
     discrete_highlight_color = c -> RGB{Float32}(1.0, 1.0, 1.0),
     line_width = 0mm));

    pmain = layer(x = expect, y = obs, Geom.point, intercept=[0], slope = [1],
        Geom.abline(color = linecolor, style = :dash),
        Theme(panel_fill = nothing, highlight_width = 0mm, point_size = 0.5mm,
        key_position = :none, grid_line_width = 0mm, panel_stroke = colorant"black",
        default_color = dotcolor); kwargs..., order = 1);

    plt1 = plot(pCI, pmain, Scale.x_continuous(minvalue = xmin, maxvalue = xmax),
        Scale.y_continuous(minvalue = ymin, maxvalue = ymax), Guide.title(titles),
        Guide.xlabel(xlabel), Guide.ylabel(ylabel), Guide.xticks(ticks = xticks),
        Guide.yticks(ticks = yticks),
        Guide.annotation(compose(context(), text(xλ, 0, "λ = " * string(round(λ, sigdigits = 4))))),
        Theme(panel_fill = nothing, highlight_width = 0mm, point_size = 0.5mm,
        key_position = :none, grid_line_width = 0mm, panel_stroke = colorant"black",
        major_label_font_size = fontsize));

    if format == "png"
        draw(PNG(outfile, 5inch, 5inch, dpi = dpi), plt1);
    elseif format == "pdf"
        draw(PDF(outfile, 5inch, 5inch, dpi = dpi), plt1);
    elseif format == "svg"
        draw(SVG(outfile, 5inch, 5inch), plt1);
    else
        throw(ArgumentError("This should not display. Please report format issue."))
    end
end

function qq(df::DataFrame; pvalvar::AbstractString = "pval", kwargs...)
    pvalvar = Symbol(pvalvar)
    if !(string(pvalvar) in names(df))
        throw(ArgumentError(string(pvalvar) * " is not in the dataframe. Please rename
        the column of pvalues to `pval` or specify the correct name using the `pvalvar` argument"))
    end
    qq(df[!, pvalvar]; kwargs...)
end

"""
    manhattan(data::DataFrame)

# Position arguments

- `data::DataFrame`: A DataFrame containing information to be used in the Manhattan Plot.
Note, DataFrame must have the following values saved under the corresponding names.
pvalues:pval, chromosome:chr. Additionally, the DataFrame must be in order of basepairs
going from first to last if there's no position arguement. Optionally, if there is
basepair information, then the position variable must be named `pos`.

    manhattan(pvalues::AbstractArray, chr::AbstractArray)

# Position arguments

- `pvalues::AbstractArray`: pvalues for the associated GWAS. Must be in the
order of the basepairs.

- `chr::AbstractArray`: An array of chromosome identifiers for each pvalue.
Must match order with pvalues.

    manhattan(pvalues::AbstractArray, chr::AbstractArray, pos::AbstractArray)

# Position arguments

- `pvalues::AbstractArray`: pvalues for the associated GWAS.
Must be in the same order of the basepairs and chromosomes.

- `chr::AbstractArray`: An array of chromosome identifiers for each pvalue.
Must match order with pvalues and positions.

- `pos::AbstractArray`: An array of basepair positions for each pvalue/chromosome.
Must match order with pvalues and chromosomes.

# Keyword arguments

- `titles::AbstractString`: Title for the plot. Default is "Manhattan Plot".
 To have blank enter "".

- `outfile::AbstractString`: output name to save the manhattan plot. Name should end in format.
Default is "manhattan.png". Supports .png, .pdf, and .svg files.

- `dpi::Int64`: dots per inch to save the png file. Higher DPI results in
larger file with higher resolution. Default dpi is 350.

- `xlabel::AbstractString`: option to replace x-label text

- `ylabel::AbstractString`: option to replace y-label text

- `ymax::Union{Float64, Int64}`: Specified maximum y value to represent
on the plot

- `ystep`: Step-size for y-axis label ticks. Default value is 2.5

- `signifline::Union{Float64, Int64}`: Line to draw significance at.
Default in Bonferonni corrected p-value for α = 0.05.

- `linecolor::AbstractString`: Color for significance line. Default
is 'deepskyblue1'.

- `fontsize` size of the axis labels. Default is "15pt".

- `pvalvar` variable indicating pvalue column name (for dataframes only). Default is "pval".

- `chrvar` variable indicating chromosome column name (for dataframes only). Default is "chr".

- `posvar` variable indicating BP position column name (for dataframes only). Default is "pos".

- `annotatevar` variable indicating annotation column name (for dataframes only). Default is "gene".

- `annotateinds` indicies of rows to include as annotation for the manhattan plot. Default is `nothing`.
"""
function manhattan(data::DataFrame;
    titles::AbstractString = "Manhattan Plot",
    outfile::AbstractString = "manhattan.png",
    dpi::Int64 = 350,
    xlabel::AbstractString = "Chromosome",
    ylabel::AbstractString = "-log10(p)",
    ymax::Union{Float64, Int64} = 0,
    ystep::Union{Float64, Int64} = 2.5,
    signifline::Union{Float64, Int64} = -1,
    linecolor::AbstractString = "deepskyblue1",
    fontsize = 15pt,
    pvalvar::Union{AbstractString, Symbol} = "pval",
    chrvar::Union{AbstractString, Symbol} = "chr",
    posvar::Union{AbstractString, Symbol} = "pos",
    annotatevar::Union{AbstractString, Symbol} = "gene",
    annotateinds::Union{Nothing, AbstractVector{<:Integer}} = nothing,
    kwargs...)

    format = lowercase(split(outfile, ".")[end])
    if !(format in ["png", "pdf", "svg"])
        throw(ArgumentError("Incorrect file format specified in outfile option.
         outfile should end in .png, .pdf, or .svg (very slow)."))
    end

    df = deepcopy(data) #so that the original dataframe isn't modified.

    pvalvar = Symbol(pvalvar)
    chrvar = Symbol(chrvar)
    posvar = Symbol(posvar)
    annotatevar = Symbol(annotatevar)

    if !(string(pvalvar) in names(df))
        throw(ArgumentError(string(pvalvar) * " is not in the dataframe. Please rename
        the column of pvalues to `pval` or specify the correct name using the `pvalvar` argument"))
    end

    if !(string(chrvar) in names(df))
        throw(ArgumentError(string(chrvar) * " is not in the dataframe. Please rename
        the column of chromsomes to `chr` or specify the correct name using the `chrvar` argument"))
    end

    if !isnothing(annotateinds)
        @assert annotatevar in propertynames(df)  "$(string(annotatevar)) is not in the dataframe. Please rename
        the column of annotations to `gene` or specify the correct name using the `annotatevar` argument"
    end

    using_basepairs = string(posvar) in names(df)
    if !using_basepairs
        println("No basepair information will be used as the basepair position variable
         was not specified properly. To use BP position info, use the `posvar` argument
         to specify the variable in the dataframe that uses basepairs.")
    end

    all(df[!, pvalvar] .>= 0.0) && all(df[!, pvalvar] .<= 1.0) ||
    @error("Values that are not between 0 and 1 detected in $pvalvar column of $data" )

    df[!, :log10pval] = -log10.(df[!, pvalvar])
    df[!, :log10pval] = clamp.(df[!, :log10pval], 0.0, -log10(floatmin(0.0)))
    if ymax == 0
        ymax = ceil(maximum(df[!, :log10pval])) + 2.5
    end
    yticks = collect(0:ystep:ymax)

    if typeof(df[!, chrvar][1]) != String
        df[!, chrvar] = string.(df[!, chrvar])
        xlabs = unique(df[!, chrvar])[1:2:end]
    else
        xlabs = unique(df[!, chrvar])[1:2:end]
    end

    if signifline == -1
        signifline = -log10(0.05 / length(df[!, chrvar]))
    end

    if using_basepairs
        a = combine(groupby(df, chrvar),  posvar => maximum => :max_pos)
        a[!, :bp_add] = cumsum(a[!, :max_pos]) - a[!, :max_pos]
        df = innerjoin(df, a, on  = chrvar)
        df[!, :adj_bp] = df[!, posvar] .+ df[!, :bp_add]
        xticks = combine(groupby(df, chrvar), :adj_bp => d -> (maximum(d) + minimum(d))/2)[!, 2][1:2:end]
        #check
        function convertlabsBP(i)
            return string.(xlabs)[findfirst(xticks .== i)]
        end
    else
        xticks = Vector{Float64}(undef, length(xlabs))
        df[!, :SNPnumber] = collect(1:size(df)[1])
        global counter = 1
        for i in parse.(Int64, unique(df[!, chrvar]))[1:2:end]
            xticks[counter] = mean(df[!, :SNPnumber][df[!, chrvar] .== string.(i)])
            global counter = counter + 1
        end
        function convertlabs(i)
            return df[!, chrvar][findfirst(df[!, :SNPnumber] .== round(i))]
        end
    end

    #make the manhattan plot
    if using_basepairs
        plt1 = layer(df, x = :adj_bp, y = :log10pval, color = chrvar, Geom.point,
                intercept=[signifline], slope = [0],
                Geom.abline(color = linecolor),
                Theme(panel_fill = nothing, highlight_width = 0mm, point_size = 0.5mm,
                key_position = :none, grid_line_width = 0mm, panel_stroke = colorant"black",
                major_label_font_size = fontsize));
        if !isnothing(annotateinds)
            plt_annotate = layer(df[annotateinds, :], x = :adj_bp, y = :log10pval,
                Geom.point, Geom.label, label = annotatevar,
                color = [colorant"black"], Theme(discrete_highlight_color = identity,
                point_size = 0.6mm, point_shapes=[Shape.circle], alphas=[0.0]))
        end
        if length(unique(df[!, chrvar])) == 22
            if isnothing(annotateinds)
                mhplot = plot(plt1,
                    Guide.xticks(ticks = xticks), Guide.xlabel(xlabel),
                    Scale.x_continuous(labels = convertlabsBP), Guide.ylabel(ylabel),
                    Guide.yticks(ticks = yticks), Guide.title(titles),
                    Scale.color_discrete_manual("#d54359", "#5ab543",
                    "#a162dc", "#a9b245", "#5861cf", "#d89c39", "#6e8be3",
                    "#c7522a", "#5ea1d5", "#dd418d", "#66b974", "#9c41a3",
                    "#5a7936", "#df6cce", "#48b5a3", "#aa4584", "#997134",
                    "#655fa0", "#dd8a69", "#bf8cd3", "#a24b5e", "#df83a3"),
                    Theme(key_position = :none,
                    grid_line_width = 0mm,
                    panel_stroke = colorant"black",
                    major_label_font_size = fontsize,
                    panel_fill = nothing); kwargs...);
            else
                mhplot = plot(plt_annotate, plt1,
                    Guide.xticks(ticks = xticks), Guide.xlabel(xlabel),
                    Scale.x_continuous(labels = convertlabsBP), Guide.ylabel(ylabel),
                    Guide.yticks(ticks = yticks), Guide.title(titles),
                    Scale.color_discrete_manual("#d54359", "#5ab543",
                    "#a162dc", "#a9b245", "#5861cf", "#d89c39", "#6e8be3",
                    "#c7522a", "#5ea1d5", "#dd418d", "#66b974", "#9c41a3",
                    "#5a7936", "#df6cce", "#48b5a3", "#aa4584", "#997134",
                    "#655fa0", "#dd8a69", "#bf8cd3", "#a24b5e", "#df83a3"),
                    Theme(key_position = :none,
                    grid_line_width = 0mm,
                    panel_stroke = colorant"black",
                    major_label_font_size = fontsize,
                    panel_fill = nothing); kwargs...);
            end
        else #more than 22 chr, use default coloring
            if isnothing(annotateinds)
                mhplot = plot(plt1,
                    Guide.xticks(ticks = xticks), Guide.xlabel(xlabel), Guide.title(titles),
                    Scale.x_continuous(labels = convertlabsBP), Guide.ylabel(ylabel),
                    Guide.yticks(ticks = yticks), Scale.color_discrete,
                    Theme(key_position = :none,
                    grid_line_width = 0mm,
                    panel_stroke = colorant"black",
                    major_label_font_size = fontsize,
                    panel_fill = nothing); kwargs...);
            else
                mhplot = plot(plt_annotate, plt1,
                    Guide.xticks(ticks = xticks), Guide.xlabel(xlabel), Guide.title(titles),
                    Scale.x_continuous(labels = convertlabsBP), Guide.ylabel(ylabel),
                    Guide.yticks(ticks = yticks), Scale.color_discrete,
                    Theme(key_position = :none,
                    grid_line_width = 0mm,
                    panel_stroke = colorant"black",
                    major_label_font_size = fontsize,
                    panel_fill = nothing); kwargs...);
            end
        end
    else
        plt1 = layer(df, x = :SNPnumber, y = :log10pval, color = chrvar, Geom.point,
                intercept=[signifline], slope = [0],
                Geom.abline(color = linecolor),
                Theme(panel_fill = nothing, highlight_width = 0mm, point_size = 0.5mm,
                key_position = :none, grid_line_width = 0mm, panel_stroke = colorant"black",
                major_label_font_size = fontsize));
        if !isnothing(annotateinds)
            plt_annotate = layer(df[annotateinds, :], x = :SNPnumber, y = :log10pval,
                Geom.point, Geom.label, label = annotatevar,
                color = [colorant"black"], Theme(point_size = 0.85mm,
                highlight_width = 0mm))
        end
        if length(unique(df[!, chrvar])) == 22
            if isnothing(annotateinds)
                mhplot = plot(plt1,
                    Guide.xticks(ticks = xticks), Guide.xlabel(xlabel),
                    Scale.x_continuous(labels = convertlabs), Guide.ylabel(ylabel),
                    Guide.yticks(ticks = yticks), Guide.title(titles),
                    Scale.color_discrete_manual("#d54359", "#5ab543",
                    "#a162dc", "#a9b245", "#5861cf", "#d89c39", "#6e8be3",
                    "#c7522a", "#5ea1d5", "#dd418d", "#66b974", "#9c41a3",
                    "#5a7936", "#df6cce", "#48b5a3", "#aa4584", "#997134",
                    "#655fa0", "#dd8a69", "#bf8cd3", "#a24b5e", "#df83a3"),
                    Theme(key_position = :none,
                    grid_line_width = 0mm,
                    panel_stroke = colorant"black",
                    major_label_font_size = fontsize,
                    panel_fill = nothing); kwargs...);
            else
                mhplot = plot(plt_annotate, plt1,
                    Guide.xticks(ticks = xticks), Guide.xlabel(xlabel),
                    Scale.x_continuous(labels = convertlabs), Guide.ylabel(ylabel),
                    Guide.yticks(ticks = yticks), Guide.title(titles),
                    Scale.color_discrete_manual("#d54359", "#5ab543",
                    "#a162dc", "#a9b245", "#5861cf", "#d89c39", "#6e8be3",
                    "#c7522a", "#5ea1d5", "#dd418d", "#66b974", "#9c41a3",
                    "#5a7936", "#df6cce", "#48b5a3", "#aa4584", "#997134",
                    "#655fa0", "#dd8a69", "#bf8cd3", "#a24b5e", "#df83a3"),
                    Theme(key_position = :none,
                    grid_line_width = 0mm,
                    panel_stroke = colorant"black",
                    major_label_font_size = fontsize,
                    panel_fill = nothing); kwargs...);
            end
        else #more than 22 chr, use default coloring
            if isnothing(annotateinds)
                mhplot = plot(plt1,
                    Guide.xticks(ticks = xticks), Guide.xlabel(xlabel), Guide.title(titles),
                    Scale.x_continuous(labels = convertlabs), Guide.ylabel(ylabel),
                    Guide.yticks(ticks = yticks), Scale.color_discrete,
                    Theme(key_position = :none,
                    grid_line_width = 0mm,
                    panel_stroke = colorant"black",
                    major_label_font_size = fontsize,
                    panel_fill = nothing); kwargs...);
            else
                mhplot = plot(plt_annotate, plt1,
                    Guide.xticks(ticks = xticks), Guide.xlabel(xlabel), Guide.title(titles),
                    Scale.x_continuous(labels = convertlabs), Guide.ylabel(ylabel),
                    Guide.yticks(ticks = yticks), Scale.color_discrete,
                    Theme(key_position = :none,
                    grid_line_width = 0mm,
                    panel_stroke = colorant"black",
                    major_label_font_size = fontsize,
                    panel_fill = nothing); kwargs...);
            end
        end
    end
    if format == "png"
        draw(PNG(outfile, 6inch, 4inch, dpi = dpi), mhplot);
    elseif format == "pdf"
        draw(PDF(outfile, 6inch, 4inch, dpi = dpi), mhplot);
    elseif format == "svg"
        draw(SVG(outfile, 6inch, 4inch), mhplot);
    else
        throw(ArgumentError("This should not display. Please report format issue."))
    end
end

function manhattan(pvalues::AbstractArray, chr::AbstractArray; kwargs...)
    df = DataFrame(pval = pvalues, chr = chr)
    manhattan(df; kwargs...)
end

function manhattan(pvalues::AbstractArray, chr::AbstractArray, pos::AbstractArray; kwargs...)
    df = DataFrame(pval = pvalues, chr = chr, pos = pos)
    manhattan(df; kwargs...)
end


text:
#=
You are given an array of non-negative integers that represents a two-dimensional elevation map where each element is unit-width wall and the integer is the height. Suppose it will rain and all spots between two walls get filled up.

Compute how many units of water remain trapped on the map in O(N) time and O(1) space.

For example, given the input [2, 1, 2], we can hold 1 unit of water in the middle.

Given the input [3, 0, 1, 3, 0, 5], we can hold 3 units in the first index, 2 in the second, and 3 in the fourth index (we cannot hold 5 since it would run off to the left), so we can trap 8 units of water.
=#

function watertrapped(elevation_map::Array{Int64, 1})
    sum_water = 0
    for map_unit in 2:length(elevation_map)-1
        max_left = maximum(elevation_map[1:map_unit-1])
        max_right = maximum(elevation_map[map_unit+1:end])
        water_trapped = maximum([minimum([max_left, max_right]) - elevation_map[map_unit], 0])
        sum_water = sum_water + water_trapped
    end
    return sum_water
end

ext: jl

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
const ObjF = Union{MOI.SingleVariable, MOI.ScalarAffineFunction{Float64}}

MOI.canget(::MosekModel, ::MOI.ObjectiveFunction{<:ObjF}) = true

# TODO get with SingleVariable
function MathOptInterface.get(m::MosekModel, ::MathOptInterface.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}})
    vid = allocatedlist(m.x_block)
    subj = getindexes(m.x_block, vid)
    coeffs = getclist(m.task, subj)
    constant = getcfix(m.task)
    MOI.ScalarAffineFunction(MOI.VariableIndex.(vid), coeffs, constant)
end

MOI.supports(::MosekModel,::MOI.ObjectiveFunction{<:ObjF})  = true
MOI.supports(::MosekModel,::MOI.ObjectiveSense) = true

function MOI.set!(m::MosekModel, ::MOI.ObjectiveFunction, func::MOI.SingleVariable)
    numvar = getnumvar(m.task)
    c = zeros(Float64,numvar)
    vid = ref2id(func.variable)
    subj = getindexes(m.x_block,vid)

    c[subj[1]] = 1.0

    putclist(m.task,Int32[1:numvar...],c)
    putcfix(m.task,0.0)
end

function MOI.set!(m::MosekModel, ::MOI.ObjectiveFunction, func::MOI.ScalarAffineFunction{Float64})
    numvar = getnumvar(m.task)
    c = zeros(Float64,numvar)
    subj = getindexes(m.x_block, ref2id.(map(t -> t.variable_index, func.terms)))
    for i in 1:length(subj)
        c[subj[i]] += func.terms[i].coefficient
    end

    putclist(m.task,Int32[1:numvar...],c)
    putcfix(m.task,func.constant)
end

function MOI.modify!(m::MosekModel, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}, change :: MOI.ScalarConstantChange)
    putcfix(m.task,change.new_constant)
end

function MOI.modify!(m::MosekModel, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}, change :: MOI.ScalarCoefficientChange)
    vid = ref2id(change.variable)
    subj = getindexes(m.x_block,vid)
    putcj(m.task,Int32(subj[1]),change.new_coefficient)
end


text:
"""
    BatchedTridiagonalSolver

A batched solver for large numbers of triadiagonal systems.
"""
struct BatchedTridiagonalSolver{A, B, C, F, T, G, P}
       a :: A
       b :: B
       c :: C
       f :: F
       t :: T
    grid :: G
  params :: P
end

"""
    BatchedTridiagonalSolver(; dl, d, du, f, grid, params=nothing)

Construct a solver for batched tridiagonal systems of the form

                               b(i, j, 1)ϕ(i, j, 1) + c(i, j,   2)ϕ(i, j,   2) = f(i, j, 1),  k = 1
    a(i, j, k-1)ϕ(i, j, k-1) + b(i, j, k)ϕ(i, j, k) + c(i, j, k+1)ϕ(i, j, k+1) = f(i, j, k),  k = 2, ..., N-1
    a(i, j, N-1)ϕ(i, j, N-1) + b(i, j, N)ϕ(i, j, N)                            = f(i, j, N),  k = N

where `dl` stores the lower diagonal coefficients `a(i, j, k)`, `d` stores the diagonal coefficients `b(i, j, k)`,
`du` stores the upper diagonal coefficients `c(i, j, k)`, and `f` stores the right-hand-side terms `f(i, j, k)`. A
`grid` must be passed in.

`dl`, `d`, `du`, and `f` can be specified in three ways to describe different batched systems:
1. A 1D array indicates that the coefficients only depend on `k` and are the same for all the tridiagonal systems, i.e.
   `a(i, j, k) = a(k)`.
2. A 3D array indicates that the coefficients are different for each tridiagonal systems and depend on `(i, j, k)`.
3. A function with the signature `b(i, j, k, grid, params)` that returns the coefficient `b(i, j, k)`.

`params` is an optional named tuple of parameters that is accessible to functions.
"""
function BatchedTridiagonalSolver(arch=CPU(), FT=Float64; dl, d, du, f, grid, params=nothing)
    ArrayType = array_type(arch)
    t = zeros(FT, grid.Nx, grid.Ny, grid.Nz) |> ArrayType

    return BatchedTridiagonalSolver(dl, d, du, f, t, grid, params)
end

@inline get_coefficient(a::AbstractArray{T, 1}, i, j, k, grid, p) where {T} = @inbounds a[k]
@inline get_coefficient(a::AbstractArray{T, 3}, i, j, k, grid, p) where {T} = @inbounds a[i, j, k]
@inline get_coefficient(a::Function, i, j, k, grid, p) = a(i, j, k, grid, p)

"""
    solve_batched_tridiagonal_system!(ϕ, arch, solver)

Solve the batched tridiagonal system of linear equations described by the
`BatchedTridiagonalSolver` `solver` using a modified TriDiagonal Matrix Algorithm (TDMA)
that is still capable of solving singular systems with a zero eigenvalue.

The result is stored in `ϕ` which must have size `(grid.Nx, grid.Ny, grid.Nz)`.

Reference implementation per Numerical Recipes, Press et. al 1992 (§ 2.4).
"""
function solve_batched_tridiagonal_system!(ϕ, arch, solver)
    a, b, c, f, t, grid, params = solver.a, solver.b, solver.c, solver.f, solver.t, solver.grid, solver.params

    event = launch!(arch, grid, :xy,
                    solve_batched_tridiagonal_system_kernel!, ϕ, a, b, c, f, t, grid, params,
                    dependencies=Event(device(arch)))

    wait(device(arch), event)

    return nothing
end

@kernel function solve_batched_tridiagonal_system_kernel!(ϕ, a, b, c, f, t, grid, p)
    Nx, Ny, Nz = grid.Nx, grid.Ny, grid.Nz

    i, j = @index(Global, NTuple)

    @inbounds begin
        β  = get_coefficient(b, i, j, 1, grid, p)
        f₁ = get_coefficient(f, i, j, 1, grid, p)
        ϕ[i, j, 1] = f₁ / β

        @unroll for k = 2:Nz
            cₖ₋₁ = get_coefficient(c, i, j, k-1, grid, p)
            bₖ   = get_coefficient(b, i, j, k,   grid, p)
            aₖ₋₁ = get_coefficient(a, i, j, k-1, grid, p)

            t[i, j, k] = cₖ₋₁ / β
            β    = bₖ - aₖ₋₁ * t[i, j, k]

            # This should only happen on last element of forward pass for problems
            # with zero eigenvalue. In that case the algorithmn is still stable.
            abs(β) < 1e-12 && break

            fₖ = get_coefficient(f, i, j, k, grid, p)

            ϕ[i, j, k] = (fₖ - aₖ₋₁ * ϕ[i, j, k-1]) / β
        end

        @unroll for k = Nz-1:-1:1
            ϕ[i, j, k] = ϕ[i, j, k] - t[i, j, k+1] * ϕ[i, j, k+1]
        end
    end
end


text:
module GeometricalPredicates

# Fast, robust 2D and 3D geometrical predicates on generic point types.
# Implementation follows algorithms described in http://arxiv.org/abs/0901.4107
# and used (for e.g.) in the Illustris Simulation
# http://www.illustris-project.org/

# Author: Ariel Keselman (skariel@gmail.com)
# License: MIT
# Bug reports welcome!

export
    min_coord, max_coord,

    AbstractPoint,
        AbstractPoint2D,
        AbstractPoint3D,

    AbstractPositivelyOrientedPrimitive,
        AbstractPositivelyOrientedTriangle,
        AbstractPositivelyOrientedTetrahedron,
    AbstractNegativelyOrientedPrimitive,
        AbstractNegativelyOrientedTriangle,
        AbstractNegativelyOrientedTetrahedron,
    AbstractTriangleUnOriented,
    AbstractTetrahedronUnOriented,

    TriangleTypes, TetrahedronTypes,

    positivelyoriented, negativelyoriented, unoriented, orientation,

    Point, Point2D, Point3D, Line, Polygon, getx, gety, getz, geta, getb, getc, getd,
    seta, setb, setc, setd, setabc, setabcd,
    setab, setbc, setcd, setac, setad, setbd,
    setabd, setacd, setbcd,
    getlines, getpoints,

    Line2D, Polygon2D, Primitive, Triangle, Tetrahedron,

    length2, area, volume, centroid, circumcenter, circumradius2, incircle, intriangle, inpolygon,

    peanokey, hilbertsort!, mssort!, clean!

const _float_err = eps(Float64)
const _abs_err_incircle_2d = 12*_float_err
const _abs_err_incircle_3d = 48*_float_err
const _abs_err_orientation_2d = 2*_float_err
const _abs_err_orientation_3d = 6*_float_err
const _abs_err_intriangle = 6*_float_err
const _abs_err_intriangle_zero = 2*_float_err
const _abs_err_intetra = 24*_float_err
const _abs_err_intetra_zero = 6*_float_err

const min_coord = 1.0
const max_coord = 2.0 - eps(Float64)

abstract type AbstractPoint end
abstract type AbstractPoint2D <: AbstractPoint end
abstract type AbstractPoint3D <: AbstractPoint end

abstract type AbstractLine2D end
abstract type AbstractPolygon2D end

abstract type AbstractPrimitive end
abstract type AbstractUnOrientedPrimitive <: AbstractPrimitive end
abstract type AbstractOrientedPrimitive <: AbstractPrimitive end
abstract type AbstractPositivelyOrientedPrimitive <: AbstractOrientedPrimitive end
abstract type AbstractNegativelyOrientedPrimitive <: AbstractOrientedPrimitive end

abstract type AbstractTriangleUnOriented <: AbstractUnOrientedPrimitive end
abstract type AbstractTetrahedronUnOriented <: AbstractUnOrientedPrimitive end
abstract type AbstractPositivelyOrientedTriangle <: AbstractPositivelyOrientedPrimitive end
abstract type AbstractNegativelyOrientedTriangle <: AbstractNegativelyOrientedPrimitive end
abstract type AbstractPositivelyOrientedTetrahedron <: AbstractPositivelyOrientedPrimitive end
abstract type AbstractNegativelyOrientedTetrahedron <: AbstractNegativelyOrientedPrimitive end

const TriangleTypes = Union{AbstractTriangleUnOriented, AbstractPositivelyOrientedTriangle, AbstractNegativelyOrientedTriangle}
const TetrahedronTypes = Union{AbstractTetrahedronUnOriented, AbstractPositivelyOrientedTetrahedron, AbstractNegativelyOrientedTetrahedron}

# standard 2D point
struct Point2D <: AbstractPoint2D
    _x::Float64
    _y::Float64
end
Point2D() = Point2D(0., 0.)

getx(p::Point2D) = p._x
gety(p::Point2D) = p._y

# standard 3D point
struct Point3D <: AbstractPoint3D
    _x::Float64
    _y::Float64
    _z::Float64
end
Point3D() = Point3D(0., 0., 0.)

getx(p::Point3D) = p._x
gety(p::Point3D) = p._y
getz(p::Point3D) = p._z

Point(x::Real, y::Real) = Point2D(x, y)
Point(x::Real, y::Real, z::Real) = Point3D(x, y, z)

struct Line2D{T<:AbstractPoint2D} <: AbstractLine2D
    _a::T
    _b::T
    _bx::Float64
    _by::Float64
end

function Line2D(a::T, b::T) where {T<:AbstractPoint2D}
    bx = getx(b) - getx(a)
    by = gety(b) - gety(a)
    Line2D(a, b, bx, by)
end

Line(a::T, b::T) where {T<:AbstractPoint2D} = Line2D(a, b)

geta(l::Line2D) = l._a
getb(l::Line2D) = l._b

length2(l::Line2D) = l._bx*l._bx + l._by*l._by

# fine filtered orientation
function _sz_orientation(l::Line2D, p::AbstractPoint2D)
    cx = getx(p) - getx(geta(l))
    cy = gety(p) - gety(geta(l))
    _pr2 = -l._bx*cy + l._by*cx
    sz = abs(cx) + abs(cy) + abs(l._bx) + abs(l._by)
    if _pr2 < -_abs_err_orientation_2d*sz
        1
    elseif _pr2 > _abs_err_orientation_2d*sz
        -1
    else
        _exact_sign_orientation_determinant!(
            _extract_bigint(getx(geta(l))), _extract_bigint(gety(geta(l))),
            _extract_bigint(getx(getb(l))), _extract_bigint(gety(getb(l))),
            _extract_bigint(getx(p)), _extract_bigint(gety(p)))
    end
end

# gross filtered orientation, asumming maximal line size (=1.0)
function orientation(l::Line2D, p::AbstractPoint2D)
    cx = getx(p) - getx(geta(l))
    cy = gety(p) - gety(geta(l))
    _pr2 = -l._bx*cy + l._by*cx
    if _pr2 < -_abs_err_orientation_2d
        1
    elseif _pr2 > _abs_err_orientation_2d
        -1
    else
        _sz_orientation(l, p)
    end
end

"""
    struct Polygon2D <: AbstractPolygon2D

Two-dimensional polygon type.
"""
struct Polygon2D{T<:AbstractPoint2D} <: AbstractPolygon2D
    _p::Vector{T}
    _l::Vector{AbstractLine2D}
    function Polygon2D{T}(p::T...) where {T<:AbstractPoint2D}
        l = map(1:length(p)-1) do i
            Line(p[i], p[i+1])
        end
        push!(l, Line(p[end], p[1]))
        new([p...;], l)
    end
end

Polygon2D(p::T...) where {T<:AbstractPoint2D} = Polygon2D{T}(p...)

Polygon(p::T...) where {T<:AbstractPoint2D} = Polygon2D(p...)

"""
    getpoints(polygon)

Return the points of a polygon.
"""
getpoints(polygon::Polygon2D) = polygon._p

"""
    getlines(polygon)

Return the lines of a polygon.
"""
getlines(polygon::Polygon2D) = polygon._l

"""
    inpolygon(polygon, point)

Return true if `point` is inside `polygon`, which is assumed to be convex.
"""
function inpolygon(polygon::Polygon2D, point::AbstractPoint2D)
    lines = getlines(polygon)
    side = orientation(lines[1], point)
    for i = 2:length(lines)
        orientation(lines[i], point) == side || return false
    end
    true
end

macro _define_triangle_type(name, abstracttype)
    oriented = !occursin("UnOriented", string(name))
    esc(quote
        mutable struct $name{T <: AbstractPoint2D} <: $abstracttype
            _a::T; _b::T; _c::T
            _bx::Float64; _by::Float64
            _cx::Float64; _cy::Float64
            _px::Float64; _py::Float64
            _pr2::Float64
            $((oriented ? tuple() : (:(_o::Int8),))...)
            function $name{T}(a::T, b::T, c::T) where T
                t = new(a, b, c, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, $((oriented ? tuple() : (0,))...))
                clean!(t)
                t
            end
        end
    end)
end

@_define_triangle_type(UnOrientedTriangle, AbstractTriangleUnOriented)
@_define_triangle_type(PositivelyOrientedTriangle, AbstractPositivelyOrientedTriangle)
@_define_triangle_type(NegativelyOrientedTriangle, AbstractNegativelyOrientedTriangle)

abstract type AbstractOrientation end
mutable struct PositivelyOriented <: AbstractOrientation; end
mutable struct NegativelyOriented <: AbstractOrientation; end
mutable struct UnOriented <: AbstractOrientation; end

const positivelyoriented = PositivelyOriented()
const negativelyoriented = NegativelyOriented()
const unoriented = UnOriented()

Triangle(a::T, b::T, c::T, ::PositivelyOriented) where {T<:AbstractPoint2D} = PositivelyOrientedTriangle{T}(a, b, c)
Triangle(a::T, b::T, c::T, ::NegativelyOriented) where {T<:AbstractPoint2D} = NegativelyOrientedTriangle{T}(a, b, c)
Triangle(a::T, b::T, c::T, ::UnOriented) where {T<:AbstractPoint2D} = UnOrientedTriangle{T}(a, b, c)
Triangle(a::T, b::T, c::T) where {T<:AbstractPoint2D} = Triangle(a, b, c, unoriented)
Triangle(ax::Float64, ay::Float64, bx::Float64, by::Float64, cx::Float64, cy::Float64, orientation::AbstractOrientation=unoriented) =
    Triangle(Point2D(ax, ay), Point2D(bx, by), Point2D(cx, cy), orientation)
Primitive(ax::Float64, ay::Float64, bx::Float64, by::Float64, cx::Float64, cy::Float64, orientation::AbstractOrientation=unoriented) =
    Triangle(Point2D(ax, ay), Point2D(bx, by), Point2D(cx, cy), orientation)
Primitive(a::T, b::T, c::T, orientation::AbstractOrientation=unoriented) where {T<:AbstractPoint2D} =
    Triangle(a, b, c, orientation)

area(tr::TriangleTypes) = abs(tr._pr2)/2

centroid(tr::TriangleTypes) =
    Point2D(
        (getx(geta(tr)) + getx(getb(tr)) + getx(getc(tr))) / 3.0,
        (gety(geta(tr)) + gety(getb(tr)) + gety(getc(tr))) / 3.0)

function circumcenter(tr::TriangleTypes)
    d = -2.0 * tr._pr2
    Point2D(
        tr._px/d + getx(geta(tr)),
        tr._py/d + gety(geta(tr))
        )
end

function circumradius2(tr::TriangleTypes)
    c = circumcenter(tr)
    x = getx(c) - getx(geta(tr))
    y = gety(c) - gety(geta(tr))
    x*x + y*y
end

macro _define_tetrahedron_type(name, abstracttype)
    oriented = !occursin("UnOriented", string(name))
    esc(quote
        mutable struct $name{T<:AbstractPoint3D} <: $abstracttype
            _a::T; _b::T; _c::T; _d::T
            _bx::Float64; _by::Float64; _bz::Float64
            _cx::Float64; _cy::Float64; _cz::Float64
            _dx::Float64; _dy::Float64; _dz::Float64
            _px::Float64; _py::Float64; _pz::Float64
            _pr2::Float64
            $((oriented ? tuple() : (:(_o::Int8),))...)
            function $name{T}(a::AbstractPoint3D, b::AbstractPoint3D, c::AbstractPoint3D, d::AbstractPoint3D) where T
                t = new(a, b, c, d, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, $((oriented ? tuple() : (0,))...))
                clean!(t)
                t
            end
        end
    end)
end

@_define_tetrahedron_type(UnOrientedTetrahedron, AbstractTetrahedronUnOriented)
@_define_tetrahedron_type(PositivelyOrientedTetrahedron, AbstractPositivelyOrientedTetrahedron)
@_define_tetrahedron_type(NegativelyOrientedTetrahedron, AbstractNegativelyOrientedTetrahedron)

Tetrahedron(a::T, b::T, c::T, d::T, ::PositivelyOriented) where {T<:AbstractPoint3D} = PositivelyOrientedTetrahedron{T}(a, b, c, d)
Tetrahedron(a::T, b::T, c::T, d::T, ::NegativelyOriented) where {T<:AbstractPoint3D} = NegativelyOrientedTetrahedron{T}(a, b, c, d)
Tetrahedron(a::T, b::T, c::T, d::T, ::UnOriented) where {T<:AbstractPoint3D} = UnOrientedTetrahedron{T}(a, b, c, d)
Tetrahedron(a::T, b::T, c::T, d::T) where {T<:AbstractPoint3D} = Tetrahedron(a, b, c, d, unoriented)
Tetrahedron(ax::Float64, ay::Float64, az::Float64, bx::Float64, by::Float64, bz::Float64,
    cx::Float64, cy::Float64, cz::Float64, dx::Float64, dy::Float64, dz::Float64, orientation::AbstractOrientation=unoriented) =
        Tetrahedron(Point3D(ax,ay,az), Point3D(bx,by,bz), Point3D(cx,cy,cz), Point3D(dx,dy,dz), orientation)
Primitive(ax::Float64, ay::Float64, az::Float64, bx::Float64, by::Float64, bz::Float64,
    cx::Float64, cy::Float64, cz::Float64, dx::Float64, dy::Float64, dz::Float64, orientation::AbstractOrientation=unoriented) =
        Tetrahedron(Point3D(ax,ay,az), Point3D(bx,by,bz), Point3D(cx,cy,cz), Point3D(dx,dy,dz), orientation)
Primitive(a::T, b::T, c::T, d::T, orientation::AbstractOrientation=unoriented) where {T<:AbstractPoint3D} =
        Tetrahedron(a, b, c, d, orientation)

volume(tr::TetrahedronTypes) = abs(tr._pr2)/2

centroid(tr::TetrahedronTypes) =
    Point3D(
        (getx(geta(tr)) + getx(getb(tr)) + getx(getc(tr)) + getx(getd(tr))) / 4.0,
        (gety(geta(tr)) + gety(getb(tr)) + gety(getc(tr)) + gety(getd(tr))) / 4.0,
        (getz(geta(tr)) + getz(getb(tr)) + getz(getc(tr)) + getz(getd(tr))) / 4.0)

function circumcenter(tr::TetrahedronTypes)
    d = -2.0 * tr._pr2
    Point3D(
        tr._px/d + getx(geta(tr)),
        tr._py/d + gety(geta(tr)),
        tr._pz/d + getz(geta(tr))
        )
end

function circumradius2(tr::TetrahedronTypes)
    c = circumcenter(tr)
    x = getx(c) - getx(geta(tr))
    y = gety(c) - gety(geta(tr))
    z = getz(c) - getz(geta(tr))
    x*x + y*y + z*z
end

# extract exact integer representation of float to be used in exact calculations when needed
_extract_int(n::Float64) = reinterpret(UInt64, n) & 0x000fffffffffffff
_extract_bigint(n::Float64) = BigInt(_extract_int(n))

# functions to re-validate cached pre-calculations in primitives
function _clean!(t::TriangleTypes)
    t._bx = getx(getb(t))-getx(geta(t)); t._by = gety(getb(t))-gety(geta(t));
    t._cx = getx(getc(t))-getx(geta(t)); t._cy = gety(getc(t))-gety(geta(t));

    br2 = t._bx*t._bx+t._by*t._by
    cr2 = t._cx*t._cx+t._cy*t._cy

    t._px  =  br2*t._cy - t._by*cr2
    t._py  = -br2*t._cx + t._bx*cr2
    t._pr2 = -t._bx *t._cy + t._by*t._cx
end

function _sz_clean!(t::AbstractTriangleUnOriented)
    sz = abs(t._bx)+abs(t._by)+abs(t._cx)+abs(t._cy)
    if t._pr2 < -_abs_err_orientation_2d*sz
        t._o = 1
    elseif t._pr2 > _abs_err_orientation_2d*sz
        t._o = -1
    else
        t._o = _exact_sign_orientation_determinant!(
            _extract_bigint(getx(geta(t))), _extract_bigint(gety(geta(t))),
            _extract_bigint(getx(getb(t))), _extract_bigint(gety(getb(t))),
            _extract_bigint(getx(getc(t))), _extract_bigint(gety(getc(t))))
    end
end
function clean!(t::AbstractTriangleUnOriented)
    _clean!(t)
    if t._pr2 < -_abs_err_orientation_2d
        t._o = 1
    elseif t._pr2 > _abs_err_orientation_2d
        t._o = -1
    else
        _sz_clean!(t)
    end
end

function _clean!(t::TetrahedronTypes)
    t._bx = getx(getb(t))-getx(geta(t)); t._by = gety(getb(t))-gety(geta(t)); t._bz = getz(getb(t))-getz(geta(t))
    t._cx = getx(getc(t))-getx(geta(t)); t._cy = gety(getc(t))-gety(geta(t)); t._cz = getz(getc(t))-getz(geta(t))
    t._dx = getx(getd(t))-getx(geta(t)); t._dy = gety(getd(t))-gety(geta(t)); t._dz = getz(getd(t))-getz(geta(t))

    br2 = t._bx*t._bx+t._by*t._by+t._bz*t._bz
    cr2 = t._cx*t._cx+t._cy*t._cy+t._cz*t._cz
    dr2 = t._dx*t._dx+t._dy*t._dy+t._dz*t._dz

    t._px  =  t._by*t._cz*dr2 - br2*t._cz*t._dy - t._by*cr2*t._dz - t._bz*t._cy*dr2 + br2*t._cy*t._dz + t._bz*cr2*t._dy
    t._py  =  br2*t._cz*t._dx + t._bz*t._cx*dr2 - br2*t._cx*t._dz - t._bx*t._cz*dr2 - t._bz*cr2*t._dx + t._bx*cr2*t._dz
    t._pz  =  br2*t._cx*t._dy + t._bx*t._cy*dr2 + t._by*cr2*t._dx - br2*t._cy*t._dx - t._bx*cr2*t._dy - t._by*t._cx*dr2
    t._pr2 = -t._bx*t._cy*t._dz  + t._bx*t._cz*t._dy  + t._by*t._cx*t._dz  - t._by*t._cz*t._dx  - t._bz*t._cx*t._dy  + t._bz*t._cy*t._dx
end

function _sz_clean!(t::AbstractTetrahedronUnOriented)
    sz = abs(t._bx) + abs(t._by) + abs(t._bz) +
                    abs(t._cx) + abs(t._cy) + abs(t._cz) +
                    abs(t._dx) + abs(t._dy) + abs(t._dz)

    # calculate the orientation
    if t._pr2 < -_abs_err_orientation_3d*sz
        t._o = 1
    elseif t._pr2 > _abs_err_orientation_3d*sz
        t._o = -1
    else
        # exact calculation is required
        t._o = _exact_sign_orientation_determinant!(
            _extract_bigint(getx(geta(t))), _extract_bigint(gety(geta(t))), _extract_bigint(getz(geta(t))),
            _extract_bigint(getx(getb(t))), _extract_bigint(gety(getb(t))), _extract_bigint(getz(getb(t))),
            _extract_bigint(getx(getc(t))), _extract_bigint(gety(getc(t))), _extract_bigint(getz(getc(t))),
            _extract_bigint(getx(getd(t))), _extract_bigint(gety(getd(t))), _extract_bigint(getz(getd(t))))
    end
end

function clean!(t::AbstractTetrahedronUnOriented)
    _clean!(t)
    # calculate the orientation
    if t._pr2 < -_abs_err_orientation_3d
        t._o = 1
    elseif t._pr2 > _abs_err_orientation_3d
        t._o = -1
    else
        _sz_clean!(t)
    end
end

clean!(t::AbstractOrientedPrimitive) = _clean!(t)

# getting points of the primitive
geta(t::AbstractPrimitive) = t._a
getb(t::AbstractPrimitive) = t._b
getc(t::AbstractPrimitive) = t._c
getd(t::TetrahedronTypes) = t._d

# changing points in the primitive
seta(t::AbstractPrimitive, p::AbstractPoint) = (t._a=p; clean!(t))
setb(t::AbstractPrimitive, p::AbstractPoint) = (t._b=p; clean!(t))
setc(t::AbstractPrimitive, p::AbstractPoint) = (t._c=p; clean!(t))
setd(t::TetrahedronTypes, p::AbstractPoint3D) = (t._d=p; clean!(t))
setabc(t::AbstractPrimitive, pa::AbstractPoint, pb::AbstractPoint, pc::AbstractPoint) = (t._a=pa; t._b=pb; t._c=pc; clean!(t))

setab(t::AbstractPrimitive, pa::AbstractPoint, pb::AbstractPoint) = (t._a=pa; t._b=pb; clean!(t))
setbc(t::AbstractPrimitive, pb::AbstractPoint, pc::AbstractPoint) = (t._b=pb; t._c=pc; clean!(t))
setac(t::AbstractPrimitive, pa::AbstractPoint, pc::AbstractPoint) = (t._a=pa; t._c=pc; clean!(t))

setabcd(t::TetrahedronTypes, pa::AbstractPoint3D, pb::AbstractPoint3D, pc::AbstractPoint3D, pd::AbstractPoint3D) = (t._a=pa; t._b=pb; t._c=pc; t._d=pd; clean!(t))

setabc(t::TetrahedronTypes, pa::AbstractPoint3D, pb::AbstractPoint3D, pc::AbstractPoint3D) = (t._a=pa; t._b=pb; t._c=pc; clean!(t))
setabd(t::TetrahedronTypes, pa::AbstractPoint3D, pb::AbstractPoint3D, pd::AbstractPoint3D) = (t._a=pa; t._b=pb; t._d=pd; clean!(t))
setacd(t::TetrahedronTypes, pa::AbstractPoint3D, pc::AbstractPoint3D, pd::AbstractPoint3D) = (t._a=pa; t._c=pc; t._d=pd; clean!(t))
setbcd(t::TetrahedronTypes, pb::AbstractPoint3D, pc::AbstractPoint3D, pd::AbstractPoint3D) = (t._b=pb; t._c=pc; t._d=pd; clean!(t))

setad(t::TetrahedronTypes, pa::AbstractPoint3D, pd::AbstractPoint3D) = (t._a=pa; t._d=pd; clean!(t))
setbd(t::TetrahedronTypes, pb::AbstractPoint3D, pd::AbstractPoint3D) = (t._b=pb; t._d=pd; clean!(t))
setcd(t::TetrahedronTypes, pc::AbstractPoint3D, pd::AbstractPoint3D) = (t._c=pc; t._d=pd; clean!(t))

orientation(p::AbstractUnOrientedPrimitive) = p._o
orientation(p::AbstractPositivelyOrientedPrimitive) = 1
orientation(p::AbstractNegativelyOrientedPrimitive) = -1
orientation(ax::Float64, ay::Float64, bx::Float64, by::Float64, cx::Float64, cy::Float64) =
    orientation(Triangle(ax, ay, bx, by, cx, cy))
orientation(ax::Float64, ay::Float64, az::Float64, bx::Float64, by::Float64, bz::Float64, cx::Float64, cy::Float64, cz::Float64, dx::Float64, dy::Float64, dz::Float64) =
    orientation(Tetrahedron(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz))

# exact orientation for triangle
function _exact_sign_orientation_determinant!(ax::BigInt, ay::BigInt, bx::BigInt, by::BigInt, cx::BigInt, cy::BigInt)
    bx -= ax; by -= ay
    cx -= ax; cy -= ay
    Int(sign(bx*cy - by*cx))
end

# exact orientation for tetrahedron
function _exact_sign_orientation_determinant!(ax::BigInt, ay::BigInt, az::BigInt, bx::BigInt, by::BigInt, bz::BigInt, cx::BigInt, cy::BigInt, cz::BigInt, dx::BigInt, dy::BigInt, dz::BigInt)
    bx -= ax; by -= ay; bz -= az
    cx -= ax; cy -= ay; cz -= az
    dx -= ax; dy -= ay; dz -= az
    Int(sign(+bx*cy*dz - bx*cz*dy - by*cx*dz + by*cz*dx + bz*cx*dy - bz*cy*dx))
end

# exact incircle for triangle
function _exact_sign_incircle_determinant!(ax::BigInt, ay::BigInt, bx::BigInt, by::BigInt, cx::BigInt, cy::BigInt, px::BigInt, py::BigInt)
    bx -= ax; by -= ay;
    cx -= ax; cy -= ay;
    px -= ax; py -= ay;
    br2 = bx*bx+by*by
    cr2 = cx*cx+cy*cy
    pr2 = px*px+py*py
    Int(sign(-br2*cx*py + br2*cy*px + bx*cr2*py - bx*cy*pr2 - by*cr2*px + by*cx*pr2))
end

# exact incircle for tetrahedron
function _exact_sign_incircle_determinant!(ax::BigInt, ay::BigInt, az::BigInt, bx::BigInt, by::BigInt, bz::BigInt, cx::BigInt, cy::BigInt, cz::BigInt, dx::BigInt, dy::BigInt, dz::BigInt, px::BigInt, py::BigInt, pz::BigInt)
    bx -= ax; by -= ay; bz -= az
    cx -= ax; cy -= ay; cz -= az
    dx -= ax; dy -= ay; dz -= az
    px -= ax; py -= ay; pz -= az
    br2 = bx*bx+by*by+bz*bz
    cr2 = cx*cx+cy*cy+cz*cz
    dr2 = dx*dx+dy*dy+dz*dz
    pr2 = px*px+py*py+pz*pz
    Int(sign(
        +br2*cx*dy*pz - br2*cx*dz*py - br2*cy*dx*pz + br2*cy*dz*px +
         br2*cz*dx*py - br2*cz*dy*px - bx*cr2*dy*pz + bx*cr2*dz*py +
         bx*cy*dr2*pz - bx*cy*dz*pr2 - bx*cz*dr2*py + bx*cz*dy*pr2 +
         by*cr2*dx*pz - by*cr2*dz*px - by*cx*dr2*pz + by*cx*dz*pr2 +
         by*cz*dr2*px - by*cz*dx*pr2 - bz*cr2*dx*py + bz*cr2*dy*px +
         bz*cx*dr2*py - bz*cx*dy*pr2 - bz*cy*dr2*px + bz*cy*dx*pr2))
end

# finer filtered incircle
function _sz_incircle(t::TriangleTypes, p::AbstractPoint2D, px::Float64, py::Float64, pr2::Float64)
    if orientation(t) != 0
        sz = abs(px)+abs(py)+abs(t._bx)+abs(t._by)+abs(t._cx)+abs(t._cy)
        d = t._px*px + t._py*py + t._pr2*pr2
        if d < -_abs_err_incircle_2d*sz
            return -orientation(t)
        elseif d > _abs_err_incircle_2d*sz
            return orientation(t)
        end
    end

    exact_in = _exact_sign_incircle_determinant!(
        _extract_bigint(getx(geta(t))), _extract_bigint(gety(geta(t))),
        _extract_bigint(getx(getb(t))), _extract_bigint(gety(getb(t))),
        _extract_bigint(getx(getc(t))), _extract_bigint(gety(getc(t))),
        _extract_bigint(getx(p))  , _extract_bigint(gety(p)))

    if orientation(t) != 0
        return orientation(t)*exact_in
    elseif exact_in == 0
        return 1
    else
        return 2
    end
end

# gross filtered incircle, asumming maximal triangle size (=1.0)
function incircle(t::TriangleTypes, p::AbstractPoint2D)
    px  = getx(p) - getx(geta(t))
    py  = gety(p) - gety(geta(t))
    pr2 = px*px + py*py
    if orientation(t) != 0
        d = t._px*px + t._py*py + t._pr2*pr2
        if d < -_abs_err_incircle_2d
            return -orientation(t)
        elseif d > _abs_err_incircle_2d
            return orientation(t)
        end
    end
    _sz_incircle(t, p, px, py, pr2)
end

# finer filtered incircle
function _sz_incircle(t::TetrahedronTypes, p::AbstractPoint3D)
    if orientation(t) != 0
        px  = getx(p) - getx(geta(t))
        py  = gety(p) - gety(geta(t))
        pz  = getz(p) - getz(geta(t))

        sz = abs(px) + abs(py) + abs(pz) +
                        abs(t._bx) + abs(t._by) + abs(t._bz) +
                        abs(t._cx) + abs(t._cy) + abs(t._cz) +
                        abs(t._dx) + abs(t._dy) + abs(t._dz)

        pr2 = px*px + py*py + pz*pz
        d = t._px*px + t._py*py + t._pz*pz + t._pr2*pr2
        if d < -_abs_err_incircle_3d*sz
            return -orientation(t)
        elseif d > _abs_err_incircle_3d*sz
            return orientation(t)
        end
    end

    exact_in = _exact_sign_incircle_determinant!(
        _extract_bigint(getx(geta(t))), _extract_bigint(gety(geta(t))), _extract_bigint(getz(geta(t))),
        _extract_bigint(getx(getb(t))), _extract_bigint(gety(getb(t))), _extract_bigint(getz(getb(t))),
        _extract_bigint(getx(getc(t))), _extract_bigint(gety(getc(t))), _extract_bigint(getz(getc(t))),
        _extract_bigint(getx(getd(t))), _extract_bigint(gety(getd(t))), _extract_bigint(getz(getd(t))),
        _extract_bigint(getx(p))  , _extract_bigint(gety(p))  , _extract_bigint(getz(p)))

    if orientation(t) != 0
        return orientation(t)*exact_in
    elseif exact_in == 0
        return 1
    else
        return 2
    end
end

# gross filtered incircle, asumming maximal possible tetra size (=1.0)
function incircle(t::TetrahedronTypes, p::AbstractPoint3D)
    if orientation(t) != 0
        px  = getx(p) - getx(geta(t))
        py  = gety(p) - gety(geta(t))
        pz  = getz(p) - getz(geta(t))
        pr2 = px*px + py*py + pz*pz
        d = t._px*px + t._py*py + t._pz*pz + t._pr2*pr2
        if d < -_abs_err_incircle_3d
            return -orientation(t)
        elseif d > _abs_err_incircle_3d
            return orientation(t)
        end
    end
    _sz_incircle(t, p)
end

# helper methods to use incircle directly with coordinates
incircle(ax::Float64, ay::Float64, bx::Float64, by::Float64, cx::Float64, cy::Float64, dx::Float64, dy::Float64) =
    incircle(Triangle(ax, ay, bx, by, cx, cy), Point2D(dx, dy))
incircle(ax::Float64, ay::Float64, az::Float64, bx::Float64, by::Float64, bz::Float64, cx::Float64, cy::Float64, cz::Float64, dx::Float64, dy::Float64, dz::Float64, ex::Float64, ey::Float64, ez::Float64) =
    incircle(Tetrahedron(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz), Point3D(ex, ey, ez))

# exact intriangle (slow!)
function _exact_intriangle!(ax::BigInt, ay::BigInt, bx::BigInt, by::BigInt, cx::BigInt, cy::BigInt, px::BigInt, py::BigInt)
    bx -= ax; by -= ay;
    cx -= ax; cy -= ay;
    px -= ax; py -= ay;

    nb = -cx*py + cy*px
    nc = bx*py - by*px
    denom = bx*cy - by*cx

    sdenom = Int(sign(denom))
    if Int(sign(nb)) * sdenom < 0
        return -2
    end
    if Int(sign(nc)) * sdenom < 0
        return -3
    end
    l = nb+nc - denom
    sl = Int(sign(l)) * sdenom
    if sl > 0
        return -1
    end

    if Int(sign(nb)) == 0
        return 3
    end
    if Int(sign(nc)) == 0
        return 4
    end
    if sl == 0
        return 2
    end
    -sl
end

function _exact_intriangle!(ax::BigInt, ay::BigInt, az::BigInt, bx::BigInt, by::BigInt, bz::BigInt, cx::BigInt, cy::BigInt, cz::BigInt, dx::BigInt, dy::BigInt, dz::BigInt, px::BigInt, py::BigInt, pz::BigInt)
    bx -= ax; by -= ay; bz -= az
    cx -= ax; cy -= ay; cz -= az
    dx -= ax; dy -= ay; dz -= az
    px -= ax; py -= ay; pz -= az

    denom = bx*cy*dz-bx*cz*dy-by*cx*dz+by*cz*dx+bz*cx*dy-bz*cy*dx

    nb = cx*dy*pz-cx*dz*py-cy*dx*pz+cy*dz*px+cz*dx*py-cz*dy*px
    sdenom = Int(sign(denom))
    if Int(sign(nb)) * sdenom < 0
        return -2
    end

    nc = -bx*dy*pz+bx*dz*py+by*dx*pz-by*dz*px-bz*dx*py+bz*dy*px
    if Int(sign(nc)) * sdenom < 0
        return -3
    end

    nd = bx*cy*pz-bx*cz*py-by*cx*pz+by*cz*px+bz*cx*py-bz*cy*px
    if Int(sign(nd)) * sdenom < 0
        return -4
    end

    l = (nb+nc+nd - denom) * sdenom
    sl = Int(sign(l))
    if sl > 0
        return -1
    end

    if Int(sign(nb)) == 0
        return 3
    end
    if Int(sign(nc)) == 0
        return 4
    end
    if Int(sign(nd)) == 0
        return 5
    end
    if sl == 0
        return 2
    end
    -sl
end

# helper methods to use the exact intriangle directly with coordinates
_exact_intriangle(t::TriangleTypes, p::AbstractPoint2D) =
    _exact_intriangle!(
        _extract_bigint(getx(geta(t))), _extract_bigint(gety(geta(t))),
        _extract_bigint(getx(getb(t))), _extract_bigint(gety(getb(t))),
        _extract_bigint(getx(getc(t))), _extract_bigint(gety(getc(t))),
        _extract_bigint(getx(p))  , _extract_bigint(gety(p)))

_exact_intriangle(t::TetrahedronTypes, p::AbstractPoint3D) =
    _exact_intriangle!(
        _extract_bigint(getx(geta(t))), _extract_bigint(gety(geta(t))), _extract_bigint(getz(geta(t))),
        _extract_bigint(getx(getb(t))), _extract_bigint(gety(getb(t))), _extract_bigint(getz(getb(t))),
        _extract_bigint(getx(getc(t))), _extract_bigint(gety(getc(t))), _extract_bigint(getz(getc(t))),
        _extract_bigint(getx(getd(t))), _extract_bigint(gety(getd(t))), _extract_bigint(getz(getd(t))),
        _extract_bigint(getx(p)), _extract_bigint(gety(p)), _extract_bigint(getz(p)))

# finer filter for intriangle, using triangle actual size
function _sz_intriangle(t::TriangleTypes, p::AbstractPoint2D, px::Float64, py::Float64)
    sz = abs(px)+abs(py)+abs(t._bx)+abs(t._by)+abs(t._cx)+abs(t._cy)

    nb = (-t._cx*py + t._cy*px) * sign(-t._pr2)
    if nb < -_abs_err_intriangle_zero*sz
        return -2
    elseif nb < _abs_err_intriangle_zero*sz
        # we need an exact calculation
        return _exact_intriangle(t, p)
    end

    nc = (t._bx*py - t._by*px) * sign(-t._pr2)
    if nc < -_abs_err_intriangle_zero*sz
        return -3
    elseif nc < _abs_err_intriangle_zero*sz
        # we need an exact calculation
        return _exact_intriangle(t, p)
    end

    l = nb+nc + t._pr2*sign(-t._pr2)
    if l < -_abs_err_intriangle*sz
        return 1
    elseif l > _abs_err_intriangle*sz
        return -1
    end
    # we need an exact calculation
    _exact_intriangle(t, p)
end

# gross filter, assuming maximal possible triangle size (=1.0)
function intriangle(t::TriangleTypes, p::AbstractPoint2D)
    px = getx(p) - getx(geta(t)); py = gety(p) - gety(geta(t))

    nb = (-t._cx*py + t._cy*px) * sign(-t._pr2)
    if nb < -_abs_err_intriangle_zero
        return -2
    elseif nb < _abs_err_intriangle_zero
        return _sz_intriangle(t, p, px, py)
    end

    nc = (t._bx*py - t._by*px) * sign(-t._pr2)
    if nc < -_abs_err_intriangle_zero
        return -3
    elseif nc < _abs_err_intriangle_zero
        return _sz_intriangle(t, p, px, py)
    end

    l = nb+nc + t._pr2*sign(-t._pr2)
    if l < -_abs_err_intriangle
        return 1
    elseif l > _abs_err_intriangle
        return -1
    end
    _sz_intriangle(t, p, px, py)
end

# fine filtered intriangle
function _sz_intriangle(t::TetrahedronTypes, p::AbstractPoint3D)
    px = getx(p) - getx(geta(t)); py = gety(p) - gety(geta(t)); pz = getz(p) - getz(geta(t))
    sz = abs(px) + abs(py) + abs(pz) +
                    abs(t._bx) + abs(t._by) + abs(t._bz) +
                    abs(t._cx) + abs(t._cy) + abs(t._cz) +
                    abs(t._dx) + abs(t._dy) + abs(t._dz)

    nb = (t._cx*t._dy*pz-t._cx*t._dz*py-t._cy*t._dx*pz+t._cy*t._dz*px+t._cz*t._dx*py-t._cz*t._dy*px) * sign(-t._pr2)
    if nb < -_abs_err_intetra_zero*sz
        return -2
    elseif nb < _abs_err_intetra_zero*sz
        # we need an exact calculation
        return _exact_intriangle(t, p)
    end

    nc = (-t._bx*t._dy*pz+t._bx*t._dz*py+t._by*t._dx*pz-t._by*t._dz*px-t._bz*t._dx*py+t._bz*t._dy*px) * sign(-t._pr2)
    if nc < -_abs_err_intetra_zero*sz
        return -3
    elseif nc < _abs_err_intetra_zero*sz
        # we need an exact calculation
        return _exact_intriangle(t, p)
    end

    nd = (t._bx*t._cy*pz-t._bx*t._cz*py-t._by*t._cx*pz+t._by*t._cz*px+t._bz*t._cx*py-t._bz*t._cy*px) * sign(-t._pr2)
    if nd < -_abs_err_intetra_zero*sz
        return -4
    elseif nd < _abs_err_intetra_zero*sz
        # we need an exact calculation
        return _exact_intriangle(t, p)
    end

    l = nb+nc+nd + t._pr2*sign(-t._pr2)
    if l < -_abs_err_intetra*sz
        return 1
    elseif l > _abs_err_intetra*sz
        return -1
    end
    # we need an exact calculation
    _exact_intriangle(t, p)
end

# gross filtered intriangle
function intriangle(t::TetrahedronTypes, p::AbstractPoint3D)
    px = getx(p) - getx(geta(t)); py = gety(p) - gety(geta(t)); pz = getz(p) - getz(geta(t))

    nb = (t._cx*t._dy*pz-t._cx*t._dz*py-t._cy*t._dx*pz+t._cy*t._dz*px+t._cz*t._dx*py-t._cz*t._dy*px) * sign(-t._pr2)
    if nb < -_abs_err_intetra_zero
        return -2
    elseif nb < _abs_err_intetra_zero
        return _sz_intriangle(t, p)
    end

    nc = (-t._bx*t._dy*pz+t._bx*t._dz*py+t._by*t._dx*pz-t._by*t._dz*px-t._bz*t._dx*py+t._bz*t._dy*px) * sign(-t._pr2)
    if nc < -_abs_err_intetra_zero
        return -3
    elseif nc < _abs_err_intetra_zero
        return _sz_intriangle(t, p)
    end

    nd = (t._bx*t._cy*pz-t._bx*t._cz*py-t._by*t._cx*pz+t._by*t._cz*px+t._bz*t._cx*py-t._bz*t._cy*px) * sign(-t._pr2)
    if nd < -_abs_err_intetra_zero
        return -4
    elseif nd < _abs_err_intetra_zero
        return _sz_intriangle(t, p)
    end

    l = nb+nc+nd + t._pr2*sign(-t._pr2)
    if l < -_abs_err_intetra
        return 1
    elseif l > _abs_err_intetra
        return -1
    end
    _sz_intriangle(t, p)
end

# helper  methods to use the filtered (fast, exact) intriangle directly with raw coordinates
intriangle(ax::Float64, ay::Float64, bx::Float64, by::Float64, cx::Float64, cy::Float64, px::Float64, py::Float64) =
    intriangle(Triangle(ax, ay, bx, by, cx, cy), Point2D(px, py))
intriangle(ax::Float64, ay::Float64, az::Float64, bx::Float64, by::Float64, bz::Float64, cx::Float64, cy::Float64, cz::Float64, dx::Float64, dy::Float64, dz::Float64, px::Float64, py::Float64, pz::Float64) =
    intriangle(Tetrahedron(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz), Point3D(px, py, pz))

###################################################################################
#
#    Hilbert stuff
#
#

# number of bits to use per dimension in calculating the peano-key
const peano_2D_bits = 31
const peano_3D_bits = 21

# implementing 2D scale dependent Peano-Hilbert indexing

 _extract_peano_bin_num(nbins::Int64, n::Float64) = trunc(Integer, (n-1)*nbins )

# calculate peano key for given point
function peanokey(p::AbstractPoint2D, bits::Int64=peano_2D_bits)
    n = 1 << bits
    s = n >> 1; d = 0
    x = _extract_peano_bin_num(n, getx(p))
    y = _extract_peano_bin_num(n, gety(p))
    while true
        rx = (x & s) > 0
        ry = (y & s) > 0
        d += s * s * xor(3 * rx, ry)
        s = s >> 1
        (s == 0) && break
        if ry == 0
            if rx == 1
                x = n - 1 - x;
                y = n - 1 - y;
            end
            x, y = y, x
        end
    end
    d
end

# Inverse calculation. I.e. calculate the point that given given peano key
function Point2D(peanokey::Int64, bits::Int64=peano_2D_bits)
    n = 1 << bits
    x = 0; y = 0; s=1
    while true

        rx = 1 & (peanokey >> 1)
        ry = 1 & xor(peanokey, rx)

        if ry == 0
            if rx == 1
                x = s - 1 - x;
                y = s - 1 - y;
            end
            x, y = y, x
        end

        x += s * rx
        y += s * ry

        s = s << 1
        (s >= n) && break

        peanokey = peanokey >> 2
    end
    Point2D(1+x/n, 1+y/n)
end

# implementing 3D scaleful Peano-Hilbert indexing

const quadrants_arr = [
  0, 7, 1, 6, 3, 4, 2, 5,
  7, 4, 6, 5, 0, 3, 1, 2,
  4, 3, 5, 2, 7, 0, 6, 1,
  3, 0, 2, 1, 4, 7, 5, 6,
  1, 0, 6, 7, 2, 3, 5, 4,
  0, 3, 7, 4, 1, 2, 6, 5,
  3, 2, 4, 5, 0, 1, 7, 6,
  2, 1, 5, 6, 3, 0, 4, 7,
  6, 1, 7, 0, 5, 2, 4, 3,
  1, 2, 0, 3, 6, 5, 7, 4,
  2, 5, 3, 4, 1, 6, 0, 7,
  5, 6, 4, 7, 2, 1, 3, 0,
  7, 6, 0, 1, 4, 5, 3, 2,
  6, 5, 1, 2, 7, 4, 0, 3,
  5, 4, 2, 3, 6, 7, 1, 0,
  4, 7, 3, 0, 5, 6, 2, 1,
  6, 7, 5, 4, 1, 0, 2, 3,
  7, 0, 4, 3, 6, 1, 5, 2,
  0, 1, 3, 2, 7, 6, 4, 5,
  1, 6, 2, 5, 0, 7, 3, 4,
  2, 3, 1, 0, 5, 4, 6, 7,
  3, 4, 0, 7, 2, 5, 1, 6,
  4, 5, 7, 6, 3, 2, 0, 1,
  5, 2, 6, 1, 4, 3, 7, 0]
quadrants(a::Int64, b::Int64, c::Int64, d::Int64) = (@inbounds x = quadrants_arr[1+a<<3+b<<2+c<<1+d]; x)
rotxmap_table = [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 0, 1, 2, 3, 17, 18, 19, 16, 23, 20, 21, 22]
rotymap_table = [1, 2, 3, 0, 16, 17, 18, 19, 11, 8, 9, 10, 22, 23, 20, 21, 14, 15, 12, 13, 4, 5, 6, 7]
rotx_table = [3, 0, 0, 2, 2, 0, 0, 1]
roty_table = [0, 1, 1, 2, 2, 3, 3, 0]
sense_table = [-1, -1, -1, +1, +1, -1, -1, -1]

function peanokey(p::AbstractPoint3D, bits::Int64=peano_3D_bits)
    n = 1 << bits
    x = _extract_peano_bin_num(n, getx(p))
    y = _extract_peano_bin_num(n, gety(p))
    z = _extract_peano_bin_num(n, getz(p))
    mask = 1 << (bits - 1)
    key = 0
    rotation = 0
    sense = 1
    for i in 1:bits
        bitx = (x & mask > 0) ? 1 : 0
        bity = (y & mask > 0) ? 1 : 0
        bitz = (z & mask > 0) ? 1 : 0

        quad = quadrants(rotation, bitx, bity, bitz)

        key <<= 3
        key += sense == 1 ? quad : 7-quad

        @inbounds rotx = rotx_table[quad+1]
        @inbounds roty = roty_table[quad+1]
        @inbounds sense *= sense_table[quad+1]

        while rotx > 0
            @inbounds rotation = rotxmap_table[rotation+1]
            rotx -= 1
        end

        while(roty > 0)
            @inbounds rotation = rotymap_table[rotation+1]
            roty -= 1
        end
        mask >>= 1
    end

    key
end

const quadrants_inverse_x = Array{Int64}(undef,24,8)
const quadrants_inverse_y = Array{Int64}(undef,24,8)
const quadrants_inverse_z = Array{Int64}(undef,24,8)

function _init_inv_peano_3d()
    for rotation in 0:23
        for bitx in 0:1
            for bity in 0:1
                for bitz in 0:1
                    quad = quadrants(rotation, bitx, bity, bitz)
                    quadrants_inverse_x[rotation+1, quad+1] = bitx;
                    quadrants_inverse_y[rotation+1, quad+1] = bity;
                    quadrants_inverse_z[rotation+1, quad+1] = bitz;
                end
            end
        end
    end
end

_init_inv_peano_3d()

# inverse transformation, give a key get a point
function Point3D(key::Int64, bits::Int64=peano_3D_bits)

    n = 1 << bits

    shift = 3*(bits - 1)
    mask = 7 << shift

    rotation = 0
    sense = 1

    x = 0
    y = 0
    z = 0

    for i in 1:bits
        keypart = (key & mask) >> shift

        quad = sense == 1 ? keypart : 7 - keypart

        x = (x << 1) + quadrants_inverse_x[rotation+1, quad+1]
        y = (y << 1) + quadrants_inverse_y[rotation+1, quad+1]
        z = (z << 1) + quadrants_inverse_z[rotation+1, quad+1]

        rotx = rotx_table[quad+1]
        roty = roty_table[quad+1]
        sense *= sense_table[quad+1]

        while rotx > 0
            rotation = rotxmap_table[rotation+1]
            rotx -= 1
        end

        while roty > 0
            rotation = rotymap_table[rotation+1]
            roty -= 1
        end

        mask >>= 3
        shift -= 3
    end
    Point3D(1+x/n, 1+y/n, 1+z/n)
end

# implementing scale-free Hilbert ordering. Real all about it here:
# http://doc.cgal.org/latest/Spatial_sorting/index.html

abstract type AbstractCoordinate end
mutable struct CoordinateX <: AbstractCoordinate end
mutable struct CoordinateY <: AbstractCoordinate end
mutable struct CoordinateZ <: AbstractCoordinate end
const coordinatex = CoordinateX()
const coordinatey = CoordinateY()
const coordinatez = CoordinateZ()
next2d(::CoordinateX) = coordinatey
next2d(::CoordinateY) = coordinatex
next3d(::CoordinateX) = coordinatey
next3d(::CoordinateY) = coordinatez
next3d(::CoordinateZ) = coordinatex
nextnext3d(::CoordinateX) = coordinatez
nextnext3d(::CoordinateY) = coordinatex
nextnext3d(::CoordinateZ) = coordinatey

abstract type AbstractDirection end
mutable struct Forward <: AbstractDirection end
mutable struct Backward <: AbstractDirection end
const forward = Forward()
const backward = Backward()
Base.:!(::Forward) = backward
Base.:!(::Backward) = forward

compare(::Forward, ::CoordinateX, p1::AbstractPoint, p2::AbstractPoint) = getx(p1) < getx(p2)
compare(::Backward, ::CoordinateX, p1::AbstractPoint, p2::AbstractPoint) = getx(p1) > getx(p2)
compare(::Forward, ::CoordinateY, p1::AbstractPoint, p2::AbstractPoint) = gety(p1) < gety(p2)
compare(::Backward, ::CoordinateY, p1::AbstractPoint, p2::AbstractPoint) = gety(p1) > gety(p2)
compare(::Forward, ::CoordinateZ, p1::AbstractPoint, p2::AbstractPoint) = getz(p1) < getz(p2)
compare(::Backward, ::CoordinateZ, p1::AbstractPoint, p2::AbstractPoint) = getz(p1) > getz(p2)

function select!(direction::AbstractDirection, coordinate::AbstractCoordinate, v::Array{T,1}, k::Int, lo::Int, hi::Int) where T<:AbstractPoint
    lo <= k <= hi || error("select index $k is out of range $lo:$hi")
    @inbounds while lo < hi
        if hi-lo == 1
            if compare(direction, coordinate, v[hi], v[lo])
                v[lo], v[hi] = v[hi], v[lo]
            end
            return v[k]
        end
        pivot = v[(lo+hi)>>>1]
        i, j = lo, hi
        while true
            while compare(direction, coordinate, v[i], pivot); i += 1; end
            while compare(direction, coordinate, pivot, v[j]); j -= 1; end
            i <= j || break
            v[i], v[j] = v[j], v[i]
            i += 1; j -= 1
        end
        if k <= j
            hi = j
        elseif i <= k
            lo = i
        else
            return pivot
        end
    end
    return v[lo]
end

function hilbertsort!(directionx::AbstractDirection, directiony::AbstractDirection, coordinate::AbstractCoordinate, a::Array{T,1}, lo::Int64, hi::Int64, lim::Int64=4) where T<:AbstractPoint2D
    hi-lo <= lim && return a

    i2 = (lo+hi)>>>1
    i1 = (lo+i2)>>>1
    i3 = (i2+hi)>>>1

    select!(directionx, coordinate, a, i2, lo, hi)
    select!(directiony, next2d(coordinate), a, i1, lo, i2)
    select!(!directiony, next2d(coordinate), a, i3, i2, hi)

    hilbertsort!(directiony, directionx, next2d(coordinate), a, lo, i1, lim)
    hilbertsort!(directionx, directiony, coordinate, a, i1, i2, lim)
    hilbertsort!(directionx, directiony, coordinate, a, i2, i3, lim)
    hilbertsort!(!directiony, !directionx, next2d(coordinate), a, i3, hi, lim)

    return a
end

function hilbertsort!(directionx::AbstractDirection, directiony::AbstractDirection, directionz::AbstractDirection, coordinate::AbstractCoordinate, a::Array{T,1}, lo::Int64, hi::Int64, lim::Int64=8) where T<:AbstractPoint3D
    hi-lo <= lim && return a

    i4 = (lo+hi)>>>1
    i2 = (lo+i4)>>>1
    i1 = (lo+i2)>>>1
    i3 = (i2+i4)>>>1
    i6 = (i4+hi)>>>1
    i5 = (i4+i6)>>>1
    i7 = (i6+hi)>>>1

    select!(directionx, coordinate, a, i4, lo, hi)
    select!(directiony, next3d(coordinate), a, i2, lo, i4)
    select!(directionz, nextnext3d(coordinate), a, i1, lo, i2)
    select!(!directionz, nextnext3d(coordinate), a, i3, i2, i4)
    select!(!directiony, next3d(coordinate), a, i6, i4, hi)
    select!(directionz, nextnext3d(coordinate), a, i5, i4, i6)
    select!(!directionz, nextnext3d(coordinate), a, i7, i6, hi)

    hilbertsort!( directionz,  directionx,  directiony, nextnext3d(coordinate), a, lo, i1, lim)
    hilbertsort!( directiony,  directionz,  directionx, next3d(coordinate),     a, i1, i2, lim)
    hilbertsort!( directiony,  directionz,  directionx, next3d(coordinate),     a, i2, i3, lim)
    hilbertsort!( directionx, !directiony, !directionz, coordinate,             a, i3, i4, lim)
    hilbertsort!( directionx, !directiony, !directionz, coordinate,             a, i4, i5, lim)
    hilbertsort!(!directiony,  directionz, !directionx, next3d(coordinate),     a, i5, i6, lim)
    hilbertsort!(!directiony,  directionz, !directionx, next3d(coordinate),     a, i6, i7, lim)
    hilbertsort!(!directionz, !directionx,  directiony, nextnext3d(coordinate), a, i7, hi, lim)

    return a
end

hilbertsort!(a::Array{T,1}) where {T<:AbstractPoint2D} = hilbertsort!(backward, backward, coordinatey, a, 1, length(a))
hilbertsort!(a::Array{T,1}, lo::Int64, hi::Int64, lim::Int64) where {T<:AbstractPoint2D} = hilbertsort!(backward, backward, coordinatey, a, lo, hi, lim)
hilbertsort!(a::Array{T,1}) where {T<:AbstractPoint3D} = hilbertsort!(backward, backward, backward, coordinatez, a, 1, length(a))
hilbertsort!(a::Array{T,1}, lo::Int64, hi::Int64, lim::Int64) where {T<:AbstractPoint3D} = hilbertsort!(backward, backward, backward, coordinatey, a, lo, hi, lim)

# multi-scale sort. Read all about it here:
# http://doc.cgal.org/latest/Spatial_sorting/classCGAL_1_1Multiscale__sort.html
function _mssort!(a::Array{T,1}, lim_ms::Int64, lim_hl::Int64, rat::Float64) where T<:AbstractPoint
    hi = length(a)
    lo = 1
    while true
        lo = hi - round(Int, (1-rat)*hi)
        hi-lo <= lim_ms && return a
        hilbertsort!(a, lo, hi, lim_hl)
        hi = lo-1
    end
    return a
end

# Utility methods, setting some different defaults for 2D and 3D. These are exported
mssort!(a::Array{T,1}; lim_ms::Int64=16, lim_hl::Int64=4, rat::Float64=0.25) where {T<:AbstractPoint2D} =
    _mssort!(a, lim_ms, lim_hl, rat)
mssort!(a::Array{T,1}; lim_ms::Int64=64, lim_hl::Int64=8, rat::Float64=0.125) where {T<:AbstractPoint3D} =
    _mssort!(a, lim_ms, lim_hl, rat)

end # module GeometicalPredicates


text:
###########################################################################
#                        Key Steps Implementation                         #
###########################################################################

# This function finds the line index (i,j) that the corresponding 2-by-2 W matrix attains the maximum λmin
function find_min_eigen(node::BB.AbstractNode)::Tuple{Int64, Int64}
    root = find_root(node)
    pm = root.auxiliary_data["PM"]
    bus_ids = PM.ids(pm, :bus)
    lookup_w_index = Dict((bi,i) for (i,bi) in enumerate(bus_ids))
    wr = PM.var(pm, :WR)
    wi = PM.var(pm, :WI)
    max_lambda = -Inf
    max_id = ()
    node_solution = node.solution
    eigenvalues = Dict()
    # for raw_k in bus_ids, raw_l in bus_ids
    #     k = lookup_w_index[raw_k]
    #     l = lookup_w_index[raw_l]
    #     if k < l
    #         lambda = 0.5 * (node_solution[wr[k,k]] + node_solution[wr[l,l]] - norm( [node_solution[wr[k,k]] - node_solution[wr[l,l]], 2 * node_solution[wr[k,l]], 2 * node_solution[wi[k,l]]] ) )
    #         println("Smallest eigenvalue for index pair ($(k), $(l)): ", lambda)
    #     end
    # end
    for (raw_i,raw_j) in PM.keys(root.auxiliary_data["Cuts"]["angle_lb"])
        i = lookup_w_index[raw_i]
        j = lookup_w_index[raw_j]
        lambda = 0.5 * (node_solution[wr[i,i]] + node_solution[wr[j,j]] - norm( [node_solution[wr[i,i]] - node_solution[wr[j,j]], 2 * node_solution[wr[i,j]], 2 * node_solution[wi[i,j]]] ) )
        eigenvalues[(raw_i, raw_j)] = lambda
        # println(lambda)
        if lambda > max_lambda
            max_id = (raw_i,raw_j)
            max_lambda = lambda
        end
    end
    node.auxiliary_data["eigenvalues"] = eigenvalues
    # println()
    # println(max_lambda)
    return max_id
end

# this is the interface for finding the complex entry to branch
function find_branching_entry(branch::SpatialBCBranch, node::BB.AbstractNode)::Tuple{Int64, Int64}
    root = find_root(node)
    model = root.model
    opt_idx = solve_candidate_nodes(branch, node)
    return opt_idx
end

# this implements weak branching to find the best, and returns the index of the complex entry to branch on
function solve_candidate_nodes(branch::SpatialBCBranch, node::BB.AbstractNode)
    function _make_fake_branch(branch::SpatialBCBranch, mode::String)
        (Lii, Uii, Ljj, Ujj, Lij, Uij) = branch.bounds
        if mode == "Wii_up" Lii = (Lii + Uii) / 2
        elseif mode == "Wii_down" Uii = (Lii + Uii) / 2
        elseif mode == "Wjj_up" Ljj = (Ljj + Ujj) / 2
        elseif mode == "Wjj_down" Ujj = (Ljj + Ujj) / 2
        elseif mode == "Wij_up" Lij = (Lij + Uij) / 2
        elseif mode == "Wij_down" Uij = (Lij + Uij) / 2
        end
        new_bounds = (Lii, Uii, Ljj, Ujj, Lij, Uij)
        new_π = compute_π(new_bounds)
        return SpatialBCBranch(branch.i, branch.j, branch.wii, branch.wjj, branch.wr, branch.wi, nothing,
                                new_bounds, new_π )
    end

    function _add_constraints_from_fake_branch!(model::JuMP.Model, branch::SpatialBCBranch)
        i = branch.i
        j = branch.j
        (Lii, Uii, Ljj, Ujj, Lij, Uij) = branch.bounds
        πs = branch.valid_ineq_coeffs
        wii = branch.wii
        wjj = branch.wjj
        wr = branch.wr
        wi = branch.wi
        new_branches = []
        JuMP.set_lower_bound(wii, Lii)
        JuMP.set_upper_bound(wii, Uii)
        JuMP.set_lower_bound(wjj, Ljj)
        JuMP.set_upper_bound(wjj, Ujj)
        JuMP.@constraint(model, Lij * wr <= wi)
        JuMP.@constraint(model, Uij * wr >= wi)
        JuMP.@constraint(model, πs[1] + πs[2] * wii + πs[3] * wjj + πs[4] * wr + πs[5] * wi >= Ujj * wii + Uii * wjj - Uii * Ujj)
        JuMP.@constraint(model, πs[1] + πs[2] * wii + πs[3] * wjj + πs[4] * wr + πs[5] * wi >= Ljj * wii + Lii * wjj - Lii * Ljj)
    end
    i = branch.i
    j = branch.j
    root = find_root(node)
    modes = ["Wii", "Wjj", "Wij"]
    best_score = -Inf
    best_mode = ""
    for mode in modes
        fake_branch = _make_fake_branch(branch, mode * "_up")
        model = JuMP.Model(MOSEK_OPTIMIZER)
        # model = JuMP.Model(SCS_OPTIMIZER)
        fake_branch.wii = JuMP.@variable(model, wii)
        fake_branch.wjj = JuMP.@variable(model, wjj)
        fake_branch.wr = JuMP.@variable(model, wr)
        fake_branch.wi = JuMP.@variable(model, wi)
        JuMP.@variable(model, λ)
        JuMP.@constraint(model, [wii + wjj - 2 * λ, wii - wjj, 2 * wr, 2 * wi] in JuMP.SecondOrderCone())
        JuMP.@objective(model, Max, λ)
        _add_constraints_from_fake_branch!(model, fake_branch)
        JuMP.optimize!(model)
        λ_up = JuMP.value(λ)

        fake_branch = _make_fake_branch(branch, mode * "_down")
        model = JuMP.Model(MOSEK_OPTIMIZER)
        fake_branch.wii = JuMP.@variable(model, wii)
        fake_branch.wjj = JuMP.@variable(model, wjj)
        fake_branch.wr = JuMP.@variable(model, wr)
        fake_branch.wi = JuMP.@variable(model, wi)
        JuMP.@variable(model, λ)
        JuMP.@constraint(model, [wii + wjj - 2 * λ, wii - wjj, 2 * wr, 2 * wi] in JuMP.SecondOrderCone())
        JuMP.@objective(model, Max, λ)
        _add_constraints_from_fake_branch!(model, fake_branch)
        JuMP.optimize!(model)
        λ_down = JuMP.value(λ)

        # compute score, update if better
        μ = 0.15
        score = μ * max(-λ_up, -λ_down) + (1-μ) * min(-λ_up, -λ_down)
        if score > best_score
            best_score = score
            best_mode = mode
        end
    end

    if best_mode == "Wii"
        return (i,i)
    elseif best_mode == "Wjj"
        return (j,j)
    else # "Wij"
        return (i,j)
    end
end

# this tracks all the branches above node and add constraints to model accordingly
# if there are multiple branches on the same complex entry, only add constraints for the one with largest depth (closest to node)
function backtracking!(model::JuMP.Model, node::BB.AbstractNode)
    node.auxiliary_data["prev_branch_crefs"] = JuMP.ConstraintRef[]
    root = find_root(node)
    pnode = node
    already_modified = []
    while !isnothing(pnode.parent)
        (i,j) = (pnode.branch.i, pnode.branch.j)
        if !((i,j) in already_modified)
            crefs = add_constraints_from_branch!(model, pnode.branch, root)
            for cref in crefs
                push!(node.auxiliary_data["prev_branch_crefs"], cref)
            end
            push!(already_modified, (i,j))
        end
        pnode = pnode.parent
    end
end

# this modifies constraints in model based on branch
function add_constraints_from_branch!(model::JuMP.Model, branch::SpatialBCBranch, root::BB.AbstractNode)::Vector{JuMP.ConstraintRef}
    i = branch.i
    j = branch.j
    (Lii, Uii, Ljj, Ujj, Lij, Uij) = branch.bounds
    πs = branch.valid_ineq_coeffs
    wii = branch.wii
    wjj = branch.wjj
    wr = branch.wr
    wi = branch.wi
    new_branches = []
    JuMP.set_lower_bound(wii, Lii)
    JuMP.set_upper_bound(wii, Uii)
    JuMP.set_lower_bound(wjj, Ljj)
    JuMP.set_upper_bound(wjj, Ujj)
    JuMP.set_normalized_coefficient(root.auxiliary_data["Cuts"]["angle_lb"][(i,j)], wr, Lij)
    JuMP.set_normalized_coefficient(root.auxiliary_data["Cuts"]["angle_ub"][(i,j)], wr, Uij)
    push!(new_branches, JuMP.@constraint(model, πs[1] + πs[2] * wii + πs[3] * wjj + πs[4] * wr + πs[5] * wi >= Ujj * wii + Uii * wjj - Uii * Ujj))
    push!(new_branches, JuMP.@constraint(model, πs[1] + πs[2] * wii + πs[3] * wjj + πs[4] * wr + πs[5] * wi >= Ljj * wii + Lii * wjj - Lii * Ljj))
    return new_branches
end

# this resets constraints (2) and deletes constraints (3) in model
function delete_prev_branch_constr!(model::JuMP.Model, node::BB.AbstractNode)
    root = find_root(node)
    pm = root.auxiliary_data["PM"]
    w = PM.var(pm, :w)
    wr = PM.var(pm, :wr)
    for (i,_) in PM.ref(pm, :bus)
        JuMP.set_lower_bound(w[i], root.auxiliary_data["Lii"][i])
        JuMP.set_upper_bound(w[i], root.auxiliary_data["Uii"][i])
    end
    for (pair,_) in PM.ref(pm, :buspairs)
        JuMP.set_normalized_coefficient(root.auxiliary_data["Cuts"]["angle_lb"][pair], wr[pair], root.auxiliary_data["Lij"][pair])
        JuMP.set_normalized_coefficient(root.auxiliary_data["Cuts"]["angle_ub"][pair], wr[pair], root.auxiliary_data["Uij"][pair])
    end
    while !isempty(node.auxiliary_data["prev_branch_crefs"])
        cref = Base.pop!(node.auxiliary_data["prev_branch_crefs"])
        JuMP.delete(model, cref)
    end
end

max_stars_repo_name: kibaekkim/ComplexBC

max_issues_repo_name: kibaekkim/ComplexBC

max_forks_repo_name: kibaekkim/ComplexBC


text:
using Plots
using LaTeXStrings

# backend(:pgfplotsx)
backend(:pyplot)
# backend(:gr)

#reset past font scaling
Plots.scalefontsizes()
#scale fonts
Plots.scalefontsizes(5.00)
#Extracting stored varibles

#Add colors
lightblue = colorant"rgb(63%,74%,78%)"
orange = colorant"rgb(85%,55%,13%)"
silver = colorant"rgb(69%,67%,66%)"
rust = colorant"rgb(72%,26%,6%)"

f(x,y) = 2 * x * y - 4 * sum((1 - y).^2) + 100 * Float64(x < 0)
flocal(x,y) = f(0.0, 1.0) + 2 * (x - 0.0) + 2 * (x - 0.0) * (y - 1.0)

N = 32
u = range(0, stop=2π, length=N)
v = range(0, stop=π, length=N)
x = cos.(u) * sin.(v)'
y = sin.(u) * sin.(v)'
z = repeat(cos.(v)',outer=[N, 1])

x_t = -0.75:0.15:0.75
y_t = -0.75:0.15:0.75
z_t = one.(x_t)

#
ν = 0.92 * 2 * π
v_e = range(0, stop=1.25 * π/2, length=100)
x_e = cos(ν) * sin.(v_e)
y_e = sin(ν) * sin.(v_e)
z_e = cos.(v_e)

srf = wireframe(x, y, z .-1.0, camera=(25,40), linecolor=silver, linealpha=0.2,linewidth=0.50, xticks=[], yticks=[], zticks=[], xlims=[-1.5,1.5], ylims=[-1.5,1.5], zlims=[-2.0,0.0], legend=false)

wireframe!(srf, x_t, y_t, (x,y) -> 0.0, linecolor=lightblue, label=L"\mathrm{approximate} \ f)", legend=:topright)

scatter!(srf, [0.0], [0.0], [0.0], color=orange, label="")
quiver!(srf, [0.0], [0.0], [0.0], color=orange, label="", quiver=(1.50 * [cos(ν)], 1.50 * [sin(ν)]), arrow = (0.5, 0.2), linewidth=2.0)
plot!(srf, x_e, y_e, z_e .- 1.00, color=rust, label="", linewidth=2.0)

savefig(srf, "./out/figures/manifold.pdf")
run(`pdfcrop --margins '-65 -35 -65 -55' ./out/figures/manifold.pdf ./out/figures/manifold.pdf`)

max_stars_repo_name: f-t-s/CMD

max_issues_repo_name: f-t-s/CMD

max_forks_repo_name: f-t-s/CMD


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule tmux_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("tmux")
JLLWrappers.@generate_main_file("tmux", UUID("558ad2b7-f92e-50c4-8d35-b0432ac1efd4"))
end  # module tmux_jll


text:
module GDBM

export DBM

include("types.jl")
include("low_level.jl")
include("dict.jl")

end


text:
square(x) = x*x
cube(x) = x*x*x
pow(x::T1, y::T2) where {T1,T2} = x^T1(y)
root(x::T1, y::T2) where {T1,T2} = x^inv(T1(y))
sqrt(x) = root(x, 2.0)
cbrt(x) = root(x, 3.0)
rsqrt(x) = sqrt(inv(x))
rcbrt(x) = cbrt(inv(x))


text:
#------------------------------------------------------------------------------------------------------------------------
# Implements solution to the diffusion equation in an infinite medium as given in [1].
# Solutions are given in the steady-state, frequency, and time domains for an isotroptic source.
#
# [1] Patterson et. al., "Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties,"
#     Appl. Opt. 28, 2331-2336 (1989)
#-------------------------------------------------------------------------------------------------------------------------
""" Structure containing inputs for simulating the fluence under the diffusion approximation in the infinite space."""
struct DAinf_params{T <: AbstractFloat} <: DiffusionParameters
    ρ::T                                 # distance away from isotropic point source (cm)
    μa::T                                # absorption coefficient (cm⁻¹)
    μsp::T                               # reduced scattering coefficient (cm⁻¹)
    n_med::T                             # medium's index of refraction
    ω::T                                 # modulation frequency (1/ns)
    t::Union{T, AbstractVector{T}}       # time vector (ns)

    # do not need to provide these arguments (calculated from previous inputs)
    D::T                                 # Diffusion coefficient
    ν::T                                 # Speed of light in medium (cm/ns)
    function DAinf_params{T}(ρ::T, μa::T, μsp::T, n_med::T, ω::T, t::Union{T, AbstractVector{T}}) where {T <: AbstractFloat}
        @assert ρ > zero(T) "ρ must be positive"
        @assert μa >= zero(T) "μa must greater than or equal to 0"
        @assert μsp > zero(T) "μsp must greater than 0"
        @assert n_med > zero(T) "n_med must be positive"
        @assert all(t .> zero(T)) "t must be positive"
        return new{T}(ρ, μa, μsp, n_med, ω, t, D_coeff(μsp), ν_coeff(n_med))
    end
end

"""
    Generator function for DAinf_params structure.

Provides parameters to use in the infinite space fluence calculation in either the CW, FD, or TD.
Arguments are given as keyword arguments (key = value)

# Keyword Arguments
- `ρ`: source-detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)
- `n_med`: medium's index of refraction
- `ω`: modulation frequency (1/ns)
- `t`: the time vector (ns).

# Examples
```
julia> data = DAinf_params() # return default parameters
julia> data = DAinf_params(ρ = 1.5) # return ρ = 1.5 with the rest of the parameters given by defaults
julia> fluence_DA_inf_CW(data) # can then be used with corresponding functions
```
"""
function DAinf_params(;
    ρ::T = 1.0,
    μa::T = 0.1,
    μsp::T = 10.0,
    n_med::T = 1.0,
    ω::T = 0.0,
    t::Union{T, AbstractVector{T}} = 1.0
) where {T<:AbstractFloat}
    return DAinf_params{T}(ρ, μa, μsp, n_med, ω, t)
end

#-------------------------------------
# Steady-State Fluence
#-------------------------------------
"""
    fluence_DA_inf_CW(ρ, μa, μsp)

Compute the steady-state fluence in an infinite medium.
This is an unsafe implementation and will not check parameters. ρ and μsp should be > 0.0 and μa >= 0.0.

# Arguments
- `ρ`: source-detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)

# Examples
```
julia> fluence_DA_inf_CW(1.0, 0.1, 10.0)
```
"""
function fluence_DA_inf_CW(ρ, μa, μsp)
    params = DiffusionKernelParams(μsp)
    return _kernel_fluence_DA_inf_CW(ρ, μa, μsp, params.D)
end
"""
    fluence_DA_inf_CW(data::DiffusionParameters)

Wrapper to `fluence_DA_inf_CW(ρ, μa, μsp)`

# Examples
```
julia> data = DAinf_params(ρ = 1.0) # use structure to generate inputs
julia> fluence_DA_inf_CW(data) # then call the function
```
"""
function fluence_DA_inf_CW(data::DiffusionParameters)
    return _kernel_fluence_DA_inf_CW(data.ρ, data.μa, data.μsp, data.D)
end

_kernel_fluence_DA_inf_CW(ρ, μa, μsp, D) = exp(-sqrt(3 * μsp * μa) * ρ) / (4 * π * ρ * D)

#-------------------------------------
# Frequency-Domain Fluence
#-------------------------------------
"""
    fluence_DA_inf_FD(ρ, μa, μsp; ω = 0.0, n_med = 1.0)

Compute the fluence for a frequency modulated source in an infinite medium.

# Arguments
- `ρ`: source-detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)

# Keyword Arguments
- `ω`: modulation frequency (1/ns)
- `n_med`: medium's index of refraction

# Examples
```
julia> fluence_DA_inf_FD(1.0, 0.1, 10.0, ω = 1.0, n_med = 1.4)
```
"""
function fluence_DA_inf_FD(ρ, μa, μsp; ω = 0.0, n_med = 1.0)
    params = DiffusionKernelParams(μsp, n_med)
    μa_complex = μa + ω * im / params.ν
    return _kernel_fluence_DA_inf_CW(ρ, μa_complex, μsp, params.D)
end
"""
    fluence_DA_inf_FD(data::DiffusionParameters)

Wrapper to `fluence_DA_inf_FD(ρ, μa, μsp; ω = 0.0, n_med = 1.0)`

# Examples
```
julia> data = DAinf_params(ω = 1.0) # use structure to generate inputs
julia> fluence_DA_inf_FD(data) # then call the function
```
"""
function fluence_DA_inf_FD(data::DiffusionParameters)
    μa_complex = data.μa + data.ω * im / data.ν
    return _kernel_fluence_DA_inf_CW(data.ρ, μa_complex, data.μsp, data.D)
end

#-------------------------------------
# Time-Domain Fluence
#-------------------------------------
"""
    fluence_DA_inf_TD(t, ρ, μa, μsp; n_med = 1.0)

Compute the time-domain fluence in an infinite medium with Eqn. 3 of Patterson. et al. 1989.
This is an unsafe implementation and will not check parameters. t, ρ and μsp should be > 0.0 and μa >= 0.0.

# Arguments
- `t`: the time vector (ns).
- `ρ`: source-detector separation (cm⁻¹)
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)

# Keyword Arguments
- `n_med`= 1.0: medium's index of refraction

# Examples
```
julia> fluence_DA_inf_TD(0.1:0.5:5.0, 1.0, 0.1, 10.0, n_med = 1.4)
```
"""
function fluence_DA_inf_TD(t, ρ, μa, μsp; n_med = 1.0)
    params = DiffusionKernelParams(μsp, n_med)
    return map(t -> _kernel_fluence_DA_inf_TD(t, ρ, μa, params.ν, params.D), t)
end
"""
    fluence_DA_inf_TD(data::DiffusionParameters)

Wrapper to `fluence_DA_inf_TD(t, ρ, μa, μsp; n_med = 1.0)`.

# Examples
```
julia> data = DAinf_params(t = 0.1:0.05:2.0) # use structure to generate inputs
julia> fluence_DA_inf_TD(data) # then call the function
```
"""
function fluence_DA_inf_TD(data::DiffusionParameters)
    return map(t -> _kernel_fluence_DA_inf_TD(t, data.ρ, data.μa, data.ν, data.D), data.t)
end

@inline function _kernel_fluence_DA_inf_TD(t, ρ, μa, ν, D)
    tmp1 = 4 * D * ν * t
    ϕ = exp(-(ρ^2 / tmp1) - μa * ν * t)
    ϕ *= ν / ((tmp1 * π )^(3/2))
end

max_issues_repo_name: heltonmc/LightPropagation


text:
using MPSKit,MPSKitModels,TensorKit

mpo = nonsym_ising_mpo();

state = InfiniteMPS([ℂ^2],[ℂ^10]);

(state,envs,_) = leading_boundary(state,mpo,Vumps(tol_galerkin=1e-10));


text:
module BasicBSpline

using LinearAlgebra
using IntervalSets
using FastGaussQuadrature
using GeometryBasics

export Point
export Knots, ♯, BSplineSpace, dim, ⊑, ⊒, ⋢, ⋣, ≃
export bsplinebasis₊₀, bsplinebasis₋₀, bsplinebasis
export bsplinebasis′₊₀, bsplinebasis′₋₀, bsplinebasis′
export bsplinesupport, bsplineunity, changebasis
export BSplineManifold, refinement, bsplinespaces, controlpoints
export isproper, properdim
export degree, knots
export 𝔫
export FastBSplineSpace, FastBSplineManifold
export BSplineCurve, BSplineSurface, BSplineSolid
export AbstractBSplineManifold, AbstractBSplineSpace
export fittingcontrolpoints

const MAX_DEGREE = 5

include("_Knots.jl")
include("_BSplineSpace.jl")
include("_FastBSplineSpace.jl")
include("_BSplineBasis.jl")
include("_FastBSplineBasis.jl")
include("_BSplineManifold.jl")
include("_FastBSplineManifold.jl")
include("_Refinement.jl")
include("_Integral.jl")
include("_Fitting.jl")

end # module


text:
module CellFieldsTests

using Test
using BlockArrays
using FillArrays
using Gridap.Arrays
using Gridap.Fields
using Gridap.TensorValues
using Gridap.Fields: MockField, MockBasis, OtherMockBasis
using Gridap.CellData

np = 3
ndofs = 4

p = Point(1,2)
x = fill(p,np)
z = 2.0

v = VectorValue(3.0,1.5)
w = VectorValue(3.4,3.5)
a = MockBasis{2}(v,ndofs)
b = MockBasis{2}(w,ndofs)
c = fill(1.0,ndofs)
f = OtherMockBasis{2}(ndofs)

g = MockField{2}(v)

l = 10
xl = Fill(x,l)
zl = [ z for  i in 1:l]
cl = fill(c,l)
fl = Fill(f,l)
ϕl = lincomb(fl,cl)
gl = fill(g,l)
al = Fill(a,l)
bl = fill(b,l)

gf = GenericCellField(gl,ϕl,Val(true))
gf_x = evaluate(gf,xl)
∇gf_x = evaluate(∇(gf),xl)
test_cell_field(gf,xl,gf_x,grad=∇gf_x)

af = GenericCellField(al,ϕl,Val(true),Fill((Base.OneTo(ndofs),),l),Val((:,)))
bf = GenericCellField(bl,ϕl,Val(true),Fill((Base.OneTo(ndofs),),l),Val((:,)))
zf = convert_to_cell_field(zl,ϕl)
df = af*zf
dft = trialize_cell_basis(df)

zfr = reindex(zf,[1,4,3])
@test length(zfr) == 3
@test length(get_array(zfr)) == 3
@test length(get_cell_map(zfr)) == 3
@test length(get_cell_axes(zfr)) == 3

# Check memoization
df_x1 = evaluate(df,xl)
df_x2 = evaluate(df,xl)
@test df_x1 === df_x2
∇gf1 = ∇(gf)
∇gf2 = ∇(gf)
@test ∇gf1 === ∇gf2
@test evaluate(∇gf1,xl) === evaluate(∇gf2,xl)
εgf1 = ε(gf)
εgf2 = ε(gf)
@test εgf1 === εgf2
@test ∇×gf === ∇×gf
@test evaluate(∇×gf,xl) === evaluate(∇×gf,xl)

@test is_test(af)
@test is_trial(dft)
@test is_basis(af)
@test is_basis(dft)
mf = af⋅dft
@test get_metasize(mf) == (:,:)
mf_x = evaluate(mf,xl)
@test size(mf_x[1]) == (np,ndofs,ndofs)

@test get_cell_axes(mf) == Fill((Base.OneTo(ndofs),Base.OneTo(ndofs)),l)

idsL = [ i*collect(1:ndofs)  for i in 1:l]
idsR = [ 2*i*collect(1:ndofs)  for i in 1:l]
axesL = Fill((Base.OneTo(ndofs),),l)
axesR = Fill((Base.OneTo(ndofs),),l)

idsS = merge_cell_dofs_at_skeleton(idsL,idsR,axesL,axesR)
@test isa(idsS,VectorOfBlockArrayCoo)

afS = merge_cell_fields_at_skeleton(af,2*af)
@test isa(afS,SkeletonCellField)

afL_x = evaluate(afS.left,xl)
afR_x = evaluate(afS.right,xl)
@test isa(afL_x,VectorOfBlockArrayCoo)
@test isa(afR_x,VectorOfBlockArrayCoo)

@test isa(afS*2,SkeletonCellField)
@test isa(afS+afS,SkeletonCellField)

# Checks associated with trial bases
df = bf
dft = trialize_cell_basis(df)
@test is_trial(dft) == true
cell_vals = [rand(ndofs) for  i in 1:l]
cell_field = lincomb(df,cell_vals)
@test isa(cell_field,CellField)
@test get_metasize(cell_field) == ()
cell_field_x = evaluate(cell_field,xl)
@test isa(cell_field_x[1],AbstractVector)

end # module


text:
using TemporalGPs: storage_type, is_of_storage_type

println("gauss_markov:")
@testset "gauss_markov" begin

    Dlats = [1, 7]
    Ns = [1, 11]
    tvs = [true, false]
    storages = [
        (name="Array{Float64}", val=ArrayStorage(Float64)),
        (name="SArray{Float64}", val=SArrayStorage(Float64)),
    ]

    # Dlatss = [7]
    # Ns = [11]
    # tvs = [true]

    @testset "time_varying=$tv, Dlat=$Dlat, N=$N, storage=$(storage.name)" for
        tv in tvs,
        Dlat in Dlats,
        N in Ns,
        storage in storages

        rng = MersenneTwister(123456)
        gmm = tv == true ?
            random_tv_gmm(rng, Forward(), Dlat, N, storage.val) :
            random_ti_gmm(rng, Forward(), Dlat, N, storage.val)

        @test eltype(gmm) == eltype(storage.val)
        @test storage_type(gmm) == storage.val

        @test length(gmm) == N
        @test getindex(gmm, N) isa TemporalGPs.SmallOutputLGC

        @test is_of_storage_type(gmm, storage.val)

        @testset "==" begin
            gmm_other = tv == true ?
                random_tv_gmm(rng, Forward(), Dlat, N, storage.val) :
                random_ti_gmm(rng, Forward(), Dlat, N, storage.val)
            @test gmm == gmm
            @test gmm != gmm_other
        end
    end
end