text:
"""

    fit(d::Type{SGS}, x::AbstractArray)

Fit an SGS distribution to the time series, x.

# Examples
```jldoctest
julia> fit(SGS, air())
SGS{Float64}(E=0.36623331624359673, b=0.7133796999284913, g=-1.1648873051087967)

```

"""
function fit(d::Type{SGS}, x::AbstractArray)

    # Calculate the moments.
    variance = var(x, corrected=false)
    skew = skewness(x)
    kurt = kurtosis(x)

    # Test for the skewness-kurtosis relationship constraint.
    skewbound = 1.5 * skew^2
    @debug "K-S inequality constraint: $kurt > $skewbound"
    kurt > skewbound || @warn "K-S inequality constraint violated."

    # Fit the SGS distribution using the method of moments.
    return fit_mm(d, variance, skew, kurt)

end

"""

    fit(d::Type{SGS}, variance::Real, skew::Real, kurt::Real)

Fit an SGS distribution to the desired variance, skewness and kurtosis.

# Examples
```jldoctest
julia> fit(SGS, 1.0, 1.0, 5.0)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

```

"""
function fit(d::Type{SGS}, variance::Real, skew::Real, kurt::Real)

    # Test for the skewness-kurtosis relationship constraint.
    skewbound = 1.5 * skew^2
    kurt > skewbound || @warn "K-S inequality constraint violated: $kurt > $skewbound"

    # Fit the SGS distribution using the method of moments.
    return fit_mm(d, variance, skew, kurt)

end

"""

    fit(d::Type{SGS}, variance::Real, skew::Real, kurt::Real)

Use a method of moments to fit an SGS distribution to the desired variance,
skewness and kurtosis.

# Examples
```jldoctest
julia> fit_mm(SGS, 1.0, 1.0, 5.0)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

```

"""
function fit_mm(d::Type{SGS}, variance::Real, skew::Real, kurt::Real)

    # Calculate the SGS parameters using the calculated moments.
    E = param_E(skew, kurt)
    b = param_b(variance, skew, E)
    g = param_g(variance, skew, E)

    @debug "E = $E"
    @debug "b = $b"
    @debug "g = $g"
    @debug "variance = $variance"
    @debug "skewness = $skew"
    @debug "kurtosis = $kurt"

    return d(E, b, g)

end

"""

    param_E(skew::Real, kurt::Real)

Solve for the parameter E in Equation 8(a) of Sardeshmukh et al. (2015).

# Examples
```jldoctest
julia> param_E(1.0, 5.0)
0.6236095644623235

```

"""
function param_E(skew::Real, kurt::Real)

    E = _param_E(skew, kurt)

    # TODO: Should the value of the kurtosis changed globally or just for the
    # purpose of calculating E?

    # Must implement the epsilon constraint.
    epsilon_constraint(skew, E) || @warn "Epsilon constraint for parameter E not satisfied, incrementing kurtosis"
    # epsilon_constraint(skew, E) || @debug "Epsilon constraint for parameter E not satisfied, incrementing kurtosis"
    while epsilon_constraint(skew, E) == false
        # println("epsilon constraint unmet, incrementing kurtosis")
        kurt += 0.1
        E = _param_E(skew, kurt)
        @debug "E = $E, kurt = $kurt"
        if kurt > 15
            @warn "E = $E  skew=$skew  kurt = $kurt"
            return E   # not great, but needs to be investigated
        end
    end

    return E

end

"""

    param_E(skew::AbstractArray, kurt::AbstractArray)

Solve for the parameter E in Equation 8(a) of Sardeshmukh et al. (2015).

# Examples
```jldoctest
julia> param_E([1.0, 1.2, 1.4], [3.0, 4.0, 5.0])
3-element Array{Float64,1}:
 0.5
 0.5186577368103327
 0.5220026556319158

```

"""
function param_E(skew::AbstractArray, kurt::AbstractArray)

    E = _param_E.(skew, kurt)

    # TODO: Must implement the epsilon constraint.

    return E

end

function _param_E(skew::Real, kurt::Real)

    E2 = (2.0 / 3.0) * ((kurt - (1.5 * skew^2)) / (kurt - skew^2 + 2.0))

    if E2 ≤ 0.0
        @debug "Constraint E ≥ 0 not satisfied, setting E to NaN"
        return NaN
    else
        return sqrt(E2)
    end

end

function _param_E(skew::AbstractArray, kurt::AbstractArray)

    E = sqrt.((2.0 / 3.0) * ((kurt .- (1.5 .* skew .^2)) / (kurt .- skew .^2 + 2.0)))
    return E

    # TODO: Implement the positive definite constraint on E here.

end

"""

    param_g(variance::Real, skew::Real, E::Real)

Solve for the parameter g in Equation 8(b) of Sardeshmukh et al. (2015).

# Examples
```jldoctest
julia> E = 0.6236095644623235
0.6236095644623235

julia> param_g(1.0, 1.0, E)
0.48997894350611143

```

"""
function param_g(variance::Real, skew::Real, E::Real)

    return skew * sqrt(variance) * ((1.0 - E^2) / (2.0 * E))

end

"""

    param_g(variance::AbstractArray, skew::AbstractArray, E::AbstractArray)

Solve for the parameter g in Equation 8(b) of Sardeshmukh et al. (2015).

# Examples
```jldoctest
julia> E = 0.6236095644623235
0.6236095644623235

julia> param_g(1.0, 1.0, E)
0.48997894350611143

```

"""
function param_g(variance::AbstractArray, skew::AbstractArray, E::AbstractArray)

    return skew .* sqrt.(variance) .* ((1.0 .- E .^2) / (2.0 .* E))

end

"""

    param_b(variance::Real, skew::Real, E::Real)

Solve for the parameter b in Equation 8(c) of Sardeshmukh et al. (2015).

# Examples
```jldoctest
julia> E = 0.6236095644623235
0.6236095644623235

julia> param_b(1.0, 1.0, E)
1.1709106481844573

```

"""
function param_b(variance::Real, skew::Real, E::Real)

    b2 = 2.0 * variance * (1.0 - (E^2 / 2.0) - (skew^2 * ((1.0 - E^2)^2 / (8.0 * E^2))))
    if b2 < 0.0
        @warn "Constraint b > 0 not satisfied, setting b to NaN"
        return NaN
    else
        return sqrt(b2)
    end

end

"""

    param_b(variance::AbstractArray, skew::AbstractArray, E::AbstractArray)

Solve for the parameter b in Equation 8(c) of Sardeshmukh et al. (2015).

# Examples
```jldoctest
julia> E = 0.6236095644623235
0.6236095644623235

julia> param_b(1.0, 1.0, E)
1.1709106481844573

```

"""
function param_b(variance::AbstractArray, skew::AbstractArray, E::AbstractArray)

    b2 .= 2.0 .* variance .* (1.0 .- (E .^2 / 2.0) .- (skew .^2 .* (1.0 .- E .^2) .^2 ./ (8.0 .* E .^2)))
    return b2 .< 0.0 ? NaN : sqrt.(b2)

end


text:
using IntervalArithmetic
using IntervalRootFinding

##
# Available strategies
##

f(x) = sin(x)
contractor = Newton(f, cos)

@info "Bread first search."
# A search takes a starting region, a contractor and a tolerance as argument
search = BreadthFirstSearch(-10..10, contractor, 1e-10)

# This is needed to avoid the data being discarded at the end of the loop use
# `local tree` instead inside a function.
global endtree
for (k, tree) in enumerate(search)
    println("Tree at iteration $k")
    println(tree)  # The tree has custom printing
    global endtree = tree
end

rts = data(endtree)  # Use `data` to only get the leaves marked as `:final`
println("Final $(length(rts)) roots: $rts")
println()

@info "Depth first search."
# The second available type of root search is depth first
search = DepthFirstSearch(-10..10, contractor, 1e-10)
for tree in search
    global endtree = tree
end  # Go to the end of the iteration

# Since the index of each node is printed and indices are attributed in the
# order nodes are processed, comparing this tree with the last iteration of the
# preceding show the difference between breadth first and depth first searches.
println(endtree)

##
# Custom strategy
##

# Define custom BBSearch that store empty intervals rather than discarding it
# Functions of the interface must be explicitely imported
import IntervalRootFinding

struct MySearch{R <: Region, C <: Contractor, T <: Real} <: BreadthFirstBBSearch{Root{R}}
    initial::Root{R}
    contractor::C
    tol::T
end

# This must be implemented, other functions needed are already implemented for
# BBSearch using Root as data (bisect function) or the general fallback match
# this case (root_element function)
function IntervalRootFinding.process(search::MySearch, r::Root)
    contracted_root = search.contractor(r, search.tol)
    status = root_status(contracted_root)

    # We only store the unrefined intervals so the final intervals cover all
    # the initial region
    unrefined_root = Root(r.interval, status)
    status == :unique && return :store, unrefined_root
    status == :empty && return :store, unrefined_root # Store largest known empty intervals
    status == :unknown && diam(contracted_root) < search.tol && return :store, unrefined_root
    return :bisect, unrefined_root  # Always bisect the original interval to bypass [NaN, NaN]
end

@info "Search with no interval discarded."

search = MySearch(Root(-10..10, :unknown), contractor, 1e-10)
for tree in search
    global endtree = tree
end  # Go to the end of the iteration

# In these tree we can see that the union of the resulting intervals (including
# those containing no solution) covers the starting interval.
println(endtree)

ext: jl


text:
using Random

struct CyclicGroup <: GroupsCore.Group
    order::UInt
end

struct CyclicGroupElement <: GroupsCore.GroupElement
    residual::UInt
    parent::CyclicGroup
    CyclicGroupElement(n::Integer, C::CyclicGroup) = new(n % C.order, C)
end

Base.one(C::CyclicGroup) = CyclicGroupElement(0, C)

Base.eltype(::Type{CyclicGroup}) = CyclicGroupElement
Base.iterate(C::CyclicGroup) = one(C), 1
Base.iterate(C::CyclicGroup, state) =
    (state < C.order ? (CyclicGroupElement(state, C), state + 1) : nothing)
Base.IteratorSize(::Type{CyclicGroup}) = Base.HasLength()

GroupsCore.order(::Type{T}, C::CyclicGroup) where {T<:Integer} = T(C.order)
GroupsCore.gens(C::CyclicGroup) = [CyclicGroupElement(1, C)]

function Base.rand(
    rng::Random.AbstractRNG,
    rs::Random.SamplerTrivial{<:CyclicGroup},
)
    C = rs[]
    return CyclicGroupElement(rand(0:C.order-1), C)
end

GroupsCore.parent(c::CyclicGroupElement) = c.parent
GroupsCore.parent_type(::Type{CyclicGroupElement}) = CyclicGroup
Base.:(==)(g::CyclicGroupElement, h::CyclicGroupElement) =
    parent(g) === parent(h) && g.residual == h.residual

# since CyclicGroupElement is isbits, there is no need to define for julia>=1.3.0
if VERSION < v"1.3.0"
    Base.deepcopy_internal(g::CyclicGroupElement, ::IdDict) =
        CyclicGroupElement(deepcopy(g.residual), parent(g))
end

Base.inv(g::CyclicGroupElement) =
    (C = parent(g); CyclicGroupElement(order(UInt, C) - g.residual, C))

function Base.:(*)(g::CyclicGroupElement, h::CyclicGroupElement)
    @assert parent(g) === parent(h)
    C = parent(g)
    return CyclicGroupElement(g.residual + h.residual, C)
end

################## Implementing Group Interface Done!

#=
### Possible performance modifications:

NOTE: Since CyclicGroupElement is immutable there is no need to implement in-place mutable arithmetic.

Base.isone(g::CyclicGroupElement) = iszero(g.residual)

function GroupsCore.order(::Type{I}, g::CyclicGroupElement) where {I<:Integer}
    isone(g) && return one(I)
    o = order(I, parent(g))
    return div(o, gcd(I(g.residual), o))
end

Base.hash(g::CyclicGroupElement, h::UInt) = hash(g.residual, hash(parent(g), h))

=#

### end of Group[Element] methods

# Some eye-candy if you please
Base.show(io::IO, C::CyclicGroup) =
    print(io, "Group of residues modulo $(order(Int, C))")
Base.show(io::IO, c::CyclicGroupElement) =
    print(io, Int(c.residual), " (mod ", order(Int, parent(c)), ")")


text:
"""
Dense layer
"""
struct Dense <: StaticLayer
    w
    b
    f::Function
end

(d::Dense)(x) = d.f.(d.w * Knet.mat(x) .+ d.b)

"""
Create dense layer by input and output size.
"""
function Dense(i_dim::Integer, o_dim::Integer; f::Function=identity, winit=Knet.xavier, b0::Bool=true)
    w = Knet.param(o_dim, i_dim; atype=AT{FT,2}, init=winit)
    b_init = b0 ? Knet.param0 : Knet.param
    b = b_init(o_dim; atype=AT{FT,1})
    return Dense(w, b, f)
end

"""
Lazy dense layer ...zZzZzZ
"""
struct LazyDense <: StaticLayer
    w
    b
    f::Function
    b0::Bool
end

function (d::LazyDense)(x, k::Int=length(d.b); b0::Bool=d.b0)
    # TODO: use `Knet.mat` here as well
    sx1 = size(x, 1)
    if sx1 > size(d.w, 2)
        d.w.value = hcat([d.w.value, Knet.param(size(d.w, 1), sx1 - size(d.w, 2); atype=AT{FT,2}).value]...)
    end
    if k > size(d.w, 1)
        k_diff = k - size(d.w, 1)
        d.w.value = vcat([d.w.value, Knet.param(k_diff, size(d.w, 2); atype=AT{FT,2}).value]...)
        d.b.value = vcat([d.b.value, (b0 ? Knet.param0 : Knet.param)(k_diff; atype=AT{FT,1}).value]...)
    end
    return d.f.(d.w[1:k,1:sx1] * x .+ d.b[1:k])
end

function LazyDense(i_dim::Integer, o_init::Integer; f::Function=identity, b0::Bool=true)
    d = Dense(i_dim, o_init; f=f, b0=b0)
    return LazyDense(d.w, d.b, d.f, b0)
end

"""
Dynamic output using RNN
"""
struct DynamicOut <: StaticLayer
    rnn::Knet.RNN
    mlp::StaticLayer
end

function DynamicOut(i_dim::Integer, r_dim::Integer;
                    skipInput=false, rnnType=:relu, bidirectional=false, f=identity)
    rnn = Knet.RNN(i_dim, r_dim; skipInput=skipInput, rnnType=rnnType, bidirectional=bidirectional, dataType=FT)
    mlp = Dense(r_dim * (bidirectional ? 2 : 1), 1; f=f)
    return DynamicOut(rnn, mlp)
end

function (dy::DynamicOut)(x, d::Integer=dy.r_dim)
    (x_dim, batch_size) = size(x)
    h = dy.rnn(reshape(hcat([x for _ = 1:d]...), x_dim, batch_size, d))
    # The `reshape` below was double-checked - don't waste time on debugging it
    y = dy.mlp(reshape(h, size(h, 1), batch_size * d))
    return reshape(y, batch_size, d)'
end

struct DynamicIn <: StaticLayer
    rnn::Knet.RNN
    mlp::StaticLayer
end

"""
Dynamic input using RNN
"""
function DynamicIn(r_dim::Integer, o_dim::Integer;
                   skipInput=false, rnnType=:relu, bidirectional=false, f=identity)
    rnn = Knet.RNN(1, r_dim; skipInput=skipInput, rnnType=rnnType, bidirectional=bidirectional, dataType=FT)
    mlp = Dense(r_dim * (bidirectional ? 2 : 1), o_dim; f=f)
    return DynamicIn(rnn, mlp)
end

function (dy::DynamicIn)(x)
    (x_dim, batch_size) = size(x)
    y = dy.rnn(reshape(x', 1, batch_size, x_dim))
    # Only the last hidden state is used
    h = if dy.rnn.direction == 0
        y[:,:,end]
    else
        h_dim = dy.rnn.hiddenSize
        vcat(y[:,:,end][1:h_dim,:], y[:,:,1][h_dim+1:end,:])
    end
    return dy.mlp(h)
end

"""
Chaining multiple layers.

NOTE: only chainning layers are allowed but not models. As models are assumed to output loss when being called.
"""
struct Chain <: AbstractLayer
    layers::Tuple
    function Chain(layers::Tuple)
        n = length(layers)
        for i = 1:n-1
            @assert layers[i] isa StaticLayer "The layers in middle should be `StaticLayer`."
        end
        @assert layers[n] isa AbstractLayer "The last layer should be a `AbstractLayer`"
        return new(layers)
    end
end
Chain(layers::AbstractLayer...) = Chain(layers)
Chain(layers) = Chain(layers...)

"""
Run chained layers.
"""
function (c::Chain)(x)
    for l in c.layers
        x = l(x)
    end
    return x
end

"""
Run chained layers with `args...` applied to the last one.
"""
function (c::Chain)(x, args...)
    n = length(c.layers)
    for i = 1:n-1
        x = c.layers[i](x)
    end
    return c.layers[n](x, args...)
end


text:
function f(x)
    return x^2 - 2
end

function f1(x)
    return 2x
end

method = "Newton"
method = "Steffensen"
niter  = 10
x      = 1.0
fv     = zeros(niter)
for iter = 1:niter
    println( f(x) )
    fv[iter] = abs( f(x) )
    if method == "Newton"
        x -= f(x)/f1(x)
    elseif method == "Steffensen"
        h  =  f(x)
        g  =  (f(x + h) - f(x)) / h
        x -= f(x)/g
    end
end

display( plot(1:niter, log10.(fv), markershape=:cross) )


text:
mutable struct LUFactorization{Tv,Ti} <: AbstractLUFactorization{Tv,Ti}
    A::Union{Nothing,ExtendableSparseMatrix{Tv,Ti}}
    fact::Union{Nothing,SuiteSparse.UMFPACK.UmfpackLU{Tv,Ti}}
    phash::UInt64
end

"""
```
LUFactorization(;valuetype=Float64, indextype=Int64)
LUFactorization(matrix)
```

Default Julia LU Factorization based on umfpack.
"""
LUFactorization(;valuetype::Type=Float64,indextype::Type=Int64)=LUFactorization{valuetype,indextype}(nothing,nothing,0)

function update!(lufact::LUFactorization)
    flush!(lufact.A)
    if lufact.A.phash!=lufact.phash
        lufact.fact=lu(lufact.A.cscmatrix)
        lufact.phash=lufact.A.phash
    else
        lufact.fact=lu!(lufact.fact,lufact.A.cscmatrix)
    end
    lufact
end


text:
export margin

function margin(G::TransferFunction; plot = true, tol = sqrt(eps()))
    # Gain Margin
    # ===========

    # Compute the phase crossover frequency, which is used to compute the gain
    # margin.
    jw  = Poly([0;im])
    p   = G.num(jw)*conj(G.den(jw))
    imp = Poly(imag(p.a))
    wv  = roots(imp)

    # Search for the real, positive frequency that yields the lowest gain
    # margin.
    wcg = NaN
    Gm  = Inf

    for w in wv
        # Discard imaginary solutions.
        (abs(imag(w)) > tol) && continue

        rw = real(w)

        # Discard negative solutions.
        (rw < 0.0) && continue

        Gw = G.num(rw*im)/G.den(rw*im)

        # The real part of the transfer function at `w` must be negative.
        ( real(Gw) > 0.0 ) && continue

        # Do not consider undefined solutions.
        (isnan(Gw)) && continue

        # Compute the new gain margin candidate.
        Gm_candidate = -real(G.den(rw*im)/G.num(rw*im))

        # Check if it is lower than the last one.
        (Gm_candidate < Gm) && (Gm = Gm_candidate; wcg = rw)
    end

    # Phase Margin
    # ============

    # Compute the gain crossover frequency, which is used to compute the phase
    # margin.
    p = polyval(G.num,jw)*conj(polyval(G.num,jw)) -
        polyval(G.den,jw)*conj(polyval(G.den,jw))

    wv = roots(p)

    # Search for the real, positive frequency that yields the lowest phase
    # margin.
    wcp  = NaN
    Pm   = Inf

    for w in wv
        # Discard imaginary solutions.
        (abs(imag(w)) > tol) && continue

        rw = real(w)

        # Discard negative solutions.
        (rw < 0.0) && continue

        # Compute the new phase margin candidate.
        Pm_candidate = 180.0 + angle(G.num(rw*im)/G.den(rw*im))*180.0/pi

        # If it is greater than 180 deg, then we must subtract 360 deg.
        (Pm_candidate > 180.0) && (Pm_candidate -= 360.0)

        # Check if it is lower than the last one.
        (Pm_candidate < Pm) && (Pm = Pm_candidate; wcp = rw)
    end

    if plot
        # Compute bode diagram and plot.
        gain, phase, w = bode(G; plot = false)
        plot_margin(gain, phase, w, Gm, Pm, wcg, wcp)
    end

    Gm, Pm, wcg, wcp
end


text:
module Visualization

using PyPlot
using Images

using DeformableRegistration: Transformation, Interpolation, ImageProcessing

export showImage, plotGrid, visualizeResults

# plotting properties
PyPlot.rc("legend",fontsize=20)
PyPlot.rc("xtick", labelsize=10, color="black", direction="out")
PyPlot.rc("ytick", labelsize=10, color="black", direction="out")

function showImage(image::regImage;name="")
    extent = [image.shift[2],
    image.shift[2] + image.voxelsize[2]*size(image.data,2),
    image.shift[1],
    image.shift[1] + image.voxelsize[1]*size(image.data,1)]
    imshow(image.data,
           interpolation="none",
           cmap="gray",
           extent=extent)
    xlabel("x",fontsize=10)
    ylabel("y",fontsize=10)
    title(name,fontsize=12)
end

function plotGrid(transformedGrid::Array{Float64,1}, shape::Tuple{Int64,Int64}; numberOfGridLines=20)
    tX = reshape(transformedGrid[1:Int(length(transformedGrid)/2)],shape)
    tY = reshape(transformedGrid[Int(length(transformedGrid)/2+1):end],shape)

    for i = convert(Array{Int}, round.(linspace(1,size(tX,2), numberOfGridLines)))
        plot(tY[:,i], tX[:,i], "k", linewidth=1.0)
    end
    for i = convert(Array{Int}, round.(linspace(1,size(tX,1), numberOfGridLines)))
        plot(tY[i,:], tX[i,:],"k", linewidth=1.0)
    end
end

function plotGrid(g::scaledArray; numberOfGridLines=20)
    plotGrid(g.data, g.dimensions, numberOfGridLines=numberOfGridLines)
end

function visualizeResults(referenceImage,templateImage;
                          displacement=scaledArray(
                            zeros(2*prod(size(referenceImage.data))),
                            size(referenceImage.data),
                            referenceImage.voxelsize,
                            referenceImage.shift),
                          affineParameters=[1,0,0,0,1,0.0],
                          numberOfGridLines = 20)

    subplot(2,3,1)
    showImage(referenceImage)
    title("Reference")

    subplot(2,3,2)
    showImage(templateImage)
    title("Template")

    subplot(2,3,3)
    showImage(createImage(abs.(Array(referenceImage.data)-Array(templateImage.data))))
    title("Template-Reference")

    subplot(2,3,4)
    centeredGrid = getCellCenteredGrid(referenceImage)
    transformedGrid = transformGridAffine(centeredGrid,affineParameters) + displacement
    plotGrid(transformedGrid, numberOfGridLines=numberOfGridLines)
    title("y (def. Grid)")

    subplot(2,3,5)
    transformedTemplate = interpolateImage(templateImage,transformedGrid)
    showImage(createImage(transformedTemplate))
    title("Template[y]")

    subplot(2,3,6)
    differenceImage = abs.(Array(referenceImage.data)-Array(transformedTemplate))
    showImage(createImage(differenceImage))
    xlabel("max: $(maximum(differenceImage[:]))")
    title("Template[y]-Reference")

end

end


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

"""
    GitIndex(repo::GitRepo)

Load the index file for the repository `repo`.
"""
function GitIndex(repo::GitRepo)
    idx_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_repository_index, :libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}), idx_ptr_ptr, repo.ptr)
    return GitIndex(repo, idx_ptr_ptr[])
end

"""
    read!(idx::GitIndex, force::Bool = false) -> GitIndex

Update the contents of `idx` by reading changes made on disk. For example, `idx`
might be updated if a file has been added to the repository since it was created.
If `force` is `true`, any changes in memory (any changes in `idx` since its last
[`write!`](@ref), or since its creation if no writes have occurred) are discarded.
If `force` is `false`, the index data is only updated from disk if the data on disk
has changed since the last time it was loaded into `idx`.
"""
function read!(idx::GitIndex, force::Bool = false)
    @check ccall((:git_index_read, :libgit2), Cint, (Ptr{Void}, Cint), idx.ptr, Cint(force))
    return idx
end

"""
    write!(idx::GitIndex) -> GitIndex

Write the state of index `idx` to disk using a file lock.
"""
function write!(idx::GitIndex)
    @check ccall((:git_index_write, :libgit2), Cint, (Ptr{Void},), idx.ptr)
    return idx
end

"""
    write_tree!(idx::GitIndex) -> GitHash

Write the index `idx` as a [`GitTree`](@ref) on disk. Trees will be recursively
created for each subtree in `idx`. The returned [`GitHash`](@ref) can be used to
create a [`GitCommit`](@ref). `idx` must have a parent repository and this
repository cannot be bare. `idx` must not contain any files with conflicts.
"""
function write_tree!(idx::GitIndex)
    oid_ptr = Ref(GitHash())
    @check ccall((:git_index_write_tree, :libgit2), Cint,
                 (Ptr{GitHash}, Ptr{Void}), oid_ptr, idx.ptr)
    return oid_ptr[]
end

function repository(idx::GitIndex)
    if isnull(idx.owner)
        throw(GitError(Error.Index, Error.ENOTFOUND, "Index does not have an owning repository."))
    else
        return Base.get(idx.owner)
    end
end

"""
    LibGit2.read_tree!(idx::GitIndex, tree::GitTree)
    LibGit2.read_tree!(idx::GitIndex, treehash::AbstractGitHash)

Read the tree `tree` (or the tree pointed to by `treehash` in the repository owned by
`idx`) into the index `idx`. The current index contents will be replaced.
"""
function read_tree!(idx::GitIndex, tree::GitTree)
    @check ccall((:git_index_read_tree, :libgit2), Cint,
                 (Ptr{Void}, Ptr{Void}), idx.ptr, tree.ptr)
end
read_tree!(idx::GitIndex, hash::AbstractGitHash) =
    read_tree!(idx, GitTree(repository(idx), hash))

"""
    add!(repo::GitRepo, files::AbstractString...; flags::Cuint = Consts.INDEX_ADD_DEFAULT)
    add!(idx::GitIndex, files::AbstractString...; flags::Cuint = Consts.INDEX_ADD_DEFAULT)

Add all the files with paths specified by `files` to the index `idx` (or the index
of the `repo`). If the file already exists, the index entry will be updated.
If the file does not exist already, it will be newly added into the index.
`files` may contain glob patterns which will be expanded and any matching files will
be added (unless `INDEX_ADD_DISABLE_PATHSPEC_MATCH` is set, see below).
If a file has been ignored (in `.gitignore` or in the config), it *will not* be
added, *unless* it is already being tracked in the index, in which case it *will* be
updated. The keyword argument `flags` is a set of bit-flags which control the behavior
with respect to ignored files:
  * `Consts.INDEX_ADD_DEFAULT` - default, described above.
  * `Consts.INDEX_ADD_FORCE` - disregard the existing ignore rules and force addition of
    the file to the index even if it is already ignored.
  * `Consts.INDEX_ADD_CHECK_PATHSPEC` - cannot be used at the same time as `INDEX_ADD_FORCE`.
    Check that each file in `files` which exists on disk is not in the ignore list. If one
    of the files *is* ignored, the function will return `EINVALIDSPEC`.
  * `Consts.INDEX_ADD_DISABLE_PATHSPEC_MATCH` - turn off glob matching, and only add files
    to the index which exactly match the paths specified in `files`.
"""
function add!(idx::GitIndex, files::AbstractString...;
              flags::Cuint = Consts.INDEX_ADD_DEFAULT)
    @check ccall((:git_index_add_all, :libgit2), Cint,
                 (Ptr{Void}, Ptr{StrArrayStruct}, Cuint, Ptr{Void}, Ptr{Void}),
                 idx.ptr, collect(files), flags, C_NULL, C_NULL)
end

"""
    update!(repo::GitRepo, files::AbstractString...)
    update!(idx::GitIndex, files::AbstractString...)

Update all the files with paths specified by `files` in the index `idx` (or the index
of the `repo`). Match the state of each file in the index with the current state on
disk, removing it if it has been removed on disk, or updating its entry in the object
database.
"""
function update!(idx::GitIndex, files::AbstractString...)
    @check ccall((:git_index_update_all, :libgit2), Cint,
                 (Ptr{Void}, Ptr{StrArrayStruct}, Ptr{Void}, Ptr{Void}),
                 idx.ptr, collect(files), C_NULL, C_NULL)
end

"""
    remove!(repo::GitRepo, files::AbstractString...)
    remove!(idx::GitIndex, files::AbstractString...)

Remove all the files with paths specified by `files` in the index `idx` (or the index
of the `repo`).
"""
function remove!(idx::GitIndex, files::AbstractString...)
    @check ccall((:git_index_remove_all, :libgit2), Cint,
                 (Ptr{Void}, Ptr{StrArrayStruct}, Ptr{Void}, Ptr{Void}),
                 idx.ptr, collect(files), C_NULL, C_NULL)
end

function add!(repo::GitRepo, files::AbstractString...;
              flags::Cuint = Consts.INDEX_ADD_DEFAULT)
    with(GitIndex, repo) do idx
        add!(idx, files..., flags = flags)
        write!(idx)
    end
    return
end

function update!(repo::GitRepo, files::AbstractString...)
    with(GitIndex, repo) do idx
        update!(idx, files...)
        write!(idx)
    end
    return
end

function remove!(repo::GitRepo, files::AbstractString...)
    with(GitIndex, repo) do idx
        remove!(idx, files...)
        write!(idx)
    end
    return
end

function read!(repo::GitRepo, force::Bool = false)
    with(GitIndex, repo) do idx
        read!(idx, force)
    end
    return
end

function Base.count(idx::GitIndex)
    return ccall((:git_index_entrycount, :libgit2), Csize_t, (Ptr{Void},), idx.ptr)
end

function Base.getindex(idx::GitIndex, i::Integer)
    Base.@gc_preserve idx begin
        ie_ptr = ccall((:git_index_get_byindex, :libgit2),
                       Ptr{IndexEntry},
                       (Ptr{Void}, Csize_t), idx.ptr, i-1)
        ie_ptr == C_NULL && return nothing
        elem = unsafe_load(ie_ptr)
    end
    return elem
end

function Base.find(path::String, idx::GitIndex)
    pos_ref = Ref{Csize_t}(0)
    ret = ccall((:git_index_find, :libgit2), Cint,
                  (Ref{Csize_t}, Ptr{Void}, Cstring), pos_ref, idx.ptr, path)
    ret == Cint(Error.ENOTFOUND) && return Nullable{Csize_t}()
    return Nullable(pos_ref[]+1)
end

"""
    stage(ie::IndexEntry) -> Cint

Get the stage number of `ie`. The stage number `0` represents the current state
of the working tree, but other numbers can be used in the case of a merge conflict.
In such a case, the various stage numbers on an `IndexEntry` describe which side(s)
of the conflict the current state of the file belongs to. Stage `0` is the state
before the attempted merge, stage `1` is the changes which have been made locally,
stages `2` and larger are for changes from other branches (for instance, in the case
of a multi-branch "octopus" merge, stages `2`, `3`, and `4` might be used).
"""
stage(ie::IndexEntry) = ccall((:git_index_entry_stage, :libgit2), Cint, (Ptr{IndexEntry},), Ref(ie))

function Base.show(io::IO, idx::GitIndex)
    println(io, "GitIndex:\nRepository: ", repository(idx), "\nNumber of elements: ", count(idx))
end

max_stars_repo_name: ihnorton/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
using Kash
using MinHash
using Test
using FASTX
using BioSequences

TESTPATH = joinpath(dirname(@__FILE__), "data", "test.fna")

## FastaIterator
@testset "FastaIterator" begin
    function test_fastaiterator(one::FastaIterator{T}, two) where T
        # FastaIterators mutate their output, so must copy
        A = [copy(i) for i in one]
        B = [FASTA.sequence(LongSequence{T}, i) for i in two]
        @test length(A) == length(B)
        @test all(i == j for (i, j) in zip(A, B))
        close(one.reader)
        close(two)
    end

    path = TESTPATH
    for T in [DNAAlphabet{2}, DNAAlphabet{4}]
        # Instantiate from FASTA Reader
        reader_one = FastaIterator{T}(FASTA.Reader(open(path)))
        reader_two = FASTA.Reader(open(path))
        test_fastaiterator(reader_one, reader_two)

        # Instantiate from IO
        reader_one = FastaIterator{T}(open(path))
        reader_two = FASTA.Reader(open(path))
        test_fastaiterator(reader_one, reader_two)

        # Instantiate from other IO
        io = open(path) do file
            IOBuffer(read(file))
        end
        reader_one = FastaIterator{T}(io)
        reader_two = FASTA.Reader(open(path))
        test_fastaiterator(reader_one, reader_two)

        # Instantiate from IO
        reader_one = FastaIterator{T}(path)
        reader_two = FASTA.Reader(open(path))
        test_fastaiterator(reader_one, reader_two)
    end
end

## CanonicalKmerIterator
@testset "CanonicalKmerIterator" begin
    function test_can_kmerit(it1::CanonicalKmerIterator{M}, it2::M) where M
        A = collect(it1)
        B = [canonical(i) for i in it2]
        @test A == B
    end

    it2 = each(DNAMer{11}, dna"TAGTAGGCGCGCGGCCGTATAGATGATGCTAGAAAGGC")
    it1 = CanonicalKmerIterator(it2)
    test_can_kmerit(it1, it2)

    it2 = each(DNAMer{3}, dna"TAGTAGNGCGCGCGGCNCGTATAGATGATGCTAGAAAGNGC")
    it1 = CanonicalKmerIterator(it2)
    test_can_kmerit(it1, it2)

    it2 = each(DNAMer{4}, dna"TAGTCGAAGTGCTGAGAGATCTCTAAAGAGAGCGCTCTGAAAAA", 3, 5)
    it1 = CanonicalKmerIterator(it2)
    test_can_kmerit(it1, it2)
end

## Kmersketcher
@testset "KmerSketcher" begin
    @testset "Instantiation" begin
        h = MinHasher{identity}(100)
        @test isa(KmerSketcher{DNAMer{4}}(h), KmerSketcher{DNAMer{4}, identity})

        @test_throws MethodError KmerSketcher{DNAMer{3},hash}(h)
        @test_throws MethodError KmerSketcher{4,identity}(h)

        sk = KmerSketcher{DNAMer{5}}(100)
        @test isa(sk.hasher, MinHasher{hash})
    end

    @testset "empty!" begin
        sk = KmerSketcher{DNAMer{3}}(10)
        seq = dna"TAGGCGTAGTGCGTATATAGCGAAAGAGCTCTA"
        Kash.update!(sk, seq)

        @test sk.bases == length(seq)
        @test sk.hasher.filled ≤ length(seq) - 2
        empty!(sk)
        @test sk.bases == 0
        @test sk.hasher.filled == 0
    end

    @testset "update! seq" begin
        function test_update(sk::KmerSketcher{M,F}, seq::BioSequence) where {M,F}
            hashes_before = KmerHashes(sk)
            Kash.update!(sk, seq)
            hashes_after = KmerHashes(sk)

            # Test bases is correct
            @test hashes_before.bases + length(seq) == hashes_after.bases

            # Test content of hashes
            new_hashes = [F(canonical(i)) for i in each(M, seq)]
            should_be = sort!(collect(Set(append!(hashes_before.sketch.hashes, new_hashes))))

            # at most max hashes
            new_length = min(length(should_be), length(sk.hasher.heap))

            should_be = should_be[1:new_length]
            @test hashes_after.sketch.hashes == should_be
        end

        sk = KmerSketcher{DNAMer{3}}(10)
        seq = dna"TAGGCGTAGTNGCGTATATAGCGAAAGAGCTCNTA"
        test_update(sk, seq)

        empty!(sk)
        Kash.update!(sk, randdnaseq(10))
        test_update(sk, randdnaseq(100))

        empty!(sk)
        test_update(sk, dna"")
    end
end

## High-level API
@testset "High-level API" begin
    function manual_kmerhashes(io::IO, A, M::Type{<:Mer}, N::Integer)
        hashes = Set{UInt}()
        reader = FASTA.Reader(io)
        for entry in reader
            seq = FASTA.sequence(LongSequence{A}, entry)
            for kmer in each(M, seq)
                push!(hashes, hash(canonical(kmer)))
            end
        end
        return sort!(collect(hashes))[1:N]
    end

    manual = open(TESTPATH) do file
        manual_kmerhashes(file, DNAAlphabet{4}, DNAMer{3}, 20)
    end
    kash = open(TESTPATH) do file
        kmer_minhash(file, 20, Val(3)).sketch.hashes
    end
    @test manual == kash

    function manual_kmerhashes_each(io::IO, A, M::Type{<:Mer}, N::Integer)
        hashvectors = []
        reader = FASTA.Reader(io)
        for entry in reader
            seq = FASTA.sequence(LongSequence{A}, entry)
            hashes = Set{UInt}()
            for kmer in each(M, seq)
                push!(hashes, hash(canonical(kmer)))
            end
            push!(hashvectors, sort!(collect(hashes))[1:min(N, length(hashes))])
        end
        return hashvectors
    end

    manual = open(TESTPATH) do file
        manual_kmerhashes_each(file, DNAAlphabet{4}, DNAMer{3}, 20)
    end
    kash = open(TESTPATH) do file
        [x.sketch.hashes for x in kmer_minhash_each(file, 20, Val(3))]
    end
    @test manual == kash

end

# TODO:
# Serialization

# Then add an AA and a IUPAC DNA file
# Run tests with all 3 alphabets
# Then check codecov
# Add documentation


text:
struct QQPair
	qx::Vector{Float64}
	qy::Vector{Float64}
end

function qqbuild(x::Vector, y::Vector)
	n = min(length(x), length(y))
	grid = [0.0:(1 / (n - 1)):1.0;]
	qx = quantile(x, grid)
	qy = quantile(y, grid)
	return QQPair(qx, qy)
end

function qqbuild(x::Vector, d::UnivariateDistribution)
	n = length(x)
	grid = [(1 / (n - 1)):(1 / (n - 1)):(1.0 - (1 / (n - 1)));]
	qx = quantile(x, grid)
	qd = quantile.(Ref(d), grid)
	return QQPair(qx, qd)
end

function qqbuild(d::UnivariateDistribution, x::Vector)
	n = length(x)
	grid = [(1 / (n - 1)):(1 / (n - 1)):(1.0 - (1 / (n - 1)));]
	qd = quantile.(Ref(d), grid)
	qx = quantile(x, grid)
	return QQPair(qd, qx)
end


text:
abstract type QGate end

function QGate(kernel:: GateKernel)
	if paramCount(kernel)==0
		return ConstGate(kernel)
	else
		return VarGate(kernel)
	end
end

#  visible to state
ITensor(gate::QGate, inds::IndexSet) = ITensor(data(gate), inds)
ITensor(gate::QGate, ind::Index...) = ITensor(data(gate), ind...)


text:
"""
`calc_depth_and_expected_rv_precission(spectrum, pixels_index, order_index; smooth_factor)`
Calculate expected RV uncertainty for one portion of spectrum given by pixels and order indicies.
Assumes only photon noise with given variances.
Returns (depth, `exp_σ_rv`)
"""
function calc_depth_and_expected_rv_precission(spectrum::ST, pixels::AR, order::Integer; smooth_factor::Real = 2.0) where { ST<:AbstractSpectra2D, AR<:AbstractRange{Int64}, AA1<:AbstractArray{Int64,1} }
   logλ = log.(view(spectrum.λ,pixels,order))
   flux = view(spectrum.flux,pixels,order)
   var = view(spectrum.var,pixels,order)
   (f_mean, f_deriv) = RvSpectML.TemporalGPInterpolation.predict_mean_and_deriv(logλ, flux, logλ;sigmasq_obs=var, use_logx=false, use_logy=false, smooth_factor=smooth_factor)
   depth = 1.0 .- minimum(f_mean)/maximum(f_mean)
   exp_sigma_rv = RvSpectMLBase.speed_of_light_mps / sqrt( sum(f_deriv.^2 ./ var) )
   return (depth=depth, exp_σ_rv=exp_sigma_rv)
 end

""" `calc_formal_rv_precission(spectrum, chunklist; smooth_factor)`
Calculate expected RV uncertainty for one spectrum and corresponding chunklist,
Assumes only photon noise with given variances.
"""
function calc_formal_rv_precission(spectrum::ST, chunklist::ACL; smooth_factor::Real = 2.0)  where {  ST<:AbstractSpectra2D, ACL<:AbstractChunkList }
  σ_rv_combo = 0
  for chid in 1:length(chunklist)
   (depth, σ_rv) = calc_depth_and_expected_rv_precission(spectrum,chunklist[chid].flux.indices[1],chunklist[chid].flux.indices[2], smooth_factor=smooth_factor)
   σ_rv_combo += 1/σ_rv^2
  end
  σ_rv_combo = sqrt(1/σ_rv_combo)
end

""" `calc_formal_rv_precission(spectra, chunklist_timeseries; smooth_factor)`
Calculate expected RV uncertainty for each spectrum in an array of spectra and corresponding chunklist_timeseries,
Assumies only photon noise with given variances.
"""
function calc_formal_rv_precission(spectra::AS, chunk_list_timeseries::ACLT; smooth_factor::Real = 2.0, verbose::Bool = false) where { ST<:AbstractSpectra, AS<:AbstractArray{ST,1}, ACLT<:AbstractChunkListTimeseries }
  @assert length(spectra) == length(chunk_list_timeseries)
  map(obs->calc_formal_rv_precission(spectra[obs], chunk_list_timeseries[obs], smooth_factor=smooth_factor), 1:length(spectra))
end


text:
module TenXData

using CodecZlib,TranscodingStreams
using DataFrames,CSV
using SparseArrays # needed for mmread
using SparseArrays # needed for mmread
using LinearAlgebra # needed for mmread
using MatrixMarket: _parseint,find_splits,skewsymmetric!,symmetric!,hermitian! # needed for mmread

import MatrixMarket.mmread

include("MM.jl")

function loadData(countFn::String,geneFn::String,barcodeFn::String)
    extn = split(countFn,".")[end]
    if extn == "gz" || extn == "gzip"
        #fh = GZip.open(countFn)
        fh = GzipDecompressorStream(open(countFn))
        #M = TenX.mmread(fh)
        M = TenX.mmread(fh)
        gh = GzipDecompressorStream(open(geneFn))
        G = map(r->r[1],CSV.File(gh,header=0))
        bh = GzipDecompressorStream(open(barcodeFn))
        B = map(r->r[1],CSV.File(bh,header=0))
    else
        M = mmread(countFn)
        G = map(r->r[1],CSV.File(geneFn,header=0))
        B = map(r->r[1],CSV.File(barcodeFn,header=0))
    end
    df = rename!(DataFrame(M),B)
    insertcols!(df,1,"gene"=>G)
    df
end

# process files with standard names
function loadFolder(dirPath::String)
    fm = Dict(:gn1=>"features.tsv.gz",
              :gn2=>"features.tsv",
              :gn3=>"genes.tsv.gz",
              :gn4=>"genes.tsv",
              :bc=>"barcodes.tsv.gz",
              :bc2=>"barcodes.tsv",
              :cnt=>"matrix.mtx.gz",
              :cnt2=>"matrix.mtx")

    if isfile(joinpath(dirPath,get(fm,:gn1,"")))
        df = TenX.loadData(joinpath(dirPath,get(fm,:cnt,"")),
                           joinpath(dirPath,get(fm,:bc,"")),
                           joinpath(dirPath,get(fm,:gn1,"")))
    elseif isfile(joinpath(dirPath,get(fm,:gn2,"")))
        df = TenX.loadData(joinpath(dirPath,get(fm,:cnt2,"")),
                           joinpath(dirPath,get(fm,:bc2,"")),
                           joinpath(dirPath,get(fm,:gn2,"")))
    elseif isfile(joinpath(dirPath,get(fm,:gn3,"")))
        df = TenX.loadData(joinpath(dirPath,get(fm,:cnt,"")),
                           joinpath(dirPath,get(fm,:bc,"")),
                           joinpath(dirPath,get(fm,:gn3,"")))
    elseif isfile(joinpath(dirPath,get(fm,:gn4,"")))
        df = TenX.loadData(joinpath(dirPath,get(fm,:cnt2,"")),
                           joinpath(dirPath,get(fm,:bc2,"")),
                           joinpath(dirPath,get(fm,:gn4,"")))
    end
    df
end
end


text:
# Heavy tailed proposal
function randh()
	return 10.0^(1.5 - 6.0*rand())*randn()
end

function logsumexp(x::Array{Float64, 1})
	top = maximum(x)
	y = exp(x - top)
	return log(sum(y)) + top
end

function logdiffexp(a::Float64, b::Float64)
	return b + log(exp(a - b) - 1.)
end

max_stars_repo_name: eggplantbren/MogTrack

max_issues_repo_name: eggplantbren/MogTrack

max_forks_repo_name: eggplantbren/MogTrack


text:
function compute()
	ans=1
	for i=1:20
		a=gcd(i,ans)
		ans*=Int(i/a);

	end
	return ans
end

print(compute())


text:
## K-core decomposition ##

function _binsortperm(deg::Vector, n::Integer)
    vert = Array{Int}(n)    # `vert` is the order of vertices
    pos  = Array{Int}(n)    # `pos` is a reverse index to `vert`
    md   = maximum(deg)
    bin  = fill(0, 1 + md)  # `bin` captures offsets for each degree
    @inbounds for i = deg
        bin[1 + i] += 1
    end
    s = 1
    @inbounds for i = 1:(1 + md)
        t = bin[i]
        bin[i] = s
        s += t
    end
    @inbounds for i = 1:n
        pos[i] = bin[1 + deg[i]]
        vert[pos[i]] = i
        bin[1 + deg[i]] += 1
    end
    return vert, pos, bin
end

"""
    kcores(g::Graph)

Compute k-core secomposition of a graph.
Returns a vector with maximum core number for each node.

Implementation based on Vladimir Batagelj and Matjaž Zaveršnik's
"An O(m) Algorithm for Cores Decomposition of Networks"
[link](http://vlado.fmf.uni-lj.si/pub/networks/doc/cores/cores.pdf)
"""
function kcores(g::Graph)
    n   = nodecount(g)
    deg = degree(g)
    vert, pos, bin = _binsortperm(deg, n)
    @inbounds for i = vert
        for j = neighbors(g, i)
            if deg[j] > deg[i]
                dj = deg[j]
                pj = pos[j]
                pk = bin[dj]
                k = vert[pk]
                if k != j
                    pos[j] = pk
                    pos[k] = pj
                    vert[pj] = k
                    vert[pk] = j
                end
                bin[dj] += 1
                deg[j]  -= 1
            end
        end
    end
    return deg
end


text:
using Statistics
using ProgressMeter
using MATLAB
using MAT
using Dates
include("ATN_Coupled_assholes.jl")
@everywhere begin
    using Statistics
    using ProgressMeter
    using MATLAB
    using MAT
    using Dates
    include("ATN_Coupled_assholes.jl")
    N = 10_000          # number of nodes
    m = 5               # number of active links per active node
    gamma = -2.1        # exponent of the power law
    alpha_size = 100
    alpha_min = 0.0
    alpha_max = 1
    n_trials = 10
    Alpha = range(alpha_min, stop = alpha_max, length = alpha_size)
    proportion_size = 100
    proportion_min = 0.4
    proportion_max = 0.9
    Proportion = range(proportion_min, stop = proportion_max, length = proportion_size)
    ActMin = 0.001      # minimum activation probability
    ActMax = 1
    agent_act = ((ActMax.^(gamma + 1) - ActMin.^(gamma + 1)) .* rand(N) .+
                 ActMin.^(gamma + 1)).^(1.0 / (gamma + 1))
    activity = [agent_act]
    N_iters = proportion_size  * alpha_size
    for k in 1:(N_iters - 1)
        push!(activity, agent_act)
    end
    activity = reshape(activity, alpha_size, proportion_size)
end

grid = [(proportion, alpha) for alpha in Alpha, proportion in Proportion]

# Dead = @showprogress pmap(ATN_Coupled_SIRD.parameter_sweep,
#                           n_trials * ones(Int, proportion_size), activity, Proportion)

Dead = @showprogress pmap(ATN_Coupled_SIRD.parameter_sweep,
                          n_trials * ones(Int, alpha_size, proportion_size), activity, grid)

location_mat = "/mnt/storage/MAT/"
right_now = replace(replace(replace(string(Dates.now()), "." => "_"), ":" => "_"),
                    "-" => "_")

file_name = "ADN_AssHats_SISW" * right_now
mat_name = location_mat * file_name * ".mat"
Dead = Float64.(Dead)
matwrite(mat_name, Dict("Dead" => Dead); compress = true)

mat"figure();set(gcf, 'Position',  [0, 0, 2500, 1500]); hold on;"
# mat"axis([ 0,15,0,10 ])"
# mat"axis([ 0,1,0,1 ])"
mat"[X,Y]=meshgrid($Proportion,$Alpha);contourf(X,Y,$Dead,15)"

# mat"colorbar;colormap('inferno')"
# mat"plot($Proportion, $Dead)"


text:
@testset "Assessing local identifiability" begin
    test_cases = []

    # 2-compartiment model

    ode = @ODEmodel(
        x0'(t) = -(a01 + a21) * x0(t) + a12 * x1(t),
        x1'(t) = a21 * x0(t) - a12 * x1(t),
        y(t) = x0(t)
    )
    funcs_to_test = [a01, a21, a12, a01 * a12, a01 + a12 + a21, (a01 + a12 + a21) // (a01 * a12), x0, x1]
    correct = [false, false, false, true, true, true, true, false]
    push!(test_cases, Dict(
        :ode => ode,
        :funcs => funcs_to_test,
        :correct => correct
    ))

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

    ode = @ODEmodel(
        x0'(t) = a * x0(t) - b * x0(t) * x1(t) + u(t),
        x1'(t) = c * x1(t) + d * x0(t) * x1(t),
        y(t) = x0(t)
    )
    funcs_to_test = [a, b, c, d, x0, x1]
    correct = [true, false, true, true, true, false]
    push!(test_cases, Dict(
        :ode => ode,
        :funcs => funcs_to_test,
        :correct => correct
    ))

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

    ode = @ODEmodel(
        S'(t) = mu - bi * S(t) * I(t) - bw * S(t) * W(t) - mu * S(t) + a * R(t),
        I'(t) = bw * S(t) * W(t) + bi * S(t) * I(t) - (gam + mu) * I(t),
        W'(t) = xi * (I(t) - W(t)),
        R'(t) = gam * I(t) - (mu + a) * R(t),
        y(t) = k * I(t)
    )
    funcs_to_test = [mu, bi, bw, a, xi, gam, mu, gam + mu, k, S, I, W, R]
    correct = [true for _ in funcs_to_test]
    push!(test_cases, Dict(
        :ode => ode,
        :funcs => funcs_to_test,
        :correct => correct
    ))

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

    ode = @ODEmodel(
        x1'(t) = -b * x1(t) + 1 / (c + x4(t)),
        x2'(t) = alpha * x1(t) - beta * x2(t),
        x3'(t) = gama * x2(t) - delta * x3(t),
        x4'(t) = sigma * x4(t) * (gama * x2(t) - delta * x3(t)) / x3(t),
        y(t) = x1(t)
    )
    funcs_to_test = [b, c, alpha, beta, delta, gama, beta + delta, beta * delta]
    correct = [true, true, false, true, true, false, true, true]
    push!(test_cases, Dict(
        :ode => ode,
        :funcs => funcs_to_test,
        :correct => correct
    ))

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

    for case in test_cases
        result = assess_local_identifiability(case[:ode], case[:funcs])
        @test result == case[:correct]
    end
end


text:
using GeoStatsBase

"""Interpolate missing values in GeoArray."""
function interpolate!(ga::GeoArray, solver::T, band=1) where T <: EstimationSolver
    # Irregular grid
    # TODO Use unstructured GeoStats method
    if is_rotated(ga)
        error("Can't interpolate warped grid yet.")

    # Regular grid
    else
        data = @view ga.A[:,:,band]
        problemdata = georef(
            (;band=data,),
            origin=Tuple(ga.f.translation),
            spacing=(abs(ga.f.linear[1]), abs(ga.f.linear[4]))
        )
    end
    domain = LinearIndices(size(data))[ismissing.(data)]
    problem = EstimationProblem(problemdata, view(problemdata.domain, domain), :band)
    solution = solve(problem, solver)

    data[domain] .= solution[:band]
    ga
end


text:
module ImagineWorker

using ImagineInterface, Unitful, NIDAQ, Distributed

@static if !Sys.iswindows()
    @warn("This package only works on Windows")
end

@static if Sys.iswindows() && !isfile("C:\\Windows\\System32\\nicaiu.dll")
    error("nicaiu.dll not found. Is NIDAQmx installed?")
end

@static if Sys.iswindows() && isfile("C:\\Windows\\System32\\nicaiu.dll")
    using NIDAQ
    global DEVICE_PREFIX = ""

    export _set_device,
		_record_analog_signals,
		_output_analog_signals

    _set_device(dev::T) where {T<:AbstractString} = global DEVICE_PREFIX = String(dev)

    include("inputs.jl")
    include("outputs.jl")
end

end #module


text:
using Base.Iterators: eachindex
using OffsetArrays: Iterators
# ---
# title: 274. H-Index
# id: problem274
# author: Indigo
# date: 2021-06-05
# difficulty: Medium
# categories: Hash Table, Sort
# link: <https://leetcode.com/problems/h-index/description/>
# hidden: true
# ---
#
# Given an array of citations (each citation is a non-negative integer) of a
# researcher, write a function to compute the researcher's h-index.
#
# According to the [definition of h-index on
# Wikipedia](https://en.wikipedia.org/wiki/H-index): "A scientist has index _h_
# if _h_ of his/her _N_ papers have **at least** _h_ citations each, and the
# other _N − h_ papers have **no more than** _h_ citations each."
#
# **Example:**
#
#
#
#     Input: citations = [3,0,6,1,5]
#     Output: 3
#     Explanation:[3,0,6,1,5] means the researcher has 5 papers in total and each of them had
#                  received 3, 0, 6, 1, 5 citations respectively.
#                  Since the researcher has 3 papers with **at least** 3 citations each and the remaining
#                  two with **no more than** 3 citations each, her h-index is 3.
#
# **Note:  **If there are several possible values for _h_ , the maximum one is
# taken as the h-index.
#
#
## @lc code=start
using LeetCode

using OffsetArrays
function h_index(citations::Vector{Int})::Int
    len = length(citations)
    cnt = OffsetArray(fill(0, len + 1), -1)
    for c in citations
        cnt[min(c, len)] += 1
    end
    c = 0
    for i in Iterators.reverse(eachindex(cnt))
        c += cnt[i]
        c == i && return i
    end
end
## @lc code=end


text:
# Cross Polar Cap Potential (CPCP) calculation at z=2.0 cut plane.
# Combination of surface integral and line integral.
# Part 1: parallel ParaView launching on clusters for stream tracing and point
# interpolation.
#
# Using Channels for parallel tasks.
#
# Hongyang Zhou, hyzhou@umich.edu 03/19/2020

using Distributed, Glob
@everywhere using DelimitedFiles, ParallelDataTransfer

# The first parameter of a channel is the Type (e.g.: Int) and the second one,
# the maximum number of elements allowed in the channel.
const jobs    = RemoteChannel(()->Channel{Int}(32))
const results = RemoteChannel(()->Channel{Tuple}(32))

"Create surface seeds for interpolation."
function createSurfaceSeeds(filename="seeds.txt", dir="./")
   xMin, xMax = -2.0, 1.3 # xrange
   yMin, yMax = -3.0, 3.1 # yrange
   z = 2.0 # z cut plane

   dx, dy = 1/32, 1/32 # resolutions

   xRange = xMin:dx:xMax
   yRange = yMin:dy:yMax

   X = [x for x in xRange, y in yRange if 3.5^2 > (x-1.2)^2 + (y-0.05)^2]
   Y = [y for x in xRange, y in yRange if 3.5^2 > (x-1.2)^2 + (y-0.05)^2]
   Z = fill(z, size(X))

   open(joinpath(dir,filename), "w") do io
      write(io, "\"X\",\"Y\",\"Z\"\n")
      writedlm(io, [X Y Z], ',')
   end
end

@everywhere function findStatus(filename, fileDir, outname)

   fieldline = readdlm(filename, ',', Float32, header=true)
   data = fieldline[1]

   # count the total number of field lines
   nlinetotal = 0
   start = Int32[]
   for i = 1:size(data,1)
      if data[i,4] == 0.0 # IntegrationTime starts at 0.0
         nlinetotal += 1
         push!(start, i)
      end
   end

   line = Array{Float32, 2}(undef, 6, nlinetotal)

   for i = 1:nlinetotal-1
      iStart, iEnd = start[i], start[i+1]-1
      line[1,i] = data[iStart,5]
      line[2,i] = data[iStart,6]
      line[3,i] = data[iStart,7]
      line[4,i] = data[iEnd,5]
      line[5,i] = data[iEnd,6]
      line[6,i] = data[iEnd,7]
   end

   line[1,end] = data[start[end],5]
   line[2,end] = data[start[end],6]
   line[3,end] = data[start[end],7]
   line[4,end] = data[end,5]
   line[5,end] = data[end,6]
   line[6,end] = data[end,7]

   # status
   # (0: open, 1: closed along B, 2: closed along -B 3: fully closed)
   # (-1: cells inside body, -2: loop ray within block, -3: strange)
   status = zeros(Int8, nlinetotal)

   # Cluster seed points based on their end points
   for i = 1:nlinetotal
      if line[4,i]^2 + line[5,i]^2 + line[6,i]^2 < 1.2^2 ||
         (line[4,i] > 0.0 && (-2.0 < line[5,i] < 2.0)) # no holes inside!
         status[i] = 1
      end
   end

   pts_closed = line[1:2, status .== 1]
   x = @view pts_closed[1,:]
   y = @view pts_closed[2,:]
   z = fill(2.0f0, length(x))

   open(fileDir*outname, "w") do io
      write(io, "\"X\",\"Y\",\"Z\"\n")
      writedlm(io, [x y z], ',')
   end
end

@everywhere function process_file(i)

   file = files[i]
   @info "processing file $file ..."
   # Create Python script with the correct input filename for tracing
   txt = read(filePy_orig[1], String)
   tracePy = filePy_orig[1][1:end-3]*string(i)*".py"
   fileTrace = "streamline"*lpad(i,4,"0")*".txt"
   fileStream = joinpath(fileDir, fileTrace)
   fileSurfaceXYZ = "surface"*lpad(i,4,"0")*".txt"
   interpPy = filePy_orig[2][1:end-3]*string(i)*".py"
   fileSurface = "surface_value"*lpad(i,4,"0")*".txt"

   if isfile(fileDir*fileSurface)
      return # skip the finished snapshots
   end

   open(tracePy, "w") do f
      txt = replace(txt, "test.vtu"=>file)
      txt = replace(txt, "streamline.txt"=>fileTrace)
      write(f, txt)
   end

   # Step 2: field line tracing (in Paraview)
   run(`$run_pvbatch $tracePy`)
   rm(tracePy, force=true)

   # Step 3: find open/closed status for each point
   findStatus(fileStream, fileDir, fileSurfaceXYZ)
   rm(fileDir*fileTrace, force=true) # hang if file not found with force=false

   # Create Python script with the correct input filename for interpolation
   txt = read(filePy_orig[2], String)
   open(interpPy, "w") do f
      txt = replace(txt, "test.vtu" => file)
      txt = replace(txt, "surface.txt" => fileSurfaceXYZ)
      txt = replace(txt, "surface_value.txt" => fileSurface)
      write(f, txt)
   end

   # Step 4: interpolate onto the boundary points (in Paraview)
   run(`$run_pvbatch $interpPy`)
   rm(interpPy, force=true)
end

@everywhere function do_work(jobs, results) # Define work function everywhere.
   while true
      job_id = take!(jobs)
      process_file(job_id)
      put!(results, (job_id, myid()))
   end
end

function make_jobs(n)
   for i in 1:n
      put!(jobs, i)
   end
end

## Main
if Sys.isapple() # Sys.MACHINE might be better
   @everywhere run_pvbatch = "pvbatch"
else
   @everywhere run_pvbatch = "run_pvbatch"
end
##
@everywhere dataDir = "../../GM/"
@everywhere fileDir = "data/"

if !isdir(fileDir)
   mkdir(fileDir)
end
# Step 1: create seeding points (same for all timesteps)
createSurfaceSeeds("seeds.txt", fileDir)
# Loop over multiple snapshots
files = basename.(glob("cut*.vtu", dataDir))

sendto(workers(), files=files)

@everywhere filePy_orig = ["streamtrace_save.py", "Interpolate_On_Surface_batch.py"]

nFile = length(files)

@async make_jobs(nFile) # Feed the jobs channel with "n" jobs.

# Start tasks on the workers to process requests in parallel.
for p in workers()
   @async remote_do(do_work, p, jobs, results) # Similar to remotecall.
end

@elapsed while nFile > 0 # Print out results.
   job_id, wid = take!(results)
   println("$job_id finished on worker $wid")
   global nFile -= 1
end

##
# Integration after the parallel work using ParaView is finished.
# See runSurfaceIntegral_part2.jl


text:
## Framework for setting up an iterative problem for finding a zero
## TODO
## * a graphic of trace when verbose=true?

#
# In McNamee & Pan (DOI:10.1016/j.camwa.2011.11.015 there are a number of
# results on efficiencies of a solution, (1/d) log_10(q)
# Those implemented here are:
# quadratic cut (Muller) .265 (in a42)
# Newton() newton = .1505   or 1/2log(2)
# Order1() secant method .20 (1/1 * log(1.6)
# FalsePostion(12) Anderson-Bjork [.226, .233]
# FalsePostion(3) (King?) .264
# A42() 0.191 but convergence guaranteed
# Order8() 8th order 4 steps: .225 (log10(8)/4
# Order16() 16th order 5 steps .240
# Order5(): 5th order, 4 steps. 0.1747

# A zero is found by specifying:
# the method to use <: AbstractUnivariateZeroMethod
# the function(s) <: CallableFunction
# the initial state through a value for x either x, [a,b], or (a,b) <: AbstractUnivariateZeroState
# the options (e.g., tolerances) <: UnivariateZeroOptions

# The minimal amount needed to add a method, is to define a Method and an update_state method.

### Methods
abstract type AbstractUnivariateZeroMethod end
abstract type AbstractBracketing <: AbstractUnivariateZeroMethod end
abstract type AbstractNonBracketing <: AbstractUnivariateZeroMethod end
abstract type AbstractSecant <: AbstractNonBracketing end

### States
abstract type  AbstractUnivariateZeroState end
mutable struct UnivariateZeroState{T,S} <: AbstractUnivariateZeroState where {T,S}
    xn1::T
    xn0::T
    xstar::T
    m::Vector{T}
    fxn1::S
    fxn0::S
    fxstar::S
    fm::Vector{S}
    steps::Int
    fnevals::Int
    stopped::Bool             # stopped, butmay not have converged
    x_converged::Bool         # converged via |x_n - x_{n-1}| < ϵ
    f_converged::Bool         # converged via |f(x_n)| < ϵ
    convergence_failed::Bool
    message::String
end

incfn(o::AbstractUnivariateZeroState, k=1)    = o.fnevals += k
incsteps(o::AbstractUnivariateZeroState, k=1) = o.steps += k

# initialize state for most methods
function init_state(method::Any, fs, x)

    x1 = float(x)
    fx1 = fs(x1); fnevals = 1
    T, S = eltype(x1), eltype(fx1)
    zT, zS =  oneunit(x1) * (0*x1)/(0*x1), oneunit(fx1) * (0*fx1)/(0*fx1)
    state = UnivariateZeroState(x1, zT, zT/Zt*oneunit(x1), T[],
                                fx1, zS, zS, S[],
                                0, fnevals,
                                false, false, false, false,
                                "")
    state
end

## This function is used to reset the state to an initial value
## As initializing a state is somewhat costly, this can be useful when many
## function calls would be used.
## An example usage, might be:
# M = Order1()
# state = Roots.init_state(M, f1, 0.9)
# options = Roots.init_options(M, state)
# out = zeros(Float64, n)
# for (i,x0) in enumerate(linspace(0.9, 1.0, n))
#    Roots.init_state!(state, M, f1, x0)
#    out[i] = find_zero(M, f1, options, state)
# end
# init_state has a method call variant
function init_state!(state::UnivariateZeroState{T,S}, x1::T, x0::T, m::Vector{T}, y1::S, y0::S, fm::Vector{S}) where {T,S}
    state.xn1 = x1
    state.xn0 = x0
    empty!(state.m); append!(state.m,  m)
    state.m = m
    state.fxn1 = y1
    state.fxn0 = y0
    empty!(state.fm); append!(state.fm,fm)
    state.steps = 0
    state.fnevals = 2
    state.stopped = false
    state.x_converged = false
    state.f_converged = false
    state.convergence_failed = false
    state.message = ""

    nothing
end

function init_state!(state::UnivariateZeroState{T,S}, ::AbstractUnivariateZeroMethod, fs, x) where {T, S}
    x1 = float(x)
    fx1 = fs(x1)

    init_state!(state, x1, oneunit(x1) * (0*x1)/(0*x1), T[],
                fx1, oneunit(fx1) * (0*fx1)/(0*fx1), S[])
end

function Base.copy(state::UnivariateZeroState{T,S}) where {T, S}
    UnivariateZeroState(state.xn1, state.xn0, state.xstar, copy(state.m),
                        state.fxn1, state.fxn0, state.fxstar, copy(state.fm),
                        state.steps, state.fnevals,
                        state.stopped, state.x_converged,
                        state.f_converged, state.convergence_failed,
                        state.message)
end

function Base.copy!(dst::UnivariateZeroState{T,S}, src::UnivariateZeroState{T,S}) where {T,S}
    dst.xn1 = src.xn1
    dst.xn0 = src.xn0
    empty!(dst.m); append!(dst.m, src.m)
    dst.fxn1 = src.fxn1
    dst.fxn0 = src.fxn0
    empty!(dst.fm); append!(dst.fm, src.fm)
    dst.steps = src.steps
    dst.fnevals = src.fnevals
    dst.stopped = src.stopped
    dst.x_converged = src.x_converged
    dst.f_converged = src.f_converged
    dst.convergence_failed = src.convergence_failed
    dst.message = src.message
    nothing
end

### Options
mutable struct UnivariateZeroOptions{Q,R,S,T}
    xabstol::Q
    xreltol::R
    abstol::S
    reltol::T
    maxevals::Int
    maxfnevals::Int
    strict::Bool
end

"""
    default_tolerances(Method, [T], [S])

The default tolerances for most methods are `xatol=eps(T)`,
`xrtol=eps(T)`, `atol=4eps(S)`, and `rtol=4eps(S)`, with the proper
units (absolute tolerances have the units of `x` and `f(x)`; relative
tolerances are unitless). For `Complex{T}` values, `T` is used.

The number of iterations is limited by `maxevals=40`, the number of
function evaluations is not capped, due to `maxfnevals=typemax(Int)`,
and `strict=false`.

"""
default_tolerances(M::AbstractUnivariateZeroMethod) = default_tolerances(M, Float64, Float64)
function default_tolerances(::AbstractUnivariateZeroMethod, ::Type{T}, ::Type{S}) where {T, S}
    xatol = eps(real(T)) * oneunit(real(T))
    xrtol = eps(real(T))  # unitless
    atol = 4.0 * eps(real(float(S))) * oneunit(real(S))
    rtol = 4.0 * eps(real(float(S))) * one(real(S))
    maxevals = 40
    maxfnevals = typemax(Int)
    strict = false
    (xatol, xrtol, atol, rtol, maxevals, maxfnevals, strict)
end

function init_options(M::AbstractUnivariateZeroMethod,
                      state::UnivariateZeroState{T,S};
                      kwargs...
                      ) where {T, S}

    d = kwargs

    defs = default_tolerances(M, T, S)
    options = UnivariateZeroOptions(get(d, :xatol, get(d, :xabstol, defs[1])),
                                    get(d, :xrtol, get(d, :xreltol, defs[2])),
                                    get(d, :atol,  get(d, :abstol,  defs[3])),
                                    get(d, :rtol,  get(d, :reltol,  defs[4])),
                                    get(d, :maxevals,   get(d, :maxsteps, defs[5])),
                                    get(d, :maxfnevals, defs[6]),
                                    get(d, :strict,     defs[7]))
    options
end

# reset options to default values
function init_options!(options::UnivariateZeroOptions{Q,R,S,T}, M::AbstractUnivariateZeroMethod) where {Q, R, S, T}

    defs = default_tolerances(M, Q, S)
    options.xabstol = defs[1]
    options.xreltol = defs[2]
    options.abstol = defs[3]
    options.reltol = defs[4]
    options.maxevals = defs[5]
    options.maxfnevals = defs[6]
    options.strict = defs[7]

    nothing
end

## Tracks (for logging actual steps)
## when no logging this should get optimized out to avoid a branch
abstract type AbstractTracks end
struct NullTracks <: AbstractTracks end
# api
log_step(s::NullTracks, M, x, init=false) = nothing

mutable struct Tracks{T,S} <: AbstractTracks
xs::Vector{T}
fs::Vector{S}
end
Tracks(s::UnivariateZeroState{T,S}) where {T, S} = Tracks(T[],S[])

log_step(s::Tracks, M::Any, o, ::Any) = log_step(s, M, o)

function log_step(s::Tracks, M::Any, o)
    push!(s.xs, o.xn1)
    push!(s.fs, o.fxn1)
    nothing
end
function show_tracks(s::Tracks, M::AbstractUnivariateZeroMethod)
    for (i, (xi, fxi)) in enumerate(zip(s.xs, s.fs))
        println(@sprintf("%s = % 18.16f,\t %s = % 18.16f", "x_$(i-1)", float(xi), "fx_$(i-1)", float(fxi)))
    end
    println("")
end

### Functions
abstract type CallableFunction end

struct FnWrapper <: CallableFunction
   f
end

struct DerivativeFree{F} <: CallableFunction
    f::F
end

struct FirstDerivative{F,FP} <: CallableFunction
    f::F
    fp::FP
end

struct SecondDerivative{F,FP,FPP} <: CallableFunction
    f::F
    fp::FP
    fpp::FPP
end

(F::FnWrapper)(x::Number) = first(F.f(x))
(F::DerivativeFree)(x::Number) = first(F.f(x))
(F::FirstDerivative)(x::Number) = first(F.f(x))
(F::SecondDerivative)(x::Number) = first(F.f(x))

## Return f, f/f'
function fΔx(F::DerivativeFree, x)
    F.f(x)
end

function fΔx(F::Union{FirstDerivative, SecondDerivative},x)
    fx, fpx = F.f(x), F.fp(x)
    fx, fx/fpx
end
function fΔx(F, x)
    F(x)
end

# return f, f/f', f'/f'' (types T, S, S)
fΔxΔΔx(F::DerivativeFree, x) = F.f(x)
fΔxΔΔx(F::FirstDerivative, x) = error("no second derivative defined")
function fΔxΔΔx(F::SecondDerivative, x)
    fx, fp, fpp = F.f(x), F.fp(x), F.fpp(x)
    (fx, fx/fp, fp/fpp)
end
fΔxΔΔx(F, x) = F(x)

# allows override for function, if desired
# the default for this specializes on the function passed
# in. When specialization occurs there is overhead due to compilation
# costs which can be amortized over subsequent calls to `find_zero`.
# However, if a function is only used once, then using `fzero` will be
# faster. (The difference is clear between `@time` and `@btime` to measure
# execution time)
#                 first call    subsequent calls  default for
# specialize       slower         faster          find_zero
# no specialize    faster         slower          fzero
#
callable_function(fs::Any) = _callable_function(fs)
function _callable_function(fs)
    if isa(fs, Tuple)
        length(fs)==1 && return DerivativeFree(fs[1])
        length(fs)==2 && return FirstDerivative(fs[1],fs[2])
        return SecondDerivative(fs[1],fs[2],fs[3])
    end
    DerivativeFree(fs)
end

## Assess convergence
@inline function _is_f_approx_0(fa, a, atol, rtol, relaxed::Any)
    aa, afa = abs(a), abs(fa)
    tol = max(_unitless(atol), _unitless(aa) * rtol)
    tol = abs(_unitless(tol))^(1/3)  # relax test
    afa < tol * oneunit(afa)
end
@inline function _is_f_approx_0(fa, a, atol, rtol)
    aa, afa = abs(a), abs(fa)
    tol = max(_unitless(atol), _unitless(aa) * rtol)
    afa < tol * oneunit(afa)
end

"""
   Roots.assess_convergence(method, state, options)

Assess if algorithm has converged.

If alogrithm hasn't converged returns `false`.

If algorithm has stopped or converged, return `true` and sets one of `state.stopped`, `state.x_converged`,  `state.f_converged`, or `state.convergence_failed`; as well, a message may be set.

* `state.x_converged = true` if `abs(xn1 - xn0) < max(xatol, max(abs(xn1), abs(xn0)) * xrtol)`

* `state.f_converged = true` if  `|f(xn1)| < max(atol, |xn1|*rtol)`

* `state.convergence_failed = true` if xn1 or fxn1 is `NaN` or an infinity

* `state.stopped = true` if the number of steps exceed `maxevals` or the number of function calls exceeds `maxfnevals`.

In `find_zero`, stopped values (and x_converged) are checked for convergence with a relaxed tolerance.

"""
function assess_convergence(method::Any, state::UnivariateZeroState{T,S}, options) where {T,S}

    xn0, xn1 = state.xn0, state.xn1
    fxn1 = state.fxn1

    if (state.x_converged || state.f_converged || state.stopped)
        if isnan(state.xstar)
            state.xstar, state.fxstar =  xn1, fxn1
        end
        return true
    end

    if isnan(xn1) || isnan(fxn1)
        state.convergence_failed = true
        state.message *= "NaN produced by algorithm. "
        return true
    end

    if isinf(xn1) || isinf(fxn1)
        state.convergence_failed = true
        state.message *= "Inf produced by algorithm. "
        return true
    end

    # f(xstar) ≈ xstar * f'(xstar)*eps(), so we pass in lambda
    if _is_f_approx_0(fxn1, xn1, options.abstol, options.reltol)
        state.xstar, state.fxstar = xn1, fxn1
        state.f_converged = true
        return true
    end

    # stop when xn1 ~ xn.
    # in find_zeros there is a check that f could be a zero with a relaxed tolerance
    if abs(xn1 - xn0) < max(options.xabstol, max(abs(xn1), abs(xn0)) * options.xreltol)
        state.xstar, state.fxstar = xn1, fxn1
        state.message *= "x_n ≈ x_{n-1}. "
        state.x_converged = true
        return true
    end

    if state.steps > options.maxevals
        state.stopped = true
        state.message *= "Too many steps taken. "
        return true
    end

    if state.fnevals > options.maxfnevals
        state.stopped = true
        state.message *= "Too many function evaluations taken. "
        return true
    end

    return false
end

function show_trace(method, N, state, tracks)
    converged = state.x_converged || state.f_converged
    println("Results of univariate zero finding:\n")
    if converged
        println("* Converged to: $(state.xn1)")
        if N == nothing || isa(method, AbstractBracketing)
            println("* Algorithm: $(method)")
        else
            println("* Algorithm: $(method), with possible bracketing with $N")
        end
        println("* iterations: $(state.steps)")
        println("* function evaluations: $(state.fnevals)")
        state.x_converged && println("* stopped as x_n ≈ x_{n-1} using atol=xatol, rtol=xrtol")
        state.f_converged && state.message == "" && println("* stopped as |f(x_n)| ≤ max(δ, max(1,|x|)⋅ϵ) using δ = atol, ϵ = rtol")
        state.message != "" && println("* Note: $(state.message)")
    else
        println("* Convergence failed: $(state.message)")
        println("* Algorithm $(method)")
    end
    println("")
    println("Trace:")
    show_tracks(tracks, method)
end

"""

    find_zero(fs, x0, method; kwargs...)

Interface to one of several methods for find zeros of a univariate function.

# Initial starting value

For most methods, `x0` is a scalar value indicating the initial value
in the iterative procedure. (Secant methods can have a tuple specify
their initial values.) Values must be a subtype of `Number` and have
methods for `float`, `real`, and `oneunit` defined.

For bracketing intervals, `x0` is specified as a tuple or a vector. A bracketing interval, (a,b), is one where f(a) and f(b) have different signs.

# Specifying a method

A method is specified to indicate which algorithm to employ:

* There are methods for bisection where a bracket is specified: `Bisection`, `Roots.A42`, `FalsePosition`

* There are several derivative-free methods: cf. `Order0`, `Order1` (secant method), `Order2` (Steffensen), `Order5`, `Order8`, and `Order16`, where the number indicates the order of the convergence. Methods `Roots.Order1B` and `Roots.Order2B` implement methods useful when the desired zero has a multiplicity.

* There are some classical methods where derivatives are required: `Roots.Newton`, `Roots.Halley`, `Roots.Schroder`. (The are not exported.)

For more detail, see the help page for each method (e.g., `?Order1`).

If no method is specified, the default method depends on `x0`:

* If `x0` is a scalar, the default is the slower, but more robust `Order0` method.

* If `x0` is a tuple or vector indicating a *bracketing* interval, the `Bisection` method is used. (The exact algorithm depends on the number type, the tolerances, and `verbose`.)

# Specifying the function

The function(s) are passed as the first argument.

For the few methods that use a derivative (`Newton`, `Halley`, `Shroder`, and
optionally `Order5`) a tuple of functions is used.

# Optional arguments (tolerances, limit evaluations, tracing)

* `xatol` - absolute tolerance for `x` values. Passed to `isapprox(x_n, x_{n-1})`
* `xrtol` - relative tolerance for `x` values. Passed to `isapprox(x_n, x_{n-1})`
* `atol`  - absolute tolerance for `f(x)` values.
* `rtol`  - relative tolerance for `f(x)` values.
* `maxevals`   - limit on maximum number of iterations
* `maxfnevals` - limit on maximum number of function evaluations
* `strict` - if `false` (the default), when the algorithm stops, possible zeros are checked with a relaxed tolerance
* `verbose` - if `true` a trace of the algorithm will be shown on successful completion.

See the help string for `Roots.assess_convergence` for details on
convergence. See the help page for `Roots.default_tolerances(method)`
for details on the default tolerances.

In general, with floating point numbers, convergence must be
understood as not an absolute statement. Even if mathematically x is
an answer and xstar the floating point realization, it may be that
f(xstar) - f(x) = f(xstar) ≈ f'(x) ⋅ eps(x), so tolerances must be
appreciated, and at times specified.

For the `Bisection` methods, convergence is guaranteed, so the tolerances are set to be 0 by default.

If a bracketing method is passed in after the method specification, if a bracket is found, the bracketing method will used to identify the zero. This is what `Order0` does by default, with a secant step initially.

Note: The order of the method is hinted at in the naming scheme. A
scheme is order `q` if, with `e_n = x_n - alpha`, `e_{n+1} = C
e_n^q`. If the error `e_n` is small enough, then essentially the error
will gain `q` times as many leading zeros each step. However, if the
error is not small, this will not be the case. Without good initial
guesses, a high order method may still converge slowly, if at all. The
`OrderN` methods have some heuristics employed to ensure a wider range
for convergence at the cost of not faithfully implementing the method,
though those are available through unexported methods.

# Examples:

```
# default methods
find_zero(sin, 3)  # use Order0()
find_zero(sin, (3,4)) # use Bisection()

# specifying a method
find_zero(sin, (3,4), Order1())            # can specify two starting points for secant method
find_zero(sin, 3.0, Order2())              # Use Steffensen method
find_zero(sin, big(3.0), Order16())        # rapid convergence
find_zero(sin, (3, 4), Roots.A42()())      # fewer function calls than Bisection(), in this case
find_zero(sin, (3, 4), FalsePosition(8))   # 1 of 12 possible algorithms for false position
find_zero((sin,cos), 3.0, Roots.Newton())  # use Newton's method
find_zero((sin, cos, x->-sin(x)), 3.0, Roots.Halley())  # use Halley's method

# changing tolerances
fn, x0, xstar = (x -> (2x*cos(x) + x^2 - 3)^10/(x^2 + 1), 3.0,  2.9806452794385368)
find_zero(fn, x0, Order2()) - xstar        # 0.014079847201995843
find_zero(fn, x0, Order2(), atol=0.0, rtol=0.0) # error: x_n ≉ x_{n-1}; just f(x_n) ≈ 0
fn, x0, xstar = (x -> (sin(x)*cos(x) - x^3 + 1)^9,        1.0,  1.117078770687451)
find_zero(fn, x0, Order2())                # 1.1122461983100858
find_zero(fn, x0, Order2(), maxevals=3)    # Roots.ConvergenceFailed: 26 iterations needed

# tracing output
find_zero(x->sin(x), 3.0, Order2(), verbose=true)   # 3 iterations
find_zero(x->sin(x)^5, 3.0, Order2(), verbose=true) # 22 iterations
find_zero(x->sin(x)^5, 3.0, Roots.Order2B(), verbose=true) # 2 iterations
```
"""
function find_zero(fs, x0, method::AbstractUnivariateZeroMethod,
                   N::Union{Nothing, AbstractBracketing}=nothing;
                   tracks::AbstractTracks=NullTracks(),
                   verbose=false,
                   kwargs...)

    F = callable_function(fs)
    state = init_state(method, F, x0)
    options = init_options(method, state; kwargs...)

    l = (verbose && isa(tracks, NullTracks)) ? Tracks(eltype(state.xn1)[], eltype(state.fxn1)[]) : tracks

    if N == nothing  || isa(method, AbstractBracketing)
        xstar = find_zero(method, F, options, state, l)
    else
        xstar = find_zero(method, N, F, options, state, l)
    end

    if verbose
        show_trace(method, N, state, l)
    end

    if isnan(xstar)
        throw(ConvergenceFailed("Stopped at: xn = $(state.xn1)"))
    else
        return xstar
    end

end

find_zero(f, x0::T; kwargs...)  where {T <: Number}= find_zero(f, x0, Order0(); kwargs...)
find_zero(f, x0::Vector; kwargs...) = find_zero(f, x0, Bisection(); kwargs...)
find_zero(f, x0::Tuple; kwargs...) = find_zero(f, x0, Bisection(); kwargs...)

# Main method
# return a zero or NaN.
## Updates state, could be `find_zero!(state, M, F, options, l)...
function find_zero(M::AbstractUnivariateZeroMethod,
                   F,
                   options::UnivariateZeroOptions,
                   state::AbstractUnivariateZeroState,
                   l::AbstractTracks=NullTracks()
                   )  #where {T<:Number, S<:Number}

    log_step(l, M, state, :init)

    while true
        val = assess_convergence(M, state, options)
        val && break
        update_state(M, F, state, options)
        log_step(l, M, state)
        incsteps(state)
    end

    return decide_convergence(M, F, state, options)
end

function find_zero(M::AbstractUnivariateZeroMethod,
                   F,
                   state::AbstractUnivariateZeroState,
                   l::AbstractTracks=NullTracks()
                   )  #where {T<:Number, S<:Number}

    options = init_options(M, state)
    find_zero(M, F, options, state, l)
end

# state has stopped, this identifies if it has converged
function decide_convergence(M::AbstractUnivariateZeroMethod,  F, state::UnivariateZeroState{T,S}, options) where {T,S}
    xn1 = state.xstar
    fxn1 = state.fxstar

    if (state.stopped || state.x_converged) && !(state.f_converged)
        ## stopped is a heuristic, x_converged can mask issues
        ## if strict == false, this will also check f(xn) ~ - with a relaxed
        ## tolerance

        ## are we at a crossing values?
        ## seems worth a check for 2 fn evals.
        if T <: Real && S <: Real
            for u in (prevfloat(xn1), nextfloat(xn1))
                fu = first(F(u))
                if iszero(fu) || _unitless(fu * fxn1) < 0
                    state.message *= "Change of sign at xn identified. "
                    state.f_converged = true
                end
            end
        end

        if options.strict
            if state.x_converged
                state.f_converged = true
            else
                state.convergence_failed = true
            end
        else
            xstar, fxstar = state.xn1, state.fxn1
            if _is_f_approx_0(fxstar, xstar, options.abstol, options.reltol, :relaxed)
                state.xstar, state.fxstar = xstar, fxstar
                msg = "Algorithm stopped early, but |f(xn)| < ϵ^(1/3), where ϵ depends on xn, rtol, and atol. "
                state.message = state.message == "" ? msg : state.message * "\n\t" * msg
                state.f_converged = true
            else
                state.convergence_failed = true
            end
        end
    end

    if state.f_converged
        return state.xstar
    end

    nan = NaN * xn1
    if state.convergence_failed
        return nan
    end
    return nan

end

# Switch to bracketing method
#M = Bisection()  # exact for floating point
#M = AlefeldPotraShi() # *usually* exact
#M = Brent()          # a bit faster, but not always convergent, as implemented (cf. RootTesting)
run_bisection(f, xs, state, options) = run_bisection(AlefeldPotraShi(), f, xs, state, options)
function run_bisection(N::AbstractBracketing, f, xs, state, options)
    steps, fnevals = state.steps, state.fnevals
    init_state!(state, N, f, xs)
    state.steps += steps
    state.fnevals += fnevals # could avoid 2 fn calls, with fxs
    init_options!(options, N)
    find_zero(N, f, options, state)
    a, b = xs
    u,v = a > b ? (b, a) : (a, b)
    state.message *= "Bracketing used over ($u, $v), those steps not shown. "
    return nothing
end

# Robust version using some tricks: idea from algorithm described in
# [The SOLVE button from the
# HP-34]C(http://www.hpl.hp.com/hpjournal/pdfs/IssuePDFs/1979-12.pdf).
# * use bracketing method if one identifed
# * limit steps so as not too far or too near the previous one
# * if not decreasing, use a quad step upto 4 times to bounce out of trap, if possible
# First uses M, then N if bracket is identified
function find_zero(M::AbstractUnivariateZeroMethod,
                   N::AbstractBracketing,
                   F,
                   options::UnivariateZeroOptions,
                   state::AbstractUnivariateZeroState,
                   l::AbstractTracks=NullTracks()
                   )

    log_step(l, M, state, :init)
    state0 = copy(state)
    quad_ctr = 0

    while true

        if assess_convergence(M, state, options)
            break
        end

        copy!(state0, state)
        update_state(M, F, state0, options) # state0 is proposed step

        adj = false
        ## did we find a zero or a bracketing interval?
        if iszero(state0.fxn1)
            copy!(state, state0)
            state.xstar, state.fxstar = state.xn1, state.fxn1
            state.f_converged = true
            break
        elseif sign(state0.fxn0) * sign(state0.fxn1) < 0
            copy!(state, state0)
            a, b = state0.xn0, state0.xn1 # could save some fn calls here
            run_bisection(N, F, (a, b), state, options)
            break
        end

        ## did we move too far?
        r, a, b = state0.xn1, state.xn0, state.xn1
        Δr = abs(r - b)
        Δx = abs(b - a)
        ts,TB = 1e-3, 1e2 # too small, too big
        if  Δr >= TB * Δx
            adj = true
            r = b + sign(r-b) * TB * Δx  ## too big
            state0.xn1 = r
            state0.fxn1 = F(r)
            incfn(state)
        elseif Δr  <= ts *  Δx
            adj = true
            r = b + sign(r - b) * ts * Δx
            fr = F(r)
            incfn(state)
            state0.xn1 = r
            state0.fxn1 = fr
        end

        # a sign change after shortening?
        if sign(state.fxn1) * sign(state0.fxn1) < 0
            state.xn0, state.fxn0 = state.xn1, state.fxn1
            state.xn1, state.fxn1 = state0.xn1, state0.fxn1
            a, b = state.xn0, state.xn1
            run_bisection(N, F, (a, b), state, options)
            break
        end

        ## did we improve?
        if adj || abs(state0.fxn1) < abs(state.fxn1)
            if isnan(state0.xn1) || isnan(state0.fxn1) || isinf(state0.xn1) || isinf(state0.fxn1)
                break
            end
            copy!(state, state0)
            log_step(l, M, state)
            incsteps(state)
            quad_ctr = 0
            continue
        end

        ## try quad_vertex, unless that has gotten old
        if quad_ctr > 4
            copy!(state, state0)
            state.stopped = true
            break
        else
            quad_ctr += 1
            r = quad_vertex(state0.xn1, state0.fxn1, state.xn1, state.fxn1, state.xn0, state.fxn0)

            if isnan(r) || isinf(r)
                copy!(state, state0)
            else

                fr = F(r)
                incfn(state)

                state0.xn1 = r
                state0.fxn1 = fr
                incfn(state)
                copy!(state, state0)
            end
        end

        log_step(l, M, state)
        incsteps(state)
    end

    decide_convergence(M, F, state, options)
end


text:
"""
    TimeSeriesHandler{T} <: AbstractDH{T}

Type for handling time series data.  As with DataHandler it is intended taht most of
the reformatting of the dataframe is done before passing it to an instance of this type.

The parameter T specifies the datatype of the input, output arrays.
"""
mutable struct TimeSeriesHandler{T} <: AbstractDH{T}

    # dataframe where all data is kept
    df::DataTable

    colsInput::Array{Symbol, 1}
    colsOutput::Array{Symbol, 1}

    colsNormalize::Array{Symbol, 1}

    mu::Array{T, 1}
    norm::Array{T, 1}
    userange::Bool

    # training data, is a subset of df
    dfTrain::DataTable
    # testing data, is a subset of df
    dfTest::DataTable

    # these are the arrays of training, testing data
    X_train::Array{T}
    y_train::Array{T}
    X_test::Array{T}
    y_test::Array{T}

    # this is where predictions are kept
    yhat::Array{T}
    yhat_train::Array{T}

    # the symbol of the timeindex
    timeindex::Symbol
    # the length of sequences (only fixed for now)
    seq_length::Int64

    """
        TimeSeriesHandler(df::DataTable, timeindex::Symbol, seq_length::Integer;
                          shuffle::Bool=false, n_test_sequences::Integer=0,
                          input_cols::Array{Symbol}=Symbol[],
                          output_cols::Array{Symbol}=Symbol[],
                          normalize_cols::Array{Symbol}=Symbol[],
                          assign::Bool=false, userange::Bool=false)

    Creates an object converting the data contained in `df` to properly formatted
    time series data with fixed-length sequences.  It is required that the input dataframe
    have a time index, the name of which should be passed as the second argument.
    The keyword arguments are otherwise the same as for `DataHandler` but note that in
    this case it makes sense for `input_cols == output_cols` (in fact that is probably the
    most common use case).
    """
    function TimeSeriesHandler{T}(df::DataTable, timeindex::Symbol, seq_length::Integer;
                                  shuffle::Bool=false,
                                  n_test_sequences::Integer=0,
                                  input_cols::Array{Symbol}=Symbol[],
                                  output_cols::Array{Symbol}=Symbol[],
                                  normalize_cols::Array{Symbol}=Symbol[],
                                  assign::Bool=false,
                                  userange::Bool=false) where T
        ndf = copy(df)
        o = new(ndf, input_cols, output_cols, normalize_cols)
        o.userange = userange
        o.timeindex = timeindex
        o.seq_length = seq_length
        splitByNSequences!(o, n_test_sequences, assign=assign)
        computeNormalizeParameters!(o, dataset=:dfTrain)
        if canNormalize(o)
            normalizeTrain!(o)
            if size(o.dfTest)[1] > 0 normalizeTest!(o) end
        end
        return o
    end
end
export TimeSeriesHandler

"""
Used by assignTrain!, and assignTest!
"""
function _get_assign_data{T}(dataframe::Symbol, dh::TimeSeriesHandler{T}; sort::Bool=true)
    df = getfield(dh, dataframe)
    if isempty(df) return end
    if sort sort!(df, cols=[dh.timeindex]) end
    npoints = size(df)[1] - dh.seq_length
    X = Array{T}(npoints, dh.seq_length, length(dh.colsInput))
    y = Array{T}(npoints, length(dh.colsOutput))
    # the loop can only go this far
    for i in 1:npoints
        nextx = reshape(convert(Array, df[i:(i+dh.seq_length-1), dh.colsInput]),
                        (1, dh.seq_length, length(dh.colsInput)))
        X[i, :, :] = nextx
        y[i, :] = convert(Array, df[i+dh.seq_length, dh.colsOutput])
    end
    return X, y
end

"""
    assignTrain!(dh[, df; sort=true, parallel=false])

Assigns the training data.  X output will be of shape (samples, seq_length, seq_width).
If `sort` is true, will sort the dataframe first.  One should be extremely careful if
`sort` is false.  If `parallel` is true the data will be generated in parallel (using
workers, not threads).  This is useful because this data manipulation is complicated and
potentially slow.

If a dataframe is provided, data will be assigned from it.  Alternatively, one can provide
a vector of dataframes.  This is useful because sequences which cross the boundaries
of the dataframes will *not* be created.

**TODO** I'm pretty sure the parallel version isn't working right because it doesn't
use shared arrays.  Revisit in v0.5 with threads.
"""
function assignTrain!(dh::TimeSeriesHandler; sort::Bool=true, parallel::Bool=false)
    if isempty(dh.dfTrain)
        throw(ErrorException("Trying to build variables from empty training dataframe."))
    end
    # TODO, again parallel not currently working
    if parallel
        gad = _get_assign_data_parallel
    else
        gad = _get_assign_data
    end
    dh.X_train, dh.y_train = gad(:dfTrain, dh, sort=sort)
end

function assignTrain!(dh::TimeSeriesHandler, df::AbstractDataTable;
                      sort::Bool=true, parallel::Bool=false)
    dh.dfTrain = df
    assignTrain!(dh, sort=sort, parallel=parallel)
end

function assignTrain!{DFT<:AbstractDataTable}(dh::TimeSeriesHandler, dfs::Vector{DFT};
                                              sort::Bool=true)
    vcat_tuples(x,y) = vcat(x[1], y[1]), vcat(x[2], y[2])
    dh.X_train, dh.y_train = mapreduce(vcat_tuples, dfs) do df
        dh.dfTrain = df
        _get_assign_data(:dfTrain, dh, sort=sort)
    end
end
export assignTrain!

"""
    assignTest!(dh[, df; sort=true])

Assigns the test data.  X output will be of shape (samples, seq_length, seq_width).
One should be extremely careful if not sorting.

If a dataframe is provided, data will be assigned from it.

Note that in the time series case this isn't very useful.  One should instead use
one of the assigned prediction functions.
"""
function assignTest!(dh::TimeSeriesHandler; sort::Bool=true)
    if isempty(dh.dfTest)
        throw(ErrorException("Trying to build variables from empty test dataframe."))
    end
    dh.X_test, dh.y_test = _get_assign_data(:dfTest, dh, sort=sort)
end

function assignTest!(dh::TimeSeriesHandler, df::AbstractDataTable; sort::Bool=true)
    dh.dfTest = df
    assignTest!(dh, sort=sort)
end
export assignTest!

"""
    assign!(dh::TimeSeriesHandler; sort::Bool=true)

Assigns both training and testing data for the `TimeSeriesHandler`.
"""
function assign!(dh::TimeSeriesHandler; sort::Bool=true)
    assignTrain!(dh, sort=sort)
    assignTest!(dh, sort=sort)
    return
end
export assign!

"""
## TimeSeriesHandler
If the `DataHandler` argument is a `TimeSeriesHandler`, one can optionall pass the flag
`sort` which specifies that the data should be sorted first.
"""
getData(dh::TimeSeriesHandler; sort::Bool=true) = _get_assign_data(:df, dh, sort=sort)
export getData

"""
    split!(dh::TimeSeriesHandler, τ₀::Integer; assign::Bool=true, sort::Bool=true)

Splits the data by time-index.  All datapoints with τ up to and including the timeindex
given (τ₀) will be in the training set, while all those with τ > τ₀ will be in
the test set.
"""
function split!(dh::TimeSeriesHandler, τ₀::Integer; assign::Bool=true, sort::Bool=true)
    constr = dropnull(dh.df[dh.timeindex]) .≤ τ₀
    dfTrain = dh.df[constr, :]
    dfTest = dh.df[!constr, :]
    if !isempty(dfTrain)
        dh.dfTrain = dfTrain
    end
    if !isempty(dfTest)
        dh.dfTest = dfTest
    end
    if assign assign!(dh, sort=sort) end
    return
end
export split!

"""
    splitByNSequences!(dh::TimeSeriesHandler, n_sequences::Integer;
                       assign::Bool=true, sort::Bool=true)

Splits the dataframe by the number of sequences in the test set.  This does nothing to
account for the possibility of missing data.

Note that the actual number of usable test sequences in the resulting test set is of
course greater than n_sequences.
"""
function splitByNSequences!(dh::TimeSeriesHandler, n_sequences::Integer;
                            assign::Bool=true, sort::Bool=true)
    τ_max = maximum(dropnull(dh.df[:, dh.timeindex]))
    nτ = n_sequences*dh.seq_length
    τ₀ = τ_max - nτ
    split!(dh, τ₀, assign=assign, sort=sort)
    return
end
export splitByNSequences!

"""
Squashes the last two indices of a rank-3 tensor into a matrix.  For internal use.
"""
function _squashXTensor(X::Array, dh::TimeSeriesHandler)
    return reshape(X, (size(X)[1], dh.seq_length*length(dh.colsInput)))
end

"""
    getSquashedTrainMatrix(dh::TimeSeriesHandler)

Gets a training input tensor in which all the inputs are arranged along a single axis
(i.e. in a matrix).

Assumes the handler's X_train is defined.
"""
function getSquashedTrainMatrix(dh::TimeSeriesHandler)
    if !isdefined(dh, :X_train)
        error("Training input not defined.  Should call assign!.")
    end
    return _squashXTensor(dh.X_train, dh)
end
export getSquashedTrainMatrix

"""
    getSquashedTestMatrix(dh::TimeSeriesHandler)

Gets a test input tensor in which all the inputs are arranged along a single axis
(i.e. in a matrix).

Assumes the handler's X_test is defined.
"""
function getSquashedTestMatrix(dh::TimeSeriesHandler)
    if !isdefined(dh, :X_test)
        error("Test input not defined.  Should call assign!.")
    end
    return _squashXTensor(dh.X_test, dh)
end
export getSquashedTestMatrix

"""
    getSquashedTrainData(dh::TimeSeriesHandler; flatten::Bool=false)

Gets the training X, y pair where X is squashed using `getSquahdedTrainMatrix`.
If `flatten`, also flatten `y`.
"""
function getSquashedTrainData(dh::TimeSeriesHandler; flatten::Bool=false)
    if flatten
        @assert size(dh.y_train)[2] == 1 "Attempting to flatten mis-shaped matrix."
        y_train = squeeze(dh.y_train, 2)
    else
        y_train = dh.y_train
    end
    return getSquashedTrainMatrix(dh), y_train
end
export getSquashedTrainData

"""
    trainOnMatrix(train::Function, dh::TimeSeriesHandler; flatten::Bool=true)

Trains an object designed to take a matrix (as opposed to a rank-3 tensor) as input.
The first argument is the function used to train.

If flatten, the training output is converted to a vector.  Most methods that take matrix
input take vector target data.

Assumes the function takes input of the form train(X, y).
"""
function trainOnMatrix(train::Function, dh::TimeSeriesHandler; flatten::Bool=true)
    # turn off flatten if y_train is a matrix
    if length(dh.colsOutput) > 1
        flatten = false
    end
    X = getSquashedTrainMatrix(dh)
    if flatten
        y_train = squeeze(dh.y_train, 2)
    else
        y_train = dh.y_train
    end
    train(X, y_train)
end
export trainOnMatrix

"""
    _check_and_reshape(lastX, on_matrix)

Private method used by generateSequence.
"""
function _check_and_reshape(lastX::Array, on_matrix::Bool)
    if on_matrix
        lastX_reshaped = reshape(lastX, (1, prod(size(lastX))))
    else
        lastX_reshaped = lastX
    end
    lastX_reshaped
end

"""
    _get_next_X(indim, outin_idx, yhat, newcols_funcs)

Private method used by `generateSequence`.
"""
function _get_next_X!{T}(lastX::Array, yhat::Array{T}, outin_idx::Vector)
    for (i, idx) ∈ enumerate(outin_idx)
        lastX[1, end, idx] = yhat[i]
    end
end

function _get_next_X!{T}(lastX::Array, yhat::Array{T}, outin_idx::Vector, newcol_funcs::Dict)
    _get_next_X!(lastX, yhat, outin_idx)
    for (n, f) ∈ newcol_funcs
        lastX[1, end, n] = f(lastX[1, end-1, :])
    end
end

function _get_next_X!{T}(lastX::Array, yhat::Array{T}, outin_idx::Vector,
                         newcol_func::Function)
    _get_next_X!(lastX, yhat, outin_idx)
    dict = newcol_func(lastX[1, end-1, :])
    for (n, x) ∈ dict
        lastX[1, end, n] = x
    end
end

"""
    generateSequence(predict, dh, seq_length[, newcol_func; on_matrix=false])

Uses the supplied prediction function `predict` to generate a sequence of length `seq_length`.
The sequence uses the end of the training dataset as initial input.

Note that this only makes sense when the output columns are a subset of the input columns.

If a function returning a dictionary or a dictionary of functions `newcol_func` is supplied,
every time a new row of the input is generated, it will have columns specified by
`newcol_func`.  The dictionary should have keys equal to the column numbers of columns
in the input matrix and values equal to functions that take a `Vector` (the previous
input row) and output a new value for the column number given by the key.  The column
numbers correspond to the index of the column in the specified input columns.

If `on_matrix` is true, the prediction function will take a matrix as input rather than
a rank-3 tensor.
"""
# we keep this method separate for a very tiny amount of efficiency
function generateSequence{T}(predict::Function, dh::TimeSeriesHandler{T},
                             seq_length::Integer;
                             on_matrix::Bool=false)
    outin_idx = indexin(dh.colsOutput, dh.colsInput)
    testsize = (seq_length, length(dh.colsInput))
    lastX = reshape(convert(Array{T}, dh.dfTrain[(end-dh.seq_length+1):end, dh.colsInput]),
                    (1, dh.seq_length, length(dh.colsInput)))
    lastX_reshaped = _check_and_reshape(lastX, on_matrix)
    yhat = Array{T}(testsize)
    for i in 1:testsize[1]
        pred = predict(lastX_reshaped)
        lastX = circshift(lastX, (0, -1, 0))
        _get_next_X!(lastX, pred, outin_idx)
        yhat[i, :] = lastX[1, end, :]
        lastX_reshaped = _check_and_reshape(lastX, on_matrix)
    end
    yhat
end

function generateSequence{T}(predict::Function, dh::TimeSeriesHandler{T},
                             seq_length::Integer, newcol_funcs::Union{Dict,Function};
                             on_matrix::Bool=false)
    outin_idx = indexin(dh.colsOutput, dh.colsInput)
    testsize = (seq_length, length(dh.colsInput))
    lastX = reshape(convert(Array{T}, dh.dfTrain[(end-dh.seq_length+1):end, dh.colsInput]),
                    (1, dh.seq_length, length(dh.colsInput)))
    lastX_reshaped = _check_and_reshape(lastX, on_matrix)
    yhat = Array{T}(testsize)
    for i ∈ 1:testsize[1]
        pred = predict(lastX_reshaped)
        lastX = circshift(lastX, (0, -1, 0))
        _get_next_X!(lastX, pred, outin_idx, newcol_funcs)
        yhat[i, :] = lastX[1, end, :]
        lastX_reshaped = _check_and_reshape(lastX, on_matrix)
    end
    yhat
end
export generateSequence

"""
    generateTest(predict::Function, dh::TimeSeriesHandler; on_matrix::Bool=true)

Uses the supplied prediction function to attempt to predict the entire test set.
Note that this assumes that the test set is ordered, sequential and immediately follows
the training set.

See the documentation for `generateSequence`.
"""
function generateTest(predict::Function, dh::TimeSeriesHandler;
                      on_matrix::Bool=true)
    return generateSequence(predict, dh, size(dh.dfTest)[1], on_matrix=on_matrix)
end
export generateTest

"""
    getRawTestTarget(dh::TimeSeriesHandler)

Returns `y_test` directly from the dataframe for comparison with the output of generateTest.
"""
function getRawTestTarget{T}(dh::TimeSeriesHandler{T})
    return convert(Array{T}, dh.dfTest[:, dh.colsOutput])
end

function getRawTestTarget{T}(dh::TimeSeriesHandler{T}, limit::Integer)
    return convert(Array{T}, dh.dfTest[1:limit, dh.colsOutput])
end

function getRawTestTarget{T}(dh::TimeSeriesHandler{T}, slice::UnitRange)
    return convert(Array{T}, dh.dfTest[slice, dh.colsOutput])
end
export getRawTestTarget