text:
module Jchemo  # Start-Module

using LinearAlgebra, SparseArrays
using Statistics, Random
using Clustering
using DataInterpolations
using Distributions
using DataFrames
using Distances
using GLMakie
using HypothesisTests
using ImageFiltering
using Interpolations
using LIBSVM
using NearestNeighbors
using StatsBase    # sample
using XGBoost

include("utility.jl")
include("preprocessing.jl")
include("rmgap.jl")
include("aov1.jl")
include("plots.jl")
include("matW.jl")

include("fda.jl")
include("pcasvd.jl") ; include("pcaeigen.jl")
include("kpca.jl")
include("rp.jl")

# Regression
include("mlr.jl")
include("rr.jl")
include("plskern.jl") ; include("plsrosa.jl")
include("plsnipals.jl") ; include("plssimp.jl")
include("plsr_agg.jl")
include("cglsr.jl")
include("pcr.jl")
include("krr.jl") ; include("kplsr.jl") ; include("dkplsr.jl")
include("aicplsr.jl")
include("wshenk.jl")

# SVM
include("svmr.jl")

# Bagging
include("baggr.jl") ; include("baggr_util.jl")

# Trees
include("treer_xgb.jl")
include("treeda_xgb.jl")

include("xfit.jl")
include("scordis.jl")

# Discrimination
include("dmnorm.jl")
include("rrda.jl")
include("lda.jl") ; include("qda.jl")
include("mlrda.jl")
include("plsrda.jl")
include("plslda.jl")
include("plsrda_agg.jl")

# SVM
include("svmda.jl")

# Local regression
include("locw.jl")
include("knnr.jl")
include("lwplsr.jl")
include("lwplsr_agg.jl")
include("cplsr_agg.jl")  # Use structure PlsrDa

# Multi-block regresssion
include("angles.jl")
include("blockscal.jl")
include("mbplsr.jl")
include("mbplsr_avg.jl")
include("mbplsr_mid.jl")
include("mbplsr_rosa.jl")
include("mbplsr_so.jl")

# Local discrimination
include("lwplsrda.jl") ; include("lwplslda.jl") ; include("lwplsqda.jl")
include("lwplsrda_agg.jl") ; include("lwplslda_agg.jl") ; include("lwplsqda_agg.jl")
include("knnda.jl")

# Variable selection/importance (direct methods)
include("covsel.jl")
include("iplsr.jl")

# Validation
include("mpars.jl")
include("scores.jl")
include("gridscore.jl")
include("segm.jl") ; include("gridcv.jl")

# Sampling
include("sampling.jl")

include("distances.jl")
include("getknn.jl")
include("wdist.jl")
include("kernels.jl")

export
    # Utilities
    aggstat,
    aov1,
    checkdupl, checkmiss,
    center, center!,
    colmeans, colvars, colstds,
    dummy,
    ensure_df, ensure_mat,
    findmax_cla,
    list,
    mweights,
    matB, matcov, matW,
    pnames, psize,
    recodcat2num, recodnum2cla,
    replacebylev,
    rmcol, rmrow,
    scale, scale!,
    sourcedir,
    summ,
    tab, tabnum,
    vcol, vrow,
   # Pre-processing
    snv, snv!, detrend, detrend!, fdif, fdif!,
    mavg, mavg!, mavg_runmean, mavg_runmean!,
    savgk, savgol, savgol!,
    interpl, interpl_mon,
    linear_int, quadratic_int, quadratic_spline, cubic_spline,
    rmgap, rmgap!,
    eposvd,
    # Pca
    pcasvd, pcasvd!, pcaeigen, pcaeigen!, pcaeigenk, pcaeigenk!,
    kpca,
    scordis, odis,
    rpmat_gauss, rpmat_li, rp,
    # Regression
    mlr, mlr!, mlrchol, mlrchol!, mlrpinv, mlrpinv!, mlrpinv_n, mlrpinv_n!,
    mlrvec!, mlrvec,
    plskern, plskern!, plsnipals, plsnipals!,
    plsrosa, plsrosa!, plssimp, plssimp!,
    cglsr, cglsr!,
    pcr,
    rr, rr!, rrchol, rrchol!,
    krr, kplsr, kplsr!, dkplsr, dkplsr!,
    plsr_agg, plsr_agg!,
    dfplsr_cg, aicplsr,
    svmr,
    treer_xgb, rfr_xgb, xgboostr, vimp_xgb,
    baggr, baggr_vi, baggr_oob,
    # Multi-block
    mblocks, blockscal, rv,
    mbplsr,
    mbplsr_avg,
    mbplsr_mid,
    mbplsr_mid_seq, mbplsr_mid_seq!,
    mbplsr_rosa, mbplsr_rosa!,
    mbplsr_so,
    # Variable selection/importance (direct methods)
    covsel,
    iplsr,
    #
    xfit, xfit!, xresid, xresid!,
    # Local regression
    locw, locwlv,
    knnr, lwplsr, lwplsr_agg,
    cplsr_agg,
    # Discrimination
    dmnorm, dmnorm!,
    fda, fda!, fdasvd, fdasvd!,
    mlrda,
    rrda, krrda,
    lda, qda,
    plsrda, kplsrda,
    plslda, plsqda,
    plsrda_agg, plslda_agg, plsqda_agg,
    svmda,
    treeda_xgb, rfda_xgb, xgboostda,
    # Local Discrimination
    lwplsrda, lwplslda, lwplsqda,
    lwplsrda_agg, lwplslda_agg, lwplsqda_agg,
    knnda,
    #
    transform, coef, predict,
    # Validation
    residreg, residcla,
    ssr, msep, rmsep, bias, sep, cor2, r2, rpd, rpdr, mse, err,
    mpars,
    gridscore, gridscorelv, gridscorelb,
    segmts, segmkf,
    gridcv, gridcvlv, gridcvlb,
    gridcv_mb, gridcvlv_mb,
    # Sampling
    sampks, sampdp, sampsys, sampclas,
    # Distances
    getknn, wdist, wdist!,
    euclsq, mahsq, mahsqchol,
    krbf, kpol,
    # Graphics
    plotsp
    # Not exported since surchage:
    # - summary => Base.summary

end # End-Module


text:
abstract type WilsonFermion_2D{NC} <: AbstractFermionfields_2D{NC}
end

function Wx!(xout::T,U::Array{G,1},x::T,A)  where  {T <: WilsonFermion_2D,G <: AbstractGaugefields}
    Wx!(xout,U,x,A,2)
    return
end

function Dx!(xout::T1,U::Array{G,1},x::T2,A) where  {T1 <: WilsonFermion_2D,T2 <: WilsonFermion_2D,G <: AbstractGaugefields}
    Dx!(xout,U,x,A,2)
    return
end

function Ddagx!(xout::T1,U::Array{G,1},x::T2,A) where  {T1 <: WilsonFermion_2D,T2 <: WilsonFermion_2D,G <: AbstractGaugefields}
    Ddagx!(xout,U,x,A,2)
    return
end

function Tx!(xout::T,U::Array{G,1},x::T,A)  where  {T <: WilsonFermion_2D,G <: AbstractGaugefields} # Tx, W = (1 - T)x
    Tx!(xout,U,x,A,2)
    return
end

function Wdagx!(xout::T,U::Array{G,1},
    x::T,A) where  {T <: WilsonFermion_2D,G <: AbstractGaugefields}
    #,temps::Array{T,1},boundarycondition) where  {T <: WilsonFermion_2D,G <: AbstractGaugefields}
    Wdagx!(xout,U,x,A,2)
    return
end

ext: jl


text:
using LightGraphs
import LightGraphs: induced_subgraph
using Base.Test

g = smallgraph(:bull)
n = 3
h = g[1:n]
@test nv(h) == n
@test ne(h) == 3

h = g[[1,2,4]]
@test nv(h) == n
@test ne(h) == 2

h = g[[1,5]]
@test nv(h) == 2
@test ne(h) == 0
@test typeof(h) == typeof(g)

g = DiGraph(100,200)
h = g[5:26]
@test nv(h) == 22
@test typeof(h) == typeof(g)

@test_throws ErrorException g[[1,1]]


text:
using PyPlot
using JLD
using DelimitedFiles
include("get_srb.jl")
function plot_smooth_indicator()
	fig = figure()
	ax = fig.add_subplot(111, projection="3d")
	eps = 1.e-8
	n_xpts, n_ypts = 100,100
	x1 = LinRange(0., 2*pi-eps, n_xpts)
	x2 = LinRange(0., 2*pi-eps, n_ypts)
	h1, h2 = 1.0, 1.0
	x1_gr = repeat(x1, 1, n_ypts)
	x2_gr = repeat(x2', n_xpts, 1)
	n_ctrs = 10
	a = 2*pi*rand(n_ctrs)
	b = 2*pi*rand(n_ctrs)
	x3_gr = zeros(n_xpts, n_ypts, n_ctrs)
	for (i, a1) = enumerate(a)
		x3_gr[:,:,i] = smooth_indicator.(x1_gr,
						x2_gr, a1, b[i], h1, h2)
		ax.plot_surface(x1_gr, x2_gr, x3_gr[:,:,i],
				label="\$ a = $a1 \$",cmap="Blues",
				vmin=0.,vmax=1.0)
	end
	ax.xaxis.set_tick_params(labelsize=28)
	ax.yaxis.set_tick_params(labelsize=28)
	ax.zaxis.set_tick_params(labelsize=28)
	ax.set_xlabel("x",fontsize=28)
	ax.set_ylabel("y",fontsize=28)
	ax.set_zlabel("z",fontsize=28)
	ax.set_title("\$ a = $(a[1]) \$",fontsize=28)
	#lgnd = fig.legend(fontsize=28)
end
function plot_density()
	#assumes get_dist() has been run
	fp = "../data/SRB_dist/"
	files = readdlm(string(fp, "files"), '\n')
	s = Array{String,1}(undef,4)
	for f in files
		s_str = split(string(f), ", ")
		s[1] = split(s_str[1],r"[a-z]|_|\[")[end]
		s[2],s[3] = s_str[2],s_str[3]
		s[4] = split(s_str[4],r"\]|_")[1]
		for (pind, pv) in enumerate(s)
			s[pind] = (length(pv) > 3) ?
					  pv[1:5] : pv[1:3]
		end
		filename = string(fp, string(f))
		X = load(filename)
		rho = X["rho"]
		n_x, n_y = size(rho)
		x = LinRange(0.,2*pi,n_x)
		y = LinRange(0.,2*pi,n_y)
		x_g = repeat(x, 1, n_y)
		y_g = repeat(y', n_x, 1)
		fig, ax = subplots(1,1)
		cplot = ax.contourf(x_g, y_g, rho,
					   cmap="Blues")
		cbar = fig.colorbar(cplot,ax=ax)
		ax.xaxis.set_tick_params(labelsize=28)
		ax.yaxis.set_tick_params(labelsize=28)
		cbar.ax.tick_params(labelsize=28)
		ax.set_title("s = $s ",fontsize=28)
	end
end
function plot_SRB_verify()
	#assumes get_dist() has been run
	fp = "../data/SRB_dist/verify/"
	files = readdlm(string(fp, "files"), '\n')
	for (i, f) = enumerate(files)
		filename = string(fp, string(f))
		X = load(filename)
		rho = X["rho"]
		n_x, n_y = size(rho)
		x = LinRange(0.,2*pi,n_x)
		y = LinRange(0.,2*pi,n_y)
		x_g = repeat(x, 1, n_y)
		y_g = repeat(y', n_x, 1)
		fig, ax = subplots(1,1)
		cplot = ax.contourf(x_g, y_g, rho,
					   cmap="Blues")
		cbar = fig.colorbar(cplot,ax=ax)
		ax.xaxis.set_tick_params(labelsize=28)
		ax.yaxis.set_tick_params(labelsize=28)
		cbar.ax.tick_params(labelsize=28)
		ax.set_title(string("bins =  ", n_x*n_y),fontsize=28)
	end
end

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
# Constructs a state-space representation of an ARIMA model given
# vectors of AR and MA coefficients and a number of differences.
# This implementation is based on Section 8.3 in Brockwell and Davis 2002,
# "Introduction to Time Series and Forecasting," 2nd ed.
function arima_statespace{T, I<:Integer}(ar::Vector{T}, d::I, ma::Vector{T},
		sigma::T)
	p = length(ar)
	q = length(ma)
	r = max(p, q + 1)

	F = [zeros(r - 1) diagm(ones(r - 1)); [zeros(max(0, r - p)), reverse(ar)]']
	V = diagm([zeros(r - 1), sigma])
	G = [zeros(max(0, r - q - 1)), reverse(ma), 1]'
	W = zeros(1, 1)
	# differencing
	if d > 0
		A = [zeros(d - 1), 1]
		B = [zeros(d - 1) diagm(ones(d - 1));
			 T[(-1.0)^(d + 1 - i) * binomial(d, d-i) for i in 0:(d - 1)]']
		# redefine matrices to incorporate differencing into state
		F = [F zeros(r, d); A*G B]
		V = [V zeros(r, d); zeros(d, r + d)]
		G = [G B[end, :]]
		x0 = zeros(r + d)
		P0 = diagm(ones(r + d) * 1e7)
	else
		x0 = zeros()
		P0 = diagm(ones(r) * 1e7)
	end
	StateSpaceModel(F, V, G, W, x0, P0)
end

function arima{T, I<:Integer}(y::Vector{T}, p::I, d::I, q::I)
	build(par::Vector{T}) = arima_statespace(par[1:p], par[p+1], par[p+2:end])
	par0 = zeros(p + d + q)
	fit(y, build, par0)
end


text:
"""
`numaxes(dev)`

Return the number of degrees of freedom of the actuator.
"""
numaxes(dev::AbstractActuator) = error("Not implemented!")
axesnames(dev::AbstractActuator) = error("Not implemented")

"""
`move(dev)`

Move a Robot.
"""
move(dev::AbstractRobot) = error("Not implemented!")

"""
`moveto(move, x)`

Move to an arbitrary point. The point is specified by vector `x`.

"""
function moveto end

"""
`moveX(dev)`

Move a Robot along x axis.
"""
moveX(dev::AbstractCartesianRobot) = error("Not implemented!")

"""
`moveY(dev)`

Move a Robot along y axis.
"""
moveY(dev::AbstractCartesianRobot) = error("Not implemented!")

"""
`moveZ(dev)`

Move a Robot along z axis.
"""
moveZ(dev::AbstractCartesianRobot) = error("Not implemented!")

"""
`rmove(dev)`

Move a Robot relative to the current position.
"""
rmove(dev::AbstractRobot) = error("Not implemented!")

"""
`rmoveX(dev)`

Move a Robot relative to the current position along x axis.
"""
rmoveX(dev::AbstractCartesianRobot) = error("Not implemented!")

"""
`rmoveY(dev)`

Move a Robot relative to the current position along y axis.
"""
rmoveY(dev::AbstractCartesianRobot) = error("Not implemented!")

"""
`rmoveZ(dev)`

Move a Robot relative to the current position along z axis.
"""
rmoveZ(dev::AbstractCartesianRobot) = error("Not implemented!")

"""
`absposition(dev)`

Position of the robot with respect to global reference frame.
"""
absposition(dev::AbstractRobot) = error("Not implemented")

"""
`abspositionX(dev)`

Position of the robot along x axis with respect to global reference frame.
"""
abspositionX(dev::AbstractCartesianRobot) = error("Not implemented")

"""
`abspositionY(dev)`

Position of the robot along y axis with respect to global reference frame.
"""
abspositionY(dev::AbstractCartesianRobot) = error("Not implemented")

"""
`abspositionZ(dev)`

Position of the robot along z axis with respect to global reference frame.
"""
abspositionZ(dev::AbstractCartesianRobot) = error("Not implemented")

import Base

"""
`position(dev)`

Position of the robot with respect to current frame of reference.
"""
Base.position(dev::AbstractRobot) = error("Not implemented")

"""
`positionX(dev)`

Position of the robot along x axis with respect to current frame of reference.
"""
positionX(dev::AbstractCartesianRobot) = error("Not implemented")

"""
`positionY(dev)`

Position of the robot along y axis with respect to current frame of reference.
"""
positionY(dev::AbstractCartesianRobot) = error("Not implemented")

"""
`positionZ(dev)`

Position of the robot along z axis with respect to current frame of reference.
"""
positionZ(dev::AbstractCartesianRobot) = error("Not implemented")

"""
`setreference(dev)`

Set current position as reference.
"""
setreference(dev::AbstractRobot) = error("Not implemented")

"""
`setreferenceX(dev)`

Set current position as reference along x axis.
"""
setreferenceX(dev::AbstractCartesianRobot) = error("Not implemented")

"""
`setreferenceY(dev)`

Set current position as reference along y axis.
"""
setreferenceY(dev::AbstractCartesianRobot) = error("Not implemented")

"""
`setreferenceZ(dev)`

Set current position as reference along z axis.
"""
setreferenceZ(dev::AbstractCartesianRobot) = error("Not implemented")

"""
`setabsreference(dev)`

Return to absolute frame of reference.
"""
setabsreference(dev::AbstractRobot) = error("Not implemented")

"""
`setabsreferenceX(dev)`

Return to absolute frame of reference along x axis.
"""
setabsreferenceX(dev::AbstractCartesianRobot) = error("Not implemented")

"""
`setabsreferenceY(dev)`

Return to absolute frame of reference along y axis.
"""
setabsreferenceY(dev::AbstractCartesianRobot) = error("Not implemented")

"""
`setabsreferenceZ(dev)`

Return to absolute frame of reference along z axis.
"""
setabsreferenceZ(dev::AbstractCartesianRobot) = error("Not implemented")

"""
`home(dev)`

Return home and set the absolute frame of reference.
"""
home(dev::AbstractRobot) = error("Not implemented!")

"""
`homeX(dev)`

Return home on x axis and set the absolute frame of reference on x axis.
"""
homeX(dev::AbstractCartesianRobot) = error("Not implemented!")

"""
`homeY(dev)`

Return home on y axis and set the absolute frame of reference on y axis.
"""
homeY(dev::AbstractCartesianRobot) = error("Not implemented!")

"""
`homeZ(dev)`

Return home on z axis and set the absolute frame of reference on z axis.
"""
homeZ(dev::AbstractCartesianRobot) = error("Not implemented!")

"""
`stopmotion(dev)`

Stop all motion of the robot.
"""
stopmotion(dev::AbstractRobot) = error("Not implemented!")

"""
`waituntildone(dev)`

Wait until all asynchronous motions of the robot is done.
"""
waituntildone(dev::AbstractRobot) = error("Not implemented!")


text:
function read_slist!(S::GphysData, fname::String, lennartz::Bool, memmap::Bool, strict::Bool, v::Integer)
  # file read
  io = memmap ? IOBuffer(Mmap.mmap(fname)) : open(fname, "r")
  hdr = readline(io)
  mark(io)
  nx = countlines(io)
  reset(io)
  X = Array{Float32,1}(undef, nx)
  i = 0
  while i < nx
    i += 1
    y = stream_float(io, 0x00)
    setindex!(X, y, i)
  end
  close(io)

  # header
  if lennartz
    id_sep = "."
    h = split(hdr)
    sta = replace(h[3], "\'" => "")
    cmp = last(split(fname, id_sep))
    id = *(id_sep, sta, id_sep, id_sep, cmp)
    ts = (Date(h[8]).instant.periods.value)*86400000000 +
          div(Time(h[9]).instant.value, 1000) -
          dtconst
    fs = 1000.0 / parse(Float64, h[5])
  else
    h = split(hdr, ',')
    id = join(split_id(split(h[1])[2], c="_"), ".")
    ts = 1000*DateTime(lstrip(h[4])).instant.periods.value - dtconst
    fs = parse(Float64, split(h[3])[1])
  end

  # Check for existing channel with same fs
  i = findid(S, id)
  if strict
    i = channel_match(S, i, fs)
  end
  if (i > 0)
    check_for_gap!(S, i, ts, nx, v)
    append!(S.x[i], X)
  else
    # New channel
    push!(S, SeisChannel(id = id, fs = fs, t = mk_t(nx, ts), x = X))
  end
  return nothing
end


text:
push!(LOAD_PATH, pwd())
include("startup.jl")

using HDF5
#using PyCall
#pygui(:tk) # :tk, :gtk3, :gtk, :qt5, :qt4, :qt, or :wx
using PyPlot
import .InterpolationFunctions
import  .MasslistFunctions
import .ResultFileFunctions

file = "ExampleFiles/TOFDATA/results/_result.hdf5"
plotHighTimeRes = true # Plot every datapoint or only file averages
plotFittedInsteadOfSummed = true # Use multi peak fitted data instead of raw
smoothing = 1 # Average n samples, 1 for raw
plotsymbol = ".-"
isobarToPlot = 0
#timedelay = Dates.Hour(0) # CLOUD12, ...
timedelay = Dates.Hour(1) # CLOUDX, CLOUD11

backgroundSubstractionMode = 0 # 0=no correction, 1=minimum found in all data,  2=mean of range given below
backgroundStart = DateTime(2016,10,02,17,14)
backgroundEnd = DateTime(2016,10,02,17,20)

include("manualMassLibrary.jl")

massesToPlot = [
# Examples for selecting what to plot:
MasslistFunctions.massFromComposition(C=10,H=16,O=2)
APINENE[1]
]

measResult = ResultFileFunctions.loadResults(file, massesToLoad=massesToPlot, useAveragesOnly=!plotHighTimeRes, raw=!plotFittedInsteadOfSummed, massMatchTolerance=0.01)
if isobarToPlot != 0
  isobarResult = ResultFileFunctions.loadResults(file, massesToLoad=[isobarToPlot+0.3], massMatchTolerance=0.5, useAveragesOnly=!plotHighTimeRes, raw=!plotFittedInsteadOfSummed)
  measResult=joinResultsMasses(measResult, isobarResult)
end

measResult.Times = measResult.Times .- timedelay

fig=figure()
ax = subplot(111)

if (backgroundSubstractionMode == 0)
  background=0
elseif (backgroundSubstractionMode == 1)
  background = minimum(InterpolationFunctions.averageSamples(measResult.Traces,smoothing),dims=1)
elseif backgroundSubstractionMode == 2
  import Statistics
  background = Statistics.mean(measResult.Traces[(measResult.Times.>backgroundStart) .& (measResult.Times.<backgroundEnd),:],dims=1)
end

bgCorrectedTraces = measResult.Traces .- background

semilogy(Dates.unix2datetime.(InterpolationFunctions.averageSamples(Dates.datetime2unix.(measResult.Times), smoothing)),InterpolationFunctions.averageSamples(bgCorrectedTraces,smoothing), plotsymbol) #linewidth=1)
#semilogy(Dates.unix2datetime(InterpolationFunctions.averageSamples(Dates.datetime2unix(measResult.Times), smoothing)),InterpolationFunctions.averageSamples(bgCorrectedTraces,smoothing), plotsymbol) #linewidth=1)

startTimeString = Dates.format(measResult.Times[1],"yyyy/mm/dd")
endTimeString = Dates.format(measResult.Times[end],"yyyy/mm/dd")
title("$startTimeString - $endTimeString")
xlabel("Time [UTC]")
ylabel("Signal [CPS]")

legStrings = []
for i = 1:length(measResult.MasslistMasses)
  push!(legStrings,"m/z $(round(measResult.MasslistMasses[i],digits=3)) - $(MasslistFunctions.sumFormulaStringFromCompositionArray(measResult.MasslistCompositions[:,i]))")
end

box = ax.get_position()
cols = 1

majorformatter = matplotlib.dates.DateFormatter("%m/%d %H:%M")
ax.xaxis.set_major_formatter(majorformatter)
legend(legStrings)
grid()

tight_layout()


text:
"""
An `Interval{L,R}(left, right)` where L,R are :open or :closed
is an interval set containg `x` such that
1. `left ≤ x ≤ right` if `L == R == :closed`
2. `left < x ≤ right` if `L == :open` and `R == :closed`
3. `left ≤ x < right` if `L == :closed` and `R == :open`, or
4. `left < x < right` if `L == R == :open`
"""
struct Interval{L,R,T}  <: TypedEndpointsInterval{L,R,T}
    left::T
    right::T

    Interval{L,R,T}(l, r) where {L,R,T} = ((a, b) = checked_conversion(T, l, r); new{L,R,T}(a, b))
end

"""
A `ClosedInterval(left, right)` is an interval set that includes both its upper and lower bounds. In
mathematical notation, the constructed range is `[left, right]`.
"""
const ClosedInterval{T} = Interval{:closed,:closed,T}

"""
An `TypedEndpointsInterval{:open,:open}(left, right)` is an interval set that includes both its upper and lower bounds. In
mathematical notation, the constructed range is `(left, right)`.
"""
const OpenInterval{T} = Interval{:open,:open,T}

Interval{L,R,T}(i::AbstractInterval) where {L,R,T} = Interval{L,R,T}(endpoints(i)...)
Interval{L,R}(left, right) where {L,R} = Interval{L,R,promote_type(typeof(left), typeof(right))}(left,right)
Interval{L,R}(left::T, right::T) where {L,R,T} = Interval{L,R,T}(left, right)
Interval(left, right) = ClosedInterval(left, right)

# Interval(::AbstractInterval) allows open/closed intervals to be changed
Interval{L,R}(i::AbstractInterval) where {L,R} = Interval{L,R}(endpoints(i)...)
Interval(i::AbstractInterval) = Interval{isleftclosed(i) ? (:closed) : (:open),
                                         isrightclosed(i) ? (:closed) : (:open)}(i)
Interval(i::TypedEndpointsInterval{L,R}) where {L,R} = Interval{L,R}(i)

endpoints(i::Interval) = (i.left, i.right)

for L in (:(:open),:(:closed)), R in (:(:open),:(:closed))
    @eval begin
        convert(::Type{Interval{$L,$R}}, i::Interval{$L,$R}) = i
        convert(::Type{Interval{$L,$R,T}}, i::Interval{$L,$R,T}) where T = i
    end
end
convert(::Type{Interval}, i::Interval) = i

function convert(::Type{II}, i::AbstractInterval) where II<:ClosedInterval
    isclosedset(i) ||  throw(InexactError(:convert,II,i))
    II(i)
end
function convert(::Type{II}, i::AbstractInterval) where II<:OpenInterval
    isopenset(i) ||  throw(InexactError(:convert,II,i))
    II(i)
end
function convert(::Type{II}, i::AbstractInterval) where II<:Interval{:open,:closed}
    (isleftopen(i) && isrightclosed(i)) ||  throw(InexactError(:convert,II,i))
    II(i)
end
function convert(::Type{II}, i::AbstractInterval) where II<:Interval{:closed,:open}
    (isleftclosed(i) && isrightopen(i)) ||  throw(InexactError(:convert,II,i))
    II(i)
end

convert(::Type{Interval}, i::AbstractInterval) = Interval(i)

convert(::Type{Domain{T}}, d::Interval{L,R}) where {L,R,T} = convert(Interval{L,R,T}, d)
convert(::Type{AbstractInterval{T}}, d::Interval{L,R}) where {L,R,T} = convert(Interval{L,R,T}, d)
convert(::Type{AbstractInterval{T}}, d::Interval{L,R,T}) where {L,R,T} = convert(Interval{L,R,T}, d)
convert(::Type{TypedEndpointsInterval{L,R,T}}, d::Interval{L,R}) where {L,R,T} = convert(Interval{L,R,T}, d)
convert(::Type{TypedEndpointsInterval{L,R,T}}, d::Interval{L,R,T}) where {L,R,T} = convert(Interval{L,R,T}, d)
convert(::Type{Domain}, d::Interval{L,R}) where {L,R} = d
convert(::Type{AbstractInterval}, d::Interval{L,R}) where {L,R} = d
convert(::Type{TypedEndpointsInterval{L,R}}, d::Interval{L,R}) where {L,R} = d

"""
    iv = l..r

Construct a ClosedInterval `iv` spanning the region from `l` to `r`.

(The symbol `..` is the same as in the package EllipsisNotation.jl.)
"""
..

(::Ellipsis)(x, y) = ClosedInterval(x, y)

"""
    iv = center±halfwidth

Construct a ClosedInterval `iv` spanning the region from
`center - halfwidth` to `center + halfwidth`.
"""
±(x, y) = ClosedInterval(x - y, x + y)
±(x::CartesianIndex, y) = (xy = y * one(x); map(ClosedInterval, (x - xy).I, (x + xy).I))

show(io::IO, I::ClosedInterval) = print(io, leftendpoint(I), "..", rightendpoint(I))
show(io::IO, I::OpenInterval) = print(io, leftendpoint(I), "..", rightendpoint(I), " (open)")
show(io::IO, I::Interval{:open,:closed}) = print(io, leftendpoint(I), "..", rightendpoint(I), " (open–closed)")
show(io::IO, I::Interval{:closed,:open}) = print(io, leftendpoint(I), "..", rightendpoint(I), " (closed–open)")

# The following are not typestable for mixed endpoint types
_left_intersect_type(::Type{Val{:open}}, ::Type{Val{L2}}, a1, a2) where L2 = a1 < a2 ? (a2,L2) : (a1,:open)
_left_intersect_type(::Type{Val{:closed}}, ::Type{Val{L2}}, a1, a2) where L2 = a1 ≤ a2 ? (a2,L2) : (a1,:closed)
_right_intersect_type(::Type{Val{:open}}, ::Type{Val{R2}}, b1, b2) where R2 = b1 > b2 ? (b2,R2) : (b1,:open)
_right_intersect_type(::Type{Val{:closed}}, ::Type{Val{R2}}, b1, b2) where R2 = b1 ≥ b2 ? (b2,R2) : (b1,:closed)

function intersect(d1::TypedEndpointsInterval{L1,R1,T}, d2::TypedEndpointsInterval{L2,R2,T}) where {L1,R1,L2,R2,T}
    a1, b1 = endpoints(d1); a2, b2 = endpoints(d2)
    a,L = _left_intersect_type(Val{L1}, Val{L2}, a1, a2)
    b,R = _right_intersect_type(Val{R1}, Val{R2}, b1, b2)
    Interval{L,R}(a,b)
end

function intersect(d1::TypedEndpointsInterval{L,R,T}, d2::TypedEndpointsInterval{L,R,T}) where {L,R,T}
    a1, b1 = endpoints(d1); a2, b2 = endpoints(d2)
    Interval{L,R}(max(a1,a2),min(b1,b2))
end

intersect(d1::TypedEndpointsInterval{L1,R1,T}, d2::TypedEndpointsInterval{L2,R2,V}) where {L1,R1,L2,R2,T,V} =
    intersect(promote(d1, d2)...)

intersect(d1::AbstractInterval, d2::AbstractInterval) = intersect(Interval(d1), Interval(d2))

union(d1::TypedEndpointsInterval{:closed,:closed},  d2::TypedEndpointsInterval{:closed,:closed}) =
    _leq_union(d1,d2)
union(d1::Interval{:open,:closed},                  d2::TypedEndpointsInterval{:closed,:closed}) =
    _leq_union(d1,d2)
union(d1::Interval{:closed,:open},                  d2::TypedEndpointsInterval{:closed,:closed}) =
    _leq_union(d1,d2)
union(d1::TypedEndpointsInterval{:closed,:closed},  d2::TypedEndpointsInterval{:open,:closed}) =
    _leq_union(d1,d2)
union(d1::TypedEndpointsInterval{:closed,:closed},  d2::TypedEndpointsInterval{:closed,:open}) =
    _leq_union(d1,d2)
union(d1::TypedEndpointsInterval{:open,:closed},    d2::TypedEndpointsInterval{:open,:closed}) =
    _leq_union(d1,d2)
union(d1::TypedEndpointsInterval{:closed,:open},    d2::TypedEndpointsInterval{:closed,:open}) =
    _leq_union(d1,d2)
union(d1::TypedEndpointsInterval{:closed,:closed},  d2::TypedEndpointsInterval{:open,:open}) =
    _leq_union(d1,d2)
union(d1::TypedEndpointsInterval{:open,:open},      d2::TypedEndpointsInterval{:closed,:closed}) =
    _leq_union(d1,d2)

function _leq_union(d1, d2)
    isempty(d1) && return d2
    isempty(d2) && return d1
    leftendpoint(d1) ≤ leftendpoint(d2) ≤ rightendpoint(d1)  || leftendpoint(d1) ≤ rightendpoint(d2) ≤ rightendpoint(d1) ||
        leftendpoint(d2) ≤ leftendpoint(d1) ≤ rightendpoint(d2) || leftendpoint(d2) ≤ rightendpoint(d1) ≤ rightendpoint(d2) ||
        throw(ArgumentError("Cannot construct union of disjoint sets."))
    _union(d1, d2)
end

function union(d1::TypedEndpointsInterval{:open,:open}, d2::TypedEndpointsInterval{:open,:open})
    isempty(d1) && return d2
    isempty(d2) && return d1
    leftendpoint(d1) ≤ leftendpoint(d2) < rightendpoint(d1)  || leftendpoint(d1) < rightendpoint(d2) ≤ rightendpoint(d1) ||
        leftendpoint(d2) ≤ leftendpoint(d1) < rightendpoint(d2)  || leftendpoint(d2) < rightendpoint(d1) ≤ rightendpoint(d2) ||
        throw(ArgumentError("Cannot construct union of disjoint sets."))
    _union(d1, d2)
end

function union(d1::TypedEndpointsInterval{:open,:open}, d2::TypedEndpointsInterval{:open,:closed})
    isempty(d1) && return d2
    isempty(d2) && return d1
    leftendpoint(d1) ≤ leftendpoint(d2) < rightendpoint(d1)  || leftendpoint(d1) ≤ rightendpoint(d2) ≤ rightendpoint(d1) ||
        leftendpoint(d2) ≤ leftendpoint(d1) ≤ rightendpoint(d2)  || leftendpoint(d2) < rightendpoint(d1) ≤ rightendpoint(d2) ||
        throw(ArgumentError("Cannot construct union of disjoint sets."))
    _union(d1, d2)
end

function union(d1::TypedEndpointsInterval{:open,:open}, d2::TypedEndpointsInterval{:closed,:open})
    isempty(d1) && return d2
    isempty(d2) && return d1
    leftendpoint(d1) ≤ leftendpoint(d2) ≤ rightendpoint(d1)  || leftendpoint(d1) < rightendpoint(d2) ≤ rightendpoint(d1) ||
        leftendpoint(d2) ≤ leftendpoint(d1) < rightendpoint(d2)  || leftendpoint(d2) ≤ rightendpoint(d1) ≤ rightendpoint(d2) ||
        throw(ArgumentError("Cannot construct union of disjoint sets."))
    _union(d1, d2)
end

union(d1::TypedEndpointsInterval, d2::TypedEndpointsInterval{:open,:open}) = union(d2, d1)

function union(d1::TypedEndpointsInterval{:closed,:open}, d2::TypedEndpointsInterval{:open,:closed})
    isempty(d1) && return d2
    isempty(d2) && return d1
    leftendpoint(d1) ≤ leftendpoint(d2) < rightendpoint(d1)  || leftendpoint(d1) ≤ rightendpoint(d2) ≤ rightendpoint(d1) ||
        leftendpoint(d2) ≤ leftendpoint(d1) ≤ rightendpoint(d2)  || leftendpoint(d2) < rightendpoint(d1) ≤ rightendpoint(d2) ||
        throw(ArgumentError("Cannot construct union of disjoint sets."))
    _union(d1, d2)
end

union(d1::TypedEndpointsInterval{:open,:closed}, d2::TypedEndpointsInterval{:closed,:open}) = union(d2, d1)

# these assume overlap
function _union(A::TypedEndpointsInterval{L,R}, B::TypedEndpointsInterval{L,R}) where {L,R}
    left = min(leftendpoint(A), leftendpoint(B))
    right = max(rightendpoint(A), rightendpoint(B))
    Interval{L,R}(left, right)
end

# this is not typestable
function _union(A::TypedEndpointsInterval{L1,R1}, B::TypedEndpointsInterval{L2,R2}) where {L1,R1,L2,R2}
    if leftendpoint(A) == leftendpoint(B)
        L = L1 == :closed ? :closed : L2
    elseif leftendpoint(A) < leftendpoint(B)
        L = L1
    else
        L = L2
    end
    if rightendpoint(A) == rightendpoint(B)
        R = R1 == :closed ? :closed : R2
    elseif rightendpoint(A) > rightendpoint(B)
        R = R1
    else
        R = R2
    end
    left = min(leftendpoint(A), leftendpoint(B))
    right = max(rightendpoint(A), rightendpoint(B))

    Interval{L,R}(left, right)
end

ClosedInterval{T}(i::AbstractUnitRange{I}) where {T,I<:Integer} = ClosedInterval{T}(minimum(i), maximum(i))
ClosedInterval(i::AbstractUnitRange{I}) where {I<:Integer} = ClosedInterval{I}(minimum(i), maximum(i))

Base.promote_rule(::Type{Interval{L,R,T1}}, ::Type{Interval{L,R,T2}}) where {L,R,T1,T2} = Interval{L,R,promote_type(T1, T2)}

# convert should only work if they represent the same thing.
@deprecate convert(::Type{R}, i::ClosedInterval{I}) where {R<:AbstractUnitRange,I<:Integer} R(i)
@deprecate length(i::ClosedInterval) IntervalSets.duration(i)


text:
"""
    TypeI(var_1::TypeJ)

Lorem ipsum dolor sit amet.
"""
struct TypeI

  var_1::TypeJ

end


text:
"""
    tab(df::DataFrame,vars::Symbol...; label_dict::Union{Nothing,Dict} = nothing, rmna::Bool = true, weights::AbstractVector = UnitWeights())

Produce an one-way or two-way frequency table from a DataFrame or any type of arrays.
Any weights defined under AbstractWeights (e.g., fweights, aweights, or pweights)
can be used to obtain weighted frequencies.

`tab` is mainly a wrapper for the excellent `FreqTables` package.
Use `skipmissing = true` to obtain frequencies that include `missing` values.
The returned table is a `NamedArray`. Frequencies are in an n-dimensional array `na.array`
where `na` is the returned NamedArray. `na.dimnames` contain row and column values.
"""
function tab(A::AbstractArray,B::AbstractArray; skipmissing::Bool=false,weights::AbstractWeights=fweights(ones(Int8,length(A))))

    # check lengths of arrays
    alen = length(A)
    blen = length(B)
    wlen = length(weights)

    if alen != blen
        error("Lengths of two vectors are not equal")
    end

    if alen != wlen
        error("Input vectors and `weights` vector are not of the same lengths")
    end

    # from FreqTables.jl/src/freqtable.jl
    d = Dict{Tuple{eltype(A),eltype(B)},Int64}()
    @inbounds @simd for i in 1:alen

        el = (A[i],B[i])

        index = Base.ht_keyindex(d, el)

        if index > 0
            @inbounds d.vals[index] += weights[i]
        else
            @inbounds d[el] = weights[i]
        end
    end

    # construct labels
    alev = Set{eltype(A)}()
    blev = Set{eltype(B)}()
    for (k1,k2) in keys(d)
        push!(alev,k1)
        push!(blev,k2)
    end

    alev = sort!(collect(alev))
    blev = sort!(collect(blev))
    if skipmissing
        alev = collect(skipmissing(alev))
        blev = collect(skipmissing(blev))
    end

    cnt = Array{Int64,2}(undef,length(alev),length(blev))
    @inbounds @simd for i=1:length(alev)
        @inbounds @simd for j=1:length(blev)
            el = (alev[i],blev[j])
            if haskey(d,el)
                cnt[i,j] = d[el]
            else
                cnt[i,j] = 0
            end
        end
    end

    return NamedArray(cnt,(alev,blev),("A","B"))
end
function tab(A::AbstractArray; skipmissing::Bool=true,weights::AbstractWeights=fweights(ones(Int8,length(A))))

    alen = length(A)
    wlen = length(weights)

    if alen != wlen
        error("The `weights` vector does not have the same length as the input vector `A`")
    end

    # from FreqTables.jl/src/freqtable.jl
    d = Dict{eltype(A),Int64}()
    @simd for i in 1:alen
        index = Base.ht_keyindex(d, A[i])

        if index > 0
            @inbounds d.vals[index] += weights[i]
        else
            @inbounds d[A[i]] = weights[i]
        end
    end

    # construct labels
    alev = sort(collect(keys(d)))

    cnt = Vector{Int64}(undef,length(alev))
    @inbounds @simd for i=1:length(alev)
        cnt[i] = d[alev[i]]
    end

    if skipmissing == true && alev[end] == missing
        return NamedArray(reshape(cnt[1:end-1],length(alev)-1,1),(alev[1:end-1],["Frequency"]),("A","Stat"))
    end

    return NamedArray(reshape(cnt,length(alev),1),(alev,["Frequency"]),("A","Stat"))
end
function tab(df::DataFrame,vars::Symbol...;  skipmissing::Bool=false, weights::AbstractWeights=fweights(ones(Int8,size(df,1))))
    nvec = length(vars)
    if nvec == 0
        error("At least one variable is required")
    elseif nvec > 2
        error("More than two variables are not supported")
    end

    # oneway table
    if nvec == 1
        na = tab(df[vars[1]], skipmissing = skipmissing, weights = weights)
        setdimnames!(na,(vars[1],"Stat"))
        return na
    end

    # two-way table
    na = tab(df[vars[1]], df[vars[2]], skipmissing = skipmissing, weights = weights)

    setdimnames!(na,(vars[1],vars[2]))
    return na
end

struct XsqResult
    chisq::Float64
    dof::Int
    p::Float64
end
function chi2test(a::Array{Float64,2})

    if ndims(t) != 2
        error("Only two dimensional arrays are supported")
    end

    rowsum = Array(sum(t,dims=2))
    colsum = Array(sum(t,dims=1))
    total = sum(t)

    ncol = length(colsum)
    nrow = length(rowsum)

    chisq = 0.
    for i = 1:nrow
        for j = 1:ncol
            expected = rowsum[i]*colsum[j]/total
            chisq += ((t[i,j] - expected)^2)/expected
        end
    end

    # degress of freedom
    df = (ncol-1)*(nrow-1)

    # return a tuple of chisq, df, p-value
    return XsqResult(chisq,df,Distributions.ccdf(Distributions.Chisq(df),chisq))
end

function rowpct(na::NamedArray)
    return 100. * na ./ sum(na.array,dims=2)
end

function colpct(na::NamedArray)
    return 100. * na ./ sum(na.array,dims=1)
end

function cellpct(na::NamedArray)
    return 100. * na ./ sum(na)
end

function chi2(t::NamedArray)
    if ndims(t) != 2
        error("Only two dimensional arrays are supported")
    end

    rowsum = sum(t,dims=2)
    colsum = sum(t,dims=1)
    total = sum(t)

    ncol = length(colsum)
    nrow = length(rowsum)

    chisq = 0.
    for i = 1:nrow
        for j = 1:ncol
            expected = rowsum[i]*colsum[j]/total
            chisq += ((t[i,j] - expected)^2)/expected
        end
    end

    # degress of freedom
    df = (ncol-1)*(nrow-1)

    # return a tuple of chisq, df, p-value
    return (chisq,df,Distributions.ccdf(Distributions.Chisq(df),chisq))
end

function rowpercent(tab::NamedArray)

    totrow = sum(tab,dims=2)

    array = 100 * tab.array ./ totrow

    return NamedArray(array,tuple(names(tab)...),tuple(dimnames(tab)...))

end

function colpercent(tab::NamedArray)

    totcol = sum(tab,dims=1)

    array = 100 * tab.array ./ tocol

    return NamedArray(array,tuple(names(tab)...),tuple(dimnames(tab)...))
end

function cellpercent(tab::NamedArray)

    total = sum(tab.array)

    array = 100 * tab.array ./ total

    return NamedArray(array,tuple(names(tab)...),tuple(dimnames(tab)...))
end

function tabprint(na::NamedArray; precision=2, labels::Union{Nothing,Label} = nothing, chisq=true, row=true, col=true, cell=false, all=false, pagewidth=78)

    ndim = ndims(na)
    if ndim == 2 && length(names(na,2)) == 1
        return tabprint1(na, precision = precision, labels = labels)
    elseif ndim == 3
        for (i,v) in enumerate(names(na,3))
            vlab = labels == nothing ? "" : string(" (",varlab(labels,dimnames(na,3)), " = ",vallab(labels,dimnames(na,3), v),")")
            print("\n",dimnames(na,3)," = ",v,vlab,"\n\n")
            tabprint(NamedArray(na.array[:,:,i],(names(na,1),names(na,2)),tuple(dimnames(na,1),dimnames(na,2))),
                precision=precision,labels=labels,chisq=chisq,row=row,col=col,cell=cell,all=all,pagewidth=pagewidth)
        end
        return
    elseif ndim > 2
        error("Only up to three dimensional arrays are currently supported")
    end

    if labels == nothing
        print(dimnames(na,1), " \\ ", dimnames(na,2),"\n")
    else
        print(varlab(labels,dimnames(na,1)), " \\ ", varlab(labels,dimnames(na,2)),"\n")
    end

    # all is true when row, col, or cell is true
    if all
        row = col = cell = true
    end

    # row and column names
    if labels == nothing
        rownames = string.(names(na,1))
        colnames = string.(names(na,2))
    else # labels defined
        rownames = findvlab(labels,dimnames(na,1),names(na,1))
        colnames = findvlab(labels,dimnames(na,2),names(na,2))
    end

    # maximum widths
    maxrowname = maximum(vcat(5,length.(rownames))) # minimum width is 5
    maxcolname = maximum(vcat(3,length.(colnames))) # minimum width is 3

    # width of data columns - the same as the greater of the length of the grand total
    # and 4 + precision
    tot = sum(na.array) # grand total
    colwidth = max(length(digits(Int(floor(tot)))),4+precision)

    # number of columns
    ncols = length(colnames)

    # number of rows
    nrows = length(rownames)

    # floating point numbers with three digits after decimal point
    if eltype(na.array) <: AbstractFloat
        colwidth += 3
    end

    # determine column widths
    colwidth = max(maxcolname,colwidth)

    #---------------------------------------------------
    # column totals
    colsum = sum(na.array,dims=1)

    # row totals
    rowsum = sum(na.array,dims=2)

    #---------------------------------------------------
    # determine how many columns can be printed within
    # the pagewidth
    maxcol = Int(floor((pagewidth - maxrowname) / colwidth))
    niter = Int(ceil(ncols / maxcol))
    nlastrow = ncols % maxcol

    for i in 1:niter

        #---------------------------------------------------
        # print header
        print(repeat(" ",maxrowname)," |")

        for j = 1:maxcol
            c = Int(maxcol*(i-1))+j
            if c > ncols
                break
            end
            print(" ",lpad(string(colnames[c]),colwidth))
        end

        print(" | ",lpad("Total",colwidth),"\n")

        print(repeat("-",maxrowname),"-+-",repeat("-",(colwidth+1)*(i == niter ? nlastrow : maxcol)),"+-",repeat("-",colwidth),"\n")

        #----------------------------------------------------
        # print values
        for r = 1:nrows

            # row name
            print(rpad(string(rownames[r]),maxrowname)," |")

            for j = 1:maxcol
                c = Int(maxcol*(i-1))+j
                if c > ncols
                    break
                end
                val = na.array[r,c]
                print(" ",lpad(val,colwidth))
            end

            # row total
            print(" |")

            val = rowsum[r]
            print(" ",lpad(val,colwidth),"\n")

            # row percentages
            if row
                print(repeat(" ",maxrowname)," |")
                for j = 1:maxcol
                    c = Int(maxcol*(i-1)) + j
                    if c > ncols
                        break
                    end
                    val = strval(100 * na.array[r,c] / rowsum[r],precision)
                    print(" ",lpad(val,colwidth," "))
                end

                # row percentage
                print(" |")

                val = strval(100.0,precision)
                print(" ",lpad(val,colwidth),"\n")
            end

            # column percentages
            if col
                print(repeat(" ",maxrowname)," |")
                for j = 1:maxcol
                    c = Int(maxcol*(i-1)) + j
                    if c > ncols
                        break
                    end
                    val = strval(100 * na.array[r,c] / colsum[c],precision)
                    print(" ",lpad(val,colwidth))
                end

                # column percent
                print(" |")

                val = strval(100 * rowsum[r] / tot,precision)
                print(" ",lpad(val,colwidth),"\n")
            end

            # cell percentages
            if cell
                print(repeat(" ",maxrowname)," |")
                for j = 1:maxcol
                    c = Int(maxcol*(i-1)) + j
                    if c > ncols
                        break
                    end
                    val = strval(100 * na.array[r,c] / tot,precision)
                    print(" ",lpad(val,colwidth))
                end

                # column percent
                print(" |")

                val = strval(100 * rowsum[r] / tot,precision)
                print(" ",lpad(val,colwidth),"\n")
            end

            if row || col || cell
                print(repeat("-",maxrowname),"-+-",repeat("-",(colwidth+1)*(i == niter ? nlastrow : maxcol)),"+-",repeat("-",colwidth),"\n")
            end

        end

        #----------------------------------------------------
        # Total
        if !(row || col || cell)
            print(repeat("-",maxrowname+1),"+",repeat("-",(colwidth+1)*(i == niter ? nlastrow : maxcol)),"-+-",repeat("-",colwidth),"\n")
        end
        print(rpad("Total",maxrowname," ")," |")

        for j = 1:maxcol
            c = Int(maxcol*(i-1)) + j
            if c > ncols
                break
            end
            val=colsum[c]
            print(" ",lpad(val,colwidth," "))
        end

        # Grand total
        val = strval(tot)
        print(" | ",lpad(val,colwidth," "),"\n")

        #----------------------------------------------------
        # row percentages
        if row
            print(repeat(" ",maxrowname)," |")
            for j = 1:maxcol
                c = Int(maxcol*(i-1)) + j
                if c > ncols
                    break
                end
                val = strval(100 * colsum[c] / tot,precision)
                print(" ",lpad(val,colwidth," "))
            end

            # column percent for the total column
            val = strval(100.0,precision)
            print(" | ",lpad(val,colwidth," "),"\n")
        end

        # column percentages
        if col
            print(repeat(" ",maxrowname)," |")
            for j = 1:maxcol
                c = Int(maxcol*(i-1)) + j
                if c > ncols
                    break
                end
                val = strval(100.0,precision)
                print(" ",lpad(val,colwidth," "))
            end

            # column percent for the total column
            val = strval(100.0,precision)
            print(" | ",lpad(val,colwidth," "),"\n")

        end

        # cell percentages
        if cell
            print(repeat(" ",maxrowname)," |")
            for j = 1:maxcol
                c = Int(maxcol*(i-1)) + j
                if c > ncols
                    break
                end
                val = strval(100.0 * colsum[c] / tot,precision)
                print(" ",lpad(val,colwidth," "))
            end

            # column percent
            val = strval(100.0,precision)
            print(" | ",lpad(val,colwidth," "),"\n")
        end

        # separator
        print("\n\n")
    end

    if chisq
        # chisq test
        dof = (nrows-1)*(ncols-1)
        chisqval = 0.
        for i = 1:nrows
            for j = 1:ncols
                expected = rowsum[i]*colsum[j]/tot
                chisqval += ((na.array[i,j] - expected)^2)/expected
            end
        end

        # degress of freedom
        dof = (ncols-1)*(nrows-1)

        # return a tuple of chisq, df, p-value
        pval = Distributions.ccdf(Distributions.Chisq(dof),chisqval)

        # chisquare output
        println("Pearson χ² (",dof,") = ",@sprintf("%.5f",chisqval)," Pr = ",@sprintf("%.5f",pval))
    end
end

function tabprint1(na::NamedArray; precision=2, labels::Union{Nothing,Label} = nothing)

    if ndims(na) == 2 && names(na,2)[1] != "Frequency"
        error("Cannot print 2 or higher dimensions")
    end

    # variable name
    dname = dimnames(na,1)

    # value names
    if labels == nothing
        rownames = string.(names(na,1))
    else
        rownames = findvlab(labels,dname,names(na,1)) # [ ismissing(x) ? "missing" : vallab(labels,dname,x) for x in names(na,1) ]
    end

    # maximum width for row labels
    if skipmissing == true && rownames[end] == "missing"
        maxrowname = maximum(vcat(5,length(string(dname)),length.(rownames[1:end-1]))) # minimum width is 5
    else
        maxrowname = maximum(vcat(5,length(string(dname)),length.(rownames))) # minimum width is 5
    end

    # number of rows
    nrows = length(rownames)

    # column names
    colnames = ["Frequency","Percent","Cum Percent"]

    # grand total
    if skipmissing == true && rownames[end] == "missing"
        tot = sum(na.array[1:end-1])
    else
        tot = sum(na)
    end

    # maximum column width
    maxcolwidth = maximum(vcat(length(string(tot)),length.(colnames)))

    # header
    print("\n")
    print(lpad(dname,maxrowname)," |")
    for i = 1:length(colnames)
        print(lpad(colnames[i],maxcolwidth+1))
    end
    print("\n")

    # separator
    print(repeat("-",maxrowname),"-+")
    for i = 1:length(colnames)
        print(repeat("-",maxcolwidth + 1))
    end
    print("\n")

    # rows
    cumpct = 0.
    for i = 1:nrows

        if skipmissing == true && rownames[i] == "missing"
            continue
        end

        print(lpad(rownames[i],maxrowname)," |")

        # frequency
        print(lpad(na.array[i],maxcolwidth + 1))

        # percent
        pct = na.array[i] / tot
        print(lpad(strval(100*pct,precision),maxcolwidth + 1))

        # cumulative percent
        cumpct += pct
        println(lpad(strval(100*cumpct,precision),maxcolwidth + 1))
    end

    # separator
    print(repeat("-",maxrowname),"-+")
    for i = 1:length(colnames)
        print(repeat("-",maxcolwidth + 1))
    end
    print("\n")

    # total
    print(lpad("Total",maxrowname)," |")

    # total frequency
    print(lpad(tot,maxcolwidth + 1))

    # percent
    print(lpad(strval(100.0,precision),maxcolwidth + 1))

    # cumulative percent
    println(lpad(strval(100.0,precision),maxcolwidth + 1))

    print("\n")
end

function findvlab(labels::Label,v::Symbol,val::Vector)

    lname = lblname(labels,v)

    if lname == nothing
        return val
    end

    return [ismissing(x) ? "missing" : vallab(labels,v,x) for x in val]
end

# function tab(x::AbstractVector...; allowmissing=false,weights::AbstractWeights = fweights(ones(Int8,length(x[1]))))
#
#     ncols = length(x)
#     nrows = length(x[1])
#     for i = 2:ncols
#         if nrows != length(x[i])
#             error("Columns do not have the same rows.")
#         end
#     end
#
#     vectypes = Missings.T.(eltype.(x))
#     println(vectypes)
#
#     # # create output arrays
#     # vdims = Vector{Vector}(ncols)
#     # vnums = zeros(Int64,ncols)
#     # naidx = falses(ncols)
#     # for i = 1:ncols
#     #     if sum(x[i].na) > 0 # there are NA values in this vector
#     #         vdims[i] = sort(levels(dropna(x[i])))
#     #         vnums[i] = size(vdims[i],1) + 1 # one for the NA
#     #         naidx[i] = true
#     #     else
#     #         vdims[i] = sort(levels(x[i]))
#     #         vnums[i] = size(vdims[i],1)
#     #     end
#     # end
#     # from FreqTables.jl/src/freqtable.jl
#     d = Dict{Tuple{Missings.T(eltype(a)),Missings.T(eltype(b))},Int64}()
#     @inbounds @simd for i in 1:alen
#         el = (a[i],b[i])
#         index = Base.ht_keyindex(d, el)
#
#         if index > 0
#             @inbounds d.vals[index] += weights[i]
#         else
#             @inbounds d[el] = weights[i]
#         end
#     end
#
#     # output array
#     cnt = Array{Int64,2}(length(alev),length(blev))
#     @inbounds @simd for i=1:length(alev)
#         @inbounds @simd for j=1:length(blev)
#             cnt[i,j] = d[alev[i],blev[j]]
#         end
#     end
#
#
#     return NamedArray(a, tuple(dimnames...), ntuple(i -> "Dim$i", ncols))
# end


text:
function g(n)
    a = rand(n,n)
    out = 0.0
    for _ =1:n
        for j=1:n
            for i=1:n
                out += sin(a[i,j])
            end
        end
    end
    out
end


text:
using JuliaEuler
ex24(1:3,4)

max_stars_repo_name: binnisb/JuliaEuler

max_issues_repo_name: binnisb/JuliaEuler

max_forks_repo_name: binnisb/JuliaEuler


text:
"""
nsolsc(f,x0, fp=difffp; rtol=1.e-6, atol=1.e-12, maxit=10,
        solver="newton", sham=1, armmax=10, resdec=.1, dx=1.e-7,
        armfix=false, pdata=nothing,
        printerr=true, keepsolhist=true, stagnationok=false)

C. T. Kelley, 2021

Newton's method for scalar equations. Has most of the features a
code for systems of equations needs. This is a wrapper for a call
to nsol.jl, the real code for systems.

Input:\n
f: function\n
x0: initial iterate\n
fp: derivative. If your derivative function is fp, you give me
its name. For example fp=foobar tells me that foobar is your
function for the derivative. The default is a forward difference
Jacobian that I provide.\n

Keyword Arguments (kwargs):\n
rtol, atol: real and absolute error tolerances\n

maxit: upper bound on number of nonlinear iterations\n

solver:\n
Your choices are "newton"(default) or "chord". However,
you have sham at your disposal only if you chose newton. "chord"
will keep using the initial derivative until the iterate converges,
uses the iteration budget, or the line search fails. It is not the
same as sham=Inf, which is smarter.\n

If you use secant and your initial iterate is poor, you have made
a mistake. I will help you by driving the line search with a finite
difference derivative.\n

sham:\n
This is the Shamanskii method. If sham=1, you have Newton.
The iteration updates the derivative every sham iterations.
The convergence rate has local q-order sham+1 if you only count
iterations where you update the derivative. You need not
provide your own derivative function to use this option. sham=Inf
is chord only if chord is converging well.\n

armmax: upper bound on stepsize reductions in linesearch

resdec: target value for residual reduction. \n
The default value is .1. In the old MATLAB codes it was .5.
I only turn Shamanskii on if the residuals are decreasing
rapidly, at least a factor of resdec, and the line search is quiescent.
If you want to eliminate resdec from the method ( you don't ) then set
resdec = 1.0 and you will never hear from it again.

dx:\n
This is the increment for forward difference, default = 1.e-7.
dx should be roughly the square root of the noise in the function.

armfix:\n
The default is a parabolic line search (ie false). Set to true and
the stepsize will be fixed at .5. Don't do this unless you are doing
experiments for research.

pdata:\n
precomputed data for the function/derivative.
Things will go better if you use this rather than hide the data
in global variables within the module for your function/derivative
If you use this option your function and derivative must take pdata
as a second argument. eg f(x,pdata) and fp(x,pdata)

printerr:\n
I print a helpful message when the solver fails. To supress that
message set printerr to false.

keepsolhist:\n
Set this to true to get the history of the iteration in the output
tuple. This is on by default for scalar equations and off for systems.
Only turn it on if you have use for the data, which can get REALLY LARGE.

stagnationok:\n
Set this to true if you want to disable the line search and either
observe divergence or stagnation. This is only useful for research
or writing a book.

Output:\n
A named tuple (solution, functionval, history, stats, idid,
               errcode, solhist)
where

solution = converged result
functionval = F(solution)
history = the vector of residual norms (||F(x)||) for the iteration
stats = named tuple of the history of (ifun, ijac, iarm), the number
of functions/derivatives/steplength reductions at each iteration.

I do not count the function values for a finite-difference derivative
because they count toward a Jacobian evaluation. I do count them for
the secant method model.

idid=true if the iteration succeeded and false if not.

errcode = 0 if if the iteration succeeded
        = -1 if the initial iterate satisifies the termination criteria
        = 10 if no convergence after maxit iterations
        = 1  if the line search failed

solhist:\n
This is the entire history of the iteration if you've set
keepsolhist=true\n
nsolsc builds solhist with a function from the Tools directory. For
systems, solhist is an N x K array where N is the length of x and K
is the number of iteration + 1. So, for scalar equations (N=1), solhist
is a row vector. Hence the use of solhist' in the example below.

### Examples for nsolsc.jl

```jldoctest
julia> nsolout=nsolsc(atan,1.0;maxit=5,atol=1.e-12,rtol=1.e-12);

julia> nsolout.history
6-element Array{Float64,1}:
 7.85398e-01
 5.18669e-01
 1.16332e-01
 1.06102e-03
 7.96200e-10
 2.79173e-24
```

# If you have an analytic derivative, I will use it.

```jldoctest
julia> fs(x)=x^2-4.0; fsp(x)=2x;

julia> nsolout=nsolsc(fs,1.0,fsp; maxit=5,atol=1.e-9,rtol=1.e-9);

julia> [nsolout.solhist'.-2 nsolout.history]
6×2 Array{Float64,2}:
 -1.00000e+00  3.00000e+00
  5.00000e-01  2.25000e+00
  5.00000e-02  2.02500e-01
  6.09756e-04  2.43940e-03
  9.29223e-08  3.71689e-07
  2.22045e-15  8.88178e-15

```

"""
function nsolsc(
    f,
    x0,
    fp = difffp;
    rtol = 1.e-6,
    atol = 1.e-12,
    maxit = 10,
    solver = "newton",
    sham = 1,
    armmax = 5,
    resdec = 0.1,
    dx = 1.e-7,
    armfix = false,
    pdata = nothing,
    printerr = true,
    keepsolhist = true,
    stagnationok = false,
)
    #
    # The scalar code is a simple wrapper for the real code (nsol). The
    # wrapper puts placeholders for the memory allocations and the precomputed
    # data.
    #
    fp0 = copy(x0)
    fpp0 = copy(x0)
    newtonout = nsol(
        f,
        x0,
        fp0,
        fpp0,
        fp;
        rtol = rtol,
        atol = atol,
        maxit = maxit,
        solver = solver,
        sham = sham,
        armmax = armmax,
        resdec = resdec,
        dx = dx,
        armfix = armfix,
        pdata = pdata,
        printerr = printerr,
        keepsolhist = keepsolhist,
        stagnationok = stagnationok,
    )
    return newtonout
end


text:
using PrettyPlots, Plots, DataFrames
gr();
#pgfplots();
#pyplot();

description = string(
"In this test: \n",c.m.name,"\n
* m.Ni=",c.m.Ni," \n
* m.Nck=",c.m.Nck,"\n
* m.tp=",c.m.tp," \n
* m.tex=",c.m.tex,"\n
* m.max_cpu_time=",c.m.max_cpu_time," \n
* m.Nck=",c.m.Nck,"\n
* m.PredictX0=",c.m.PredictX0," \n
* m.FixedTp=",c.m.FixedTp,"\n
")
results_dir=string("testE_",c.m.name,"_",c.m.solver,"/")
resultsDir!(r,results_dir;description=description);

println("Plotting the Final Results!")
s=Settings(;reset=false,save=true,simulate=true,MPC=false,format=:png);
#pmain(n,r,s,c)
#pp(n,r,s,c)
mainSimPath(n,r,s,c,pa)
if r.dfs_opt[r.eval_num-1][:status]==:Infeasible
  s=Settings(;evalConstraints=true,save=true,MPC=false,simulate=false,format=:png);
  postProcess(n,r,s)
  # trouble getting a feasible solution? -> look at the constraints
  print(r.constraint.value)
end # then consider relaxing tolerances etc. to make it work


text:
@muladd function ode_determine_initdt(u0,t,tdir,dtmax,abstol,reltol,internalnorm,prob::DiffEqBase.AbstractODEProblem{uType,tType,true},integrator) where {tType,uType}
  _tType = eltype(tType)
  f = prob.f
  p = integrator.p
  oneunit_tType = oneunit(_tType)
  dtmax_tdir = tdir*dtmax

  if eltype(u0) <: Number && !(typeof(integrator.alg) <: CompositeAlgorithm)
    cache = get_tmp_cache(integrator)
    sk = first(cache)
    @.. sk = abstol+internalnorm(u0,t)*reltol
  else
    sk = @.. abstol+internalnorm(u0,t)*reltol
  end

  if get_current_isfsal(integrator.alg, integrator.cache) && typeof(integrator) <: ODEIntegrator
    # Right now DelayDiffEq has issues with fsallast not being initialized
    f₀ = integrator.fsallast
    f(f₀,u0,p,t)
  else
    # TODO: use more caches
    f₀ = zero.(u0./t)
    f(f₀,u0,p,t)
  end

  # TODO: use more caches
  #tmp = cache[2]
  tmp = @.. u0/sk

  d₀ = internalnorm(tmp,t)

  #=
  Try/catch around the linear solving. This will catch singular matrices defined
  by DAEs and thus we use the _tType(1//10^(6)) default from Hairer. Note that
  this will not always catch singular matrices, an example from Andreas:

  julia> A = fill(rand(), 2, 2)
  2×2 Array{Float64,2}:
   0.637947  0.637947
   0.637947  0.637947

  julia> inv(A)
  2×2 Array{Float64,2}:
    9.0072e15  -9.0072e15
   -9.0072e15   9.0072e15

  The only way to make this more correct is to check

  issingular(A) = rank(A) < min(size(A)...)

  but that would introduce another svdfact in rank (which may not be possible
  anyways if the mass_matrix is not actually an array). Instead we stick to the
  user-chosen factorization. Sometimes this will cause `ftmp` to be absurdly
  large like shown there, but that later gets caught in the quick estimates
  below which then makes it spit out the default

  dt₀ = _tType(1//10^(6))

  so that is a a cheaper way to get to the same place than an svdfact and
  still works for matrix-free definitions of the mass matrix.
  =#

  dtmin = nextfloat(integrator.opts.dtmin)
  smalldt = convert(_tType,oneunit_tType*1//10^(6))
  if prob.f.mass_matrix != I && (!(typeof(prob.f)<:DynamicalODEFunction) || any(mm != I for mm in prob.f.mass_matrix))
    ftmp = zero(f₀)
    try
      integrator.alg.linsolve(ftmp, copy(prob.f.mass_matrix), f₀, true)
      f₀ .= ftmp
    catch
      return tdir*max(smalldt, dtmin)
    end
  end

  if integrator.opts.verbose && any(x->any(isnan, x), f₀)
    @warn("First function call produced NaNs. Exiting.")
  end

  @.. tmp = f₀/sk*oneunit_tType
  d₁ = internalnorm(tmp,t)

  if d₀ < 1//10^(5) || d₁ < 1//10^(5)
    dt₀ = smalldt
  else
    dt₀ = convert(_tType,oneunit_tType*(d₀/d₁)/100)
  end
  dt₀ = min(dt₀,dtmax_tdir)

  if typeof(one(_tType)) <: AbstractFloat && dt₀ < 10eps(_tType)*oneunit(_tType)
    # This catches Andreas' non-singular example
    # should act like it's singular
    return tdir*max(smalldt, dtmin)
  end

  dt₀_tdir = tdir*dt₀

  u₁ = zero(u0) # required by DEDataArray
  @.. u₁ = u0 + dt₀_tdir*f₀
  f₁ = zero(f₀)
  f(f₁,u₁,p,t+dt₀_tdir)

  if prob.f.mass_matrix != I && (!(typeof(prob.f)<:DynamicalODEFunction) || any(mm != I for mm in prob.f.mass_matrix))
    integrator.alg.linsolve(ftmp, prob.f.mass_matrix, f₁, false)
    f₁ .= ftmp
  end

  # Constant zone before callback
  # Just return first guess
  # Avoids AD issues
  f₀ == f₁ && return tdir*max(dtmin, 100dt₀)

  @.. tmp = (f₁-f₀)/sk*oneunit_tType
  d₂ = internalnorm(tmp,t)/dt₀*oneunit_tType
  # Hairer has d₂ = sqrt(sum(abs2,tmp))/dt₀, note the lack of norm correction

  max_d₁d₂ = max(d₁,d₂)
  if max_d₁d₂ <= 1//Int64(10)^(15)
    dt₁ = max(convert(_tType,oneunit_tType*1//10^(6)),dt₀*1//10^(3))
  else
    dt₁ = convert(_tType,oneunit_tType*10.0^(-(2+log10(max_d₁d₂))/get_current_alg_order(integrator.alg,integrator.cache)))
  end
  return tdir*max(dtmin, min(100dt₀,dt₁,dtmax_tdir))
end

@muladd function ode_determine_initdt(u0,t,tdir,dtmax,abstol,reltol,internalnorm,prob::DiffEqBase.AbstractODEProblem{uType,tType,false},integrator) where {uType,tType}
  _tType = eltype(tType)
  f = prob.f
  p = prob.p
  oneunit_tType = oneunit(_tType)
  dtmax_tdir = tdir*dtmax

  sk = @.. abstol + internalnorm(u0,t) * reltol
  d₀ = internalnorm(u0 ./ sk,t)

  f₀ = f(u0,p,t)
  if integrator.opts.verbose && any(x->any(isnan, x), f₀)
    @warn("First function call produced NaNs. Exiting.")
  end

  d₁ = internalnorm(f₀ ./ sk .* oneunit_tType,t)

  dtmin = nextfloat(integrator.opts.dtmin)
  smalldt = convert(_tType,oneunit_tType*1//10^(6))

  if d₀ < 1//10^(5) || d₁ < 1//10^(5)
    dt₀ = smalldt
  else
    dt₀ = convert(_tType,oneunit_tType*(d₀/d₁)/100)
  end
  dt₀ = min(dt₀,dtmax_tdir)
  dt₀_tdir = tdir*dt₀

  u₁ = @.. u0 + dt₀_tdir * f₀
  f₁ = f(u₁,p,t+dt₀_tdir)

  # Constant zone before callback
  # Just return first guess
  # Avoids AD issues
  f₀ == f₁ && return tdir*max(dtmin, 100dt₀)

  d₂ = internalnorm((f₁ .- f₀) ./ sk .* oneunit_tType,t) / dt₀ * oneunit_tType

  max_d₁d₂ = max(d₁, d₂)
  if max_d₁d₂ <= 1//Int64(10)^(15)
    dt₁ = max(smalldt,dt₀*1//10^(3))
  else
    dt₁ = _tType(oneunit_tType*10.0^(-(2+log10(max_d₁d₂))/get_current_alg_order(integrator.alg,integrator.cache)))
  end
  return tdir*max(dtmin, min(100dt₀,dt₁,dtmax_tdir))
end

@inline function ode_determine_initdt(u0,t,tdir,dtmax,abstol,reltol,internalnorm,prob::DiffEqBase.AbstractDAEProblem{duType,uType,tType},integrator) where {duType,uType,tType}
  _tType = eltype(tType)
  oneunit_tType = oneunit(_tType)
  return convert(_tType,oneunit_tType*1//10^(6))
end


text:
"""
    RootSolvers

Contains functions for solving roots of non-linear
equations. See [`find_zero`](@ref).

## Example

```julia
julia> using RootSolvers

julia> sol = find_zero(x -> x^2 - 100^2,
                       SecantMethod{Float64}(0.0, 1000.0),
                       CompactSolution());

julia> sol
RootSolvers.CompactSolutionResults{Float64}(99.99999999994358, true)
```

"""
module RootSolvers

using DocStringExtensions: FIELDS

export find_zero,
    SecantMethod, RegulaFalsiMethod, NewtonsMethodAD, NewtonsMethod
export CompactSolution, VerboseSolution
export AbstractTolerance, ResidualTolerance, SolutionTolerance

import ForwardDiff

# Input types
const FTypes = Union{AbstractFloat, AbstractArray}

abstract type RootSolvingMethod{FT <: FTypes} end
Base.broadcastable(method::RootSolvingMethod) = Ref(method)

"""
    SecantMethod

# Fields
$(FIELDS)
"""
struct SecantMethod{FT} <: RootSolvingMethod{FT}
    "lower"
    x0::FT
    "upper bound"
    x1::FT
end

"""
    RegulaFalsiMethod

# Fields
$(FIELDS)
"""
struct RegulaFalsiMethod{FT} <: RootSolvingMethod{FT}
    "lower bound"
    x0::FT
    "upper bound"
    x1::FT
end

"""
    NewtonsMethodAD

# Fields
$(FIELDS)
"""
struct NewtonsMethodAD{FT} <: RootSolvingMethod{FT}
    "initial guess"
    x0::FT
end

"""
    NewtonsMethod

# Fields
$(FIELDS)
"""
struct NewtonsMethod{FT, F′ <: Function} <: RootSolvingMethod{FT}
    "initial guess"
    x0::FT
    "`f′` derivative of function `f` whose zero is sought"
    f′::F′
end

abstract type SolutionType end
Base.broadcastable(soltype::SolutionType) = Ref(soltype)

"""
    VerboseSolution <: SolutionType

Used to return a [`VerboseSolutionResults`](@ref)
"""
struct VerboseSolution <: SolutionType end

abstract type AbstractSolutionResults{AbstractFloat} end

"""
    VerboseSolutionResults{FT} <: AbstractSolutionResults{FT}

Result returned from `find_zero` when
`VerboseSolution` is passed as the `soltype`.
"""
struct VerboseSolutionResults{FT} <: AbstractSolutionResults{FT}
    "solution ``x^*`` of the root of the equation ``f(x^*) = 0``"
    root::FT
    "indicates convergence"
    converged::Bool
    "error of the root of the equation ``f(x^*) = 0``"
    err::FT
    "number of iterations performed"
    iter_performed::Int
    "solution per iteration"
    root_history::Vector{FT}
    "error of the root of the equation ``f(x^*) = 0`` per iteration"
    err_history::Vector{FT}
end
SolutionResults(soltype::VerboseSolution, args...) =
    VerboseSolutionResults(args...)

"""
    CompactSolution <: SolutionType

Used to return a [`CompactSolutionResults`](@ref)
"""
struct CompactSolution <: SolutionType end

"""
    CompactSolutionResults{FT} <: AbstractSolutionResults{FT}

Result returned from `find_zero` when
`CompactSolution` is passed as the `soltype`.
"""
struct CompactSolutionResults{FT} <: AbstractSolutionResults{FT}
    "solution ``x^*`` of the root of the equation ``f(x^*) = 0``"
    root::FT
    "indicates convergence"
    converged::Bool
end
SolutionResults(soltype::CompactSolution, root, converged, args...) =
    CompactSolutionResults(root, converged)

init_history(::VerboseSolution, x::FT) where {FT <: AbstractFloat} = FT[x]
init_history(::CompactSolution, x) = nothing
init_history(::VerboseSolution, ::Type{FT}) where {FT <: AbstractFloat} = FT[]
init_history(::CompactSolution, ::Type{FT}) where {FT <: AbstractFloat} =
    nothing

function push_history!(
    history::Vector{FT},
    x::FT,
    ::VerboseSolution,
) where {FT <: AbstractFloat}
    push!(history, x)
end
function push_history!(
    history::Nothing,
    x::FT,
    ::CompactSolution,
) where {FT <: AbstractFloat}
    nothing
end

abstract type AbstractTolerance{FT <: FTypes} end
Base.broadcastable(tol::AbstractTolerance) = Ref(tol)

"""
    ResidualTolerance

A tolerance type based on the residual of the equation ``f(x) = 0``
"""
struct ResidualTolerance{FT} <: AbstractTolerance{FT}
    tol::FT
end

"""
    (tol::ResidualTolerance)(x1, x2, y)

Evaluates residual tolerance, based on ``|f(x)|``
"""
(tol::ResidualTolerance)(x1, x2, y) = abs(y) < tol.tol

"""
    SolutionTolerance

A tolerance type based on the solution ``x`` of the equation ``f(x) = 0``
"""
struct SolutionTolerance{FT} <: AbstractTolerance{FT}
    tol::FT
end

"""
    (tol::SolutionTolerance)(x1, x2, y)

Evaluates solution tolerance, based on ``|x2-x1|``
"""
(tol::SolutionTolerance)(x1, x2, y) = abs(x2 - x1) < tol.tol

# TODO: CuArrays.jl has trouble with isapprox on 1.1
# we use simple checks for now, will switch to relative checks later.

"""
    sol = find_zero(
            f::F,
            method::RootSolvingMethod{FT},
            soltype::SolutionType,
            tol::Union{Nothing, AbstractTolerance} = nothing,
            maxiters::Union{Nothing, Int} = 10_000,
            )

Finds the nearest root of `f`. Returns a the value of the root `x` such
that `f(x) ≈ 0`, and a Boolean value `converged` indicating convergence.

 - `f` function of the equation to find the root
 - `method` can be one of:
    - `SecantMethod()`: [Secant method](https://en.wikipedia.org/wiki/Secant_method)
    - `RegulaFalsiMethod()`: [Regula Falsi method](https://en.wikipedia.org/wiki/False_position_method#The_regula_falsi_(false_position)_method)
    - `NewtonsMethodAD()`: [Newton's method](https://en.wikipedia.org/wiki/Newton%27s_method) using Automatic Differentiation
    - `NewtonsMethod()`: [Newton's method](https://en.wikipedia.org/wiki/Newton%27s_method)
- `soltype` is a solution type which may be one of:
      `CompactSolution` GPU-capable. Solution has `converged` and `root` only, see [`CompactSolutionResults`](@ref)
      `VerboseSolution` CPU-only. Solution has additional diagnostics, see [`VerboseSolutionResults`](@ref)
- `tol` is a tolerance type to determine when to stop iterations.
- `maxiters` is the maximum number of iterations.
"""
function find_zero end

# Main entry point: Dispatch to specific method
function find_zero(
    f::F,
    method::RootSolvingMethod{FT},
    soltype::SolutionType,
    tol::Union{Nothing, AbstractTolerance} = nothing,
    maxiters::Union{Nothing, Int} = 10_000,
) where {FT <: FTypes, F <: Function}
    if tol === nothing
        tol = SolutionTolerance{eltype(FT)}(1e-3)
    end
    return find_zero(f, method, method_args(method)..., soltype, tol, maxiters)
end

# Allow broadcast:
function Broadcast.broadcasted(
    ::typeof(find_zero),
    f::F,
    method::RootSolvingMethod{FT},
    soltype::SolutionType,
    tol::Union{Nothing, AbstractTolerance} = nothing,
    maxiters::Union{Nothing, Int} = 10_000,
) where {FT <: FTypes, F}
    if tol === nothing
        tol = SolutionTolerance{eltype(FT)}(1e-3)
    end
    return broadcast(
        find_zero,
        f,
        method,
        method_args(method)...,
        soltype,
        tol,
        maxiters,
    )
end

####
#### Numerical methods
####

"""
    method_args(method::RootSolvingMethod)

Return tuple of positional args for `RootSolvingMethod`.
"""
function method_args end

method_args(method::SecantMethod) = (method.x0, method.x1)

function find_zero(
    f::F,
    ::SecantMethod,
    x0::FT,
    x1::FT,
    soltype::SolutionType,
    tol::AbstractTolerance{FT},
    maxiters::Int,
) where {F <: Function, FT <: FTypes}
    y0 = f(x0)
    y1 = f(x1)
    x_history = init_history(soltype, x0)
    y_history = init_history(soltype, y0)
    for i in 1:maxiters
        Δx = x1 - x0
        Δy = y1 - y0
        x0, y0 = x1, y1
        push_history!(x_history, x0, soltype)
        push_history!(y_history, y0, soltype)
        x1 -= y1 * Δx / Δy
        y1 = f(x1)
        if tol(x0, x1, y1)
            return SolutionResults(
                soltype,
                x1,
                true,
                y1,
                i,
                x_history,
                y_history,
            )
        end
    end
    return SolutionResults(
        soltype,
        x1,
        false,
        y1,
        maxiters,
        x_history,
        y_history,
    )
end

method_args(method::RegulaFalsiMethod) = (method.x0, method.x1)

function find_zero(
    f::F,
    ::RegulaFalsiMethod,
    x0::FT,
    x1::FT,
    soltype::SolutionType,
    tol::AbstractTolerance{FT},
    maxiters::Int,
) where {F <: Function, FT}
    y0 = f(x0)
    y1 = f(x1)
    @assert y0 * y1 < 0
    x_history = init_history(soltype, x0)
    y_history = init_history(soltype, y0)
    lastside = 0
    local x, y
    for i in 1:maxiters
        x = (x0 * y1 - x1 * y0) / (y1 - y0)
        y = f(x)
        push_history!(x_history, x, soltype)
        push_history!(y_history, y, soltype)
        if y * y0 < 0
            if tol(x, x1, y)
                return SolutionResults(
                    soltype,
                    x,
                    true,
                    y,
                    i,
                    x_history,
                    y_history,
                )
            end
            x1, y1 = x, y
            if lastside == +1
                y0 /= 2
            end
            lastside = +1
        else
            if tol(x0, x, y)
                return SolutionResults(
                    soltype,
                    x,
                    true,
                    y,
                    i,
                    x_history,
                    y_history,
                )
            end
            x0, y0 = x, y
            if lastside == -1
                y1 /= 2
            end
            lastside = -1
        end
    end
    return SolutionResults(soltype, x, false, y, maxiters, x_history, y_history)
end

"""
    value_deriv(f, x)

Compute the value and derivative `f(x)` using ForwardDiff.jl.
"""
function value_deriv(f, x::FT) where {FT}
    tag = typeof(ForwardDiff.Tag(f, FT))
    y = f(ForwardDiff.Dual{tag}(x, one(x)))
    ForwardDiff.value(tag, y), ForwardDiff.partials(tag, y, 1)
end

method_args(method::NewtonsMethodAD) = (method.x0,)

function find_zero(
    f::F,
    ::NewtonsMethodAD,
    x0::FT,
    soltype::SolutionType,
    tol::AbstractTolerance{FT},
    maxiters::Int,
) where {F <: Function, FT}
    local y
    x_history = init_history(soltype, FT)
    y_history = init_history(soltype, FT)
    if soltype isa VerboseSolution
        y, y′ = value_deriv(f, x0)
        push_history!(x_history, x0, soltype)
        push_history!(y_history, y, soltype)
    end
    for i in 1:maxiters
        y, y′ = value_deriv(f, x0)
        x1 = x0 - y / y′
        push_history!(x_history, x1, soltype)
        push_history!(y_history, y, soltype)
        if tol(x0, x1, y)
            return SolutionResults(
                soltype,
                x1,
                true,
                y,
                i,
                x_history,
                y_history,
            )
        end
        x0 = x1
    end
    return SolutionResults(
        soltype,
        x0,
        false,
        y,
        maxiters,
        x_history,
        y_history,
    )
end

method_args(method::NewtonsMethod) = (method.x0, method.f′)

function find_zero(
    f::F,
    ::NewtonsMethod,
    x0::FT,
    f′::F′,
    soltype::SolutionType,
    tol::AbstractTolerance{FT},
    maxiters::Int,
) where {F <: Function, F′ <: Function, FT, IT <: Int}
    x_history = init_history(soltype, FT)
    y_history = init_history(soltype, FT)
    if soltype isa VerboseSolution
        y, y′ = f(x0), f′(x0)
        push_history!(x_history, x0, soltype)
        push_history!(y_history, y, soltype)
    end
    for i in 1:maxiters
        y, y′ = f(x0), f′(x0)
        x1 = x0 - y / y′
        push_history!(x_history, x1, soltype)
        push_history!(y_history, y, soltype)
        if tol(x0, x1, y)
            return SolutionResults(
                soltype,
                x1,
                true,
                y,
                i,
                x_history,
                y_history,
            )
        end
        x0 = x1
    end
    return SolutionResults(
        soltype,
        x0,
        false,
        y,
        maxiters,
        x_history,
        y_history,
    )
end

end


text:
module MouseButtons
export MouseButton
@enum MouseButton::UInt8 begin
    LEFT = 1
    MIDDLE = 2
    RIGHT = 3
    WHEEL_UP = 4
    WHEEL_DOWN = 5
end
end

using .MouseButtons

# This somewhat convoluted code below is to map julia code
# like `Keys.A` to the SDLK_* enums, without actually enumerating
# all the values in Julia. We used to do that, but that duplication
# was somewhat problematic.

# A Keys.X invocation returns the keyscan code (an UInt32) corresponding
# to the SDLK_X enum defined in SDL2.
# In contrast, a keyboard.X returns true or false depending of whether the X
# key is currently pressed.
# Thus user code to check if X is pressed is written as
#       ` if game.keyboard.X`
const singleCharStrings = string.(collect('a':'z'))
struct KeyHolder; end
const Keys = KeyHolder()
function Base.getproperty(k::KeyHolder, s::Symbol)
    st=string(s)
    if findfirst(x->x==lowercase(st), singleCharStrings) !== nothing
        st = lowercase(st)
    end
    st = "SDLK_" * st
    s=Symbol(st)
    try
        return getproperty(GameZero.SimpleDirectMediaLayer.LibSDL2, s)
    catch
        @error "Unknown key: $s"
        return nothing
    end
end

struct Keyboard
    pressed::Array{UInt32, 1}
end

Keyboard() = Keyboard(Array{Any, 1}())

function Base.getproperty(k::Keyboard, s::Symbol)
    return getproperty(GameZero.Keys, s) in getfield(k, :pressed)
end

Base.push!(k::Keyboard, item) = push!(getfield(k, :pressed), item)
function Base.delete!(k::Keyboard, item)
    a = getfield(k, :pressed)
    deleteat!(a, a.==item)
end


text:
module RDF

using URIParser

export
    # types
    Graph,
    # functions (graph manipulation)
    push!,
    pop!,
    # functions (graph serialization)
    ntriples,
    nquads,
    turtle

type Graph
    name::URI
    size::Int64
    statements::Dict{String,Dict{String,Set{Union(Number,String,URI)}}}

    Graph(name::URI) = new(name, 0, Dict{String,Dict{String,Set{Union(Number,String,URI)}}}())
end

function push!(graph::Graph,
               subject::URI,
               predicate::URI,
               object::Union(Number,String,URI))
    # Get dict mappings:
    subject_dict = get(graph.statements, string(subject), Dict{String,Set{Union(Number,String,URI)}}())
    predicate_set = get(subject_dict, string(predicate), Set{Union(Number,String,URI)}())

    # Add statement:
    objects_num = length(predicate_set)
    Base.push!(predicate_set, object)
    if objects_num != length(predicate_set)
        graph.size += 1
    end

    # Store newly created dicts/sets:
    if !haskey(subject_dict, object)
        subject_dict[string(predicate)] = predicate_set
    end
    if !haskey(graph.statements, subject)
        graph.statements[string(subject)] = subject_dict
    end
end

function pop!(graph::Graph,
              subject::URI)
    # Get dict mappings;
    if !haskey(graph.statements, string(subject))
        return 0
    end
    subject_dict = Base.pop!(graph.statements, string(subject))

    # Remove statements:
    removed_statements = 0
    for predicate_objects = subject_dict
        removed_statements += length(predicate_objects[2])
    end
    graph.size -= removed_statements

    # Final number of removed statements:
    return removed_statements
end

function pop!(graph::Graph,
              subject::URI,
              predicate::URI)
    # Get dict mappings:
    if !haskey(graph.statements, string(subject))
        return 0
    end
    subject_dict = graph.statements[string(subject)]
    if !haskey(subject_dict, string(predicate))
        return 0
    end

    # Remove statements:
    objects = Base.pop!(subject_dict, string(predicate))
    removed_statements = length(objects)
    graph.size -= removed_statements

    # Remove empty dict:
    if length(subject_dict) == 0
        Base.pop!(graph.statements, string(subject))
    end

    # Final number of removed statements:
    return removed_statements
end

function pop!(graph::Graph,
              subject::URI,
              predicate::URI,
              object::Union(Number,String,URI))
    # Get dict mappings:
    if !haskey(graph.statements, string(subject))
        return 0
    end
    subject_dict = get(graph.statements, string(subject))
    if !haskey(subject_dict, string(predicate))
        return 0
    end
    predicate_set = get(subject_dict, string(predicate))

    # Remove statement:
    Base.pop!(predicate_set, object)
    graph.size -= 1

    # Remove empty dicts/sets:
    if length(predicate_set) == 0
        Base.pop!(subject_dict, string(predicate))
    end
    if length(subject_dict) == 0
        Base.pop!(graph.statements, string(subject))
    end

    # Final number of removed statements:
    return 1
end

function ntriples(graph::Graph,
                  out::Any)
    for statement = graph.statements
        for predicate_object = statement[2]
            for object = predicate_object[2]
                write(out, join(triple(statement[1], predicate_object[1], object)), " .\n")
            end
        end
    end
end

function nquads(graph::Graph,
                out::Any)
    for statement = graph.statements
        for predicate_object = statement[2]
            for object = predicate_object[2]
                write(out, join(triple(statement[1], predicate_object[1], object)), " <", string(graph.name), "> .\n")
            end
        end
    end
end

function turtle(graph::Graph,
                out::Any)
    subject_same = false
    predicate_same = false
    last_subject = nothing
    last_predicate = nothing
    last_object = nothing
    for statement = graph.statements
        for predicate_object = statement[2]
            for object = predicate_object[2]
                # Statement:
                subject = statement[1]
                predicate = predicate_object[1]

                # Serialize:
                prefix = ""
                if subject_same
                    prefix *= "    "
                elseif last_subject != nothing
                    prefix *= uri_rdf(last_subject) * " "
                end
                if predicate_same
                    prefix *= "    "
                elseif last_predicate != nothing
                    prefix *= uri_rdf(last_predicate) * " "
                end
                subject_same = last_subject == subject
                predicate_same = subject_same & (last_predicate == predicate)
                if subject_same
                    if predicate_same
                        write(out, prefix, object_rdf(last_object), " ,\n")
                    else
                        write(out, prefix, object_rdf(last_object), " ;\n")
                    end
                elseif last_subject != nothing
                    write(out, prefix, object_rdf(last_object), " .\n")
                end

                # Remember this statement, which only determined `separator`, but has not been serialized yet:
                last_subject = subject
                last_predicate = predicate
                last_object = object
            end
        end
    end

    # Serialize saved statement, if exists:
    if last_subject != nothing
        if subject_same
            if predicate_same
                write(out, "        ", object_rdf(last_object), " .\n")
            else
                write(out, "    ", uri_rdf(last_predicate), " ", object_rdf(last_object), " .\n")
            end
        else
            write(out, join(triple(last_subject, last_predicate, last_object)), " .\n")
        end
    end
end

function triple(subject::String,
                predicate::String,
                object::Union(String,Number,URI))
    return uri_rdf(subject), " ", uri_rdf(predicate), " ", object_rdf(object)
end

function uri_rdf(uri::String)
    # TODO Not complete yet; needs to handle URI encoding/escaping.
    return string("<", uri, ">")
end

function object_rdf(value::String)
    # TODO I wonder whether there is a simpler way to do this...
    # TODO Also, this function is not complete yet.
    previous = nothing
    escaped_string = utf8("")
    for character = value
        if previous == nothing && character == '"'
            escaped_string *= "\\\""
        elseif character == '"' && previous != '"'
            escaped_string *= "\\\""
        else
            escaped_string *= string(character)
        end
        previous = character
    end
    return "\"" * escaped_string * "\""
end

function object_rdf(value::Number)
    # TODO Add RDF/RDFS type based on the particular number we're seeing.
    return string(value)
end

function object_rdf(value::URI)
    return uri_rdf(string(value))
end

end # module RDF


text:
function bubble_sort!(array::Vector{T}) where {T}
    for i in 1:length(array)
        flag = false
        for j in 1:(length(array)-i)
            if array[j] > array[j+1]
                array[j], array[j+1] = array[j+1], array[j]
                flag = true
            end
        end

        if !flag
            return nothing
        end
    end
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
"""
    ast_is(ast::VPtr, what::Symbol)::Bool

Helper for quickly recognizing kinds of ASTs
"""
ast_is(ast::VPtr, what::Symbol)::Bool = ccall(sbml(what), Cint, (VPtr,), ast) != 0

"""
    parse_math_children(ast::VPtr)::Vector{Math}

Recursively parse all children of an AST node.
"""
parse_math_children(ast::VPtr)::Vector{Math} = [
    parse_math(ccall(sbml(:ASTNode_getChild), VPtr, (VPtr, Cuint), ast, i - 1)) for
    i = 1:ccall(sbml(:ASTNode_getNumChildren), Cuint, (VPtr,), ast)
]

"""
    parse_math(ast::VPtr)::Math

This attempts to parse out a decent Julia-esque ([`Math`](@ref) AST from a
pointer to `ASTNode_t`.
"""
function parse_math(ast::VPtr)::Math
    if ast_is(ast, :ASTNode_isName)
        if ccall(sbml(:ASTNode_getType), Cint, (VPtr,), ast) == 262
            # This is a special case checking for the value of "simulation
            # time" as defined by SBML. The constant `262` is the value of the
            # enum AST_NAME_TIME in `libsbml/src/sbml/math/ASTNodeType.h`,
            # needs to be kept up to date with the library (otherwise this
            # breaks).
            return MathTime(get_string(ast, :ASTNode_getName))
        else
            return MathIdent(get_string(ast, :ASTNode_getName))
        end
    elseif ast_is(ast, :ASTNode_isConstant)
        return MathConst(get_string(ast, :ASTNode_getName))
    elseif ast_is(ast, :ASTNode_isInteger)
        return MathVal(ccall(sbml(:ASTNode_getInteger), Cint, (VPtr,), ast))
    elseif ast_is(ast, :ASTNode_isReal)
        return MathVal(ccall(sbml(:ASTNode_getReal), Cdouble, (VPtr,), ast))
    elseif ast_is(ast, :ASTNode_isFunction)
        return MathApply(get_string(ast, :ASTNode_getName), parse_math_children(ast))
    elseif ast_is(ast, :ASTNode_isOperator) || ast_is(ast, :ASTNode_isRelational)
        return MathApply(
            string(Char(ccall(sbml(:ASTNode_getCharacter), Cchar, (VPtr,), ast))),
            parse_math_children(ast),
        )
    elseif ast_is(ast, :ASTNode_isLambda)
        children = parse_math_children(ast)
        if !isempty(children)
            body = pop!(children)
            return MathLambda(broadcast((x::MathIdent) -> x.id, children), body)
        else
            @warn "invalid function definition found"
            return MathIdent("?invalid?")
        end
    else
        @warn "unsupported math element found"
        return MathIdent("?unsupported?")
    end
end


text:
@testset "notequal.jl" begin
    @testset "NotEqualConstant()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(2, 6, "x", trailer)

        constraint = SeaPearl.NotEqualConstant(x, 3, trailer)

        @test constraint in x.onDomainChange
        @test constraint.active.value
    end
    @testset "propagate!(::NotEqualConstant)" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(2, 6, "x", trailer)

        constraint = SeaPearl.NotEqualConstant(x, 3, trailer)

        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test length(x.domain) == 4
        @test 2 in x.domain
        @test !(3 in x.domain)
        @test prunedDomains == SeaPearl.CPModification("x" => [3])

        y = SeaPearl.IntVar(2, 2, "y", trailer)

        cons2 = SeaPearl.NotEqualConstant(y, 2, trailer)

        @test !SeaPearl.propagate!(cons2, toPropagate, prunedDomains)

        @test isempty(y.domain)

        z = SeaPearl.IntVar(2, 6, "z", trailer)
        constraint1 = SeaPearl.NotEqualConstant(z, 3, trailer)
        constraint2 = SeaPearl.NotEqualConstant(z, 4, trailer)

        toPropagate2 = Set{SeaPearl.Constraint}()
        SeaPearl.propagate!(constraint1, toPropagate2, prunedDomains)

        @test constraint2 in toPropagate2

    end

    @testset "propagate!(::NotEqual)" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(2, 6, "x", trailer)
        y = SeaPearl.IntVar(5, 8, "y", trailer)

        constraint = SeaPearl.NotEqual(x, y, trailer)
        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test length(x.domain) == 5
        @test length(y.domain) == 4
        @test prunedDomains == SeaPearl.CPModification()

        # Propagation test
        z = SeaPearl.IntVar(5, 6, "z", trailer)
        constraint2 = SeaPearl.NotEqual(y, z, trailer)
        @test SeaPearl.propagate!(constraint2, toPropagate, prunedDomains)

        # Domain not reduced => not propagation
        @test !(constraint in toPropagate)
        @test !(constraint2 in toPropagate)

        # Domain reduced => propagation
        SeaPearl.remove!(z.domain, 6)
        @test SeaPearl.propagate!(constraint2, toPropagate, prunedDomains)
        @test constraint in toPropagate
        @test !(constraint2 in toPropagate)
        @test prunedDomains == SeaPearl.CPModification("y" => [5])
        @test length(y.domain) == 3

        #Unfeasible test
        t = SeaPearl.IntVar(5, 5, "t", trailer)
        constraint3 = SeaPearl.NotEqual(z, t, trailer)
        @test !SeaPearl.propagate!(constraint3, toPropagate, prunedDomains)

        # Same with IntVarViewOpposite
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(2, 6, "x", trailer)
        y = SeaPearl.IntVar(-8, -5, "y", trailer)
        minusY = SeaPearl.IntVarViewOpposite(y, "-y")

        constraint = SeaPearl.NotEqual(x, minusY, trailer)
        toPropagate = Set{SeaPearl.Constraint}()
        prunedDomains = SeaPearl.CPModification()

        @test SeaPearl.propagate!(constraint, toPropagate, prunedDomains)

        @test length(x.domain) == 5
        @test length(y.domain) == 4
        @test prunedDomains == SeaPearl.CPModification()

        # Propagation test
        z = SeaPearl.IntVar(5, 6, "z", trailer)
        constraint2 = SeaPearl.NotEqual(minusY, z, trailer)
        @test SeaPearl.propagate!(constraint2, toPropagate, prunedDomains)

        # Domain not reduced => not propagation
        @test !(constraint in toPropagate)
        @test !(constraint2 in toPropagate)

        # Domain reduced => propagation
        SeaPearl.remove!(z.domain, 6)
        @test SeaPearl.propagate!(constraint2, toPropagate, prunedDomains)
        @test constraint in toPropagate
        @test !(constraint2 in toPropagate)
        @test prunedDomains == SeaPearl.CPModification(
            "-y" => [5],
            "y"  => [-5]
        )
        @test length(y.domain) == 3

        #Unfeasible test
        t = SeaPearl.IntVar(5, 5, "t", trailer)
        constraint3 = SeaPearl.NotEqual(z, t, trailer)
        @test !SeaPearl.propagate!(constraint3, toPropagate, prunedDomains)
    end
end


text:
using Images
using ImageView
using Plots
using Printf
using StatsBase

include("KMeans.jl")
using .KMeans

const MIN_BRIGHTNESS = 0.1
const MAX_BRIGHTNESS = 0.85
const SHOW_PLOT = false
const SHOW_COLORS = false
const NITER = 25
const COLOR_SPACE = "lab"
const TEST_IMG = "/home/kubouch/pictures/kodim/raw/kodim23.png"

function select_base_colors(
    img::Base.ReinterpretArray{T,3,Lab{T}}
) where {T<:Real}
    nch = size(img)[1]

    if size(img)[1] > 3
        error("Error: Wrong number of channels: ", size(img)[1])
    end

    # Select the base color
    base = reshape(mean(img; dims=(2, 3)), nch)

    # Scale the base color between min and max brightness
    # rgb_base = Lab(base[1], base[2], base[3])
    # hsl_base = HSL(rgb_base)

    bases = zeros(nch, 8)

    for (i, l) in enumerate(LinRange(MIN_BRIGHTNESS*100, MAX_BRIGHTNESS*100, 8))
        # res = Lab(HSL(hsl_base.h, hsl_base.s, l))
        res = Lab(l, base[2], base[3])
        bases[1, i] = res.l
        bases[2, i] = res.a
        bases[3, i] = res.b
    end

    return bases
end

function select_base_colors(
    img::Base.ReinterpretArray{T,3,RGB{T}}
) where {T<:Real}
    nch = size(img)[1]

    if size(img)[1] > 3
        error("Error: Wrong number of channels: ", size(img)[1])
    end

    # Select the base color
    base = reshape(mean(img; dims=(2, 3)), nch)

    # Scale the base color between min and max brightness
    rgb_base = RGB(base[1], base[2], base[3])
    hsl_base = HSL(rgb_base)

    bases = zeros(nch, 8)

    for (i, l) in enumerate(LinRange(MIN_BRIGHTNESS, MAX_BRIGHTNESS, 8))
        res = RGB(HSL(hsl_base.h, hsl_base.s, l))
        bases[1, i] = res.r
        bases[2, i] = res.g
        bases[3, i] = res.b
    end

    return bases
end

function show_colors(means_rgb)
    display(colorview(RGB, means_rgb))
end

function show_plot(img, bases, means)
    # Plotly backend opens a browser at the end
    plotly()

    # Plot the image pixels
    # Limit max. number of plotted pixels, otherwise it's too slow
    sz = (min(100, size(img)[2]), min(100, size(img)[2]))
    plot_img = channelview(sample(img, sz))

    scatter(
        plot_img[1, :, :],
        plot_img[2, :, :],
        plot_img[3, :, :];
        # Too slow:
        # marker_z=plot_img[3, :, :],
        markersize=0.5,
        markerstrokewidth=0,
        markeralpha=0.5,
        markercolor="black",
        legend=false,
        size=(1000, 1000),
    )

    # Plot the fixed base colors first
    nbases = if ndims(bases) == 1
        size(bases)[1]
    else
        size(bases)[2]
    end

    scatter!(
        means[1, 1:nbases],
        means[2, 1:nbases],
        means[3, 1:nbases];
        # marker_z=means[3, :, :],
        markersize=2,
        # markerstrokewidth=0,
        markercolor="blue",
        legend=false,
        size=(1000, 1000),
    )

    # Plot the rest of the base colors and show the plot
    gui(
        scatter!(
            means[1, (nbases + 1):end],
            means[2, (nbases + 1):end],
            means[3, (nbases + 1):end];
            # marker_z=means[3, :, :],
            markersize=2,
            # markerstrokewidth=0,
            legend=false,
            size=(1000, 1000),
        ),
    )
end

function print_colors(means_rgb)
    for i in 1:length(means_rgb)
        col = means_rgb[i]

        r = reinterpret(N0f8(col.r))
        g = reinterpret(N0f8(col.g))
        b = reinterpret(N0f8(col.b))

        @printf("#%02x%02x%02x\n", r, g, b)
    end
end

function main()

    prog = Base.basename(Base.source_path())
    img_name = if length(ARGS) == 0
        TEST_IMG
    else
        ARGS[1]
    end

    img = load(img_name)

    color_space = if COLOR_SPACE == "lab"
        Lab
    elseif COLOR_SPACE == "rgb"
        RGB
    else
        error("Unknown color space: ", COLOR_SPACE)
    end

    # Color transform
    img = color_space.(img)

    # Image types other than RGB do not support math ops => need raw pixels
    inp_img = channelview(img)

    # Select base colors
    bases = select_base_colors(inp_img)

    # Calculate k-means
    means = KMeans.kmeans(inp_img, bases, 8, NITER)

    # Convert the result back to RGB
    means_rgb = RGB.(colorview(color_space, means))

    if SHOW_COLORS
        show_colors(means_rgb)
    end

    if SHOW_PLOT
        show_plot(img, bases, means)
    end

    # Print out the final colors
    print_colors(means_rgb)

    return nothing
end

main()


text:
using TheNaNTrap
using Test

@testset "TheNaNTrap.jl" begin
    # Write your tests here.
end


text:
export binhist

"""
    _non0hist(points, binning_scheme::RectangularBinning, dims)

Determine which bins are visited by `points` given the rectangular `binning_scheme`,
considering only the marginal along dimensions `dims`. Bins are referenced
relative to the axis minima.

Returns the unordered (sum-normalized, if `normalize==true`) histogram (visitation
frequency) over the array of bin visits.

# Example
```julia
using using Entropies, DelayEmbeddings
pts = Dataset([rand(5) for i = 1:100]);

# Histograms directly from points given a rectangular binning scheme
h1 = _non0hist(pts, RectangularBinning(0.2), 1:3)
h2 = _non0hist(pts, RectangularBinning(0.2), [1, 2])

# Test that we're actually getting normalised histograms
sum(h1) ≈ 1.0, sum(h2) ≈ 1.0
```
"""
function _non0hist(points, binning_scheme::RectangularBinning, dims)
    bin_visits = marginal_visits(points, binning_scheme, dims)
    _non0hist(bin_visits)
end

function _non0hist(points::AbstractDataset{N, T}, binning_scheme::RectangularBinning) where {N, T}
    _non0hist(points, binning_scheme, 1:N)
end

# TODO: ϵ::RectangularBinning(Float64) allocates slightly more memory than the method below
# because of the call to `minima_edgelengths`. Can be optimized, but keep both versions for
# now and merge docstrings.
probabilities(data::AbstractDataset, ϵ::RectangularBinning) = _non0hist(data, ϵ)[1]
function _non0hist(data::AbstractDataset{D, T}, ϵ::RectangularBinning) where {D, T<:Real}

    # TODO: this allocates a lot, but is not performance critical
    mini, edgelengths = minima_edgelengths(data, ϵ)

    # Map each datapoint to its bin edge and sort the resulting list:
    bins = map(point -> floor.(Int, (point .- mini) ./ edgelengths), data)
    sort!(bins, alg=QuickSort)

    # Reserve enough space for histogram:
    L = length(data)
    hist = Vector{Float64}()
    sizehint!(hist, L)

    # Fill the histogram by counting consecutive equal bins:
    prev_bin, count = bins[1], 0
    for bin in bins
        if bin == prev_bin
            count += 1
        else
            push!(hist, count)
            prev_bin = bin
            count = 1
        end
    end
    push!(hist, count)

    # Shrink histogram capacity to fit its size:
    sizehint!(hist, length(hist))
    return Probabilities(hist ./ L), bins, mini, edgelengths
end

function binhist(x::AbstractDataset{D, T}, ϵ::RectangularBinning) where {D, T<:Real}
    hist, bins, mini, edgelengths = _non0hist(x, ϵ)
    unique!(bins)
    b = [β .* edgelengths .+ mini for β in bins]
    return hist, b
end

# Proper dispatch for when an integer is given.
probabilities(data, n::Integer) =
probabilities(Dataset(data), VisitationFrequency(RectangularBinning(n)))

# The following is originally from ChaosTools.jl
probabilities(data::AbstractDataset, ε::AbstractFloat) = _non0hist(data, ε)[1]
function _non0hist(data::AbstractDataset{D, T}, ε::AbstractFloat) where {D, T<:Real}
    mini = minima(data)
    L = length(data)
    hist = Vector{Float64}()
    # Reserve enough space for histogram:
    sizehint!(hist, L)

    # Map each datapoint to its bin edge and sort the resulting list:
    bins = map(point -> floor.(Int, (point - mini)/ε), data)
    sort!(bins, alg=QuickSort)

    # Fill the histogram by counting consecutive equal bins:
    prev_bin, count = bins[1], 0
    for bin in bins
        if bin == prev_bin
            count += 1
        else
            push!(hist, count)
            prev_bin = bin
            count = 1
        end
    end
    push!(hist, count)

    # Shrink histogram capacity to fit its size:
    sizehint!(hist, length(hist))
    return Probabilities(hist ./ L), bins, mini
end

# TODO: This needs to be expanded to allow `est::RectangularBinning(vector)`
"""
    binhist(x::AbstractDataset, ε::Real) → p, bins
    binhist(x::AbstractDataset, ε::RectangularBinning) → p, bins

Hyper-optimized histogram calculation for `x` with rectangular binning `ε`.
Returns the probabilities `p` of each bin of the histogram as well as the bins.
Notice that `bins` are the starting corners of each bin. If `ε isa Real`, then the actual
bin size is `ε` across each dimension. If `ε isa RectangularBinning`, then the bin
size for each dimension will depend on the binning scheme.

See also: [`RectangularBinning`](@ref).
"""
function binhist(data, ε)
    hist, bins, mini = _non0hist(data, ε)
    unique!(bins)
    b = [β .* ε .+ mini for β in bins]
    return hist, b
end