text:
__precompile__()
module KeywordDispatch

export KeywordDispatchFunction
export KeywordDispathFunctionNotFound
export KeywordDispatchFunctionNotImplemented
export add_specializer!, remove_specializer!
export @defkwspec

"""
Data structure for a function specializing on keywrod arguments.

Field `specializers` is a dictionary with keys being (highest priority)
keyword arguments and every entry of the form

    (function, dict)

where `function` is a function to be called for current specialization and
`dict` provides further specializations on lower priority keyword arguments
"""
struct KeywordDispatchFunction
	specializers :: Dict
	KeywordDispatchFunction() = new(Dict())
end

"""
`KeywordDispatchFunction` is funcallable (functor) object that accepts
keyword arguments only
"""
function (d :: KeywordDispatchFunction)(; kwargs...)
    dks = Dict(kwargs)
	find_function(dks, d.specializers, notfound(d, collect(keys(dks))))
end

"""
General error
"""
abstract type KeywordDispatchError <: Exception end

"""
Error to be thrown if specializers are not found
"""
struct KeywordDispathFunctionNotFound <: KeywordDispatchError
    kwdfunction :: KeywordDispatchFunction
    keywords :: Vector{Symbol}
end

function Base.showerror(io :: IO, e :: KeywordDispathFunctionNotFound)
    println(io, "Keywords ", e.keywords,
        " a not found in specializers of ", e.kwdfunction)
end

"""
Error to be thrown if the function is requested with
incomplete list of specializers
"""
struct KeywordDispatchFunctionNotImplemented <: KeywordDispatchError
    kwdfunction :: KeywordDispatchFunction
    keywords :: Vector{Symbol}
end

function Base.showerror(io :: IO, e :: KeywordDispatchFunctionNotImplemented)
    println(io,
        "Incomplete keywords specializers ", e.keywords,
        " for ", e.kwdfunction)
end

"""
Function of one argument throwing an error if specializer is not found.
It is used a placeholder when searching for an effective function
"""
function notfound(f, keys)
    x -> throw(KeywordDispathFunctionNotFound(f, keys))
end

function notimplemented(f, keys)
    x -> throw(KeywordDispatchFunctionNotImplemented(f, keys))
end

"""
Finds function applicable to `args` in `dict` which has a form of
`specializers` field of `KeywordDispatchFunction` **and** applies it
to `args`
"""
function find_function(args, dict, default)
	for key in keys(dict)
		if haskey(args, key)
			subdefault, subdict = dict[key]
			return find_function(args, subdict, subdefault)
		end
	end
	default(args)
end

"""
Adds specialized function to a keyword generic function.

- `f` is a new specialized function
- `dispatch` is a `KeywordDispatchFunction` generic function
- `ks` is a sequence of symbols to specialize `dispatch` on. Note
  that the order of symbols matter: they must be ordered from higher priority
  to lower priority
"""
function add_specializer!(f, dispatch :: KeywordDispatchFunction, ks)
    n = length(ks)
    dict = dispatch.specializers
    # Loop over all but last keys:
    # just walking down the dictionary and putting new things on the way
    for i in 1:n-1
        if !haskey(dict, ks[i])
            # No entry: create a placeholder
            ff = notimplemented(dispatch, ks[1:i])
            dict[ks[i]] = (ff, Dict())
        end
        # Simply go down
        default, dict = dict[ks[i]]
    end
    # By now, `dict` is the last dictionary in the structure
    if !haskey(dict, ks[n])
        # no further entries: just create a new one
        dict[ks[n]] = (f, Dict())
    else
        # there is an entry: put new function, keep sub-dictionary
        default, subdict = dict[ks[n]]
        dict[ks[n]] = (f, subdict)
    end
	dispatch
end

"""
Removes a specializer defined by the sequence of keys `ks` from
the keyword-dispatch function `dispatch`.
If the specializer with required keys is not found, the warning is issued
"""
function remove_specializer!(dispatch :: KeywordDispatchFunction, ks)
    n = length(ks)
    dict = dispatch.specializers
    for i in 1:(n-1)
        if !haskey(dict, ks[i])
            warn("remove_specializer!: Specializer for keys ",
                ks[1:i], " not found")
            return dispatch
        end
        f, dict = dict[ks[i]]
    end
    if !haskey(dict, ks[n])
        warn("remove_specializer!: Specializer for keys ",
            ks[1:i], " not found")
        return dispatch
    end
    f, subdict = dict[ks[n]]
    dict[ks[n]] = (notimplemented(dispatch, ks), subdict)
    dispatch
end

"""
Defines new keyword specialization of the keyword-generic function `name`
"""
macro defkwspec(name, kw, args, body)
	return quote
		let f = $(esc(args)) -> ($(esc(body)))
		    add_specializer!(f, $(esc(name)), $(esc(kw)))
        end
    end
end

# find_function(Dict(:x=>2, :y=> 3), add.specializers)

# Example
# add = KeywordDispatchFunction()
# macroexpand(:(@defkwspec add [:x :y] (xx) xx[:x] + xx[:y]))
#
# ff = (:x, :y)
#
# macro foo(x)
#     return quote
#         f() = $(esc(x))
#     end
# end
#
# macroexpand(:(@foo (z,w)))

# add_specializer(args -> args[:x], add, [:x])
# add(x = 2, y = 3)
# @defkwspec add [:x :y] args begin
# 	println("Yay!")
# 	args[:x] + args[:y]
# end

end # module


text:
using Logging
using PortfolioOpt
using JuMP

"""
    worst_case_return(decision::Array{Float64,1}, formulation::PortfolioOpt.RobustBertsimas, solver)

Returns worst case return (WCR) in Bertsimas's uncertainty set ([`RobustBertsimas`](@ref)) for a defined decision:

$(PortfolioOpt._portfolio_return_latex_RobustBertsimas_primal())

"""
function worst_case_return(decision::Array{Float64,1}, formulation::PortfolioOpt.RobustBertsimas, solver)
    r̄ = formulation.predicted_mean
    numA = formulation.number_of_assets
    Δ = formulation.uncertainty_delta
    Λ = formulation.bertsimas_budget

    model = Model(solver)
    @variable(model, μ[i=1:numA])
    @variable(model, 0.0 <= z[i=1:numA] <= 1)
    @constraint(model, [i=1:numA], μ[i] <= r̄[i] + z[i] * Δ[i])
    @constraint(model, [i=1:numA], μ[i] >= r̄[i] - z[i] * Δ[i])
    @constraint(model, sum(z) <= Λ)
    @objective(model, Min, decision'μ)

    optimize!(model)
    status = termination_status(model)
    status !== MOI.OPTIMAL && @warn "Did not find an optimal solution: status=$status"

    return objective_value(model)
end

function best_case_return(decision::Array{Float64,1}, formulation::PortfolioOpt.RobustBertsimas, solver)
    r̄ = formulation.predicted_mean
    numA = formulation.number_of_assets
    Δ = formulation.uncertainty_delta
    Λ = formulation.bertsimas_budget

    model = Model(solver)
    @variable(model, μ[i=1:numA])
    @variable(model, 0.0 <= z[i=1:numA] <= 1)
    @constraint(model, [i=1:numA], μ[i] <= r̄[i] + z[i] * Δ[i])
    @constraint(model, [i=1:numA], μ[i] >= r̄[i] - z[i] * Δ[i])
    @constraint(model, sum(z) <= Λ)
    @objective(model, Max, decision'μ)

    optimize!(model)
    status = termination_status(model)
    status !== MOI.OPTIMAL && @warn "Did not find an optimal solution: status=$status"

    return objective_value(model)
end

"""
    worst_case_return(decision::Array{Float64,1}, formulation::PortfolioOpt.RobustBenTal, solver)

Returns worst case return (WCR) in BenTal's uncertainty set ([`RobustBenTal`](@ref)) for a defined decision:

$(PortfolioOpt._portfolio_return_latex_RobustBenTal_primal())

"""
function worst_case_return(decision::Array{Float64,1}, formulation::PortfolioOpt.RobustBenTal, solver)
    r̄ = formulation.predicted_mean
    sqrt_Σ_inv = inv(sqrt(formulation.predicted_covariance))
    numA = formulation.number_of_assets
    δ = formulation.uncertainty_delta

    model = Model(solver)
    @variable(model, μ[i=1:numA])
    @constraint(model, [δ; sqrt_Σ_inv * (μ - r̄)] in JuMP.SecondOrderCone())
    @objective(model, Min, decision'μ)

    optimize!(model)
    status = termination_status(model)
    status !== MOI.OPTIMAL && @warn "Did not find an optimal solution: status=$status"

    return objective_value(model)
end

function best_case_return(decision::Array{Float64,1}, formulation::PortfolioOpt.RobustBenTal, solver)
    r̄ = formulation.predicted_mean
    sqrt_Σ_inv = inv(sqrt(formulation.predicted_covariance))
    numA = formulation.number_of_assets
    δ = formulation.uncertainty_delta

    model = Model(solver)
    @variable(model, μ[i=1:numA])
    @constraint(model, [δ; sqrt_Σ_inv * (μ - r̄)] in JuMP.SecondOrderCone())
    @objective(model, Max, decision'μ)

    optimize!(model)
    status = termination_status(model)
    status !== MOI.OPTIMAL && @warn "Did not find an optimal solution: status=$status"

    return objective_value(model)
end

################################### To be Updated #####################################
function returns_montecarlo(Σ, r̄, numS)
    d = MvNormal(r̄, Σ)
    r = rand(d, numS)
    P = pdf(d, r)
    return r', P / sum(P)
end

function compute_solution_stoc(model::JuMP.Model, w; solver=DEFAULT_SOLVER)
    set_optimizer(model, solver)
    optimize!(model)
    status = termination_status(model)
    status !== MOI.OPTIMAL && @warn "Did not find an optimal solution: status=$status"

    w_values = value.(w)
    if sum(w_values) > 1.0
        w_values = w_values / sum(w_values)
    end
    Cvar = -objective_value(model)
    r = sum(r̄'w_values)
    q1_α = value.(model[:z])
    return w_values, r, Cvar, q1_α
end

function compute_solution_stoc_2(model::JuMP.Model, w; solver=DEFAULT_SOLVER)
    set_optimizer(model, solver)
    optimize!(model)
    status = termination_status(model)
    status !== MOI.OPTIMAL && @warn "Did not find an optimal solution: status=$status"

    w_values = value.(w)
    if sum(w_values) > 1.0
        w_values = w_values / sum(w_values)
    end
    r = objective_value(model)
    q1_α = value.(model[:z])
    return w_values, r, q1_α
end

ext: jl

max_stars_repo_name: andrewrosemberg/PortfolioOpt

max_issues_repo_name: andrewrosemberg/PortfolioOpt

max_forks_repo_name: andrewrosemberg/PortfolioOpt


text:
#filepath cross system compatability code
if Sys.KERNEL === :NT || Sys.KERNEL === :Windows
    pathmark = "\\";
else
    pathmark = "/";
end

include("concordance.jl")

numSectors = 20
# for the 20 sector model
IOIGAs20=Array{Union{Nothing, String}}(nothing, length(IOIG));
for i in eachindex(IOIG);
    IOIGAs20[i] = IOIGTo20[IOIG[i]]
end

# reducing sectors by column
IOHalf = DataFrame(copy(transpose(IOSource[[4:117; 119; 121:126; 128], 3:116])),:auto);
insertcols!(IOHalf ,1, :Industry => IOIGAs20);
IOHalfSplitIndustry = groupby(IOHalf, :Industry);
IOHalf = combine(IOHalfSplitIndustry, valuecols(IOHalfSplitIndustry) .=> sum);

sort!(IOHalf)
IOHalf = select!(IOHalf, Not(:Industry));
IOHalf = transpose(Matrix(IOHalf[:,:]));

# reducing sectors by row
IO20 = DataFrame(copy([IOHalf[1:114, :] IOSource[4:117,117:126]]),:auto);
insertcols!(IO20 ,1, :Industry => IOIGAs20);
IO20SplitIndustry = groupby(IO20, :Industry);
IO20 = combine(IO20SplitIndustry, valuecols(IO20SplitIndustry) .=> sum);
IO20 = select!(IO20, Not(:Industry));
IO20 = Matrix(IO20[:,:])

reducedIO = [IO20; [IOHalf[115:122,:] IOSource[[119; 121:126; 128], 117:126]]];

reducedCapFlows = Matrix(select!(CSV.read("data"*pathmark*"capitalFlowsDraft.csv", DataFrame), Not(:Titles)));

max_stars_repo_name: uq-aibe/aim

max_issues_repo_name: uq-aibe/aim

max_forks_repo_name: uq-aibe/aim


text:
# BFS algorithm using cvdlab SuiteSparseGraphBLAS Julia wrapper
# (can be found at https://github.com/cvdlab/SuiteSparseGraphBLAS.jl)
# The following function is very similar to the one at
# https://github.com/cvdlab/SuiteSparseGraphBLAS.jl/blob/master/examples/bfs.jl
# All credits go to the original author.
using SuiteSparseGraphBLAS

"""
    GraphBLAS_bfs_cvd(A::GBMatrix, source::Integer)

Compute BFS on the graph represented by `A` starting from `source`
using GraphBLAS CVD version.

# Examples
```julia-repl
julia> nrows, ncols = 100, 100;

julia> A = from_type(Int64, nrows, ncols);

julia> GraphBLAS_sparse_square_matrix_cvd!(A);

julia> source = 1;

julia> GraphBLAS_bfs_cvd(A, source);
```
"""
function GraphBLAS_bfs_cvd(
    A::GBMatrix,
    source::Integer
    )

    source > 0 || throw(ArgumentError("source must be positive"))

    nrows = size(A,1)
    size(A, 2) == nrows || throw(DimensionMismatch("Input matrix must be square"))

    q = from_lists([source], [true], n = nrows)
    desc = descriptor(outp=>SuiteSparseGraphBLAS.replace, mask=>scmp)
    v = from_type(Int64, nrows)

    for i in 1:nrows-1
        @with q v[:] = i
        vxm(q, A, out = q, semiring = Semirings.LOR_LAND, mask = v, desc = desc)
        if !reduce(q, monoid = Monoids.LOR)
            break
        end
    end

    return v
end


text:
# Geometric Utilities

export show_cs

show_cs(p, scale=1) =
    let rcyl = scale/10
        rcon = scale/5
        lcyl = scale
        lcon = scale/5
        px = add_x(p, 3*lcyl)
        py = add_y(p, 2*lcyl)
        pz = add_z(p, 1*lcyl)
        union(cylinder(p, rcyl, px),
              cone(px, rcon, add_x(px, lcon)),
              cylinder(p, rcyl, py),
              cone(py, rcon, add_y(py, lcon)),
              cylinder(p, rcyl, pz))
    end

project_to_world(surf) =
    transform(surf, inverse_transformation(frame_at(surf, 0, 0)))

#project_to_world(surface_polygon(xyz(1,1,1), xyz(10,1,1), xyz(10,1,5), xyz(1,1,5)))

#=

Given a poligonal line described by its vertices, we need to compute another
polygonal line that is parallel to the first one.

=#

v_in_v(v0, v1) =
  let v = v0 + v1
    v*dot(v0, v0)/dot(v, v0)
  end

rotated_v(v, alpha) =
  vpol(pol_rho(v), pol_phi(v) + alpha)

centered_rectangle(p0, w, p1) =
  let v0 = p1 - p0,
      v1 = rotated_v(v0, pi/2),
      c = loc_from_o_vx_vy(p0, v0, v1)
    rectangle(c-vy(w/2, c.cs), distance(p0, p1), w)
  end

offset_vertices(ps::Locs, d::Real, closed) =
  let qs = closed ? [ps[end], ps..., ps[1]] : ps,
      vs = map((p0, p1) -> rotated_v(unitized(p1 - p0)*d, pi/2), qs[1:end-1], qs[2:end]),
      ws = map(v_in_v, vs[1:end-1], vs[2:end])
    map(+, ps, closed ? ws : [vs[1], ws..., vs[end]])
  end

offset(path::Union{Path,Shape}, d::Real) = d == 0 ? path : nonzero_offset(path, d)
nonzero_offset(path::RectangularPath, d::Real) =
  rectangular_path(add_xy(path.corner, d, d), path.dx - 2d, path.dy - 2d)
nonzero_offset(path::OpenPolygonalPath, d::Real) =
  open_polygonal_path(offset_vertices(path.vertices, d, false))
nonzero_offset(path::ClosedPolygonalPath, d::Real) =
  closed_polygonal_path(offset_vertices(path.vertices, d, true))
nonzero_offset(l::Line, d::Real) =
  line(offset(l.vertices, d, false))
nonzero_offset(path::CircularPath, d::Real) =
  circular_path(path.center, path.radius - d)
nonzero_offset(path::ArcPath, d::Real) =
  arc_path(path.center, path.radius - d, path.start_angle, path.amplitude)

export offset

# Polygon combination

closest_vertices_indexes(pts1, pts2) =
  # This is a brute force method. There are better algorithms to do this.
  let min_dist = Inf,
      min_i = nothing,
      min_j = nothing
    for (i, pt1) in enumerate(pts1)
      for (j, pt2) in enumerate(pts2)
        let dist = distance(pt1, pt2)
          if dist < min_dist
            min_dist = dist
            min_i = i
            min_j = j
          end
        end
      end
    end
    min_i, min_j
  end

#=
using Test
pts1 = regular_polygon_vertices(4, xy(1,2), 4)
pts2 = regular_polygon_vertices(4, xy(2,2), 1)
pts3 = regular_polygon_vertices(5, xy(1,4), 1)
pts4 = regular_polygon_vertices(6, xy(0,1), 1)

polygon(pts1)
polygon(pts2)
polygon(pts3)
polygon(pts4)

@test closest_vertices_indexes(pts1, pts2) == (1,1)
@test closest_vertices_indexes(pts1, pts3) == (2,2)
@test closest_vertices_indexes(pts1, pts4) == (4,6)
@test closest_vertices_indexes(pts2, pts3) == (2,5)
=#

point_in_segment(r, p, q) =
  let pr = r-p,
      pq = q-p,
      rx = cx(pr)/cx(pq),
      ry = cy(pr)/cy(pq)
    isapprox(rx, ry) && isapprox(ry, cz(pr)/cz(pq))
  end

collinear_segments(p1, p2, q1, q2) =
  point_in_segment(q1, p1, p2) &&
  point_in_segment(q2, p1, p2)

collinear_vertices_indexes(pts1, pts2) =
  for (i1, p1) in enumerate(pts1)
    let i2 = (i1+1)%length(pts1),
        p2 = pts1[i2]
      for (j1, q1) in enumerate(pts2)
        let j2 = (j1+1)%length(pts2),
            q2 = pts2[j2]
          if collinear_segments(p1, p2, q1, q2)
            return (i1, j1)
          end
        end
      end
    end
  end

subtract_polygon_vertices(pts1, pts2) =
  let ij = collinear_vertices_indexes(pts1, pts2)
    isnothing(ij) ?
      inject_polygon_vertices_at_indexes(pts1, pts2, closest_vertices_indexes(pts1, pts2)) :
      splice_polygon_vertices_at_indexes(pts1, pts2, ij)
  end

inject_polygon_vertices_at_indexes(pts1, pts2, (i, j)) =
  [pts1[1:i]..., reverse([pts2[j:end]..., pts2[1:j]...])..., pts1[i:end]...]

export closest_vertices_indexes, inject_polygon_vertices_at_indexes, subtract_polygon_vertices

#=
pts1 = subtract_polygon_vertices(pts1, pts2)
pts1 = subtract_polygon_vertices(pts1, pts3)
pts1 = subtract_polygon_vertices(pts1, pts4)
polygon(pts1)
for (i,p) in enumerate(pts1)
  text(string(i), p, 0.1)
  sleep(1)
end
=#

# Intersection

segments_intersection(p0, p1, p2, p3) =
  let denom = (p3.y - p2.y)*(p1.x - p0.x) - (p3.x - p2.x)*(p1.y - p0.y)
    if denom == 0
      nothing
    else
      let u = ((p3.x - p2.x)*(p0.y - p2.y) - (p3.y - p2.y)*(p0.x - p2.x))/denom,
          v = ((p1.x - p0.x)*(p0.y - p2.y) - (p1.y - p0.y)*(p0.x - p2.x))/denom
        if 0 <= u <= 1 && 0 <= v <= 1
          xy(p0.x + u*(p1.x - p0.x), p0.y + u*(p1.y - p0.y))
        else
          nothing
        end
      end
    end
  end

lines_intersection(p0, p1, p2, p3) =
  let denom = (p3.y - p2.y)*(p1.x - p0.x) - (p3.x - p2.x)*(p1.y - p0.y)
    if denom == 0
      nothing
    else
      let u = ((p3.x - p2.x)*(p0.y - p2.y) - (p3.y - p2.y)*(p0.x - p2.x))/denom,
          v = ((p1.x - p0.x)*(p0.y - p2.y) - (p1.y - p0.y)*(p0.x - p2.x))/denom
        xy(p0.x + u*(p1.x - p0.x), p0.y + u*(p1.y - p0.y))
      end
    end
  end

const collinearity_tolerance = Parameter(1e-2)
# are the three points sufficiently collinear?
collinear_points(p0, pm, p1, epsilon=collinearity_tolerance()) =
  let a = distance(p0, pm),
      b = distance(pm, p1),
      c = distance(p1, p0)
    triangle_area(a, b, c) < epsilon
  end

triangle_area(a, b, c) =
  let s = (a + b + c)/2.0
    sqrt(s*(s - a)*(s - b)*(s - c))
  end


text:
# Makie

# without consistency bounds
function Makie.convert_arguments(
    ::Type{<:Reliability},
    probabilities::AbstractVector{<:Real},
    frequencies::AbstractVector{<:Real},
)
    points = map(probabilities, frequencies) do x, y
        return Makie.Point4f0(x, y, NaN32, NaN32)
    end
    return (points,)
end

# with separate consistency bounds
function Makie.convert_arguments(
    ::Type{<:Reliability},
    probabilities::AbstractVector{<:Real},
    frequencies::AbstractVector{<:Real},
    low::AbstractVector{<:Real},
    high::AbstractVector{<:Real},
)
    points = map(probabilities, frequencies, low, high) do x, y, low, high
        return Makie.Point4f0(x, y, low, high)
    end
    return (points,)
end

# with tuples of consistency bounds
function Makie.convert_arguments(
    ::Type{<:Reliability},
    probabilities::AbstractVector{<:Real},
    frequencies::AbstractVector{<:Real},
    low_high::AbstractVector,
)
    points = map(probabilities, frequencies, low_high) do x, y, (low, high)
        return Makie.Point4f0(x, y, low, high)
    end
    return (points,)
end

# with outcomes instead of frequencies
function Makie.used_attributes(
    ::Type{<:Reliability}, ::AbstractVector{<:Real}, ::AbstractVector{<:Bool}
)
    return (:binning, :consistencybars)
end
function Makie.convert_arguments(
    ::Type{P},
    probabilities::AbstractVector{<:Real},
    outcomes::AbstractVector{Bool};
    binning=EqualMass(),
    consistencybars=ConsistencyBars(),
) where {P<:Reliability}
    meanprobabilities, meanfrequencies, low_high = means_and_bars(
        probabilities, outcomes; binning=binning, consistencybars=consistencybars
    )
    return Makie.convert_arguments(P, meanprobabilities, meanfrequencies, low_high)
end

# Plots

# default constructor (recipe forwards arguments as tuple)
function ReliabilityPlot(
    (
        probabilities, frequencies, low_high
    )::Tuple{
        <:AbstractVector{<:Real},
        <:AbstractVector{<:Real},
        <:AbstractVector{<:Tuple{<:Real,<:Real}},
    },
)
    return ReliabilityPlot(probabilities, frequencies, low_high)
end
function ReliabilityPlot(
    (
        probabilities, frequencies, low, high
    )::Tuple{
        <:AbstractVector{<:Real},
        <:AbstractVector{<:Real},
        <:AbstractVector{<:Real},
        <:AbstractVector{<:Real},
    },
)
    return ReliabilityPlot((probabilities, frequencies, map(tuple, low, high)))
end

# without consistency bars
function ReliabilityPlot(
    (
        probabilities, frequencies_or_outcomes
    )::Tuple{<:AbstractVector{<:Real},<:AbstractVector{<:Real}},
)
    return ReliabilityPlot(probabilities, frequencies_or_outcomes, nothing)
end


text:
module Underlying

export xunderlying

include("xunderlying.jl")

end


text:
struct IncrementalSUniward{B,T<:AbstractFloat} <: AbstractImageDistortion
    image::Array{UInt8}
    invrs::Array{B,1}
    waveletplanes::Array{B,1}
    distortion::T
    height::Int
    width::Int
end

Base.similar(i::IncrementalSUniward) = IncrementalSUniward(deepcopy(i.image),i.invrs,deepcopy(i.waveletplanes),i.distortion,i.height,i.width)
IncrementalSUniward(cover::Array{UInt8};sigma=1.0,T::Type=Float64) = IncrementalSUniward(cover,1,sigma=sigma,T=T)[1]

function IncrementalSUniward(cover::Array{UInt8},n::Int;sigma=1.0,T::Type=Float64)
    paddedcover = padarray(T.(cover),Pad(:symmetric,[16,16],[16,16]))
    invrs=Array{typeof(paddedcover),1}(3);
    for f in 1:length(invrs)
        invrs[f] = padarray(zeros(T,size(cover)),Pad(:symmetric,[16,16],[16,16]))
        invrs[f].parent .= 1 ./ (abs.(sameconv2(paddedcover.parent, flts[f])) + T(sigma))
    end
    mn = min(map(minimum,invrs)...)
    mx = max(map(maximum,invrs)...)
    [IncrementalSUniward(deepcopy(cover),
        invrs,
        [padarray(zeros(T,size(cover)),Pad(:symmetric,[16,16],[16,16])) for i in 1:3],
        0.0,size(cover,1),size(cover,2)) for i in 1:n]
end

function setindex!(image::IncrementalSUniward,v,i)
    setindex!(image,v,ind2sub((image.height,image.width),i)...)
end

function setindex!(image::IncrementalSUniward,v,i,j)
    (v>255 || v<0) && error("Value $v outside the allowed range [0,255]")
    idxs=(max(1,i-8):min(image.height,i+7),max(1,j-8):min(image.width,j+7))
    delta=eltype(image.waveletplanes[1])(v)-eltype(image.waveletplanes[1])(image.image[i,j])
    for f in 1:3
        A=image.waveletplanes[f]
        invr=image.invrs[f]
        flt=rotflts[f]
        for l in max(1,10-j):min(image.width+9-j,16)
            @simd for k in max(1,10-i):min(image.height+9-i,16)
                @inbounds image.distortion-=abs.(A[i-9+k,j-9+l])
                δ = delta*flt[k,l]*invr[i-9+k,j-9+l]
                @inbounds A[i-9+k,j-9+l]+=δ
                @inbounds image.distortion+=abs.(A[i-9+k,j-9+l])
            end
        end
    end
    image.image[i,j]=v
end

function tryvalue(image::IncrementalSUniward,v,i)
    tryvalue(image,v,ind2sub((image.height,image.width),i)...)
end

function tryvalue(image::IncrementalSUniward{B,T},v,i,j) where {B,T}
    if v>255 || v<0
        return(typemax(T))
    end
    distortion=image.distortion
    delta=v-image.image[i,j]
    for f in 1:3
        A=image.waveletplanes[f]
        invr=image.invrs[f]
        flt=rotflts[f]
        for l in max(1,10-j):min(image.width+9-j,16)
            @simd for k in max(1,10-i):min(image.height+9-i,16)
                @inbounds distortion-=abs.(A[i-9+k,j-9+l])
                @inbounds distortion+=abs.(A[i-9+k,j-9+l]+delta*flt[k,l]*invr[i-9+k,j-9+l])
            end
        end
    end
    return(distortion)
end

function trypmone(image::IncrementalSUniward,i)
    trypmone(image,ind2sub((image.height,image.width),i)...)
end

function trypmone(image::IncrementalSUniward,i,j)
    dp=0.0
    dm=0.0
    d=image.distortion
    for f in 1:3
        A=image.waveletplanes[f]
        invr=image.invrs[f]
        flt=rotflts[f]
        for l in max(1,10-j):min(image.width+9-j,16)
            @simd for k in max(1,10-i):min(image.height+9-i,16)
                @inbounds d-=abs.(A[i-9+k,j-9+l])
                @inbounds dp+=abs.(A[i-9+k,j-9+l]+flt[k,l]*invr[i-9+k,j-9+l])
                @inbounds dm+=abs.(A[i-9+k,j-9+l]-flt[k,l]*invr[i-9+k,j-9+l])
            end
        end
    end
    dp = (image[i,j] == 255) ? typemax(Float64) : dp
    dm = (image[i,j] == 0) ? typemax(Float64) : dm
    return(d+dp,d+dm)
end

function copy!(dest::IncrementalSUniward,src::IncrementalSUniward)
    if dest.width!=src.width || dest.height!=src.height
        error("source and destination image has to have the same height")
    end
    copy!(dest.image,src.image)
    for f in 1:3
        dest.invrs[f]=src.invrs[f]
        copy!(dest.waveletplanes[f].parent,src.waveletplanes[f].parent)
    end
    dest.distortion=src.distortion
    return(dest)
end


text:
using BinDeps

@BinDeps.setup

windllname = "IpOpt-vc10"
libipopt = library_dependency("libipopt", aliases=[windllname])

ipoptname = "Ipopt-3.11.7"

provides(Sources, URI("http://www.coin-or.org/download/source/Ipopt/$ipoptname.tgz"),
    libipopt, os = :Unix)

prefix=joinpath(BinDeps.depsdir(libipopt),"usr")
patchdir=BinDeps.depsdir(libipopt)
srcdir = joinpath(BinDeps.depsdir(libipopt),"src",ipoptname)

provides(SimpleBuild,
    (@build_steps begin
        GetSources(libipopt)
        @build_steps begin
            ChangeDirectory(srcdir)
            @build_steps begin
                ChangeDirectory(joinpath(srcdir,"ThirdParty","Blas"))
                CreateDirectory("build", true)
                `./get.Blas`
            end
            @build_steps begin
                ChangeDirectory(joinpath(srcdir,"ThirdParty","Blas","build"))
                `../configure --prefix=$prefix --disable-shared --with-pic`
                `make install`
            end
            @build_steps begin
                ChangeDirectory(joinpath(srcdir,"ThirdParty","Lapack"))
                CreateDirectory("build", true)
                `./get.Lapack`
            end
            @build_steps begin
                ChangeDirectory(joinpath(srcdir,"ThirdParty","Lapack","build"))
                `../configure --prefix=$prefix --disable-shared --with-pic`
                `make install`
            end
            @build_steps begin
                ChangeDirectory(joinpath(srcdir,"ThirdParty","Mumps"))
                `./get.Mumps`
            end
            `./configure --prefix=$prefix --enable-dependency-linking --with-blas=$prefix/lib/libcoinblas.a --with-lapack=$prefix/lib/libcoinlapack.a`
            `make install`
        end
    end),libipopt, os = :Unix)

# OS X
@osx_only begin
    using Homebrew
    provides(Homebrew.HB, "ipopt", libipopt, os = :Darwin)
end

# Windows
downloadsdir = BinDeps.downloadsdir(libipopt)
libdir = BinDeps.libdir(libipopt)
downloadname = "Ipopt-3.11.0-Win32-Win64-dll.7z"
windir = WORD_SIZE == 32 ? "Win32" : "x64"

# BinDeps complains about the .7z file on other platforms...
@windows_only provides(BuildProcess,
    (@build_steps begin
        FileDownloader("http://www.coin-or.org/download/binary/Ipopt/$downloadname", joinpath(downloadsdir, downloadname))
	CreateDirectory(BinDeps.srcdir(libipopt), true)
	FileUnpacker(joinpath(downloadsdir, downloadname), BinDeps.srcdir(libipopt), joinpath(BinDeps.srcdir(libipopt),"lib","Win32"))
	CreateDirectory(libdir, true)
	@build_steps begin
	    ChangeDirectory(joinpath(BinDeps.srcdir(libipopt),"lib",windir,"ReleaseMKL"))
	    FileRule(joinpath(libdir,"$(windllname).dll"), @build_steps begin
	        `cp *.dll $(libdir)`
	    end)
	end
     end), libipopt, os = :Windows)

@windows_only push!(BinDeps.defaults, BuildProcess)

@BinDeps.install [:libipopt => :libipopt]

@windows_only pop!(BinDeps.defaults)


text:
using ModelingToolkit, SparseArrays, BenchmarkTools,
    LinearAlgebra

# Generate RAM matrices of factor models
function gen_model(nfact, nitem)
    nfact = Int64(nfact)
    nitem = Int64(nitem)

    ## Model definition
    nobs = Int64(nfact*nitem)
    nnod = nfact+nobs
    @ModelingToolkit.variables λ[1:nobs], σ[1:nobs]

    #F
    Ind = collect(1:nobs)
    Jnd = collect(1:nobs)
    V = fill(1,nobs)
    F = sparse(Ind, Jnd, V, nobs, nnod)

    #S
    Ind = collect(1:nnod)
    Jnd = collect(1:nnod)
    V = [σ; fill(1.0, nfact)]
    S = sparse(Ind, Jnd, V, nnod, nnod)

    #A
    Ind = collect(1:nobs)
    Jnd = vcat([fill(nobs+i, nitem) for i in 1:nfact]...)
    V = λ
    A = sparse(Ind, Jnd, V, nnod, nnod)
    return(A, S, F)
end

function neumann_series(mat::SparseMatrixCSC)
    inverse = I + mat
    next_term = mat^2

    while nnz(next_term) != 0
        inverse += next_term
        next_term *= mat
    end

    return inverse
end

function implied_cov(A, S, F)
    invia = neumann_series(A)

    imp_cov_sym = F*invia*S*permutedims(invia)*permutedims(F)

    imp_cov_sym = Array(imp_cov_sym)
    imp_cov_sym = ModelingToolkit.simplify.(imp_cov_sym)
    return imp_cov_sym
end

function symbolic_ULS(Sigma, S)
    diff = S - Sigma
    F = sum(x -> ^(x, 2), diff)
    #F = diff*permutedims(diff)
    #F = F[diagind(F)]
    F = simplify.(F)
    return diff
end

# small factor models
A, S, F = gen_model(1, 4)

Sigma = implied_cov(A, S, F)

nobs = size(Sigma, 1)

@ModelingToolkit.variables s[1:nobs, 1:nobs]

F = symbolic_ULS(Sigma, s)


text:
@testset "unbiased.jl" begin
    @testset "Unbiased: Two-dimensional example" begin
        # categorical distributions
        skce = UnbiasedSKCE(SqExponentialKernel() ⊗ WhiteKernel())
        for predictions in ([[1, 0], [0, 1]], ColVecs([1 0; 0 1]), RowVecs([1 0; 0 1]))
            @test iszero(@inferred(skce(predictions, [1, 2])))
            @test iszero(@inferred(skce(predictions, [1, 1])))
            @test @inferred(skce(predictions, [2, 1])) ≈ -2 * exp(-1)
            @test iszero(@inferred(skce(predictions, [2, 2])))
        end

        # probabilities
        skce = UnbiasedSKCE(
            (SqExponentialKernel() ∘ ScaleTransform(sqrt(2))) ⊗ WhiteKernel()
        )
        @test iszero(@inferred(skce([1, 0], [true, false])))
        @test iszero(@inferred(skce([1, 0], [true, true])))
        @test @inferred(skce([1, 0], [false, true])) ≈ -2 * exp(-1)
        @test iszero(@inferred(skce([1, 0], [false, false])))
    end

    @testset "Unbiased: Basic properties" begin
        skce = UnbiasedSKCE((ExponentialKernel() ∘ ScaleTransform(0.1)) ⊗ WhiteKernel())
        estimates = Vector{Float64}(undef, 1_000)

        # categorical distributions
        for nclasses in (2, 10, 100)
            dist = Dirichlet(nclasses, 1.0)

            predictions = [Vector{Float64}(undef, nclasses) for _ in 1:20]
            targets = Vector{Int}(undef, 20)

            for i in 1:length(estimates)
                rand!.(Ref(dist), predictions)
                targets .= rand.(Categorical.(predictions))

                estimates[i] = skce(predictions, targets)
            end

            @test any(x -> x > zero(x), estimates)
            @test any(x -> x < zero(x), estimates)
            @test mean(estimates) ≈ 0 atol = 1e-3
        end

        # probabilities
        predictions = Vector{Float64}(undef, 20)
        targets = Vector{Bool}(undef, 20)
        for i in 1:length(estimates)
            rand!(predictions)
            map!(targets, predictions) do p
                rand() < p
            end
            estimates[i] = skce(predictions, targets)
        end

        @test any(x -> x > zero(x), estimates)
        @test any(x -> x < zero(x), estimates)
        @test mean(estimates) ≈ 0 atol = 1e-3
    end

    @testset "Block: Two-dimensional example" begin
        # categorical distributions
        skce = BlockUnbiasedSKCE(SqExponentialKernel() ⊗ WhiteKernel())
        for predictions in ([[1, 0], [0, 1]], ColVecs([1 0; 0 1]), RowVecs([1 0; 0 1]))
            @test iszero(@inferred(skce(predictions, [1, 2])))
            @test iszero(@inferred(skce(predictions, [1, 1])))
            @test @inferred(skce(predictions, [2, 1])) ≈ -2 * exp(-1)
            @test iszero(@inferred(skce(predictions, [2, 2])))
        end

        # two predictions, ten times replicated
        for predictions in (
            repeat([[1, 0], [0, 1]], 10),
            ColVecs(repeat([1 0; 0 1], 1, 10)),
            RowVecs(repeat([1 0; 0 1], 10, 1)),
        )
            @test iszero(@inferred(skce(predictions, repeat([1, 2], 10))))
            @test iszero(@inferred(skce(predictions, repeat([1, 1], 10))))
            @test @inferred(skce(predictions, repeat([2, 1], 10))) ≈ -2 * exp(-1)
            @test iszero(@inferred(skce(predictions, repeat([2, 2], 10))))
        end

        # probabilities
        skce = BlockUnbiasedSKCE(
            (SqExponentialKernel() ∘ ScaleTransform(sqrt(2))) ⊗ WhiteKernel()
        )
        @test iszero(@inferred(skce([1, 0], [true, false])))
        @test iszero(@inferred(skce([1, 0], [true, true])))
        @test @inferred(skce([1, 0], [false, true])) ≈ -2 * exp(-1)
        @test iszero(@inferred(skce([1, 0], [false, false])))

        # two predictions, ten times replicated
        @test iszero(@inferred(skce(repeat([1, 0], 10), repeat([true, false], 10))))
        @test iszero(@inferred(skce(repeat([1, 0], 10), repeat([true, true], 10))))
        @test @inferred(skce(repeat([1, 0], 10), repeat([false, true], 10))) ≈ -2 * exp(-1)
        @test iszero(@inferred(skce(repeat([1, 0], 10), repeat([false, false], 10))))
    end

    @testset "Block: Basic properties" begin
        nsamples = 20
        kernel = (ExponentialKernel() ∘ ScaleTransform(0.1)) ⊗ WhiteKernel()
        skce = UnbiasedSKCE(kernel)
        blockskce = BlockUnbiasedSKCE(kernel)
        blockskce_all = BlockUnbiasedSKCE(kernel, nsamples)
        estimates = Vector{Float64}(undef, 1_000)

        # categorical distributions
        for nclasses in (2, 10, 100)
            dist = Dirichlet(nclasses, 1.0)

            predictions = [Vector{Float64}(undef, nclasses) for _ in 1:nsamples]
            targets = Vector{Int}(undef, nsamples)

            for i in 1:length(estimates)
                rand!.(Ref(dist), predictions)
                targets .= rand.(Categorical.(predictions))

                estimates[i] = blockskce(predictions, targets)

                # consistency checks
                @test estimates[i] ≈ mean(
                    skce(predictions[(2 * i - 1):(2 * i)], targets[(2 * i - 1):(2 * i)]) for
                    i in 1:(nsamples ÷ 2)
                )
                @test skce(predictions, targets) == blockskce_all(predictions, targets)
            end

            @test any(x -> x > zero(x), estimates)
            @test any(x -> x < zero(x), estimates)
            @test mean(estimates) ≈ 0 atol = 5e-3
        end

        # probabilities
        predictions = Vector{Float64}(undef, nsamples)
        targets = Vector{Bool}(undef, nsamples)

        for i in 1:length(estimates)
            rand!(predictions)
            map!(targets, predictions) do p
                return rand() < p
            end
            estimates[i] = blockskce(predictions, targets)

            # consistency checks
            @test estimates[i] ≈ mean(
                skce(predictions[(2 * i - 1):(2 * i)], targets[(2 * i - 1):(2 * i)]) for
                i in 1:(nsamples ÷ 2)
            )
            @test skce(predictions, targets) == blockskce_all(predictions, targets)
        end

        @test any(x -> x > zero(x), estimates)
        @test any(x -> x < zero(x), estimates)
        @test mean(estimates) ≈ 0 atol = 5e-3
    end
end


text:
using GaussQuadrature, Plots
using ClimateMachine.Mesh.Elements
using ClimateMachine
gr(size = (400,400))
n = 4
r,w = GaussQuadrature.legendre(n+1, both)
T = eltype(r)
a, b = GaussQuadrature.legendre_coefs(T, n)
x = collect(range(-1,1,length = 100))
V = GaussQuadrature.orthonormal_poly(x, a, b)
p = []
push!(p,plot(x, V[:,1], ylims = (-2,2),
    title = string(1-1) * "'th",
    legend = false)
)
for i in 2:n+1
    push!(p,
    plot(x, V[:,i], ylims = (-2,2),
    title = string(i-1) * "'th",
    legend = false)
    )
end
display(plot(p...))
## Show decomposition
f(x) = sin(π*x) * exp(x) / 2
ry = @. f(r)
xy = @. f(x)
plot(r,ry)
plot!(x, xy)
𝒯 = GaussQuadrature.orthonormal_poly(r, a, b)
spectrum = inv(𝒯) * f.(r)
index = collect(range(0, n,length = n+1))
scatter(index, log10.(abs.(spectrum)),
xlabel = "Legendre Index",
ylabel = "Log10(Abs(Amplitude))",
label = false,
title = "Legendre Modal Amplitudes")
##
function plot_approximation(V, spectrum, i; exact = true)
    p1 = plot(x, V[:,1:i] * spectrum[1:i], ylims = (-1.1,1.1),
        title = "Modes 0 through " * string(i-1),
        label = "approximation", legend = :bottomright)
    if exact
        plot!(x, xy, color = :red, label = "truth")
    end
    return p1
end

p = []
for i in 1:n+1
    push!(p,
        plot_approximation(V, spectrum, i)
    )
end
display(plot(p...))
for i in 1:n+1
    display(p[i])
    sleep(0.3)
end

## Mass Matrix and Stiffness Matrix
Mi =  𝒯 * 𝒯'
M = inv(Mi)
D = ClimateMachine.Mesh.Elements.spectralderivative(r)
Sᵀ = (M * D)'

## Show that derivatives are not accurate for polynomials of order n+1
# but accurate for polynomials below that order
L = hcat([r .^ i for i in 1:n+1]...)
DL = hcat([i * r .^(i-1) for i in 1:n+1]...)
D * L - DL
## Compute Differentiation Matrix (numerically unstable for large N)
L  = hcat([r .^ i for i in 0:n]...)
DL = hcat([i * r .^ (i-1) for i in 0:n]...)
D*L - DL
D - DL / L # check computation of differentiation matrix

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
using Test

@testset "Quiqbox tests" begin

    unit1 = "Support Functions"
    println("Testing $(unit1)...")
    t1 = @elapsed @testset "$(unit1)" begin
        include("unit-tests/Tools-test.jl")
        include("unit-tests/FileIO-test.jl")
        include("unit-tests/Overload-test.jl")
    end
    println("$(unit1) test finished in $t1 seconds.\n")

    @testset "Module Functions" begin

        unit2_1 = "Main Functions"
        println("Testing $(unit2_1)...")
        t2_1 = @elapsed @testset "$(unit2_1)" begin
            include("unit-tests/Integration/OneBody-test.jl")
            include("unit-tests/Integration/TwoBody-test.jl")
            include("unit-tests/Basis-test.jl")
            include("unit-tests/Differentiation-test.jl")
        end
        println("$(unit2_1) test finished in $t2_1 seconds.\n")

        unit2_2 = "Data Structure"
        println("Testing $(unit2_2)...")
        t2_2 = @elapsed @testset "$(unit2_2)" begin
            include("unit-tests/Box-test.jl")
            include("unit-tests/Library-test.jl")
            include("unit-tests/Molecule-test.jl")
        end
        println("$(unit2_2) test finished in $t2_2 seconds.\n")

        unit2_3 = "Applications"
        println("Testing $(unit2_3)...")
        t2_3 = @elapsed @testset "$(unit2_3)" begin
            include("unit-tests/HartreeFock-test.jl")
            include("unit-tests/Optimization-test.jl")
        end
        println("$(unit2_3) test finished in $t2_3 seconds.\n")

    end

    @testset "Submodule Functions" begin
        unit3_1 = "Submodule: Molden"
        println("Testing $(unit3_1)...")
        t3_1 = @elapsed @testset "$(unit3_1)" begin
            include("unit-tests/SubModule/Molden-test.jl")
        end
        println("$(unit3_1) test finished in $t3_1 seconds.\n")
    end

end


text:
ToDataVec(M::AbstractVector) = floatify(M)
ToDataVec(M::AbstractMatrix) = reduce(vcat, collect(eachrow(ToArray(M))))
ToDataVec(M::DataFrame) = M |> ToArray |> ToDataVec

ToArray(df::AbstractVector{<:Real}) = df |> floatify
ToArray(df::AbstractMatrix) = size(df,2) == 1 ? floatify(df[:,1]) : floatify(df)
ToArray(df::DataFrame) = size(df,2) == 1 ? convert(Vector{suff(df)}, floatify(df[:,1])) : convert(Matrix{suff(df)}, floatify(df))
function ToArray(df::AbstractVector{<:Union{Missing, AbstractFloat}})
    !all(x-> !ismissing(x), df) && throw("Input contains missing values.")
    convert(Vector{suff(df[1])},floatify(df))
end

function ReadIn(df::DataFrame, xdims::Int=1, ydims::Int=Int((size(df,2)-1)/2); xerrs::Bool=false, stripedXs::Bool=true, stripedYs::Bool=true)
    if xerrs
        (size(df,2) != 2xdims + 2ydims) && throw("Inconsistent no. of columns on DataFrame: got $(size(df,2))")
        Xcols = stripedXs ? (1:2:2xdims) : (1:xdims)
        Xerrs = stripedXs ? (2:2:2xdims) : (xdims+1:2xdims)
        Ycols = stripedYs ? ((2xdims+1):2:(2xdims + 2ydims)) : ((2xdims+1):(2xdims + ydims))
        Yerrs = stripedYs ? ((2xdims+2):2:(2xdims + 2ydims)) : ((2xdims+ydims+1):(2xdims + 2ydims))
    else
        (size(df,2) != xdims + 2ydims) && throw("Inconsistent no. of columns on DataFrame: got $(size(df,2))")
        Xcols = 1:xdims;        Xerrs = Val(false)
        Ycols = stripedYs ? ((xdims+1):2:(xdims + 2ydims)) : ((xdims+1):(xdims + ydims))
        Yerrs = stripedYs ? ((xdims+2):2:(xdims + 2ydims)) : ((xdims+ydims+1):(xdims + 2ydims))
    end
    xnames = names(df[:,Xcols]);    ynames = names(df[:,Ycols])
    println("Order of DataSets: $xnames $ynames.")
    DSs = _ReadIn(df, Xcols, Xerrs, Ycols, Yerrs)
    InformNames(DSs, xnames, ynames)
end

function _ReadIn(df::DataFrame, xcols::AbstractVector{<:Int}, xerrs::AbstractVector{<:Int}, ycols::AbstractVector{<:Int}, yerrs::AbstractVector{<:Int})
    X = df[:, xcols];    Xerr = df[:, xerrs];   Y = df[:, ycols];    Yerr = df[:, yerrs]
    DSs = Array{AbstractDataSet}(undef, size(Y,2))
    for (i,Col) in enumerate(eachcol(Y))
        inds = DitchMissingRows(X, Col)
        DSs[i] = DataSetExact(ToDataVec(X[inds,:]), ToArray(Xerr[inds,:]), ToArray(Y[inds,:]), ToArray(Yerr[inds,:]), (sum(inds), size(X,2), 1))
    end;    DSs
end

function _ReadIn(df::DataFrame, xcols::AbstractVector{<:Int}, xerrs::Val{false}, ycols::AbstractVector{<:Int}, yerrs::AbstractVector{<:Int})
    X = df[:, xcols];    Y = df[:, ycols];    Yerr = df[:, yerrs]
    DSs = Array{AbstractDataSet}(undef, size(Y,2))
    for (i,Col) in enumerate(eachcol(Y))
        inds = DitchMissingRows(X, Col)
        DSs[i] = DataSet(ToDataVec(X[inds,:]), ToArray(Y[inds,i]), ToArray(Yerr[inds,i]), (sum(inds), size(X,2), 1))
    end;    DSs
end

DitchMissingRows(df1, df2) = hcat(df1, df2) |> DitchMissingRows
DitchMissingRows(df) = DitchMissingRows(DataFrame(df, :auto))
function DitchMissingRows(df::Union{DataFrame, AbstractArray{<:Union{Missing,AbstractFloat}}})
    inds = falses(size(df,1))
    for i in eachindex(inds)
        if !any(ismissing, df[i,:])
            inds[i] = true
        end
    end;    inds
end

function SplitDS(DS::DataSet)
    if typeof(ysigma(DS)) <: AbstractVector
        return [InformNames(DataSet(xdata(DS), ydata(DS)[i:ydim(DS):end], ysigma(DS)[i:ydim(DS):end], (Npoints(DS), xdim(DS), 1)), xnames(DS), ynames(DS)[i:ydim(DS):end]) for i in 1:ydim(DS)]
    else
        return [InformNames(DataSet(xdata(DS), ydata(DS)[i:ydim(DS):end], ysigma(DS)[i:ydim(DS):end,i:ydim(DS):end], (Npoints(DS), xdim(DS), 1)), xnames(DS), ynames(DS)[i:ydim(DS):end]) for i in 1:ydim(DS)]
    end
end
function SplitDS(DS::DataSetExact)
    if typeof(ysigma(DS)) <: AbstractVector
        return [InformNames(DataSetExact(xdata(DS), xsigma(DS), ydata(DS)[i:ydim(DS):end], ysigma(DS)[i:ydim(DS):end], (Npoints(DS), xdim(DS), 1)), xnames(DS), ynames(DS)[i:ydim(DS):end]) for i in 1:ydim(DS)]
    else

        return [InformNames(DataSetExact(xdata(DS), xsigma(DS), ydata(DS)[i:ydim(DS):end], ysigma(DS)[i:ydim(DS):end,i:ydim(DS):end], (Npoints(DS), xdim(DS), 1)), xnames(DS), ynames(DS)[i:ydim(DS):end]) for i in 1:ydim(DS)]
    end
end

# Add Namelist of y-components to be able to look up order of columns
"""
The `CompositeDataSet` type is a more elaborate (and typically less performant) container for storing data.
Essentially, it splits observed data which has multiple `y`-components into separate data containers (e.g. of type `DataSet`), each of which corresponds to one of the components of the `y`-data.
Crucially, each of the smaller data containers still shares the same "kind" of `x`-data, that is, the same `xdim`, units and so on, although they do **not** need to share the exact same particular `x`-data.

The main advantage of this approach is that it can be applied when there are `missing` `y`-components in some observations.
A typical use case for `CompositeDataSet`s are time series where multiple quantities are tracked but not every quantity is necessarily recorded at each time step.
Example:
```julia
using DataFrames
t = [1,2,3,4]
y₁ = [2.5, 6, missing, 9];      y₂ = [missing, 5, 3.1, 1.4]
σ₁ = 0.3*ones(4);               σ₂ = [missing, 0.2, 0.1, 0.5]
df = DataFrame([t y₁ σ₁ y₂ σ])

xdim = 1;   ydim = 2
CompositeDataSet(df, xdim, ydim; xerrs=false, stripedYs=true)
```
The boolean-valued keywords `stripedXs` and `stripedYs` can be used to indicate to the constructor whether the values and corresponding ``1\\sigma`` uncertainties are given in alternating order, or whether the initial block of `ydim` many columns are the values and the second `ydim` many columns are the corresponding uncertainties.
Also, `xerrs=true` can be used to indicate that the `x`-values also carry uncertainties.
Basically all functions which can be called on other data containers such as `DataSet` have been specialized to also work with `CompositeDataSet`s.
"""
struct CompositeDataSet <: AbstractDataSet
    DSs::AbstractVector{<:AbstractDataSet}
    InvCov::AbstractMatrix{<:Number}
    logdetInvCov::Real
    WoundX::AbstractVector
    SharedYdim::Val
    function CompositeDataSet(pDSs::AbstractVector{<:AbstractDataSet})
        !all(DS->xdim(DS)==xdim(pDSs[1]), pDSs) && throw("Inconsistent dimensionality of x-data between data containers.")
        DSs = reduce(vcat, map(SplitDS, pDSs))
        InvCov = mapreduce(yInvCov, BlockMatrix, DSs) |> HealthyCovariance
        CompositeDataSet(DSs, InvCov, logdet(InvCov), unique(mapreduce(WoundX, vcat, DSs)), Val(all(DS->ydim(DS)==ydim(DSs[1]), DSs)))
    end
    function CompositeDataSet(DSs::AbstractVector{<:AbstractDataSet}, InvCov::AbstractMatrix, logdetInvCov::Real, WoundX::AbstractVector, SharedYdim::Val)
        new(DSs, InvCov, logdetInvCov, WoundX, SharedYdim)
    end
end
CompositeDataSet(DS::AbstractDataSet) = CompositeDataSet([DS])
function CompositeDataSet(df::DataFrame, xdims::Int=1, ydims::Int=Int((size(df,2)-1)/2); xerrs::Bool=false, stripedXs::Bool=true, stripedYs::Bool=true)
    CompositeDataSet(ReadIn(floatify(df), xdims, ydims; xerrs=xerrs, stripedXs=stripedXs, stripedYs=stripedYs))
end

# For SciMLBase.remake
CompositeDataSet(;
DSs::AbstractVector{<:AbstractDataSet}=[DataSet([0.],[0.],[1.])],
InvCov::AbstractMatrix=Diagonal([1,2.]),
logdetInvCov::Real=-Inf,
WoundX::AbstractVector=[0.],
SharedYdim::Val=Val(true)) = CompositeDataSet(DSs, logdetInvCov, WoundX, SharedYdim)

Data(CDS::CompositeDataSet) = CDS.DSs
xdata(CDS::CompositeDataSet) = mapreduce(xdata, vcat, Data(CDS))
ydata(CDS::CompositeDataSet) = mapreduce(ydata, vcat, Data(CDS))

# BlockReduce(X::AbstractVector{<:AbstractVector{<:Number}}) = reduce(vcat, X)
# BlockReduce(X::AbstractVector{<:AbstractMatrix{<:Number}}) = reduce(BlockMatrix, X)
BlockReduce(X::AbstractVector{<:AbstractArray{<:Number}}) = reduce(BlockMatrix, [(typeof(x) <: AbstractVector ? Diagonal(x.^2) : x) for x in X])

ysigma(CDS::CompositeDataSet) = map(ysigma, Data(CDS)) |> BlockReduce
xsigma(CDS::CompositeDataSet) = map(xsigma, Data(CDS)) |> BlockReduce
yInvCov(CDS::CompositeDataSet) = mapreduce(yInvCov, BlockMatrix, Data(CDS))

Npoints(CDS::CompositeDataSet) = mapreduce(Npoints, +, Data(CDS))
ydim(CDS::CompositeDataSet) = mapreduce(ydim, +, Data(CDS))
xdim(CDS::CompositeDataSet) = xdim(Data(CDS)[1])

WoundX(CDS::CompositeDataSet) = CDS.WoundX
logdetInvCov(CDS::CompositeDataSet) = CDS.logdetInvCov

DataspaceDim(CDS::CompositeDataSet) = mapreduce(DS->Npoints(DS)*ydim(DS), +, Data(CDS))

xnames(CDS::CompositeDataSet) = xnames(Data(CDS)[1])
ynames(CDS::CompositeDataSet) = mapreduce(ynames, vcat, Data(CDS))

function InformNames(CDS::CompositeDataSet, xnames::AbstractVector{String}, ynames::AbstractVector{String})
    CompositeDataSet(InformNames(Data(CDS), xnames, ynames))
end
function InformNames(DSs::AbstractVector{<:AbstractDataSet}, xnames::AbstractVector{String}, ynames::AbstractVector{String})
    # Use InformNames for single DataSet recursively
    @assert length(ynames) == sum(ydim.(DSs)) && all(x->xdim(x)==length(xnames), DSs)
    Res = Vector{AbstractDataSet}(undef, length(DSs))
    j = 1   # Use this to infer how many elements have been popped from ynames
    for i in 1:length(DSs)
        Res[i] = InformNames(DSs[i], xnames, ynames[j:j-1+ydim(DSs[i])])
        j += ydim(DSs[i])
    end;    Res
end

function _CustomOrNot(CDS::CompositeDataSet, model::Function, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{false}; kwargs...)
    @assert CDS.SharedYdim isa Val{true} && ydim(Data(CDS)[1]) == 1
    # reduce(vcat, transpose) faster than Unpack?
    X = unique(woundX)
    _FillResVector(CDS, X, reduce(vcat, map(z->transpose(model(z, θ; kwargs...)), X)))
end

# Apparently reduce(vcat, map(z->transpose(G(z)), X))  just as fast as   transpose(reshape(reduce(vcat, map(z->transpose(G(z)), X)), ydim, :))
function _CustomOrNot(CDS::CompositeDataSet, model::Function, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{true}; kwargs...)
    @assert CDS.SharedYdim isa Val{true} && ydim(Data(CDS)[1]) == 1
    # reduce(vcat, transpose) faster than Unpack?
    X = unique(woundX)
    _FillResVector(CDS, X, (ydim(CDS) == 1 ? reshape(model(X, θ; kwargs...), :, 1) : transpose(reshape(model(X, θ; kwargs...), ydim(CDS), :))))
end

function _FillResVector(CDS::CompositeDataSet, X::AbstractVector, Mapped::AbstractMatrix{<:Number})
    Res = Vector{suff(Mapped)}(undef, DataspaceDim(CDS));      i = 1
    for SetInd in 1:length(Data(CDS))
        for xval in WoundX(Data(CDS)[SetInd])
            # Res[i] = view(Mapped, findfirst(isequal(xval),X), SetInd]
            Res[i] = Mapped[findfirst(isequal(xval),X), SetInd]
            i += 1
        end
    end;    return Res
end

function _CustomOrNotdM(CDS::CompositeDataSet, dmodel::Function, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{false}; kwargs...)
    @assert CDS.SharedYdim isa Val{true} && ydim(Data(CDS)[1]) == 1
    X = unique(woundX)
    _FillResMatrix(CDS, X, map(z->dmodel(z,θ; kwargs...), X))
end

function _CustomOrNotdM(CDS::CompositeDataSet, dmodel::Function, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{true}; kwargs...)
    @assert CDS.SharedYdim isa Val{true} && ydim(Data(CDS)[1]) == 1
    X = unique(woundX);    Mapped = dmodel(X, θ; kwargs...)
    [view(Mapped, (1 + (i-1)*ydim(CDS)):(i*ydim(CDS)) , :) for i in 1:length(X)]
    _FillResMatrix(CDS, X, map(z->dmodel(z,θ; kwargs...), X))
end

_FillResMatrix(CDS::CompositeDataSet, X::AbstractVector, Mapped::AbstractVector{<:AbstractMatrix{<:Number}}) = reduce(vcat, map(i->_getViewDmod(CDS,i,X,Mapped), 1:length(Data(CDS))))

@inline function _getViewDmod(CDS::CompositeDataSet, SetInd::Int, X::AbstractVector, Mapped::AbstractVector{<:AbstractMatrix{<:Number}})
    subXs = WoundX(Data(CDS)[SetInd])
    Res = Mapped[findfirst(isequal(subXs[1]), X)][SetInd,:] |> transpose
    for i in 2:length(subXs)
        Res = vcat(Res, transpose(Mapped[findfirst(isequal(subXs[i]), X)][SetInd,:]))
    end;    Res
end

RecipesBase.@recipe function f(CDS::CompositeDataSet)
    xdim(CDS) != 1 && throw("Not programmed for plotting xdim != 1 yet.")
    !all(x->ydim(x)==1, Data(CDS)) && throw("Not programmed for plotting ydim > 1 yet.")
    xguide -->  xnames(CDS)[1]
    yguide -->  "Observations"
    for (i,DS) in enumerate(Data(CDS))
        @series begin
            label --> "Data: " * ynames(DS)[1]
            markercolor --> [:red,:blue,:green,:orange,:grey][i]
            DS
        end
    end
end

function ResidualStandardError(CDS::CompositeDataSet, model::ModelOrFunction, MLE::AbstractVector{<:Number})
    ydiff = ydata(CDS) - EmbeddingMap(CDS, model, MLE);    Res = zeros(length(Data(CDS)));    startind = 1
    for i in eachindex(Res)
        ypred = view(ydiff, startind:startind+Npoints(Data(CDS)[i])*ydim(Data(CDS)[i])-1)
        startind += DataspaceDim(Data(CDS)[i])
        Res[i] = sqrt(sum(abs2, ypred) / (DataspaceDim(Data(CDS)[i]) - length(MLE)))
    end
    @assert (startind - 1) == DataspaceDim(CDS)
    ydim(CDS) == 1 ? Res[1] : Res
end


text:
abstract type KNNAbstract{T} <: SurpriseModel{T} end

struct KNNBasic{T} <: KNNAbstract{T}
    object::PyObject
    k::Int
    min_k::Int
    preference::Persa.Preference{T}
    users::Int
    items::Int
end

"""
    KNNBasic(dataset::Persa.Dataset; k = 40, min_k = 1)

A basic KNN algorithm.

# Arguments
- `k::Int = 40`: Maximum Number of neighbors used.
- `min_k::Int = 1`: Minimum number of neighbors used.
"""
function KNNBasic(dataset::Persa.Dataset; k = 40, min_k = 1)
    return KNNBasic(
        surprise.KNNBasic(k = k, min_k = min_k),
        k,
        min_k,
        dataset.preference,
        Persa.users(dataset),
        Persa.items(dataset)
    )
end

struct KNNBaseline{T} <: KNNAbstract{T}
    object::PyObject
    k::Int
    min_k::Int
    preference::Persa.Preference{T}
    users::Int
    items::Int
end

"""
    KNNBaseline(dataset::Persa.Dataset; k = 40, min_k = 1)
A basic KNN algorithm but using a baseline factor.
# Arguments
- `k::Int = 40`: Maximum Number of neighbors used.
- `min_k::Int = 1`: Minimum number of neighbors used.
"""
function KNNBaseline(dataset::Persa.Dataset; k = 40, min_k = 1)
    return KNNBaseline(
        surprise.KNNBaseline(k = k, min_k = min_k),
        k,
        min_k,
        dataset.preference,
        Persa.users(dataset),
        Persa.items(dataset)
    )
end

struct KNNWithMeans{T} <: KNNAbstract{T}
    object::PyObject
    k::Int
    min_k::Int
    preference::Persa.Preference{T}
    users::Int
    items::Int
end

"""
    KNNWithMeans(dataset::Persa.Dataset; k = 40, min_k = 1)
A basic KNN algorithm but using user or item mean.
# Arguments
- `k::Int = 40`: Maximum Number of neighbors used.
- `min_k::Int = 1`: Minimum number of neighbors used.
"""
function KNNWithMeans(dataset::Persa.Dataset; k::Int = 40, min_k::Int = 1)
    sim_options = Dict()
    sim_options["name"] = "cosine"
    return KNNWithMeans(
        surprise.KNNWithMeans(k = k, min_k = min_k, sim_options = sim_options),
        k,
        min_k,
        dataset.preference,
        Persa.users(dataset),
        Persa.items(dataset)
    )
end

struct SlopeOne{T} <: SurpriseModel{T}
    object::PyObject
    preference::Persa.Preference{T}
    users::Int
    items::Int
end

"""
    SlopeOne(dataset::Persa.Dataset)
SlopeOne algorithm.
"""
function SlopeOne(dataset::Persa.Dataset)
    return SlopeOne(
        surprise.SlopeOne(),
        dataset.preference,
        Persa.users(dataset),
        Persa.items(dataset)
    )
end

function Persa.predict(model::SurpriseModel, user::Int, item::Int)
    uid, vid = rawtoid(model.object, user, item)

    if isnan(uid) || isnan(vid)
        return NaN
    end

    return model.object.estimate(uid, vid)
end

function Persa.predict(model::KNNAbstract, user::Int, item::Int)
    uid, vid = rawtoid(model.object, user, item)

    if isnan(uid) || isnan(vid)
        return missing
    end

    return model.object.estimate(uid, vid)[1]
end


text:
module TestIR

import PerformanceTestTools
PerformanceTestTools.@include("__test_ir.jl")

end  # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule wrfuser_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("wrfuser")
JLLWrappers.@generate_main_file("wrfuser", UUID("66307a7f-b63f-5cd2-ae25-ee9b20ebbd65"))
end  # module wrfuser_jll


text:
"""
Tableau of Crank-Nicolson two-stage, 2nd order method

```julia
TableauCrankNicolson(::Type{T}=Float64) where {T}
```

The constructor takes one optional argument, that is the element type of the tableau.

Reference:

    J. Crank and P. Nicolson.
    A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type.
    Mathematical Proceedings of the Cambridge Philosophical Society, Volume 43, Issue 1, Pages 50-67, 1947.
    doi: 10.1017/S0305004100023197

"""
function TableauCrankNicolson(::Type{T}=Float64) where {T}
    a = @big [[ 0     0   ]
              [ 1/2   1/2 ]]
    b = @big  [ 1/2,  1/2 ]
    c = @big  [ 0,    1   ]
    o = 2

    Tableau{T}(:cranknicolson, o, a, b, c)
end

"""
Tableau of Kraaijevanger and Spijker's two-stage, 2nd order method

```julia
TableauKraaijevangerSpijker(::Type{T}=Float64) where {T}
```

The constructor takes one optional argument, that is the element type of the tableau.

Reference:

    J. F. B. M. Kraaijevanger and M. N. Spijker.
    Algebraic stability and error propagation in Runge-Kutta methods.
    Applied Numerical Mathematics, Volume 5, Issues 1-2, Pages 71-87, 1989.
    doi: 10.1016/0168-9274(89)90025-1

"""
function TableauKraaijevangerSpijker(::Type{T}=Float64) where {T}
    a = @big [[ 1/2   0   ]
              [-1/2   2   ]]
    b = @big  [-1/2,  3/2 ]
    c = @big  [ 1/2,  3/2 ]
    o = 2

    Tableau{T}(:kraaijevangerspijker, o, a, b, c)
end

"""
Tableau of Qin and Zhang's symplectic two-stage, 2nd order method

```julia
TableauQinZhang(::Type{T}=Float64) where {T}
```

The constructor takes one optional argument, that is the element type of the tableau.

Reference:

    M.-Z. Qin and M.-Q. Zhang.
    Symplectic Runge-Kutta algorithms for Hamilton systems.
    Journal of Computational Mathematics, Supplementary Issue, Pages 205-215, 1992.

"""
function TableauQinZhang(::Type{T}=Float64) where {T}
    a = @big [[ 1/4   0   ]
              [ 1/2   1/4 ]]
    b = @big  [ 1/2,  1/2 ]
    c = @big  [ 1/4,  3/4 ]
    o = 2

    Tableau{T}(:qinzhang, o, a, b, c)
end

"""
Tableau of Crouzeix's two-stage, 3rd order method

```julia
TableauCrouzeix(::Type{T}=Float64) where {T}
```

The constructor takes one optional argument, that is the element type of the tableau.

Reference:

    M.Crouzeix.
    Sur L'approximation des équations différentielles opérationelles linéaires par des méthodes de Runge-Kutta.
    Thesis. Université de Paris, 1975.

"""
function TableauCrouzeix(::Type{T}=Float64) where {T}
    fac = @big 1/2 / √3
    a = @big [[ 1/2+fac  0       ]
              [  -2*fac  1/2+fac ]]
    b = @big  [ 1/2,     1/2     ]
    c = @big  [ 1/2+fac, 1/2-fac ]
    o = 3

    Tableau{T}(:crouzeix, o, a, b, c)
end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

# convert binary into HDF5 data
using HDF5

base_dir = "./"

if length(ARGS) > 0
  base_dir = ARGS[1]
end

datasets = Dict("train" => ["$base_dir/train-labels-idx1-ubyte","$base_dir/train-images-idx3-ubyte"],
                "test" => ["$base_dir/t10k-labels-idx1-ubyte","$base_dir/t10k-images-idx3-ubyte"])

for key in keys(datasets)
  label_fn, data_fn = datasets[key]
  label_f = open(label_fn)
  data_f  = open(data_fn)

  label_header = read(label_f, Int32, 2)
  @assert ntoh(label_header[1]) == 2049
  n_label = round(Int, ntoh(label_header[2]))
  data_header = read(data_f, Int32, 4)
  @assert ntoh(data_header[1]) == 2051
  n_data = round(Int, ntoh(data_header[2]))
  @assert n_label == n_data
  h = round(Int, ntoh(data_header[3]))
  w = round(Int, ntoh(data_header[4]))

  println("Exporting $n_data digits of size $h x $w")

  h5open("$base_dir/$key.hdf5", "w") do h5
    dset_data = d_create(h5, "data", datatype(Float32), dataspace(w, h, 1, n_data))
    dset_label = d_create(h5, "label", datatype(Float32), dataspace(1, n_data))

    img = readbytes(data_f, n_data * h*w)
    img = convert(Array{Float32},img) / 256 # scale into [0,1)
    class = readbytes(label_f, n_data)
    class = convert(Array{Float32},class)

    idx = 1:n_data
    println("  $idx...")

    idx = collect(idx)
    rp = randperm(length(idx))
    for j = 1:length(idx)
      r_idx = rp[j]
      dset_data[:,:,1,idx[j]] = img[(r_idx-1)*h*w+1:r_idx*h*w]
      dset_label[1,idx[j]] = class[r_idx]
    end
  end

  close(label_f)
  close(data_f)
end


text:
# Generic concepts and functions

# a premetric is a function d that satisfies:
#
#   d(x, y) >= 0
#   d(x, x) = 0
#
abstract type PreMetric end

# a semimetric is a function d that satisfies:
#
#   d(x, y) >= 0
#   d(x, x) = 0
#   d(x, y) = d(y, x)
#
abstract type SemiMetric <: PreMetric end

# a metric is a semimetric that satisfies triangle inequality:
#
#   d(x, y) + d(y, z) >= d(x, z)
#
abstract type Metric <: SemiMetric end

# Generic functions

result_type(::PreMetric, ::AbstractArray, ::AbstractArray) = Float64

# Generic column-wise evaluation

function colwise!(r::AbstractArray, metric::PreMetric, a::AbstractVector, b::AbstractMatrix)
    n = size(b, 2)
    length(r) == n || throw(DimensionMismatch("Incorrect size of r."))
    @inbounds for j = 1:n
        r[j] = evaluate(metric, a, view(b, :, j))
    end
    r
end

function colwise!(r::AbstractArray, metric::PreMetric, a::AbstractMatrix, b::AbstractVector)
    n = size(a, 2)
    length(r) == n || throw(DimensionMismatch("Incorrect size of r."))
    @inbounds for j = 1:n
        r[j] = evaluate(metric, view(a, :, j), b)
    end
    r
end

function colwise!(r::AbstractArray, metric::PreMetric, a::AbstractMatrix, b::AbstractMatrix)
    n = get_common_ncols(a, b)
    length(r) == n || throw(DimensionMismatch("Incorrect size of r."))
    @inbounds for j = 1:n
        r[j] = evaluate(metric, view(a, :, j), view(b, :, j))
    end
    r
end

function colwise!(r::AbstractArray, metric::SemiMetric, a::AbstractMatrix, b::AbstractVector)
    colwise!(r, metric, b, a)
end

function colwise(metric::PreMetric, a::AbstractMatrix, b::AbstractMatrix)
    n = get_common_ncols(a, b)
    r = Vector{result_type(metric, a, b)}(undef, n)
    colwise!(r, metric, a, b)
end

function colwise(metric::PreMetric, a::AbstractVector, b::AbstractMatrix)
    n = size(b, 2)
    r = Vector{result_type(metric, a, b)}(undef, n)
    colwise!(r, metric, a, b)
end

function colwise(metric::PreMetric, a::AbstractMatrix, b::AbstractVector)
    n = size(a, 2)
    r = Vector{result_type(metric, a, b)}(undef, n)
    colwise!(r, metric, a, b)
end

# Generic pairwise evaluation

function pairwise!(r::AbstractMatrix, metric::PreMetric, a::AbstractMatrix, b::AbstractMatrix)
    na = size(a, 2)
    nb = size(b, 2)
    size(r) == (na, nb) || throw(DimensionMismatch("Incorrect size of r."))
    @inbounds for j = 1:size(b, 2)
        bj = view(b, :, j)
        for i = 1:size(a, 2)
            r[i, j] = evaluate(metric, view(a, :, i), bj)
        end
    end
    r
end

function pairwise!(r::AbstractMatrix, metric::PreMetric, a::AbstractMatrix)
    pairwise!(r, metric, a, a)
end

function pairwise!(r::AbstractMatrix, metric::SemiMetric, a::AbstractMatrix)
    n = size(a, 2)
    size(r) == (n, n) || throw(DimensionMismatch("Incorrect size of r."))
    @inbounds for j = 1:n
        aj = view(a, :, j)
        for i = (j + 1):n
            r[i, j] = evaluate(metric, view(a, :, i), aj)
        end
        r[j, j] = 0
        for i = 1:(j - 1)
            r[i, j] = r[j, i]   # leveraging the symmetry of SemiMetric
        end
    end
    r
end

function pairwise(metric::PreMetric, a::AbstractMatrix, b::AbstractMatrix)
    m = size(a, 2)
    n = size(b, 2)
    r = Matrix{result_type(metric, a, b)}(undef, m, n)
    pairwise!(r, metric, a, b)
end

function pairwise(metric::PreMetric, a::AbstractMatrix)
    n = size(a, 2)
    r = Matrix{result_type(metric, a, a)}(undef, n, n)
    pairwise!(r, metric, a)
end


text:
"""
    abstract type Stage end

Base type that represents a step of iterated staging.

Iterated staging allows the program to repeatedly self-recompile its parts.
The first iterations are

`Stage` is the root of a layered type hierarchy:

- Direct subtypes of it (ContextStage, EvalStage) represent staging
techniques available in Julia.
- Downstream subtypes represent means of staging:
Initialization, Extension, Optimization, Configuration.
"""
abstract type Stage end

"""
    abstract type ContextStage <: Stage end

Stage technique that generates a new stage context type, which can be used
for dispatching or in `@generated` functions.
Assembled types may be regenerated during a `ContextStage`, depending on TODO.

A context stage does not create a new world but run in
the same world than the previous stage.
"""
abstract type ContextStage <: Stage end

"""
    abstract type EvalStage <: Stage end

Stage technique that runs in a new world and generates a new stage context type.
"""
abstract type EvalStage <: Stage end

"""
    abstract type Optimization <: ContextStage end

Stage that does not change the functional behavior of the program,
only its performance characteristics.
"""
abstract type Optimization <: ContextStage end

"""
    abstract type Configuration <: ContextStage end

Stage that potentially changes the behavior of the program
without evaluating previously unknown code.
"""
abstract type Configuration <: ContextStage end

"""
    abstract type Initialization <: EvalStage end

A classical Stage that runs before the normal operation of the program.

Multiple initialization stages may run, but only before the first
non-Initialization stage.
"""
abstract type Initialization <: EvalStage end

"""
    abstract type Extension <: EvalStage end

Stage that evaluates previously unknown code, e.g. loads
new plugins.
"""
abstract type Extension <: EvalStage end


text:
using Documenter
using OnlinePCA

makedocs()

makedocs(
    format = :html,
    sitename = "OnlinePCA.jl",
    modules = [OnlinePCA],
    pages = [
        "Home" => "index.md",
        "Julia API" => "juliaapi.md",
        "Command line Tool" => "commandline.md"
    ])

deploydocs(
    repo = "github.com/rikenbit/OnlinePCA.jl.git",
    julia = "1.0",
    target = "build",
    deps = nothing,
    make = nothing)


text:
"""
    struct BoolVar <: AbstractVar

A "simple" boolean variable.
The constraints that affect this variable are stored in the `onDomainChange` array.
"""
struct BoolVar <: AbstractBoolVar
    onDomainChange      ::Array{Constraint}
    domain              ::SeaPearl.BoolDomain
    id                  ::String
end

"""
    function BoolVar(id::String, trailer::Trailer)

Create a `BoolVar` with a domain being equal to [false, true] with the `id` string identifier
and that will be backtracked by `trailer`.
"""
function BoolVar(id::String, trailer::Trailer)
    dom = BoolDomain(trailer)

    return BoolVar(Constraint[], dom, id)
end

function Base.show(io::IO, var::BoolVar)
    print(io, var.id, " = ", var.domain)
end

function Base.show(io::IO, ::MIME"text/plain", var::BoolVar)
    print(io, typeof(var), ": ", var.id, " = ", var.domain)
end

"""
    isbound(x::AbstractBoolVar)

Check whether x has an assigned value.
"""
isbound(x::AbstractBoolVar) = length(x.domain) == 1

"""
    assign!(x::BoolVar, value::Bool)

Remove everything from the domain of `x` but `value`.
"""
assign!(x::BoolVar, value::Bool) = assign!(x.domain, value)

"""
    assignedValue(x::BoolVar)

Return the assigned value of `x`. Throw an error if `x` is not bound.
"""
function assignedValue(x::BoolVar)
    @assert isbound(x)

    return convert(Bool, minimum(x.domain.inner))
end

"""
    rootVariable(x::BoolVar)

Return the "true" variable behind `x`. For a `BoolVar`, it simply returns `x`.
"""
rootVariable(x::BoolVar) = x


text:
using RunMetrics, Test

cd(@__DIR__)

sample_data = joinpath("..", "sample_data", "sample.tcx")

rs = parse_run(sample_data)
n = length(rs)

@test n == 946

lap_df = parse_laps(sample_data)
@test lap_df.seconds[end] == 21.0
@test lap_df.meters[end] == 66.24