text:
const _allAxes = [:auto, :left, :right]
@compat const _axesAliases = Dict(
    :a => :auto,
    :l => :left,
    :r => :right
  )

const _3dTypes = [:path3d, :scatter3d, :surface, :wireframe]
const _allTypes = vcat([
                        :none, :line, :path, :steppre, :steppost, :sticks, :scatter,
                        :heatmap, :hexbin, :hist, :density, :bar, :hline, :vline, :ohlc,
                        :contour, :pie
                       ], _3dTypes)
@compat const _typeAliases = Dict(
    :n             => :none,
    :no            => :none,
    :l             => :line,
    :p             => :path,
    :stepinv       => :steppre,
    :stepsinv      => :steppre,
    :stepinverted  => :steppre,
    :stepsinverted => :steppre,
    :step          => :steppost,
    :steps         => :steppost,
    :stair         => :steppost,
    :stairs        => :steppost,
    :stem          => :sticks,
    :stems         => :sticks,
    :dots          => :scatter,
    :histogram     => :hist,
    :pdf           => :density,
    :contours      => :contour,
    :line3d        => :path3d,
    :surf          => :surface,
    :wire          => :wireframe,
  )

ishistlike(lt::Symbol) = lt in (:hist, :density)
islinelike(lt::Symbol) = lt in (:line, :path, :steppre, :steppost)
isheatmaplike(lt::Symbol) = lt in (:heatmap, :hexbin)

const _allStyles = [:auto, :solid, :dash, :dot, :dashdot, :dashdotdot]
@compat const _styleAliases = Dict(
    :a    => :auto,
    :s    => :solid,
    :d    => :dash,
    :dd   => :dashdot,
    :ddd  => :dashdotdot,
  )

# const _allMarkers = [:none, :auto, :ellipse, :rect, :diamond, :utriangle, :dtriangle,
#                      :cross, :xcross, :star5, :star8, :hexagon, :octagon, Shape]
const _allMarkers = vcat(:none, :auto, sort(collect(keys(_shapes))))
@compat const _markerAliases = Dict(
    :n            => :none,
    :no           => :none,
    :a            => :auto,
    :circle       => :ellipse,
    :c            => :ellipse,
    :square       => :rect,
    :sq           => :rect,
    :r            => :rect,
    :d            => :diamond,
    :^            => :utriangle,
    :ut           => :utriangle,
    :utri         => :utriangle,
    :uptri        => :utriangle,
    :uptriangle   => :utriangle,
    :v            => :dtriangle,
    :V            => :dtriangle,
    :dt           => :dtriangle,
    :dtri         => :dtriangle,
    :downtri      => :dtriangle,
    :downtriangle => :dtriangle,
    :+            => :cross,
    :plus         => :cross,
    :x            => :xcross,
    :X            => :xcross,
    :star         => :star5,
    :s            => :star5,
    :star1        => :star5,
    :s2           => :star8,
    :star2        => :star8,
    :p            => :pentagon,
    :pent         => :pentagon,
    :h            => :hexagon,
    :hex          => :hexagon,
    :hep          => :heptagon,
    :o            => :octagon,
    :oct          => :octagon,
    :spike        => :vline,
  )

const _allScales = [:identity, :ln, :log2, :log10, :asinh, :sqrt]
@compat const _scaleAliases = Dict(
    :none => :identity,
    :log  => :log10,
  )

# -----------------------------------------------------------------------------

const _seriesDefaults = Dict{Symbol, Any}()

# series-specific
_seriesDefaults[:axis]            = :left
_seriesDefaults[:label]           = "AUTO"
_seriesDefaults[:linetype]        = :path
_seriesDefaults[:linestyle]       = :solid
_seriesDefaults[:linewidth]       = 1
_seriesDefaults[:linecolor]       = :auto
_seriesDefaults[:linealpha]       = nothing
# _seriesDefaults[:linestroke]      = Stroke(1, :auto, nothing, :solid)  # linewidth, linecolor, linealpha, linestyle
# _seriesDefaults[:fillbrush]       = Brush(nothing, :match, nothing)  # fillrange, fillcolor, fillalpha
_seriesDefaults[:fillrange]       = nothing   # ribbons, areas, etc
_seriesDefaults[:fillcolor]       = :match
_seriesDefaults[:fillalpha]       = nothing
_seriesDefaults[:markershape]     = :none
# _seriesDefaults[:markerstroke]    = Stroke(1, :match_foreground, nothing, :solid)
# _seriesDefaults[:markerbrush]     = Brush(6, :match, nothing)
_seriesDefaults[:markercolor]     = :match
_seriesDefaults[:markeralpha]     = nothing
_seriesDefaults[:markersize]      = 6
_seriesDefaults[:markerstrokestyle] = :solid
_seriesDefaults[:markerstrokewidth] = 1
_seriesDefaults[:markerstrokecolor] = :match
_seriesDefaults[:markerstrokealpha] = nothing
# _seriesDefaults[:ribbon]          = nothing
# _seriesDefaults[:ribboncolor]     = :match
_seriesDefaults[:nbins]           = 30               # number of bins for heatmaps and hists
_seriesDefaults[:smooth]          = false               # regression line?
_seriesDefaults[:group]           = nothing           # groupby vector
# _seriesDefaults[:annotation]      = nothing           # annotation tuple(s)... (x,y,annotation)
_seriesDefaults[:x]               = nothing
_seriesDefaults[:y]               = nothing
_seriesDefaults[:z]               = nothing           # depth for contour, surface, etc
_seriesDefaults[:zcolor]          = nothing           # value for color scale
# _seriesDefaults[:surface]         = nothing
# _seriesDefaults[:nlevels]         = 15
_seriesDefaults[:levels]          = 15
_seriesDefaults[:orientation]     = :vertical

const _plotDefaults = Dict{Symbol, Any}()

# plot globals
_plotDefaults[:title]             = ""
_plotDefaults[:xlabel]            = ""
_plotDefaults[:ylabel]            = ""
_plotDefaults[:yrightlabel]       = ""
_plotDefaults[:legend]            = true
_plotDefaults[:background_color]  = colorant"white"
_plotDefaults[:foreground_color]  = :auto
_plotDefaults[:xlims]             = :auto
_plotDefaults[:ylims]             = :auto
_plotDefaults[:xticks]            = :auto
_plotDefaults[:yticks]            = :auto
_plotDefaults[:xscale]            = :identity
_plotDefaults[:yscale]            = :identity
_plotDefaults[:xflip]             = false
_plotDefaults[:yflip]             = false
_plotDefaults[:size]              = (600,400)
_plotDefaults[:pos]               = (0,0)
_plotDefaults[:windowtitle]       = "Plots.jl"
_plotDefaults[:show]              = false
_plotDefaults[:layout]            = nothing
_plotDefaults[:n]                 = -1
_plotDefaults[:nr]                = -1
_plotDefaults[:nc]                = -1
_plotDefaults[:color_palette]     = :auto
_plotDefaults[:link]              = false
_plotDefaults[:linkx]             = false
_plotDefaults[:linky]             = false
_plotDefaults[:linkfunc]          = nothing
_plotDefaults[:tickfont]          = font(8)
_plotDefaults[:guidefont]         = font(11)
_plotDefaults[:legendfont]        = font(8)
_plotDefaults[:grid]              = true
_plotDefaults[:annotation]      = nothing           # annotation tuple(s)... (x,y,annotation)

# TODO: x/y scales

const _allArgs = sort(collect(union(keys(_seriesDefaults), keys(_plotDefaults))))
supportedArgs(::PlottingPackage) = _allArgs
supportedArgs() = supportedArgs(backend())

# -----------------------------------------------------------------------------

makeplural(s::Symbol) = symbol(string(s,"s"))

autopick(arr::AVec, idx::Integer) = arr[mod1(idx,length(arr))]
autopick(notarr, idx::Integer) = notarr

autopick_ignore_none_auto(arr::AVec, idx::Integer) = autopick(setdiff(arr, [:none, :auto]), idx)
autopick_ignore_none_auto(notarr, idx::Integer) = notarr

function aliasesAndAutopick(d::Dict, sym::Symbol, aliases::Dict, options::AVec, plotIndex::Int)
  if d[sym] == :auto
    d[sym] = autopick_ignore_none_auto(options, plotIndex)
  elseif haskey(aliases, d[sym])
    d[sym] = aliases[d[sym]]
  end
end

function aliases(aliasMap::Dict, val)
  # sort(vcat(val, collect(keys(filter((k,v)-> v==val, aliasMap)))))
  sortedkeys(filter((k,v)-> v==val, aliasMap))
end

# -----------------------------------------------------------------------------

# Alternate args

@compat const _keyAliases = Dict(
    :c                  => :linecolor,
    :color              => :linecolor,
    :colour             => :linecolor,
    :lab                => :label,
    :l                  => :line,
    :w                  => :linewidth,
    :width              => :linewidth,
    :lw                 => :linewidth,
    :la                 => :linealpha,
    :lineopacity        => :linealpha,
    :type               => :linetype,
    :lt                 => :linetype,
    :t                  => :linetype,
    :style              => :linestyle,
    :s                  => :linestyle,
    :ls                 => :linestyle,
    :m                  => :marker,
    :mark               => :marker,
    :shape              => :markershape,
    :mc                 => :markercolor,
    :mcolor             => :markercolor,
    :markercolour       => :markercolor,
    :ms                 => :markersize,
    :msize              => :markersize,
    :ma                 => :markeralpha,
    :alpha              => :markeralpha,
    :opacity            => :markeralpha,
    :markeropacity      => :markeralpha,
    :f                  => :fill,
    :area               => :fill,
    :fillrng            => :fillrange,
    :fc                 => :fillcolor,
    :fcolor             => :fillcolor,
    :fillcolour         => :fillcolor,
    :fa                 => :fillalpha,
    :fillopacity        => :fillalpha,
    :g                  => :group,
    :nb                 => :nbins,
    :nbin               => :nbins,
    :rib                => :ribbon,
    :ann                => :annotation,
    :anns               => :annotation,
    :annotate           => :annotation,
    :annotations        => :annotation,
    :xlab               => :xlabel,
    :ylab               => :ylabel,
    :yrlab              => :yrightlabel,
    :ylabr              => :yrightlabel,
    :y2lab              => :yrightlabel,
    :ylab2              => :yrightlabel,
    :ylabelright        => :yrightlabel,
    :ylabel2            => :yrightlabel,
    :y2label            => :yrightlabel,
    :leg                => :legend,
    :bg                 => :background_color,
    :bgcolor            => :background_color,
    :bg_color           => :background_color,
    :background         => :background_color,
    :background_colour  => :background_color,
    :fg                 => :foreground_color,
    :fgcolor            => :foreground_color,
    :fg_color           => :foreground_color,
    :foreground         => :foreground_color,
    :foreground_colour  => :foreground_color,
    :regression         => :smooth,
    :reg                => :smooth,
    :nlevels            => :levels,
    :nlev               => :levels,
    :levs               => :levels,
    :xlim               => :xlims,
    :xlimit             => :xlims,
    :xlimits            => :xlims,
    :ylim               => :ylims,
    :ylimit             => :ylims,
    :ylimits            => :ylims,
    :xtick              => :xticks,
    :ytick              => :yticks,
    :windowsize         => :size,
    :wsize              => :size,
    :wtitle             => :windowtitle,
    :gui                => :show,
    :display            => :show,
    :palette            => :color_palette,
    :xlink              => :linkx,
    :ylink              => :linky,
  )

# add all pluralized forms to the _keyAliases dict
for arg in keys(_seriesDefaults)
  _keyAliases[makeplural(arg)] = arg
end

# -----------------------------------------------------------------------------

# update the defaults globally

"""
`default(key)` returns the current default value for that key
`default(key, value)` sets the current default value for that key
`default(; kw...)` will set the current default value for each key/value pair
"""

function default(k::Symbol)
  k = get(_keyAliases, k, k)
  if haskey(_seriesDefaults, k)
    return _seriesDefaults[k]
  elseif haskey(_plotDefaults, k)
    return _plotDefaults[k]
  else
    error("Unknown key: ", k)
  end
end

function default(k::Symbol, v)
  k = get(_keyAliases, k, k)
  if haskey(_seriesDefaults, k)
    _seriesDefaults[k] = v
  elseif haskey(_plotDefaults, k)
    _plotDefaults[k] = v
  else
    error("Unknown key: ", k)
  end
end

function default(; kw...)
  for (k,v) in kw
    default(k, v)
  end
end

# -----------------------------------------------------------------------------

wraptuple(x::@compat(Tuple)) = x
wraptuple(x) = (x,)

trueOrAllTrue(f::Function, x::AbstractArray) = all(f, x)
trueOrAllTrue(f::Function, x) = f(x)

function handleColors!(d::Dict, arg, csym::Symbol)
  try
    if arg == :auto
      d[csym] = :auto
    else
      c = colorscheme(arg)
      d[csym] = c
    end
    return true
  end
  false
end

# given one value (:log, or :flip, or (-1,1), etc), set the appropriate arg
# TODO: use trueOrAllTrue for subplots which can pass vectors for these
function processAxisArg(d::Dict, axisletter::@compat(AbstractString), arg)
  T = typeof(arg)
  # if T <: Symbol

  arg = get(_scaleAliases, arg, arg)

  if arg in _allScales
    d[symbol(axisletter * "scale")] = arg

  elseif arg in (:flip, :invert, :inverted)
    d[symbol(axisletter * "flip")] = true

  elseif T <: @compat(AbstractString)
    d[symbol(axisletter * "label")] = arg

  # xlims/ylims
  elseif (T <: Tuple || T <: AVec) && length(arg) == 2
    d[symbol(axisletter * "lims")] = arg

  # xticks/yticks
  elseif T <: AVec
    d[symbol(axisletter * "ticks")] = arg

  elseif arg == nothing
    d[symbol(axisletter * "ticks")] = []

  else
    warn("Skipped $(axisletter)axis arg $arg")

  end
end

function processLineArg(d::Dict, arg)

  # linetype
  if trueOrAllTrue(a -> get(_typeAliases, a, a) in _allTypes, arg)
    d[:linetype] = arg

  # linestyle
  elseif trueOrAllTrue(a -> get(_styleAliases, a, a) in _allStyles, arg)
    d[:linestyle] = arg

  elseif typeof(arg) <: Stroke
    arg.width == nothing || (d[:linewidth] = arg.width)
    arg.color == nothing || (d[:linecolor] = arg.color == :auto ? :auto : colorscheme(arg.color))
    arg.alpha == nothing || (d[:linealpha] = arg.alpha)
    arg.style == nothing || (d[:linestyle] = arg.style)

  elseif typeof(arg) <: Brush
    arg.size  == nothing || (d[:fillrange] = arg.size)
    arg.color == nothing || (d[:fillcolor] = arg.color == :auto ? :auto : colorscheme(arg.color))
    arg.alpha == nothing || (d[:fillalpha] = arg.alpha)

  # linealpha
  elseif trueOrAllTrue(a -> typeof(a) <: Real && a > 0 && a < 1, arg)
    d[:linealpha] = arg

  # linewidth
  elseif trueOrAllTrue(a -> typeof(a) <: Real, arg)
    d[:linewidth] = arg

  # color
  elseif !handleColors!(d, arg, :linecolor)
    warn("Skipped line arg $arg.")

  end
end

function processMarkerArg(d::Dict, arg)

  # markershape
  if trueOrAllTrue(a -> get(_markerAliases, a, a) in _allMarkers, arg)
    d[:markershape] = arg
  elseif trueOrAllTrue(a -> isa(a, Shape), arg)
    d[:markershape] = arg

  # stroke style
  elseif trueOrAllTrue(a -> get(_styleAliases, a, a) in _allStyles, arg)
    d[:markerstrokestyle] = arg

  elseif typeof(arg) <: Stroke
    arg.width == nothing || (d[:markerstrokewidth] = arg.width)
    arg.color == nothing || (d[:markerstrokecolor] = arg.color == :auto ? :auto : colorscheme(arg.color))
    arg.alpha == nothing || (d[:markerstrokealpha] = arg.alpha)
    arg.style == nothing || (d[:markerstrokestyle] = arg.style)

  elseif typeof(arg) <: Brush
    arg.size  == nothing || (d[:markersize]  = arg.size)
    arg.color == nothing || (d[:markercolor] = arg.color == :auto ? :auto : colorscheme(arg.color))
    arg.alpha == nothing || (d[:markeralpha] = arg.alpha)

  # linealpha
  elseif trueOrAllTrue(a -> typeof(a) <: Real && a > 0 && a < 1, arg)
    d[:markeralpha] = arg

  # markersize
  elseif trueOrAllTrue(a -> typeof(a) <: Real, arg)
    d[:markersize] = arg

  # markercolor
  elseif !handleColors!(d, arg, :markercolor)
    warn("Skipped marker arg $arg.")

  end
end

function processFillArg(d::Dict, arg)

  if typeof(arg) <: Brush
    arg.size  == nothing || (d[:fillrange] = arg.size)
    arg.color == nothing || (d[:fillcolor] = arg.color == :auto ? :auto : colorscheme(arg.color))
    arg.alpha == nothing || (d[:fillalpha] = arg.alpha)

  elseif !handleColors!(d, arg, :fillcolor)
    d[:fillrange] = arg
  end
end

"Handle all preprocessing of args... break out colors/sizes/etc and replace aliases."
function preprocessArgs!(d::Dict)
  replaceAliases!(d, _keyAliases)

  # handle axis args
  for axisletter in ("x", "y")
    asym = symbol(axisletter * "axis")
    for arg in wraptuple(get(d, asym, ()))
      processAxisArg(d, axisletter, arg)
    end
    delete!(d, asym)
  end

  # handle line args
  for arg in wraptuple(get(d, :line, ()))
    processLineArg(d, arg)
  end
  delete!(d, :line)

  # handle marker args... default to ellipse if shape not set
  anymarker = false
  for arg in wraptuple(get(d, :marker, ()))
    processMarkerArg(d, arg)
    anymarker = true
  end
  delete!(d, :marker)
  if anymarker && !haskey(d, :markershape)
    d[:markershape] = :ellipse
  end

  # handle fill
  for arg in wraptuple(get(d, :fill, ()))
    processFillArg(d, arg)
  end
  delete!(d, :fill)

  # convert into strokes and brushes

  # handle subplot links
  if haskey(d, :link)
    l = d[:link]
    if isa(l, Bool)
      d[:linkx] = l
      d[:linky] = l
    elseif isa(l, Function)
      d[:linkx] = true
      d[:linky] = true
      d[:linkfunc] = l
    else
      warn("Unhandled/invalid link $l.  Should be a Bool or a function mapping (row,column) -> (linkx, linky), where linkx/y can be Bool or Void (nothing)")
    end
    delete!(d, :link)
  end

  return
end

# -----------------------------------------------------------------------------

"A special type that will break up incoming data into groups, and allow for easier creation of grouped plots"
type GroupBy
  groupLabels::Vector{UTF8String}   # length == numGroups
  groupIds::Vector{Vector{Int}}     # list of indices for each group
end

# this is when given a vector-type of values to group by
function extractGroupArgs(v::AVec, args...)
  groupLabels = sort(collect(unique(v)))
  n = length(groupLabels)
  if n > 20
    error("Too many group labels. n=$n  Is that intended?")
  end
  groupIds = Vector{Int}[filter(i -> v[i] == glab, 1:length(v)) for glab in groupLabels]
  GroupBy(map(string, groupLabels), groupIds)
end

# expecting a mapping of "group label" to "group indices"
function extractGroupArgs{T, V<:AVec{Int}}(idxmap::Dict{T,V}, args...)
  groupLabels = sortedkeys(idxmap)
  groupIds = VecI[collect(idxmap[k]) for k in groupLabels]
  GroupBy(groupLabels, groupIds)
end

# -----------------------------------------------------------------------------

function warnOnUnsupportedArgs(pkg::PlottingPackage, d::Dict)
  for k in sortedkeys(d)
    if (!(k in supportedArgs(pkg))
        && k != :subplot
        && d[k] != default(k))
      warn("Keyword argument $k not supported with $pkg.  Choose from: $(supportedArgs(pkg))")
    end
  end
end

function warnOnUnsupported(pkg::PlottingPackage, d::Dict)
  (d[:axis] in supportedAxes(pkg)
    || warn("axis $(d[:axis]) is unsupported with $pkg.  Choose from: $(supportedAxes(pkg))"))
  (d[:linetype] == :none
    || d[:linetype] in supportedTypes(pkg)
    || warn("linetype $(d[:linetype]) is unsupported with $pkg.  Choose from: $(supportedTypes(pkg))"))
  (d[:linestyle] in supportedStyles(pkg)
    || warn("linestyle $(d[:linestyle]) is unsupported with $pkg.  Choose from: $(supportedStyles(pkg))"))
  (d[:markershape] == :none
    || d[:markershape] in supportedMarkers(pkg)
    || (Shape in supportedMarkers(pkg) && typeof(d[:markershape]) <: Shape)
    || warn("markershape $(d[:markershape]) is unsupported with $pkg.  Choose from: $(supportedMarkers(pkg))"))
end

function warnOnUnsupportedScales(pkg::PlottingPackage, d::Dict)
  for k in (:xscale, :yscale)
    if haskey(d, k)
      d[k] in supportedScales(pkg) || warn("scale $(d[k]) is unsupported with $pkg.  Choose from: $(supportedScales(pkg))")
    end
  end
end

# -----------------------------------------------------------------------------

# 1-row matrices will give an element
# multi-row matrices will give a column
# anything else is returned as-is
# getArgValue(v::Tuple, idx::Int) = v[mod1(idx, length(v))]
function getArgValue(v::AMat, idx::Int)
  c = mod1(idx, size(v,2))
  size(v,1) == 1 ? v[1,c] : v[:,c]
end
getArgValue(v, idx) = v

# given an argument key (k), we want to extract the argument value for this index.
# if nothing is set (or container is empty), return the default.
function setDictValue(d_in::Dict, d_out::Dict, k::Symbol, idx::Int, defaults::Dict)
  if haskey(d_in, k) && !(typeof(d_in[k]) <: @compat(Union{AbstractArray, Tuple}) && isempty(d_in[k]))
    d_out[k] = getArgValue(d_in[k], idx)
  else
    d_out[k] = defaults[k]
  end
end

# -----------------------------------------------------------------------------

# build the argument dictionary for the plot
function getPlotArgs(pkg::PlottingPackage, kw, idx::Int; set_defaults = true)
  kwdict = Dict(kw)
  d = Dict()

  # add defaults?
  if set_defaults
    for k in keys(_plotDefaults)
      setDictValue(kwdict, d, k, idx, _plotDefaults)
    end
  end

  for k in (:xscale, :yscale)
    if haskey(_scaleAliases, d[k])
      d[k] = _scaleAliases[d[k]]
    end
  end

  # convert color
  handlePlotColors(pkg, d)

  # no need for these
  delete!(d, :x)
  delete!(d, :y)

  d
end

# build the argument dictionary for a series
function getSeriesArgs(pkg::PlottingPackage, plotargs::Dict, kw, commandIndex::Int, plotIndex::Int, globalIndex::Int)  # TODO, pass in plotargs, not plt
  kwdict = Dict(kw)
  d = Dict()

  # add defaults?
  for k in keys(_seriesDefaults)
    setDictValue(kwdict, d, k, commandIndex, _seriesDefaults)
  end

  # groupby args?
  for k in (:idxfilter, :numUncounted, :dataframe)
    if haskey(kwdict, k)
      d[k] = kwdict[k]
    end
  end

  if haskey(_typeAliases, d[:linetype])
    d[:linetype] = _typeAliases[d[:linetype]]
  end

  aliasesAndAutopick(d, :axis, _axesAliases, supportedAxes(pkg), plotIndex)
  aliasesAndAutopick(d, :linestyle, _styleAliases, supportedStyles(pkg), plotIndex)
  aliasesAndAutopick(d, :markershape, _markerAliases, supportedMarkers(pkg), plotIndex)

  # update color
  d[:linecolor] = getSeriesRGBColor(d[:linecolor], plotargs, plotIndex)

  # update markercolor
  c = d[:markercolor]
  c = (c == :match ? d[:linecolor] : getSeriesRGBColor(c, plotargs, plotIndex))
  d[:markercolor] = c

  # update markerstrokecolor
  c = d[:markerstrokecolor]
  c = (c == :match ? plotargs[:foreground_color] : getSeriesRGBColor(c, plotargs, plotIndex))
  d[:markerstrokecolor] = c

  # update fillcolor
  c = d[:fillcolor]
  c = (c == :match ? d[:linecolor] : getSeriesRGBColor(c, plotargs, plotIndex))
  d[:fillcolor] = c

  # TODO: rebuild

  # set label
  label = d[:label]
  label = (label == "AUTO" ? "y$globalIndex" : label)
  if d[:axis] == :right && !(length(label) >= 4 && label[end-3:end] != " (R)")
    label = string(label, " (R)")
  end
  d[:label] = label

  warnOnUnsupported(pkg, d)

  d
end


text:
using Catalyst, LinearAlgebra, OrdinaryDiffEq, Test, NonlinearSolve

# Repressilator model
@parameters t α₀ α K n δ β μ
@variables m(t) P(t) R(t)
rxs = [
       Reaction(α₀, nothing, [m]),
       Reaction(α / (1 + (R/K)^n), nothing, [m]),
       Reaction(δ, [m], nothing),
       Reaction(β, [m], [m,P]),
       Reaction(μ, [P], nothing)
      ]

specs = [m,P,R]
pars  = [α₀,α,K,n,δ,β,μ]
@named rs = ReactionSystem(rxs, t, specs, pars)

# using ODESystem components
@named sys₁ = convert(ODESystem, rs; include_zero_odes=false)
@named sys₂ = convert(ODESystem, rs; include_zero_odes=false)
@named sys₃ = convert(ODESystem, rs; include_zero_odes=false)
connections = [sys₁.R ~ sys₃.P,
               sys₂.R ~ sys₁.P,
               sys₃.R ~ sys₂.P]
@named connected = ODESystem(connections, t, [], [], systems=[sys₁,sys₂,sys₃])
oderepressilator = structural_simplify(connected)

pvals = [sys₁.α₀ => 5e-4,
         sys₁.α => .5,
         sys₁.K => 40.0,
         sys₁.n => 2,
         sys₁.δ => (log(2)/120),
         sys₁.β => (20*log(2)/120),
         sys₁.μ => (log(2)/600),
         sys₂.α₀ => 5e-4,
         sys₂.α => .5,
         sys₂.K => 40.0,
         sys₂.n => 2,
         sys₂.δ => (log(2)/120),
         sys₂.β => (20*log(2)/120),
         sys₂.μ => (log(2)/600),
         sys₃.α₀ => 5e-4,
         sys₃.α => .5,
         sys₃.K => 40.0,
         sys₃.n => 2,
         sys₃.δ => (log(2)/120),
         sys₃.β => (20*log(2)/120),
         sys₃.μ => (log(2)/600)]
u₀    = [sys₁.m => 0.0, sys₁.P => 20.0, sys₁.R => 0.0, sys₂.m => 0.0, sys₂.P => 0.0, sys₂.R => 0.0, sys₃.m => 0.0, sys₃.P => 0.0, sys₃.R => 0.0]
tspan = (0.0, 100000.0)
oprob = ODEProblem(oderepressilator, u₀, tspan, pvals)
sol = solve(oprob, Tsit5())

# hardcoded network
function repress!(f, y, p, t)
    α = p.α; α₀ = p.α₀; β = p.β; δ = p.δ; μ = p.μ; K = p.K; n = p.n
    f[1] = α / (1 + (y[6] / K)^n) - δ * y[1] + α₀
    f[2] = α / (1 + (y[4] / K)^n) - δ * y[2] + α₀
    f[3] = α / (1 + (y[5] / K)^n) - δ * y[3] + α₀
    f[4] = β * y[1] - μ * y[4]
    f[5] = β * y[2] - μ * y[5]
    f[6] = β * y[3] - μ * y[6]
    nothing
end
ps = (α₀=5e-4, α=.5, K=40.0, n=2, δ=(log(2)/120), β=(20*log(2)/120), μ=(log(2)/600))
u0 = [0.0,0.0,0.0,20.0,0.0,0.0]
oprob2 = ODEProblem(repress!, u0, tspan, ps)
sol2 = solve(oprob2, Tsit5())
tvs = 0:1:tspan[end]
@test all(isapprox.(sol(tvs,idxs=sys₁.P),sol2(tvs, idxs=4),atol=1e-4))

# using ReactionSystem components
@named sys₁ = ReactionSystem(rxs, t, specs, pars)
@named sys₂ = ReactionSystem(rxs, t, specs, pars)
@named sys₃ = ReactionSystem(rxs, t, specs, pars)
connections = [ParentScope(sys₁.R) ~ ParentScope(sys₃.P),
               ParentScope(sys₂.R) ~ ParentScope(sys₁.P),
               ParentScope(sys₃.R) ~ ParentScope(sys₂.P)]
@named csys = ODESystem(connections, t, [], [])
@named repressilator = ReactionSystem(t; systems=[csys,sys₁,sys₂,sys₃])
@named oderepressilator2 = convert(ODESystem, repressilator, include_zero_odes=false)
sys2 = structural_simplify(oderepressilator2)  # FAILS currently
oprob = ODEProblem(sys2, u₀, tspan, pvals)
sol = solve(oprob, Tsit5())
@test all(isapprox.(sol(tvs,idxs=sys₁.P),sol2(tvs, idxs=4),atol=1e-4))

# test conversion to nonlinear system
@named nsys = NonlinearSystem(connections, [], [])
@named ssrepressilator = ReactionSystem(t; systems=[nsys,sys₁,sys₂,sys₃])
@named nlrepressilator = convert(NonlinearSystem, ssrepressilator, include_zero_odes=false)
sys2 = structural_simplify(nlrepressilator)
@test length(equations(sys2)) == 6
nlprob = NonlinearProblem(sys2, u₀, pvals)
sol = solve(nlprob, NewtonRaphson(), tol=1e-9)
@test sol[sys₁.P] ≈ sol[sys₂.P] ≈ sol[sys₃.P]
@test sol[sys₁.m] ≈ sol[sys₂.m] ≈ sol[sys₃.m]
@test sol[sys₁.R] ≈ sol[sys₂.R] ≈ sol[sys₃.R]

# flattening
fsys = Catalyst.flatten(ssrepressilator)
@named nlrepressilator = convert(NonlinearSystem, fsys, include_zero_odes=false)
sys2 = structural_simplify(nlrepressilator)
@test length(equations(sys2)) == 6
nlprob = NonlinearProblem(sys2, u₀, pvals)
sol = solve(nlprob, NewtonRaphson(), tol=1e-9)
@test sol[sys₁.P] ≈ sol[sys₂.P] ≈ sol[sys₃.P]
@test sol[sys₁.m] ≈ sol[sys₂.m] ≈ sol[sys₃.m]
@test sol[sys₁.R] ≈ sol[sys₂.R] ≈ sol[sys₃.R]

# test constraints
connections = [sys₁.R ~ sys₃.P,
               sys₂.R ~ sys₁.P,
               sys₃.R ~ sys₂.P]
@named csys = NonlinearSystem(connections, [sys₁.R,sys₃.P,sys₂.R,sys₁.P,sys₃.R,sys₂.P],[])
@named repressilator2 = ReactionSystem(t; constraints=csys, systems=[sys₁,sys₂,sys₃])
@named nlrepressilator = convert(NonlinearSystem, repressilator2, include_zero_odes=false)
sys2 = structural_simplify(nlrepressilator)
@test length(equations(sys2)) == 6
nlprob = NonlinearProblem(sys2, u₀, pvals)
sol = solve(nlprob, NewtonRaphson(), tol=1e-9)
@test sol[sys₁.P] ≈ sol[sys₂.P] ≈ sol[sys₃.P]
@test sol[sys₁.m] ≈ sol[sys₂.m] ≈ sol[sys₃.m]
@test sol[sys₁.R] ≈ sol[sys₂.R] ≈ sol[sys₃.R]

# test can make ODESystem
@named oderepressilator = convert(ODESystem, repressilator2, include_zero_odes=false)
sys2 = structural_simplify(oderepressilator)  # FAILS currently
oprob = ODEProblem(sys2, u₀, tspan, pvals)
sol = solve(oprob, Tsit5())
@test all(isapprox.(sol(tvs,idxs=sys₁.P),sol2(tvs, idxs=4),atol=1e-4))

# test extending with NonlinearSystem
@named repressilator2 = ReactionSystem(t; systems=[sys₁,sys₂,sys₃])
repressilator2 = Catalyst.flatten(repressilator2)
repressilator2 = extend(csys, repressilator2)
@named nlrepressilator = convert(NonlinearSystem, repressilator2, include_zero_odes=false)
sys2 = structural_simplify(nlrepressilator)
@test length(equations(sys2)) == 6
nlprob = NonlinearProblem(sys2, u₀, pvals)
sol = solve(nlprob, NewtonRaphson(), tol=1e-9)
@test sol[sys₁.P] ≈ sol[sys₂.P] ≈ sol[sys₃.P]
@test sol[sys₁.m] ≈ sol[sys₂.m] ≈ sol[sys₃.m]
@test sol[sys₁.R] ≈ sol[sys₂.R] ≈ sol[sys₃.R]

# TODO add conversion to SDE and JumpSystems once supported

# adding algebraic constraints
@parameters t, r₊, r₋, β
@variables A(t), B(t), C(t), D(t)
rxs1 = [Reaction(r₊, [A,B], [C])]
rxs2 = [Reaction(r₋, [C], [A,B])]
@named rs1 = ReactionSystem(rxs1, t, [A,B,C], [r₊])
@named rs2 = ReactionSystem(rxs2, t, [A,B,C], [r₋])
@named rs  = extend(rs1, rs2)
@test issetequal(states(rs), [A,B,C])
@test issetequal(parameters(rs), [r₊,r₋])
@test issetequal(equations(rs), union(rxs1,rxs2))
A2 = ModelingToolkit.ParentScope(A)
B2 = ModelingToolkit.ParentScope(B)
nseqs = [D ~ 2*A2 + β*B2]
@named ns = ODESystem(nseqs, t, [A2,B2,D], [β])
rs = compose(rs, [ns])
osys = convert(ODESystem, rs; include_zero_odes=false)
p  = [r₊ => 1.0, r₋ => 2.0, ns.β => 3.0]
u₀ = [A => 1.0, B => 2.0, C => 0.0]
oprob = ODEProblem(structural_simplify(osys), u₀, (0.0,10.0), p)
sol = solve(oprob, Tsit5())
@test isapprox(0, norm(sol[ns.D] .- 2*sol[A] - 3*sol[B]), atol=(100*eps()))

# test API functions for composed model
@test issetequal(species(rs), [A,B,C])
@test issetequal(states(rs), [A,B,C,ns.D])
@test issetequal(reactionparams(rs), [r₊,r₋])
@test issetequal(parameters(rs), [r₊,r₋,ns.β])
@test issetequal(reactions(rs), union(rxs1,rxs2))
@test issetequal(filter(eq -> eq isa Reaction, equations(rs)), union(rxs1,rxs2))
@test issetequal(filter(eq -> eq isa Equation, equations(rs)), [ModelingToolkit.namespace_equation(nseqs[1],ns)])

# check several levels of nesting namespace and filter ok for the API functions
@parameters t, p1, p2a, p2b, p3a, p3b
@variables A1(t), A2a(t), A2b(t), A3a(t), A3b(t)
rxs1 = [Reaction(p1, [A1], nothing)]
rxs2 = [Reaction(p2a, [A2a], nothing), Reaction(p2b, [ParentScope(A1)], nothing)]
eqs2 = [ParentScope(A1) ~ ParentScope(p1)*A2b]
rxs3 = [Reaction(p3a, [A3a], nothing), Reaction(ParentScope(p2a), nothing, [ParentScope(A2a)])]
eqs3 = [ParentScope(A2a) ~ p3b*A3b]

@named rs3 = ReactionSystem(rxs3,t,[A3a,ParentScope(A2a)],[p3a,ParentScope(p2a)])
@named ns3 = NonlinearSystem(eqs3,[ParentScope(A2a),A3b],[p3b])
@named rs2 = ReactionSystem(rxs2,t,[A2a,ParentScope(A1)],[p2a,p2b], systems=[rs3,ns3])
@named ns2 = NonlinearSystem(eqs2,[ParentScope(A1),A2b],[ParentScope(p1)])
@named rs1 = ReactionSystem(rxs1,t,[A1],[p1],systems=[rs2,ns2])

# namespaced reactions
nrxs1 = [Reaction(p1, [A1], nothing)]
nrxs2 = [Reaction(rs2.p2a, [rs2.A2a], nothing), Reaction(rs2.p2b, [A1], nothing)]
neqs2 = [0 ~ p1*ns2.A2b - A1]
nrxs3 = [Reaction(rs2.rs3.p3a, [rs2.rs3.A3a], nothing), Reaction(rs2.p2a, nothing, [rs2.A2a])]
neqs3 = [0 ~ rs2.ns3.p3b*rs2.ns3.A3b - rs2.A2a]
rxs = vcat(nrxs1, nrxs2, nrxs3)
eqs = vcat(nrxs1, nrxs2, neqs2, nrxs3, neqs3)

@test issetequal(states(rs1), [A1,rs2.A2a,ns2.A2b,rs2.rs3.A3a,rs2.ns3.A3b])
@test issetequal(species(rs1), [A1,rs2.A2a,rs2.rs3.A3a])
@test issetequal(parameters(rs1), [p1,rs2.p2a,rs2.p2b,rs2.rs3.p3a,rs2.ns3.p3b])
@test issetequal(reactionparams(rs1), [p1, rs2.p2a, rs2.p2b, rs2.rs3.p3a])
@test issetequal(rxs, reactions(rs1))
@test issetequal(eqs, equations(rs1))

# test throw error if there are ODE constraints and convert to NonlinearSystem
rn = @reaction_network rn begin
    (k1,k2), A <--> B
    end k1 k2
@parameters a,b
@variables t,A(t),C(t)
D = Differential(t)
eqs = [D(A) ~ -a*A + C, D(C) ~ -b*C + a*A]
@named osys = ODESystem(eqs,t,[A,C],[a,b])
rn2 = extend(osys,rn)
rnodes = convert(ODESystem,rn2)
@test_throws ErrorException convert(NonlinearSystem,rn2)

# ensure right number of equations are generated
@variables G(t)
eqs = [D(G) ~ -G]
@named osys2 = ODESystem(eqs,t)
rn3 = compose(rn2, osys2)
@test length(equations(rn3)) == 5

# check conversions work with algebraic constraints
eqs = [0 ~ -a*A + C, 0 ~ -b*C + a*A]
@named nlsys = NonlinearSystem(eqs,[A,C],[a,b])
rn2 = extend(nlsys,rn)
rnodes = convert(ODESystem,rn2)
rnnlsys = convert(NonlinearSystem,rn2)
@named nlsys = ODESystem(eqs,t,[A,C],[a,b])
rn2 = extend(nlsys,rn)
rnodes = convert(ODESystem,rn2)
rnnlsys = convert(NonlinearSystem,rn2)

# https://github.com/SciML/ModelingToolkit.jl/issues/1274
@parameters p1 p2
@variables t A(t)
rxs1 = [Reaction(p1, [A], nothing)]
rxs2 = [Reaction(p2, [ParentScope(A)], nothing)]
@named rs1 = ReactionSystem(rxs1, t)
@named rs2 = ReactionSystem(rxs2, t)
rsc = compose(rs1,[rs2])
orsc = convert(ODESystem, rsc)
@test length(equations(orsc)) == 1

ext: jl


text:
using CutPruners
using Compat, Compat.Test

include("redund.jl")
include("start.jl")
include("oldnew.jl")
include("avg.jl")
include("decay.jl")
include("dematos.jl")
include("exact.jl")


text:
struct FullTimeGF{T, scalar} <: AbstractTimeGF{T, scalar}
  grid::FullTimeGrid
  gtr::AntiHermitianStorage{T,scalar}
  les::AntiHermitianStorage{T,scalar}
  rm::GenericStorage{T,scalar}
  mat::PeriodicStorage{T,scalar}
  ξ::GFSignEnum
end

norbitals(G::FullTimeGF) = G.gtr.norb

function FullTimeGF(::Type{T}, grid::FullTimeGrid, norb=1, ξ::GFSignEnum=fermionic, scalar=false) where T <: Number
  nt = grid.nt
  ntau = grid.ntau

  gtr = AntiHermitianStorage(T, nt, norb, scalar)
  les = AntiHermitianStorage(T, nt, norb, scalar)
  rm = GenericStorage(T, ntau, nt, norb, scalar)
  mat = PeriodicStorage(T, ntau, norb, scalar)

  FullTimeGF(grid, gtr, les, rm, mat, ξ)
end
FullTimeGF(grid::FullTimeGrid, norb=1, ξ::GFSignEnum=fermionic, scalar=false) = FullTimeGF(ComplexF64, grid, norb, ξ, scalar)

@inline function Base.getindex(G::FullTimeGF, k, l, t1::TimeGridPoint, t2::TimeGridPoint, greater=true)
  greater = t1 == t2 ? greater : heaviside(t1.bpoint, t2.bpoint)

  i = t1.ridx
  j = t2.ridx
  ξ = Int(G.ξ)

  real_branches = (forward_branch, backward_branch)
  if (t1.bpoint.domain ∈ real_branches && t2.bpoint.domain ∈ real_branches)
    return greater ? G.gtr[k,l,i,j] : G.les[k,l,i,j]
  elseif (t1.bpoint.domain == imaginary_branch && (t2.bpoint.domain == forward_branch || t2.bpoint.domain == backward_branch))
    return G.rm[k,l,i,j]
  elseif ((t1.bpoint.domain == forward_branch || t1.bpoint.domain == backward_branch) && t2.bpoint.domain == imaginary_branch)
    ntau = G.grid.ntau
    return -ξ * adjoint(G.rm[l,k,ntau+1-j,i]) # Aoki 19c
  else
    return greater ? G.mat[k,l,i,j,greater] : ξ * G.mat[k,l,i,j,greater]
  end
end

function Base.setindex!(G::FullTimeGF, v, k, l, t1::TimeGridPoint, t2::TimeGridPoint)
  greater = heaviside(t1.bpoint, t2.bpoint)

  i = t1.ridx
  j = t2.ridx
  ξ = Int(G.ξ)

  real_branches = (forward_branch, backward_branch)
  if (t1.bpoint.domain ∈ real_branches && t2.bpoint.domain ∈ real_branches)
    return greater ? G.gtr[k,l,i,j] = v : G.les[k,l,i,j] = v
  elseif (t1.bpoint.domain == imaginary_branch && (t2.bpoint.domain == forward_branch || t2.bpoint.domain == backward_branch))
    return G.rm[k,l,i,j] = v
  elseif ((t1.bpoint.domain == forward_branch || t1.bpoint.domain == backward_branch) && t2.bpoint.domain == imaginary_branch)
    ntau = G.grid.ntau
    return G.rm[l,k,ntau+1-j,i] = -ξ * adjoint(v) # Aoki 19c
  else
    if greater
      return G.mat[k,l,i,j] = v
    else
      return G.mat[k,l,i,j] = ξ * v
    end
  end
end

function Base.similar(G::T, ξ::GFSignEnum) where T <: FullTimeGF
  T(G.grid, similar(G.gtr), similar(G.les), similar(G.rm), similar(G.mat), ξ)
end

function Base.zero(G::T, ξ::GFSignEnum) where T <: FullTimeGF
  T(G.grid, zero(G.gtr), zero(G.les), zero(G.rm), zero(G.mat), ξ)
end

function Base.similar(G::T) where T <: FullTimeGF
  similar(G, G.ξ)
end

function Base.zero(G::T) where T <: FullTimeGF
  zero(G, G.ξ)
end

function TimeDomain(G::FullTimeGF)
  grid = G.grid
  tplus = grid[forward_branch]
  tminus = reverse(grid[backward_branch])
  tau = grid[imaginary_branch]

  nt = G.grid.nt
  ntau = G.grid.ntau

  points = Tuple{TimeGridPoint, TimeGridPoint}[]

  for i in 1:nt
    for j in 1:i
      push!(points, (tminus[j], tminus[i])) # greater
      push!(points, (tplus[j], tminus[i])) # lesser
    end
  end

  for i in 1:nt
    for j in 1:ntau
      push!(points, (tau[j], tminus[i])) # right mixing
    end
  end

  for i in 1:ntau
    push!(points, (tau[i], tau[1])) # matsubara
  end

  return TimeDomain(points)
end

function FullTimeGF(f::Function, ::Type{T}, grid::FullTimeGrid, norb=1, ξ::GFSignEnum=fermionic, scalar=false) where T <: Number
  G = FullTimeGF(T, grid, norb, ξ, scalar)
  D = TimeDomain(G)

  for (t1, t2) in D.points
    G[t1, t2] = f(t1, t2)
  end

  return G
end

FullTimeGF(f::Function, grid::FullTimeGrid, norb=1, ξ::GFSignEnum=fermionic, scalar=false) = FullTimeGF(f, ComplexF64, grid, norb, ξ, scalar)

function FullTimeGF(dos::AbstractDOS, grid::FullTimeGrid)
  # result is time invariant so do calculation for time invariant case and then copy over
  G = TimeInvariantFullTimeGF(dos, grid)
  FullTimeGF(grid,1,fermionic,true) do t1, t2
    G[t1,t2]
  end
end


text:
@testset "Regular Expression" begin
    @test isa(biore"A+"d, BioSequences.RE.Regex{DNA})
    @test isa(biore"A+"r, BioSequences.RE.Regex{RNA})
    @test isa(biore"A+"a, BioSequences.RE.Regex{AminoAcid})
    @test isa(biore"A+"dna, BioSequences.RE.Regex{DNA})
    @test isa(biore"A+"rna, BioSequences.RE.Regex{RNA})
    @test isa(biore"A+"aa, BioSequences.RE.Regex{AminoAcid})
    @test_throws Exception eval(:(biore"A+"))
    @test_throws Exception eval(:(biore"A+"foo))
    @test string(biore"A+"dna) == "biore\"A+\"dna"
    @test string(biore"A+"rna) == "biore\"A+\"rna"
    @test string(biore"A+"aa) == "biore\"A+\"aa"

    @test  occursin(biore"A"d, dna"A")
    @test !occursin(biore"A"d, dna"C")
    @test !occursin(biore"A"d, dna"G")
    @test !occursin(biore"A"d, dna"T")
    @test  occursin(biore"N"d, dna"A")
    @test  occursin(biore"N"d, dna"C")
    @test  occursin(biore"N"d, dna"G")
    @test  occursin(biore"N"d, dna"T")
    @test  occursin(biore"[AT]"d, dna"A")
    @test !occursin(biore"[AT]"d, dna"C")
    @test !occursin(biore"[AT]"d, dna"G")
    @test  occursin(biore"[AT]"d, dna"T")
    @test !occursin(biore"[^AT]"d, dna"A")
    @test  occursin(biore"[^AT]"d, dna"C")
    @test  occursin(biore"[^AT]"d, dna"G")
    @test !occursin(biore"[^AT]"d, dna"T")

    re = biore"^A(C+G*)(T{2,})N$"d
    @test !occursin(re, dna"AC")
    @test !occursin(re, dna"AGTT")
    @test !occursin(re, dna"CCGTT")
    @test !occursin(re, dna"ACTT")
    @test !occursin(re, dna"ACTTGT")
    @test  occursin(re, dna"ACGTTA")
    @test  occursin(re, dna"ACGTTT")
    @test  occursin(re, dna"ACCGGTTT")
    @test  occursin(re, dna"ACCGGTTT")
    @test  occursin(re, dna"ACCGGTTTA")
    @test  occursin(re, dna"ACCGGTTTG")

    @test matched(match(re, dna"ACCGTTTTA")) == dna"ACCGTTTTA"
    @test captured(match(re, dna"ACCGTTTTA"))[1] == dna"CCG"
    @test captured(match(re, dna"ACCGTTTTA"))[2] == dna"TTTT"

    # greedy
    @test matched(match(biore"A*"d, dna"AAA")) == dna"AAA"
    @test matched(match(biore"A+"d, dna"AAA")) == dna"AAA"
    @test matched(match(biore"A?"d, dna"AAA")) == dna"A"
    @test matched(match(biore"A{2}"d, dna"AAA")) == dna"AA"
    @test matched(match(biore"A{2,}"d, dna"AAA")) == dna"AAA"
    @test matched(match(biore"A{2,4}"d, dna"AAA")) == dna"AAA"

    # lazy
    @test matched(match(biore"A*?"d, dna"AAA")) == dna""
    @test matched(match(biore"A+?"d, dna"AAA")) == dna"A"
    @test matched(match(biore"A??"d, dna"AAA")) == dna""
    @test matched(match(biore"A{2}?"d, dna"AAA")) == dna"AA"
    @test matched(match(biore"A{2,}?"d, dna"AAA")) == dna"AA"
    @test matched(match(biore"A{2,4}?"d, dna"AAA")) == dna"AA"

    # search
    @test findfirst(biore"A+"d, dna"ACGTAAT") == 1:1
    @test findfirst(biore"A+"d, dna"ACGTAAT", 1) == 1:1
    @test findfirst(biore"A+"d, dna"ACGTAAT", 2) == 5:6
    @test findfirst(biore"A+"d, dna"ACGTAAT", 7) == nothing

    # eachmatch
    matches = [dna"CG", dna"GC", dna"GC", dna"CG"]
    for (i, m) in enumerate(eachmatch(biore"GC|CG"d, dna"ACGTTATGCATGGCG"))
        @test matched(m) == matches[i]
    end
    matches = [dna"CG", dna"GC", dna"GC"]
    for (i, m) in enumerate(collect(eachmatch(biore"GC|CG"d, dna"ACGTTATGCATGGCG", false)))
        @test matched(m) == matches[i]
    end

    # matchall
    matchall(pat, seq, overlap=true) =
        collect(map(matched, eachmatch(pat, seq, overlap)))

    @test matchall(biore"A*"d, dna"") == [dna""]
    @test matchall(biore"A*"d, dna"AAA") == [
        dna"AAA", dna"AA", dna"A", dna"",
        dna"AA",  dna"A",  dna"",
        dna"A",   dna""]
    @test matchall(biore"AC*G*T"d, dna"ACCGGGT") == [dna"ACCGGGT"]

    @test matchall(biore"A*"d, dna"", false) == [dna""]
    @test matchall(biore"A*"d, dna"AAA", false) == [dna"AAA"]
    @test matchall(biore"AC*G*T"d, dna"ACCGGGT", false) == [dna"ACCGGGT"]

    # RNA and Amino acid
    @test  occursin(biore"U(A[AG]|GA)$"r, rna"AUUGUAUGA")
    @test !occursin(biore"U(A[AG]|GA)$"r, rna"AUUGUAUGG")
    @test  occursin(biore"T+[NQ]A?P"a, aa"MTTQAPMFTQPL")
    @test  occursin(biore"T+[NQ]A?P"a, aa"MTTAAPMFTQPL")
    @test !occursin(biore"T+[NQ]A?P"a, aa"MTTAAPMFSQPL")

    # PROSITE
    @test  occursin(prosite"[AC]-x-V-x(4)-{ED}", aa"ADVAARRK")
    @test  occursin(prosite"[AC]-x-V-x(4)-{ED}", aa"CPVAARRK")
    @test !occursin(prosite"[AC]-x-V-x(4)-{ED}", aa"ADVAARRE")
    @test !occursin(prosite"[AC]-x-V-x(4)-{ED}", aa"CPVAARK")
    @test  occursin(prosite"<[AC]-x-V-x(4)-{ED}>", aa"ADVAARRK")
    @test !occursin(prosite"<[AC]-x-V-x(4)-{ED}>", aa"AADVAARRK")
    @test !occursin(prosite"<[AC]-x-V-x(4)-{ED}>", aa"ADVAARRKA")
end


text:
################################################################################
#
#	ABSTRACT TYPE
#
#   ReciprocalUnitcell{P,B} <: Unitcell{P,B}
#   --> P is the reciprocal site type (<: AbstractReciprocalPoint{D})
#       --> D is the dimension of embedding space
#   --> B is the bond type (<: AbstractBond{L,N})
#       --> N is the dimension of the Bravais lattice that the bond is located in
#       --> L is the label type of bonds
#
#   FILE CONTAINS
#       - abstract type definition
#       - interface definition
#       - TODO interface testing
#
################################################################################

################################################################################
#
#   ABSTRACT TYPE DEFINITION
#
################################################################################
abstract type AbstractReciprocalUnitcell{
        P <: AbstractReciprocalPoint{D} where {D},
        B <: AbstractBond{L,N} where {L,N}
    } <: AbstractUnitcell{P,B} end

# export the type
export AbstractReciprocalUnitcell

################################################################################
#
#	INTERFACING / ACCESSING UNITCELLS
#	(functions that are overwritten by abstract reciprocal type, NOT concrete)
#
################################################################################

# default constructor interface
# used for creation of new unitcells
function newUnitcell(
            ::Type{U},
            lattice_vectors :: Vector{<:Vector{<:Real}},
            sites           :: Vector{S},
            bonds           :: Vector{B}
        ) :: U where {D,N,L,S<:AbstractReciprocalPoint{D},B<:AbstractBond{L,N},U<:AbstractReciprocalUnitcell{S,B}}

    # check if only a single site is given
    @assert length(sites) == 1
    # check if the site is located at 0
    @assert sum(abs.(point(sites[1]))) < 1e-8

    # pass to the reciprocal unitcell constructor interface
    return newReciprocalUnitcell(U, lattice_vectors, bonds)
end

# export the newUnitcell interface
export newUnitcell

# accessing a list of sites
function sites(
            unitcell :: U
        ) :: Vector{S} where {D,N,L,S<:AbstractReciprocalPoint{D},B<:AbstractBond{L,N},U<:AbstractReciprocalUnitcell{S,B}}

    # return the Gamma point inside an array
    return S[gammaPoint(unitcell),]
end
# setting a list of sites
function sites!(
            unitcell :: U,
            sites    :: Vector{S}
        ) where {D,N,L,S<:AbstractReciprocalPoint{D},B<:AbstractBond{L,N},U<:AbstractReciprocalUnitcell{S,B}}

    # check if only a single site is given
    @assert length(sites) == 1
    # check if the site is located at 0
    @assert sum(abs.(point(sites[1]))) < 1e-8

    # set the Gamma point
    gammaPoint!(unitcell, sites[1])
end

# export the sites interface
export sites, sites!

################################################################################
#
#	INTERFACING / ACCESSING RECIPROCAL UNITCELLS
#	(functions that have to be overwritten by concrete types)
#
################################################################################

# constructor
function newReciprocalUnitcell(
            ::Type{R},
            lattice_vectors :: Vector{<:Vector{<:Real}},
            bonds           :: Vector{B}
        ) :: R where {D,N,L,P<:AbstractReciprocalPoint{D}, B<:AbstractBond{N,L}, R<:AbstractReciprocalUnitcell{P,B}}

    # throw an error because implementation for concrete type is missing
    error(  "not implemented function 'newReciprocalUnitcell' for concrete reciprocal unitcell type " *
            string(R) * " with bond type " * string(B) )
end

# export the constructor
export newReciprocalUnitcell

# accessing a list of lattice vectors
function latticeVectors(
            unitcell :: U
        ) :: Vector{Vector{Float64}} where {P,B,U<:AbstractReciprocalUnitcell{P,B}}

    # print an error because implementation for concrete type is missing
    error(  "not implemented function 'latticeVectors' for concrete reciprocal unitcell type " * string(U) )
end
# setting a list of lattice vectors
function latticeVectors!(
            unitcell        :: U,
            lattice_vectors :: Vector{<:Vector{<:Real}}
        ) where {P,B,U<:AbstractReciprocalUnitcell{P,B}}

    # print an error because implementation for concrete type is missing
    error(  "not implemented function 'latticeVectors!' for concrete reciprocal unitcell type " * string(U) )
end

# export the latticeVectors interface
export latticeVectors, latticeVectors!

# accessing the gamma point
function gammaPoint(
            unitcell :: U
        ) :: P where {D,N,L,P<:AbstractReciprocalPoint{D},B<:AbstractBond{L,N},U<:AbstractReciprocalUnitcell{P,B}}

    # print an error because implementation for concrete type is missing
    error(  "not implemented function 'gammaPoint' for concrete reciprocal unitcell type " * string(U) )
end

# setting the gamma point
function gammaPoint!(
            unitcell :: U,
            site     :: P
        ) :: P where {D,N,L,P<:AbstractReciprocalPoint{D},B<:AbstractBond{L,N},U<:AbstractReciprocalUnitcell{P,B}}

    # print an error because implementation for concrete type is missing
    error(  "not implemented function 'gammaPoint!' for concrete reciprocal unitcell type " * string(U) )
end

# export gamma point interface
export gammaPoint, gammaPoint!

# accessing a list of bonds
function bonds(
            unitcell :: U
        ) :: Vector{B} where {N,L,P,B<:AbstractBond{L,N},U<:AbstractReciprocalUnitcell{P,B}}

    # print an error because implementation for concrete type is missing
    error(  "not implemented function 'bonds' for concrete reciprocal unitcell type " *
            string(U) * " with bond type " * string(B)   )
end
# setting a list of bonds
function bonds!(
            unitcell :: U,
            bonds    :: Vector{B}
        ) where {N,L,P,B<:AbstractBond{L,N},U<:AbstractReciprocalUnitcell{P,B}}

    # print an error because implementation for concrete type is missing
    error(  "not implemented function 'bonds!' for concrete reciprocal unitcell type " *
            string(U) * " with bond type " * string(B)   )
end

# export the bonds interface
export bonds, bonds!

################################################################################
#
#   ADDITIONAL SYNTAX FOR RECIPROCAL UNITCELLS
#
################################################################################

# specific to reciprocal Bravais lattices (just rename a1 etc to b1 etc)

# b1
function b1(
            unitcell :: U
        ) :: Vector{Float64} where {D,N,L,P<:AbstractReciprocalPoint{D},B<:AbstractBond{L,N},U<:AbstractReciprocalUnitcell{P,B}}

    # return the first entry in the lattice vectors list
    return latticeVectors(unitcell)[1]
end

# b2
function b2(
            unitcell :: U
        ) :: Vector{Float64} where {D,N,L,P<:AbstractReciprocalPoint{D},B<:AbstractBond{L,N},U<:AbstractReciprocalUnitcell{P,B}}

    # return the second entry in the lattice vectors list
    return latticeVectors(unitcell)[2]
end

# b3
function b3(
            unitcell :: U
        ) :: Vector{Float64} where {D,N,L,P<:AbstractReciprocalPoint{D},B<:AbstractBond{L,N},U<:AbstractReciprocalUnitcell{P,B}}

    # return the third entry in the lattice vectors list
    return latticeVectors(unitcell)[3]
end

# export the Bravais lattice interfaces
export b1, b2, b3


text:
const grh = Genie.Renderer.Html

include("table/table_tables.jl")
include("table/table_csv.jl")
include("table/table_dataframes.jl")


text:
using Test

include("run-length-encoding.jl")

# Tests adapted from `problem-specifications//canonical-data.json` @ v1.0.0
# Encode and decode the strings under the given specifications.

@testset "encode strings" begin
    @test encode("") == ""
    @test encode("XYZ") == "XYZ"
    @test encode("AABBBCCCC") == "2A3B4C"
    @test encode("WWWWWWWWWWWWBWWWWWWWWWWWWBBBWWWWWWWWWWWWWWWWWWWWWWWWB") == "12WB12W3B24WB"
    @test encode("aabbbcccc") == "2a3b4c"
    @test encode("  hsqq qww  ") == "2 hs2q q2w2 "
end

@testset "decode strings" begin
    @test decode("") == ""
    @test decode("XYZ") == "XYZ"
    @test decode("2A3B4C") == "AABBBCCCC"
    @test decode("12WB12W3B24WB") == "WWWWWWWWWWWWBWWWWWWWWWWWWBBBWWWWWWWWWWWWWWWWWWWWWWWWB"
    @test decode("2a3b4c") == "aabbbcccc"
    @test decode("2 hs2q q2w2 ") == "  hsqq qww  "
end

@testset "encode then decode" begin
    @test decode(encode("zzz ZZ  zZ")) == "zzz ZZ  zZ"
end


text:
doc"""
    Gumbel(μ, θ)

The *Gumbel distribution*  with location `μ` and scale `θ` has probability density function

$f(x; \mu, \theta) = \frac{1}{\theta} e^{-(z + e^z)},
\quad \text{ with } z = \frac{x - \mu}{\theta}$

```julia
Gumbel()            # Gumbel distribution with zero location and unit scale, i.e. Gumbel(0, 1)
Gumbel(u)           # Gumbel distribution with location u and unit scale, i.e. Gumbel(u, 1)
Gumbel(u, b)        # Gumbel distribution with location u and scale b

params(d)        # Get the parameters, i.e. (u, b)
location(d)      # Get the location parameter, i.e. u
scale(d)         # Get the scale parameter, i.e. b
```

External links

* [Gumbel distribution on Wikipedia](http://en.wikipedia.org/wiki/Gumbel_distribution)
"""
immutable Gumbel{T<:Real} <: ContinuousUnivariateDistribution
    μ::T  # location
    θ::T  # scale

    Gumbel(μ::T, θ::T) = (@check_args(Gumbel, θ > zero(θ)); new(μ, θ))
end

Gumbel{T<:Real}(μ::T, θ::T) = Gumbel{T}(μ, θ)
Gumbel(μ::Real, θ::Real) = Gumbel(promote(μ, θ)...)
Gumbel(μ::Integer, θ::Integer) = Gumbel(Float64(μ), Float64(θ))
Gumbel(μ::Real) = Gumbel(μ, 1.0)
Gumbel() = Gumbel(0.0, 1.0)

@distr_support Gumbel -Inf Inf

const DoubleExponential = Gumbel

#### Conversions

convert{T <: Real, S <: Real}(::Type{Gumbel{T}}, μ::S, θ::S) = Gumbel(T(μ), T(θ))
convert{T <: Real, S <: Real}(::Type{Gumbel{T}}, d::Gumbel{S}) = Gumbel(T(d.μ), T(d.θ))

#### Parameters

location(d::Gumbel) = d.μ
scale(d::Gumbel) = d.θ
params(d::Gumbel) = (d.μ, d.θ)
@inline partype{T<:Real}(d::Gumbel{T}) = T

#### Statistics

mean(d::Gumbel) = d.μ + d.θ * γ

median{T<:Real}(d::Gumbel{T}) = d.μ - d.θ * log(T(logtwo))

mode(d::Gumbel) = d.μ

var{T<:Real}(d::Gumbel{T}) = T(π)^2/6 * d.θ^2

skewness{T<:Real}(d::Gumbel{T}) = 112*sqrt(T(6))*zeta(T(3))/π/π/π

kurtosis{T<:Real}(d::Gumbel{T}) = T(12)/5

entropy(d::Gumbel) = log(d.θ) + 1 + γ

#### Evaluation

zval(d::Gumbel, x::Real) = (x - d.μ) / d.θ
xval(d::Gumbel, z::Real) = x * d.θ + d.μ

function pdf(d::Gumbel, x::Real)
    z = zval(d, x)
    exp(-z - exp(-z)) / d.θ
end

function logpdf(d::Gumbel, x::Real)
    z = zval(d, x)
    - (z + exp(-z) + log(d.θ))
end

cdf(d::Gumbel, x::Real) = exp(-exp(-zval(d, x)))
logcdf(d::Gumbel, x::Real) = -exp(-zval(d, x))

quantile(d::Gumbel, p::Real) = d.μ - d.θ * log(-log(p))

gradlogpdf(d::Gumbel, x::Real) = - (1 + exp((d.μ - x) / d.θ)) / d.θ

#### Sampling

rand(d::Gumbel) = quantile(d, rand())


text:
function ensure_package_compiler()
    Pkg.activate()  # root environment
    Pkg.add("PackageCompiler")
    return nothing
end

function output_dir()
    dir = joinpath(DEPOT_PATH[1], "images", string(VERSION))
    mkpath(dir)  # ensure path exists
    return dir
end

function activate_build_directory(c::SysImageConfig)
    build_dir = mktempdir()
    packages = string.(getproperty.(c.dependencies, :package))
    Pkg.activate(build_dir)
    Pkg.add.(packages)
    return build_dir
end

function Base.pathof(dep::Dependency)
    return Base.locate_package(Base.PkgId(String(dep.package)))
end

function deactivate_build_directory()
    Pkg.activate()
end

function find_precompile_files(c::SysImageConfig)
    files = String[]

    # get pacakge-provided precompile files
    for pkg in c.dependencies
        if pkg.precompile_file !== nothing
            push!(files, pathof(pkg) * "../../../" * pkg.precompile_file)
        end
    end

    # add custom precompile files
    custom_dir = joinpath(@__DIR__, "../precompile/$(c.name)")
    if isdir(custom_dir)
        append!(files, readdir(custom_dir; join = true))
    end
    return files
end

function make_sysimage(
    c::SysImageConfig,
    output_file::String,
    precompile_files::Vector{<:AbstractString},
)
    modules = [dep.package for dep in c.dependencies]
    return create_sysimage(
        modules,
        sysimage_path = output_file,
        precompile_execution_file = precompile_files,
    )
end

function build(c::SysImageConfig)
    output_file = joinpath(output_dir(), "$(c.name).so")
    try
        activate_build_directory(c)
        precompile_files = find_precompile_files(c)
        @info "files" precompile_files
        return make_sysimage(c, output_file, precompile_files)
    finally
        deactivate_build_directory()
    end
end

function build(image_name::String)
    idx = findfirst(c -> c.name == image_name, CONFIGS)
    if idx !== nothing
        build(CONFIGS[idx])
    else
        println("ERROR: image config not found for `$image_name`")
        println("Please use one of the following supported images:")
        foreach(c -> println("  ", c.name), CONFIGS)
    end
end


text:
# # Linear rotating shallow water dynamics
#
#md # This example can be run online via [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/literated/OneDShallowWaterGeostrophicAdjustment.ipynb).
#md # Also, it can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/literated/OneDShallowWaterGeostrophicAdjustment.ipynb).
#
#
# This example solves the linear 1D rotating shallow water equations
# for the ``u(x, t)``, ``v(x, t)`` and the surface surface elevation ``\eta(x, t)``,
# for a fluid with constant rest-depth ``H``. That is, the total fluid's depth
# is ``H + \eta(x, t)`` with ``|\eta| \ll H``.
#
# The linearized equations for the evolution of ``u``, ``v``, ``\eta`` are:
#
# ```math
# \begin{aligned}
# \partial_t u - f v & = - g \partial_x \eta - \mathrm{D} u, \\
# \partial_t v + f u & = - \mathrm{D} v, \\
# \partial_t \eta + H \partial_x u & = - \mathrm{D} \eta.
# \end{aligned}
# ```
#
# Above, ``g`` is the gravitational acceleration, ``f`` is the  Coriolis parameter, and
# ``\mathrm{D}`` indicates a hyperviscous linear operator of the form ``(-1)^{n_ν} ν \nabla^{2 n_ν}``,
# with ``ν`` the viscosity coefficient and ``n_ν`` the order of the operator.
#
# Rotation introduces the deformation length scale, ``L_d = \sqrt{g H} / f``. Disturbances with
# length scales much smaller than ``L_d`` don't "feel" the rotation and propagate as inertia-gravity
# waves. Disturbances with length scales comparable or larger than ``L_d`` should be approximately
# in geostrophic balance, i.e., the Coriolis acceleration ``f \widehat{\bm{z}} \times \bm{u}``
# should be in approximate balance with the pressure gradient ``-g \bm{\nabla} \eta``.

using FourierFlows, Plots, Printf, Random
using LinearAlgebra: mul!, ldiv!

# ## Coding up the equations
# ### A demonstration of FourierFlows.jl framework
#
# What follows is a step-by-step tutorial demonstrating how you can create your own solver
# for an equation of your liking.

# The basic building blocks for a `FourierFlows.Problem` are:
# - `Grid` struct containining the physical and wavenumber grid for the problem,
# - `Params` struct containining all the parameters of the problem,
# - `Vars` struct containining arrays with the variables used in the problem,
# - `Equation` struct containining the coefficients of the linear operator ``L`` and the function that computes the nonlinear terms, usually named `calcN!()`.
#
# The `Grid` structure is provided by FourierFlows.jl. We simply have to call one of either
# `OneDGrid()`, `TwoDGrid()`, or `ThreeDGrid()` constructors, depending on the dimensionality
# of the problem. All other structs mentioned above are problem-specific and need to be constructed
# for every set of equations we want to solve.

# First let's construct the `Params` struct that contains all parameters of the problem.

struct Params{T} <: AbstractParams
   ν :: T         # Hyperviscosity coefficient
  nν :: Int       # Order of the hyperviscous operator
   g :: T         # Gravitational acceleration
   H :: T         # Fluid depth
   f :: T         # Coriolis parameter
end
nothing #hide

# Now the `Vars` struct that contains all variables used in this problem. For this
# problem `Vars` includes the representations of the flow fields in physical space
# `u`, `v` and `η` and their Fourier transforms `uh`, `vh`, and `ηh`.

struct Vars{Aphys, Atrans} <: AbstractVars
   u :: Aphys
   v :: Aphys
   η :: Aphys
  uh :: Atrans
  vh :: Atrans
  ηh :: Atrans
end
nothing #hide

# A constructor populates empty arrays based on the dimension of the `grid`
# and then creates `Vars` struct.
"""
    Vars(dev, grid)
Constructs Vars for 1D shallow water based on the dimensions of arrays of the `grid`.
"""
function Vars(::Dev, grid) where Dev
  T = eltype(grid)
  @devzeros Dev T grid.nx u v η
  @devzeros Dev Complex{T} grid.nkr uh vh ηh

  return Vars(u, v, η, uh, vh, ηh)
end
nothing #hide

# In Fourier space, the 1D linear shallow water dynamics read:
#
# ```math
# \begin{aligned}
# \frac{\partial \hat{u}}{\partial t} & = \underbrace{ f \hat{v} - i k g \hat{\eta} }_{N_u} \; \underbrace{- \nu k^2 }_{L_u} \hat{u} , \\
# \frac{\partial \hat{v}}{\partial t} & = \underbrace{ - f \hat{u} }_{N_v} \; \underbrace{- \nu k^2 }_{L_v} \hat{v} , \\
# \frac{\partial \hat{\eta}}{\partial t} & = \underbrace{ - i k H \hat{u} }_{N_{\eta}} \; \underbrace{- \nu k^2 }_{L_{\eta}} \hat{\eta} .
# \end{aligned}
# ```
# Although, e.g., terms involving the Coriolis accelaration are, in principle, linear we include
# them in the nonlinear term ``N`` because they render the linear operator ``L`` non-diagonal.
#
# With these in mind, we construct function `calcN!` that computes the nonlinear terms.
#
"""
    calcN!(N, sol, t, clock, vars, params, grid)
Compute the nonlinear terms for 1D linear shallow water dynamics.
"""
function calcN!(N, sol, t, clock, vars, params, grid)
  @. vars.uh = sol[:, 1]
  @. vars.vh = sol[:, 2]
  @. vars.ηh = sol[:, 3]

  @. N[:, 1] =   params.f * vars.vh - im * grid.kr * params.g * vars.ηh    #  + f v - g ∂η/∂x
  @. N[:, 2] = - params.f * vars.uh                                        #  - f u
  @. N[:, 3] = - im * grid.kr * params.H * vars.uh                         #  - H ∂u/∂x

  dealias!(N, grid)

  return nothing
end
nothing #hide

# Next we construct the `Equation` struct:

"""
    Equation(dev, params, grid)
Construct the equation: the linear part, in this case the hyperviscous dissipation,
and the nonlinear part, which is computed by `calcN!` function.
"""
function Equation(dev, params, grid)
  T = eltype(grid)
  L = zeros(dev, T, (grid.nkr, 3))
  D = @. - params.ν * grid.kr^(2*params.nν)

  L[:, 1] .= D # for u equation
  L[:, 2] .= D # for v equation
  L[:, 3] .= D # for η equation

  return FourierFlows.Equation(L, calcN!, grid)
end
nothing #hide

# We now have all necessary building blocks to construct a `FourierFlows.Problem`.
# It would be useful, however, to define some more "helper functions". For example,
# a function that updates all variables given the solution `sol` which comprises ``\hat{u}``,
# ``\hat{v}`` and ``\hat{\eta}``:

"""
    updatevars!(prob)
Update the variables in `prob.vars` using the solution in `prob.sol`.
"""
function updatevars!(prob)
  vars, grid, sol = prob.vars, prob.grid, prob.sol

  @. vars.uh = sol[:, 1]
  @. vars.vh = sol[:, 2]
  @. vars.ηh = sol[:, 3]

  ldiv!(vars.u, grid.rfftplan, deepcopy(sol[:, 1])) # use deepcopy() because irfft destroys its input
  ldiv!(vars.v, grid.rfftplan, deepcopy(sol[:, 2])) # use deepcopy() because irfft destroys its input
  ldiv!(vars.η, grid.rfftplan, deepcopy(sol[:, 3])) # use deepcopy() because irfft destroys its input

  return nothing
end
nothing #hide

# Another useful function is one that prescribes an initial condition to the state variable `sol`.

"""
    set_uvη!(prob, u0, v0, η0)
Sets the state variable `prob.sol` as the Fourier transforms of `u0`, `v0`, and `η0`
and update all variables in `prob.vars`.
"""
function set_uvη!(prob, u0, v0, η0)
  vars, grid, sol = prob.vars, prob.grid, prob.sol

  A = typeof(vars.u) # determine the type of vars.u

  mul!(vars.uh, grid.rfftplan, A(u0)) # A(u0) converts u0 to the same type as vars expects (useful if u0 is a CPU array while working on the GPU)
  mul!(vars.vh, grid.rfftplan, A(v0)) # A(v0) converts u0 to the same type as vars expects (useful if v0 is a CPU array while working on the GPU)
  mul!(vars.ηh, grid.rfftplan, A(η0)) # A(η0) converts u0 to the same type as vars expects (useful if η0 is a CPU array while working on the GPU)

  @. sol[:, 1] = vars.uh
  @. sol[:, 2] = vars.vh
  @. sol[:, 3] = vars.ηh

  updatevars!(prob)

  return nothing
end
nothing #hide

# ## Let's prescibe parameter values and solve the PDE
#
# We are now ready to write up a program that sets up parameter values, constructs
# the problem `prob`, # time steps the solutions `prob.sol` and plots it.

# ## Choosing a device: CPU or GPU

dev = CPU()    # Device (CPU/GPU)
nothing # hide

# ## Numerical parameters and time-stepping parameters

     nx = 512            # grid resolution
stepper = "FilteredRK4"  # timestepper
     dt = 20.0           # timestep (s)
 nsteps = 320            # total number of time-steps
nothing # hide

# ## Physical parameters

Lx = 500e3      # Domain length (m)
g  = 9.8        # Gravitational acceleration (m s⁻²)
H  = 200.0      # Fluid depth (m)
f  = 1e-2       # Coriolis parameter (s⁻¹)
ν  = 100.0      # Viscosity (m² s⁻¹)
nν = 1          # Viscosity order (nν = 1 means Laplacian ∇²)
nothing # hide

# ## Construct the `struct`s and you are ready to go!
# Create a `grid` and also `params`, `vars`, and the `equation` structs. Then
# give them all as input to the `FourierFlows.Problem()` constructor to get a
# problem struct, `prob`, that contains all of the above.

    grid = OneDGrid(dev, nx, Lx)
  params = Params(ν, nν, g, H, f)
    vars = Vars(dev, grid)
equation = Equation(dev, params, grid)

    prob = FourierFlows.Problem(equation, stepper, dt, grid, vars, params, dev)
nothing #hide

# ## Setting initial conditions

# For initial condition we take the fluid at rest (``u = v = 0``). The free surface elevation
# is perturbed from its rest position (``\eta=0``); the disturbance we impose a Gaussian
# bump with half-width greater than the deformation radius and on top of that we
# superimpose some random noise with scales smaller than the deformation radius.
# We mask the small-scale perturbations so that it only applies in the central part
# of the domain by applying
#
# The system develops geostrophically-balanced jets around the Gaussian bump,
# while the smaller-scale noise propagates away as inertia-gravity waves.

# First let's construct the Gaussian bump.

gaussian_width = 6e3
gaussian_amplitude = 3.0
gaussian_bump = @. gaussian_amplitude * exp( - grid.x^2 / (2*gaussian_width^2) )

plot(grid.x/1e3, gaussian_bump,    # divide with 1e3 to convert m -> km
     color = :black,
    legend = false,
 linewidth = 2,
     alpha = 0.7,
     xlims = (-Lx/2e3, Lx/2e3),
    xlabel = "x [km]",
    ylabel = "η [m]",
     title = "A gaussian bump with half-width ≈ "*string(gaussian_width/1e3)*" km",
      size = (600, 260))

# Next the noisy perturbation. The `mask` is simply a product of hyperbolic tangent functions.
mask = @. 1/4 * (1 + tanh( -(grid.x - 100e3) / 10e3)) * (1 + tanh( (grid.x + 100e3) / 10e3))

noise_amplitude = 0.1 # the amplitude of the noise for η(x,t=0) (m)
η_noise = noise_amplitude * Random.randn(size(grid.x))
@. η_noise *= mask    # mask the noise

plot_noise = plot(grid.x/1e3, η_noise,      # divide with 1e3 to convert m -> km
                 color = :black,
                legend = :false,
             linewidth = [3 2],
                 alpha = 0.7,
                 xlims = (-Lx/2e3, Lx/2e3), # divide with 1e3 to convert m -> km
                 ylims = (-0.3, 0.3),
                xlabel = "x [km]",
                ylabel = "η [m]")

plot_mask = plot(grid.x/1e3, mask,          # divide with 1e3 to convert m -> km
                 color = :gray,
                legend = :false,
             linewidth = [3 2],
                 alpha = 0.7,
                 xlims = (-Lx/2e3, Lx/2e3), # divide with 1e3 to convert m -> km
                xlabel = "x [km]",
                ylabel = "mask")

title = plot(title = "Small-scale noise",
              grid = false,
          showaxis = false,
            xticks = [],
            yticks = [],
     bottom_margin = -20Plots.px)

plot(title, plot_noise, plot_mask,
           layout = @layout([A{0.01h}; [B; C]]),
             size = (600, 400))

# Sum the Gaussian bump and the noise and then call `set_uvη!()` to set the initial condition to the problem `prob`.

η0 = @. gaussian_bump + η_noise
u0 = zeros(grid.nx)
v0 = zeros(grid.nx)

set_uvη!(prob, u0, v0, η0)

plot(grid.x/1e3, η0,    # divide with 1e3 to convert m -> km
     color = :black,
    legend = false,
 linewidth = 2,
     alpha = 0.7,
     xlims = (-Lx/2e3, Lx/2e3),
    xlabel = "x [km]",
    ylabel = "η [m]",
     title = "initial surface elevation, η(x, t=0)",
      size = (600, 260))

# ## Visualizing the simulation

# We define a function that plots the surface elevation ``\eta`` and the
# depth-integrated velocities ``u`` and ``v``.

function plot_output(prob)
  plot_η = plot(grid.x/1e3, vars.η,         # divide with 1e3 to convert m -> km
                 color = :blue,
                legend = false,
             linewidth = 2,
                 alpha = 0.7,
                 xlims = (-Lx/2e3, Lx/2e3), # divide with 1e3 to convert m -> km
                xlabel = "x [km]",
                ylabel = "η [m]")

  plot_u = plot(grid.x/1e3, vars.u,         # divide with 1e3 to convert m -> km
                 color = :red,
                legend = false,
             linewidth = 2,
                 alpha = 0.7,
                 xlims = (-Lx/2e3, Lx/2e3), # divide with 1e3 to convert m -> km
                 ylims = (-0.3, 0.3),
                xlabel = "x [km]",
                ylabel = "u [m s⁻¹]")

  plot_v = plot(grid.x/1e3, vars.v,         # divide with 1e3 to convert m -> km
                 color = :green,
                legend = false,
             linewidth = 2,
                 alpha = 0.7,
                 xlims = (-Lx/2e3, Lx/2e3), # divide with 1e3 to convert m -> km
                 ylims = (-0.3, 0.3),
                xlabel = "x [km]",
                ylabel = "v [m s⁻¹]")

  Ld = @sprintf "%.2f" sqrt(g*H)/f /1e3     # divide with 1e3 to convert m -> km
  plottitle = "Deformation radius √(gh) / f = "*string(Ld)*" km"

  title = plot(title = plottitle,
                grid = false,
            showaxis = false,
              xticks = [],
              yticks = [],
       bottom_margin = -30Plots.px)

  return plot(title, plot_η, plot_u, plot_v,
           layout = @layout([A{0.01h}; [B; C; D]]),
             size = (600, 800))
end
nothing # hide

# ## Time-stepping the `Problem` forward

# We time-step the `Problem` forward in time. We update variables by calling
# `updatevars!()` and we also update the plot. We enclose the `for` loop in
# an `@animate` macro to produce an animation of the solution.

p = plot_output(prob)

anim = @animate for j = 0:nsteps
  updatevars!(prob)

  p[2][1][:y] = vars.η    # updates the plot for η
  p[2][:title] = "t = " * @sprintf("%.1f", prob.clock.t/60) * " min" # updates time in the title
  p[3][1][:y] = vars.u    # updates the plot for u
  p[4][1][:y] = vars.v    # updates the plot for v

  stepforward!(prob)
end

mp4(anim, "onedshallowwater.mp4", fps=18)

# ## Geostrophic balance

# It is instructive to compare the solution for ``v`` with its geostrophically balanced approximation, ``f \widehat{\bm{z}} \times \bm{u}_{\rm geostrophic} = - g \bm{\nabla} \eta``, i.e.,
#
# ```math
# \begin{aligned}
# v_{\rm geostrophic} & =   \frac{g}{f} \frac{\partial \eta}{\partial x} \ , \\
# u_{\rm geostrophic} & = - \frac{g}{f} \frac{\partial \eta}{\partial y} = 0 \ .
# \end{aligned}
# ```
# The geostrophic solution should capture well the the behavior of the flow in
# the center of the domain, after small-scale disturbances propagate away.

u_geostrophic = zeros(grid.nx)  # -g/f ∂η/∂y = 0
v_geostrophic = params.g / params.f * irfft(im * grid.kr .* vars.ηh, grid.nx)  #g/f ∂η/∂x

plot_u = plot(grid.x/1e3, [vars.u u_geostrophic], # divide with 1e3 to convert m -> km
                 color = [:red :purple],
                labels = ["u" "- g/f ∂η/∂y"],
             linewidth = [3 2],
                 alpha = 0.7,
                 xlims = (-Lx/2e3, Lx/2e3),       # divide with 1e3 to convert m -> km
                 ylims = (-0.3, 0.3),
                xlabel = "x [km]",
                ylabel = "u [m s⁻¹]")

plot_v = plot(grid.x/1e3, [vars.v v_geostrophic], # divide with 1e3 to convert m -> km
                 color = [:green :purple],
                labels = ["v" "g/f ∂η/∂x"],
             linewidth = [3 2],
                 alpha = 0.7,
                 xlims = (-Lx/2e3, Lx/2e3),       # divide with 1e3 to convert m -> km
                 ylims = (-0.3, 0.3),
                xlabel = "x [km]",
                ylabel = "v [m s⁻¹]")

title = plot(title = "Geostrophic balance",
              grid = false,
          showaxis = false,
            xticks = [],
            yticks = [],
     bottom_margin = -20Plots.px)

plot(title, plot_u, plot_v,
           layout = @layout([A{0.01h}; [B; C]]),
             size = (600, 400))


text:
@doc raw"""Legendre basis

In the original paper [L. Boehnke et al., PRB 84, 075145 (2011)],
they used:

    G(\tau) = \sum_{l=0} \sqrt{2l+1} P_l[x(\tau)] G_l/beta,

where P_l[x] is the $l$-th Legendre polynomial.

In this type, the basis functions are defined by

    U_l(\tau) \equiv c_l (\sqrt{2l+1}/beta) * P_l[x(\tau)],

where c_l are additional l-dependent constant factors.
By default, we take c_l = 1, which reduces to the original definition.
"""
struct LegendreBasis{T<:AbstractFloat} <: AbstractBasis
    statistics::Statistics
    β::Float64
    cl::Vector{T}
    u::PiecewiseLegendrePolyVector{T}
    uhat::PiecewiseLegendreFTVector{T}
end

function LegendreBasis(
    statistics::Statistics,
    beta::Float64,
    size::Int;
    cl::Vector{Float64}=ones(Float64, size),
)
    beta > 0 || throw(DomainError(beta, "inverse temperature beta must be positive"))
    size > 0 || throw(DomainError(size, "size of basis must be positive"))

    # u
    knots = Float64[0, beta]
    data = zeros(Float64, size, length(knots) - 1, size)
    symm = (-1) .^ collect(0:(size - 1))
    for l in 1:size
        data[l, 1, l] = sqrt(((l - 1) + 0.5) / beta) * cl[l]
    end
    u = PiecewiseLegendrePolyVector(data, knots; symm)

    # uhat
    uhat_base = PiecewiseLegendrePolyVector(sqrt(beta) .* data, Float64[-1, 1]; symm)
    even_odd = Dict(fermion => :odd, boson => :even)[statistics]
    uhat = hat.(uhat_base, even_odd)

    return LegendreBasis(statistics, beta, cl, u, uhat)
end

function Base.getproperty(obj::LegendreBasis, d::Symbol)
    if d === :v
        return nothing
    else
        return getfield(obj, d)
    end
end

iswellconditioned(basis::LegendreBasis) = true

function default_tau_sampling_points(basis::LegendreBasis)
    x = gauss(length(basis.u))[1]
    return (getbeta(basis) / 2) .* (x .+ 1)
end

struct _ConstTerm{T<:Number}
    value::T
end

(ct::_ConstTerm)(n::AbstractVector{<:Integer}) = fill(ct.value, (1, length(n)))
(ct::_ConstTerm)(n::Integer) = ct([n])

"""
Constant term in matsubara-frequency domain
"""
struct MatsubaraConstBasis{T<:AbstractFloat} <: AbstractBasis
    statistics::Statistics
    β::Float64
    uhat::_ConstTerm{T}
end

function MatsubaraConstBasis(statistics::Statistics, beta::Float64; value=1)
    beta > 0 || throw(DomainError(beta, "inverse temperature beta must be positive"))
    return MatsubaraConstBasis(statistics, beta, _ConstTerm(value))
end

Base.size(::MatsubaraConstBasis) = (1,)


text:
function solve()
    x, k, d = [parse(Int, x) for x in split(readline())]
    x = abs(x)
    t = x ÷ d
    if t > k
        x - k * d
    else
        k = k - t
        x = x - t * d
        if k % 2 == 0
            x
        else
            d - x
        end
    end
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
using Documenter, Gadfly, Miletus, GraphPlot

load_dir(x) = map(file -> joinpath("lib", x, file), readdir(joinpath(Base.source_dir(), "src", "lib", x)))

makedocs(
    modules = [Miletus],
    clean = false,
    format = [:html, :latex],
    sitename = "JuliaFin∕Miletus",
    pages = Any[
        "Introduction" => "index.md",
        "Tutorial" => "tutorial.md",
        "Examples" => "examples.md"
    ],
    assets = ["assets/jc.css"]
)

deploydocs(
    repo   = "github.com/JuliaComputing/Miletus.jl.git",
    julia  = "0.5",
    osname = "linux",
    deps = nothing,
    make = nothing,
    target = "build",
)


text:
function sieve(n :: Int)
    isprime = trues(n)
    isprime[1] = false
    for p in 2:n
        if isprime[p]
            j = p * p
            if j > n
                return findall(isprime)
            else
                for k in j:p:n
                  isprime[k] = false
                end
            end
        end
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# module for estimating pose by extended kalman filter
# initial pose is decided randomly
# global localization problem

using Distributions, LinearAlgebra, StatsBase

include(joinpath(split(@__FILE__, "src")[1], "src/model/map/map.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/common/covariance_ellipse/covariance_ellipse.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/common/state_transition/state_transition.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/common/observation_function/observation_function.jl"))

mutable struct GlobalKf
  belief
  motion_noise_stds
  estimated_pose
  estimated_cov
  map
  dist_dev
  dir_dev

  function GlobalKf(init_pose::Array;
                    motion_noise_stds::Dict=Dict("nn"=>0.20, "no"=>0.001, "on"=>0.11, "oo"=>0.20),
                    env_map=nothing, dist_dev_rate=0.14, dir_dev=0.05)
    self = new()
    self.belief = MvNormal([rand(Uniform(-5.0, 5.0)), rand(Uniform(-5.0, 5.0)), rand(Uniform(-pi, pi))],
                           diagm(0 => [1e+4, 1e+4, 1e+4]))
    self.motion_noise_stds = motion_noise_stds
    self.estimated_pose = self.belief.μ
    self.estimated_cov = self.belief.Σ
    self.map = env_map
    self.dist_dev = dist_dev_rate
    self.dir_dev = dir_dev
    return self
  end
end

function mat_M(speed, yaw_rate, time, stds)
  return diagm(0 => [stds["nn"]^2*abs(speed)/time + stds["no"]^2*abs(yaw_rate)/time,
                     stds["on"]^2*abs(speed)/time + stds["oo"]^2*abs(yaw_rate)/time])
end

function mat_A(speed, yaw_rate, time, theta)
  st, ct = sin(theta), cos(theta)
  stw, ctw = sin(theta + yaw_rate * time), cos(theta + yaw_rate * time)
  return [(stw - st)/yaw_rate  -speed/(yaw_rate^2)*(stw - st) + speed/yaw_rate*time*ctw;
          (-ctw + ct)/yaw_rate -speed/(yaw_rate^2)*(-ctw + ct) + speed/yaw_rate*time*stw;
          0                    time]
end

function mat_F(speed, yaw_rate, time, theta)
  F = diagm(0 => [1.0, 1.0, 1.0])
  F[1, 3] = speed / yaw_rate * (cos(theta + yaw_rate * time) - cos(theta))
  F[2, 3] = speed / yaw_rate * (sin(theta + yaw_rate * time) - sin(theta))
  return F
end

function mat_H(mu_pose, obj_pose)
  obj_x, obj_y = obj_pose[1], obj_pose[2]
  mu_x, mu_y = mu_pose[1], mu_pose[2]
  mu_l = sqrt((mu_x - obj_x)^2 + (mu_y - obj_y)^2)
  return [(mu_x - obj_x)/mu_l (mu_y - obj_y)/mu_l 0.0;
          (obj_y - mu_y)/(mu_l^2) (mu_x - obj_x)/(mu_l^2) -1.0]
end

function mat_Q(dist_dev, dir_dev)
  return [dist_dev^2 0.0;
          0.0 dir_dev^2]
end

function motion_update(self::GlobalKf, speed,
                       yaw_rate, time)
  if abs(yaw_rate) < 1e-5 # to prevent division by zero
    yaw_rate = 1e-5
  end
  M = mat_M(speed, yaw_rate, time, self.motion_noise_stds)
  A = mat_A(speed, yaw_rate, time, self.estimated_pose[3])
  F = mat_F(speed, yaw_rate, time, self.estimated_pose[3])
  self.estimated_cov = F * self.estimated_cov * F' + A * M * A'
  self.estimated_pose = state_transition(speed, yaw_rate, time, self.estimated_pose)
end

function observation_update(self::GlobalKf, observation)
  for obs in observation
    z = obs[1] # [distance, direction]
    id = obs[2]

    H = mat_H(self.estimated_pose, self.map.objects[id].pose)
    estimated_z = observation_function(self.estimated_pose,
                                       self.map.objects[id].pose)

    # observation noise
    Q = mat_Q(estimated_z[1]*self.dist_dev, self.dir_dev)

    # kalman gain
    K = self.estimated_cov * H' * inv(Q + H*self.estimated_cov*H')
    self.estimated_pose += K * (z - estimated_z)
    self.estimated_cov = (Matrix{Float64}(I, 3, 3) - K*H) * self.estimated_cov
  end
end

function draw!(self::GlobalKf)
  draw_covariance_ellipse!(self.estimated_pose, self.estimated_cov, 3)
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
#MDDatasets DataRS (Recursive Sweep) defninitions
#-------------------------------------------------------------------------------

#==Main types
===============================================================================#

@doc """
    DataRS{T} <: DataMD

Recursive data structure used to store results from a parametric sweep.

See also: [fill(::DataRS, ...)](@ref) [ndims(::DataRS)](@ref) [paramlist(::DataRS)](@ref)
""" DataRS

#Linked-list representation of multi-dimensional datasets:
#-------------------------------------------------------------------------------
mutable struct DataRS{T} <: DataMD
	sweep::PSweep
	elem::Vector{T}

	function DataRS{T}(sweep::PSweep, elem::Vector{T}) where {T}
		if !elemallowed(DataRS, eltype(elem))
			msg = "Can only create DataRS{T} for T ∈ {DataRS, DataF1, DataFloat, DataInt, DataComplex}"
			throw(ArgumentError(msg))
		elseif length(sweep) != length(elem)
			throw(ArgumentError("sweep length does not match number of elem"))
		end
		return new(sweep, elem)
	end
end

#Shorthand (because default (non-parameterized) constructor was overwritten):
DataRS(sweep::PSweep, elem::Vector{T}) where T = DataRS{T}(sweep, elem)

elemallowed(::Type{DataRS}, t::Type{T}) where T = elemallowed(DataMD, t) #Allow basic types
elemallowed(::Type{DataRS}, ::Type{DataRS}) = true #Also allow recursive structures

#Generate empty DataRS structure:
(::Type{DataRS{T}})(sweep::PSweep) where T = DataRS{T}(sweep, Array{T}(undef, length(sweep)))

#==Type promotions
===============================================================================#
Base.promote_rule(::Type{T1}, ::Type{T2}) where {T1<:DataRS, T2<:Number} = DataRS

#==Unchecked accessors
===============================================================================#
#Get expected # of dimensions by walking across first path of ::DataRS
function _ndims_p1(d::DataRS)
	depth = 0

	#Count depth:
	w = d #walks across data structure
	while isa(w, DataRS)
		depth += 1
		w = w.elem[1]
	end
	return depth
end

#Get expected leaf element by walking across first path of ::DataRS
function _leaftype_p1(d::DataRS)
	w = d #walks across data structure
	leaftype = DataRS
	while isa(w, DataRS)
		leaftype = eltype(w)
		w = w.elem[1]
	end
	return leaftype
end

#Get a list of sweeps along the first path of a ::DataRS.
function _getsweeplist_p1(d::DataRS)
	sweeplist = PSweep[]
	w = d #walks across data structure
	while isa(w, DataRS)
		push!(sweeplist, w.sweep)
		w = w.elem[1]
	end
	return sweeplist
end

#==Data integrity
===============================================================================#
"""
    sweepsmatch(s1, s2)

Return true if two sweeps match.

TODO: Extend Base.(=)() instead??
"""
sweepsmatch(s1::PSweep, s2::PSweep) = (s1.id==s2.id && s1.v==s2.v)
#NOTE: fill!() creates different sweep objects for each pass (even if values are the same)

function ensure_validsweep(d::DataRS)
	if length(d.elem) != length(d.sweep)
		id = d.sweep.id
		msg = "Sweep length : data length mismatch: $id"
		throw(ArgumentError(msg))
	end
end

function _validate_dimensionality(d::DataRS, depth::Int)
	function validate_depth(d, depth) #Leaf element
		ensure(0 == depth, ArgumentError("DataRS: Inconsistent dimensionality"))
	end

	function validate_depth(d::DataRS, depth)
		ensure(depth > 0, ArgumentError("DataRS: Inconsistent dimensionality")) #Recursion safety
		ensure_validsweep(d)
		nextdepth = depth - 1
		for e in d.elem
			validate_depth(e, nextdepth)
		end
		return
	end
	return validate_depth(d, depth)
end
validate_dimensionality(d::DataRS) = _validate_dimensionality(d, _ndims_p1(d))

#Ensure that all sweeps @depth are equal to psweep:
function _validate_paramsweep(d::DataRS, sweeplist::Vector{PSweep})
	function validate_sweep(d::DataRS{T}, sweeplist) where T
		ensure(length(sweeplist)>0,
			ArgumentError("More sweeps than expected: $(d.sweep)")
		)
		ensure(sweepsmatch(d.sweep, sweeplist[1]),
			ArgumentError("Unexpected sweep: $(d.sweep). Expected: $(sweeplist[1]).")
		)
		if !(T<:DataRS) #No more sweeps to check
			ensure(length(sweeplist)<2,
				ArgumentError("Missing sweeps:\n $(sweeplist[2:end])")
			)
			return #Stop recursion
		end
		ensure(length(sweeplist)>1,
			ArgumentError("More sweeps than expected. Should stop at: $(sweeplist[1]).")
		)
		for e in d.elem
			validate_sweep(e, sweeplist[2:end])
		end
		return
	end

	return validate_sweep(d, sweeplist)
end

#==Accessor functions
===============================================================================#
Base.eltype(d::DataRS{T}) where T = T
Base.length(d::DataRS) = length(d.elem)

"""
    ndims(d::DataRS)

Return number of dimensions for the parametric sweep in d.
""" Base.ndims(::DataRS)
function Base.ndims(d::DataRS)
	depth = _ndims_p1(d)
	_validate_dimensionality(d, depth)
	return depth
end

"""
    sweeps(d::DataRS)

Return a list of parameter sweeps for d.
"""
function sweeps(d::DataRS)
	ref = _getsweeplist_p1(d)
	_validate_paramsweep(d, ref)
	return ref
end

"""
    paramlist(d::DataRS)

Return a list of parameter values being swept.
"""
function paramlist(d::DataRS)
	ref = _getsweeplist_p1(d)
	_validate_paramsweep(d, ref)
	return paramlist(ref)
end

#==Help with construction
===============================================================================#

#Implement "fill(DataRS, ...) do sweepval" syntax:
function Base.fill!(fn::Function, d::DataRS)
	for i in 1:length(d.sweep)
		d.elem[i] = fn(d.sweep.v[i])
	end
	return d
end
Base.fill(fn::Function, ::Type{DataRS{T}}, sweep::PSweep) where T =
	fill!(fn, DataRS(sweep, Array{T}(undef, length(sweep))))
Base.fill(fn::Function, ::Type{DataRS}, sweep::PSweep) = fill(fn, DataRS{DataRS}, sweep)

#Fill with a specified value (leverage DataHR algorithms for now):
Base.fill(elem_0::DF1_Num, ::Type{DataRS}, swlist::Vector{PSweep}) =
	convert(DataRS, fill(elem_0, DataHR, swlist))
Base.fill(elem_0::DF1_Num, ::Type{DataRS}, sweep::PSweep) = fill(elem_0, DataRS, PSweep[sweep])

#zeros & ones:
Base.zeros(::Type{DataRS}, swlist::Vector{PSweep}) = fill(0.0, DataRS, swlist)
Base.zeros(::Type{DataRS}, sweep::PSweep) = zeros(DataRS, PSweep[sweep])
Base.ones(::Type{DataRS}, swlist::Vector{PSweep}) = fill(1.0, DataRS, swlist)
Base.ones(::Type{DataRS}, sweep::PSweep) = ones(DataRS, PSweep[sweep])

#==Data generation
===============================================================================#
function _ensuresweepunique(d::DataRS, sweepid::String)
	if sweepid == d.sweep.id
		msg = "Sweep occurs multiple times in DataRS: $sweepid"
		throw(ArgumentError(msg))
	end
end

#Define "parameter".
#(Generates a DataRS object containing the value of a given swept parameter)
#-------------------------------------------------------------------------------

@doc """
    parameter(d::DataRS, sweepid::String)

Get parameter values for a particular sweep.
""" parameter

#Deal with non-leaf elements, once the sweep value is found:
function _parameter(d::DataRS{DataRS}, sweepid::String, sweepval::T) where T
	_ensuresweepunique(d, sweepid)
	elem = DataRS[_parameter(d.elem[i], sweepid, sweepval) for i in 1:length(d.sweep)]
	return DataRS(d.sweep, elem)
end

#Deal with leaf elements, once the sweep value is found:
function _parameter(d::DataRS, sweepid::String, sweepval::T) where T
	_ensuresweepunique(d, sweepid)
	elem = T[sweepval for i in 1:length(d.sweep)]
	return DataRS(d.sweep, elem)
end

#Main "parameter" algorithm (non-leaf elements):
function parameter(d::DataRS{DataRS}, sweepid::String)
	if sweepid == d.sweep.id #Sweep found
		elem = DataRS[_parameter(d.elem[i], sweepid, d.sweep.v[i]) for i in 1:length(d.sweep)]
	else
		elem = DataRS[parameter(d.elem[i], sweepid) for i in 1:length(d.sweep)]
	end
	return DataRS(d.sweep, elem)
end
#Main "parameter" algorithm (leaf elements):
function parameter(d::DataRS, sweepid::String)
	T = eltype(d.sweep.v)
	if sweepid == d.sweep.id #Sweep found
		return DataRS(d.sweep, d.sweep.v)
	else
		msg = "Sweep not found in DataRS: $sweepid"
		throw(ArgumentError(msg))
	end
end

#Generate DataRS from DataHR.
#-------------------------------------------------------------------------------
function _buildDataRS(d::DataHR, firstinds::Vector{Int})
	curidx = length(firstinds) + 1
	sweep = d.sweeps[curidx]
	if curidx < length(d.sweeps)
		result = DataRS{DataRS}(sweep)
		for i in 1:length(sweep.v)
			result.elem[i] = _buildDataRS(d, vcat(firstinds, i))
		end
	else #Last index.  Copy data over:
		result = DataRS{eltype(d.elem)}(sweep)
		for i in 1:length(sweep.v)
			result.elem[i] = d.elem[firstinds..., i]
		end
	end
	return result
end
_buildDataRS(d::DataHR) = _buildDataRS(d, Int[]) #Start of recursive algorithm

#Generate DataHR from DataRS.
#-------------------------------------------------------------------------------
function _buildDataHR(src::DataRS)
	function _build(data, idxlist, depth, src::DataRS{DataHR}) #Corrupt structure
		idxlist = String[v in idxlist]
		for i in depth:length(idxlist)
			stridx = "?"
		end
		idxstr = join(idxlist, " ")
		throw(ArgumentError("Found DataHR @ index [$idxstr]"))
	end
	function _build(data, idxlist, depth, src::DataRS) #Leaf element
		for (i, elem) in enumerate(src.elem)
			idxlist[depth] = i
			data[idxlist...] = elem #TODO: copy???
		end
	end
	function _build(data, idxlist, depth, src::DataRS{DataRS})
		for (i, elem) in enumerate(src.elem)
			idxlist[depth] = i
			_build(data, idxlist, depth+1, elem)
		end
	end

	LT = _leaftype_p1(src)
	sweeplist = sweeps(src) #Validates that structure is compatible with DataHR
	idxlist = ones(Int, length(sweeplist))
	asize = size(DataHR, sweeplist)
	data = Array{LT}(undef, asize)
		_build(data, idxlist, 1, src)
	return DataHR(sweeplist, data)
end

Base.convert(::Type{DataRS}, d::DataHR) = _buildDataRS(d)
Base.convert(::Type{DataHR}, d::DataRS) = _buildDataHR(d)

@deprecate DataRS(d::DataHR) convert(DataRS, d)

#==User-friendly show functions
===============================================================================#

#Print leaf element:
function printDataRSelem(io::IO, ds::DataRS, idx::Int, indent::String)
	if isassigned(ds.elem, idx)
		println(io, ds.elem[idx])
	else
		println(io, indent, "UNDEFINED")
	end
end
#Print next level of recursive DataRS:
function printDataRSelem(io::IO, ds::DataRS{T}, idx::Int, indent::String) where T<:DataRS
	println(io)
	if isassigned(ds.elem, idx)
		printDataRS(io, ds.elem[idx], indent)
	else
		println(io, indent, "UNDEFINED")
	end
end
#Print DataRS structure:
function printDataRS(io::IO, ds::DataRS, indent::String)
	for i in 1:length(ds.elem)
		print(io, "$indent", ds.sweep.id, "=", ds.sweep.v[i], ": ")
		printDataRSelem(io, ds, i, "$indent  ")
	end
end

function Base.show(io::IO, ds::DataRS)
	print(io, "DataRS[\n")
	printDataRS(io, ds, "  ")
	print(io, "]\n")
end


text:
module VisualParticipationAnalytics

using Clustering
using DataFrames
using Distances
using JSON3
using Languages
using LinearAlgebra
using SparseArrays
using SQLite
using Tables
using TextAnalysis

export clustering,
        corpus,
        preprocess,
        readdb,
        save,
        similarities,
        slatetext,
        topicmodel,
        topkwords

include("io.jl")
include("preprocessing.jl")

""""""
function docterm(crps::Corpus, tfidf::Bool)
    update_lexicon!(crps)
    dt = DocumentTermMatrix(crps)
    m = dtm(dt)
    if tfidf
        m = tf_idf(m)
    end
    return m
end

""""""
function similarities(crps::Corpus, tfidf=true)
    dt = docterm(crps, tfidf)
    return sparse(1 .- pairwise(CosineDist(), Matrix(dt'), dims=2))
end

""""""
function clustering(crps::Corpus, k::Int, tfidf=true)
    dt = docterm(crps, tfidf)
    return kmeans(Matrix(dt'), k; init=:kmpp)
end

""""""
function topicmodel(crps::Corpus, k::Int, iter::Int, alpha::Float64,
                beta::Float64)
    update_lexicon!(crps)
    return lda(DocumentTermMatrix(crps), k, iter, alpha, beta)
end

""""""
function topkwords(topicword::SparseMatrixCSC{Float64,Int}, i::Int,
                    crps::Corpus, k::Int=10)
    terms = collect(keys(lexicon(crps)))
    topicword = permutedims(topicword)
    rows = rowvals(topicword)
    vals = nonzeros(topicword)
    wordprobs = Array{Tuple{String,Float64},1}()
    for j in nzrange(topicword, i)
        push!(wordprobs, (terms[rows[j]], vals[j]))
    end
    wordprobs = sort(wordprobs; by=w -> w[2], rev=true)
    n = k < length(wordprobs) ? k : length(wordprobs)
    return wordprobs[1:n]
end

end # module VisualParticipationAnalytics


text:
using Manopt, Manifolds, Plots
#
# Settings
resultsFolder = "examples/Total_Variation/S1_TV/"
experimentName = "S1_TV12"
plotAndExportResult = true

n = 500
σ = 0.2
α = 0.5
β = 1.0

dataColor = RGBA{Float64}(colorant"#BBBBBB")
s2dColor = RGBA{Float64}(colorant"#EE7733") # data Color: Tol Vibrant Orange
s1Color = RGBA{Float64}(colorant"#0077BB") # control point data color: Tol Virbant Blue
nColor = RGBA{Float64}(colorant"#33BBEE") # tangent vector: Tol Vibrant Teal

if !isdir(resultsFolder)
    mkdir(resultsFolder)
end
#
# Manifolds and Data
M = Circle()
N = PowerManifold(M, n)
f = artificial_S1_signal(n)
xCompare = f
fn = exp.(Ref(M),f, random_tangent.(Ref(M),f, Val(:Gaussian), σ))
data = fn
t = range(0.0, 1.0, length = n)

if plotAndExportResult
    scene = scatter(
        t,
        f,
        markersize = 2,
        markercolor = dataColor,
        markerstrokecolor = dataColor,
        lab = "original",
    )
    scatter!(
        scene,
        t,
        fn,
        markersize = 2,
        markercolor = nColor,
        markerstrokecolor = nColor,
        lab = "noisy",
    )
    yticks!(
        [-π, -π/2, 0, π/2, π],
        [raw"$-\pi$",raw"$-\frac{\pi}{2}$",raw"$0$",raw"$\frac{\pi}{2}$",raw"$\pi$"]
    )
    png(scene, "$(resultsFolder)$(experimentName)-original.png")
end
#
# Setup and Optimize
F = x -> costL2TVTV2(N, data, α, β, x)
proxes = [
    (λ, x) -> prox_distance(N, λ, data, x,2),
    (λ, x) -> prox_TV(N, α * λ, x),
    (λ, x) -> prox_TV2(N, β * λ, x),
]

o = cyclic_proximal_point(
    N,
    F,
    proxes,
    data;
    λ = i -> π / (2*i),
    debug = Dict(
        :Stop => DebugStoppingCriterion(),
        :Step => DebugEvery(
            DebugGroup([
                DebugIteration(),
                DebugDivider(),
                DebugProximalParameter(),
                DebugDivider(),
                DebugCost(),
                DebugDivider(),
                DebugChange(),
                DebugDivider("\n"),
            ]),
            1000,
        ),
        :Start => DebugDivider("Starting the solver\n"),
    ),
    record = [:Iteration, :Cost, :Change, :Iterate],
    return_options = true,
)
fR = get_solver_result(o)
r = get_record(o)
#
# Result
if plotAndExportResult
    scene = scatter(
        t,
        f,
        markersize = 2,
        markercolor = dataColor,
        markerstrokecolor = dataColor,
        lab = "original",
    )
    scatter!(
        scene,
        t,
        fR,
        markersize = 2,
        markercolor = nColor,
        markerstrokecolor = nColor,
        lab = "reconstruction",
    )
    yticks!(
        [-π, -π/2, 0, π/2, π],
        [raw"$-\pi$",raw"$-\frac{\pi}{2}$",raw"$0$",raw"$\frac{\pi}{2}$",raw"$\pi$"]
    )
    png(scene, "$(resultsFolder)$(experimentName)-result.png")
end

print("MSE (input):  ", 1 / n * distance(N, xCompare, data)^2, "\n")
print("MSE (result): ", 1 / n * distance(N, xCompare, fR)^2, "\n")


text:
"""
    mutable struct CategoricalSummarizer{T} <: Summarizer

Get unique values of a single categorical variable of type T.

Results will be stored with the key of `:unique` in the schema metadata.
"""
mutable struct CategoricalSummarizer{T} <: Summarizer
    name::Symbol
    seen::Set{T}
    uniq::Vector{T}
end

CategoricalSummarizer{T}(name::Symbol) where T = CategoricalSummarizer(name, Set{T}(), T[])

function update!(cs::CategoricalSummarizer{T}, nt) where T
    x = nt[cs.name]::T
    if x ∉ cs.seen
        push!(cs.seen, x)
        push!(cs.uniq, x)
    end
    cs
end

"""
    get(cs::CategoricalSummarizer{T}) = :unique => T[values...]
"""
Base.get(cs::CategoricalSummarizer) = :unique => cs.uniq

Base.hash(cs::CategoricalSummarizer, h::UInt64) = hash(typeof(cs), hash(cs.name, h))
Base.:(==)(c1::CategoricalSummarizer, c2::CategoricalSummarizer) = c1.name == c2.name

is_categorical(s, sch::Data.Schema) = is_categorical(string(s), sch)
is_categorical(s::String, sch::Data.Schema) = is_categorical(Data.types(sch)[sch[s]])
is_categorical(::Type{<:Real}) = false
is_categorical(::Type{Nullable{<:Real}}) = false
is_categorical(::Type) = true


text:
using MicrostructureNoise, Distributions
using Plots

using Statistics, DelimitedFiles, Dates
# uncomment if you do not mind to create this large file
# Base.download("https://www.truefx.com/dev/data//2015/MARCH-2015/EURUSD-2015-03.zip","./data/EURUSD-2015-03.zip")
# run(`unzip ./data/EURUSD-2015-03.zip -d ./data`)
dat = readdlm("./data/EURUSD-2015-03.csv", ',')
times = map(a -> DateTime(a, "yyyymmdd H:M:S.s"), dat[1:10:130260,2])
tt = Float64[1/1000*(times[i].instant.periods.value-times[1].instant.periods.value) for i in 1:length(times)]
n = length(tt)-1
T = tt[end]
y = Float64.(dat[1:10:130260, 3])

prior = MicrostructureNoise.Prior(
N = 40,

α1 = 0.0,
β1 = 0.0,

αη = 0.3,
βη = 0.3,

Πα = LogNormal(1., 0.5),
μ0 = 0.0,
C0 = 5.0
)

α = 0.3
σα = 0.1
td, θs, ηs, αs, p = MicrostructureNoise.MCMC(prior, tt, y, α, σα, 1500)

posterior = MicrostructureNoise.posterior_volatility(td, θs)

plot(MicrostructureNoise.piecewise(posterior.post_t, posterior.post_mean_root)..., label="Volatility estimate")


text:
module TokamakIterCylindrical

    import ElectromagneticFields.AxisymmetricTokamakCylindrical

    export charged_particle_3d_pode, charged_particle_3d_iode,
           hamiltonian, toroidal_momentum

    AxisymmetricTokamakCylindrical.@code_iter() # inject magnetic field code

    include("charged_particle_3d_canonical.jl")

    const qᵢ = [7.0, 0.0, 0.0]
    const vᵢ = [3.43E-3, 6.75, -3.41E-1]
    const pᵢ = charged_particle_3d_pᵢ(qᵢ, vᵢ)

end