text:
module Tables

using LinearAlgebra, DataValueInterfaces, DataAPI, TableTraits, IteratorInterfaceExtensions

export rowtable, columntable

if !hasmethod(getproperty, Tuple{Tuple, Int})
    Base.getproperty(t::Tuple, i::Int) = t[i]
end

"Abstract row type with a simple required interface: row values are accessible via `getproperty(row, field)`; for example, a NamedTuple like `nt = (a=1, b=2, c=3)` can access its value for `a` like `nt.a` which turns into a call to the function `getproperty(nt, :a)`"
abstract type Row end

"""
The Tables.jl package provides simple, yet powerful interface functions for working with all kinds tabular data through predictable access patterns.

```julia
    Tables.rows(table) => Rows
    Tables.columns(table) => Columns
```
Where `Rows` and `Columns` are the duals of each other:
* `Rows` is an iterator of property-accessible objects (any type that supports `propertynames(row)` and `getproperty(row, nm::Symbol`)
* `Columns` is a property-accessible object of iterators (i.e. each column is an iterator)

In addition to these `Rows` and `Columns` objects, it's useful to be able to query properties of these objects:
* `Tables.schema(x::Union{Rows, Columns}) => Union{Tables.Schema, Nothing}`: returns a `Tables.Schema` object, or `nothing` if the table's schema is unknown
* For the `Tables.Schema` object:
  * column names can be accessed as a tuple of Symbols like `sch.names`
  * column types can be accessed as a tuple of types like `sch.types`
  * See `?Tables.Schema` for more details on this type

A big part of the power in these simple interface functions is that each (`Tables.rows` & `Tables.columns`) is defined for any table type, even if the table type only explicitly implements one interface function or the other.
This is accomplished by providing performant, generic fallback definitions in Tables.jl itself (though obviously nothing prevents a table type from implementing each interface function directly).

With these simple definitions, powerful workflows are enabled:
* A package providing data cleansing, manipulation, visualization, or analysis can automatically handle any number of decoupled input table types
* A tabular file format can have automatic integration with in-memory structures and translation to other file formats

So how does one go about satisfying the Tables.jl interface functions? It mainly depends on what you've already defined and the natural access patterns of your table:

First:
* `Tables.istable(::Type{<:MyTable}) = true`: this provides an explicit affirmation that your type implements the Tables interface

To support `Rows`:
* Define `Tables.rowaccess(::Type{<:MyTable}) = true`: this signals to other types that `MyTable` supports valid `Row`-iteration
* Define `Tables.rows(x::MyTable)`: return a `Row`-iterator object (perhaps the table itself if already defined)
* Define `Tables.schema(Tables.rows(x::MyTable))` to either return a `Tables.Schema` object, or `nothing` if the schema is unknown or non-inferrable for some reason

To support `Columns`:
* Define `Tables.columnaccess(::Type{<:MyTable}) = true`: this signals to other types that `MyTable` supports returning a valid `Columns` object
* Define `Tables.columns(x::MyTable)`: return a `Columns`, property-accessible object (perhaps the table itself if it naturally supports property-access to columns)
* Define `Tables.schema(Tables.columns(x::MyTable))` to either return a `Tables.Schema` object, or `nothing` if the schema is unknown or non-inferrable for some reason

The final question is how `MyTable` can be a "sink" for any other table type. The answer is quite simple: use the interface functions!

* Define a function or constructor that takes, at a minimum, a single, untyped argument and then calls `Tables.rows` or `Tables.columns` on that argument to construct an instance of `MyTable`

For example, if `MyTable` is a row-oriented format, I might define my "sink" function like:
```julia
function MyTable(x)
    Tables.istable(x) || throw(ArgumentError("MyTable requires a table input"))
    rows = Tables.rows(x)
    sch = Tables.schema(rows)
    names = sch.names
    types = sch.types
    # custom constructor that creates an "empty" MyTable according to given column names & types
    # note that the "unknown" schema case should be considered, i.e. when `sch.types => nothing`
    mytbl = MyTable(names, types)
    for row in rows
        # a convenience function provided in Tables.jl for "unrolling" access to each column/property of a `Row`
        # it works by applying a provided function to each value; see `?Tables.eachcolumn` for more details
        Tables.eachcolumn(sch, row) do val, col, name
            push!(mytbl[col], val)
        end
    end
    return mytbl
end
```

Alternatively, if `MyTable` is column-oriented, perhaps my definition would be more like:
```julia
function MyTable(x)
    Tables.istable(x) || throw(ArgumentError("MyTable requires a table input"))
    cols = Tables.columns(x)
    # here we use Tables.eachcolumn to iterate over each column in a `Columns` object
    return MyTable(collect(propertynames(cols)), [collect(col) for col in Tables.eachcolumn(cols)])
end
```

Obviously every table type is different, but via a combination of `Tables.rows` and `Tables.columns` each table type should be able to construct an instance of itself.
"""
abstract type Table end

# default definitions
istable(x::T) where {T} = istable(T) || TableTraits.isiterabletable(x) === true
istable(::Type{T}) where {T} = false
rowaccess(x::T) where {T} = rowaccess(T)
rowaccess(::Type{T}) where {T} = false
columnaccess(x::T) where {T} = columnaccess(T)
columnaccess(::Type{T}) where {T} = false
schema(x) = nothing
materializer(x) = columntable

# Schema implementation
"""
    Tables.Schema(names, types)

Create a `Tables.Schema` object that holds the column names and types for a tabular data object.
`Tables.Schema` is dual-purposed: provide an easy interface for users to query these properties,
as well as provide a convenient "structural" type for code generation.

To get a table's schema, one can call `Tables.schema(tbl)`, but also note that a table may return `nothing`,
indicating that it's column names and/or column types are unknown (usually not inferrable). This is similar
to the `Base.EltypeUnknown()` trait for iterators when `Base.IteratorEltype` is called. Users should account
for the `Tables.schema(tbl) => nothing` case by using the properties of the results of `Tables.rows(x)` and `Tables.columns(x)`
directly.

To access the names, one can simply call `sch.names` to return the tuple of Symbols.
To access column types, one can similarly call `sch.types`, which will return a tuple of types (like `(Int64, Float64, String)`).

The actual type definition is
```julia
struct Schema{names, types} end
```
Where `names` is a tuple of Symbols, and `types` is a tuple _type_ of types (like `Tuple{Int64, Float64, String}`).
Encoding the names & types as type parameters allows convenient use of the type in generated functions
and other optimization use-cases.
"""
struct Schema{names, types} end
Schema(names::Tuple{Vararg{Symbol}}, types::Type{T}) where {T <: Tuple} = Schema{names, T}()
Schema(::Type{NamedTuple{names, types}}) where {names, types} = Schema{names, types}()

# pass through Ints to allow Tuples to act as rows
sym(x) = Symbol(x)
sym(x::Int) = x

Schema(names, ::Nothing) = Schema{Tuple(Base.map(sym, names)), nothing}()
Schema(names, types) = Schema{Tuple(Base.map(sym, names)), Tuple{types...}}()

function Base.show(io::IO, sch::Schema{names, types}) where {names, types}
    println(io, "Tables.Schema:")
    Base.print_matrix(io, hcat(collect(names), types === nothing ? fill(nothing, length(names)) : collect(fieldtype(types, i) for i = 1:fieldcount(types))))
end

function Base.getproperty(sch::Schema{names, types}, field::Symbol) where {names, types}
    if field === :names
        return names
    elseif field === :types
        return types === nothing ? nothing : Tuple(fieldtype(types, i) for i = 1:fieldcount(types))
    else
        throw(ArgumentError("unsupported property for Tables.Schema"))
    end
end

Base.propertynames(sch::Schema) = (:names, :types)

# helper functions
include("utils.jl")

# reference implementations: Vector of NamedTuples and NamedTuple of Vectors
include("namedtuples.jl")

# generic fallback definitions
include("fallbacks.jl")

# allow any valid iterator to be a table
include("tofromdatavalues.jl")

# simple table operations on table inputs
include("operations.jl")

# matrix integration
include("matrix.jl")

"Return the column index (1-based) of a `colname` in a table with a known schema; returns 0 if `colname` doesn't exist in table"
columnindex(table, colname) = columnindex(schema(table).names, colname)

"Return the column type of a `colname` in a table with a known schema; returns Union{} if `colname` doesn't exist in table"
columntype(table, colname) = columntype(schema(table), colname)

end # module


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

# Expr head => argument count bounds
const VALID_EXPR_HEADS = IdDict{Any,Any}(
    :call => 1:typemax(Int),
    :invoke => 2:typemax(Int),
    :static_parameter => 1:1,
    :gotoifnot => 2:2,
    :(&) => 1:1,
    :(=) => 2:2,
    :method => 1:4,
    :const => 1:1,
    :new => 1:typemax(Int),
    :splatnew => 2:2,
    :return => 1:1,
    :unreachable => 0:0,
    :the_exception => 0:0,
    :enter => 1:1,
    :leave => 1:1,
    :pop_exception => 1:1,
    :inbounds => 1:1,
    :boundscheck => 0:0,
    :copyast => 1:1,
    :meta => 0:typemax(Int),
    :global => 1:1,
    :foreigncall => 5:typemax(Int), # name, RT, AT, nreq, cconv, args..., roots...
    :cfunction => 5:5,
    :isdefined => 1:1,
    :code_coverage_effect => 0:0,
    :loopinfo => 0:typemax(Int),
    :gc_preserve_begin => 0:typemax(Int),
    :gc_preserve_end => 0:typemax(Int),
    :thunk => 1:1,
    :throw_undef_if_not => 2:2
)

# @enum isn't defined yet, otherwise I'd use it for this
const INVALID_EXPR_HEAD = "invalid expression head"
const INVALID_EXPR_NARGS = "invalid number of expression args"
const INVALID_LVALUE = "invalid LHS value"
const INVALID_RVALUE = "invalid RHS value"
const INVALID_RETURN = "invalid argument to :return"
const INVALID_CALL_ARG = "invalid :call argument"
const EMPTY_SLOTNAMES = "slotnames field is empty"
const SLOTFLAGS_MISMATCH = "length(slotnames) < length(slotflags)"
const SSAVALUETYPES_MISMATCH = "not all SSAValues in AST have a type in ssavaluetypes"
const SSAVALUETYPES_MISMATCH_UNINFERRED = "uninferred CodeInfo ssavaluetypes field does not equal the number of present SSAValues"
const NON_TOP_LEVEL_METHOD = "encountered `Expr` head `:method` in non-top-level code (i.e. `nargs` > 0)"
const NON_TOP_LEVEL_GLOBAL = "encountered `Expr` head `:global` in non-top-level code (i.e. `nargs` > 0)"
const SIGNATURE_NARGS_MISMATCH = "method signature does not match number of method arguments"
const SLOTNAMES_NARGS_MISMATCH = "CodeInfo for method contains fewer slotnames than the number of method arguments"

struct InvalidCodeError <: Exception
    kind::AbstractString
    meta::Any
end
InvalidCodeError(kind::AbstractString) = InvalidCodeError(kind, nothing)

function validate_code_in_debug_mode(linfo::MethodInstance, src::CodeInfo, kind::String)
    if JLOptions().debug_level == 2
        # this is a debug build of julia, so let's validate linfo
        errors = validate_code(linfo, src)
        if !isempty(errors)
            for e in errors
                if linfo.def isa Method
                    println(stderr, "WARNING: Encountered invalid ", kind, " code for method ",
                            linfo.def, ": ", e)
                else
                    println(stderr, "WARNING: Encountered invalid ", kind, " code for top level expression in ",
                            linfo.def, ": ", e)
                end
            end
        end
    end
end

"""
    validate_code!(errors::Vector{>:InvalidCodeError}, c::CodeInfo)

Validate `c`, logging any violation by pushing an `InvalidCodeError` into `errors`.
"""
function validate_code!(errors::Vector{>:InvalidCodeError}, c::CodeInfo, is_top_level::Bool = false)
    function validate_val!(@nospecialize(x))
        if isa(x, Expr)
            if x.head === :call || x.head === :invoke
                f = x.args[1]
                if f isa GlobalRef && (f.name === :cglobal) && x.head === :call
                    # TODO: these are not yet linearized
                else
                    for arg in x.args
                        if !is_valid_argument(arg)
                            push!(errors, InvalidCodeError(INVALID_CALL_ARG, arg))
                        else
                            validate_val!(arg)
                        end
                    end
                end
            end
        elseif isa(x, SSAValue)
            id = x.id
            !in(id, ssavals) && push!(ssavals, id)
        end
    end

    ssavals = BitSet()
    lhs_slotnums = BitSet()
    for x in c.code
        if isa(x, Expr)
            head = x.head
            if !is_top_level
                head === :method && push!(errors, InvalidCodeError(NON_TOP_LEVEL_METHOD))
                head === :global && push!(errors, InvalidCodeError(NON_TOP_LEVEL_GLOBAL))
            end
            narg_bounds = get(VALID_EXPR_HEADS, head, -1:-1)
            nargs = length(x.args)
            if narg_bounds == -1:-1
                push!(errors, InvalidCodeError(INVALID_EXPR_HEAD, (head, x)))
            elseif !in(nargs, narg_bounds)
                push!(errors, InvalidCodeError(INVALID_EXPR_NARGS, (head, nargs, x)))
            elseif head === :(=)
                lhs, rhs = x.args
                if !is_valid_lvalue(lhs)
                    push!(errors, InvalidCodeError(INVALID_LVALUE, lhs))
                elseif isa(lhs, SlotNumber) && !in(lhs.id, lhs_slotnums)
                    n = lhs.id
                    push!(lhs_slotnums, n)
                end
                if !is_valid_rvalue(rhs)
                    push!(errors, InvalidCodeError(INVALID_RVALUE, rhs))
                end
                validate_val!(lhs)
                validate_val!(rhs)
            elseif head === :gotoifnot
                if !is_valid_argument(x.args[1])
                    push!(errors, InvalidCodeError(INVALID_CALL_ARG, x.args[1]))
                end
                validate_val!(x.args[1])
            elseif head === :return
                if !is_valid_return(x.args[1])
                    push!(errors, InvalidCodeError(INVALID_RETURN, x.args[1]))
                end
                validate_val!(x.args[1])
            elseif head === :call || head === :invoke || head === :gc_preserve_end || head === :meta ||
                head === :inbounds || head === :foreigncall || head === :cfunction ||
                head === :const || head === :enter || head === :leave || head === :pop_exception ||
                head === :method || head === :global || head === :static_parameter ||
                head === :new || head === :splatnew || head === :thunk || head === :loopinfo ||
                head === :throw_undef_if_not || head === :unreachable || head === :code_coverage_effect
                validate_val!(x)
            else
                # TODO: nothing is actually in statement position anymore
                #push!(errors, InvalidCodeError("invalid statement", x))
            end
        elseif isa(x, NewvarNode)
        elseif isa(x, GotoNode)
        elseif x === nothing
        elseif isa(x, SlotNumber)
        elseif isa(x, GlobalRef)
        elseif isa(x, LineNumberNode)
        elseif isa(x, PiNode)
        elseif isa(x, PhiCNode)
        elseif isa(x, PhiNode)
        elseif isa(x, UpsilonNode)
        else
            #push!(errors, InvalidCodeError("invalid statement", x))
        end
    end
    nslotnames = length(c.slotnames)
    nslotflags = length(c.slotflags)
    nssavals = length(c.code)
    !is_top_level && nslotnames == 0 && push!(errors, InvalidCodeError(EMPTY_SLOTNAMES))
    nslotnames < nslotflags && push!(errors, InvalidCodeError(SLOTFLAGS_MISMATCH, (nslotnames, nslotflags)))
    if c.inferred
        nssavaluetypes = length(c.ssavaluetypes)
        nssavaluetypes < nssavals && push!(errors, InvalidCodeError(SSAVALUETYPES_MISMATCH, (nssavals, nssavaluetypes)))
    else
        c.ssavaluetypes != nssavals && push!(errors, InvalidCodeError(SSAVALUETYPES_MISMATCH_UNINFERRED, (nssavals, c.ssavaluetypes)))
    end
    return errors
end

"""
    validate_code!(errors::Vector{>:InvalidCodeError}, mi::MethodInstance,
                   c::Union{Nothing,CodeInfo} = Core.Compiler.retrieve_code_info(mi))

Validate `mi`, logging any violation by pushing an `InvalidCodeError` into `errors`.

If `isa(c, CodeInfo)`, also call `validate_code!(errors, c)`. It is assumed that `c` is
the `CodeInfo` instance associated with `mi`.
"""
function validate_code!(errors::Vector{>:InvalidCodeError}, mi::Core.MethodInstance,
                        c::Union{Nothing,CodeInfo} = Core.Compiler.retrieve_code_info(mi))
    is_top_level = mi.def isa Module
    if is_top_level
        mnargs = 0
    else
        m = mi.def::Method
        mnargs = m.nargs
        n_sig_params = length(Core.Compiler.unwrap_unionall(m.sig).parameters)
        if (m.isva ? (n_sig_params < (mnargs - 1)) : (n_sig_params != mnargs))
            push!(errors, InvalidCodeError(SIGNATURE_NARGS_MISMATCH, (m.isva, n_sig_params, mnargs)))
        end
    end
    if isa(c, CodeInfo)
        mnargs > length(c.slotnames) && push!(errors, InvalidCodeError(SLOTNAMES_NARGS_MISMATCH))
        validate_code!(errors, c, is_top_level)
    end
    return errors
end

validate_code(args...) = validate_code!(Vector{InvalidCodeError}(), args...)

is_valid_lvalue(@nospecialize(x)) = isa(x, Slot) || isa(x, GlobalRef)

function is_valid_argument(@nospecialize(x))
    if isa(x, Slot) || isa(x, SSAValue) || isa(x, GlobalRef) || isa(x, QuoteNode) ||
        (isa(x,Expr) && (x.head in (:static_parameter, :boundscheck))) ||
        isa(x, Number) || isa(x, AbstractString) || isa(x, AbstractChar) || isa(x, Tuple) ||
        isa(x, Type) || isa(x, Core.Box) || isa(x, Module) || x === nothing
        return true
    end
    # TODO: consider being stricter about what needs to be wrapped with QuoteNode
    return !(isa(x,Expr) || isa(x,Symbol) || isa(x,GotoNode) ||
             isa(x,LineNumberNode) || isa(x,NewvarNode))
end

function is_valid_rvalue(@nospecialize(x))
    is_valid_argument(x) && return true
    if isa(x, Expr) && x.head in (:new, :splatnew, :the_exception, :isdefined, :call, :invoke, :foreigncall, :cfunction, :gc_preserve_begin, :copyast)
        return true
    end
    return false
end

is_valid_return(@nospecialize(x)) = is_valid_argument(x) || (isa(x, Expr) && x.head === :lambda)

is_flag_set(byte::UInt8, flag::UInt8) = (byte & flag) == flag

ext: jl

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
using WignerFamilies
using PyPlot
using WignerSymbols

j₂, j₃, m₂, m₃ = 10, 10, 0, 0
w = WignerF(Float64, j₂, j₃, m₂, m₃)
println((w.nₘᵢₙ, w.nₘₐₓ))
w3j = get_wigner_array(w)
WignerFamilies.wigner3j_f!(w, w3j)
j_array = collect(eachindex(w3j))
w3j

##

plt.clf()
plot(j_array, Float64.(w3j.symbols))
plt.plot(j_array,
    [Float64(wigner3j(BigFloat, n, j₂, j₃, 0, 0)) for n in j_array], alpha=0.5)
plt.gcf()

##

for j₂ in 1:100
    j₃, m₂, m₃ = 10, 0, 0
    w3j = wigner3j_f(Float64, j₂, j₃, m₂, m₃)
    j_array = collect(eachindex(w3j))
    for j in j_array
        @assert w3j[j] ≈ Float64(WignerSymbols.wigner3j(BigFloat, j, j₂, j₃, -m₂-m₃, m₂))
    end
end


text:
using Documenter, Wavefronts

makedocs(
    sitename="Wavefronts.jl documentation",
    pages = [
        "index.md",
        "example.md"
    ],
    modules = [Wavefronts]
)

deploydocs(
    repo = "github.com/Klafyvel/Wavefronts.jl.git",
)


text:
module GridFit

using SparseArrays
using Statistics
using StatsBase
using LinearAlgebra

import C2PO: nan, meshgrid, histc
export gridfit

# define parameter type
mutable struct Parameters
    smoothness::Array{Float64}
    interp::String
    regularizer::String
    solver::String
    maxiter::Int64
    tilesize::Float64
    overlap::Float64
    mask::AbstractArray
    autoscale::String
    xscale::Float64
    yscale::Float64
    maskflag::Bool
    extend::String
end

function gridfit(x::Array{Float64}, y::Array{Float64}, z::Array{Float64}, xnodes::Array{Float64}, ynodes::Array{Float64}; smoothness::Array{Float64}=[1.0 1.0], interp::String="triangle", regularizer::String="gradient", solver::String="normal", maxiter::Int64=10000, extend::String="warning", tilesize::Float64=Inf, overlap::Float64=0.2, autoscale::String="on")

    # gridfit: estimates a surface on a 2d grid, based on scattered data
    #          Replicates are allowed. All methods extrapolate to the grid
    #          boundaries. Gridfit uses a modified ridge estimator to
    #          generate the surface, where the bias is toward smoothness.
    #
    #          Gridfit is not an interpolant. Its goal is a smooth surface
    #          that approximates your data, but allows you to control the
    #          amount of smoothing.
    #
    # usage #1: zgrid = gridfit(x,y,z,xnodes,ynodes);
    # usage #2: [zgrid,xgrid,ygrid] = gridfit(x,y,z,xnodes,ynodes);
    # usage #3: zgrid = gridfit(x,y,z,xnodes,ynodes,prop,val,prop,val,...);
    #
    # Arguments: (input)
    #  x,y,z - vectors of equal lengths, containing arbitrary scattered data
    #          The only constraint on x and y is they cannot ALL fall on a
    #          single line in the x-y plane. Replicate points will be treated
    #          in a least squares sense.
    #
    #          ANY points containing a NaN are ignored in the estimation
    #
    #  xnodes - vector defining the nodes in the grid in the independent
    #          variable (x). xnodes need not be equally spaced. xnodes
    #          must completely span the data. If they do not, then the
    #          'extend' property is applied, adjusting the first and last
    #          nodes to be extended as necessary. See below for a complete
    #          description of the 'extend' property.
    #
    #          If xnodes is a scalar integer, then it specifies the number
    #          of equally spaced nodes between the min and max of the data.
    #
    #  ynodes - vector defining the nodes in the grid in the independent
    #          variable (y). ynodes need not be equally spaced.
    #
    #          If ynodes is a scalar integer, then it specifies the number
    #          of equally spaced nodes between the min and max of the data.
    #
    #          Also see the extend property.
    #
    #  Additional arguments follow in the form of property/value pairs.
    #  Valid properties are:
    #    'smoothness', 'interp', 'regularizer', 'solver', 'maxiter'
    #    'extend', 'tilesize', 'overlap'
    #
    #  Any UNAMBIGUOUS shortening (even down to a single letter) is
    #  valid for property names. All properties have default values,
    #  chosen (I hope) to give a reasonable result out of the box.
    #
    #   'smoothness' - scalar or vector of length 2 - determines the
    #          eventual smoothness of the estimated surface. A larger
    #          value here means the surface will be smoother. Smoothness
    #          must be a non-negative real number.
    #
    #          If this parameter is a vector of length 2, then it defines
    #          the relative smoothing to be associated with the x and y
    #          variables. This allows the user to apply a different amount
    #          of smoothing in the x dimension compared to the y dimension.
    #
    #          Note: the problem is normalized in advance so that a
    #          smoothness of 1 MAY generate reasonable results. If you
    #          find the result is too smooth, then use a smaller value
    #          for this parameter. Likewise, bumpy surfaces suggest use
    #          of a larger value. (Sometimes, use of an iterative solver
    #          with too small a limit on the maximum number of iterations
    #          will result in non-convergence.)
    #
    #          DEFAULT: 1
    #
    #
    #   'interp' - character, denotes the interpolation scheme used
    #          to interpolate the data.
    #
    #          DEFAULT: 'triangle'
    #
    #          'bilinear' - use bilinear interpolation within the grid
    #                     (also known as tensor product linear interpolation)
    #
    #          'triangle' - split each cell in the grid into a triangle,
    #                     then linear interpolation inside each triangle
    #
    #          'nearest' - nearest neighbor interpolation. This will
    #                     rarely be a good choice, but I included it
    #                     as an option for completeness.
    #
    #
    #   'regularizer' - character flag, denotes the regularization
    #          paradignm to be used. There are currently three options.
    #
    #          DEFAULT: 'gradient'
    #
    #          'diffusion' or 'laplacian' - uses a finite difference
    #              approximation to the Laplacian operator (i.e, del^2).
    #
    #              We can think of the surface as a plate, wherein the
    #              bending rigidity of the plate is specified by the user
    #              as a number relative to the importance of fidelity to
    #              the data. A stiffer plate will result in a smoother
    #              surface overall, but fit the data less well. I've
    #              modeled a simple plate using the Laplacian, del^2. (A
    #              projected enhancement is to do a better job with the
    #              plate equations.)
    #
    #              We can also view the regularizer as a diffusion problem,
    #              where the relative thermal conductivity is supplied.
    #              Here interpolation is seen as a problem of finding the
    #              steady temperature profile in an object, given a set of
    #              points held at a fixed temperature. Extrapolation will
    #              be linear. Both paradigms are appropriate for a Laplacian
    #              regularizer.
    #
    #          'gradient' - attempts to ensure the gradient is as smooth
    #              as possible everywhere. Its subtly different from the
    #              'diffusion' option, in that here the directional
    #              derivatives are biased to be smooth across cell
    #              boundaries in the grid.
    #
    #              The gradient option uncouples the terms in the Laplacian.
    #              Think of it as two coupled PDEs instead of one PDE. Why
    #              are they different at all? The terms in the Laplacian
    #              can balance each other.
    #
    #          'springs' - uses a spring model connecting nodes to each
    #              other, as well as connecting data points to the nodes
    #              in the grid. This choice will cause any extrapolation
    #              to be as constant as possible.
    #
    #              Here the smoothing parameter is the relative stiffness
    #              of the springs connecting the nodes to each other compared
    #              to the stiffness of a spting connecting the lattice to
    #              each data point. Since all springs have a rest length
    #              (length at which the spring has zero potential energy)
    #              of zero, any extrapolation will be minimized.
    #
    #          Note: The 'springs' regularizer tends to drag the surface
    #          towards the mean of all the data, so too large a smoothing
    #          parameter may be a problem.
    #
    #
    #   'solver' - character flag - denotes the solver used for the
    #          resulting linear system. Different solvers will have
    #          different solution times depending upon the specific
    #          problem to be solved. Up to a certain size grid, the
    #          direct \ solver will often be speedy, until memory
    #          swaps causes problems.
    #
    #          What solver should you use? Problems with a significant
    #          amount of extrapolation should avoid lsqr. \ may be
    #          best numerically for small smoothnesss parameters and
    #          high extents of extrapolation.
    #
    #          Large numbers of points will slow down the direct
    #          \, but when applied to the normal equations, \ can be
    #          quite fast. Since the equations generated by these
    #          methods will tend to be well conditioned, the normal
    #          equations are not a bad choice of method to use. Beware
    #          when a small smoothing parameter is used, since this will
    #          make the equations less well conditioned.
    #
    #          DEFAULT: 'normal'
    #
    #          '\' - uses matlab's backslash operator to solve the sparse
    #                     system. 'backslash' is an alternate name.
    #
    #          'symmlq' - uses matlab's iterative symmlq solver
    #
    #          'lsqr' - uses matlab's iterative lsqr solver
    #
    #          'normal' - uses \ to solve the normal equations.
    #
    #
    #   'maxiter' - only applies to iterative solvers - defines the
    #          maximum number of iterations for an iterative solver
    #
    #          DEFAULT: min(10000,length(xnodes)*length(ynodes))
    #
    #
    #   'extend' - character flag - controls whether the first and last
    #          nodes in each dimension are allowed to be adjusted to
    #          bound the data, and whether the user will be warned if
    #          this was deemed necessary to happen.
    #
    #          DEFAULT: 'warning'
    #
    #          'warning' - Adjust the first and/or last node in
    #                     x or y if the nodes do not FULLY contain
    #                     the data. Issue a warning message to this
    #                     effect, telling the amount of adjustment
    #                     applied.
    #
    #          'never'  - Issue an error message when the nodes do
    #                     not absolutely contain the data.
    #
    #          'always' - automatically adjust the first and last
    #                     nodes in each dimension if necessary.
    #                     No warning is given when this option is set.
    #
    #
    #   'tilesize' - grids which are simply too large to solve for
    #          in one single estimation step can be built as a set
    #          of tiles. For example, a 1000x1000 grid will require
    #          the estimation of 1e6 unknowns. This is likely to
    #          require more memory (and time) than you have available.
    #          But if your data is dense enough, then you can model
    #          it locally using smaller tiles of the grid.
    #
    #          My recommendation for a reasonable tilesize is
    #          roughly 100 to 200. Tiles of this size take only
    #          a few seconds to solve normally, so the entire grid
    #          can be modeled in a finite amount of time. The minimum
    #          tilesize can never be less than 3, although even this
    #          size tile is so small as to be ridiculous.
    #
    #          If your data is so sparse than some tiles contain
    #          insufficient data to model, then those tiles will
    #          be left as NaNs.
    #
    #          DEFAULT: inf
    #
    #
    #   'overlap' - Tiles in a grid have some overlap, so they
    #          can minimize any problems along the edge of a tile.
    #          In this overlapped region, the grid is built using a
    #          bi-linear combination of the overlapping tiles.
    #
    #          The overlap is specified as a fraction of the tile
    #          size, so an overlap of 0.20 means there will be a 20%
    #          overlap of successive tiles. I do allow a zero overlap,
    #          but it must be no more than 1/2.
    #
    #          0 <= overlap <= 0.5
    #
    #          Overlap is ignored if the tilesize is greater than the
    #          number of nodes in both directions.
    #
    #          DEFAULT: 0.20
    #
    #
    #   'autoscale' - Some data may have widely different scales on
    #          the respective x and y axes. If this happens, then
    #          the regularization may experience difficulties.
    #
    #          autoscale = 'on' will cause gridfit to scale the x
    #          and y node intervals to a unit length. This should
    #          improve the regularization procedure. The scaling is
    #          purely internal.
    #
    #          autoscale = 'off' will disable automatic scaling
    #
    #          DEFAULT: 'on'
    #
    #
    # Arguments: (output)
    #  zgrid   - (nx,ny) array containing the fitted surface
    #
    #  xgrid, ygrid - as returned by meshgrid(xnodes,ynodes)
    #
    #
    # Speed considerations:
    #  Remember that gridfit must solve a LARGE system of linear
    #  equations. There will be as many unknowns as the total
    #  number of nodes in the final lattice. While these equations
    #  may be sparse, solving a system of 10000 equations may take
    #  a second or so. Very large problems may benefit from the
    #  iterative solvers or from tiling.
    #
    #
    # Example usage:
    #
    #  x = rand(100,1);
    #  y = rand(100,1);
    #  z = exp(x+2*y);
    #  xnodes = 0:.1:1;
    #  ynodes = 0:.1:1;
    #
    #  g = gridfit(x,y,z,xnodes,ynodes);
    #
    # Note: this is equivalent to the following call:
    #
    #  g = gridfit(x,y,z,xnodes,ynodes, ...
    #              'smooth',1, ...
    #              'interp','triangle', ...
    #              'solver','normal', ...
    #              'regularizer','gradient', ...
    #              'extend','warning', ...
    #              'tilesize',inf);
    #
    #
    # Author: John D'Errico
    # e-mail address: woodchips@rochester.rr.com
    # Release: 2.0
    # Release date: 5/23/06

    testflag = 0
    if testflag == 1
        # create test dataset
        #x = [collect(1:0.5:100.0); 200.0];
        #y = deepcopy(x);
        #z = deepcopy(x);
        #x[10] = NaN;
        #y[30] = NaN;
        #z[70] = NaN;
        #xnodes = collect(1:100.0);
        #ynodes = collect(1:150.0);

        n1 = 15;
        n2 = 20;
        theta = rand(n1,1) .* pi/2;
        r = rand(1,n2);

        x = cos.(theta) * r;
        y = sin.(theta) * r;
        x=x[:];
        y=y[:];

        x = [Array{Float64,1}([0;0;1;1]); x; x; 1 .- x; 1 .- x];
        y = [Array{Float64,1}([0;1;0;1]); y; 1 .- y; y; 1 .- y];

        z = sin.(4 .* x .+ 5 .* y) .* cos.(7 .* (x .- y)) .+ exp.(x .+ y);

        xi = collect(0.0 : 1.0/200.0 : 1.0);
        yi = collect(0.0 : 1.0/100.0 : 1.0);
    end

    # ensure all of x,y,z,xnodes,ynodes are column vectors, also drop any NaN data
    x = x[:]; y = y[:]; z = z[:];
    k = findall((.!isnan.(x) .& .!isnan.(y) .& .!isnan.(z)) .== true);
    x = x[k]; y = y[k]; z = z[k];
    xmin = minimum(x);
    xmax = maximum(x);
    ymin = minimum(y);
    ymax = maximum(y);

    # did they supply a scalar for the nodes?
    if length(xnodes)==1
        xnodes = collect(xmin : (xmax-xmin)/xnodes : xmax);
        xnodes[end] = xmax; # make sure it hits the max
    end
    if length(ynodes)==1
        ynodes = collect(ymin : (ymax-ymin)/ynodes : ymax);
        ynodes[end] = ymax; # make sure it hits the max
    end

    xnodes=xnodes[:];
    ynodes=ynodes[:];
    dx = diff(xnodes);
    dy = diff(ynodes);
    nx = length(xnodes);
    ny = length(ynodes);
    ngrid = nx*ny;

    # set default parameters
    #params = Parameters([1.0 1.5], "triangle", "gradient", "backslash", 10000, Inf, 0.2, Array{Bool,2}(undef,nx,ny), "on", 1.0, 1.0, Bool(0), "always");
    #params.mask .= 1;
    params = Parameters([1.0 1.5], "triangle", "gradient", "backslash", 10000, Inf, 0.2, [], "on", 1.0, 1.0, Bool(0), "always");

    # set the scaling if autoscale was on
    if params.autoscale == "on"
        params.xscale = mean(dx);
        params.yscale = mean(dy);
        params.autoscale = "off";
    end

    if params.tilesize < max(nx,ny)
        # split it into smaller tiles. compute zgrid and ygrid
        # at the very end if requested
        #zgrid = tiled_gridfit(x,y,z,xnodes,ynodes,params);
        zgrid = [];
        return zgrid;
    else
        #it's a single tile.

        # mask must be either an empty array, or a boolean array of the same size as the final grid.
        nmask = size(params.mask);
        if !isempty(params.mask)
            if ((nmask[2] != nx) | (nmask[1] != ny))
                if ((nmask[2]==ny) | (nmask[1]==nx))
                    ErrorException("Mask array is probably transposed from proper orientation.")
                else
                    ErrorException("Mask array must be the same size as the final grid.")
                end
            end
        end
        if !isempty(params.mask)
            params.maskflag = 1;
        else
            params.maskflag = 0;
        end

        # default for maxiter?
        if isempty(params.maxiter)
            params.maxiter = min(10000,nx*ny);
        end

        # check lengths of the data
        n = length(x);
        if (length(y) != n) | (length(z) != n)
            ErrorException("Data vectors are incompatible in size.")
        end
        if n < 3
            ErrorException("Insufficient data for surface estimation.")
        end

        # verify the nodes are distinct
        if any(diff(xnodes) .<= 0) .| any(diff(ynodes) .<= 0)
            ErrorException("xnodes and ynodes must be monotone increasing")
        end

        # do we need to tweak the first or last node in x or y?
        if xmin<xnodes[1]
            if params.extend == "always"
                xnodes[1] = xmin;
            elseif param.extend == "warning"
                display("WARNING: GRIDFIT:extend xnodes[1] was decreased by: " * string(xnodes[1]-xmin) * ", new node = " * string(xmin))
                xnodes[1] = xmin;
            elseif param.extend == "never"
                ErrorException("Some x (" * string(xmin) * ") falls below xnodes[1] by: " * string(xnodes[1]-xmin))
            end
        end
        if xmax>xnodes[end]
            if params.extend == "always"
                xnodes[end] = xmax;
            elseif params.extend == "warning"
                display("WARNING: GRIDFIT:extend xnodes(end) was increased by: " * string(xmax-xnodes[end]) * ", new node = " * string(xmax))
                xnodes[end] = xmax;
            elseif params.extend == "never"
                ErrorException("Some x (" * string(xmax) * ") falls above xnodes[end] by: " * string(xmax-xnodes[end]))
            end
        end
        if ymin<ynodes[1]
            if params.extend == "always"
                ynodes[1] = ymin;
            elseif params.extend == "warning"
                display("WARNING: GRIDFIT:extend ynodes[1] was decreased by: " * string(ynodes[1]-ymin) * ", new node = " * string(ymin))
                ynodes[1] = ymin;
            elseif params.extend == "never"
                ErrorException("Some y (" * string(ymin) * ") falls below ynodes[1] by: " * string(ynodes[1]-ymin))
            end
        end
        if ymax>ynodes[end]
            if params.extend == "always"
                ynodes[end] = ymax;
            elseif params.extend == "warning"
                display("WARNING: GRIDFIT:extend ynodes[end] was increased by: " * string(ymax-ynodes[end]) * ", new node = " * string(ymax))
                ynodes[end] = ymax;
            elseif params.extend == "never"
                ErrorException("Some y (" * string(ymax) * ") falls above ynodes[end] by: " * string(ymax-ynodes[end]))
            end
        end

        # determine which cell in the array each point lies in
        #hx = fit(Histogram, x, xnodes);
        #hy = fit(Histogram, y, ynodes);
        indx = histc(x, xnodes);
        indy = histc(y, ynodes);

        # any point falling at the last node is taken to be inside the last cell in x or y.
        k = findall(indx .== nx);
        indx[k]=indx[k] .- 1;
        k = findall(indy .== ny);
        indy[k]=indy[k] .- 1;
        ind = indy .+ ny .* (indx .- 1);

        # Do we have a mask to apply?
        if params.maskflag
            # if we do, then we need to ensure that every
            # cell with at least one data point also has at
            # least all of its corners unmasked.
            params.mask[ind] .= 1;
            params.mask[ind .+ 1] .= 1;
            params.mask[ind .+ ny] .= 1;
            params.mask[ind .+ ny .+ 1] .= 1;
        end

        # interpolation equations for each point
        tx = min.(1.0,max.(0.0,(x .- xnodes[indx])./dx[indx]));
        ty = min.(1.0,max.(0.0,(y .- ynodes[indy])./dy[indy]));

        if params.interp == "triangle"
            # linear interpolation inside each triangle
            k = findall(tx .> ty);
            L = ones(n);
            L[k] .= ny;

            t1 = min(tx,ty);
            t2 = max(tx,ty);
            #A = sparse(repeat(collect(1:n),[1,3]...), [ind ind .+ ny .+ 1 ind .+ L],[1 .- t2 t1 t2 .- t1], n, ngrid);
            A = sparse(repeat(collect(1:n),[3]...), [ind; ind .+ ny .+ 1; ind .+ L],[1 .- t2; t1; t2 .- t1], n, ngrid);

        elseif params.interp == "nearest"
            # nearest neighbor interpolation in a cell
            k = round.(1 .- ty) + round.(1 .- tx) .* ny;
            A = sparse(collect(1:n), ind .+ k, ones(n), n, ngrid);

        elseif params.interp == "bilinear"
            # bilinear interpolation in a cell
            A = sparse(repeat(collect(1:n),[4]...), [ind; ind .+ 1; ind .+ ny .+ 1; ind .+ ny .+ 1], [(1 .- tx) .* (1 .- ty); (1 .- tx) .* ty; tx .* (1 .- ty); tx .* ty], n, ngrid);
            #A = sparse(repmat((1:n)',1,4),[ind,ind+1,ind+ny,ind+ny+1],[(1-tx).*(1-ty), (1-tx).*ty, tx.*(1-ty), tx.*ty], n,ngrid);
        end
        rhs = z;

        # do we have relative smoothing parameters?
        if length(params.smoothness) == 1
            # it was scalar, so treat both dimensions equally
            smoothparam = params.smoothness;
            xyRelativeStiffness = [1.0; 1.0];
        else
            # It was a vector, so anisotropy reigns.
            # I've already checked that the vector was of length 2
            smoothparam = sqrt(prod(params.smoothness));
            xyRelativeStiffness = params.smoothness[:] ./ smoothparam;
        end

        # Build regularizer. Add del^4 regularizer one day.
        if params.regularizer == "springs"
            # zero "rest length" springs
            #i = [i for i in 1:nx, j in 1:(ny-1)];
            #j = [j for i in 1:nx, j in 1:(ny-1)];
            (i,j) = meshgrid(1:nx,1:(ny-1));
            ind = j[:] .+ ny .* (i[:] .- 1);
            m = nx*(ny-1);
            stiffness = 1 ./ (dy ./ params.yscale);
            Areg = sparse(repeat(collect(1:m),[2]...), [ind; ind .+ 1], (xyRelativeStiffness[2] * stiffness[j[:]] .* [-1 1])[:], m, ngrid);
            #Areg = sparse(repmat((1:m)',1,2),[ind,ind+1], xyRelativeStiffness(2)*stiffness(j(:))*[-1 1], m,ngrid); #matlab version

            #i = [i for i in 1:(nx-1), j in 1:ny];
            #j = [j for i in 1:(nx-1), j in 1:ny];
            (i,j) = meshgrid(1:(nx-1),1:ny);
            ind = j[:] .+ ny .* (i[:] .- 1);
            m = (nx-1)*ny;
            stiffness = 1 ./ (dx ./ params.xscale);
            Areg = [Areg; sparse(repeat(collect(1:m),[2]...), [ind; ind .+ ny], (xyRelativeStiffness[1] .* stiffness[i[:]] .* [-1 1])[:], m, ngrid)];
            #Areg = [Areg;sparse(repmat((1:m)',1,2),[ind,ind+ny], xyRelativeStiffness(1)*stiffness(i(:))*[-1 1],m,ngrid)]; #matlab version

            #i = [i for i in 1:(nx-1), j in 1:(ny-1)];
            #j = [j for i in 1:(nx-1), j in 1:(ny-1)];
            (i,j) = meshgrid(1:(nx-1),1:(ny-1));
            ind = j[:] .+ ny .* (i[:] .- 1);
            m = (nx-1)*(ny-1);
            stiffness = 1 ./ sqrt.((dx[i[:]] ./ params.xscale ./ xyRelativeStiffness[1]) .^ 2 + (dy[j[:]] ./ params.yscale ./ xyRelativeStiffness[2]) .^ 2);
            Areg = [Areg; sparse(repeat(collect(1:m),[2]...), [ind; ind .+ ny .+ 1], (stiffness .* [-1 1])[:], m, ngrid)];
            #Areg = [Areg;sparse(repmat((1:m)',1,2),[ind,ind+ny+1], stiffness*[-1 1],m,ngrid)];
            Areg = [Areg; sparse(repeat(collect(1:m),[2]...), [ind .+ 1; ind .+ ny], (stiffness .* [-1 1])[:], m, ngrid)];
            #Areg = [Areg;sparse(repmat((1:m)',1,2),[ind+1,ind+ny], stiffness*[-1 1],m,ngrid)];

        elseif (params.regularizer == "diffusion") | (params.regularizer == "laplacian")
            # thermal diffusion using Laplacian (del^2)
            #i = [i for i in 1:nx, j in 2:(ny-1)];
            #j = [j for i in 1:nx, j in 2:(ny-1)];
            (i,j) = meshgrid(1:nx,2:(ny-1));
            ind = j[:] .+ ny .* (i[:] .- 1);
            dy1 = dy[j[:] .- 1] ./ params.yscale;
            dy2 = dy[j[:]] ./ params.yscale;
            Areg = sparse(repeat(ind,[3]...), [ind .- 1; ind; ind .+ 1], xyRelativeStiffness[2] .* [-2 ./ (dy1 .* (dy1 .+ dy2)); 2 ./ (dy1 .* dy2); -2 ./ (dy2 .* (dy1 .+ dy2))], ngrid, ngrid);
            #Areg = sparse(repmat(ind,1,3),[ind-1,ind,ind+1], xyRelativeStiffness(2)*[-2./(dy1.*(dy1+dy2)), 2./(dy1.*dy2), -2./(dy2.*(dy1+dy2))],ngrid,ngrid);

            #i = [i for i in 2:(nx-1), j in 1:ny];
            #j = [j for i in 2:(nx-1), j in 1:ny];
            (i,j) = meshgrid(2:(nx-1),1:ny);
            ind = j[:] .+ ny .* (i[:] .- 1);
            dx1 = dx[i[:] .- 1] ./ params.xscale;
            dx2 = dx[i[:]] ./ params.xscale;
            Areg .= Areg .+ sparse(repeat(ind,[3]...), [ind .- ny; ind; ind .+ ny], xyRelativeStiffness[1] .* [-2 ./ (dx1 .* (dx1 .+ dx2)); 2 ./ (dx1 .* dx2); -2 ./ (dx2 .* (dx1 .+ dx2))], ngrid, ngrid);
            #Areg = Areg + sparse(repmat(ind,1,3),[ind-ny,ind,ind+ny], xyRelativeStiffness(1)*[-2./(dx1.*(dx1+dx2)), 2./(dx1.*dx2), -2./(dx2.*(dx1+dx2))],ngrid,ngrid);

        elseif params.regularizer == "gradient"
            # Subtly different from the Laplacian. A point for future
            # enhancement is to do it better for the triangle interpolation case.
            #i = [i for i in 1:nx, j in 2:(ny-1)];
            #j = [j for i in 1:nx, j in 2:(ny-1)];
            (i,j) = meshgrid(1:nx,2:(ny-1));
            ind = j[:] .+ ny .* (i[:] .- 1);
            dy1 = dy[j[:] .- 1] ./ params.yscale;
            dy2 = dy[j[:]] ./ params.yscale;

            Areg = sparse(repeat(ind,[3]...), [ind .- 1; ind; ind .+ 1], xyRelativeStiffness[2] .* [-2 ./ (dy1 .* (dy1 .+ dy2)); 2 ./ (dy1 .* dy2); -2 ./ (dy2 .* (dy1 .+ dy2))], ngrid, ngrid);
            #Areg = sparse(repmat(ind,1,3),[ind-1,ind,ind+1], xyRelativeStiffness(2)*[-2./(dy1.*(dy1+dy2)), 2./(dy1.*dy2), -2./(dy2.*(dy1+dy2))],ngrid,ngrid);

            #i = [i for i in 2:(nx-1), j in 1:ny];
            #j = [j for i in 2:(nx-1), j in 1:ny];
            (i,j) = meshgrid(2:(nx-1),1:ny);
            ind = j[:] .+ ny .* (i[:] .- 1);
            dx1 = dx[i[:] .- 1] ./ params.xscale;
            dx2 = dx[i[:]] ./ params.xscale;
            Areg = [Areg; sparse(repeat(ind,[3]...), [ind .- ny; ind; ind .+ ny], xyRelativeStiffness[1] .* [-2 ./ (dx1 .* (dx1 .+ dx2)); 2 ./ (dx1 .* dx2); -2 ./ (dx2 .* (dx1 .+ dx2))], ngrid, ngrid)];
            #Areg = [Areg;sparse(repmat(ind,1,3),[ind-ny,ind,ind+ny], xyRelativeStiffness(1)*[-2./(dx1.*(dx1+dx2)), 2./(dx1.*dx2), -2./(dx2.*(dx1+dx2))],ngrid,ngrid)];
        end
        nreg = size(Areg,1);

        # Append the regularizer to the interpolation equations,
        # scaling the problem first. Use the 1-norm for speed.
        NA = norm(A,1);
        NR = norm(Areg,1);
        A = [A; Areg .* (smoothparam*NA/NR)];
        rhs = [rhs; zeros(nreg,1)];
        # do we have a mask to apply?
        if params.maskflag == true
            unmasked = findall(params.mask[:]);
            unmaskedCI = findall(params.mask);
        end
        # solve the full system, with regularizer attached
        if (params.solver == "backslash") | (params.solver == "\\")
            if params.maskflag
                # there is a mask to use
                zgrid=Array{Float64,2}(undef,ny,nx);
                zgrid .= NaN;
                zgrid[unmaskedCI] = A[:,unmasked] \ rhs;
            else
                # no mask
                zgrid = reshape(A \ rhs, (ny, nx));
            end
        elseif params.solver == "normal"
            # The normal equations, solved with \. Can be faster
            # for huge numbers of data points, but reasonably
            # sized grids. The regularizer makes A well conditioned
            # so the normal equations are not a terribly bad thing
            # here.
            if params.maskflag
                # there is a mask to use
                Aunmasked = A[:,unmasked];
                zgrid=Array{Float64,2}(undef,ny,nx);
                zgrid .= NaN;
                zgrid[unmaskedCI] = (Aunmasked' * Aunmasked) \ (Aunmasked' * rhs);
            else
                zgrid = reshape((A' * A) \ (A' * rhs), (ny, nx));
            end
        end  # if params.solver
        return zgrid;
    end

end

# define a histogram index function similar to matlab histc
function histc(x, xnodes)
    bin = Array{Number}(undef,length(x));
    for i = 1:length(x)
        a = findlast(x[i] .>= xnodes);
        if !isnothing(a)
            bin[i] = a;
        end
    end
    return bin;
end

# emulate the behavior of Matlab's nan function
function nan(m::Unsigned,n::Unsigned)
    nanarray=Array{Float64,2}(undef,m,n);
    nanarray .= NaN;
end

# emulate the behavior of Matlab's meshgrid
function meshgrid(xgrid::Array{<:AbstractFloat,1},ygrid::Array{<:AbstractFloat,1})
    nx = length(xgrid);
    ny = length(ygrid);
    minx = minimum(xgrid);
    maxx = maximum(xgrid);
    miny = minimum(ygrid);
    maxy = maximum(ygrid);
    dx = (maxx .- minx) ./ (nx-1);
    dy = (maxy .- miny) ./ (ny-1);
    i = [i for j in miny:dy:maxy, i in minx:dx:maxx];
    j = [j for j in miny:dy:maxy, i in minx:dx:maxx];
    return (i,j)
end

function meshgrid(xgrid::Array{<:Integer,1},ygrid::Array{<:Integer,1})
    minx = minimum(xgrid);
    maxx = maximum(xgrid);
    miny = minimum(ygrid);
    maxy = maximum(ygrid);
    i = [i for j in miny:maxy, i in minx:maxx];
    j = [j for j in miny:maxy, i in minx:maxx];
    return (i,j)
end

function meshgrid(xgrid::UnitRange{<:Integer},ygrid::UnitRange{<:Integer})
    #minx = minimum(xgrid);
    #maxx = maximum(xgrid);
    #miny = minimum(ygrid);
    #maxy = maximum(ygrid);
    i = [i for j in ygrid[1]:ygrid[end], i in xgrid[1]:xgrid[end]];
    j = [j for j in ygrid[1]:ygrid[end], i in xgrid[1]:xgrid[end]];
    return (i,j)
end

end # module GridFit


text:
struct CustomLikelihood{T} <: Likelihood{T}
    pdf::Function
    nLatent::Int
end

function treat_labels!(y,likelihood::CustomLikelihood)
    return y,likelihood.nLatent,likelihood
end

pdf(d::CustomLikelihood,y::Real,fs::AbstractVector)=d.pdf(y,fs)


text:
export Zeros

"""
`Zeros(domainType::Type, dim_in::Tuple, [codomainType::Type,] dim_out::Tuple)`

Create a `LinearOperator` which, when multiplied with an array `x` of size `dim_in`, returns an array `y` of size `dim_out` filled with zeros.

For convenience `Zeros` can be constructed from any `AbstractOperator`.

```julia
julia> Zeros(Eye(10,20))
0  ℝ^(10, 20) -> ℝ^(10, 20)

julia> Zeros([Eye(10,20) Eye(10,20)])
[0,0]  ℝ^(10, 20)  ℝ^(10, 20) -> ℝ^(10, 20)
```

"""
struct Zeros{C,N,D,M} <: LinearOperator
	dim_out::NTuple{N, Int}
	dim_in::NTuple{M, Int}
end

# Constructors
#default
Zeros(domainType::Type, dim_in::NTuple{M,Int},
      codomainType::Type, dim_out::NTuple{N,Int}) where {N,M} =
Zeros{codomainType,N,domainType,M}(dim_out,dim_in)

Zeros(domainType::Type, dim_in::NTuple{M,Int}, dim_out::NTuple{N,Int}) where {N,M} =
Zeros{domainType,N,domainType,M}(dim_out,dim_in)

function Zeros(domainType::NTuple{NN,Type}, dim_in::NTuple{NN,Tuple},
	       codomainType::Type, dim_out::Tuple) where {NN}
	HCAT([Zeros(domainType[i], dim_in[i], codomainType, dim_out) for i =1:NN]...)
end

function Zeros(domainType::Type, dim_in::Tuple,
	       codomainType::NTuple{NN,Type}, dim_out::NTuple{NN,Tuple}) where {NN}
	VCAT([Zeros(domainType, dim_in, codomainType[i], dim_out[i]) for i =1:NN]...)
end

Zeros(A::AbstractOperator) = Zeros(domainType(A),size(A,2),codomainType(A),size(A,1))

# Mappings

mul!(y::AbstractArray{C,N}, A::Zeros{C,N,D,M}, b::AbstractArray{D,M}) where {C,N,D,M}= fill!(y,zero(C))
mul!(y::AbstractArray{D,M}, A::AdjointOperator{Zeros{C,N,D,M}}, b::AbstractArray{C,N}) where {C,N,D,M}= fill!(y,zero(D))

# Properties

  domainType(L::Zeros{C,N,D,M}) where {C,N,D,M} = D
codomainType(L::Zeros{C,N,D,M}) where {C,N,D,M} = C

size(L::Zeros) = (L.dim_out, L.dim_in)

fun_name(A::Zeros)  = "0"

is_null(L::Zeros) = true


text:
module Console
export Instruction, parseInstructions, Program, step!, hasFinished

struct Instruction
    name::String
    arg::Int
end

Instruction(data::AbstractArray{<: AbstractString}) = Instruction(data[1], parse(Int, data[2]))

function parseInstructions(iter)
    [Instruction(split(line,' ')) for line in iter]
end

mutable struct Program
    instructions::Vector{Instruction}
    ip::Int
    acc::Int
end

Program(instructions::Vector{Instruction}) = Program(instructions, 1, 0)

Base.copy(prog::Program) = Program(prog.instructions, prog.ip, prog.acc)

hasFinished(prog::Program) = !(1 <= prog.ip <= length(prog.instructions))

function step!(prog::Program)
    if hasFinished(prog)
        return false
    end
    instruction = prog.instructions[prog.ip]
    if instruction.name == "nop"
        prog.ip += 1
    elseif instruction.name == "acc"
        prog.acc += instruction.arg
        prog.ip += 1
    elseif instruction.name == "jmp"
        prog.ip += instruction.arg
    else
        error("What is this instruction: $(instruction.name) ???")
    end
    return true
end

end #module


text:
module WaveformCommunications

using QuadGK, Statistics

export cosinepulse, halfsinepulse, rcpulse, srrcpulse, gaussianpulse,
       Constellation, pam, qam, psk,
       Pulse, pulseshaper, eyediag

include("pulses.jl")
include("constellations.jl")
include("utils.jl")

"""
    pulseshaper(c, pulse, nsyms = 20)

Return a function of time that generates a pulse-shaped waveform.

Symbols are drawn at random from constellation `c`. The pulse shape is specified by
`pulse`. The waveform starts at time 0 and ends after `nsyms` symbols have been generated.
"""
function pulseshaper(c::Constellation, pulse::Pulse, nsyms = 20)
    syms = rand(c, nsyms)
    return t -> begin
        s = zero(ComplexF64)
        for k = 1:nsyms
            s += syms[k].*pulse.(t.-(k-1)*pulse.Tp)
        end
        return s
    end
end

end


text:
module FastPrimeSieve

const ps = (1, 7, 11, 13, 17, 19, 23, 29)

"""
Population count of a vector of UInt8s for counting prime numbers.
See https://github.com/JuliaLang/julia/issues/34059
"""
function vec_count_ones(xs::Vector{UInt8}, n = length(xs))
    count = 0
    chunks = n ÷ sizeof(UInt)
    GC.@preserve xs begin
        ptr = Ptr{UInt}(pointer(xs))
        for i in 1:chunks
            count += count_ones(unsafe_load(ptr, i))
        end
    end

    @inbounds for i in 8chunks+1:n
        count += count_ones(xs[i])
    end

    count
end

function to_idx(x)
    x ==  1 && return 1
    x ==  7 && return 2
    x == 11 && return 3
    x == 13 && return 4
    x == 17 && return 5
    x == 19 && return 6
    x == 23 && return 7
    return 8
end

include("generate_sieving_loop.jl")
include("sieve_small.jl")
include("presieve.jl")
include("siever.jl")
include("sieve.jl")
include("parallel.jl")

end # module


text:
using Documenter, LikelihoodProfiler

makedocs(
    modules = [LikelihoodProfiler],
    sitename = "LikelihoodProfiler.jl",
    pages = [
    "Home" => "index.md",
    "Methods" => "methods.md",
    "Visualization" => "visualization.md",
    "API" => "api.md",
    ],
)

deploydocs(
    repo   = "github.com/insysbio/LikelihoodProfiler.jl.git"
)


text:
# iterator
struct ComponentIterator{C<:Colorant}
    c::C
end

eltype(::Type{ComponentIterator{C}}) where {T, C <: Colorant{T}} = T
length(::ComponentIterator{C}) where {N, C <: ColorantN{N}} = N

function Base.iterate(itr::ComponentIterator{C}, state::Int=0) where {N, C <: ColorantN{N}}
    state < 0 && return nothing
    state >= N && return nothing
    state += 1
    return (itr[state], state)
end

@inline function Base.getindex(itr::ComponentIterator{C}, i::Integer) where {N, C <: ColorantN{N}}
    N > 0 && i == 1 && return comp1(itr.c)
    N > 1 && i == 2 && return comp2(itr.c)
    N > 2 && i == 3 && return comp3(itr.c)
    N > 3 && i == 4 && return comp4(itr.c)
    N > 4 && i == 5 && return comp5(itr.c)
    throw(BoundsError(itr, i))
end
Base.getindex(itr::ComponentIterator, r::AbstractRange) = Tuple(itr[i] for i in r)
Base.getindex(itr::ComponentIterator, ::Colon) = itr

Base.firstindex(::ComponentIterator) = 1
Base.lastindex(itr::ComponentIterator) = length(itr)

function Base.BroadcastStyle(::Type{<:ComponentIterator{C}}) where {T, N, C <: Colorant{T, N}}
    Base.BroadcastStyle(NTuple{N, T})
end
Base.axes(::ComponentIterator{C}) where {N, C <: ColorantN{N}} = (Base.OneTo(N),)
Base.ndims(::Type{ComponentIterator{C}}) where {C} = 1
function Base.broadcastable(itr::ComponentIterator{C}) where {T, N, C <: Colorant{T, N}}
    (itr...,)::NTuple{N, T}
end

comps(c::Colorant) = ComponentIterator(c) # TODO: design public APIs

# comparison
_is_same_colorspace(a, b) = base_colorant_type(a) === base_colorant_type(b)
_is_same_colorspace(a::TransparentColor, b::TransparentColor) = base_color_type(a) === base_color_type(b)
_is_same_colorspace(a::AbstractGray, b::AbstractGray) = true
_is_same_colorspace(a::TransparentGray, b::TransparentGray) = true
_is_same_colorspace(a::AbstractRGB, b::AbstractRGB) = true
_is_same_colorspace(a::TransparentRGB, b::TransparentRGB) = true

function ==(a::ColorantN{N}, b::ColorantN{N}) where {N}
    _is_same_colorspace(a, b) || return false
    all(comps(a) .== comps(b))
end
==(x::Number, y::AbstractGray) = x == gray(y)
==(x::AbstractGray, y::Number) = ==(y, x)

function isapprox(a::ColorantN{N}, b::ColorantN{N}; kwargs...) where {N}
    _is_same_colorspace(a, b) || return false
    componentapprox(x, y) = isapprox(x, y; kwargs...)
    all(componentapprox.(comps(a), comps(b)))
end
isapprox(a::Colorant, b::Colorant; kwargs...) = false
isapprox(a::Number, b::AbstractGray; kwargs...) = isapprox(a, gray(b); kwargs...)
isapprox(a::AbstractGray, b::Number; kwargs...) = isapprox(b, a; kwargs...)

isless(a::AbstractGray, b::AbstractGray) = isless(gray(a), gray(b))
isless(a::AbstractGray, b::Real)         = isless(gray(a), b)
isless(a::Real,         b::AbstractGray) = isless(a,       gray(b))

<(a::AbstractGray, b::AbstractGray) = gray(a) < gray(b)
<(a::AbstractGray, b::Real)         = gray(a) < b
<(a::Real,         b::AbstractGray) = a       < gray(b)

# hash
hash(c::AbstractGray, hx::UInt) = hash(gray(c), hx)
hash(c::TransparentGray, hx::UInt) = hash(alpha(c), hash(gray(c), hx))
hash(c::AbstractRGB, hx::UInt) = hash(blue(c), hash(green(c), hash(red(c), hx)))
hash(c::TransparentRGB, hx::UInt) = hash(alpha(c), hash(blue(c), hash(green(c), hash(red(c), hx))))

Base.adjoint(c::Colorant) = c

# gamut{min,max}
# the following values are based on the bounding boxes of sRGB gamut (D65).
gamutmin(::Type{<:ColorantN{N}}) where {N} = ntuple(_ -> 0, N) # fallback

gamutmax(::Type{<:AbstractGray}) = (1,)
gamutmax(::Type{<:AbstractRGB}) = (1, 1, 1)

gamutmax(::Type{<:Union{HSV, HSL, HSI}}) = (360, 1, 1)

gamutmax(::Type{<:XYZ}) = (0.9505, 1, 1.089)

gamutmax(::Type{<:xyY}) = (0.3127, 0.3290, 1)

gamutmax(::Type{<:Lab}) = (100, 98.2343, 94.4779)
gamutmin(::Type{<:Lab}) = (0, -86.1827, -107.8602)

gamutmax(::Type{<:Luv}) = (100, 175.0150, 107.3985)
gamutmin(::Type{<:Luv}) = (0, -83.077, -134.1030)

gamutmax(::Type{<:LCHab}) = (100, 133.8076, 360)

gamutmax(::Type{<:LCHuv}) = (100, 179.0414, 360)

gamutmax(::Type{<:DIN99}) = (100.0013, 36.1753, 31.1551)
gamutmin(::Type{<:DIN99}) = (0, -27.4504, -33.3955)

gamutmax(::Type{<:DIN99d}) = (100.0002, 43.6867, 43.6318)
gamutmin(::Type{<:DIN99d}) = (0, -38.1436, -45.8498)

gamutmax(::Type{<:DIN99o}) = (99.9997, 45.5242, 44.3662)
gamutmin(::Type{<:DIN99o}) = (0, -40.1101, -40.4900)

gamutmax(::Type{<:LMS}) = (0.9493, 1.0354, 1.0872)

gamutmax(::Type{<:YIQ}) = (1, 0.5957, 0.5226)
gamutmin(::Type{<:YIQ}) = (0, -0.5957, -0.5226)

gamutmax(::Type{<:YCbCr}) = (235, 240, 240)
gamutmin(::Type{<:YCbCr}) = (16, 16, 16)

gamutmax(::Type{C}) where {C<:TransparentColor} = (gamutmax(color_type(C))..., 1)
gamutmin(::Type{C}) where {C<:TransparentColor} = (gamutmin(color_type(C))..., 0)

# rand
const Rand01Normd = Union{N0f8, N0f16, N0f32, N0f64}
const Rand01Type = Union{AbstractFloat, Rand01Normd}

# TODO: Remove the following once it is guaranteed to be implemented in FixedPointNumbers.
if which(rand, Tuple{AbstractRNG, SamplerType{<:FixedPoint}}).module === Random
    function rand(r::AbstractRNG, ::SamplerType{X}) where X <: FixedPoint
        reinterpret(X, rand(r, FixedPointNumbers.rawtype(X)))
    end
end

function rand(r::AbstractRNG, ::SamplerType{C}) where {C<:Colorant}
    rand(r, base_colorant_type(C){Float64})
end
function rand(r::AbstractRNG, ::SamplerType{C}) where {T, C<:Colorant{T}}
    Cmax = C(gamutmax(C)...)
    Cmin = C(gamutmin(C)...)
    mapc((m, n) -> T((m - n) * rand(r, floattype(T)) + n), Cmax, Cmin)
end
function rand(r::AbstractRNG, ::SamplerType{C}) where {T<:Rand01Type, C0<:AbstractGray{T},
                                                       C<:Union{C0, TransparentGray{C0, T}}}
    mapc(_ -> rand(r, T), C())
end
function rand(r::AbstractRNG, ::SamplerType{C}) where {T<:Rand01Type, C0<:AbstractRGB{T},
                                                       C<:Union{C0, TransparentRGB{C0, T}}}
    mapc(_ -> rand(r, T), C())
end
function rand(r::AbstractRNG, ::SamplerType{AGray32}) # Gray24 has little benefit of specialization.
    reinterpret(AGray32, (rand(r, UInt32) & 0xff0000ff) * 0x010101)
end
rand(r::AbstractRNG, ::SamplerType{RGB24}) = reinterpret(RGB24, rand(r, UInt32) & 0xffffff)
rand(r::AbstractRNG, ::SamplerType{ARGB32}) = reinterpret(ARGB32, rand(r, UInt32))

function rand(r::AbstractRNG, ::Type{C}, dims::Dims) where {C <: Colorant}
    CC = isconcretetype(C) ? C : base_colorant_type(C){Float64}
    rand!(r, Array{CC}(undef, dims), CC)
end

# rand!
function rand!(r::AbstractRNG, A::Array{C}, ::SamplerType{C}) where {C<:Colorant}
    rand!(r, A, SamplerType{base_colorant_type(C){Float64}}())
end
function rand!(r::AbstractRNG, A::Array{C}, ::SamplerType{C}) where {T, C<:Colorant{T}}
    A .= rand.((r,), C)
end
function _rand01!(r::AbstractRNG, A::Array{C},
               ::SamplerType{C}) where {T<:Rand01Type, C<:Colorant{T}}
    N = sizeof(C) ÷ sizeof(T)
    T0 = T <: FixedPoint ? FixedPointNumbers.rawtype(T) : T
    At = unsafe_wrap(Array, reinterpret(Ptr{T0}, pointer(A)), (N, size(A)...))
    rand!(r, At, T0)
    A
end
function rand!(r::AbstractRNG, A::Array{C},
               s::SamplerType{C}) where {T<:AbstractFloat, C<:Colorant{T}}
    _rand01!(r, A, s)
    Cmin = C(gamutmin(C)...)
    Cs = C((gamutmax(C) .- gamutmin(C))...)
    f(c) = mapc((a, b) -> T(a + b), mapc(*, c, Cs), Cmin)
    A .= f.(A)
end
function rand!(r::AbstractRNG, A::Array{C},
               s::SamplerType{C}) where {T<:Rand01Type, C<:Union{Gray{T}, AGray{T}, GrayA{T}}}
    _rand01!(r, A, s)
end
function rand!(r::AbstractRNG, A::Array{C},
               s::SamplerType{C}) where {T<:Rand01Type,
                                         C<:Union{RGB{T}, ARGB{T}, RGBA{T}, XRGB{T}, RGBX{T},
                                                  BGR{T}, ABGR{T}, BGRA{T}}}
    _rand01!(r, A, s)
end
function rand!(r::AbstractRNG, A::Array{C},
               ::SamplerType{C}) where {C<:Union{Gray24, AGray32, RGB24, ARGB32}}
    At = unsafe_wrap(Array, reinterpret(Ptr{UInt32}, pointer(A)), size(A))
    rand!(r, At, UInt32)
    if C === Gray24
        At .= At .& 0xff .* 0x010101
    elseif C === AGray32
        At .= At .& 0xff0000ff .* 0x010101
    elseif C === RGB24
        At .&= 0xffffff
    end
    A
end

# broadcast

# Mapping a function over color channels
"""
    mapc(f, rgb) -> rgbf
    mapc(f, rgb1, rgb2) -> rgbf

`mapc` applies the function `f` to each color channel of the input
color(s), returning an output color in the same colorspace.

# Examples:

    julia> mapc(x->clamp(x,0,1), RGB(-0.2,0.3,1.2))
    RGB{Float64}(0.0,0.3,1.0)

    julia> mapc(max, RGB(0.1,0.8,0.3), RGB(0.5,0.5,0.5))
    RGB{Float64}(0.5,0.8,0.5)

    julia> mapc(+, RGB(0.1,0.8,0.3), RGB(0.5,0.5,0.5))
    RGB{Float64}(0.6,1.3,0.8)
"""
@inline mapc(f, c::C) where {C<:Colorant} = base_colorant_type(C)(f.(comps(c))...)
@inline mapc(f, x::Number) = f(x)

@inline mapc(f::F, x, y) where {F} = _same_colorspace(x, y)(f.(comps(x), comps(y))...)

_same_colorspace(x::Colorant, y::Colorant) = _same_colorspace(base_colorant_type(x),
                                                              base_colorant_type(y))
_same_colorspace(::Type{C}, ::Type{C}) where {C<:Colorant} = C
@noinline _same_colorspace(::Type{C1}, ::Type{C2}) where {C1<:Colorant,C2<:Colorant} =
    throw(ArgumentError("$C1 and $C2 are from different colorspaces"))

"""
    reducec(op, v0, c)
    reducec(op, c; [init])

Reduce across color channels of `c` with the binary operator `op`.
`v0` or `init` is the neutral element used to initiate the reduction.
For grayscale,

    reducec(op, v0, c::Gray) = op(v0, comp1(c))

whereas for RGB

    reducec(op, v0, c::RGB) = op(comp3(c), op(comp2(c), op(v0, comp1(c))))

If `c` has an alpha channel, it is always the last one to be folded into the reduction.

!!! compat "ColorTypes 0.12"
    The keyword argument `init` and its omission using the default value require
    ColorTypes v0.12 or later.
"""
@inline reducec(op, v0, c::C) where {C <: Colorant} = reduce(op, Tuple(c); init=v0)
@inline reducec(op, v0, x::Number) = op(v0, x)
@inline reducec(op, c::C; kw...) where {C <: Colorant} = reduce(op, Tuple(c); kw...)
@inline reducec(op, x::Number; kw...) = reduce(op, x; kw...)

"""
    mapreducec(f, op, v0, c)
    mapreducec(f, op, c; [init])

Reduce across color channels of `c` with the binary operator `op`,
first applying `f` to each channel.
`v0` or `init` is the neutral element used to initiate the reduction.
For grayscale,

    mapreducec(f, op, v0, c::Gray) = op(v0, f(comp1(c)))

whereas for RGB

    mapreducec(f, op, v0, c::RGB) = op(f(comp3(c)), op(f(comp2(c)), op(v0, f(comp1(c)))))

If `c` has an alpha channel, it is always the last one to be folded into the reduction.

!!! compat "ColorTypes 0.12"
    The keyword argument `init` and its omission using the default value require
    ColorTypes v0.12 or later.
"""
@inline mapreducec(f, op, v0, c::C) where {C <: Colorant} = reduce(op, f.(comps(c)); init=v0)
@inline mapreducec(f, op, v0, x::Number) = op(v0, f(x))
@inline mapreducec(f, op, c::C; kw...) where {C <: Colorant} = reduce(op, f.(comps(c)); kw...)
@inline mapreducec(f, op, x::Number; kw...) = reduce(op, f(x); kw...)


text:
import GeoStatsBase: solve_single

# define a dummy solver for testing
@simsolver Dummy begin end
function solve_single(problem::SimulationProblem, var::Symbol,
                      solver::Dummy, preproc)
  pdomain = domain(problem)
  npts = npoints(pdomain)
  V = variables(problem)[var]
  fill(one(V), npts)
end


text:
using DelimitedFiles

function ReadData(fname::AbstractString)
    dataMtx = readdlm(fname, '-', String, '\n')

    data = Dict([(x, String[]) for x in dataMtx])
    for i=1:size(dataMtx)[1]
        push!(data[dataMtx[i,1]], dataMtx[i,2])
        push!(data[dataMtx[i,2]], dataMtx[i,1])
    end

    return data
end

function FindPaths1(data::Dict, start::String, path::Array{String,1})
    paths = Array{String,1}[]
    for node in data[start]
        if all(islowercase, node) && node in path
            continue
        elseif node == "end"
            push!(paths, cat(path, node, dims=1))
        else
            foreach(x -> push!(paths, x), FindPaths1(data, node, cat(path, node, dims=1)))
        end
    end

    return paths
end

# count the maximum number of times any single lower case node appears in path
function CountMaxLower(path::Array{String,1})
    return max( map(x -> count(y -> y==x, path), filter(z -> all(islowercase, z), path))...)
end

function FindPaths2(data::Dict, start::String, path::Array{String,1})
    paths = Array{String,1}[]
    for node in data[start]
        if node == "start" && length(path) > 0
            continue
        elseif node == "end"
            push!(paths, cat(path, node, dims=1))
        elseif all(islowercase, node) && CountMaxLower(path) == 2 && node in path
            continue
        else
            foreach(x -> push!(paths, x), FindPaths2(data, node, cat(path, node, dims=1)))
        end
    end

    return paths
end

function Part1(data)
    return length(FindPaths1(data, "start", ["start"]))
end

function Part2(data)
    return length(FindPaths2(data, "start", ["start"]))
end

#dataFile = "../data/test12.txt"
dataFile = "../data/day12.txt"
data = ReadData(dataFile)

println(Part1(data))
println(Part2(data))


text:
function pi(n)
	numerator = 4.0
	denominator = 1.0
	operation = 1.0
	pi = 0.0
	for i in 1:n
		pi += operation * (numerator / denominator)
		denominator += 2.0
		operation *= -1.0
	end
	return pi
end

println(@time pi(1000000000))


text:
# # Predicting credit approval by applicants' characteristics

# ### Instructions

# In this exercise we will predict the approval (yes/no) for credit applications from the applicant characteristics.

# As the data comes from a real-world log from a financial institution, both fields' names and values have been replaced with meaningless symbols to preserve anonymity.

# In detail, the attributes of this dataset are:
# - A1: b, a.
# - A2: continuous.
# - A3: continuous.
# - A4: u, y, l, t.
# - A5: g, p, gg.
# - A6: c, d, cc, i, j, k, m, r, q, w, x, e, aa, ff.
# - A7: v, h, bb, j, n, z, dd, ff, o.
# - A8: continuous.
# - A9: t, f.
# - A10: t, f.
# - A11: continuous.
# - A12: t, f.
# - A13: g, p, s.
# - A14: continuous.
# - A15: continuous.
# - A16: +,- (class attribute) - what we want to predict

# Further information concerning this dataset can be found online on the [UCI Machine Learning Repository dedicated page](https://archive.ics.uci.edu/ml/datasets/Credit+Approval)

# Our prediction concerns the positive or negative outcome of the credit application.

# While you can use any supervised ML algorithm, I suggest the [`Random Forests`](https://sylvaticus.github.io/BetaML.jl/dev/Trees.html) from BetaML because of their ease of use and the presence of numerous categorical data and missing data that would require additional work with most other algorithms.

# ------------------------------------------------------------------------------
# ### 1) Start by setting the working directory to the directory of this file and activate it.
# If you have the provided `Manifest.toml` file in the directory, just run `Pkg.instantiate()`, otherwise manually add the packages Pipe, HTTP, CSV, DataFrames, Plots and BetaML
# Also, seed the random seed with the integer `123`.
cd(@__DIR__)
using Pkg
Pkg.activate(".")
# If using a Julia version different than 1.7 please uncomment and run the following line (reproductibility guarantee will hower be lost)
# Pkg.resolve()
Pkg.instantiate()
using Random
Random.seed!(123)

# ------------------------------------------------------------------------------
# ### 2) Load the packages/modules Pipe, HTTP, CSV, DataFrames, Plots, BetaML
using  Pipe, HTTP, CSV, DataFrames, Plots, BetaML

# ### 3) Load from internet or from local file the input data.
# You can use a pipeline from HTTP.get() to CSV.File to finally a DataFrame.
# Use the parameter `missingstring="?"` in the `CSV.File()` call.

dataURL = "https://archive.ics.uci.edu/ml/machine-learning-databases/credit-screening/crx.data"

# [...] write your code here...

# ------------------------------------------------------------------------------
# ### 4) Now create the X matrix and Y vector
# Create the X matrix of features using the first to the second-to-last column of the data you loaded above and the Y vector by taking the last column.
# If you use the random forests algorithm suggested above, the only data preprocessing you need to do is to convert the X from a DataFrame to a Matrix and to `collect` the Y to a vector. Otherwise be sure to encode the categorical data, skip or impute the missing data and scale the feature matrix as required by the algorithm you employ.

# [...] write your code here...

# ------------------------------------------------------------------------------
# ### 5) Partition your data in (xtrain,xtest) and (ytrain,ytest)
# (e.g. using 80% for the training and 20% for testing)
# You can use the BetaML [`partition()`](https://sylvaticus.github.io/BetaML.jl/dev/Utils.html#BetaML.Api.partition-Union{Tuple{T},%20Tuple{AbstractVector{T},%20AbstractVector{Float64}}}%20where%20T%3C:AbstractArray) function.
# Be sure to shuffle your data if you didn't do it earlier! (that's done by default)

# [...] write your code here...

# ------------------------------------------------------------------------------
# ### 6) (optional but suggested) Find the best hyper-parameters for your model, i.e. the ones that lead to the highest accuracy under the records not used for training.
# You can use the [`crossValidation`](https://sylvaticus.github.io/BetaML.jl/dev/Utils.html#BetaML.Utils.crossValidation) function here.
# The idea is that for each hyper-parameter you have a range of possible values, and for each hyper-parameter, you first set `bestAcc=0.0` and then loop on each possible value, you run crossValidation with that particular value to compute the average training accuracy with that specific value under different data samples, and if it is better than the current `bestAcc` you save it as the new `bestAcc` and the parameter value as the best value for that specific hyper-parameter.
# After you have found the best hyper-parameter value for one specific hyper-parameter, you can switch to the second hyper-parameter repeating the procedure but using the best value for the first hyper-parameter that you found earlier, and you continue with the other hyper-parameters.
# Note that if you limit the hyper-parameter space sufficiently you could also directly loop over all the possible combinations of hyper-parameters.

# If you use the Random Forests from BetaML consider the following hyper-parameter ranges:

nTrees_range             = 20:5:60
splittingCriterion_range = [gini,entropy]
maxDepth_range           = [10,15,20,25,30,500]
minRecords_range         = [1,2,3,4,5]
maxFeatures_range        = [2,3,4,5,6]
β_range                  = [0.0,0.5,1,2,5,10,20,50,100]

# To train a Random Forest in BetaML use:
# `myForest = buildForest(xtrain,ytrain, nTrees; <other hyper-parameters>)`
# And then to predict and compute the accuracy use:
# ```julia
# ŷtrain=predict(myforest,xtrain)
# trainAccuracy = accuracy(ŷtrain,ytrain)
# ```
# This activity is "semi-optional" as Random Forests have very good default values, so the gain you will likely obtain with tuning the various hyper-parameters is not expected to be very high.  But it is a good exercise to arrive at this result by yourself !

# [...] write your code here...

# ------------------------------------------------------------------------------
# ### 7) Perform the final training with the best hyperparameters and compute the accuracy on the test set
# If you have chosen good hyperparameters, your accuracy should be in the 98%-99% range for training and 81%-89% range for testing

# [...] write your code here...

max_stars_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


text:
using Test, Pkg

@testset "InputOutput" begin
    all_tests = isempty(ARGS) || "all" in ARGS ? true : false
    for submodule in ["Writers"]
        if all_tests || "$submodule" in ARGS || "Utilities" in ARGS
            include_test(submodule)
        end
    end

end


text:
# Encoding support
#
# Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones
# Licensed under MIT License, see LICENSE.md
#
# Encodings inspired from collaboration with @nalimilan (Milan Bouchet-Valat) on
# [StringEncodings](https://github.com/nalimilan/StringEncodings.jl)

@api public  Encoding, Native1Byte, UTF8Encoding
@api develop encoding_types

struct Encoding{Enc} end

"""List of installed encodings"""
const encoding_types = []

Encoding(s) = Encoding{Symbol(s)}()

const Native1Byte  = Encoding{:Byte}
const UTF8Encoding = Encoding{:UTF8}

push!(encoding_types, Native1Byte, UTF8Encoding)

# Allow handling different endian encodings

for (n, l, b, s) in (("2Byte", :LE2,     :BE2,     "16-bit"),
                     ("4Byte", :LE4,     :BE4,     "32-bit"),
                     ("UTF16", :UTF16LE, :UTF16BE, "UTF-16"))
    nat, swp = BIG_ENDIAN ? (b, l) : (l, b)
    natnam = symstr("Native",  n)
    swpnam = symstr("Swapped", n)
    @eval const $natnam = Encoding{$(quotesym("N", n))}
    @eval const $swpnam = Encoding{$(quotesym("S", n))}
    @eval const $nat = $natnam
    @eval const $swp = $swpnam
    @eval push!(encoding_types, $natnam, $swpnam)
    @eval @api public $natnam, $swpnam
end

show(io::IO, ::Type{Encoding{S}}) where {S}  = print(io, "Encoding{:", string(S), "}")

encoding(::Type{T}) where {T<:AbstractString} = encoding(cse(T)) # Default unless overridden
encoding(str::AbstractString) = encoding(cse(str))


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

using Markdown
using InteractiveUtils

# ╔═╡ 0e0ff076-f409-11ea-3774-250375ce2cc0
function my_sum1(a::Int, b::Int)
	a + b
end

# ╔═╡ 20ac3636-f409-11ea-319f-05c98f069ef1
function my_sum2(a::Int, b::Int)
	a + b
end

# ╔═╡ 1a89f222-f40a-11ea-1d2c-67a8f61101b2
function my_subtract(a::Int, b::Int)
	a - b
end

# ╔═╡ Cell order:
# ╠═0e0ff076-f409-11ea-3774-250375ce2cc0
# ╠═20ac3636-f409-11ea-319f-05c98f069ef1
# ╠═1a89f222-f40a-11ea-1d2c-67a8f61101b2


text:
const Tms = typeof(1.0f0*ms)

"""
    CellsTrial(t, dFoF)

Store cellular responses for a single trial. `t` is the list of frame times (in units of time, see Unitful.jl) and
`dFoF` is a matrix with one row for each time in `t`, and one column per cell.

# Examples

```jldoctest demo; setup=:(using EcoTrialStructure)
julia> dFoF = [ 0.1 0.8;
               -0.1 0.7;
                0.2 0.6;
                0.1 0.5;
                0.0 0.4];

julia> ct = CellsTrial((100:100:500) * ms, dFoF)
2 cells with 5 timepoints
```

It's straightforward to extract the dFoF values over a particular time interval:

```jldoctest demo
julia> tframes, df = ct[175ms..310ms, :];  # frames range nearest to given start..stop times

julia> tframes
200.0f0 ms..300.0f0 ms

julia> df
2×2 Matrix{Float64}:
 -0.1  0.7
  0.2  0.6
```
This extracted the `dFoF` values recorded between 200ms and 300ms, inclusive.

In other cases, you might want a specific number of frames, starting at a particular time:

```jldoctest demo
julia> tframes, df = ct[FrameSeq(175ms, 3), :];  # start nearest to 175ms, and grab 3 frames' worth

julia> tframes
200.0f0 ms..400.0f0 ms

julia> df
3×2 OffsetArray(::Matrix{Float64}, 0:2, 1:2) with eltype Float64 with indices 0:2×1:2:
 -0.1  0.7
  0.2  0.6
  0.1  0.5
```
"""
struct CellsTrial{T<:AbstractFloat}
    t::Vector{Tms}
    dFoF::Matrix{T}

    function CellsTrial{T}(t, dFoF) where T
        issorted(t) || throw(ArgumentError("frame times must be increasing"))
        axes(dFoF, 1) == axes(t, 1) || throw(DimensionMismatch("frame times and dFoF must match, got sizes $(length(t)) and $(size(dFoF))"))
        return new{T}(t, dFoF)
    end
end
CellsTrial(t, dFoF::AbstractMatrix{T}) where T<:AbstractFloat = CellsTrial{T}(t, dFoF)
CellsTrial(t, dFoF::AbstractMatrix) = CellsTrial{Float64}(t, dFoF)

Base.:(==)(a::CellsTrial, b::CellsTrial) = a.t == b.t && a.dFoF == b.dFoF
Base.isequal(a::CellsTrial, b::CellsTrial) = isequal(a.t, b.t) && isequal(a.dFoF, b.dFoF)

const hash_seed_CellsTrial = Int === Int64 ? 0x80450cff5158f34c : 0xbcdbf5e0
Base.hash(ct::CellsTrial, h::UInt) = hash(hash_seed_CellsTrial, hash(ct.t, hash(ct.dFoF, h)))

Base.show(io::IO, ct::CellsTrial) = print(io, size(ct.dFoF, 2), " cells with ", length(ct.t), " timepoints")

# indexing (primarily, selecting frames based on time)

Base.getindex(ct::CellsTrial, ti::Union{Colon,AbstractVector{<:Integer}}, ci) = (ct.t[ti], ct.dFoF[ti, ci])

function Base.getindex(ct::CellsTrial, ti::AbstractInterval, ci)
    ibegin, iend = idxof(ct.t, minimum(ti)), idxof(ct.t, maximum(ti))
    return ct.t[ibegin] .. ct.t[iend], ct.dFoF[ibegin:iend, ci]
end

"""
    FrameSeq(tstart, nframes)
    FrameSeq(tstart, offsetrange)
    FrameSeq(eventfield::Symbol, nframes)

A sequence of `nframes` frames starting at the nearest timepoint to `tstart`.
Optionally specify a range of frames offset from `tstart`.

```jldoctest; setup=:(using EcoTrialStructure)
julia> fs = FrameSeq(100ms, 5)
FrameSeq(100.0f0 ms, 0:4)

julia> fs = FrameSeq(100ms, -2:4)
FrameSeq(100.0f0 ms, -2:4)
```

See [`CellsTrial`](@ref) for an example using this in indexing.

Alternatively, this can be constructed specifying a particular fieldname of [`EventTiming`](@ref),
in which case the concrete timing can be deferred until a later time based on a specific
trial:

```jldoctest; setup=:(using EcoTrialStructure)
julia> fs = FrameSeq(:go, 5)
FrameSeq(:go, 0:4)

julia> et = EventTiming(0ms, 100ms, 400ms, 450ms, 837ms, 1.2s)
EventTiming(trial_start=0.0f0 ms, offer_on=100.0f0 ms, offer_off=400.0f0 ms, go=450.0f0 ms, choice=837.0f0 ms, trial_end=1200.0f0 ms)

julia> fs(et)
FrameSeq(450.0f0 ms, 0:4)
```
"""
struct FrameSeq
    start::Union{Tms,Symbol}
    idx::UnitRange{Int}
end
FrameSeq(start::Unitful.Quantity, len::Integer) = FrameSeq(Tms(start), 0:len-1)
FrameSeq(start::Symbol, len::Integer) = FrameSeq(start, 0:len-1)
FrameSeq(start::Unitful.Quantity, idx::AbstractUnitRange) = FrameSeq(Tms(start), idx)

function Base.getindex(ct::CellsTrial, ti::FrameSeq, ci)
    isdeferred(ti) && throw(ArgumentError("indexing requires a concrete `FrameSeq`, use `fs(et::EventTiming)`"))
    irange = idxsof(ct.t, ti)
    return ct.t[first(irange)] .. ct.t[last(irange)], ct.dFoF[irange, ci]
end

Base.length(fs::FrameSeq) = length(fs.idx)
Base.axes(fs::FrameSeq) = (Base.IdentityUnitRange(fs.idx),)
Base.axes(fs::FrameSeq, d) = axes(fs)[d]

function Base.checkbounds(::Type{Bool}, ct::CellsTrial, ti::FrameSeq, ci)
    r1, r2 = _idxof(ct.t, ti.start)
    checkbounds(Bool, ct.t, r1) && checkbounds(Bool, ct.t, r2) || return false
    istart = _idxof(ct.t, ti.start, r1, r2)
    return checkbounds(Bool, ct.dFoF, istart .+ ti.idx, ci)
end

"""
    OfferType(nA, nB, leftA::Bool)

Encode the offers (`nA` and `nB` are the number of drops of A and B, respectively),
and whether A was on the left.

If you have a list of `OfferType`s, you can `sort` them. High-A trials will come before high-B trials,
and left before right.
"""
struct OfferType
    nA::Int8
    nB::Int8
    leftA::Bool
end

# This constructor makes it possible to copy/paste the output of the `show` method below
OfferType(; nA, nB, leftA) = OfferType(nA, nB, leftA)

OfferType(tt::OfferType) = tt

Base.show(io::IO, tt::OfferType) = print(io, "OfferType(nA=", tt.nA, ", nB=", tt.nB, ", leftA=", tt.leftA, ")")

function Base.isless(a::OfferType, b::OfferType)
    ratio(tt::OfferType) = tt.nB/tt.nA

    # Ordering puts high-A trials before high-B trials
    rA, rB = ratio(a), ratio(b)
    isless(rA, rB) && return true
    isless(rB, rA) && return false
    a.nA > b.nA && return true
    a.nA < b.nA && return false
    a.nB > b.nB && return false
    a.nB < b.nB && return true
    return a.leftA > b.leftA   # if all else are equal, put left before right
end

"""
    TrialResult(nA, nB, leftA::Bool, choseA::Union{Bool,Missing})
    TrialResult(tt::OfferType, choseA::Union{Bool,Missing})

Encode the offer configuration (see [`OfferType`](@ref)), and whether the animal chose A, B, or failed to make a choice
(`choseA = true | false | missing`, respectively).
"""
struct TrialResult
    tt::OfferType
    choseA::Union{Bool,Missing}   # missing if the animal didn't lick
end
TrialResult(nA, nB, leftA, choseA) = TrialResult(OfferType(nA, nB, leftA), choseA)
TrialResult(; nA, nB, leftA, choseA) = TrialResult(nA, nB, leftA, choseA)
TrialResult(tr::TrialResult) = tr

function Base.getproperty(tr::TrialResult, name::Symbol)
    name === :nA && return getfield(tr, :tt).nA
    name === :nB && return getfield(tr, :tt).nB
    name === :leftA && return getfield(tr, :tt).leftA
    name === :tt && return getfield(tr, :tt)
    return getfield(tr, name)
end

Base.show(io::IO, tr::TrialResult) = print(io, "TrialResult(nA=", tr.nA, ", nB=", tr.nB, ", leftA=", tr.leftA, ", choseA=", tr.choseA, ")")

OfferType(tr::TrialResult) = tr.tt

function Base.isless(a::TrialResult, b::TrialResult)
    tta, ttb = OfferType(a), OfferType(b)
    isless(tta, ttb) && return true
    isless(ttb, tta) && return false
    # This implements `true` > `missing` > `false` logic
    rank(choice) = isa(choice, Bool) ? Float64(choice) : 0.5
    return isless(rank(b.choseA), rank(a.choseA))
end

"""
    EventTiming(trial_start, offer_on, offer_off, go, choice, trial_end)

Specify the timing of events during a trial. All times should be in physical units (`s` or `ms`).

# Examples

```jldoctest; setup=:(using EcoTrialStructure)
julia> EventTiming(0ms, 100ms, 400ms, 450ms, 837ms, 1.2s)
EventTiming(trial_start=0.0f0 ms, offer_on=100.0f0 ms, offer_off=400.0f0 ms, go=450.0f0 ms, choice=837.0f0 ms, trial_end=1200.0f0 ms)
```
"""
struct EventTiming
    trial_start::Tms
    offer_on::Tms
    offer_off::Tms
    go::Tms
    choice::Tms
    trial_end::Tms
end
EventTiming(; trial_start, offer_on, offer_off, go, choice, trial_end) =
    EventTiming(trial_start, offer_on, offer_off, go, choice, trial_end)

Base.show(io::IO, et::EventTiming) = print(io, "EventTiming(trial_start=", et.trial_start, ", offer_on=", et.offer_on, ", offer_off=", et.offer_off, ", go=", et.go, ", choice=", et.choice, ", trial_end=", et.trial_end, ")")

function (fs::FrameSeq)(et::EventTiming)
    start = fs.start
    return isa(start, Symbol) ? FrameSeq(getfield(et, fs.start), fs.idx) : fs
end


text:
using TaylorSeries

const order = 20
const x, y, z, w = set_variables(Int128, "x", numvars=4, order=2order)

function fateman1(degree::Int)
    T = Int128
    oneH = HomogeneousPolynomial(one(T), 0)
    # s = 1 + x + y + z + w
    s = TaylorN( [oneH, HomogeneousPolynomial([one(T),one(T),one(T),one(T)],1)], degree )
    s = s^degree
    # s is converted to order 2*ndeg
    s = TaylorN(s, 2*degree)

    s * ( s+TaylorN(oneH, 2*degree) )
end

function fateman2(degree::Int)
    T = Int128
    oneH = HomogeneousPolynomial(one(T), 0)
    # s = 1 + x + y + z + w
    s = TaylorN( [oneH, HomogeneousPolynomial([one(T),one(T),one(T),one(T)],1)], degree )
    s = s^degree
    # s is converted to order 2*ndeg
    s = TaylorN(s, 2*degree)
    return s^2 + s
end

function fateman3(degree::Int)
    s = x + y + z + w + 1
    s = s^degree
    s * (s+1)
end

function fateman4(degree::Int)
    s = x + y + z + w + 1
    s = s^degree
    s^2 + s
end

function run_fateman(N)
    results = Any[]
    nn = 5
    for f in (fateman1, fateman2, fateman3, fateman4)
        f(0)
        println("Running $f")
        @time result = f(N)
        # push!(results, result) # This may take a lot of memory
        t = Inf
        tav = 0.0
        for i = 1:nn
            ti = @elapsed f(N)
            tav += ti
            t = min(t,ti)
        end
        println("\tAverage time of $nn runs: ", tav/nn)
        println("\tMinimum time of $nn runs: ", t)
    end
    results
end

println("Running Fateman with order $order...")

results = run_fateman(order);

println("Done.")

# @assert results[1] == results[2] == results[3] == results[4]


text:
include("common.jl")
@testset "saveat" begin
    prob = prob_dde_1delay_long
    alg = MethodOfSteps(Tsit5())

    # reference integrator and solution
    dde_int = init(prob, alg)
    sol = solve!(dde_int)

    @testset "reference" begin
        # solution equals solution of ODE integrator
        @test sol.t == dde_int.sol.t
        @test sol.u == dde_int.sol.u
    end

    # do not save every step
    @testset "not every step (save_start=$save_start)" for save_start in (false, true)
        # for time(s) as scalar (implicitly adds end point as well!) and vectors
        for saveat in (25.0, [25.0, 50.0, 75.0])
            ## minimal ODE solution
            dde_int_min = init(prob, alg; saveat=saveat, save_start=save_start,
                               minimal_solution=true)

            # solution of ODE integrator will be reduced
            @test dde_int_min.saveat !== nothing

            sol_min = solve!(dde_int_min)

            # time point of solution
            @test sol_min.t == (save_start ? [0.0, 25.0, 50.0, 75.0, 100.0] :
                                [25.0, 50.0, 75.0, 100.0])

            # solution of ODE integrator is reduced:
            # [0.0, ≈24.67, ≈25.89, ≈49.23, ≈50.47, ≈73.91, ≈75.14, 100.0]
            @test dde_int_min.sol.t ≈ [0.0, 24.67, 25.89, 49.23, 50.47, 73.91,
                                       75.14, 100.0] atol=1.1e-2

            # solution lies on interpolation of full solution above
            @test sol(sol_min.t).u == sol_min.u

            ## full ODE solution
            dde_int_full = init(prob, alg; saveat=saveat, save_start=save_start,
                                minimal_solution=false)

            # solution of ODE integrator will not be reduced
            @test dde_int_full.saveat === nothing

            sol_full = solve!(dde_int_full)

            # solution of ODE integrator equals full solution above
            @test sol.t == dde_int_full.sol.t && sol.u == dde_int_full.sol.u

            # solution equals reduced solution above
            @test sol_min.t == sol_full.t && sol_min.u == sol_full.u

            ## dense interpolation
            dde_int_dense = init(prob, alg; saveat=saveat, save_start=save_start,
                                 dense=true)

            # solution of ODE integrator will not be reduced
            @test dde_int_dense.saveat === nothing

            sol_dense = solve!(dde_int_dense)

            # time steps of solution
            @test sol_dense.t == (save_start ? [0.0, 25.0, 50.0, 75.0, 100.0] :
                                  [25.0, 50.0, 75.0, 100.0])

            # solution of ODE integrator equals full solution above
            @test sol.t == dde_int_dense.sol.t && sol.u == dde_int_dense.sol.u

            # solution lies on interpolation of full solution above
            @test sol(sol_dense.t).u == sol_dense.u

            # full solution above lies on interpolation of solution
            @test sol_dense(sol.t).u == sol.u
        end
    end

    # save every step
    @testset "every step (save_start=$save_start)" for save_start in (false, true)
        # for time(s) as scalar (implicitly adds end point as well!) and vectors
        for saveat in (25.0, [25.0, 50.0, 75.0])
            ## minimal ODE solution (is not possible to minimize ODE solution actually!)
            dde_int_min = init(prob, alg; saveat=saveat, save_everystep=true,
                               save_start=save_start, minimal_solution=true)

            # solution of ODE integrator will not be reduced
            @test dde_int_min.saveat === nothing

            ## full ODE solution
            dde_int_full = init(prob, alg; saveat=saveat, save_everystep=true,
                                save_start=save_start, minimal_solution=false)

            # solution of ODE integrator will not be reduced
            @test dde_int_full.saveat === nothing

            sol_full = solve!(dde_int_full)

            # time steps of solution
            @test symdiff(sol.t, sol_full.t) == (save_start ? [25.0, 50.0, 75.0] :
                                                 [0.0, 25.0, 50.0, 75.0])

            # solution of ODE integrator equals full solution above
            @test sol.t == dde_int_full.sol.t && sol.u == dde_int_full.sol.u

            # solution lies on interpolation of full solution above
            @test sol(sol_full.t).u == sol_full.u

            ## dense interpolation
            dde_int_dense = init(prob, alg; saveat=saveat, save_everystep=true,
                                 save_start=save_start, dense=true)

            # solution of ODE integrator will not be reduced
            @test dde_int_dense.saveat === nothing

            sol_dense = solve!(dde_int_dense)

            # time steps of solution
            @test symdiff(sol.t, sol_dense.t) == (save_start ? [25.0, 50.0, 75.0] :
                                                  [0.0, 25.0, 50.0, 75.0])

            # solution of ODE integrator equals full solution above
            @test sol.t == dde_int_dense.sol.t && sol.u == dde_int_dense.sol.u

            # solution lies on interpolation of full solution above
            @test sol(sol_dense.t).u == sol_dense.u

            # full solution above lies on interpolation of solution
            @test sol_dense(sol.t).u == sol.u
        end
    end
end


text:
## conversion between spectrum envelop and mel-cepstrum

# TODO(ryuichi) import from MelGeneralizedCepstrums or SPTK
# copied from r9y9/MelGeneralizedCepstrums.jl
function freqt!(wc::AbstractVector, c::AbstractVector, α;
                prev=Array(eltype(wc), length(wc)))
    fill!(wc, zero(eltype(wc)))
    dst_order = length(wc) - 1

    m1 = length(c)-1
    for i in -m1:0
        copy!(prev, wc)
        if dst_order >= 0
            @inbounds wc[1] = c[-i+1] + α*prev[1]
        end
        if dst_order >= 1
            wc[2] = (1.0-α*α)*prev[1] + α*prev[2]
        end
        for m=3:dst_order+1
            @inbounds wc[m] = prev[m-1] + α*(prev[m] - wc[m-1])
        end
    end

    wc
end

function freqt(c::AbstractVector, order=25, α=0.35)
    wc = Array(eltype(c), order+1)
    freqt!(wc, c, α)
end

"""
$(SIGNATURES)

sp2mc converts power spectrum envelope to mel-cepstrum

``
|X(\\omega)|^{2} -> c\_{\\alpha}(m)
``
"""
function sp2mc(powerspec::AbstractVector,
               order,
               α; # all-pass constant
               fftlen::Int=(length(powerspec)-1)*2
    )
    # |X(ω)|² -> log(|X(ω)²|)
    logperiodogram = log(powerspec)

    # transform log-periodogram to real cepstrum
    # log(|X(ω)|²) -> c(m)
    c = real(irfft(logperiodogram, fftlen))
    c[1] /= 2.0

    # c(m) -> cₐ(m)
    freqt(c, order, α)
end

"""
$(SIGNATURES)

mc2sp converts mel-cepstrum to power spectrum envelope.

``
c\_{\\alpha}(m) -> |X(\\omega)|^{2}
``

equivalent: `exp(2real(MelGeneralizedCepstrums.mgc2sp(mc, α, 0.0, fftlen)))`
Note that `MelGeneralizedCepstrums.mgc2sp` returns log magnitude spectrum.
"""
function mc2sp{T}(mc::AbstractVector{T}, α, fftlen)
    # back to cepstrum from mel-cesptrum
    # cₐ(m) -> c(m)
    c = freqt(mc, fftlen>>1, -α)
    c[1] *= 2.0

    symc = zeros(T, fftlen)
    copy!(symc, c)
    for i=1:length(c)-1
        symc[end-i+1] = symc[i+1]
    end

    # back to power spectrum
    # c(m) -> log(|X(ω)|²) -> |X(ω)|²
    exp(real(rfft(symc)))
end

# extend vector to vector transformation for matrix input
for f in [:sp2mc,
          :mc2sp,
          ]
    @eval begin
        function $f(x::AbstractMatrix, args...; kargs...)
            outbuf = $f(view(x, :, 1), args...; kargs...)
            ret = Array(eltype(outbuf), length(outbuf), size(x, 2))
            copy!(ret, 1, outbuf, 1, length(outbuf))
            for i = 2:size(x, 2)
                @inbounds ret[:, i] = $f(view(x, :, i), args...; kargs...)
            end
            ret
        end
    end
end


text:
s = """
if VERSION < v"0.4-" && VERSION > v"0.1"
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4")
else
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4")
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test occursin("Reachable by 0.3", msgs[1].message)
@test occursin("Unreachable by 0.4", msgs[2].message)
@test occursin("Unreachable by 0.3", msgs[3].message)
@test occursin("Reachable by 0.4", msgs[4].message)

s = """
test() = true
if VERSION < v"0.4-" && test()
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4.0-dev+1833")
else
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4.0-dev+1833")
end
"""
msgs = lintstr(s)
@test length(msgs) == 4
@test msgs[1].code == :I271
@test occursin("Reachable by 0.3", msgs[1].message)
@test occursin("Unreachable by 0.4", msgs[2].message)
@test occursin("Reachable by 0.3", msgs[3].message)
@test occursin("Reachable by 0.4", msgs[4].message) # we cannot prove unreachable

s = """
test() = true
if VERSION < v"0.4-" || test()
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4.0-dev+1833")
else
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4.0-dev+1833")
end
"""
msgs = lintstr(s)
@test length(msgs) == 4
@test msgs[1].code == :I271
@test occursin("Reachable by 0.3", msgs[1].message)
@test occursin("Reachable by 0.4", msgs[2].message)
@test occursin("Unreachable by 0.3", msgs[3].message)
@test occursin("Reachable by 0.4", msgs[4].message)

# testing `||`, should be rare
s = """
if VERSION >= v"0.4-" || VERSION < v"0.3"
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4.0-dev+1833")
else
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4.0-dev+1833")
end
"""
msgs = lintstr(s)
@test length(msgs) == 4
@test msgs[1].code == :I271
@test occursin("Unreachable by 0.3", msgs[1].message)
@test occursin("Reachable by 0.4", msgs[2].message)
@test occursin("Reachable by 0.3", msgs[3].message)
@test occursin("Unreachable by 0.4", msgs[4].message)

s = """
if !(VERSION >= v"0.4-")
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4.0-dev+1833")
else
    @lintpragma("Info version 0.3")
    @lintpragma("Info version 0.4.0-dev+1833")
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test occursin("Reachable by 0.3", msgs[1].message)
@test occursin("Unreachable by 0.4", msgs[2].message)
@test occursin("Unreachable by 0.3", msgs[3].message)
@test occursin("Reachable by 0.4", msgs[4].message)