text:
module Scale

using Colors
using Compat
using Compose
using DataArrays
using DataStructures
using Gadfly
using Showoff
using IndirectArrays
using CategoricalArrays

import Gadfly: element_aesthetics, isconcrete, concrete_length, discretize_make_ia
import Distributions: Distribution

include("color_misc.jl")

# Return true if var is categorical.
iscategorical(scales::Dict{Symbol, Gadfly.ScaleElement}, var::Symbol) =
        haskey(scales, var) && isa(scales[var], DiscreteScale)

# Apply some scales to data in the given order.
#
# Args:
#   scales: An iterable object of ScaleElements.
#   aess: Aesthetics (of the same length as datas) to update with scaled data.
#   datas: Zero or more data objects. (Yes, I know "datas" is not a real word.)
#
# Returns:
#   nothing
#
function apply_scales(scales, aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    for scale in scales
        apply_scale(scale, aess, datas...)
    end

    for (aes, data) in zip(aess, datas)
        aes.titles = data.titles
    end
end

# Apply some scales to data in the given order.
#
# Args:
#   scales: An iterable object of ScaleElements.
#   datas: Zero or more data objects.
#
# Returns:
#   A vector of Aesthetics of the same length as datas containing scaled data.
#
function apply_scales(scales, datas::Gadfly.Data...)
    aess = Gadfly.Aesthetics[Gadfly.Aesthetics() for _ in datas]
    apply_scales(scales, aess, datas...)
    aess
end

# Transformations on continuous scales
struct ContinuousScaleTransform
    f::Function     # transform function
    finv::Function  # f's inverse

    # A function taking one or more values and returning an array of
    # strings.
    label::Function
end

identity_formatter(xs::AbstractArray, format=:auto) = showoff(xs, format)
const identity_transform = ContinuousScaleTransform(identity, identity, identity_formatter)

log10_formatter(xs::AbstractArray, format=:plain) =
        [@sprintf("10%s", x) for x in showoff(xs, format)]
const log10_transform = ContinuousScaleTransform(log10, x -> 10^x, log10_formatter)

log2_formatter(xs::AbstractArray, format=:plain) =
        [@sprintf("2%s", x) for x in showoff(xs, format)]
const log2_transform = ContinuousScaleTransform(log2, x -> 2^x, log2_formatter)

ln_formatter(xs::AbstractArray, format=:plain) =
        [@sprintf("e%s", x) for x in showoff(xs, format)]
const ln_transform = ContinuousScaleTransform(log, exp, ln_formatter)

asinh_formatter(xs::AbstractArray, format=:plain) =
        [@sprintf("sinh(%s)", x) for x in showoff(xs, format)]
const asinh_transform = ContinuousScaleTransform(asinh, sinh, asinh_formatter)

sqrt_formatter(xs::AbstractArray, format=:plain) =
        [@sprintf("%s2", x) for x in showoff(xs, format)]
const sqrt_transform = ContinuousScaleTransform(sqrt, x -> x^2, sqrt_formatter)

# Continuous scale maps data on a continuous scale simple by calling
# `convert(Float64, ...)`.
struct ContinuousScale <: Gadfly.ScaleElement
    vars::Vector{Symbol}
    trans::ContinuousScaleTransform
    minvalue
    maxvalue
    minticks
    maxticks
    labels::Union{(Void), Function}
    format
    scalable

    function ContinuousScale(vars, trans,
                             minvalue, maxvalue, minticks, maxticks,
                             labels, format, scalable)
        minvalue != nothing && maxvalue != nothing && minvalue > maxvalue &&
                error("Cannot construct a ContinuousScale with minvalue > maxvalue")
        new(vars, trans, minvalue, maxvalue, minticks, maxticks, labels, format, scalable)
    end
end

function ContinuousScale(vars, trans;
                         labels=nothing,
                         minvalue=nothing, maxvalue=nothing,
                         minticks=2, maxticks=10, format=nothing, scalable=true)
    ContinuousScale(vars, trans, minvalue, maxvalue, minticks, maxticks, labels, format, scalable)
end

function make_labeler(scale::ContinuousScale)
    if scale.labels != nothing
        xs -> [scale.labels(x) for x in xs]
    elseif scale.format == nothing
        scale.trans.label
    else
        xs -> scale.trans.label(xs, scale.format)
    end
end

const x_vars = [:x, :xmin, :xmax, :xintercept, :intercept, :xviewmin, :xviewmax, :xend]
const y_vars = [:y, :ymin, :ymax, :yintercept, :slope, :middle, :upper_fence, :lower_fence,
                :upper_hinge, :lower_hinge, :yviewmin, :yviewmax, :yend]

function continuous_scale_partial(vars::Vector{Symbol}, trans::ContinuousScaleTransform)
    function(; minvalue=nothing, maxvalue=nothing, labels=nothing, format=nothing, minticks=2,
                 maxticks=10, scalable=true)
        ContinuousScale(vars, trans, minvalue=minvalue, maxvalue=maxvalue,
                        labels=labels, format=format, minticks=minticks,
                        maxticks=maxticks, scalable=scalable)
    end
end

# Commonly used scales.
const x_continuous = continuous_scale_partial(x_vars, identity_transform)
const y_continuous = continuous_scale_partial(y_vars, identity_transform)
const x_log10      = continuous_scale_partial(x_vars, log10_transform)
const y_log10      = continuous_scale_partial(y_vars, log10_transform)
const x_log2       = continuous_scale_partial(x_vars, log2_transform)
const y_log2       = continuous_scale_partial(y_vars, log2_transform)
const x_log        = continuous_scale_partial(x_vars, ln_transform)
const y_log        = continuous_scale_partial(y_vars, ln_transform)
const x_asinh      = continuous_scale_partial(x_vars, asinh_transform)
const y_asinh      = continuous_scale_partial(y_vars, asinh_transform)
const x_sqrt       = continuous_scale_partial(x_vars, sqrt_transform)
const y_sqrt       = continuous_scale_partial(y_vars, sqrt_transform)

const size_continuous = continuous_scale_partial([:size], identity_transform)

element_aesthetics(scale::ContinuousScale) = scale.vars

# Apply a continuous scale.
#
# Args:
#   scale: A continuous scale.
#   datas: Zero or more data objects.
#   aess: Aesthetics (of the same length as datas) to update with scaled data.
#
# Return:
#   nothing
#
function apply_scale(scale::ContinuousScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    for (aes, data) in zip(aess, datas)
        for var in scale.vars
            vals = getfield(data, var)
            if vals isa CategoricalArray
                throw(ArgumentError("continuous scale for $var aesthetic when stored as a CategoricalArray. Consider using a discrete scale or convert data to an Array."))
            end
            vals === nothing && continue

            # special case for function arrays bound to :y
            # pass the function values through and wait for the scale to
            # be reapplied by Stat.func
            if var == :y && eltype(vals) == Function
                aes.y = vals
                continue
            end

            # special case for Distribution values bound to :x or :y. wait for
            # scale to be re-applied by Stat.qq
            if in(var, [:x, :y]) && typeof(vals) <: Distribution
                setfield!(aes, var, vals)
                continue
            end

            T = Any
            for d in vals
                if isconcrete(d)
                    T = typeof(scale.trans.f(d))
                    break
                end
            end

            if T <: Measure
                T = Measure
            end

            ds = any(ismissing, vals) ? DataArray(T, length(vals)) : Array{T}(length(vals))
            apply_scale_typed!(ds, vals, scale)

            if var == :xviewmin || var == :xviewmax ||
               var == :yviewmin || var == :yviewmax
                setfield!(aes, var, ds[1])
            else
                setfield!(aes, var, ds)
            end

            if var in x_vars
                label_var = :x_label
            elseif var in y_vars
                label_var = :y_label
            else
                label_var = Symbol(var, "_label")
            end

            if in(label_var, Set(fieldnames(aes)))
                setfield!(aes, label_var, make_labeler(scale))
            end
        end

        if scale.minvalue != nothing
            if scale.vars === x_vars
                aes.xviewmin = scale.trans.f(scale.minvalue)
            elseif scale.vars === y_vars
                aes.yviewmin = scale.trans.f(scale.minvalue)
            end
        end

        if scale.maxvalue != nothing
            if scale.vars === x_vars
                aes.xviewmax = scale.trans.f(scale.maxvalue)
            elseif scale.vars === y_vars
                aes.yviewmax = scale.trans.f(scale.maxvalue)
            end
        end
    end
end

function apply_scale_typed!(ds, field, scale::ContinuousScale)
    for i in 1:length(field)
        d = field[i]
        ds[i] = isconcrete(d) ? scale.trans.f(d) : d
    end
end

function discretize(values, levels=nothing, order=nothing, preserve_order=true)
    if levels == nothing
        if preserve_order
            levels = OrderedSet()
            for value in values
                push!(levels, value)
            end
            da = discretize_make_ia(values, collect(eltype(values), levels))
        else
            da = discretize_make_ia(values)
        end
    else
        da = discretize_make_ia(values, levels)
    end

    if order != nothing
        return discretize_make_ia(da, da.values[order])
    else
        return da
    end
end

struct DiscreteScaleTransform
    f::Function
end

struct DiscreteScale <: Gadfly.ScaleElement
    vars::Vector{Symbol}

    # Labels are either a function that takes an array of values and returns
    # an array of string labels, a vector of string labels of the same length
    # as the number of unique values in the discrete data, or nothing to use
    # the default labels.
    labels::Union{(Void), Function}

    # If non-nothing, give values for the scale. Order will be respected and
    # anything in the data that's not represented in values will be set to missing.
    levels::Union{(Void), AbstractVector}

    # If non-nothing, a permutation of the pool of values.
    order::Union{(Void), AbstractVector}
end
DiscreteScale(vals; labels=nothing, levels=nothing, order=nothing) =
        DiscreteScale(vals, labels, levels, order)

const discrete = DiscreteScale

element_aesthetics(scale::DiscreteScale) = scale.vars

x_discrete(; labels=nothing, levels=nothing, order=nothing) =
        DiscreteScale(x_vars, labels=labels, levels=levels, order=order)

y_discrete(; labels=nothing, levels=nothing, order=nothing) =
        DiscreteScale(y_vars, labels=labels, levels=levels, order=order)

group_discrete(; labels=nothing, levels=nothing, order=nothing) =
        DiscreteScale([:group], labels=labels, levels=levels, order=order)

shape_discrete(; labels=nothing, levels=nothing, order=nothing) =
        DiscreteScale([:shape], labels=labels, levels=levels, order=order)

size_discrete(; labels=nothing, levels=nothing, order=nothing) =
        DiscreteScale([:size], labels=labels, levels=levels, order=order)

function apply_scale(scale::DiscreteScale, aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    for (aes, data) in zip(aess, datas)
        for var in scale.vars
            label_var = Symbol(var, "_label")
            getfield(data, var) === nothing && continue

            disc_data = discretize(getfield(data, var), scale.levels, scale.order)
            setfield!(aes, var, discretize_make_ia(Int64.(disc_data.index)))

            # The leveler for discrete scales is a closure over the discretized data.
            if scale.labels === nothing
                function default_labeler(xs)
                    lvls = filter(!ismissing, disc_data.values)
                    vals = Any[1 <= x <= length(lvls) ? lvls[x] : "" for x in xs]
                    if all([isa(val, AbstractFloat) for val in vals])
                        return showoff(vals)
                    else
                        return [string(val) for val in vals]
                    end
                end
                labeler = default_labeler
            else
                function explicit_labeler(xs)
                    lvls = filter(!ismissing, disc_data.values)
                    return [string(scale.labels(lvls[x])) for x in xs]
                end
                labeler = explicit_labeler
            end

            in(label_var, Set(fieldnames(aes))) && setfield!(aes, label_var, labeler)
        end
    end
end

struct NoneColorScale <: Gadfly.ScaleElement
end

const color_none = NoneColorScale

element_aesthetics(scale::NoneColorScale) = [:color]

function apply_scale(scale::NoneColorScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    for aes in aess
        aes.color = nothing
    end
end

struct IdentityColorScale <: Gadfly.ScaleElement
end

const color_identity = IdentityColorScale

element_aesthetics(scale::IdentityColorScale) = [:color]

function apply_scale(scale::IdentityColorScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    for (aes, data) in zip(aess, datas)
        data.color === nothing && continue
        aes.color = discretize_make_ia(data.color)
        aes.color_key_colors = Dict()
    end
end

struct DiscreteColorScale <: Gadfly.ScaleElement
    f::Function # A function f(n) that produces a vector of n colors.

    # If non-nothing, give values for the scale. Order will be respected and
    # anything in the data that's not represented in values will be set to missing.
    levels::Union{(Void), AbstractVector}

    # If non-nothing, a permutation of the pool of values.
    order::Union{(Void), AbstractVector}

    # If true, order levels as they appear in the data
    preserve_order::Bool
end
DiscreteColorScale(f; levels=nothing, order=nothing, preserve_order=true) =
        DiscreteColorScale(f, levels, order, preserve_order)

element_aesthetics(scale::DiscreteColorScale) = [:color]

function default_discrete_colors(n)
    convert(Vector{Color},
         distinguishable_colors(n, [LCHab(70, 60, 240)],
             transform=c -> deuteranopic(c, 0.5),
             lchoices=Float64[65, 70, 75, 80],
             cchoices=Float64[0, 50, 60, 70],
             hchoices=linspace(0, 330, 24),
         )
     )
end

# Common discrete color scales
function color_discrete_hue(f=default_discrete_colors;
                            levels=nothing,
                            order=nothing,
                            preserve_order=true)
    DiscreteColorScale(
        f,
        levels=levels,
        order=order,
        preserve_order=preserve_order)
end

@deprecate discrete_color_hue(; levels=nothing, order=nothing, preserve_order=true) color_discrete_hue(; levels=levels, order=order, preserve_order=preserve_order)

const color_discrete = color_discrete_hue

@deprecate discrete_color(; levels=nothing, order=nothing, preserve_order=true) color_discrete(; levels=levels, order=order, preserve_order=preserve_order)

color_discrete_manual(colors::AbstractString...; levels=nothing, order=nothing) =
        color_discrete_manual(Gadfly.parse_colorant(colors)...; levels=levels, order=order)

function color_discrete_manual(colors::Color...; levels=nothing, order=nothing)
    cs = [colors...]
    f = n -> distinguishable_colors(n, cs)
    DiscreteColorScale(f, levels=levels, order=order)
end

@deprecate discrete_color_manual(colors...; levels=nothing, order=nothing) color_discrete_manual(colors...; levels=levels, order=order)

function apply_scale(scale::DiscreteColorScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    levelset = OrderedSet()
    for (aes, data) in zip(aess, datas)
        data.color === nothing && continue
        for d in data.color
            # Remove missing values
            # FixMe! The handling of missing values shouldn't be this scattered across the source
            ismissing(d) || push!(levelset, d)
        end
    end

    if scale.levels == nothing
        scale_levels = [levelset...]
        scale.preserve_order || sort!(scale_levels)
    else
        scale_levels = scale.levels
    end
    scale.order == nothing || permute!(scale_levels, scale.order)
    colors = convert(Vector{RGB{Float32}}, scale.f(length(scale_levels)))

    color_map = Dict([(color, string(label))
                              for (label, color) in zip(scale_levels, colors)])
    labeler(xs) = [color_map[x] for x in xs]

    for (aes, data) in zip(aess, datas)
        data.color === nothing && continue
        # Remove missing values
        # FixMe! The handling of missing values shouldn't be this scattered across the source
        ds = discretize([c for c in data.color if !ismissing(c)], scale_levels)

        colorvals = colors[ds.index]

        colored_ds = discretize_make_ia(colorvals, colors)
        aes.color = colored_ds

        aes.color_label = labeler
        aes.color_key_colors = OrderedDict()
        for (i, c) in enumerate(colors)
            aes.color_key_colors[c] = i
        end
    end
end

struct ContinuousColorScale <: Gadfly.ScaleElement
    # A function of the form f(p) where 0 <= p <= 1, that returns a color.
    f::Function
    trans::ContinuousScaleTransform

    minvalue
    maxvalue
end
ContinuousColorScale(f, trans=identity_transform; minvalue=nothing, maxvalue=nothing) =
        ContinuousColorScale(f, trans, minvalue, maxvalue)

function continuous_color_scale_partial(trans::ContinuousScaleTransform)
    lch_diverge2 = function(l0=30, l1=100, c=40, h0=260, h1=10, hmid=20, power=1.5)
        lspan = l1 - l0
        hspan1 = hmid - h0
        hspan0 = h1 - hmid
        function(r)
            r2 = 2r - 1
            return LCHab(min(80, l1 - lspan * abs(r2)^power), max(10, c * abs(r2)),
                         (1-r)*h0 + r * h1)
        end
    end

    function(; minvalue=nothing, maxvalue=nothing, colormap=lch_diverge2())
        ContinuousColorScale(colormap, trans, minvalue=minvalue, maxvalue=maxvalue)
    end
end

# Commonly used scales.
const color_continuous = continuous_color_scale_partial(identity_transform)
const color_log10      = continuous_color_scale_partial(log10_transform)
const color_log2       = continuous_color_scale_partial(log2_transform)
const color_log        = continuous_color_scale_partial(ln_transform)
const color_asinh      = continuous_color_scale_partial(asinh_transform)
const color_sqrt       = continuous_color_scale_partial(sqrt_transform)

const color_continuous_gradient = color_continuous

element_aesthetics(::ContinuousColorScale) = [:color]

@deprecate continuous_color_gradient(;minvalue=nothing, maxvalue=nothing) color_continuous_gradient(;minvalue=minvalue, maxvalue=maxvalue)

@deprecate continuous_color(;minvalue=nothing, maxvalue=nothing) color_continuous(;minvalue=nothing, maxvalue=nothing)

function apply_scale(scale::ContinuousColorScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    cmin, cmax = Inf, -Inf
    for data in datas
        data.color === nothing && continue

        for c in data.color
            ismissing(c) && continue

            c = convert(Float64, c)
            if c < cmin
                cmin = c
            end

            if c > cmax
                cmax = c
            end
        end
    end

    (cmin == Inf || cmax == -Inf) && return

    if scale.minvalue != nothing
        cmin = scale.minvalue
    end

    if scale.maxvalue  != nothing
        cmax = scale.maxvalue
    end

    cmin, cmax = promote(cmin, cmax)

    cmin = scale.trans.f(cmin)
    cmax = scale.trans.f(cmax)

    ticks, viewmin, viewmax = Gadfly.optimize_ticks(cmin, cmax)
    if ticks[1] == 0 && cmin >= 1
        ticks[1] = 1
    end

    cmin, cmax = ticks[1], ticks[end]
    cspan = cmax != cmin ? cmax - cmin : 1.0

    for (aes, data) in zip(aess, datas)
        data.color === nothing && continue

        aes.color = Array{RGB{Float32}}(length(data.color))
        apply_scale_typed!(aes.color, data.color, scale, cmin, cspan)

        color_key_colors = Dict{Color, Float64}()
        color_key_labels = Dict{Color, AbstractString}()

        tick_labels = scale.trans.label(ticks)
        for (i, j, label) in zip(ticks, ticks[2:end], tick_labels[1:end-1])
            r = (i - cmin) / cspan
            c = scale.f(r)
            color_key_colors[c] = r
            color_key_labels[c] = label
        end
        c = scale.f((ticks[end] - cmin) / cspan)
        color_key_colors[c] = (ticks[end] - cmin) / cspan
        color_key_labels[c] = tick_labels[end]

        labeler(xs) = [get(color_key_labels, x, "") for x in xs]

        aes.color_function = scale.f
        aes.color_label = labeler
        aes.color_key_colors = color_key_colors
        aes.color_key_continuous = true
    end
end

function apply_scale_typed!(ds, field, scale::ContinuousColorScale,
                            cmin::Float64, cspan::Float64)
    for (i, d) in enumerate(field)
        if isconcrete(d)
            ds[i] = convert(RGB{Float32},
                        scale.f((convert(Float64, scale.trans.f(d)) - cmin) / cspan))
        else
            ds[i] = missing
        end
    end
end

# Label scale is always discrete, hence we call it 'label' rather
# 'label_discrete'.
struct LabelScale <: Gadfly.ScaleElement
end

function apply_scale(scale::LabelScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    for (aes, data) in zip(aess, datas)
        data.label === nothing && continue
        aes.label = discretize(data.label)
    end
end

element_aesthetics(::LabelScale) = [:label]

const label = LabelScale

# Scale applied to grouping aesthetics.
struct GroupingScale <: Gadfly.ScaleElement
    var::Symbol
end

xgroup(; labels=nothing, levels=nothing, order=nothing) =
        DiscreteScale([:xgroup], labels=labels, levels=levels, order=order)

ygroup(; labels=nothing, levels=nothing, order=nothing) =
        DiscreteScale([:ygroup], labels=labels, levels=levels, order=order)

# Catchall scale for when no transformation of the data is necessary
struct IdentityScale <: Gadfly.ScaleElement
    var::Symbol
end

element_aesthetics(scale::IdentityScale) = [scale.var]

function apply_scale(scale::IdentityScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    for (aes, data) in zip(aess, datas)
        getfield(data, scale.var) === nothing && continue
        setfield!(aes, scale.var, getfield(data, scale.var))
    end
end

z_func() = IdentityScale(:z)
y_func() = IdentityScale(:y)
x_distribution() = IdentityScale(:x)
y_distribution() = IdentityScale(:y)
shape_identity() = IdentityScale(:shape)
size_identity() = IdentityScale(:size)

end # module Scale


text:
abstract GoogleChart

## Concrete type for Core charts
type CoreChart <: GoogleChart
    packages::Vector
    chart_type::AbstractString
    id::AbstractString
    width::Integer
    height::Integer
    options::Dict
    data::Any
    xtra::AbstractString   ## xtra javascript as a string
end

## Charts from
## https://developers.google.com/chart/interactive/docs/gallery
## and
## https://developers.google.com/chart/interactive/docs/more_charts
charts = (
          (:area_chart,     "AreaChart",       ["corechart"], Dict(:title => "Area chart"),
           "https://developers.google.com/chart/interactive/docs/gallery/areachart")
          ,(:bar_chart,     "BarChart",        ["corechart"], Dict(:title => "Bar chart"),
            "https://developers.google.com/chart/interactive/docs/gallery/areachart")
          ,(:bubble_chart,  "BubbleChart",     ["corechart"], Dict(:title => "Bubble chart", :bubble=>Dict(:textStyle=>Dict(:fontSize=>11))),
            "https://developers.google.com/chart/interactive/docs/gallery/bubblechart")
          ,(:candlestick_chart,    "CandlestickChart",   ["corechart"], Dict(:title => "Candlestick chart"),
            "https://developers.google.com/chart/interactive/docs/gallery/candlestickchart")
          ,(:column_chart,   "ColumnChart",    ["corechart"], Dict(:title => "Column chart"),
            "https://developers.google.com/chart/interactive/docs/gallery/columnchart")
          ,(:combo_chart,    "ComboChart",     ["corechart"], Dict(:title => "Combo chart"),
            "https://developers.google.com/chart/interactive/docs/gallery/combochart")
          ,(:gauge_chart,    "Gauge",     ["gauge"], Dict(),
            "https://developers.google.com/chart/interactive/docs/gallery/gaugechart")
          ,(:geo_chart,      "GeoChart",       ["geochart"], Dict(:title => "Geo chart"),
            "https://developers.google.com/chart/interactive/docs/gallery/geochart")
          ,(:line_chart,     "LineChart",      ["corechart"], Dict(:title => "Line chart"),
            "https://developers.google.com/chart/interactive/docs/gallery/linechart")
          ,(:pie_chart,      "PieChart",       ["corechart"], Dict(:title => "Pie chart"),
            "https://developers.google.com/chart/interactive/docs/gallery/piechart")
          ,(:scatter_chart,  "ScatterChart",   ["corechart"], Dict(:title => "Scatter chart"),
            "https://developers.google.com/chart/interactive/docs/gallery/scatterchart")
          ,(:stepped_area_chart,    "SteppedAreaChart",   ["corechart"], Dict(:title => "Stepped area chart"),
            "https://developers.google.com/chart/interactive/docs/gallery/steppedareachart")
          ,(:table_chart,    "Table",     ["table"], Dict(:title => "Table chart"),
            "https://developers.google.com/chart/interactive/docs/gallery/tablechart")
          ,(:tree_chart,     "TreeMap",      ["treemap"], Dict(:title => "Tree chart"),
            "https://developers.google.com/chart/interactive/docs/gallery/treechart")
          ## See note about viewing from a file here: https://developers.google.com/chart/interactive/docs/gallery/motionchart#Overview
          ,(:annotated_time_line,    "AnnotatedTimeLine",   ["annotatedtimeline"], Dict(:title => "Annotated time line"),
            "https://developers.google.com/chart/interactive/docs/gallery/annotatedtimeline")
          ,(:intensity_map,  "IntensityMap",   ["intensitymap"], Dict(:title => "Intensity map"),
            "https://developers.google.com/chart/interactive/docs/gallery/intensitymap")
          ,(:motion_chart,   "MotionChart",    ["motionchart"], Dict(:width=>900,:height=>600),
            "https://developers.google.com/chart/interactive/docs/gallery/motionchart")
          ,(:org_chart,      "OrgChart",       ["orgchart"], Dict(:title => "Org chart"),
            "https://developers.google.com/chart/interactive/docs/gallery/orgchart")
          ,(:image_spark_line,    "ImageSparkLine",   ["imagesparkline"], Dict(),
            "https://developers.google.com/chart/interactive/docs/gallery/imagesparkline")
          )

## Make constructors
## e.g.
## line_type(data, opts, width, height)
## line_type(data, opts) ## default 900, 600
## line_type(data)       ## deafult title, 900, 600
function merge_dict(d1, d2)
    d = Dict{Symbol, Any}()
    for (k, v) in d1 d[k] = v end
    for (k, v) in d2 d[k] = v end
    d
end
for (nm, ctype, packages, defaults, url) in charts
    @eval begin
        function $(nm)(data::DataFrame, opts::MaybeDict, width::MaybeInt, height::MaybeInt)
            obj = CoreChart($packages, string($ctype), get_id(), 900, 600, $defaults, nothing, "" )
            if isa(opts, Dict) obj.options = merge_dict(obj.options, opts) end
            if isa(width, Integer) obj.width = width end
            if isa(height, Integer) obj.height = height end
            obj.data = make_data_array(obj.id, data)
            obj
        end
        $(nm)(data::DataFrame, opts::MaybeDict) = $(nm)(data, opts, nothing, nothing)
        function $(nm)(data::DataFrame, width::Int, height::Int; kwargs...)
            d = Dict()
            [d[s]=v for (s,v) in kwargs]
            $(nm)(data, d, width, height)
        end
        $(nm)(data::DataFrame; kwargs...) = $(nm)(data, 900, 600; kwargs...)
        $(nm)(data::DataFrame) = $(nm)(data, nothing, nothing, nothing)
    end
end

## get help for chart named...
function help_on_chart(nm::AbstractString)
    default_url = "https://developers.google.com/chart/interactive/docs/gallery"
    for (name, ctype, packages, defaults, url) in charts
        if nm == string(name)
            open_url(url)
            return(nothing)
        end
    end
    println("Chart $nm not found. Opening list of charts ...")
    open_url(default_url)
end

ext: jl


text:
module DaggerGPU

using Dagger, MemPool, Requires, Adapt
using Distributed
using KernelAbstractions

import Dagger: Chunk

macro gpuproc(PROC, T)
    quote
        # Assume that we can run anything
        Dagger.iscompatible_func(proc::$PROC, opts, f) = true
        Dagger.iscompatible_arg(proc::$PROC, opts, x) = true

        # CPUs shouldn't process our array type
        Dagger.iscompatible_arg(proc::Dagger.ThreadProc, opts, x::$T) = false

        # Adapt to/from the appropriate type
        function Dagger.move(from_proc::OSProc, to_proc::$PROC, x::Chunk)
            from_pid = from_proc.pid
            to_pid = Dagger.get_parent(to_proc).pid
            @assert myid() == to_pid
            adapt($T, remotecall_fetch(from_pid, x) do x
                poolget(x.handle)
            end)
        end
        function Dagger.move(from_proc::$PROC, to_proc::OSProc, x::Chunk)
            from_pid = Dagger.get_parent(from_proc).pid
            to_pid = to_proc.pid
            @assert myid() == to_pid
            remotecall_fetch(from_pid, x) do x
                adapt(Array, poolget(x.handle))
            end
        end
        function Dagger.move(from_proc::OSProc, to_proc::$PROC, x)
            adapt($T, x)
        end
        function Dagger.move(from_proc::$PROC, to_proc::OSProc, x)
            adapt(Array, x)
        end
    end
end

processor(kind::Symbol) = processor(Val(kind))
processor(::Val) = Dagger.ThreadProc
cancompute(kind::Symbol) = cancompute(Val(kind))
cancompute(::Val) = false

kernel_backend() = kernel_backend(Dagger.Sch.thunk_processor())
kernel_backend(::Dagger.ThreadProc) = CPU()

function __init__()
    @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" begin
        include("cu.jl")
    end
    @require AMDGPU="21141c5a-9bdb-4563-92ae-f87d6854732e" begin
        include("roc.jl")
    end
end

end


text:
module PartialSvdStoch

using Match

using Logging
using Base.CoreLogging

using Printf
using LinearAlgebra
using LowRankApprox
using Base.CoreLogging
using Statistics

debug_log = stdout
logger = ConsoleLogger(stdout, CoreLogging.Info)
global_logger(logger)

export VrPCA,
        SvdMode,
        LowRankApproxMc,
        reduceToRank,
        reduceToEpsil,
        iterateVrPCA,
        getindex,
        getResidualError,
        SvdApprox,
        isConverged

include("vrpca.jl")
include("lowRankApproxMc.jl")

end


text:
### A Pluto.jl notebook ###
# v0.15.0

using Markdown
using InteractiveUtils

# ╔═╡ f4cc5c93-624b-45ca-a4da-fa299c691c8c
begin
	using Pkg

	Pkg.activate()
end

# ╔═╡ 23dd09a8-d5d0-11eb-2029-c3cd330fdfa2
begin
	using SatisfiabilityInterface

	using Symbolics
	using Symbolics: Sym, Term

	using PlutoUI

	using LightGraphs, GraphPlot, Colors
	using Compose
end

# ╔═╡ a70a8996-513c-4bbd-aaf8-255d349ee739
using SimpleSATSolver

# ╔═╡ 25d094ff-bfe9-48af-b6b9-f2203ad49e47
using Latexify

# ╔═╡ e2ece52b-efbd-41de-b98e-b9e4e52bf033
using Luxor

# ╔═╡ fbcd65de-d0a4-4914-84f7-5bbbfacd0aa9
emph(s) = HTML("$s")

# ╔═╡ 533cff90-d873-4c99-8e63-3f0a77811cf6
# TableOfContents()

# ╔═╡ 1a9aad3a-bc69-4a77-bf28-c268a0bc1737
set_default_graphic_size(6*Compose.cm, 6*Compose.cm)  # size for Compose figures

# ╔═╡ 22297817-e948-467c-9a4a-8a1cee20c1e9
html"<button onclick=present()>Present</button>"

# ╔═╡ ec63ded2-3d52-4256-bf8e-638be5133902
md"""
# Solving discrete problems using $(emph("Boolean satisfiability")) with Julia
"""

# ╔═╡ 5c1ac720-70b4-495a-a66c-4ed952fbd137
md"""
### David P. Sanders

#### Facultad de Ciencias $(HTML("
")) Universidad Nacional Autónoma de México
#### Department of Mathematics & Julia Lab, MIT
"""

# ╔═╡ f4e269bd-1e67-4bb3-9dc0-fa4c51eae012
md"""
## Outline
"""

# ╔═╡ 67713a33-da7d-495e-911b-051fc34bf638
md"""

- A $(emph("symbolic")) modelling interface for discrete problems
"""

# ╔═╡ 9635873c-7362-4b38-ad9e-1e1eacb49220
md"""
- How to $(emph("encode")) problems as Boolean satisfiability problems
"""

# ╔═╡ fc64bb2e-7c4b-4467-ad68-8d141ae67484
md"""
- How to $(emph("solve")) Boolean satisfiability problems
"""

# ╔═╡ f7598c9b-2673-4b8b-9b2f-b44d648c2b66
HTML("
")

# ╔═╡ f42290e3-27f8-4a74-99f0-ed823c5dd16f
md"""
- Related work:

  - `JuliaConstraints` organisation

  - `ConstraintProgrammingExtensions.jl` for JuMP / MOI

"""

# ╔═╡ 7a4f1c53-53cf-44be-9ac0-fdc960bdae51
md"""
## Motivation: Colouring $(emph("graphs"))
"""

# ╔═╡ a1723672-5ce5-4fe8-a808-2cc1da3c8c60
md"""
- A $(emph("graph")) has a set `V` of $(emph("vertices")):
"""

# ╔═╡ d0477e2c-9be5-46de-a85f-bbf94c3f3fe7
V = [1, 2, 3]

# ╔═╡ 13a4d93a-114a-4ccd-aaee-28c80fb0606f
md"""
- Joined by $(emph("edges")) in the set `E`:
"""

# ╔═╡ 90789020-b24c-41a8-847b-e6309babd31f
E = [(1, 2), (2, 3)]

# ╔═╡ f47c8993-1d2a-4baf-88c7-fa5e605e34b5
gplot(V, E)

# ╔═╡ e5ea5a80-af37-48ce-91ca-4bc2652c1673
md"""
Graphs rendered using `LightGraphs.jl` and `GraphPlot.jl`
"""

# ╔═╡ 0f305f88-26ca-4acb-9696-ca0bc31ab979
md"""
## Graph $(emph("colouring"))
"""

# ╔═╡ 6a1342de-507b-4008-bedb-5f387f17a850
md"""
- The $(emph("graph colouring")) problem is:
"""

# ╔═╡ e6bd1f4e-e23a-45a0-be82-2f47c1f665cf
md"""

> - Assign a colour to each vertex
"""

# ╔═╡ 729aae42-b4f2-493a-ade8-71a7b790561f
md"""
>   - Such that two vertices joined by an edge have *distinct* colours
"""

# ╔═╡ d61c615e-c751-490c-9889-506c8dee30aa
md"""
## $(emph("Modelling")) the colouring problem
"""

# ╔═╡ c4cc188f-1476-4d39-b4cb-3f5632d23c79
md"""
- We want a mathematical $(emph("model")) for the colouring problem
"""

# ╔═╡ fd63036a-c674-4f30-b8cc-21b67528dd42
md"""
- We use symbolic $(emph("variables")) from `Symbolics.jl`
"""

# ╔═╡ 12628428-1bb1-4f45-8909-88a1d71b8ee7
md"""
- We call $\quad c_i$ = colour of vertex $i$
"""

# ╔═╡ 57e9e379-6ef1-4cc8-ae7b-bd17b4e8e11b
n = length(V)   # number of vertices

# ╔═╡ 4c2c7c62-3b60-48cb-96db-4020885dee1c
@variables c[1:n];    # from Symbolics.jl

# ╔═╡ 8ac35f21-1ed6-4b2e-b019-4ddfd5830851
[c[i] for i in 1:n]

# ╔═╡ 170aaa75-7a8b-4886-90ad-21f417aac55b
md"""
## Domain constraints
"""

# ╔═╡ 759127e1-0703-4fa1-b61f-f4c26c126cc6
md"""
- We need a $(emph("domain")) for each variable
"""

# ╔═╡ b5b9a579-3c53-4724-8f3c-339254fd522d
md"""
- Let's try to colour the graph with 2 colours:
"""

# ╔═╡ c54d16eb-d06a-4dbc-a3e7-c3513da710ee
k = 2

# ╔═╡ e5efb464-042e-48fc-9224-c6505d8e04a3
colours = [:red, :green, :yellow, :blue, :black]

# ╔═╡ 79370a37-d35c-42e5-b6dd-abfd60711fa9
md"""
- We can write domains $(emph("symbolically")):
"""

# ╔═╡ fec874c5-f32d-4988-a6e0-f449d362fc9c
c[1] ∈ colours[1:k]

# ╔═╡ 182e2eac-cf13-4225-a003-cc2c82553ae5
domains = [c[i] ∈ colours[1:k] for i in 1:n]

# ╔═╡ 3fb4ca38-7f53-4251-aa5e-d183b03376af
md"""
## Edge constraints
"""

# ╔═╡ 7ade9ef2-416f-438c-b126-421db2379776
md"""
- If nodes $i$ and $j$ are joined by an edge then they must have distinct colours:
"""

# ╔═╡ d1070d1c-f74d-4bd3-b0ce-862d848fd565
restrictions = [c[i] ≠ c[j] for (i, j) in E]

# ╔═╡ 6e821823-1f30-4b2e-a3e3-4de2a550dc19
md"""
- Let's collect all the constraints together:
"""

# ╔═╡ de3fd8b8-37c2-424e-98be-45824154edd3
constraints = domains ∪ restrictions

# ╔═╡ 0181743f-1a2b-4d44-85b1-77bf08edbf19
md"""
## Constraint satisfaction problem
"""

# ╔═╡ 58ce7a7f-8797-41c1-af47-314157dcd873
md"""
- A $(emph("constraint satisfaction problem")) (CSP) is a collection of variables, domains and constraints
"""

# ╔═╡ 49977007-7b86-426c-bdfa-4fc218bf1e7e
md"""
- We make a `DiscreteCSP` object from `SatisfiabilityInterface.jl` using the constraints:
"""

# ╔═╡ 8d214801-c65d-4d71-bc0c-b05852343428
prob = DiscreteCSP(constraints)

# ╔═╡ c8acfa87-a671-437d-bb9b-b776d014c491
md"""
## Solving the constraint satisfaction problem
"""

# ╔═╡ c600bd16-b655-422d-bb28-7dc124e0ba5a
md"""
- We call `solve` to solve the constraint satisfaction problem:
"""

# ╔═╡ 53a0cc3a-aee5-479e-9a56-9d6563c64c8a
status, results = solve(prob);

# ╔═╡ 760fda23-0027-412d-a3c5-aab300631bc8
results

# ╔═╡ 1df9c28b-54ef-4a1a-adbe-2610861a8559
gplot(V, E, prob, results)

# ╔═╡ d8092ba3-57b9-46aa-ae32-9f66444edb5b
md"""
## Ring graph
"""

# ╔═╡ bb584a01-1051-4b4d-936e-a64e83e58f17
md"""
- Let's make a graph in the form of a ring:
"""

# ╔═╡ 9f6ae168-e89f-4e29-bbf1-057d0dc66c1a
md"""
- Let's try to colour it with 2 colours:
"""

# ╔═╡ b8fe62d5-70e3-4bc7-9d59-9c8ac19ae56f
md"""
	- We see that this is $(emph("unsatisfiable"))!
	"""

# ╔═╡ 4bf219bd-dc77-44e8-91af-9aa67bb26c24
md"""

	## Colouring a ring again

	- But of course it can be done with *three* colours:
	"""

# ╔═╡ 7b86ec36-6e8d-45a7-bbc7-81c8eb0037e0
md"""
## Colouring a Penrose tiling -- @Cormullion
"""

# ╔═╡ f41caa13-d354-4a42-9bc4-f7ba12f9fffe
PlutoUI.LocalResource("penrose.png")

# ╔═╡ 0110824f-c664-4312-81a4-2cdb967f926b
md"""
## How do we actually $(emph("solve")) these problems?
"""

# ╔═╡ c9db9c18-d1f6-47f4-8e8f-2f158a6eedde
md"""
- We **encode** the problem to a $(emph("boolean satisfiability")) problem (SAT)
"""

# ╔═╡ ef59a7e9-c952-4081-86d0-e5091cf950e5
md"""
- Simple encoding of a variable $x$ with domain 1:$n$ :

  - Use $n$ **Boolean** variables $x_1, \ldots, x_n$

  - Restrict such that **exactly one** of the $x_i$ is true
"""

# ╔═╡ f6769053-14b4-444b-8761-a217ebf3fbf2
encoded = encode(prob)

# ╔═╡ 6a95e734-7f1c-4364-a28c-f7ac7807336a
md"""
## SAT problem in CNF form
"""

# ╔═╡ 5f8ae350-aa2b-4331-9443-af183e355529
md"""
- CNF = $(emph("Conjunctive Normal Form"))

- A conjunction (AND) of clauses

- Each clause is an OR
"""

# ╔═╡ 7442450f-8919-4fe9-ba94-36a27135997d
encoded.p

# ╔═╡ d6d38c75-d38a-4009-acbd-9195370bcbcc
md"""
## SAT solvers
"""

# ╔═╡ a59ffbb7-56a4-4723-a40a-561ed4b995a5
md"""
- The Boolean satisfiability problem is $(emph("NP-complete"))
"""

# ╔═╡ 647b6d87-62ca-43fa-b772-afa8ab9555ac
md"""
- But it can be solved efficiently in many cases, by one of many $(emph("SAT solvers")) that are available
"""

# ╔═╡ e891d669-9431-4ee2-8eb0-d4921c1b2d38
md"""
- Usually written in C or C++
"""

# ╔═╡ f2dbb4d7-5c49-46e7-9535-dc10500e7e61
md"""
- `PicoSAT.jl` wraps the popular `PicoSAT` solver
"""

# ╔═╡ ea0028ed-df33-4f05-adb0-036bc07a8eb2
md"""
- `SatisfiabilityInterface.jl` currently calls out to `cryptominisat5` by default
"""

# ╔═╡ c4e2bb24-462c-450b-a145-c1fea2bd9ee8
md"""
## `SimpleSATSolver.jl`
"""

# ╔═╡ 17f423cf-a241-48f3-a718-51b17ab6a268
md"""
- A very simple SAT solver in $(emph("pure Julia"))
"""

# ╔═╡ bb498e44-559f-414e-ab11-e055037a4b87
md"""
- Implements the $(emph("DPLL algorithm")) (1962)
"""

# ╔═╡ f7061757-93e3-4405-80d3-92f854697e5b
md"""
- Meant to be pedagogical
"""

# ╔═╡ 958cea58-82b1-43b3-975d-aade709c718d
md"""
- But seems to be performant for small systems
"""

# ╔═╡ 03f70163-6712-4e92-9a47-ff246e227813
typeof(encoded.p)

# ╔═╡ ae00c44f-b04b-4233-95cf-9c9a8c6dd799
solve(encoded.p, SimpleSAT())

# ╔═╡ 1528a881-8501-4122-be89-bf6b75203a60
md"""
## Algorithm for solving boolean satisfiability
"""

# ╔═╡ 44638a74-7b54-4ef2-8a64-c1e8da6cd553
md"""
- Start assigning variables to `true` or `false`
"""

# ╔═╡ fa9d0583-522b-43ea-90b1-70e902b9d55c
md"""
- Check if any constraints are violated and $(emph("backtrack")) (choose the opposite value) if so
	"""

# ╔═╡ 3932867d-74c0-49cc-840e-2d1421a3867a
md"""
- Apply $(emph("unit propagation")):

   - If $(x \vee y)$ and $x$ is false then $y$ must be true

"""

# ╔═╡ ddb8ef7d-6e90-4699-a635-5351622f79d6
md"""
## Discrete numerical problems
"""

# ╔═╡ 812ba4be-e5b2-42f0-be4e-7ff643284662
md"""
- We can also solve discrete numerical problems, e.g.:
"""

# ╔═╡ d1f99f24-d9fb-409d-8ef6-316797a5f88a
vars = @variables x, y, z

# ╔═╡ 955768ef-c550-47ee-b56a-2dc4bc5baad4
numerical_constraints = [
    x ∈ 1:2
    y ∈ 2:5
    z ∈ 3:7
    z == x * y
	x ≠ 1
]

# ╔═╡ 8ed7146d-ca65-4676-b0c9-f97cf67aff32
md"""
## Solving numerical problems
"""

# ╔═╡ 1edf43b4-418d-4615-82a7-00ff7b1ac41c
prob6 = DiscreteCSP(numerical_constraints);

# ╔═╡ 0b0da897-61bb-40cf-a3f9-b30d95c47563
solve(prob6)   # single solution

# ╔═╡ f03160be-b968-4b19-a369-d48cb7dc010a
all_solutions(prob6)

# ╔═╡ 601d3d5d-e613-4c5a-a164-805ccee5f73d
md"""
## Conclusions
"""

# ╔═╡ 2193ce7e-dfde-4bed-a6fb-e4b4831f36fa
md"""
- Encoding to SAT is a $(emph("powerful paradigm")) for solving discrete problems
"""

# ╔═╡ f9580d63-0bdf-4f3b-9dfb-98838f3b666f
md"""
- Our work aims to make this accessible from Julia
"""

# ╔═╡ 6be5dbf5-bee0-4935-a5c7-e5daaad7360d
md"""

- Our packages:

   - `SatisfiabilityInterface.jl`

   - `SimpleSATSolver.jl`

"""

# ╔═╡ 420da272-eeb9-431c-8bcb-4b8666746655
md"""
# Appendix
"""

# ╔═╡ 97b7b89a-0919-4f81-812a-23d3823b38ba

# ╔═╡ e556a339-a596-49dd-83c3-e0c5c72c7d4d
function convert_to_lightgraphs(V, E)

	g = Graph()
	add_vertices!(g, length(V))

	for (i, j) in E
		add_edge!(g, i, j)
	end

	return g
end

# ╔═╡ b0d6ba28-d60d-4aee-a489-79af197a29de
function GraphPlot.gplot(V, E, c, results; nodelabel=V)
	RGB_colours = Dict(col => parse(RGBA, col) for col in colours)
	colour_indices = [results[c[i]] for i in 1:length(c)]

	nodefillc = [RGB_colours[c] for c in colour_indices]

	g = convert_to_lightgraphs(V, E)

	gplot(g,
		nodelabel=nodelabel,
		nodefillc=nodefillc
	)

end

# ╔═╡ 403216ab-a405-43c3-adf6-926be6b063bc
function GraphPlot.gplot(V, E, prob::DiscreteCSP, results; nodelabel=V)
	gplot(V, E, sort(collect(keys(prob.varmap))), results; nodelabel=nodelabel)
end

# ╔═╡ b48b43ba-115c-4b69-9e9c-481e25f12581
function GraphPlot.gplot(V, E; nodelabel=V)
	p = gplot(convert_to_lightgraphs(V, E), nodelabel=nodelabel)
end

# ╔═╡ 80e4d7f7-eb31-480a-beac-abd4bd8b2798

"k is the number of colours"
function graph_colouring_problem(V, E, k=3)

    # colours = [:red, :green, :yellow, :blue, :black][1:k]
    # c = [Num(Variable(:c, i)) for i in 1:length(V)]   # colour variables

    @variables c[1:length(V)]

    constraints =
    [
        [c[i] ∈ colours[1:k] for i in 1:length(V)]

        [c[i] ≠ c[j] for (i, j) in E]  #  different_colours(E, c)
    ]

    return DiscreteCSP(constraints)
end

# ╔═╡ c1ff01b8-916c-4333-bc34-d0c04debb362
### Fix Latexify output:

# ╔═╡ f796d811-dbaa-4f5f-83a5-5911ac3dae43
begin
	Latexify.latexify(x::Sym)  = Latexify.latexify(Num(x))
	Latexify.latexify(x::Term) = Latexify.latexify(Num(x))
end

# ╔═╡ a16ee047-921f-46bb-bfda-71be201740c3
function ring_graph(n=11)
    V = 1:n
    E = [(i, mod1((i+1), n)) for i in 1:n]

    return V, E
end

# ╔═╡ 10f71d92-0575-4940-9c57-27ea52b5fbb1
V2, E2 = ring_graph(3)

# ╔═╡ a79fb891-8ab0-4c19-b801-b7881e5525ca
gplot(V2, E2)

# ╔═╡ 5ee64af1-a53b-4c66-b4b9-0c11b05cae78
prob2 = graph_colouring_problem(V2, E2, 2);

# ╔═╡ a1e31033-998c-4bfc-9d23-8295d5da742d
status2, results2 = solve(prob2)

# ╔═╡ d9b22906-e72a-43dd-81f2-2241e5dca6ee
prob3 = graph_colouring_problem(V2, E2, 3);

# ╔═╡ 26ef3f82-be05-4545-bc1e-cefe96bcddde
status3, results3 = solve(prob3);

# ╔═╡ 9bdaa197-cf4b-45e2-b231-8ec571b3fc95
gplot(V2, E2, prob3, results3)

# ╔═╡ f801d658-0562-4239-9806-e5eccf239fff
## Penrose tiling -- code by Cormullion

# ╔═╡ cd29bbff-a06c-433f-882d-e18b4987114a
struct PenroseTriangle
    red::Bool # not red, more of a flag
    pointA::Point
    pointB::Point
    pointC::Point
end

# ╔═╡ 307e96c3-a270-4b12-a0ac-f865cb20d49a
function subdivide(triangles;
        type=:P3)
    result = []
    for triangle in triangles
        A, B, C = triangle.pointA, triangle.pointB, triangle.pointC
        if triangle.red == true
            if type == :P3
                #  P3 rhombus
                P = A + (B - A) / MathConstants.golden
                push!(result, PenroseTriangle(true, C, P, B))
                push!(result, PenroseTriangle(false, P, C, A))
            else # P2 half kite
                Q = A + (B - A) / MathConstants.golden
                R = B + (C - B) / MathConstants.golden
                push!(result, PenroseTriangle(false, R, Q, B))
                push!(result, PenroseTriangle(true,  Q, A, R))
                push!(result, PenroseTriangle(true,  C, A, R))
            end
        else
            if type == :P3
                # P3 rhombus
                Q = B + (A - B) / MathConstants.golden
                R = B + (C - B) / MathConstants.golden
                push!(result, PenroseTriangle(true, R, Q, A))
                push!(result, PenroseTriangle(false, R, A, C))
                push!(result, PenroseTriangle(false, Q, R, B))
            else # P2 kite/dart
                P = C + (A - C) / MathConstants.golden
                push!(result, PenroseTriangle(false, B, P, A))
                push!(result, PenroseTriangle(true,  P, C, B))
            end
        end
    end
    return result
end

# ╔═╡ 17ae28f4-ba09-42bb-83f0-51f771d0bbc5
function PenroseTiling(centerpos::Point, radius, depth=4;
        type=:P3)
    triangles = PenroseTriangle[]
    A = centerpos
    n = 10 # a circle of triangles 360/10 -> 36°
    for i in 1:n
        phi = (i - 1) * (2π / n)
        C = A + polar(radius, phi)
        phi = (i) * (2π / n)
        B = A + polar(radius, phi)
        if type == :P3
            if i % 2 == 1
                triangle = PenroseTriangle(true, A, C, B)
            else
                triangle = PenroseTriangle(true, A, B, C)
            end
        else # P2
            if i % 2 == 1
                triangle = PenroseTriangle(true, B, A, C)
            else
                triangle = PenroseTriangle(true, C, A, B)
            end
        end
        push!(triangles, triangle)
    end
    for i in 1:depth
        triangles = subdivide(triangles, type=type)
    end
    return triangles
end

# ╔═╡ 759d23f1-360c-4a36-8cf6-0d1f6c6133cd

function incidence_graph(triangles)
    n = length(triangles)
    g = SimpleGraph(n)
    for i in 1:n
        for j in i+1:n
            ti = triangles[i]
            tj = triangles[j]

            points_i = [ti.pointA, ti.pointB, ti.pointC]
            points_j = [tj.pointA, tj.pointB, tj.pointC]

            coincidental = 0
            tol = 1e-2

            for k in 1:3, l in 1:3
                if distance(points_i[k], points_j[l]) < tol
                    coincidental += 1
                end
            end

            if coincidental == 2  # share edge
                add_edge!(g, i, j)
            end
        end
    end

    return g
end

# ╔═╡ 4a0baabf-2d98-4d1a-9906-b05c0654ce86

# ╔═╡ e4281582-9748-4801-972f-2880d6031f1b
tiling = PenroseTiling(O, 1, 4)

# ╔═╡ f32d038d-c126-4fda-bd1c-f33dc039d042
g = incidence_graph(tiling)

# ╔═╡ 0bb71243-3bb9-4ffb-b062-fd8cffa43701
VV = vertices(g)

# ╔═╡ 6f90a7ba-e4cb-4cb1-a60c-f1918642f2b6
EE = [(e.dst, e.src) for e in edges(g)]

# ╔═╡ 0fcbd390-bc8f-4284-b7d3-61428ab6ec82
prob4 = graph_colouring_problem(VV, EE, 4)

# ╔═╡ cf89ab12-e70f-4533-b7e4-52bdc253c18a
@time status4, results4 = solve(prob4)  # 0.05 seconds

# ╔═╡ 38aa8bc1-5b9a-48ae-90e6-5826c86aaff5
unique(values(results4))

# ╔═╡ ff73ed7c-ef11-4876-a345-15f265e753a8

# ╔═╡ 77645c06-1dde-4a29-86bb-5d62f06ee3b0
set_default_graphic_size(20*Compose.cm, 20*Compose.cm)  # size for Compose figures

# ╔═╡ cc04b735-3bff-4f60-a392-2c3d9e6b0220
gplot(VV, EE; nodelabel=nothing)

# ╔═╡ 70405d58-e2e5-4d55-b145-ddaf5a3ed36f
gplot(VV, EE, prob4, results4; nodelabel=nothing)

# ╔═╡ cb819540-f084-45af-bec2-6afb535ab5f9
g.fadjlist[70]

# ╔═╡ bf2a29da-52ed-44f0-b572-38cbb53bd260
results4[c[73]]

# ╔═╡ 31b3bb83-286b-4937-9e86-bed812991013
V5, E5 = ring_graph(31)

# ╔═╡ 7fcbbd88-785f-4c65-96b2-d74239de6e56
prob5 = graph_colouring_problem(V5, E5, 3)

# ╔═╡ c90239b3-d451-4379-adf5-56262ab20b4e
status5, results5 = solve(prob5)

# ╔═╡ 0f1b0462-e0af-4d01-953b-231de11fdf79
gplot(V5, E5, prob5, results5)

# ╔═╡ 8872ee87-e03f-4445-b706-a0b15d86ba07
sort(collect(keys(prob4.varmap)))

# ╔═╡ fa0c842a-b252-4421-90c9-dc713bb1a040
isless(string(c[100]),  string(c[10]))

# ╔═╡ 42429a49-0c05-4b41-96ed-67f58bb762c8
c[100] < c[10]

# ╔═╡ 875f7bd7-a1bd-44f6-a8bb-9a6464cbcc0f
typeof(c[100])

# ╔═╡ eb29d8b9-37b8-4822-a4ec-5a98bd777892
sort(collect(keys(prob4.varmap)))

# ╔═╡ ece5d0ae-8edb-42a0-a467-9effd1253af8
Base.isless(x::Term, y::Term) = isless(x.args[2], y.args[2])

# ╔═╡ 31cdbb1c-044c-4bb0-9fed-ecbb1f8b126e
Base.isless(x::Num, y::Num) = isless(x.val, y.val)

# ╔═╡ 37cb6adb-c6af-47b6-b663-2e89dc11fc0a
isless(c[1], c[100])

# ╔═╡ 1706b53e-0bab-4206-ad88-570039430839
Dump(d.val)

# ╔═╡ Cell order:
# ╠═f4cc5c93-624b-45ca-a4da-fa299c691c8c
# ╠═23dd09a8-d5d0-11eb-2029-c3cd330fdfa2
# ╠═fbcd65de-d0a4-4914-84f7-5bbbfacd0aa9
# ╠═533cff90-d873-4c99-8e63-3f0a77811cf6
# ╠═1a9aad3a-bc69-4a77-bf28-c268a0bc1737
# ╟─22297817-e948-467c-9a4a-8a1cee20c1e9
# ╟─ec63ded2-3d52-4256-bf8e-638be5133902
# ╟─5c1ac720-70b4-495a-a66c-4ed952fbd137
# ╟─f4e269bd-1e67-4bb3-9dc0-fa4c51eae012
# ╟─67713a33-da7d-495e-911b-051fc34bf638
# ╟─9635873c-7362-4b38-ad9e-1e1eacb49220
# ╟─fc64bb2e-7c4b-4467-ad68-8d141ae67484
# ╟─f7598c9b-2673-4b8b-9b2f-b44d648c2b66
# ╟─f42290e3-27f8-4a74-99f0-ed823c5dd16f
# ╟─7a4f1c53-53cf-44be-9ac0-fdc960bdae51
# ╟─a1723672-5ce5-4fe8-a808-2cc1da3c8c60
# ╠═d0477e2c-9be5-46de-a85f-bbf94c3f3fe7
# ╟─13a4d93a-114a-4ccd-aaee-28c80fb0606f
# ╠═90789020-b24c-41a8-847b-e6309babd31f
# ╠═f47c8993-1d2a-4baf-88c7-fa5e605e34b5
# ╟─e5ea5a80-af37-48ce-91ca-4bc2652c1673
# ╟─0f305f88-26ca-4acb-9696-ca0bc31ab979
# ╟─6a1342de-507b-4008-bedb-5f387f17a850
# ╟─e6bd1f4e-e23a-45a0-be82-2f47c1f665cf
# ╟─729aae42-b4f2-493a-ade8-71a7b790561f
# ╟─d61c615e-c751-490c-9889-506c8dee30aa
# ╟─c4cc188f-1476-4d39-b4cb-3f5632d23c79
# ╟─fd63036a-c674-4f30-b8cc-21b67528dd42
# ╟─12628428-1bb1-4f45-8909-88a1d71b8ee7
# ╠═57e9e379-6ef1-4cc8-ae7b-bd17b4e8e11b
# ╠═4c2c7c62-3b60-48cb-96db-4020885dee1c
# ╠═8ac35f21-1ed6-4b2e-b019-4ddfd5830851
# ╟─170aaa75-7a8b-4886-90ad-21f417aac55b
# ╟─759127e1-0703-4fa1-b61f-f4c26c126cc6
# ╟─b5b9a579-3c53-4724-8f3c-339254fd522d
# ╠═c54d16eb-d06a-4dbc-a3e7-c3513da710ee
# ╠═e5efb464-042e-48fc-9224-c6505d8e04a3
# ╟─79370a37-d35c-42e5-b6dd-abfd60711fa9
# ╠═fec874c5-f32d-4988-a6e0-f449d362fc9c
# ╠═182e2eac-cf13-4225-a003-cc2c82553ae5
# ╟─3fb4ca38-7f53-4251-aa5e-d183b03376af
# ╟─7ade9ef2-416f-438c-b126-421db2379776
# ╠═d1070d1c-f74d-4bd3-b0ce-862d848fd565
# ╟─6e821823-1f30-4b2e-a3e3-4de2a550dc19
# ╠═de3fd8b8-37c2-424e-98be-45824154edd3
# ╟─0181743f-1a2b-4d44-85b1-77bf08edbf19
# ╟─58ce7a7f-8797-41c1-af47-314157dcd873
# ╟─49977007-7b86-426c-bdfa-4fc218bf1e7e
# ╠═8d214801-c65d-4d71-bc0c-b05852343428
# ╟─c8acfa87-a671-437d-bb9b-b776d014c491
# ╟─c600bd16-b655-422d-bb28-7dc124e0ba5a
# ╠═53a0cc3a-aee5-479e-9a56-9d6563c64c8a
# ╠═760fda23-0027-412d-a3c5-aab300631bc8
# ╠═1df9c28b-54ef-4a1a-adbe-2610861a8559
# ╟─d8092ba3-57b9-46aa-ae32-9f66444edb5b
# ╟─bb584a01-1051-4b4d-936e-a64e83e58f17
# ╠═10f71d92-0575-4940-9c57-27ea52b5fbb1
# ╠═a79fb891-8ab0-4c19-b801-b7881e5525ca
# ╟─9f6ae168-e89f-4e29-bbf1-057d0dc66c1a
# ╠═5ee64af1-a53b-4c66-b4b9-0c11b05cae78
# ╠═a1e31033-998c-4bfc-9d23-8295d5da742d
# ╟─b8fe62d5-70e3-4bc7-9d59-9c8ac19ae56f
# ╟─4bf219bd-dc77-44e8-91af-9aa67bb26c24
# ╠═d9b22906-e72a-43dd-81f2-2241e5dca6ee
# ╠═26ef3f82-be05-4545-bc1e-cefe96bcddde
# ╠═9bdaa197-cf4b-45e2-b231-8ec571b3fc95
# ╟─7b86ec36-6e8d-45a7-bbc7-81c8eb0037e0
# ╟─f41caa13-d354-4a42-9bc4-f7ba12f9fffe
# ╟─0110824f-c664-4312-81a4-2cdb967f926b
# ╟─c9db9c18-d1f6-47f4-8e8f-2f158a6eedde
# ╟─ef59a7e9-c952-4081-86d0-e5091cf950e5
# ╠═f6769053-14b4-444b-8761-a217ebf3fbf2
# ╟─6a95e734-7f1c-4364-a28c-f7ac7807336a
# ╟─5f8ae350-aa2b-4331-9443-af183e355529
# ╠═7442450f-8919-4fe9-ba94-36a27135997d
# ╟─d6d38c75-d38a-4009-acbd-9195370bcbcc
# ╟─a59ffbb7-56a4-4723-a40a-561ed4b995a5
# ╟─647b6d87-62ca-43fa-b772-afa8ab9555ac
# ╟─e891d669-9431-4ee2-8eb0-d4921c1b2d38
# ╟─f2dbb4d7-5c49-46e7-9535-dc10500e7e61
# ╟─ea0028ed-df33-4f05-adb0-036bc07a8eb2
# ╟─c4e2bb24-462c-450b-a145-c1fea2bd9ee8
# ╟─17f423cf-a241-48f3-a718-51b17ab6a268
# ╟─bb498e44-559f-414e-ab11-e055037a4b87
# ╟─f7061757-93e3-4405-80d3-92f854697e5b
# ╟─958cea58-82b1-43b3-975d-aade709c718d
# ╠═a70a8996-513c-4bbd-aaf8-255d349ee739
# ╠═03f70163-6712-4e92-9a47-ff246e227813
# ╠═ae00c44f-b04b-4233-95cf-9c9a8c6dd799
# ╟─1528a881-8501-4122-be89-bf6b75203a60
# ╟─44638a74-7b54-4ef2-8a64-c1e8da6cd553
# ╟─fa9d0583-522b-43ea-90b1-70e902b9d55c
# ╟─3932867d-74c0-49cc-840e-2d1421a3867a
# ╟─ddb8ef7d-6e90-4699-a635-5351622f79d6
# ╟─812ba4be-e5b2-42f0-be4e-7ff643284662
# ╠═d1f99f24-d9fb-409d-8ef6-316797a5f88a
# ╠═955768ef-c550-47ee-b56a-2dc4bc5baad4
# ╟─8ed7146d-ca65-4676-b0c9-f97cf67aff32
# ╠═1edf43b4-418d-4615-82a7-00ff7b1ac41c
# ╠═0b0da897-61bb-40cf-a3f9-b30d95c47563
# ╠═f03160be-b968-4b19-a369-d48cb7dc010a
# ╟─601d3d5d-e613-4c5a-a164-805ccee5f73d
# ╟─2193ce7e-dfde-4bed-a6fb-e4b4831f36fa
# ╟─f9580d63-0bdf-4f3b-9dfb-98838f3b666f
# ╟─6be5dbf5-bee0-4935-a5c7-e5daaad7360d
# ╟─420da272-eeb9-431c-8bcb-4b8666746655
# ╠═97b7b89a-0919-4f81-812a-23d3823b38ba
# ╠═e556a339-a596-49dd-83c3-e0c5c72c7d4d
# ╠═b0d6ba28-d60d-4aee-a489-79af197a29de
# ╠═403216ab-a405-43c3-adf6-926be6b063bc
# ╠═b48b43ba-115c-4b69-9e9c-481e25f12581
# ╠═80e4d7f7-eb31-480a-beac-abd4bd8b2798
# ╠═c1ff01b8-916c-4333-bc34-d0c04debb362
# ╠═25d094ff-bfe9-48af-b6b9-f2203ad49e47
# ╠═f796d811-dbaa-4f5f-83a5-5911ac3dae43
# ╠═a16ee047-921f-46bb-bfda-71be201740c3
# ╠═f801d658-0562-4239-9806-e5eccf239fff
# ╠═e2ece52b-efbd-41de-b98e-b9e4e52bf033
# ╠═cd29bbff-a06c-433f-882d-e18b4987114a
# ╠═17ae28f4-ba09-42bb-83f0-51f771d0bbc5
# ╠═307e96c3-a270-4b12-a0ac-f865cb20d49a
# ╠═759d23f1-360c-4a36-8cf6-0d1f6c6133cd
# ╠═4a0baabf-2d98-4d1a-9906-b05c0654ce86
# ╠═e4281582-9748-4801-972f-2880d6031f1b
# ╠═f32d038d-c126-4fda-bd1c-f33dc039d042
# ╠═0bb71243-3bb9-4ffb-b062-fd8cffa43701
# ╠═6f90a7ba-e4cb-4cb1-a60c-f1918642f2b6
# ╠═0fcbd390-bc8f-4284-b7d3-61428ab6ec82
# ╠═cf89ab12-e70f-4533-b7e4-52bdc253c18a
# ╠═38aa8bc1-5b9a-48ae-90e6-5826c86aaff5
# ╠═ff73ed7c-ef11-4876-a345-15f265e753a8
# ╠═77645c06-1dde-4a29-86bb-5d62f06ee3b0
# ╠═cc04b735-3bff-4f60-a392-2c3d9e6b0220
# ╠═70405d58-e2e5-4d55-b145-ddaf5a3ed36f
# ╠═cb819540-f084-45af-bec2-6afb535ab5f9
# ╠═bf2a29da-52ed-44f0-b572-38cbb53bd260
# ╠═31b3bb83-286b-4937-9e86-bed812991013
# ╠═7fcbbd88-785f-4c65-96b2-d74239de6e56
# ╠═c90239b3-d451-4379-adf5-56262ab20b4e
# ╠═0f1b0462-e0af-4d01-953b-231de11fdf79
# ╠═8872ee87-e03f-4445-b706-a0b15d86ba07
# ╠═fa0c842a-b252-4421-90c9-dc713bb1a040
# ╠═42429a49-0c05-4b41-96ed-67f58bb762c8
# ╠═875f7bd7-a1bd-44f6-a8bb-9a6464cbcc0f
# ╠═eb29d8b9-37b8-4822-a4ec-5a98bd777892
# ╠═ece5d0ae-8edb-42a0-a467-9effd1253af8
# ╠═31cdbb1c-044c-4bb0-9fed-ecbb1f8b126e
# ╠═37cb6adb-c6af-47b6-b663-2e89dc11fc0a
# ╠═1706b53e-0bab-4206-ad88-570039430839


text:
using CliffordAlgebras
using Test
import LinearAlgebra.SingularException

@testset "CliffordAlgebras.jl" begin

    @testset "Algebra" begin
        @test CliffordAlgebra(1,0,0) == CliffordAlgebra(1)
        @test CliffordAlgebra(0,1,0) == CliffordAlgebra(0,1)
        @test CliffordAlgebra(0,0,1) == CliffordAlgebra(0,0,1)

        @test_throws AssertionError CliffordAlgebra(-1)
        @test_throws AssertionError CliffordAlgebra(0,-1)
        @test_throws AssertionError CliffordAlgebra(0,0,-1)

        @test_throws AssertionError CliffordAlgebra(0,0,0, (:e1,))
        @test_throws AssertionError CliffordAlgebra(1,0,0, (:e1,:e2,))

        @test typeof(CliffordAlgebra(1,0,0, (:e1,))) <: CliffordAlgebra
        @test typeof(CliffordAlgebra(1,1,0, (:e₊,:e₋))) <: CliffordAlgebra
        @test typeof(CliffordAlgebra(1,1,1, (:e₊,:e₋,:e₀))) <: CliffordAlgebra

        @test CliffordAlgebra(:Hyper) == CliffordAlgebra(:Hyperbolic)
        @test CliffordAlgebra(:Complex) == CliffordAlgebra(:ℂ)
        @test CliffordAlgebra(:Dual) == CliffordAlgebra(:Grassmann)
        @test CliffordAlgebra(:Grassmann2D) == CliffordAlgebra(:G2)
        @test CliffordAlgebra(:Grassmann3D) == CliffordAlgebra(:G3)
        @test CliffordAlgebra(:Quaternions) == CliffordAlgebra(:ℍ)
        @test CliffordAlgebra(:Quaternions) == CliffordAlgebra(:ℍ)
        @test CliffordAlgebra(:Cl2) == CliffordAlgebra(2)
        @test CliffordAlgebra(:Cl3) == CliffordAlgebra(3)
        @test CliffordAlgebra(:Spacetime) == CliffordAlgebra(:STA)
        @test CliffordAlgebra(:PGA2D) == CliffordAlgebra(:Projective2D) == CliffordAlgebra(:Plane2D)
        @test CliffordAlgebra(:PGA3D) == CliffordAlgebra(:Projective3D) == CliffordAlgebra(:Plane3D)
        @test CliffordAlgebra(:CGA2D) == CliffordAlgebra(:Conformal2D)
        @test CliffordAlgebra(:CGA3D) == CliffordAlgebra(:Conformal3D)
        @test CliffordAlgebra(:DCGA3D) == CliffordAlgebra(:DoubleConformal3D)
        @test CliffordAlgebra(:TCGA3D) == CliffordAlgebra(:TripleConformal3D)
        @test CliffordAlgebra(:DCGSTA) == CliffordAlgebra(:DoubleConformalSpacetime)
        @test CliffordAlgebra(:QCGA) == CliffordAlgebra(:QuadricConformal)

        @test_throws ArgumentError CliffordAlgebra(:UnknownAlgebra)

        cl2 = CliffordAlgebra(2)

        @test basesymbol(cl2, 1) == :𝟏
        @test basesymbol(cl2, 2) == :e1
        @test basesymbol(cl2, 3) == :e2
        @test basesymbol(cl2, 4) == :e1e2

        @test character(cl2) == -1
        @test dimension(cl2) == 4
        @test order(cl2) == 2
        @test signature(cl2) == (2,0,0)

        io = IOBuffer()
        signaturetable(io, cl2)
        @test String(take!(io)) == "┌────┬────┐\n│ e1 │ +1 │\n│ e2 │ +1 │\n└────┴────┘\n"
        cayleytable(io, cl2)
        @test String(take!(io)) == "┌───────┬──────────────┬───────┐\n│  +1   │  +e1    +e2  │ +e1e2 │\n├───────┼──────────────┼───────┤\n│  +e1  │  +1    +e1e2 │  +e2  │\n│  +e2  │ -e1e2   +1   │  -e1  │\n├───────┼──────────────┼───────┤\n│ +e1e2 │  -e2    +e1  │  -1   │\n└───────┴──────────────┴───────┘\n"
    end

    @testset "MultiVector" begin
        cl3 = CliffordAlgebra(:Cl3)

        @test scalar(MultiVector(cl3,1)) == 1
        @test MultiVector(cl3, (1,2,3)).e1 == 1
        @test MultiVector(cl3, (1,2,3)).e2 == 2
        @test MultiVector(cl3, (1,2,3)).e3 == 3
        @test MultiVector(cl3, (1,2,3)).𝟏 == 0
        @test MultiVector(cl3, (1,2,3)).e1e2 == 0
        @test MultiVector(cl3, (1,2,3)).e2e3 == 0
        @test MultiVector(cl3, (1,2,3)).e3e1 == 0
        @test MultiVector(cl3, (1,2,3)).e1e2e3 == 0

        @test_throws AssertionError MultiVector(cl3, (1,2))

        mv = 1 + cl3.e1 - cl3.e2 + cl3.e3 - cl3.e1e2 + cl3.e2e3 - cl3.e3e1 + cl3.e1e2e3

        @test MultiVector(cl3, 0) == zero(mv)
        @test MultiVector(cl3, 1) == one(mv)

        @test iszero(zero(mv))
        @test isone(one(mv))

        @test algebra(mv) == cl3

        @test eltype(mv) == typeof(scalar(mv))

        @test convert(Float64, one(mv)) == one(Float64)
        @test_throws InexactError convert(Float64, mv)

        mv2 = cl3.e1 + 0 * cl3.e2

        @test prune(mv) == mv
        @test typeof(prune(mv2)) != typeof(mv2)

        @test isgrade(mv2, 1)
        @test !isgrade(mv2, 0)
        @test !isgrade(mv2, 2)
        @test !isgrade(mv2, 3)

        @test !isgrade(mv, 0)
        @test !isgrade(mv, 1)
        @test !isgrade(mv, 2)
        @test !isgrade(mv, 3)

        @test isgrade(grade(mv,2),2)

        @test even(mv2) == zero(mv2)
        @test odd(mv2) == mv2

        @test grin(mv2) == -mv2

        @test dual(mv2) == cl3.e2e3

        @test reverse(cl3.𝟏) == cl3.𝟏
        @test reverse(cl3.e2e3) == -cl3.e2e3
        @test reverse(cl3.e1e2e3) == -cl3.e1e2e3

        @test reverse(mv) == ~mv

        @test pseudoscalar(cl3) == cl3.e1e2e3

        mva = cl3.e1 - cl3.e1e2
        mvb = 1 + cl3.e2 + cl3.e1e2e3

        @assert matrix(mva) * matrix(mvb) == matrix(mva*mvb)
        @assert matrix(mva) * vector(mvb) == vector(mva*mvb)
    end

    @testset "arithmetic for $cl" for cl in (
            CliffordAlgebra(3,0,0),
            CliffordAlgebra(2,1,0),
            CliffordAlgebra(1,2,0),
            CliffordAlgebra(2,0,1),
            CliffordAlgebra(1,0,2),
            CliffordAlgebra(0,2,1),
            CliffordAlgebra(0,1,2),
            CliffordAlgebra(1,1,1)
        )

        cl = CliffordAlgebra(1,1,1)
        e1 = cl.e1
        e2 = cl.e2
        e3 = cl.e3
        e12 = cl.e1e2
        e23 = cl.e2e3
        e31 = cl.e3e1
        e123 = cl.e1e2e3

        s1 = scalar(e1*e1)
        s2 = scalar(e2*e2)
        s3 = scalar(e3*e3)

        @testset "addition/subtraction" begin
            @test e1 + e2 == e2 + e1
            @test e1 - e2 == -(e2 - e1)
            @test e3 + 1 == 1 + e3
        end

        @testset "geometric product" begin
            @test e1 * e1 == s1
            @test e2 * e2 == s2
            @test e3 * e3 == s3

            @test e1 * e2 == e12
            @test e2 * e3 == e23
            @test e1 * e3 == -e31
            @test e1 * e2 * e3 == e123
        end

        @testset "exterior product" begin
            @test e1 ∧ e1 == 0
            @test e2 ∧ e2 == 0
            @test e3 ∧ e3 == 0

            @test e1 ∧ e2 == e12
            @test e2 ∧ e3 == e23
            @test e3 ∧ e1 == e31

            @test e12 ∧ e31 == 0
            @test e23 ∧ e31 == 0
            @test e12 ∧ e23 == 0

            @test e1 ∧ e23 == e123
            @test e12 ∧ e3 == e123

            @test 2 ∧ e1 == MultiVector(cl,2) ∧ e1
            @test 2 ∧ 3 == MultiVector(cl,2) ∧ MultiVector(cl,3)
        end

        @testset "fatdot product" begin
            @test e1 ⋅ e1 == s1
            @test e2 ⋅ e2 == s2
            @test e3 ⋅ e3 == s3

            @test e1 ⋅ e12 == s1*e2
            @test e1 ⋅ e23 == 0
            @test e1 ⋅ e31 == -s1*e3

            @test e123 ⋅ e1 == s1*e23
            @test e123 ⋅ e2 == s2*e31
            @test e123 ⋅ e3 == s3*e12
            @test e123 ⋅ e12 == -s1*s2*e3
            @test e123 ⋅ e23 == -s2*s3*e1
            @test e123 ⋅ e31 == -s3*s1*e2
            @test e123 ⋅ e123 == -s1*s2*s3

            @test 2 ⋅ e1 == MultiVector(cl,2) ⋅ e1
            @test 2 ⋅ 3 == MultiVector(cl,2) ⋅ MultiVector(cl,3)
        end

        @testset "left contraction product" begin
            @test e1 ⨼ e1 == s1
            @test e2 ⨼ e2 == s2
            @test e3 ⨼ e3 == s3

            @test e1 ⨼ e12 == s1*e2
            @test e1 ⨼ e23 == 0
            @test e1 ⨼ e31 == -s1*e3

            @test e123 ⨼ e1 == 0
            @test e123 ⨼ e2 == 0
            @test e123 ⨼ e3 == 0

            @test e123 ⨼ e12 == 0
            @test e123 ⨼ e23 == 0
            @test e123 ⨼ e31 == 0

            @test e123 ⨼ e123 == 0

            @test 2 ⨼ e1 == MultiVector(cl,2) ⨼ e1
            @test 2 ⨼ 3 == MultiVector(cl,2) ⨼ MultiVector(cl,3)
        end

        @testset "right contraction product" begin
            @test e1 ⨽ e1 == s1
            @test e2 ⨽ e2 == s2
            @test e3 ⨽ e3 == s3

            @test e1 ⨽ e12 == 0
            @test e1 ⨽ e23 == 0
            @test e1 ⨽ e31 == 0

            @test e123 ⨽ e1 == s1*e23
            @test e123 ⨽ e2 == s2*e31
            @test e123 ⨽ e3 == s3*e12

            @test e123 ⨽ e12 == -s1*s2*e3
            @test e123 ⨽ e23 == -s2*s3*e1
            @test e123 ⨽ e31 == -s3*s1*e2

            @test e123 ⨽ e123 == 0

            @test 2 ⨽ e1 == MultiVector(cl,2) ⨽ e1
            @test 2 ⨽ 3 == MultiVector(cl,2) ⨽ MultiVector(cl,3)
        end

        @testset "scalar product" begin
            @test e1 ⋆ e1 == s1
            @test e2 ⋆ e2 == s2
            @test e3 ⋆ e3 == s3

            @test e1 ⋆ e12 == 0
            @test e1 ⋆ e23 == 0
            @test e1 ⋆ e31 == 0

            @test e123 ⋆ e1 == 0
            @test e123 ⋆ e2 == 0
            @test e123 ⋆ e3 == 0

            @test e123 ⋆ e123 == s1 * s2 * s3

            @test 2 ⋆ e1 == MultiVector(cl,2) ⋆ e1
            @test 2 ⋆ 3 == MultiVector(cl,2) ⋆ MultiVector(cl,3)
        end

        @testset "regressive product" begin
            @test e1 ∨ e1 == 0
            @test e2 ∨ e2 == 0
            @test e3 ∨ e3 == 0

            @test e1 ∨ e23 == 1
            @test e2 ∨ e31 == 1
            @test e3 ∨ e12 == 1

            @test e123 ∨ e1 == e1
            @test e123 ∨ e2 == e2
            @test e123 ∨ e3 == e3

            @test e123 ∨ e12 == e12
            @test e123 ∨ e23 == e23
            @test e123 ∨ e31 == e31

            @test e123 ∨ e123 == e123

            @test 2 ∨ e1 == MultiVector(cl,2) ∨ e1
            @test 2 ∨ 3 == MultiVector(cl,2) ∨ MultiVector(cl,3)
        end

        @testset "commutator product" begin
            @test 1 ×₋ 2 == 0
            @test 2 ×₋ e1 == 0
            @test 2 ×₋ e2 == 0
            @test 2 ×₋ e3 == 0
            @test 2 ×₋ e12 == 0
            @test 2 ×₋ e23 == 0
            @test 2 ×₋ e31 == 0
            @test 2 ×₋ e123 == 0

            @test e1 ×₋ e1 == 0
            @test e2 ×₋ e2 == 0
            @test e3 ×₋ e3 == 0

            @test e1 ×₋ e2 == e12
            @test e1 ×₋ e3 == -e31
            @test e2 ×₋ e3 == e23

            @test e12 ×₋ e12 == 0
            @test e23 ×₋ e23 == 0
            @test e31 ×₋ e31 == 0

            @test e23 ×₋ e12 == s2 * e31
            @test e23 ×₋ e31 == -s3 * e12
            @test e12 ×₋ e31 == s1 * e23

            @test e123 ×₋ e123 == 0
        end

        @testset "anti-commutator product" begin
            @test 1 ×₊ 2 == 2
            @test 2 ×₊ e1 == 2*e1
            @test 2 ×₊ e2 == 2*e2
            @test 2 ×₊ e3 == 2*e3
            @test 2 ×₊ e12 == 2*e12
            @test 2 ×₊ e23 == 2*e23
            @test 2 ×₊ e31 == 2*e31
            @test 2 ×₊ e123 == 2*e123

            @test e1 ×₊ e1 == s1
            @test e2 ×₊ e2 == s2
            @test e3 ×₊ e3 == s3

            @test e1 ×₊ e2 == 0
            @test e1 ×₊ e3 == 0
            @test e2 ×₊ e3 == 0

            @test e12 ×₊ e12 == -s1 * s2
            @test e23 ×₊ e23 == -s2 * s3
            @test e31 ×₊ e31 == -s3 * s1

            @test e23 ×₊ e12 == 0
            @test e23 ×₊ e31 == 0
            @test e12 ×₊ e31 == 0

            @test e123 ×₊ e123 == -s1 * s2 * s3
        end

        @testset "sandwich product" for n = 1:100
            r = -100:100
            mva = rand(r) + rand(r) * e1 +
                  rand(r) * e2 + rand(r) * e3 +
                  rand(r) * e12 + rand(r) * e23 +
                  rand(r) * e31 + rand(r) * e123

            mvb = rand(r) + rand(r) * e1 +
                  rand(r) * e2 + rand(r) * e3 +
                  rand(r) * e12 + rand(r) * e23 +
                  rand(r) * e31 + rand(r) * e123

            @test mva ≀ mvb == mva * mvb * reverse(mva)
        end

        @testset "misc functions" begin
            mv = 1 + e1 - e2 + e3 + e12 - e23 + e31 - e123
            @test polarize(mv)*pseudoscalar(algebra(mv)) == character(algebra(mv))
            @test Λᵏ(mv,0) == grade(mv,0)
            @test Λᵏ(mv,1) == grade(mv,1)
            @test Λᵏ(mv,2) == grade(mv,2)
            @test Λᵏ(mv,3) == grade(mv,3)

            @test norm_sqr(e1) == s1
            @test norm_sqr(2*e2) == 4*s2
            @test norm_sqr(3*e3) == 9*s3
            @test norm_sqr(e12) == s1*s2
            @test norm_sqr(e23) == s2*s3
            @test norm_sqr(e31) == s3*s1
            @test norm_sqr(e123) == s1*s2*s3

            @test norm(2*e1) == 2*s1
            @test norm(-2*e1) == 2*s1

            if !iszero(s1)
                @test inv(e1) * e1 == e1 * inv(e1) == 1
            else
                @test_throws SingularException inv(e1)
                @test_throws SingularException inv(e12)
                @test_throws SingularException inv(e31)
                @test_throws SingularException inv(e123)
            end

            if !iszero(s2)
                @test inv(e2) * e2 == e2 * inv(e2) == 1
            else
                @test_throws SingularException inv(e2)
                @test_throws SingularException inv(e12)
                @test_throws SingularException inv(e23)
                @test_throws SingularException inv(e123)
            end

            if !iszero(s3)
                @test inv(e3) * e3 == e3 * inv(e3) == 1
            else
                @test_throws SingularException inv(e3)
                @test_throws SingularException inv(e31)
                @test_throws SingularException inv(e23)
                @test_throws SingularException inv(e123)
            end

            if !iszero(s1) && !iszero(s2)
                @test inv(e12) * e12 == e12 * inv(e12) == 1
            end

            if !iszero(s1) && !iszero(s3)
                @test inv(e31) * e31 == e31 * inv(e31) == 1
            end

            if !iszero(s2) && !iszero(s3)
                @test inv(e23) * e23 == e23 * inv(e23) == 1
            end

            if !iszero(s1) && !iszero(s2) && !iszero(s3)
                @test (e1+e2+e3) \ (e1+e2+e3) == 1
                @test (e12 + e23 + e31) / (e12 + e23 + e31) == 1
                @test inv(e123) * e123 == 1
            end

            for mv in (e1, e2, e3, e12, e23, e31, e123)
                if scalar(mv*mv) < 0
                    @test exp(mv) == cos(1) + sin(1) * mv
                elseif scalar(mv*mv) > 0
                    @test exp(mv) == cosh(1) + sinh(1) * mv
                else
                    @test exp(mv) == 1 + mv
                end
                @test convert(Float32,exp(mv) * exp(-mv)) ≈ 1
            end

            @test exp(1 + e1 + e12 + e123) == exp(1 + e123) * exp(e1 + e12)

            mv = 1 + e1 - e2 - e3 + e12 - e23 + e31 - e123

            M = [1 0 0; 0 1 0 ; 0 0 1]
            @test outermorphism(M, mv) == mv

            M = [0 0 0; 0 0 0 ; 0 0 0]
            @test outermorphism(M, mv) == 1

            M = [0 1 0; 1 0 0 ; 0 0 1]
            @test outermorphism(M, mv) == 1 + e2 - e1 - e3 - e12 + e31 - e23 + e123

            M = [0 0 1; 0 1 0 ; 1 0 0]
            @test outermorphism(M, mv) == 1 + e3 - e2 - e1 - e23 + e12 - e31 + e123

            M = [0 0 1; 0 1 0 ; 0 0 0]
            @test outermorphism(M, mv) == 1 - e2 - e1 + e12
        end

    end

end

max_stars_repo_name: ATell-SoundTheory/CliffordAlgebras

max_issues_repo_name: ATell-SoundTheory/CliffordAlgebras


text:
using SparseArrays

@testset "base" begin

@testset "helper functions" begin
    @testset "findallmax" begin
        @test findallmax([-Inf, -Inf, -Inf]) == (-Inf, [1, 2, 3])
        @test findallmax([Inf, Inf, Inf]) == (Inf, [1, 2, 3])
        @test findallmax([0,1,2,1,2,1,0]) == (2, [3, 5])
        @test begin
            max_val, inds = findallmax([NaN, NaN, NaN])
            isnan(max_val) && inds == [1, 2, 3]
        end
    end
    @testset "onehot" begin
        @test onehot(5, 3) == [0, 0, 1, 0, 0]
        @test onehot(5, 3, Float64) == [0., 0., 1., 0., 0.]
        @test onehot(5, 3, Float64; isdense=false) == SparseVector([0., 0., 1., 0., 0.])
    end
    @testset "deletefirst!" begin
        @test deletefirst!([1,2,3,2,1], 1) == [2,3,2,1]
        @test deletefirst!([1,2,3,2,1], 2) == [1,3,2,1]
        @test deletefirst!([1,2,3,2,1], 3) == [1,2,2,1]
    end
    @testset "discounted reward" begin
        @test discounted_reward([2., 3., 4.], 0.9) ≈ 2. + 0.9 * 3. + 0.9 * 0.9 * 4
    end
    @testset "reverse_discounted_rewards" begin
        rewards, γ = [2., 3., 4.], 0.9
        @test collect(reverse_discounted_rewards(rewards, γ)) == reverse([discounted_reward(rewards, γ), discounted_reward(rewards[2:end], γ), discounted_reward(rewards[3:end], γ)])
    end
    @testset "importance_weight" begin
        actions = [1, 2]
        π_target = DeterministicPolicy([1,2,1,2], length(actions))
        π_behavior = RandomPolicy([0.5, 0.5])

        @test importance_weight(π_target, π_behavior, [1, 2, 3, 4], [1, 2, 1, 2]) == 16.
        @test importance_weight(π_target, π_behavior, [1, 2, 3, 4], [1, 1, 1, 1]) == 0.
    end
    @testset "reverse_importance_weights" begin
        actions = [1, 2]
        π_target = DeterministicPolicy([1,2,1,2], length(actions))
        π_behavior = RandomPolicy([0.5, 0.5])

        states, actions = [1, 2, 3, 4], [2, 2, 2, 2]
        @test collect(reverse_importance_weights(π_target, π_behavior, states, actions)) == [1., 2., 0., 0., 0.]
    end
    @testset "huber_loss" begin
        @testset "all linear delta" begin
            δ = 0.5
            predictions = [1.5, -1.4, -1.0, 0.0]
            labels = [0.0, 1.0, 0.0, 1.5]
            expected = mean(δ .* [1.5, 2.4, 1.0, 1.5]) .- 0.5 .* δ ^2
            @test expected ≈ huber_loss(labels, predictions;δ=δ)
        end
        @testset "all quadratic delta" begin
            δ = 0.5
            predictions = [1.5, -1.4, -0.5, 0.0]
            labels = [1.0, -1.0, 0.0, 0.5]
            expected = mean(0.5 .* [0.5, 0.4, 0.5, 0.5] .^2)
            @test expected ≈ huber_loss(labels, predictions;δ=δ)
        end
    end
end

@testset "iterators" begin
    @testset "reductions" begin
        @test reductions(+, 1:4) == [1, 3, 6, 10]
        @test collect(reductions(+, 1:4; init=10)) == [10, 11, 13, 16, 20]
        @test collect(reductions(+, [1 3; 2 4]; init=10)) == [10, 11, 13, 16, 20]
    end

    @testset "is_first_visit" begin
        @test is_first_visit([1,2,1,3,1,2,3]) ==  [true, true, false, true, false, false, false]
        @test collect(Iterators.reverse(IsFirstVisit([1,2,1,3,1,2,3]))) == [false, false, false, true, false, true, true]
    end
end

@testset "decays" begin
    @testset "incrementer" begin
        x = incrementer()
        for i in 1:10
            @test x() == i
        end
    end
    @testset "multiplier" begin
        x = multiplier(ratio=0.99)
        for i in 1:10
            @test x() ≈ 0.99 ^ (i - 1)
        end
    end
    @testset "inverse_decay" begin
        x = inverse_decay()
        for i in 1:10
            @test x() == 1 / i
        end
    end
    @testset "cached_inverse_decay" begin
        x = cached_inverse_decay()
        for i in 1:3
            for j in 1:3
                @test x(j) == 1 / i
            end
        end
        @test x(3) == 1 / 4
        @test x(4) == 1
    end
    @testset "exp_decay" begin
        init, λ, decay_step, clip = 1.0, 0.5, 2, 0.1
        f = exp_decay(init=init, λ=λ, decay_step=decay_step, clip=clip)
        f_real(i) = max(init * exp(-λ*i/decay_step), clip)
        for i in 0:10
            @test f() == f_real(i)
        end
    end
end

@testset "callbacks" begin
    @testset "stop_at_step" begin
        env = SimpleRandomWalkEnv()
        agent = Agent(MonteCarloLearner(TabularQ(length(observationspace(env)), length(actionspace(env))), EpsilonGreedySelector(0.01)), EpisodeSARDBuffer())
        cb = stop_at_step(2, false)

        @test cb(env, agent) == false
        @test cb() == 1
        @test cb(env, agent) == true
        @test cb() == 2
    end

    @testset "stop_at_episode" begin
        env = SimpleRandomWalkEnv()
        agent = Agent(MonteCarloLearner(TabularQ(length(observationspace(env)), length(actionspace(env))), EpsilonGreedySelector(0.01)), EpisodeSARDBuffer())
        cb = stop_at_episode(2, false)

        reset!(env)
        @test isend(env) == false
        @test cb(env, agent) == false

        env.state = env.N
        @test isend(env) == true
        @test cb(env, agent) == false  # the first episode

        reset!(env)
        @test isend(env) == false
        @test cb(env, agent) == false

        env.state = env.N
        @test isend(env) == true
        @test cb(env, agent) == true  # the second episode
    end

    @testset "stop_when_done" begin
        env = SimpleRandomWalkEnv()
        agent = Agent(MonteCarloLearner(TabularQ(length(observationspace(env)), length(actionspace(env))), EpsilonGreedySelector(0.01)), EpisodeSARDBuffer())
        cb = stop_when_done

        reset!(env)
        @test isend(env) == false
        @test cb(env, agent) == false

        env.state = env.N
        @test isend(env) == true
        @test cb(env, agent) == true
    end

    @testset "steps_per_episode" begin
        env = SimpleRandomWalkEnv()
        agent = Agent(MonteCarloLearner(TabularQ(length(observationspace(env)), length(actionspace(env))), EpsilonGreedySelector(0.01)), EpisodeSARDBuffer())
        cb = steps_per_episode()
        steps = [2,3,4]

        for n in steps
            reset!(env)
            for i in 1:(n-1)
                cb(env, agent)
            end
            env.state = env.N
            cb(env, agent)
        end
        @test cb() == steps
    end

    @testset "rewards_of_each_episode" begin
        env = SimpleRandomWalkEnv()
        agent = Agent(MonteCarloLearner(TabularQ(length(observationspace(env)), length(actionspace(env))), EpsilonGreedySelector(0.01)), EpisodeSARDBuffer())
        cb = rewards_of_each_episode()
        rewards = [[1., 2., 3.], [4., 5.], [6.]]

        for reward in rewards
            reset!(env)
            @test isend(env) == false
            push!(buffer(agent), 1, 1)
            for r in reward[1:end-1]
                push!(buffer(agent), r, false, 1, 1)
                cb(env, agent)
            end
            env.state = env.N
            @test isend(env) == true
            push!(buffer(agent), reward[end], true, 1, 1)
            cb(env, agent)
        end
        @test cb() == map(sum, rewards)
    end
end

end


text:
# Example 2.8 of
# P. Parrilo and A. Jadbabaie
# Approximation of the joint spectral radius using sum of squares
# Linear Algebra and its Applications, Elsevier, 2008, 428, 2385-2402
# The JSR was conjectured to be 8.9149 = √ρ(A_1 * A_3)

using LinearAlgebra

@testset "[PJ08] Example 5.4" begin
    # values with log-accuracy 4e-7 taken from the examples/PJ08e54.ipynb
    expected_lb = [5.635326998733677, 6.777596245727604, 7.423337986847027]
    expected_ub = [9.760675006197351, 8.91982041593713,  8.914964296278484]
    A1 = [ 0  1  7  4;
           1  6 -2 -3;
          -1 -1 -2 -6;
           3  0  9  1]
    A2 = [-3  3  0 -2;
          -2  1  4  9;
           4 -3  1  1;
           1 -5 -1 -2]
    A3 = [ 1  4  5 10;
           0  5  1 -4;
           0 -1  4  6;
          -1  5  0  1]
    s = discreteswitchedsystem([A1, A2, A3])
    smp = periodicswitching(s, [1, 3])

    @testset "Gripenberg" begin
        for p in [1, 2, Inf]
            psw, ub = gripenberg(s, max_length = 1, matrix_norm = A -> opnorm(A, p), verbose=0)
            @test psw == periodicswitching(s, [2])
            @test ub ≈ opnorm(A3, p)
        end
        psw, ub = gripenberg(s, verbose=0)
        @test psw == smp
        @test ub ≈ psw.growthrate + 1e-2 # 1e-2 is the default δ
    end

    for optimizer_constructor in sdp_factories
        sosdata(s).lb = 0
        for d in 1:3
            tol = 1e-5
            lb, ub = soslyapb(s, d, optimizer_constructor=optimizer_constructor, tol=tol)
            @test abs(log(expected_lb[d]) - log(lb)) <= tol
            @test abs(log(expected_ub[d]) - log(ub)) <= tol
            seq = sosbuildsequence(s, d, p_0=:Primal)
            psw = findsmp(seq)
            @test psw !== nothing
            @test psw == smp
        end
    end
end


text:
## vanilla hyperthetical TE without failure scenarios
function HYPERTHETICAL(GRB_ENV, edges, capacity, flows, demand, T, Tf)
    printstyled("\n** solving Vanilla TE LP... ", color=:yellow)

    nedges = size(edges,1)
    nflows = size(flows,1)
    ntunnels = size(T,1)

    tunnel_num = 0
    for x in 1:length(Tf)
        if size(Tf[x],1) > tunnel_num
            tunnel_num = size(Tf[x],1)
        end
    end

    #CREATE TUNNEL EDGE MATRIX
    L = zeros(ntunnels, nedges)
    for t in 1:ntunnels
        for e in 1:nedges
            if in(e, T[t])
                L[t,e] = 1
            end
        end
    end

    model = Model(() -> Gurobi.Optimizer(GRB_ENV))
    set_optimizer_attribute(model, "OutputFlag", 0)
    set_optimizer_attribute(model, "Threads", 32)

    model_def_time = @elapsed begin

    @variable(model, b[1:nflows] >= 0)
    @variable(model, a[1:nflows,1:tunnel_num] >= 0)

    for f in 1:nflows
        @constraint(model, sum(a[f,t] for t in 1:size(Tf[f],1)) >= b[f])   #the sum of all allocated bandwidths on every flow must be >= the total bandwidth for that flow
    end

    for e in 1:nedges
        @constraint(model, sum(a[f,t] * L[Tf[f][t],e] for f in 1:nflows for t in 1:size(Tf[f],1)) <= capacity[e])   #overlapping flows cannot add up to the capacity of that link
    end

    for f in 1:nflows
        @constraint(model, b[f] >= 0)
    end

    for f in 1:nflows
        @constraint(model, b[f] <= demand[f])   #all allocated bandwidths must be less than the demand for that flow
        for t in 1:size(Tf[f],1)
            @constraint(model, a[f,t] >= 0)     #each allocated bandwidth on for flow f on tunnel t >= 0
        end
    end

    @objective(model, Max, sum((b[i] for i in 1:size(b,1))))
    end

    optimize!(model)
    solve_runtime = solve_time(model)
    opt_runtime = solve_runtime + model_def_time

    return value.(a), value.(b), solve_runtime, opt_runtime
end

max_stars_repo_name: hipersys-team/arrow

max_issues_repo_name: arrow-bow/arrow

max_forks_repo_name: arrow-bow/arrow