text:
const DOCDIR = joinpath(Pkg.dir("ExamplePlots"), "docs", "examples")
const IMGDIR = joinpath(DOCDIR, "img")

"""
Holds all data needed for a documentation example... header, description, and plotting expression (Expr)
"""
type PlotExample
  header::AbstractString
  desc::AbstractString
  exprs::Vector{Expr}
end

function pretty_print_expr(io::IO, expr::Expr)
    lines = split(string(expr), "\n")
    start_idx = 1
    if expr.head == :block
        lines = lines[2:end-1]
        start_idx = 5
    end
    lines_without_comments = map(str -> split(str,"#")[1][start_idx:end], lines)
    for line in lines_without_comments
        println(io, expr)
    end
end

# the _examples we'll run for each
const _examples = PlotExample[
  PlotExample("Lines",
              "A simple line plot of the columns.",
              [
                :(plot(Plots.fakedata(50,5), w=3))
              ]),
  PlotExample("Functions, adding data, and animations",
              "Plot multiple functions.  You can also put the function first, or use the form `plot(f, xmin, xmax)` where f is a Function or AbstractVector{Function}.\n\nGet series data: `x, y = plt[i]`.  Set series data: `plt[i] = (x,y)`. Add to the series with `push!`/`append!`.\n\nEasily build animations.  (`convert` or `ffmpeg` must be available to generate the animation.)  Use command `gif(anim, filename, fps=15)` to save the animation.",
              [
                :(p = plot([sin,cos], zeros(0), leg=false)),
                :(anim = Animation()),
                :(for x in linspace(0, 10π, 100)
                    push!(p, x, Float64[sin(x), cos(x)])
                    frame(anim)
                  end)
              ]),
  PlotExample("Parametric plots",
              "Plot function pair (x(u), y(u)).",
              [
                :(plot(sin, x->sin(2x), 0, 2π, line=4, leg=false, fill=(0,:orange)))
              ]),
  PlotExample("Colors",
              "Access predefined palettes (or build your own with the `colorscheme` method).  Line/marker colors are auto-generated from the plot's palette, unless overridden.  Set the `z` argument to turn on series gradients.",
              [
                :(y = rand(100)),
                :(plot(0:10:100,rand(11,4),lab="lines",w=3,palette=:grays,fill=(0,:auto), α=0.6)),
                :(scatter!(y, zcolor=abs(y-.5), m=(:heat,0.8,stroke(1,:green)), ms=10*abs(y-0.5)+4, lab="grad"))
              ]),
  PlotExample("Global",
              "Change the guides/background/limits/ticks.  Convenience args `xaxis` and `yaxis` allow you to pass a tuple or value which will be mapped to the relevant args automatically.  The `xaxis` below will be replaced with `xlabel` and `xlims` args automatically during the preprocessing step. You can also use shorthand functions: `title!`, `xaxis!`, `yaxis!`, `xlabel!`, `ylabel!`, `xlims!`, `ylims!`, `xticks!`, `yticks!`",
              [
                :(y = rand(20,3)),
                :(plot(y, xaxis=("XLABEL",(-5,30),0:2:20,:flip), background_color = RGB(0.2,0.2,0.2), leg=false)),
                :(hline!(mean(y,1)+rand(1,3), line=(4,:dash,0.6,[:lightgreen :green :darkgreen]))),
                :(vline!([5,10])),
                :(title!("TITLE")),
                :(yaxis!("YLABEL", :log10))
              ]),
  PlotExample("Two-axis",
              "Use the `axis` arguments.",
              [
                :(plot(Vector[randn(100), randn(100)*100], axis = [:l :r], ylabel="LEFT", yrightlabel="RIGHT", xlabel="X", title="TITLE"))
              ]),
  PlotExample("Arguments",
              "Plot multiple series with different numbers of points.  Mix arguments that apply to all series (marker/markersize) with arguments unique to each series (colors).  Special arguments `line`, `marker`, and `fill` will automatically figure out what arguments to set (for example, we are setting the `linestyle`, `linewidth`, and `color` arguments with `line`.)  Note that we pass a matrix of colors, and this applies the colors to each series.",
              [
                :(ys = Vector[rand(10), rand(20)]),
                :(plot(ys, color=[:black :orange], line=(:dot,4), marker=([:hex :d],12,0.8,stroke(3,:gray))))
              ]),
  PlotExample("Build plot in pieces",
              "Start with a base plot...",
              [
                :(plot(rand(100)/3, reg=true, fill=(0,:green)))
              ]),
  PlotExample("",
              "and add to it later.",
              [
                :(scatter!(rand(100), markersize=6, c=:orange))
              ]),
  PlotExample("Histogram2D",
              "",
              [
                :(histogram2d(randn(10000), randn(10000), nbins=20))
              ]),
  PlotExample("Line types",
              "",
              [
                # :(types = intersect(supportedTypes(), [:line, :path, :steppre, :steppost, :sticks, :scatter])'),
                :(linetypes = [:path :steppre :steppost :sticks :scatter]),
                :(n = length(linetypes)),
                :(x = Vector[sort(rand(20)) for i in 1:n]),
                :(y = rand(20,n)),
                :(plot(x, y, line=(linetypes,3), lab=map(string,linetypes), ms=15))
              ]),
  PlotExample("Line styles",
              "",
              [
                :(styles = setdiff(supportedStyles(), [:auto])'),
                :(plot(cumsum(randn(20,length(styles)),1), style=:auto, label=map(string,styles), w=5))
              ]),
  PlotExample("Marker types",
              "",
              [
                :(markers = setdiff(supportedMarkers(), [:none,:auto,Shape])'),
                :(n = length(markers)),
                :(x = linspace(0,10,n+2)[2:end-1]),
                :(y = repmat(reverse(x)', n, 1)),
                :(scatter(x, y, m=(8,:auto), lab=map(string,markers), bg=:linen, xlim=(0,10), ylim=(0,10)))
              ]),
  PlotExample("Bar",
              "x is the midpoint of the bar. (todo: allow passing of edges instead of midpoints)",
              [
                :(bar(randn(99)))
              ]),
  PlotExample("Histogram",
              "",
              [
                :(histogram(randn(1000), nbins=20))
              ]),
  PlotExample("Subplots",
              """
                subplot and subplot! are distinct commands which create many plots and add series to them in a circular fashion.
                You can define the layout with keyword params... either set the number of plots `n` (and optionally number of rows `nr` or
                number of columns `nc`), or you can set the layout directly with `layout`.
              """,
              [
                :(plot(randn(100,5), layout=@layout([a;b;grid(1,3)]), t=[:line :histogram :scatter :steppre :bar], nbins=10, leg=false))
              ]),
  PlotExample("Adding to subplots",
              "Note here the automatic grid layout, as well as the order in which new series are added to the plots.",
              [
                :(plot(Plots.fakedata(100,10), layout=4, palette=[:grays :blues :heat :lightrainbow], bg_inside=[:orange :pink :darkblue :black]))
              ]),
  PlotExample("",
              "",
              [
                :(plot!(Plots.fakedata(100,10)))
              ]),
  PlotExample("Open/High/Low/Close",
              "Create an OHLC chart.  Pass in a list of (open,high,low,close) tuples as your `y` argument.  This uses recipes to first convert the tuples to OHLC objects, and subsequently create a :path series with the appropriate line segments.",
              [
                :(n=20),
                :(hgt=rand(n)+1),
                :(bot=randn(n)),
                :(openpct=rand(n)),
                :(closepct=rand(n)),
                :(y = OHLC[(openpct[i]*hgt[i]+bot[i], bot[i]+hgt[i], bot[i], closepct[i]*hgt[i]+bot[i]) for i in 1:n]),
                :(ohlc(y))
              ]),
  PlotExample("Annotations",
              "The `annotations` keyword is used for text annotations in data-coordinates.  Pass in a tuple (x,y,text) or a vector of annotations.  `annotate!(ann)` is shorthand for `plot!(; annotation=ann)`.  Series annotations are used for annotating individual data points.  They require only the annotation... x/y values are computed.  A `PlotText` object can be build with the method `text(string, attr...)`, which wraps font and color attributes.",
              [
                :(y = rand(10)),
                :(plot(y, annotations = (3,y[3],text("this is #3",:left)), leg=false)),
                :(annotate!([(5, y[5], text("this is #5",16,:red,:center)),
                             (10, y[10], text("this is #10",:right,20,"courier"))])),
                :(scatter!(linspace(2,8,6), rand(6), marker=(50,0.2,:orange), series_annotations = ["series","annotations","map","to","series",text("data",:green)]))
              ]),
  PlotExample("Custom Markers",
              "A `Plots.Shape` is a light wrapper around vertices of a polygon.  For supported backends, pass arbitrary polygons as the marker shapes.  Note: The center is (0,0) and the size is expected to be rougly the area of the unit circle.",
              [
                :(verts = [(-1.0,1.0),(-1.28,0.6),(-0.2,-1.4),(0.2,-1.4),(1.28,0.6),(1.0,1.0),
                           (-1.0,1.0),(-0.2,-0.6),(0.0,-0.2),(-0.4,0.6),(1.28,0.6),(0.2,-1.4),
                           (-0.2,-1.4),(0.6,0.2),(-0.2,0.2),(0.0,-0.2),(0.2,0.2),(-0.2,-0.6)])
                :(plot(0.1:0.2:0.9, 0.7rand(5)+0.15,
                       l=(3,:dash,:lightblue),
                       m=(Shape(verts),30,RGBA(0,0,0,0.2)),
                       bg=:pink, fg=:darkblue,
                       xlim = (0,1), ylim=(0,1), leg=false))
              ]),

  PlotExample("Contours",
              "",
              [
                :(x = 1:0.3:20),
                :(y = x),
                :(f(x,y) = sin(x)+cos(y)),
                :(contour(x, y, f, fill=true))
              ]),

  PlotExample("Pie",
              "",
              [
                :(x = ["Nerds", "Hackers", "Scientists"]),
                :(y = [0.4, 0.35, 0.25]),
                :(pie(x, y, title="The Julia Community", l=0.5))
              ]),

  PlotExample("3D",
              "",
              [
                :(n = 100),
                :(ts = linspace(0,8π,n)),
                :(x = ts .* map(cos,ts)),
                :(y = 0.1ts .* map(sin,ts)),
                :(z = 1:n),
                :(plot(x, y, z, zcolor=reverse(z), m=(10,0.8,:blues,stroke(0)), leg=false, w=5)),
                :(plot!(zeros(n),zeros(n),1:n, w=10))
              ]),

  PlotExample("DataFrames",
              "Plot using DataFrame column symbols.",
              [
                # :(import DataFrames, RDatasets),
                :(iris = RDatasets.dataset("datasets", "iris")),
                :(scatter(iris, :SepalLength, :SepalWidth, group=:Species,
                          title = "My awesome plot", xlabel = "Length", ylabel = "Width",
                          m=(0.5, [:+ :h :star7], 12), bg=RGB(.2,.2,.2)))
              ]),

PlotExample("Groups and Subplots",
            "",
            [
              :(group = rand(map(i->"group $i",1:4),100)),
              :(plot(rand(100), layout=@layout([a b;c]), group=group, n=3, linetype=[:bar :scatter :steppre]))
            ]),

PlotExample("Polar Plots",
            "",
            [
                :(Θ = linspace(0,1.5π,100)),
                :(r = abs(0.1randn(100)+sin(3Θ))),
                :(plot(Θ, r, proj=:polar, m=2))
            ]),
PlotExample("Heatmap, categorical axes, and aspect_ratio",
            "",
            [
                :(xs = [string("x",i) for i=1:10]),
                :(ys = [string("y",i) for i=1:4]),
                :(z = float((1:4)*(1:10)')),
                :(heatmap(xs, ys, z, aspect_ratio=1))
            ]),
PlotExample("Layouts, margins, label rotation, title location",
            "",
            [
                :(plot(rand(100,6),layout=@layout([a b; c]),title=["A" "B" "C"],
                title_location=:left, left_margin=[20mm 0mm],
                bottom_margin=50px, xrotation=60))
            ]),
PlotExample("Boxplot and Violin series recipes", "", [:(begin
    using RDatasets
    singers = dataset("lattice", "singer")
    boxplot(singers, :VoicePart, :Height, marker = (0.3, :orange, stroke(2)))
    violin!(singers, :VoicePart, :Height, marker = (0.2, :blue, stroke(0)))
end)])

]

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

function createStringOfMarkDownCodeValues(arr, prefix = "")
  string("`", prefix, join(sort(map(string, arr)), "`, `$prefix"), "`")
end
createStringOfMarkDownSymbols(arr) = isempty(arr) ? "" : createStringOfMarkDownCodeValues(arr, ":")

function generate_markdown(pkgname::Symbol; skip = [])

  # set up the backend, and don't show the plots by default
  pkg = backend(pkgname)
  # default(:show, false)
  default(reuse = true)

  # mkdir if necessary
  pkgdir = joinpath(IMGDIR, string(pkgname))
  try
    mkdir(pkgdir)
  catch
  end

  # open the markdown file
  md = open("$DOCDIR/$(pkgname).md", "w")

  # write(md, "## Examples for backend: $pkgname\n\n")

  write(md, "### Initialize\n\n```julia\nusing Plots\n$(pkgname)()\n```\n\n")

  for (i,example) in enumerate(_examples)

    i in skip && continue

    try

      # we want to always produce consistent results
      srand(1234)

      # run the code
      map(eval, example.exprs)

      # # save the png
      # imgname = "$(pkgname)_example_$i.png"

      # NOTE: uncomment this to overwrite the images as well
      if i == 2
        imgname = "$(pkgname)_example_$i.gif"
        gif(anim, "$pkgdir/$imgname", fps=15)
      else
        imgname = "$(pkgname)_example_$i.png"
        png("$pkgdir/$imgname")
      end

      # write out the header, description, code block, and image link
      write(md, "### $(example.header)\n\n")
      write(md, "$(example.desc)\n\n")
      write(md, "```julia\n$(join(map(string, example.exprs), "\n"))\n```\n\n")
      write(md, "![](img/$pkgname/$imgname)\n\n")

    catch ex
      # TODO: put error info into markdown?
      warn("Example $pkgname:$i failed with: $ex")
    end

    #
  end

  write(md, "- Supported arguments: $(createStringOfMarkDownCodeValues(supportedArgs(pkg)))\n")
  write(md, "- Supported values for axis: $(createStringOfMarkDownSymbols(supportedAxes(pkg)))\n")
  write(md, "- Supported values for linetype: $(createStringOfMarkDownSymbols(supportedTypes(pkg)))\n")
  write(md, "- Supported values for linestyle: $(createStringOfMarkDownSymbols(supportedStyles(pkg)))\n")
  write(md, "- Supported values for marker: $(createStringOfMarkDownSymbols(supportedMarkers(pkg)))\n")
  write(md, "- Is `subplot`/`subplot!` supported? $(subplotSupported(pkg) ? "Yes" : "No")\n\n")

  write(md, "(Automatically generated: $(now()))")
  close(md)

end

# make and display one plot
function test_examples(pkgname::Symbol, idx::Int; debug = false, disp = true)
  Plots._debugMode.on = debug
  println("Testing plot: $pkgname:$idx:$(_examples[idx].header)")
  backend(pkgname)
  backend()
  map(eval, _examples[idx].exprs)
  plt = current()
  if disp
    gui(plt)
  end
  plt
end

# generate all plots and create a dict mapping idx --> plt
function test_examples(pkgname::Symbol; debug = false, disp = true)
  Plots._debugMode.on = debug
  plts = Dict()
  for i in 1:length(_examples)

    try
      plt = test_examples(pkgname, i, debug=debug, disp=disp)
      plts[i] = plt
    catch ex
      # TODO: put error info into markdown?
      warn("Example $pkgname:$i:$(_examples[i].header) failed with: $ex")
    end
  end
  plts
end


text:
function read_dataframe(tablegroup::HDF5.Group; separate_index=true, kwargs...)
    columns = read_attribute(tablegroup, "column-order")

    if separate_index && haskey(attributes(tablegroup), "_index")
        indexdsetname = read_attribute(tablegroup, "_index")
        rownames = read(tablegroup[indexdsetname])
    else
        rownames = nothing
    end

    df = DataFrame()

    for col in columns
        column = read_matrix(tablegroup[col])
        if sum(size(column) .> 1) > 1
            @warn "column $col has more than 1 dimension for data frame $(HDF5.name(tablegroup)), skipping"
        end
        df[!, col] = column
    end

    return df, rownames
end

function read_matrix(f::HDF5.Dataset; kwargs...)
    mat = read(f)
    if HDF5.API.h5t_get_class(datatype(f)) == HDF5.API.H5T_COMPOUND
        return StructArray(mat)
    end

    if ndims(f) == 0
        return mat
    end
    if ndims(f) > 1
        mat = PermutedDimsArray(mat, ndims(mat):-1:1) # transpose for h5py compatibility
    end
    if haskey(attributes(f), "categories")
        categories = f[read_attribute(f, "categories")]
        ordered =
            haskey(attributes(categories), "ordered") &&
            read_attribute(categories, "ordered") == true
        cats = read(categories)
        mat = mat .+ 0x1
        mat = compress(
            CategoricalArray{eltype(cats), ndims(mat)}(
                mat,
                CategoricalPool{eltype(cats), eltype(mat)}(cats, ordered),
            ),
        )
    end
    return mat
end

function read_matrix(f::HDF5.Group; kwargs...)
    enctype = read_attribute(f, "encoding-type")

    if enctype == "csc_matrix" || enctype == "csr_matrix"
        shape = read_attribute(f, "shape")
        iscsr = enctype[1:3] == "csr"

        indptr = read(f, "indptr")
        indices = read(f, "indices")
        data = read(f, "data")

        indptr .+= eltype(indptr)(1)
        indices .+= eltype(indptr)(1)

        # the row indices in every column need to be sorted
        @views for (colstart, colend) in zip(indptr[1:(end - 1)], indptr[2:end])
            sort!(indices[colstart:(colend - 1)], data[colstart:(colend - 1)])
        end

        if iscsr
            reverse!(shape)
        end
        mat = SparseMatrixCSC(shape..., indptr, indices, data)
        return iscsr ? mat' : mat
    else
        error("unknown encoding $enctype")
    end
end

function read_dict_of_matrices(f::HDF5.Group; kwargs...)
    return Dict{String, AbstractArray{<:Number}}(
        key => read_matrix(f[key]; kwargs...) for key in keys(f)
    )
end

read_auto(f::HDF5.Dataset; kwargs...) = (read_matrix(f; kwargs...), nothing)
function read_auto(f::HDF5.Group; kwargs...)
    if haskey(attributes(f), "encoding-type")
        enctype = read_attribute(f, "encoding-type")
        if enctype == "dataframe"
            return read_dataframe(f; kwargs...)
        elseif endswith(enctype, "matrix")
            return read_matrix(f; kwargs), nothing
        else
            error("unknown encoding $enctype")
        end
    else
        return read_dict_of_mixed(f; kwargs...), nothing
    end
end

function read_dict_of_mixed(f::HDF5.Group; kwargs...)
    ret = Dict{
        String,
        Union{
            DataFrame,
            StructArray,
            <:AbstractArray{<:Number},
            <:AbstractArray{<:AbstractString},
            <:AbstractString,
            <:Number,
            Dict,
        },
    }()
    for k in keys(f)
        ret[k] = read_auto(f[k]; kwargs...)[1] # assume data frames are properly aligned, so we can discard rownames
    end
    return ret
end

function write_attr(parent::Union{HDF5.File, HDF5.Group}, name::AbstractString, data; kwargs...)
    if haskey(parent, name)
        delete_object(parent, name)
    end
    write_impl(parent, name, data; kwargs...)
end

function write_impl(
    parent::Union{HDF5.File, HDF5.Group},
    name::AbstractString,
    data::Union{<:Number, <:AbstractString};
    kwargs...,
)
    parent[name] = data
end

function write_impl(
    parent::Union{HDF5.File, HDF5.Group},
    name::AbstractString,
    data::AbstractDict{<:AbstractString, <:Any};
    kwargs...,
)
    if length(data) > 0
        g = create_group(parent, name)
        for (key, val) in data
            write_impl(g, key, val; kwargs...)
        end
    end
end

function write_impl(
    parent::Union{HDF5.File, HDF5.Group},
    name::AbstractString,
    data::CategoricalArray;
    kwargs...,
)
    write_impl(parent, name, data.refs .- 0x1; kwargs...)
    write_impl(parent, "__categories/$name", levels(data); kwargs...)
    attributes(parent["__categories/$name"])["ordered"] = UInt8(isordered(data))
    attributes(parent[name])["categories"] = HDF5.Reference(parent["__categories"], name)
end

function write_impl(
    parent::Union{HDF5.File, HDF5.Group},
    name::AbstractString,
    data::AbstractDataFrame;
    index::AbstractVector{<:AbstractString}=nothing,
    kwargs...,
)
    g = create_group(parent, name)
    attrs = attributes(g)
    attrs["encoding-type"] = "dataframe"
    attrs["encoding-version"] = "0.1.0"
    attrs["column-order"] = names(data)

    for (name, column) in pairs(eachcol(data))
        write_impl(g, string(name), column; kwargs...)
    end

    idxname = "_index"
    columns = names(data)
    if !isnothing(index)
        while idxname ∈ columns
            idxname = "_" * idxname
        end
    else
        if idxname ∈ columns
            index = data[!, idxname]
            select!(data, Not(idxname))
        else
            @warn "Data frame $(HDF5.name(parent))/$name does not have an _index column, a row number index will be written"
            index = collect(1:nrow(data))
        end
    end
    g = parent[name]
    write_impl(g, idxname, values(index); kwargs...)
    attributes(g)["_index"] = idxname
end

# see https://github.com/JuliaIO/HDF5.jl/issues/827
# and https://github.com/JuliaIO/HDF5.jl/issues/826
write_impl(
    parent::Union{HDF5.File, HDF5.Group},
    name::AbstractString,
    data::Union{<:AbstractArray{Bool}, BitArray{1}};
    extensible::Bool=false,
    compress::UInt8=UInt8(9),
    kwargs...,
) = write_impl(parent, name, Int8.(data); extensible=extensible, compress=compress, kwargs...)

write_impl(parent::Union{HDF5.File, HDF5.Group}, name::AbstractString, data::SubArray; kwargs...) =
    write_impl(parent, name, copy(data); kwargs...)

function write_impl(
    parentgrp::Union{HDF5.File, HDF5.Group},
    name::AbstractString,
    data::AbstractArray,
    ;
    extensible::Bool=false,
    compress::UInt8=0x9,
    kwargs...,
)
    if ndims(data) > 1
        data =
            data isa PermutedDimsArray && typeof(data).parameters[3] == Tuple(ndims(data):-1:1) ?
            parent(data) : permutedims(data, ndims(data):-1:1) # transpose for h5py compatibility
    end                                             # copy because HDF5 apparently can't handle lazy Adjoints

    if extensible
        dims = (size(data), Tuple(-1 for _ in 1:ndims(data)))
    else
        dims = size(data)
    end
    dtype = datatype(data)
    write_impl_array(parentgrp, name, data, dtype, dims, compress)
end

function write_impl_array(
    parent::Union{HDF5.File, HDF5.Group},
    name::AbstractString,
    data::AbstractArray{<:Number},
    dtype::HDF5.Datatype,
    dims::Union{Tuple{Vararg{<:Integer}}, Tuple{Tuple{Vararg{<:Integer, N}}, Tuple{Vararg{<:Integer, N}}}},
    compress::UInt8,
) where N
    chunksize = HDF5.heuristic_chunk(data)
    if length(chunksize) == 0
        chunksize = Tuple(100 for _ in 1:ndims(data))
    end
    d = create_dataset(parent, name, dtype, dims, chunk=chunksize, deflate=compress)
    write_dataset(d, dtype, data)
end

# variable-length strings cannot be compressed in HDF5
function write_impl_array(
    parent::Union{HDF5.File, HDF5.Group},
    name::AbstractString,
    data::AbstractArray{<:AbstractString},
    dtype::HDF5.Datatype,
    dims::Union{Tuple{Vararg{<:Integer}}, Tuple{Tuple{Vararg{<:Integer, N}}, Tuple{Vararg{<:Integer, N}}}},
    compress::UInt8,
) where N
    d = create_dataset(parent, name, dtype, dims)
    write_dataset(d, dtype, data)
end

function write_impl(
    parent::Union{HDF5.File, HDF5.Group},
    name::AbstractString,
    data::SparseMatrixCSC{<:Number, <:Integer};
    transposed=false,
    kwargs...,
)
    g = create_group(parent, name)
    attrs = attributes(g)
    attrs["encoding-type"] = transposed ? "csr_matrix" : "csc_matrix"
    attrs["encoding-version"] = "0.1.0"

    shape = collect(size(data))
    transposed && reverse!(shape)
    attrs["shape"] = shape
    write_impl(g, "indptr", data.colptr .- 1, extensible=true)
    write_impl(g, "indices", data.rowval .- 1, extensible=true)
    write_impl(g, "data", data.nzval, extensible=true)
end

write_impl(
    prt::Union{HDF5.File, HDF5.Group},
    name::AbstractString,
    data::Adjoint{T, SparseMatrixCSC{T, V}} where {T <: Number, V <: Integer};
    kwargs...,
) = write_impl(prt, name, parent(data), transposed=true; kwargs...)

write_impl(
    parent::Union{HDF5.File, HDF5.Group},
    name::AbstractString,
    data::Union{HDF5.Dataset, SparseDataset, TransposedDataset};
    kwargs...,
) = copy_object(data, parent, name)

write_impl(parent::Union{HDF5.File, HDF5.Group}, name::AbstractString, ::Nothing; kwargs...) =
    nothing

function write_impl(
    parent::Union{HDF5.File, HDF5.Group},
    name::AbstractString,
    data::StructArray;
    extensible::Bool=false,
    compress::UInt8=0x9,
    kwargs...)
    ety = eltype(data)
    dtype = create_datatype(HDF5.API.H5T_COMPOUND, sizeof(ety))
    for (i, (fname, ftype)) in enumerate(zip(fieldnames(ety), fieldtypes(ety)))
        HDF5.API.h5t_insert(dtype, string(fname), fieldoffset(ety, i), _datatype(ftype))
    end
    chunksize = HDF5.heuristic_chunk(data)
    dims = size(data)
    if extensible
        dims = (size(data), Tuple(-1 for _ in 1:ndims(data)))
    end
    dset = create_dataset(parent, name, dtype, dataspace(dims), chunk=chunksize, deflate=compress, kwargs...)
    write_dataset(dset, dtype, [val for row in data for val in row])
end

_datatype(::Type{T}) where T = datatype(T)

function _datatype(::Type{T}) where {T <: AbstractString}
    strdtype = HDF5.API.h5t_copy(HDF5.API.H5T_C_S1)
    HDF5.API.h5t_set_size(strdtype, HDF5.API.H5T_VARIABLE)
    HDF5.API.h5t_set_cset(strdtype, HDF5.API.H5T_CSET_UTF8)

    HDF5.Datatype(strdtype)
end

ext: jl


text:
using TypedDelegation
using Test

import Base: string, abs,
             (+), (-), (*), div, fld, cld,
             (==), (!=), (<), (<=), (>=), (>)

struct MyInt16
    value::Int16
end

@delegate_onefield( MyInt16, value, [ string ]);
@delegate_onefield_twovars( MyInt16, value, [ (==), (!=), (<), (<=), (>=), (>) ] );
@delegate_onefield_astype( MyInt16, value, [ abs, (-) ]); # unary minus
@delegate_onefield_twovars_astype( MyInt16, value, [ (+), (-), (*), div, fld, cld ] );

# values to be used in testing
three = MyInt16(3);
seven = MyInt16(7);

# testing four of the eight exported macros, one of each kind

  @test  string(three) == "3"
  @test  three < seven
  @test -three < three
  @test  seven == div( seven*three, three )

# end tests


text: abstract type Motion end


text:
"""
    block_jackknife(Y::JArray{Float64,2}, subsample::Float64)

Generate block jackknife samples as in Kunsch (1989).

This technique subsamples a time series dataset by removing, in turn, all the blocks of consecutive observations
with a given size.

# Arguments
- `Y`: observed measurements (`nxT`), where `n` and `T` are the number of series and observations.
- `subsample`: block size as a percentage of the original sample size. It is bounded between 0 and 1.

# References
Kunsch (1989) and Pellegrino (2020)
"""
function block_jackknife(Y::JArray{Float64,2}, subsample::Float64)

    # Check inputs
    check_bounds(subsample, 0, 1);

    # Dimensions
    n, T = size(Y);

    # Block size
    block_size = Int64(ceil(subsample*T));
    if block_size == 0
        error("subsample is too small!");
    end

    # Number of block jackknifes samples - as in Kunsch (1989)
    samples = T-block_size+1;

    # Initialise jackknife_data
    jackknife_data = JArray{Float64,3}(undef, n, T, samples);

    # Loop over j=1, ..., samples
    for j=1:samples

        # Index of missings
        ind_j = collect(j:j+block_size-1);

        # Block jackknife data
        jackknife_data[:, :, j] = Y;
        jackknife_data[:, ind_j, j] .= missing;
    end

    # Return jackknife_data
    return jackknife_data;
end

"""
    optimal_d(n::Int64, T::Int64)

Select the optimal value for d.

# Arguments
- `n`: Number of series
- `T`: Number of observations
"""
function optimal_d(n::Int64, T::Int64)
    objfun_array = zeros(n*T);

    for d=1:fld(n*T,2)
        objfun_array[d] = objfun_optimal_d(n, T, d);
    end

    return argmax(objfun_array);
end

"""
    objfun_optimal_d(n::Int64, T::Int64, d::Int64)

Objective function to select the optimal value for d.

# Arguments
- `n`: Number of series
- `T`: Number of observations
- `d`: Candidate d
"""
function objfun_optimal_d(n::Int64, T::Int64, d::Int64)
    fun = no_combinations(n*T, d) - (d>=n).*no_combinations(n*T-n, d-n)*T;

    for i=2:fld(d, n)
        fun -= (-1^(i-1))*no_combinations(T, i)*no_combinations(n*T-i*n, d-i*n);
    end

    return fun;
end

"""
    artificial_jackknife(Y::JArray{Float64,2}, subsample::Float64, max_samples::Int64)

Generate artificial jackknife samples as in Pellegrino (2020).

The artificial delete-d jackknife is an extension of the delete-d jackknife for dependent data problems.
This technique replaces the actual data removal step with a fictitious deletion, which consists of
imposing `d`-dimensional (artificial) patterns of missing observations to the data. This approach
does not alter the data order nor destroy the correlation structure.

# Arguments
- `Y`: observed measurements (`nxT`), where `n` and `T` are the number of series and observations.
- `subsample`: `d` as a percentage of the original sample size. It is bounded between 0 and 1.
- `max_samples`: if `C(n*T,d)` is large, artificial_jackknife would generate `max_samples` jackknife samples.

# References
Pellegrino (2020)
"""
function artificial_jackknife(Y::JArray{Float64,2}, subsample::Float64, max_samples::Int64)

    # Dimensions
    n, T = size(Y);
    nT = n*T;

    # Check inputs
    check_bounds(subsample, 0, 1);

    # Set d using subsample
    d = Int64(ceil(subsample*nT));

    # Error management
    if d == 0
        error("subsample is too small!");
    end

    # Warning
    if subsample > 0.5
        @warn "this algorithm might be unstable for `subsample` larger than 0.5!";
    end

    # Get vec(Y)
    vec_Y = convert(JArray{Float64}, Y[:]);

    # Initialise loop (controls)
    C_nT_d = no_combinations(nT, d);
    samples = convert(Int64, min(C_nT_d, max_samples));
    zeros_vec = zeros(d);

    # Initialise loop (output)
    ind_missings = Array{Int64}(zeros(d, samples));
    jackknife_data = JArray{Float64,3}(undef, n, T, samples);

    # Loop over j=1, ..., samples
    for j=1:samples

        # First draw
        if j == 1
            if samples == C_nT_d
                ind_missings[:,j] = rand_without_replacement(nT, d);
            else
                ind_missings[:,j] = rand_without_replacement(n, T, d);
            end

        # Iterates until ind_missings[:,j] is neither a vector of zeros, nor already included in ind_missings
        else
            while ind_missings[:,j] == zeros_vec || is_vector_in_matrix(ind_missings[:,j], ind_missings[:, 1:j-1])
                if samples == C_nT_d
                    ind_missings[:,j] = rand_without_replacement(nT, d);
                else
                    ind_missings[:,j] = rand_without_replacement(n, T, d);
                end
            end
        end

        # Add (artificial) missing observations
        vec_Y_j = copy(vec_Y);
        vec_Y_j[ind_missings[:,j]] .= missing;

        # Store data
        jackknife_data[:, :, j] = reshape(vec_Y_j, n, T);
    end

    # Return jackknife_data
    return jackknife_data;
end


text:
module ScalarWaveEquation

export FunArr
export FunArr2D
export scalarwave_and_energy
export integrate_difference

mutable struct FunArr{T, U}
	xmin::T            # x-range minimum
	xmax::T            # x-range maximum
	dx::T              # step size of z-data
	ydat::Vector{U}    # range data

	function FunArr{T, U}(ydat_) where {T, U}
		size = length(ydat_)
		if size < 3
			error("The y-data Vector must have a length 3 or greater.")
		end

		xmin_ = 0.0
		xmax_ = 1.0
		dx_   = (xmax_ - xmin_)/(size - 1)
		new(xmin_, xmax_, dx_, ydat_)
	end
end

mutable struct FunArr2D{T, U, V}
	xmin::T
	xmax::T
	dx::T
	Nx::Int
	tmin::U
	tmax::U
	dt::U
	Nt::Int
	xdat::Array{T}
	tdat::Array{U}
	u2dat::Array{V, 2}

	function FunArr2D{T, U, V}(u2dat_, t_duration) where {T, U, V}
		# test for appropriate sizes
		Nx_ = size(u2dat_, 1)
		if Nx_ < 3
			error("The first dimension of the array must have 3 or more elements.")
		end
		Nt_ = size(u2dat_, 2)
		if Nt_ < 3
			error("The second dimension of the array must have 3 or more elements.")
		end

		# the duration of elapsed time must be positive
		if t_duration <= 0.0
			error("The duration of time that passes must be positive.")
		end

		xmin_ = 0.0
		xmax_ = 1.0
		tmin_ = 0.0
		tmax_ = t_duration
		dx_ = (xmax_ - xmin_)/(Nx_ - 1)
		dt_ = (tmax_ - tmin_)/(Nt_ - 1)
		xdat_ = range(0.0, stop = 1.0, length = Nx_)
		tdat_ = range(0.0, stop = t_duration, length = Nt_)
		new(xmin_, xmax_, dx_, Nx_, tmin_, tmax_, dt_, Nt_, xdat_, tdat_, u2dat_)
	end
end

function first_deriv(f::FunArr{T, U})::FunArr{T, U} where {T, U}
	N  = length(f.ydat)
	dx = U(f.dx)
	dydat = similar(f.ydat)

	# take second derivative of end-points
	dydat[1] = (f.ydat[2] - f.ydat[1]  ) / dx
	dydat[N] = (f.ydat[N] - f.ydat[N-1]) / dx

	# take the second derivative of all points in between
	for i in 2:(N-1)
		dydat[i] = (f.ydat[i+1] - f.ydat[i-1]) / (2*dx)
	end

	FunArr{T, U}(dydat)
end

function second_deriv(f::FunArr{T, U})::FunArr{T, U} where {T, U}
	N  = length(f.ydat)
	dx = U(f.dx)
	d2ydat = similar(f.ydat)

	# take second derivative of end-points
	d2ydat[1] = (f.ydat[3] - 2.0*f.ydat[2] + f.ydat[1]) / dx^2
	d2ydat[N] = (f.ydat[N] - 2.0*f.ydat[N-1] + f.ydat[N-2]) / dx^2

	# take the second derivative of all points in between
	for i in 2:(N-1)
		d2ydat[i] = (f.ydat[i+1] - 2*f.ydat[i] + f.ydat[i-1]) / dx^2
	end

	FunArr{T, U}(d2ydat)
end

function forward_euler2_disc(i0::Int, yfun::FunArr{T, U}, vfun::FunArr{T, U}, afun::FunArr{T, U}, h::T) where {T, U}
	# i    : index value between 1 and N inclusive
	# yfun : the discretized position function
	# vfun : the discretized velocity function
	# afun : the discretized acceleration function
	y0 = yfun.ydat[i0]
	v0 = vfun.ydat[i0]
	a0 = afun.ydat[i0]
	dt = U(h)
	v1 = v0 + a0*dt
	y1 = y0 + v0*dt + 0.5*a0*dt*dt

	return y1, v1
end

function get_energy(utfun::FunArr{T, U}, uxfun::FunArr{T, U}) where {T, U}
	etotal = 0.0
	N = length(utfun.ydat)
	h = utfun.dx

	for i in 1:N
		etotal += utfun.ydat[i]^2 + uxfun.ydat[i]^2
	end

	etotal *= 0.5*h

	return etotal
end

function nearest_lower_index(N::I, xmin::T, xmax::T, x::T) where {T, I<:Int}
	# takes care of the rand() = 0.0 corner case
	if x == 0.0
		x += eps(1.0)
	end
	@assert xmax > x > xmin
	ratio = (x - xmin)/(xmax - xmin)
	lower_index = ceil(typeof(N), (N-1)*ratio)
	return lower_index
end

function bilinear_interp_regular(u2fun::FunArr2D{T, U, V}, x::T, t::U) where {T, U, V}
	# performs bilinear interpolation on a function discretized on a regular grid
	# find the x's and y's around the given (x, y) point
	# they will be (x1, y1), (x1, y2), (x2, y1), (x2, y2)

	xsize = u2fun.Nx; xmin = u2fun.xmin; xmax = u2fun.xmax; dx = u2fun.dx
	tsize = u2fun.Nt; tmin = u2fun.tmin; tmax = u2fun.tmax; dt = u2fun.dt

	Nx1 = nearest_lower_index(xsize, xmin, xmax, x)
	Nt1 = nearest_lower_index(tsize, tmin, tmax, t)

	x1 = xmin + (Nx1-1)*dx
	x2 = xmin + Nx1*dx
	t1 = tmin + (Nt1-1)*dt
	t2 = tmin + Nt1*dt

	u11 = u2fun.u2dat[Nx1  , Nt1  ]
	u12 = u2fun.u2dat[Nx1  , Nt1+1]
	u21 = u2fun.u2dat[Nx1+1, Nt1  ]
	u22 = u2fun.u2dat[Nx1+1, Nt1+1]

	coeff = 1.0/((x2 - x1)*(t2 - t1))
	xV    = [x2 - x, x - x1]
	tV    = [t2 - t, t - t1]
	uM    = [u11 u12; u21 u22]

	interp_value = (xV'*uM*tV)*coeff
	return interp_value
end

function integrate_difference(u2fun1::FunArr2D{T, U, V}, u2fun2::FunArr2D{T, U, V}, Nxdiv::I, Ntdiv::I) where {T, U, V, I<:Int}
	@assert u2fun1.xmin == u2fun2.xmin
	@assert u2fun1.xmax == u2fun2.xmax
	@assert u2fun1.tmin == u2fun2.tmin
	@assert u2fun1.tmax == u2fun2.tmax

	xmin = u2fun1.xmin; xmax = u2fun1.xmax
	tmin = u2fun1.tmin; tmax = u2fun1.tmax

	dx = (xmax - xmin)/T(Nxdiv - 1)
	dt = (tmax - tmin)/U(Ntdiv - 1)

	total = 0.0
	for ix in 1:(Nxdiv-1)
		x = (ix - 0.5)*dx
		for it in 1:(Ntdiv-1)
			t = (it - 0.5)*dt
			term_fun1 = bilinear_interp_regular(u2fun1, x, t)
			term_fun2 = bilinear_interp_regular(u2fun2, x, t)
			total += dx * dt * (term_fun1 - term_fun2)^2
		end
	end

	return total
end

function scalarwave_and_energy(ufun::FunArr{T, U}, utfun::FunArr{T, U}, Nt::I, duration::T) where {T, U, I<:Int}
	# ufun  : initial positions at t = 0
	# utfun : initial first time derivatives at t = 0
	# Nt    : total number of time steps
	# duration : total amount of time elapsed
	# returns : FunArr2D instance for the u(x, t) grid, and an array of energies

	Nx    = length(ufun.ydat)
	Nx_ut = length(utfun.ydat)
	if (Nx != Nx_ut)
		error("The initial positions and initial velocities Arrays must be of the same size.")
	end

	h = duration/(Nt-1)  # time step
	u2fun    = FunArr2D{Float64, Float64, Float64}(zeros(Nx, Nt), duration)
	energies = zeros(U, Nt)

	for it in 1:Nt
		# propagate forward in time
		uxfun = first_deriv(ufun)
		afun  = second_deriv(ufun)
		for ix in 1:Nx
			u_1, ut_1 = forward_euler2_disc(ix, ufun, utfun, afun, h)
			ufun.ydat[ix] = u_1
			utfun.ydat[ix] = ut_1
		end

		# set the values in the 2D matrix
		for ix in 1:Nx
			u2fun.u2dat[ix, it] = ufun.ydat[ix]
		end

		# collect the energy
		energies[it] = get_energy(utfun, uxfun)
	end

	return u2fun, energies
end

end # module


text:
function sequencedist(x::Seq,y::Seq;sub::Int=1,ins::Int=10,del::Int=10)
	levenshtein(x[2],y[2],del,ins,sub) # high penalty for indels
end

function sequencedist(x::Reference,y::Reference;sub::Int=1,ins::Int=10,del::Int=10)
	@assert all(v->isequal(v[1][1],v[2][1]), zip(x, y)) "Cannot compare reference sequences with different names"

	# sum distance
	sum(v->sequencedist(v[1],v[2],sub=sub,ins=ins,del=del), zip(x, y))
end

# elementwise comparison of strings
function hamming(x::Seq,y::Seq)
	length(x[2])!=length(y[2]) && return typemax(Int)
	s = 0
	for (cx,cy) in zip(x[2],y[2])
		s += cx!=cy
	end
	s
end

function hamming(x::Reference,y::Reference) # not used atm.
	sum(v->hamming(v[1],v[2]), zip(x,y))
end

function hasindels(x::Seq,y::Seq)
	length(x[2]) != length(y[2]) && return true # different length, must have indels
	sequencedist(x,y) != hamming(x,y) # if levenshtein distance is different from hamming distance, there are indels
end

function hasindels(x::Reference,y::Reference)
	any(v->hasindels(v[1],v[2]), zip(x,y))
end

function reference_sanity_check(sample::Sample,
                                references::Array{Reference},
                                referenceNames::Array{String,1}, log)
	if isempty(sample.consensus)
		printiferror(log, "$(sample.name): No consensus file specified. (Alignment failed?)")
		return
	end

	c = loadfasta(sample.consensus)

	d = zeros(Int,length(references))
	for (i,r) in enumerate(references)
		# allow some name mismatching since we are normally including PhiX
		if length(c) == 1 && length(r) == 1
			d[i] = sequencedist(c[1],r[1]) # allow name mismatch
		elseif length(c) != length(r)
			d[i] = typemax(Int)
		else
			d[i] = sequencedist(c,r)
		end
	end

	i = indmin(d)
	j = findfirst(referenceNames,sample.referenceName)

	if d[j] > d[i]
		ri = splitdir(referenceNames[i])[2]
		rj = splitdir(referenceNames[j])[2]
		printifwarning(log, "$(sample.name) is closer to reference \"$ri\" than to reference \"$rj\".")
	end

	# are there any indels (when compared to the best reference)?
	hasindels(c, references[i]) && printifwarning(log, "$(sample.name) has indels.")
	return
end

function reference_sanity_check(samples::Vector{Sample},
                                refs::Vector;
                                log=devnull)
	r = Reference[loadfasta(f[4]) for f in refs] # loaded fastas
	names = String[f[2] for f in refs]

	for i=1:length(samples)
		reference_sanity_check(samples[i], r, names, log)
	end
end


text:
using Documenter
using PersistenceDiagrams
using Test

if VERSION ≥ v"1.7-DEV" || VERSION < v"1.6-DEV"
    @warn "Doctests were set up on Julia v1.6. Skipping."
else
    DocMeta.setdocmeta!(
        PersistenceDiagrams, :DocTestSetup, :(using PersistenceDiagrams); recursive=true
    )
    doctest(PersistenceDiagrams)
end


text:
using BioMASS

const model = Model("erbb_network_jl")

if abspath(PROGRAM_FILE) == @__FILE__
    optimize(
        model,
        parse(Int64, ARGS[1]),
        max_generation=50,
        allowable_error=6.0,
        popsize=3,
        local_search_method="DE",
        maxiter=200,
    )
end


text:
#=== CW cell ===#
struct Cell <: AbstractCell
    d::Int
    boundary::Vector{Cell}
    hash::UInt64
    Cell(d::Int, cells::Vector{Cell}) = new(d, cells, hash(cells))
end
Cell() = Cell(0, Cell[])
Cell(d::Int, splx::Cell...)= Cell(d, [splx...])

# Private methods

function Base.show(io::IO, c::Cell)
    vstr = ""
    print(io, "C$(c.d)")
    if length(c.boundary) > 0
        vstr = foldl((v,u)->v*"$(u), ", c.boundary, init="")[1:end-2]
        print(io, "[$vstr]")
    end
end
Base.eltype(c::Cell) = Cell

# Public methods

dim(c::Cell) = c.d

==(a::Cell, b::Cell) = a.hash == b.hash

hash(c::Cell) = c.hash

faces(c::Cell) = c.boundary

function vertecies(c::Cell)
    ret = Cell[]
    return ret
end

function union(u::Cell, v::Cell)
    dim(u) != dim(v) && throw(DimensionMismatch("Cell dimensions must match for union"))
    return Cell(dim(v)+1, u, v)
end


text:
@testset "Forcings" begin
    fx = zeros(5,5)
    fy = zeros(5,5)
    f1 = zeros(30)
    # For structs
    sys = Swalbe.SysConst(Lx=5, Ly=5)
    state = Swalbe.Sys(sys, "CPU")
    state.height .= 1.0
    # One dim model
    sys1D = Swalbe.SysConst_1D(L=30)
    state1D = Swalbe.Sys(sys1D)
    state1D.height .= 1.0
    @testset "Slippage" begin
        # No velocities
        Swalbe.slippage!(fx, fy, ones(5,5), zeros(5,5), zeros(5,5), 1.0, 1/6)
        Swalbe.slippage!(state, sys)
        @test all(fx .== 0.0)
        @test all(fy .== 0.0)
        @test all(state.slipx .== 0.0)
        @test all(state.slipy .== 0.0)

        # Velocity in x
        Swalbe.slippage!(fx, fy, ones(5,5), fill(0.1,5,5), zeros(5,5), 1.0, 1/6)
        state.velx .= 0.1
        Swalbe.slippage!(state, sys)
        @test all(fx .== 0.1/11)
        @test all(fy .== 0.0)
        # At some point I have to find out how to use all() with atol
        # Thanks to Tilman for the tip!
        @test all(isapprox.(state.slipx, 0.1/11; atol=1e-8))
        @test all(state.slipy .== 0.0)
        # Velocity in y
        Swalbe.slippage!(fx, fy, ones(5,5), zeros(5,5), fill(0.1,5,5), 1.0, 1/6)
        state.velx .= 0.0
        state.vely .= 0.1
        Swalbe.slippage!(state, sys)
        @test all(fx .== 0.0)
        @test all(fy .== 0.1/11)
        @test all(state.slipx .== 0.0)
        all(isapprox.(state.slipy, 0.1/11; atol=1e-8))
        # Velocity
        Swalbe.slippage!(fx, fy, ones(5,5), fill(-0.1,5,5), fill(0.1,5,5), 1.0, 1/6)
        state.velx .= -0.1
        state.vely .= 0.1
        Swalbe.slippage!(state, sys)
        @test all(fx .== -0.1/11)
        @test all(fy .== 0.1/11)
        @test all(isapprox.(state.slipx, -0.1/11; atol=1e-8))
        @test all(isapprox.(state.slipy, 0.1/11; atol=1e-8))
        # No slip
        Swalbe.slippage!(fx, fy, ones(5,5), fill(-0.1,5,5), fill(0.1,5,5), 0.0, 1/6)
        state.velx .= -0.1
        state.vely .= 0.1
        sys2 = Swalbe.SysConst(Lx=5,Ly=5,δ=0.0)
        Swalbe.slippage!(state, sys2)
        @test all(fx .== -0.1/2)
        @test all(fy .== 0.1/2)
        @test all(isapprox.(state.slipx, -0.1/2; atol=1e-8))
        @test all(isapprox.(state.slipy, 0.1/2; atol=1e-8))
    end

    @testset "Slippage 1D" begin
        # No velocities
        Swalbe.slippage!(f1, ones(30), zeros(30), 1.0, 1/6)
        Swalbe.slippage!(state1D, sys1D)
        @test all(f1 .== 0.0)
        @test all(state1D.slip .== 0.0)
        # Velocity in x
        Swalbe.slippage!(f1, ones(30), fill(0.1,30), 1.0, 1/6)
        state1D.vel .= 0.1
        state1D.slip .= 0.0
        Swalbe.slippage!(state1D, sys1D)
        @test all(f1 .== 0.1/11)
        @test all(isapprox.(state1D.slip, 0.1/11; atol=1e-8))

        # No slip
        Swalbe.slippage!(f1, ones(30), fill(-0.1,30), 0.0, 1/6)
        state1D.vel .= -0.1
        state1D.slip .= 0.0
        sys1D2 = Swalbe.SysConst_1D(L=30, δ=0.0)
        Swalbe.slippage!(state1D, sys1D2)
        @test all(f1 .== -0.1/2)
        @test all(isapprox.(state1D.slip, -0.1/2; atol=1e-8))
    end

    @testset "Pressure gradient" begin
        sys = Swalbe.SysConst(Lx=5, Ly=5)
        state = Swalbe.Sys(sys, "CPU")
        state.pressure .= reshape(collect(1.0:25),5,5)
        state.height .= 1.0
        Swalbe.h∇p!(state)

        solx = [-1.5 -1.5 -1.5 -1.5 -1.5;
                 1.0 1.0 1.0 1.0 1.0;
                 1.0 1.0 1.0 1.0 1.0;
                 1.0 1.0 1.0 1.0 1.0;
                -1.5 -1.5 -1.5 -1.5 -1.5]
        soly = [-7.5 5.0 5.0 5.0 -7.5;
                -7.5 5.0 5.0 5.0 -7.5;
                -7.5 5.0 5.0 5.0 -7.5;
                -7.5 5.0 5.0 5.0 -7.5;
                -7.5 5.0 5.0 5.0 -7.5]

        @test all(state.h∇px .== solx)
        @test all(state.h∇py .== soly)
    end

    @testset "Thermal" begin
        f1 = ones(50,50)
        f2 = ones(50,50)
        vartest = 2*0.01/11
        Swalbe.thermal!(f1, f2, ones(50,50), 0.01, 1/6, 1.0)
        @test mean(f1) ≈ 0.0 atol=1e-2
        @test mean(f2) ≈ 0.0 atol=1e-2
        @test var(f1) ≈ vartest atol=vartest/10
        @test var(f2) ≈ vartest atol=vartest/10
        # Using structs
        sys = Swalbe.SysConst(Lx=50, Ly=50, kbt=0.01)
        state = Swalbe.Sys(sys, "CPU", kind="thermal")
        Swalbe.thermal!(state, sys)
        @test mean(state.kbtx) ≈ 0.0 atol=1e-2
        @test mean(state.kbty) ≈ 0.0 atol=1e-2
        @test var(state.kbtx) ≈ vartest atol=vartest/10
        @test var(state.kbtx) ≈ vartest atol=vartest/10
        # More thermal energy
        vartest = 0.2/11
        Swalbe.thermal!(f1, f2, ones(50,50), 0.1, 1/6, 1.0)
        @test mean(f1) ≈ 0.0 atol=1e-2
        @test mean(f2) ≈ 0.0 atol=1e-2
        @test var(f1) ≈ vartest atol=vartest/10
        @test var(f2) ≈ vartest atol=vartest/10
        # Using structs
        sys = Swalbe.SysConst(Lx=50, Ly=50, kbt=0.1)
        state = Swalbe.Sys(sys, "CPU", kind="thermal")
        Swalbe.thermal!(state, sys)
        @test mean(state.kbtx) ≈ 0.0 atol=1e-2
        @test mean(state.kbty) ≈ 0.0 atol=1e-2
        @test var(state.kbtx) ≈ vartest atol=vartest/10
        @test var(state.kbtx) ≈ vartest atol=vartest/10

    end
    @testset "Thermal 1D" begin
        f1D = ones(100000)
        vartest = 2*0.01/11
        Swalbe.thermal!(f1D, ones(100000), 0.01, 1/6, 1.0)
        @test mean(f1D) ≈ 0.0 atol=1e-2
        @test var(f1D) ≈ vartest atol=vartest/10
        # Using structs
        sys = Swalbe.SysConst_1D(L=100000,δ=1.0,kbt=0.01)
        state = Swalbe.Sys(sys, kind="thermal")
        Swalbe.thermal!(state, sys)
        @test mean(state.kbt) ≈ 0.0 atol=1e-2
        @test var(state.kbt) ≈ vartest atol=vartest/10
        vartest = 0.2/11
        Swalbe.thermal!(f1D, ones(100000), 0.1, 1/6, 1.0)
        @test mean(f1D) ≈ 0.0 atol=1e-2
        @test var(f1D) ≈ vartest atol=vartest/10
        # Using structs
        sys = Swalbe.SysConst_1D(L=100000,δ=1.0,kbt=0.1)
        state = Swalbe.Sys(sys, kind="thermal")
        Swalbe.thermal!(state, sys)
        @test mean(state.kbt) ≈ 0.0 atol=1e-2
        @test var(state.kbt) ≈ vartest atol=vartest/10
    end

    @testset "rho update" begin
        @testset "Constant fields" begin
            rho = ones(25)
            height = ones(25)
            output = zeros(25)
            Swalbe.update_rho!(rho, output, height, zeros(25,2), zeros(25,4))
            @test all(rho .== 1)
        end
    end

    @testset "view four" begin
        dummy = reshape(collect(1:20),5,4)
        d1, d2, d3, d4 = Swalbe.view_four(dummy)
        @test all(d1 .== dummy[:,1])
        @test all(d2 .== dummy[:,2])
        @test all(d3 .== dummy[:,3])
        @test all(d4 .== dummy[:,4])
    end
end


text:
module DashMantineComponents
using Dash

const resources_path = realpath(joinpath( @__DIR__, "..", "deps"))
const version = "0.2.0"

include("jl/accordion.jl")
include("jl/accordionitem.jl")
include("jl/affix.jl")
include("jl/alert.jl")
include("jl/anchor.jl")
include("jl/badge.jl")
include("jl/blockquote.jl")
include("jl/breadcrumbs.jl")
include("jl/button.jl")
include("jl/center.jl")
include("jl/checkbox.jl")
include("jl/chips.jl")
include("jl/col.jl")
include("jl/container.jl")
include("jl/datepicker.jl")
include("jl/daterangepicker.jl")
include("jl/divider.jl")
include("jl/drawer.jl")
include("jl/grid.jl")
include("jl/group.jl")
include("jl/image.jl")
include("jl/modal.jl")
include("jl/multiselect.jl")
include("jl/notificationhandler.jl")
include("jl/notificationsprovider.jl")
include("jl/paper.jl")
include("jl/prism.jl")
include("jl/progress.jl")
include("jl/radiogroup.jl")
include("jl/segmentedcontrol.jl")
include("jl/select.jl")
include("jl/simplegrid.jl")
include("jl/skeleton.jl")
include("jl/slider.jl")
include("jl/space.jl")
include("jl/spoiler.jl")
include("jl/switch.jl")
include("jl/tab.jl")
include("jl/table.jl")
include("jl/tabs.jl")
include("jl/text.jl")
include("jl/textinput.jl")
include("jl/title.jl")
include("jl/tooltip.jl")

function __init__()
    DashBase.register_package(
        DashBase.ResourcePkg(
            "dash_mantine_components",
            resources_path,
            version = version,
            [
                DashBase.Resource(
    relative_package_path = "dash_mantine_components.min.js",
    external_url = "https://unpkg.com/dash_mantine_components@0.2.0/dash_mantine_components/dash_mantine_components.min.js",
    dynamic = nothing,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dash_mantine_components.min.js.map",
    external_url = "https://unpkg.com/dash_mantine_components@0.2.0/dash_mantine_components/dash_mantine_components.min.js.map",
    dynamic = true,
    async = nothing,
    type = :js
)
            ]
        )

    )
end
end


text:
##############################################################################
##
## Reshaping
##
## Also, see issue # ??
##
##############################################################################

##############################################################################
##
## stack()
## melt()
##
##############################################################################

"""
Stacks a DataTable; convert from a wide to long format

```julia
stack(dt::AbstractDataTable, [measure_vars], [id_vars];
      variable_name::Symbol=:variable, value_name::Symbol=:value)
melt(dt::AbstractDataTable, [id_vars], [measure_vars];
     variable_name::Symbol=:variable, value_name::Symbol=:value)
```

### Arguments

* `dt` : the AbstractDataTable to be stacked

* `measure_vars` : the columns to be stacked (the measurement
  variables), a normal column indexing type, like a Symbol,
  Vector{Symbol}, Int, etc.; for `melt`, defaults to all
  variables that are not `id_vars`. If neither `measure_vars`
  or `id_vars` are given, `measure_vars` defaults to all
  floating point columns.

* `id_vars` : the identifier columns that are repeated during
  stacking, a normal column indexing type; for `stack` defaults to all
  variables that are not `measure_vars`

* `variable_name` : the name of the new stacked column that shall hold the names
  of each of `measure_vars`

* `value_name` : the name of the new stacked column containing the values from
  each of `measure_vars`

### Result

* `::DataTable` : the long-format datatable with column `:value`
  holding the values of the stacked columns (`measure_vars`), with
  column `:variable` a Vector of Symbols with the `measure_vars` name,
  and with columns for each of the `id_vars`.

See also `stackdt` and `meltdt` for stacking methods that return a
view into the original DataTable. See `unstack` for converting from
long to wide format.

### Examples

```julia
d1 = DataTable(a = repeat([1:3;], inner = [4]),
               b = repeat([1:4;], inner = [3]),
               c = randn(12),
               d = randn(12),
               e = map(string, 'a':'l'))

d1s = stack(d1, [:c, :d])
d1s2 = stack(d1, [:c, :d], [:a])
d1m = melt(d1, [:a, :b, :e])
d1s_name = melt(d1, [:a, :b, :e], variable_name=:somemeasure)
```

"""
function stack(dt::AbstractDataTable, measure_vars::Vector{Int},
               id_vars::Vector{Int}; variable_name::Symbol=:variable,
               value_name::Symbol=:value)
    N = length(measure_vars)
    cnames = names(dt)[id_vars]
    insert!(cnames, 1, value_name)
    insert!(cnames, 1, variable_name)
    DataTable(Any[Compat.repeat(_names(dt)[measure_vars], inner=nrow(dt)),   # variable
                  vcat([dt[c] for c in measure_vars]...),                    # value
                  [Compat.repeat(dt[c], outer=N) for c in id_vars]...],      # id_var columns
              cnames)
end
function stack(dt::AbstractDataTable, measure_var::Int, id_var::Int;
               variable_name::Symbol=:variable, value_name::Symbol=:value)
    stack(dt, [measure_var], [id_var];
          variable_name=variable_name, value_name=value_name)
end
function stack(dt::AbstractDataTable, measure_vars::Vector{Int}, id_var::Int;
               variable_name::Symbol=:variable, value_name::Symbol=:value)
    stack(dt, measure_vars, [id_var];
          variable_name=variable_name, value_name=value_name)
end
function stack(dt::AbstractDataTable, measure_var::Int, id_vars::Vector{Int};
               variable_name::Symbol=:variable, value_name::Symbol=:value)
    stackdt(dt, [measure_var], id_vars;
            variable_name=variable_name, value_name=value_name)
end
function stack(dt::AbstractDataTable, measure_vars, id_vars;
               variable_name::Symbol=:variable, value_name::Symbol=:value)
    stack(dt, index(dt)[measure_vars], index(dt)[id_vars];
          variable_name=variable_name, value_name=value_name)
end
# no vars specified, by default select only numeric columns
numeric_vars(dt::AbstractDataTable) =
    [T <: AbstractFloat || (T <: Nullable && eltype(T) <: AbstractFloat)
     for T in eltypes(dt)]

function stack(dt::AbstractDataTable, measure_vars = numeric_vars(dt);
               variable_name::Symbol=:variable, value_name::Symbol=:value)
    mv_inds = index(dt)[measure_vars]
    stack(dt, mv_inds, setdiff(1:ncol(dt), mv_inds);
          variable_name=variable_name, value_name=value_name)
end

"""
Stacks a DataTable; convert from a wide to long format; see
`stack`.
"""
function melt(dt::AbstractDataTable, id_vars::@compat(Union{Int,Symbol});
              variable_name::Symbol=:variable, value_name::Symbol=:value)
    melt(dt, [id_vars]; variable_name=variable_name, value_name=value_name)
end
function melt(dt::AbstractDataTable, id_vars;
              variable_name::Symbol=:variable, value_name::Symbol=:value)
    id_inds = index(dt)[id_vars]
    stack(dt, setdiff(1:ncol(dt), id_inds), id_inds;
          variable_name=variable_name, value_name=value_name)
end
function melt(dt::AbstractDataTable, id_vars, measure_vars;
              variable_name::Symbol=:variable, value_name::Symbol=:value)
    stack(dt, measure_vars, id_vars; variable_name=variable_name,
          value_name=value_name)
end
melt(dt::AbstractDataTable; variable_name::Symbol=:variable, value_name::Symbol=:value) =
    stack(dt; variable_name=variable_name, value_name=value_name)

##############################################################################
##
## unstack()
##
##############################################################################

"""
Unstacks a DataTable; convert from a long to wide format

```julia
unstack(dt::AbstractDataTable, rowkey, colkey, value)
unstack(dt::AbstractDataTable, colkey, value)
unstack(dt::AbstractDataTable)
```

### Arguments

* `dt` : the AbstractDataTable to be unstacked

* `rowkey` : the column with a unique key for each row, if not given,
  find a key by grouping on anything not a `colkey` or `value`

* `colkey` : the column holding the column names in wide format,
  defaults to `:variable`

* `value` : the value column, defaults to `:value`

### Result

* `::DataTable` : the wide-format datatable

### Examples

```julia
wide = DataTable(id = 1:12,
                 a  = repeat([1:3;], inner = [4]),
                 b  = repeat([1:4;], inner = [3]),
                 c  = randn(12),
                 d  = randn(12))

long = stack(wide)
wide0 = unstack(long)
wide1 = unstack(long, :variable, :value)
wide2 = unstack(long, :id, :variable, :value)
```
Note that there are some differences between the widened results above.

"""
function unstack(dt::AbstractDataTable, rowkey::Int, colkey::Int, value::Int)
    # `rowkey` integer indicating which column to place along rows
    # `colkey` integer indicating which column to place along column headers
    # `value` integer indicating which column has values
    refkeycol = NullableCategoricalArray(dt[rowkey])
    valuecol = dt[value]
    keycol = NullableCategoricalArray(dt[colkey])
    Nrow = length(refkeycol.pool)
    Ncol = length(keycol.pool)
    payload = DataTable(Any[similar_nullable(valuecol, Nrow) for i in 1:Ncol], map(Symbol, levels(keycol)))
    nowarning = true
    for k in 1:nrow(dt)
        j = Int(CategoricalArrays.order(keycol.pool)[keycol.refs[k]])
        i = Int(CategoricalArrays.order(refkeycol.pool)[refkeycol.refs[k]])
        if i > 0 && j > 0
            if nowarning && !isnull(payload[j][i])
                warn("Duplicate entries in unstack.")
                nowarning = false
            end
            payload[j][i]  = valuecol[k]
        end
    end
    levs = levels(refkeycol)
    col = similar_nullable(dt[rowkey], length(levs))
    insert!(payload, 1, copy!(col, levs), _names(dt)[rowkey])
end
unstack(dt::AbstractDataTable, rowkey, colkey, value) =
    unstack(dt, index(dt)[rowkey], index(dt)[colkey], index(dt)[value])

# Version of unstack with just the colkey and value columns provided
unstack(dt::AbstractDataTable, colkey, value) =
    unstack(dt, index(dt)[colkey], index(dt)[value])

function unstack(dt::AbstractDataTable, colkey::Int, value::Int)
    # group on anything not a key or value:
    g = groupby(dt, setdiff(_names(dt), _names(dt)[[colkey, value]]), sort=true)
    groupidxs = [g.idx[g.starts[i]:g.ends[i]] for i in 1:length(g.starts)]
    rowkey = zeros(Int, size(dt, 1))
    for i in 1:length(groupidxs)
        rowkey[groupidxs[i]] = i
    end
    keycol = NullableCategoricalArray(dt[colkey])
    valuecol = dt[value]
    dt1 = nullable!(dt[g.idx[g.starts], g.cols], g.cols)
    Nrow = length(g)
    Ncol = length(levels(keycol))
    dt2 = DataTable(Any[similar_nullable(valuecol, Nrow) for i in 1:Ncol], map(Symbol, levels(keycol)))
    nowarning = true
    for k in 1:nrow(dt)
        j = Int(CategoricalArrays.order(keycol.pool)[keycol.refs[k]])
        i = rowkey[k]
        if i > 0 && j > 0
            if nowarning && !isnull(dt2[j][i])
                warn("Duplicate entries in unstack at row $k.")
                nowarning = false
            end
            dt2[j][i]  = valuecol[k]
        end
    end
    hcat(dt1, dt2)
end

unstack(dt::AbstractDataTable) = unstack(dt, :id, :variable, :value)

##############################################################################
##
## Reshaping using referencing (issue #145)
## New AbstractVector types (all read only):
##     StackedVector
##     RepeatedVector
##
##############################################################################

"""
An AbstractVector{Any} that is a linear, concatenated view into
another set of AbstractVectors

NOTE: Not exported.

### Constructor

```julia
StackedVector(d::AbstractVector...)
```

### Arguments

* `d...` : one or more AbstractVectors

### Examples

```julia
StackedVector(Any[[1,2], [9,10], [11,12]])  # [1,2,9,10,11,12]
```

"""
type StackedVector <: AbstractVector{Any}
    components::Vector{Any}
end

function Base.getindex(v::StackedVector,i::Real)
    lengths = [length(x)::Int for x in v.components]
    cumlengths = [0; cumsum(lengths)]
    j = searchsortedlast(cumlengths .+ 1, i)
    if j > length(cumlengths)
        error("indexing bounds error")
    end
    k = i - cumlengths[j]
    if k < 1 || k > length(v.components[j])
        error("indexing bounds error")
    end
    v.components[j][k]
end

function Base.getindex{I<:Real}(v::StackedVector,i::AbstractVector{I})
    result = similar(v.components[1], length(i))
    for idx in 1:length(i)
        result[idx] = v[i[idx]]
    end
    result
end

Base.size(v::StackedVector) = (length(v),)
Base.length(v::StackedVector) = sum(map(length, v.components))
Base.ndims(v::StackedVector) = 1
Base.eltype(v::StackedVector) = promote_type(map(eltype, v.components)...)
Base.similar(v::StackedVector, T, dims::Dims) = similar(v.components[1], T, dims)

CategoricalArrays.CategoricalArray(v::StackedVector) = CategoricalArray(v[:]) # could be more efficient

"""
An AbstractVector that is a view into another AbstractVector with
repeated elements

NOTE: Not exported.

### Constructor

```julia
RepeatedVector(parent::AbstractVector, inner::Int, outer::Int)
```

### Arguments

* `parent` : the AbstractVector that's repeated
* `inner` : the numer of times each element is repeated
* `outer` : the numer of times the whole vector is repeated after
  expanded by `inner`

`inner` and `outer` have the same meaning as similarly named arguments
to `repeat`.

### Examples

```julia
RepeatedVector([1,2], 3, 1)   # [1,1,1,2,2,2]
RepeatedVector([1,2], 1, 3)   # [1,2,1,2,1,2]
RepeatedVector([1,2], 2, 2)   # [1,2,1,2,1,2,1,2]
```

"""
type RepeatedVector{T} <: AbstractVector{T}
    parent::AbstractVector{T}
    inner::Int
    outer::Int
end

function Base.getindex{T,I<:Real}(v::RepeatedVector{T},i::AbstractVector{I})
    N = length(v.parent)
    idx = Int[Base.fld1(mod1(j,v.inner*N),v.inner) for j in i]
    v.parent[idx]
end
function Base.getindex{T}(v::RepeatedVector{T},i::Real)
    N = length(v.parent)
    idx = Base.fld1(mod1(i,v.inner*N),v.inner)
    v.parent[idx]
end
Base.getindex(v::RepeatedVector,i::Range) = getindex(v, [i;])

Base.size(v::RepeatedVector) = (length(v),)
Base.length(v::RepeatedVector) = v.inner * v.outer * length(v.parent)
Base.ndims(v::RepeatedVector) = 1
Base.eltype{T}(v::RepeatedVector{T}) = T
Base.reverse(v::RepeatedVector) = RepeatedVector(reverse(v.parent), v.inner, v.outer)
Base.similar(v::RepeatedVector, T, dims::Dims) = similar(v.parent, T, dims)
Base.unique(v::RepeatedVector) = unique(v.parent)

function CategoricalArrays.CategoricalArray(v::RepeatedVector)
    res = CategoricalArrays.CategoricalArray(v.parent)
    res.refs = repeat(res.refs, inner = [v.inner], outer = [v.outer])
    res
end

##############################################################################
##
## stackdt()
## meltdt()
## Reshaping using referencing (issue #145), using the above vector types
##
##############################################################################

"""
A stacked view of a DataTable (long format)

Like `stack` and `melt`, but a view is returned rather than data
copies.

```julia
stackdt(dt::AbstractDataTable, [measure_vars], [id_vars];
        variable_name::Symbol=:variable, value_name::Symbol=:value)
meltdt(dt::AbstractDataTable, [id_vars], [measure_vars];
       variable_name::Symbol=:variable, value_name::Symbol=:value)
```

### Arguments

* `dt` : the wide AbstractDataTable

* `measure_vars` : the columns to be stacked (the measurement
  variables), a normal column indexing type, like a Symbol,
  Vector{Symbol}, Int, etc.; for `melt`, defaults to all
  variables that are not `id_vars`

* `id_vars` : the identifier columns that are repeated during
  stacking, a normal column indexing type; for `stack` defaults to all
  variables that are not `measure_vars`

### Result

* `::DataTable` : the long-format datatable with column `:value`
  holding the values of the stacked columns (`measure_vars`), with
  column `:variable` a Vector of Symbols with the `measure_vars` name,
  and with columns for each of the `id_vars`.

The result is a view because the columns are special AbstractVectors
that return indexed views into the original DataTable.

### Examples

```julia
d1 = DataTable(a = repeat([1:3;], inner = [4]),
               b = repeat([1:4;], inner = [3]),
               c = randn(12),
               d = randn(12),
               e = map(string, 'a':'l'))

d1s = stackdt(d1, [:c, :d])
d1s2 = stackdt(d1, [:c, :d], [:a])
d1m = meltdt(d1, [:a, :b, :e])
```

"""
function stackdt(dt::AbstractDataTable, measure_vars::Vector{Int},
                 id_vars::Vector{Int}; variable_name::Symbol=:variable,
                 value_name::Symbol=:value)
    N = length(measure_vars)
    cnames = names(dt)[id_vars]
    insert!(cnames, 1, value_name)
    insert!(cnames, 1, variable_name)
    DataTable(Any[RepeatedVector(_names(dt)[measure_vars], nrow(dt), 1),   # variable
                  StackedVector(Any[dt[:,c] for c in measure_vars]),     # value
                  [RepeatedVector(dt[:,c], 1, N) for c in id_vars]...],     # id_var columns
              cnames)
end
function stackdt(dt::AbstractDataTable, measure_var::Int, id_var::Int;
                 variable_name::Symbol=:variable, value_name::Symbol=:value)
    stackdt(dt, [measure_var], [id_var]; variable_name=variable_name,
            value_name=value_name)
end
function stackdt(dt::AbstractDataTable, measure_vars, id_var::Int;
                 variable_name::Symbol=:variable, value_name::Symbol=:value)
    stackdt(dt, measure_vars, [id_var]; variable_name=variable_name,
            value_name=value_name)
end
function stackdt(dt::AbstractDataTable, measure_var::Int, id_vars;
                 variable_name::Symbol=:variable, value_name::Symbol=:value)
    stackdt(dt, [measure_var], id_vars; variable_name=variable_name,
            value_name=value_name)
end
function stackdt(dt::AbstractDataTable, measure_vars, id_vars;
                 variable_name::Symbol=:variable, value_name::Symbol=:value)
    stackdt(dt, index(dt)[measure_vars], index(dt)[id_vars];
            variable_name=variable_name, value_name=value_name)
end
function stackdt(dt::AbstractDataTable, measure_vars = numeric_vars(dt);
                 variable_name::Symbol=:variable, value_name::Symbol=:value)
    m_inds = index(dt)[measure_vars]
    stackdt(dt, m_inds, setdiff(1:ncol(dt), m_inds);
            variable_name=variable_name, value_name=value_name)
end

"""
A stacked view of a DataTable (long format); see `stackdt`
"""
function meltdt(dt::AbstractDataTable, id_vars; variable_name::Symbol=:variable,
                value_name::Symbol=:value)
    id_inds = index(dt)[id_vars]
    stackdt(dt, setdiff(1:ncol(dt), id_inds), id_inds;
            variable_name=variable_name, value_name=value_name)
end
function meltdt(dt::AbstractDataTable, id_vars, measure_vars;
                variable_name::Symbol=:variable, value_name::Symbol=:value)
    stackdt(dt, measure_vars, id_vars; variable_name=variable_name,
            value_name=value_name)
end
meltdt(dt::AbstractDataTable; variable_name::Symbol=:variable, value_name::Symbol=:value) =
    stackdt(dt; variable_name=variable_name, value_name=value_name)


text:
module PkgUtility

using CLIMAParameters
using CLIMAParameters: AbstractEarthParameterSet, avogad, gas_constant, h_Planck, k_Boltzmann, light_speed
using CLIMAParameters.Planet: LH_v0, MSLP, R_v, T_freeze, T_triple, cp_d, cp_l, cp_v, grav, molmass_dryair, molmass_water, press_triple, ρ_cloud_liq
using CLIMAParameters.SubgridScale: von_karman_const
using CSV: File, write
using DataFrames: DataFrame
using Dates: Date, DateTime, format, isleapyear, now
using DocStringExtensions: METHODLIST
using NCDatasets: Dataset, defDim, defVar, listVar
using Pkg.Artifacts: archive_artifact, artifact_exists, artifact_hash, artifact_meta, artifact_path, bind_artifact!, create_artifact, download_artifact
using Pkg.PlatformEngines: unpack
using Statistics: mean, median, std
using StatsBase: percentile

# global constants
const MDAYS_LEAP  = [0,31,60,91,121,152,182,213,244,274,305,335,366];
const MDAYS       = [0,31,59,90,120,151,181,212,243,273,304,334,365];
const NDAYS_LEAP  = [31,29,31,30,31,30,31,31,30,31,30,31];
const NDAYS       = [31,28,31,30,31,30,31,31,30,31,30,31];
const TIME_FORMAT = ["YYYYMMDD", "YYYYMMDDhh", "YYYYMMDDhhmm", "YYYYMMDDhhmmss"];
const TIME_OUTPUT = ["DATE", "DATETIME", "DOY", "FDOY"];

# export public functions --- artifact
export deploy_artifact!, predownload_artifact!

# export public functions --- date
export month_days, month_ind, parse_timestamp, terror, tinfo, twarn

# export public functions --- display
export pretty_display!

# export public functions --- io
export dataframe, read_csv, read_nc, save_csv!, save_nc!, size_nc, send_email!

# export land wrapper functions --- land
export AVOGADRO, CP_D, CP_D_MOL, CP_L, CP_V, GAS_R, GRAVITY, H_PLANCK, K_BOLTZMANN, K_STEFAN, K_VON_KARMAN, LH_V0, LIGHT_SPEED, M_DRYAIR, M_H₂O, P_ATM,
       PRESS_TRIPLE, R_V, RT_25, T_0, T_25, T_TRIPLE, V_H₂O, YEAR_D, ρ_H₂O, ρg_MPa

# export public functions --- math
export lower_quadratic, mae, mape, mase, nanmax, nanmean, nanmedian, nanmin, nanpercentile, nanstd, numerical∫, rmse, upper_quadratic

# export public functions --- test
export FT_test, NaN_test

include("artifact/deploy.jl"  )
include("artifact/download.jl")

include("date/doy.jl"  )
include("date/tinfo.jl")

include("display/recursive.jl")

include("io/csv.jl"      )
include("io/dataframe.jl")
include("io/email.jl"    )
include("io/nc.jl"       )

include("land/clima.jl")

include("math/integral.jl"  )
include("math/quadratic.jl" )
include("math/statistics.jl")

include("test/recursive.jl")

include("deprecated.jl")

end # module


text:
using Compat, NullableArrays

immutable CSVError <: Exception
    msg::ASCIIString
end

# used only during the type detection process
immutable NullField
end

const ZIP_FILE_EXTS = Set{UTF8String}([".zip",".gz","-gz",".z","-z","_z",".Z",".bz2", ".bz", ".tbz2", ".tbz"])

function unzip(ext,file)
    if ext in (".Z",".gz",".zip")
        return `gzip -dc $file`
    elseif ext in (".bz2", ".bz", ".tbz2", ".tbz")
        return `bzip2 -dck $file`
    else
        throw(ArgumentError("ERROR: unsupported compression file format: $ext"))
    end
    return cmd
end

const RETURN  = @compat UInt8('\r')
const NEWLINE = @compat UInt8('\n')
const COMMA   = @compat UInt8(',')
const QUOTE   = @compat UInt8('"')
const ESCAPE  = @compat UInt8('\\')
const PERIOD  = @compat UInt8('.')
const SPACE   = @compat UInt8(' ')
const TAB     = @compat UInt8('\t')
const MINUS   = @compat UInt8('-')
const PLUS    = @compat UInt8('+')
const NEG_ONE = @compat UInt8('0')-UInt8(1)
const ZERO    = @compat UInt8('0')
const TEN     = @compat UInt8('9')+UInt8(1)
Base.isascii(c::UInt8) = c < 0x80

# accept URLs?
# ignore columns
type File
    fullpath::UTF8String
    compression::UTF8String

    delim::UInt8
    quotechar::UInt8
    escapechar::UInt8
    separator::UInt8
    decimal::UInt8
    null::UTF8String

    cols::Int
    header::Vector{UTF8String}
    types::Vector{DataType}
    formats::Vector{UTF8String}

    datapos::Int
    datarow::Int
    skipblankrows::Bool
    footerskip::Int
end

function Base.show(io::IO,f::File)
    println(io,"fullpath: \"",f.fullpath,"\"")
    println(io,"compression: ",f.compression)

    println(io,"delim: '",@compat(Char(f.delim)),"'")
    println(io,"quotechar: '",@compat(Char(f.quotechar)),"'")
    println(io,"escapechar: '\\",@compat(Char(f.escapechar)),"'")
    println(io,"separator: '",@compat(Char(f.separator)),"'")
    println(io,"decimal: '",@compat(Char(f.decimal)),"'")
    println(io,"null: \"",f.null,"\"")

    println(io,"cols: ",f.cols)
    println(io,"header: ",f.header)
    println(io,"types: ",f.types)
    println(io,"formats: ",f.formats)

    println(io,"datapos: ",f.datapos)
    println(io,"skipblankrows: ",f.skipblankrows)
    println(io,"footerskip: ",f.footerskip)
end

function File{T<:AbstractString}(fullpath::AbstractString;
              compression="",

              delim=COMMA,
              quotechar=QUOTE,
              escapechar=ESCAPE,
              separator=COMMA,
              decimal=PERIOD,
              null::AbstractString="",

              header::Union(Integer,UnitRange{Int},Vector{T})=1, # header can be a row number, range of rows, or actual string vector
              datarow::Int=2, # by default, data starts immediately after header or start of file
              types::Union(Dict{Int,DataType},Vector{DataType})=DataType[],
              formats::Union(Dict{Int,T},Vector{T})=UTF8String[],

              skipblankrows::Bool=true,
              footerskip::Int=0,
              rows_for_type_detect::Int=250)
    # argument checks
    isfile(fullpath) || throw(ArgumentError("\"$fullpath\" is not a valid file"))
    isa(header,Integer) && (datarow > header || throw(ArgumentError("data row ($datarow) must come after header row ($header)")))

    # open the file for property detection; handle possible compression types
    if compression in ZIP_FILE_EXTS || splitext(fullpath)[end] in ZIP_FILE_EXTS
        temp = tempname()
        splitext(fullpath)[end] in ZIP_FILE_EXTS && (compression = splitext(fullpath)[end])
        unzipcmd = unzip(splitext(fullpath)[end],fullpath)
        run(pipe(pipe(unzipcmd,`head -$(rows_for_type_detect+datarow)`);stdout=temp))
        f = open(temp)
        rm(temp)
    elseif compression == ""
        f = open(fullpath)
    else
        throw(ArgumentError("unsupported compression type: $compression"))
    end

    # make sure character args are UInt8
    isascii(delim) || throw(ArgumentError("non-ASCII characters not supported for delim argument: $delim"))
    isascii(quotechar) || throw(ArgumentError("non-ASCII characters not supported for quotechar argument: $quotechar"))
    isascii(escapechar) || throw(ArgumentError("non-ASCII characters not supported for escapechar argument: $escapechar"))
    isascii(decimal) || throw(ArgumentError("non-ASCII characters not supported for decimal argument: $decimal"))
    isascii(separator) || throw(ArgumentError("non-ASCII characters not supported for separator argument: $separator"))
    delim = delim % UInt8
    quotechar = quotechar % UInt8
    escapechar = escapechar % UInt8
    decimal = decimal % UInt8
    separator = separator % UInt8

    # figure out # of columns and header, either an Integer, Range, or Vector{String}
    # also ensure that `f` is positioned at the start of data
    cols = 0
    if isa(header,Integer)
        CSV.skipto!(f,1,header,quotechar,escapechar)
        headerpos = position(f)
        cols = length(split(readline(f,quotechar),@compat(Char(delim))))
        seek(f,headerpos)
        columnnames = map(utf8,split(readline(f,quotechar),@compat(Char(delim))))
        datarow != 0 && CSV.skipto!(f,header+1,datarow,quotechar,escapechar)
        datapos = position(f)
    elseif isa(header,Range)
        CSV.skipto!(f,1,first(header),quotechar,escapechar)
        headerpos = position(f)
        cols = length(split(readline(f,quotechar),@compat(Char(delim))))
        seek(f,headerpos)
        columnames = fill(utf8(""),cols)
        for row in header
            for (i,c) in map(utf8,split(readline(f,quotechar),@compat(Char(delim))))
                columnames[i] *= "_" * c
            end
        end
        datarow != 0 && CSV.skipto!(f,last(header)+1,datarow,quotechar,escapechar)
        datapos = position(f)
    else
        datarow = datarow == 0 ? 1 : datarow
        CSV.skipto!(f,1,datarow,quotechar,escapechar)
        datapos = position(f)
        cols = length(split(readline(f,quotechar),@compat(Char(delim))))
        seek(f,datapos)
        if isempty(header)
            columnnames = ["Column$i" for i = 1:cols]
        else
            length(header) == cols || throw(ArgumentError("length of provided header doesn't match number of columns of data at row $datarow"))
            columnnames = header
        end
    end

    # Detect column types
    columnformats = fill(utf8(""),cols)
    if isa(formats,Vector) && length(formats) == cols
        columnformats = formats
    elseif isa(formats,Dict) || isempty(formats)
        for (k,v) in formats
            columnformats[k] = v
        end
    else
        throw(ArgumentError("$cols number of columns detected; `formats` argument has $(length(formats)) entries"))
    end
    columntypes = Array(DataType,cols)
    if isa(types,Vector) && length(types) == cols
        columntypes = types
    elseif isa(types,Dict) || isempty(types)
        poss_types = Array(DataType,rows_for_type_detect,cols)
        lineschecked = 0
        while !eof(f) && lineschecked < rows_for_type_detect
            vals = split(readline(f,quotechar),@compat(Char(delim)))
            lineschecked += 1
            for i = 1:cols
               poss_types[lineschecked,i], columnformats[i] = CSV.detecttype(vals[i],columnformats[i],@compat(Char(separator)),null)
            end
        end
        # detect most common/general type of each column of types
        d = Set{DataType}()
        for i = 1:cols
            for n = 1:lineschecked
                t = poss_types[n,i]
                t == NullField && continue
                push!(d,t)
            end
            columntypes[i] = (isempty(d) || AbstractString in d ) ? AbstractString :
                (Date     in d) ? Date :
                (DateTime in d) ? DateTime :
                (Float64  in d) ? Float64 :
                (Int      in d) ? Int : AbstractString
            empty!(d)
        end
    else
        throw(ArgumentError("$cols number of columns detected; `types` argument has $(length(types)) entries"))
    end
    if isa(types,Dict)
        for (col,typ) in types
            columptypes[col] = typ
        end
    end
    return File(utf8(fullpath),compression,
                delim,quotechar,escapechar,separator,decimal,utf8(null),
                cols,columnnames,columntypes,columnformats,
                datapos,datarow,skipblankrows,footerskip)
end

# for potentially embedded newlines in quoted fields
function readline(f::IOStream,q::UInt8)
    buf = IOBuffer()
    while !eof(f)
        b = read(f, UInt8)
        write(buf, b)
        if b == q
            while !eof(f)
                b = read(f, UInt8)
                write(buf, b)
                b == q && break
            end
        elseif b == NEWLINE
            break
        elseif b == RETURN
            mark(f)
            read(f, UInt8) != NEWLINE && reset(f)
            break
        end
    end
    return takebuf_string(buf)
end

# for skipping from `cur` line number to `to`
function skipto!(f,cur,to,q,e)
    while !eof(f)
        cur == to && break
        b = read(f, UInt8)
        if b == q
            while !eof(f)
                b = read(f, UInt8)
                if b == e
                    b = read(f, UInt8)
                elseif b == q
                    break
                end
            end
        elseif b == NEWLINE
            cur += 1
        elseif b == RETURN
            mark(f)
            read(f, UInt8) != NEWLINE && reset(f)
            cur += 1
        end
    end
    return cur
end

function detecttype(val,f::AbstractString,sep::Char,null::AbstractString)
    (val == "" || val == null) && return NullField, ""
    val2 = replace(val, sep, "")
    t = tryparse(Int,val2)
    !isnull(t) && return Int, ""
    t = tryparse(Float64,val2)
    !isnull(t) && return Float64, ""
    if f != ""
        try
            Date(val,f)
            return Date, f
        end
        try
            DateTime(val,f)
            return DateTime, f
        end
    end
    return Str, ""
end

# we need to keep a reference to our mmapped arrays for Strs
const GLOBALREF = []

type Stream <: IO
    file::CSV.File
    data::Vector{UInt8} # Mmap.Array
    pos::Int
end
function Base.open(file::CSV.File)
    m = Mmap.mmap(file.fullpath)
    push!(GLOBALREF,m)
    Stream(file,m,1)
end
# CSV.open(file::AbstractString) = Stream(CSV.File(file),mmap(file),1)

@inline function Base.read(io::CSV.Stream, ::Type{UInt8}=UInt8)
    @inbounds b = io.data[io.pos]
    io.pos += 1
    return b
end
peek(io::CSV.Stream) = (@inbounds b = io.data[io.pos]; return b)
Base.eof(io::CSV.Stream) = io.pos > length(io.data)
Base.position(io::CSV.Stream) = io.pos
Base.seek(io::CSV.Stream,i::Int) = (io.pos = i)
function skipn!(f::CSV.Stream,n,q,e)
    cur = 0
    while !eof(f)
        cur == n && break
        b = read(f)
        if b == q
            while !eof(f)
                b = read(f)
                if b == e
                    b = read(f)
                elseif b == q
                    break
                end
            end
        elseif b == NEWLINE
            cur += 1
        elseif b == RETURN
            peek(f) == NEWLINE && (f.pos+=1)
            cur += 1
        end
    end
    return cur
end

immutable Str <: AbstractString
    ptr::Ptr{UInt8}
    len::Int
end
const NULLSTRING = Str(C_NULL,0)
Base.print(io::IO, x::Str) = show(io,x)
Base.show(io::IO, x::Str) = print(io,x == NULLSTRING ? "" : "\"$(bytestring(x.ptr,x.len))\"")
Base.endof(x::Str) = x.len
Base.string(x::Str) = x == NULLSTRING ? "" : bytestring(x.ptr,x.len)

@inline function readfield{T<:Integer}(io::CSV.Stream, ::Type{T}, row, col)
    pos = io.pos
    @inbounds begin
    b = read(io)
    while !eof(io) && (b == SPACE || b == TAB || b == io.file.quotechar)
        b = read(io)
    end
    if b == io.file.delim || b == NEWLINE
        return zero(T), true
    elseif b == RETURN
        peek(io) == NEWLINE && (io.pos+=1)
        return zero(T), true
    end
    negative = false
    if b == MINUS
        negative = true
        b = read(io)
    elseif b == PLUS
        b = read(io)
    end
    v = zero(T)
    while NEG_ONE < b < TEN
        # process digits
        v *= 10
        v += b - ZERO
        eof(io) && break
        b = read(io)
    end
    end # @inbounds
    if b == io.file.delim || b == NEWLINE || eof(io)
        return ifelse(negative,-v,v), false
    elseif b == RETURN
        peek(io) == NEWLINE && (io.pos+=1)
        return ifelse(negative,-v,v), false
    else
        throw(CSV.CSVError("error parsing $T on column $col, row $row; parsed $v before encountering $(@compat(Char(b))) character"))
    end
end

const REF = Array(Ptr{UInt8},1)

@inline function readfield(io::CSV.Stream, ::Type{Float64}, row, col)
    b = read(io)
    while !eof(io) && (b == CSV.SPACE || b == CSV.TAB || b == io.file.quotechar)
        b = read(io)
    end
    if b == io.file.delim || b == NEWLINE # || b == PERIOD
        return NaN, true
    elseif b == RETURN
        peek(io) == NEWLINE && (io.pos+=1)
        return NaN, true
    end
    ptr = pointer(io.data) + UInt(io.pos-2)
    v = ccall(:strtod, Float64, (Ptr{UInt8},Ptr{Ptr{UInt8}}), ptr, REF)
    io.pos += REF[1] - ptr
    b = io.data[io.pos-1]
    if b == io.file.delim || b == NEWLINE || eof(io)
        return v, false
    elseif b == RETURN
        peek(io) == NEWLINE && (io.pos+=1)
        return v, false
    else
        throw(CSV.CSVError("error parsing Float64 on row $row, column $col; parsed '$v' before encountering $(@compat(Char(b))) character"))
    end
end

const UINT8NULL = convert(Ptr{UInt8},C_NULL)

@inline function readfield{T<:AbstractString}(io::CSV.Stream, ::Type{T}, row, col)
    pos = position(io)
    @inbounds while !eof(io)
        b = read(io)
        if b == io.file.quotechar
            while !eof(io)
                b = read(io)
                if b == io.file.escapechar
                    b = read(io)
                elseif b == io.file.quotechar
                    break
                end
            end
        elseif b == io.file.delim || b == NEWLINE
            break
        elseif b == RETURN
            peek(io) == NEWLINE && (io.pos+=1)
            break
        end
    end
    if pos == position(io)-1
        return NULLSTRING, true
    else
        return Str(pointer(io.data)+Uint(pos-1), position(io)-pos-1), false
    end
end

itr(io,n,val) = (for i = 1:n; val += 10; val += read(io) - ZERO; end; return val)

@inline function readfield(io::CSV.Stream, ::Type{Date}, row, col)
    year = itr(io,4,0)
    read(io)
    month = itr(io,2,0)
    read(io)
    day = itr(io,2,0)
    read(io)
    return Date(year,month,day), false
end

gettype(x) = x
gettype{T<:AbstractString}(::Type{T}) = Str

# add skip, limit, ignore columns
function Base.read(file::CSV.File)
    io = CSV.open(file)
    seek(io,file.datapos+1)
    rows = countlines(file.fullpath)
    N = file.datarow
    cols = file.cols
    result = Any[]
    for i = 1:cols
        push!(result,NullableArray(CSV.gettype(file.types[i]), rows - N + 1))
    end
    ind = 1
    while !eof(io)
        for i = 1:cols
            @inbounds result[i][ind] = CSV.readfield(io, file.types[i], N, i)
        end
        N += 1; ind += 1
        b = peek(io)
        empty = b == NEWLINE || b == RETURN
        if empty
            N == rows && break
            file.skipblankrows && skipn!(io,1,file.quotechar,file.escapechar)
        end
    end
    return result
end

function Base.setindex!{T}(a::NullableVector{T},x::Tuple{T,Bool},i::Int)
    a.values[i], a.isnull[i] = x
    return a
end


text:
using Documenter
#include("../src/LightGraphs.jl")
using LightGraphs

# same for contributing and license
cp(normpath(@__FILE__, "../../CONTRIBUTING.md"), normpath(@__FILE__, "../src/contributing.md"); force=true)
cp(normpath(@__FILE__, "../../LICENSE.md"), normpath(@__FILE__, "../src/license.md"); force=true)

makedocs(
    modules     = [LightGraphs],
    format      = Documenter.HTML(),
    sitename    = "LightGraphs",
    doctest     = false,
    pages       = Any[
        "Getting Started"                   => "index.md",
        "Choosing A Graph Type"             => "graphtypes.md",
        "LightGraphs Types"                 => "types.md",
        "Accessing Properties"              => "basicproperties.md",
        "Making and Modifying Graphs"       => "generators.md",
        "Reading / Writing Graphs"          => "persistence.md",
        "Operators"                         => "operators.md",
        "Plotting Graphs"                   => "plotting.md",
        "Path and Traversal"                => "pathing.md",
        "Coloring"                          => "coloring.md",
        "Distance"                          => "distance.md",
        "Centrality Measures"               => "centrality.md",
        "Linear Algebra"                    => "linalg.md",
        "Matching"                          => "matching.md",
        "Community Structures"              => "community.md",
        "Degeneracy"                        => "degeneracy.md",
        "Integration with other packages"   => "integration.md",
        "Experimental Functionality"        => "experimental.md",
        "Parallel Algorithms"               => "parallel.md",
        "Contributing"                      => "contributing.md",
        "Developer Notes"                   => "developing.md",
        "License Information"               => "license.md",
        "Citing LightGraphs"                => "citing.md"
    ]
)

deploydocs(
    deps        = nothing,
    make        = nothing,
    repo        = "github.com/JuliaGraphs/LightGraphs.jl.git",
    target      = "build",
    julia       = "nightly",
    osname      = "linux"
)

rm(normpath(@__FILE__, "../src/contributing.md"))
rm(normpath(@__FILE__, "../src/license.md"))


text:
using BenchmarkTools

list_test_str = "199
200
208
210
200
207
240
269
260
263"

list_test_str_split = split(list_test_str, '\n')
list_test = parse.(Int, list_test_str_split)

list_str = readlines("day1_input.txt")
list = parse.(Int, list_str)

# Part 1

f(l) = count(l[n+1] > l[n] for n in 1:length(l)-1)

@info "Part 1"
@info "Test successful?: $(f(list_test)== 7)"
@info "Challenge successful?: $(f(list) == 1655)"
@btime f($list_test)

# Part 2

g(l) = count(sum(@view(l[n+1:n+3])) > sum(@view(l[n:n+2])) for n in 1:length(l)-3)

@info "Part 2"
@info "Test successful?: $(g(list_test) == 5)"
@info "Challenge successful?: $(g(list) == 1683)"

@btime g($list)

# Code golfing

res = parse.(Int, readlines("day1_input.txt")) |> l -> count(l[n+3] > l[n] for n in 1:length(l)-3)
@info "Code golfing success?: $(res == 1683)"

# Performance

function h(list)
    s = 0
    for j in 1:lastindex(list)-3
        @inbounds s += ifelse(list[j+3] > list[j], 1, 0)
    end
    return s
end

@info "Performance"
@btime h($list)


text:
function eye(middlesizes::Matrix{Int64})

    cores = Vector{Array{Float64,4}}(undef,size(middlesizes,2));
    for i = 1:size(middlesizes,2)
        core = Matrix(I,(middlesizes[1,i],middlesizes[2,i]));
        core = reshape(core, (1, middlesizes[1,i],middlesizes[2,i], 1) );
        cores[i] = convert(Array{Float64}, core);
    end

    MPT(cores)
end

max_stars_repo_name: clarazen/BigMat

max_issues_repo_name: clarazen/BigMat

max_forks_repo_name: clarazen/BigMat


text:
function splat!(a::Image, x::Int, y::Int, c::Color)
    a.pic[x, y] = c
    nothing
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

struct ParameterMappedDensity{
    D<:AbstractDensity,
    M<:ParameterMapping,
    B<:AbstractParamBounds,
} <: AbstractDensity
    orig_density::D
    parmap::M
    new_bounds::B
end

ParameterMappedDensity(density::AbstractDensity, parmap::ParameterMapping) =
    ParameterMappedDensity(density, parmap, invmap_params(parmap, param_bounds(density)))

import Base.∘
∘(density::AbstractDensity, parmap::ParameterMapping) =
    ParameterMappedDensity(density, parmap)

Base.parent(density::ParameterMappedDensity) = density.orig_density

BAT.param_bounds(density::ParameterMappedDensity) = density.new_bounds

BAT.nparams(density::ParameterMappedDensity) = nparams(density.new_bounds)

function BAT.density_logval(density::ParameterMappedDensity, params::AbstractVector{<:Real})
    BAT.density_logval(density.orig_density, map_params(density.parmap, params))
end

import Base.∘
∘(density::ConstDensity{<:HyperRectBounds}, parmap::ParameterMapping) =
    ConstDensity(invmap_params(parmap, density.bounds), density.log_value)


text:
# Code for Video 2.1:

Array{Float64, 2}
Dict{Symbol, Float64}

# Array{T, 2}
# Dict{K, V}

function double(f::Float64)
    f * 2
end

local n::Int32 = 7
# (31 + 42)::Float64
# ERROR: TypeError: typeassert: expected Float64, got Int64

super(Int64) # Signed

subtypes(Signed)
# 5-element Array{Any,1}:
#  Int128
#  Int16
#  Int32
#  Int64
#  Int8

Int32 <: Signed # true
Float64 <: Number # true

type Iris
    sepal_length::Float64
    sepal_width::Float64
    petal_length::Float64
    petal_width::Float64
    species::ASCIIString
end

fieldnames(Iris)
# 5-element Array{Symbol,1}:
#  :sepal_length
#  :sepal_width
#  :petal_length
#  :petal_width
#  :species

iris1 = Iris(5.1, 3.5, 1.4, 0.2, "I. setosa")
# Iris(5.1,3.5,1.4,0.2,"I. setosa")
iris1.petal_width # 0.2
iris1.petal_width = 0.3 # 0.3
iris1 # Iris(5.1,3.5,1.4,0.3,"I. setosa")

immutable Iris2
    sepal_length::Float64
    sepal_width::Float64
    petal_length::Float64
    petal_width::Float64
    species::ASCIIString
end
iris2 = Iris2(5.1, 3.5, 1.4, 0.2, "I. setosa")
# iris1.petal_width = 0.2 # ERROR: LoadError: type Iris2 is immutable

arr_iris = Iris[] # 0-element Array{Iris,1}
push!(arr_iris, iris1)
# 1-element Array{Iris,1}:
#  Iris(5.1,3.5,1.4,0.2,"I. setosa")

type IrisG{T <: Real}
   sepal_length::T
   sepal_width::T
   petal_length::T
   petal_width::T
   species::ASCIIString
end

add{T}(x::T, y::T) = x + y
# add (generic function with 1 method)
add{T <: Number}(x::T, y::T) = x + y

# add("I. setosa", "I. versicolor")
# ERROR: MethodError: `+` has no method matching +(::ASCIIString, ::ASCIIString)
# Closest candidates are:
#   +(::Any, ::Any, ::Any, ::Any...)
#  in add at none:1

function add(x::Float64, y::Float64)
    println("Adding 2 floating point numbers:")
    println(x + y)
    x + y
end

methods(add)
# 3 methods for generic function "add":
# add(x::Float64, y::Float64) at none:2
# add{T<:Number}(x::T<:Number, y::T<:Number) at none:1
# add{T}(x::T, y::T) at none:1

add(2, 5) # 7
add(3.1, 5.4)
# Adding 2 floating point numbers:
# 8.5

methods(+)
# 171 methods for generic function "+":
# +(x::Bool) at bool.jl:33
# +(x::Bool, y::Bool) at bool.jl:36
# +(y::AbstractFloat, x::Bool) at bool.jl:46
# +(x::Int64, y::Int64) at int.jl:8
# +(x::Int8, y::Int8) at int.jl:16
# +(x::UInt8, y::UInt8) at int.jl:16
# ....

type Setosa
    # some fields
end

type Versicolor
    # some fields
end

type Virginica
    # some fields
end

IrisU = Union{Setosa, Versicolor, Virginica}
Union{Setosa,Versicolor,Virginica}

set1 = Setosa() # Setosa() value
isa(set1, IrisU) # true

max_stars_repo_path: Section 2/2_1.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 2/2_1.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 2/2_1.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
using LinearAlgebra, MixedModels, Random, SparseArrays, Test

@testset "utilities" begin
	@test MixedModels.average(1.1, 1.2) == 1.15
	@test MixedModels.densify(sparse(1:5, 1:5, ones(5))) == Diagonal(ones(5))
	rsparsev = SparseVector(float.(rand(MersenneTwister(123454321), Bool, 20)))
	@test MixedModels.densify(rsparsev) == Vector(rsparsev)
	@test MixedModels.densify(Diagonal(rsparsev)) == Diagonal(Vector(rsparsev))
end

@testset "threaded_replicate" begin
	rng = MersenneTwister(42);
	single_thread = replicate(10,use_threads=false) do; randn(rng, 1)[1] ; end
	rng = MersenneTwister(42);
	multi_thread = replicate(10,use_threads=true) do
		if Threads.threadid() % 2 == 0
			sleep(0.001)
		end
		r = randn(rng, 1)[1];
	end

	@test all(sort!(single_thread) .== sort!(multi_thread))
end

@testset "datasets" begin
	@test isa(MixedModels.datasets(), Vector{String})
	@test size(MixedModels.dataset(:dyestuff)) == (30, 2)
	@test size(MixedModels.dataset("dyestuff")) == (30, 2)
	@test_throws ArgumentError MixedModels.dataset(:foo)
end