text:
type DirectJumpAggregation{T,F1,F2} <: AbstractJumpAggregator
  next_jump::T
  end_time::T
  cur_rates::Vector{T}
  sum_rate::T
  rates::F1
  affects!::F2
  save_positions::Tuple{Bool,Bool}
end

@inline function (p::DirectJumpAggregation)(u,t,integrator) # condition
  p.next_jump == t
end

function (p::DirectJumpAggregation)(integrator) # affect!
  rng_val         = rand() * p.sum_rate
  i               = searchsortedfirst(p.cur_rates, rng_val)
  @inbounds p.affects![i](integrator)
  p.sum_rate,ttnj = time_to_next_jump(integrator.u,integrator.p,integrator.t,p.rates,p.cur_rates)
  p.next_jump     = integrator.t + ttnj
  if p.next_jump < p.end_time
    add_tstop!(integrator,p.next_jump)
  end
  nothing
end

function (p::DirectJumpAggregation)(dj,u,t,integrator) # initialize
  sum_rate,next_jump = time_to_next_jump(u,integrator.p,t,p.rates,p.cur_rates)
  p.sum_rate = sum_rate
  p.next_jump = t + next_jump
  if p.next_jump < p.end_time
    add_tstop!(integrator,p.next_jump)
  end
  nothing
end

################# Tuple based time to next jump ##################
function time_to_next_jump(u, p, t, rates::Tuple, cur_rates)
  @inbounds fill_cur_rates(u,p,t,cur_rates,1,rates...)
  @inbounds for i in 2:length(cur_rates) # cumsum
    cur_rates[i] = cur_rates[i] + cur_rates[i-1]
  end
  @inbounds sum_rate = cur_rates[end]
  sum_rate,randexp()/sum_rate
end

@inline function fill_cur_rates(u,p,t,cur_rates,idx,rate,rates...)
  @inbounds cur_rates[idx] = rate(u,p,t)
  idx += 1
  fill_cur_rates(u,p,t,cur_rates,idx,rates...)
end

@inline function fill_cur_rates(u,p,t,cur_rates,idx,rate)
  @inbounds cur_rates[idx] = rate(u,p,t)
  nothing
end

############ FunctionWrapper based time to next jump #############
function time_to_next_jump(u, p, t, rates::Vector{T}, cur_rates) where T
  prev_rate = zero(t)
  new_rate  = zero(t)
  @inbounds for i in 1:length(rates)
      new_rate     = rates[i](u, p, t)
      cur_rates[i] = new_rate + prev_rate
      prev_rate    = cur_rates[i]
  end
  @inbounds sum_rate = cur_rates[end]
  sum_rate,randexp()/sum_rate
end

@inline function aggregate(aggregator::Direct,u,p,t,end_time,constant_jumps,save_positions)

  # decide if representing rates/affects as tuples or function wrappers
  if length(constant_jumps) < TUPLE_TO_FWRAPPER_CUTOFF
    rates,affects! = get_jump_info_tuples(constant_jumps)
  else
    rates,affects! = get_jump_info_fwrappers(u,p,t,constant_jumps)
  end

  cur_rates = zeros(typeof(t),length(rates))
  sum_rate  = zero(T)
  next_jump = zero(T)
  DirectJumpAggregation(next_jump, end_time, cur_rates, sum_rate, rates, affects!, save_positions)
end

DiscreteCallback(c::DirectJumpAggregation) = DiscreteCallback(c,c,initialize=c,save_positions=c.save_positions)

max_stars_repo_name: isaacsas/DiffEqJumpExtensions

max_issues_repo_name: isaacsas/DiffEqJumpExtensions

max_forks_repo_name: isaacsas/DiffEqJumpExtensions


text:
## @lc code=start
using LeetCode

function count_operations(num1::Int, num2::Int)::Int
    if num2 > num1
        num1, num2 = num2, num1
    end
    res = 0
    while num1 > 0
        res += num2 ÷ num1
        num1, num2 = num2 % num1, num1
    end
    res
end
## @lc code=end

## @lc test=start
@testset "2169.count-operations-to-obtain-zero.jl" begin
    @test count_operations(2, 3) == 3
    @test count_operations(10, 10) == 1
    @test count_operations(2022, 1077) == 20
end
## @lc test=end

## @lc info=start
# link: [solution to 2169](https://leetcode-cn.com/problems/count-operations-to-obtain-zero/solution/6004-zhan-zhuan-xiang-chu-de-ci-shu-juli-23oc/)
## @lc info=end

ext: jl


text:
include("../src/MDDF.jl")

# The solute is a single protein molecule (infinte dilution case). In this case,
# use the option nmols=1
solute = MDDF.Selection( MDDF.VMDselect("structure.pdb","protein",vmd="/usr/local/bin/vmd"),
                      nmols=1 )

# The solvent is TMAO, which has 14 atoms. Use the natomspermol to indicate how many
# atoms each molecule has, such that there is no ambiguity on how to split the coordinates
# of the selection into individual molecules.
solvent = MDDF.Selection( MDDF.VMDselect("structure.pdb","resname TMAO",vmd="/usr/local/bin/vmd"),
                        natomspermol=14 )

# Initialize trajectroy data structure and open input stream
trajectory = MDDF.PDBTraj("./trajectory.pdb",solute,solvent)

# Run MDDF calculation
#mddf = MDDF.mddf_naive(solute,solvent,


text:
using OnlinePCA
using OnlinePCA: readcsv, writecsv
using Test
using Pkg
using DelimitedFiles
using Statistics
using Distributions

tmp = mktempdir()
println(tmp)
julia = joinpath(Sys.BINDIR, "julia")
bindir = joinpath(dirname(pathof(OnlinePCA)), "..", "bin")

function testfilesize(remove::Bool, x...)
	for n = 1:length(x)
		@test filesize(x[n]) != 0
		if remove
			rm(x[n])
		end
	end
end

#####################################
input = Int64.(ceil.(rand(NegativeBinomial(1, 0.5), 300, 99)))
input[1:50, 1:33] .= 100*input[1:50, 1:33]
input[51:100, 34:66] .= 100*input[51:100, 34:66]
input[101:150, 67:99] .= 100*input[101:150, 67:99]
writecsv(joinpath(tmp, "Data.csv"), input)
#####################################

#####################################
println("####### Binarization (Julia API) #######")
csv2bin(csvfile=joinpath(tmp, "Data.csv"),
	binfile=joinpath(tmp, "Data.zst"))

testfilesize(true, joinpath(tmp, "Data.zst"))
###################################

#####################################
println("####### Binarization (Command line) #######")
run(`$(julia) $(joinpath(bindir, "csv2bin")) --csvfile $(joinpath(tmp, "Data.csv")) --binfile $(joinpath(tmp, "Data.zst"))`)

testfilesize(false, joinpath(tmp, "Data.zst"))
#####################################

#####################################
println("####### Summarization (Julia API) #######")
sumr(binfile=joinpath(tmp, "Data.zst"), outdir=tmp)

testfilesize(true,
	joinpath(tmp, "Sample_NoCounts.csv"),
	joinpath(tmp, "Feature_Means.csv"),
	joinpath(tmp, "Feature_LogMeans.csv"),
	joinpath(tmp, "Feature_FTTMeans.csv"),
	joinpath(tmp, "Feature_CPMMeans.csv"),
	joinpath(tmp, "Feature_LogCPMMeans.csv"),
	joinpath(tmp, "Feature_FTTCPMMeans.csv"),
	joinpath(tmp, "Feature_CPTMeans.csv"),
	joinpath(tmp, "Feature_LogCPTMeans.csv"),
	joinpath(tmp, "Feature_FTTCPTMeans.csv"),
	joinpath(tmp, "Feature_CPMEDMeans.csv"),
	joinpath(tmp, "Feature_LogCPMEDMeans.csv"),
	joinpath(tmp, "Feature_FTTCPMEDMeans.csv"),
	joinpath(tmp, "Feature_Vars.csv"),
	joinpath(tmp, "Feature_LogVars.csv"),
	joinpath(tmp, "Feature_FTTVars.csv"),
	joinpath(tmp, "Feature_CPMVars.csv"),
	joinpath(tmp, "Feature_LogCPMVars.csv"),
	joinpath(tmp, "Feature_FTTCPMVars.csv"),
	joinpath(tmp, "Feature_CPTVars.csv"),
	joinpath(tmp, "Feature_LogCPTVars.csv"),
	joinpath(tmp, "Feature_FTTCPTVars.csv"),
	joinpath(tmp, "Feature_CPMEDVars.csv"),
	joinpath(tmp, "Feature_LogCPMEDVars.csv"),
	joinpath(tmp, "Feature_FTTCPMEDVars.csv"),
	joinpath(tmp, "Feature_CV2s.csv"),
	joinpath(tmp, "Feature_NoZeros.csv"))
###################################

#####################################
println("####### Summarization (Command line) #######")
run(`$(julia) $(joinpath(bindir, "sumr")) --binfile $(joinpath(tmp, "Data.zst")) --outdir $(tmp)`)

testfilesize(false,
	joinpath(tmp, "Sample_NoCounts.csv"),
	joinpath(tmp, "Feature_Means.csv"),
	joinpath(tmp, "Feature_LogMeans.csv"),
	joinpath(tmp, "Feature_FTTMeans.csv"),
	joinpath(tmp, "Feature_CPMMeans.csv"),
	joinpath(tmp, "Feature_LogCPMMeans.csv"),
	joinpath(tmp, "Feature_FTTCPMMeans.csv"),
	joinpath(tmp, "Feature_CPTMeans.csv"),
	joinpath(tmp, "Feature_LogCPTMeans.csv"),
	joinpath(tmp, "Feature_FTTCPTMeans.csv"),
	joinpath(tmp, "Feature_CPMEDMeans.csv"),
	joinpath(tmp, "Feature_LogCPMEDMeans.csv"),
	joinpath(tmp, "Feature_FTTCPMEDMeans.csv"),
	joinpath(tmp, "Feature_Vars.csv"),
	joinpath(tmp, "Feature_LogVars.csv"),
	joinpath(tmp, "Feature_FTTVars.csv"),
	joinpath(tmp, "Feature_CPMVars.csv"),
	joinpath(tmp, "Feature_LogCPMVars.csv"),
	joinpath(tmp, "Feature_FTTCPMVars.csv"),
	joinpath(tmp, "Feature_CPTVars.csv"),
	joinpath(tmp, "Feature_LogCPTVars.csv"),
	joinpath(tmp, "Feature_FTTCPTVars.csv"),
	joinpath(tmp, "Feature_CPMEDVars.csv"),
	joinpath(tmp, "Feature_LogCPMEDVars.csv"),
	joinpath(tmp, "Feature_FTTCPMEDVars.csv"),
	joinpath(tmp, "Feature_CV2s.csv"),
	joinpath(tmp, "Feature_NoZeros.csv"))
#####################################

#####################################
println("####### HVG (Julia API) #######")
hvg(binfile=joinpath(tmp, "Data.zst"),
	rowmeanlist=joinpath(tmp, "Feature_Means.csv"),
	rowvarlist=joinpath(tmp, "Feature_Vars.csv"),
	rowcv2list=joinpath(tmp, "Feature_CV2s.csv"), outdir=tmp)

testfilesize(true,
	joinpath(tmp, "HVG_pval.csv"),
	joinpath(tmp, "HVG_a0.csv"),
	joinpath(tmp, "HVG_a1.csv"),
	joinpath(tmp, "HVG_afit.csv"),
	joinpath(tmp, "HVG_useForFit.csv"),
	joinpath(tmp, "HVG_varFitRatio.csv"),
	joinpath(tmp, "HVG_df.csv"))
#####################################

#####################################
println("####### HVG (Command line) #######")
run(`$(julia) $(joinpath(bindir, "hvg")) --binfile $(joinpath(tmp, "Data.zst")) --rowmeanlist $(joinpath(tmp, "Feature_Means.csv")) --rowvarlist $(joinpath(tmp, "Feature_Vars.csv")) --rowcv2list $(joinpath(tmp, "Feature_CV2s.csv")) --outdir $(tmp)`)

testfilesize(false,
	joinpath(tmp, "HVG_pval.csv"),
	joinpath(tmp, "HVG_a0.csv"),
	joinpath(tmp, "HVG_a1.csv"),
	joinpath(tmp, "HVG_afit.csv"),
	joinpath(tmp, "HVG_useForFit.csv"),
	joinpath(tmp, "HVG_varFitRatio.csv"),
	joinpath(tmp, "HVG_df.csv"))
#####################################

#####################################
println("####### Filtering (Julia API) #######")
filtering(input=joinpath(tmp, "Data.zst"),
	featurelist=joinpath(tmp, "Feature_Means.csv"),
	samplelist=joinpath(tmp, "Sample_NoCounts.csv"),
	thr1=10, thr2=10,
	direct1="+", direct2="+",
	outdir=tmp)

filtering(input=joinpath(tmp, "Data.zst"),
	featurelist=joinpath(tmp, "Feature_Means.csv"),
	thr1=10, thr2=10,
	direct1="+", direct2="+",
	outdir=tmp)

filtering(input=joinpath(tmp, "Data.zst"),
	samplelist=joinpath(tmp, "Sample_NoCounts.csv"),
	thr1=10, thr2=10,
	direct1="+", direct2="+",
	outdir=tmp)

testfilesize(true,
	joinpath(tmp, "filtered.zst"))
testfilesize(true,
	joinpath(tmp, "filteredFeature.csv"))
testfilesize(true,
	joinpath(tmp, "filteredSample.csv"))
####################################

#####################################
println("####### Filtering (Command line) #######")
run(`$(julia) $(joinpath(bindir, "filtering")) --input $(joinpath(tmp, "Data.zst")) --featurelist $(joinpath(tmp, "Feature_Means.csv")) --samplelist $(joinpath(tmp, "Sample_NoCounts.csv")) --thr1 10 --thr2 10 --outdir $(tmp)`)

run(`$(julia) $(joinpath(bindir, "filtering")) --input $(joinpath(tmp, "Data.zst")) --featurelist $(joinpath(tmp, "Feature_Means.csv")) --thr1 10 --thr2 10 --outdir $(tmp)`)

run(`$(julia) $(joinpath(bindir, "filtering")) --input $(joinpath(tmp, "Data.zst")) --samplelist $(joinpath(tmp, "Sample_NoCounts.csv")) --thr1 10 --thr2 10 --outdir $(tmp)`)

testfilesize(false,
	joinpath(tmp, "filtered.zst"))
testfilesize(false,
	joinpath(tmp, "filteredFeature.csv"))
testfilesize(false,
	joinpath(tmp, "filteredSample.csv"))
#####################################

#####################################
println("####### GD (Julia API) #######")
out_gd1 = gd(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="robbins-monro",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

out_gd2 = gd(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="momentum",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

out_gd3 = gd(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="nag",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

out_gd4 = gd(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="adagrad",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

@test size(out_gd1[1]) == (99, 3)
@test size(out_gd1[2]) == (3, )
@test size(out_gd1[3]) == (300, 3)
@test size(out_gd1[4]) == (99, 3)
@test size(out_gd1[5]) == ()
@test size(out_gd1[6]) == ()

@test size(out_gd2[1]) == (99, 3)
@test size(out_gd2[2]) == (3, )
@test size(out_gd2[3]) == (300, 3)
@test size(out_gd2[4]) == (99, 3)
@test size(out_gd2[5]) == ()
@test size(out_gd2[6]) == ()

@test size(out_gd3[1]) == (99, 3)
@test size(out_gd3[2]) == (3, )
@test size(out_gd3[3]) == (300, 3)
@test size(out_gd3[4]) == (99, 3)
@test size(out_gd3[5]) == ()
@test size(out_gd3[6]) == ()

@test size(out_gd4[1]) == (99, 3)
@test size(out_gd4[2]) == (3, )
@test size(out_gd4[3]) == (300, 3)
@test size(out_gd4[4]) == (99, 3)
@test size(out_gd4[5]) == ()
@test size(out_gd4[6]) == ()
#####################################

#####################################
println("####### GD (Command line) #######")
run(`$(julia) $(joinpath(bindir, "gd")) --input $(joinpath(tmp, "Data.zst")) --outdir $(tmp) --dim 3 --scheduling robbins-monro --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(false,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))

run(`$(julia) $(joinpath(bindir, "gd")) --initW $(joinpath(tmp, "Eigen_vectors.csv")) --input $(joinpath(tmp, "Data.zst"))  --outdir $(tmp) --dim 3 --scheduling momentum --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(false,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))

run(`$(julia) $(joinpath(bindir, "gd")) --initV $(joinpath(tmp, "Loadings.csv")) --input $(joinpath(tmp, "Data.zst"))  --outdir $(tmp) --dim 3 --scheduling nag --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))

run(`$(julia) $(joinpath(bindir, "gd")) --input $(joinpath(tmp, "Data.zst"))  --outdir $(tmp) --dim 3 --scheduling adagrad --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))
#####################################

#####################################
println("####### Oja (Julia API) #######")
out_oja1 = oja(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="robbins-monro",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

out_oja2 = oja(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="momentum",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

out_oja3 = oja(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="nag",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

out_oja4 = oja(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="adagrad",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

@test size(out_oja1[1]) == (99, 3)
@test size(out_oja1[2]) == (3, )
@test size(out_oja1[3]) == (300, 3)
@test size(out_oja1[4]) == (99, 3)
@test size(out_oja1[5]) == ()
@test size(out_oja1[6]) == ()

@test size(out_oja2[1]) == (99, 3)
@test size(out_oja2[2]) == (3, )
@test size(out_oja2[3]) == (300, 3)
@test size(out_oja2[4]) == (99, 3)
@test size(out_oja2[5]) == ()
@test size(out_oja2[6]) == ()

@test size(out_oja3[1]) == (99, 3)
@test size(out_oja3[2]) == (3, )
@test size(out_oja3[3]) == (300, 3)
@test size(out_oja3[4]) == (99, 3)
@test size(out_oja3[5]) == ()
@test size(out_oja3[6]) == ()

@test size(out_oja4[1]) == (99, 3)
@test size(out_oja4[2]) == (3, )
@test size(out_oja4[3]) == (300, 3)
@test size(out_oja4[4]) == (99, 3)
@test size(out_oja4[5]) == ()
@test size(out_oja4[6]) == ()
#####################################

#####################################
println("####### Oja (Command line) #######")
run(`$(julia) $(joinpath(bindir, "oja")) --input $(joinpath(tmp, "Data.zst")) --outdir $(tmp) --dim 3 --scheduling robbins-monro --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))

run(`$(julia) $(joinpath(bindir, "oja")) --input $(joinpath(tmp, "Data.zst"))  --outdir $(tmp) --dim 3 --scheduling momentum --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))

run(`$(julia) $(joinpath(bindir, "oja")) --input $(joinpath(tmp, "Data.zst"))  --outdir $(tmp) --dim 3 --scheduling nag --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))

run(`$(julia) $(joinpath(bindir, "oja")) --input $(joinpath(tmp, "Data.zst"))  --outdir $(tmp) --dim 3 --scheduling adagrad --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))
#####################################

#####################################
println("####### SGD (Julia API) #######")
out_sgd1 = sgd(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="robbins-monro",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

out_sgd2 = sgd(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="momentum",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

out_sgd3 = sgd(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="nag",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

out_sgd4 = sgd(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="adagrad",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

@test size(out_sgd1[1]) == (99, 3)
@test size(out_sgd1[2]) == (3, )
@test size(out_sgd1[3]) == (300, 3)
@test size(out_sgd1[4]) == (99, 3)
@test size(out_sgd1[5]) == ()
@test size(out_sgd1[6]) == ()

@test size(out_sgd2[1]) == (99, 3)
@test size(out_sgd2[2]) == (3, )
@test size(out_sgd2[3]) == (300, 3)
@test size(out_sgd2[4]) == (99, 3)
@test size(out_sgd2[5]) == ()
@test size(out_sgd2[6]) == ()

@test size(out_sgd3[1]) == (99, 3)
@test size(out_sgd3[2]) == (3, )
@test size(out_sgd3[3]) == (300, 3)
@test size(out_sgd3[4]) == (99, 3)
@test size(out_sgd3[5]) == ()
@test size(out_sgd3[6]) == ()

@test size(out_sgd4[1]) == (99, 3)
@test size(out_sgd4[2]) == (3, )
@test size(out_sgd4[3]) == (300, 3)
@test size(out_sgd4[4]) == (99, 3)
@test size(out_sgd4[5]) == ()
@test size(out_sgd4[6]) == ()
#####################################

#####################################
println("####### SGD (Command line) #######")
run(`$(julia) $(joinpath(bindir, "sgd")) --input $(joinpath(tmp, "Data.zst")) --outdir $(tmp) --dim 3 --scheduling robbins-monro --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))

run(`$(julia) $(joinpath(bindir, "sgd")) --input $(joinpath(tmp, "Data.zst"))  --outdir $(tmp) --dim 3 --scheduling momentum --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))

run(`$(julia) $(joinpath(bindir, "sgd")) --input $(joinpath(tmp, "Data.zst"))  --outdir $(tmp) --dim 3 --scheduling nag --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))

run(`$(julia) $(joinpath(bindir, "sgd")) --input $(joinpath(tmp, "Data.zst"))  --outdir $(tmp) --dim 3 --scheduling adagrad --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))
#####################################

#####################################
println("####### CCIPCA (Julia API) #######")
out_ccipca1 = ccipca(input=joinpath(tmp, "Data.zst"),
	dim=3, stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

@test size(out_ccipca1[1]) == (99, 3)
@test size(out_ccipca1[2]) == (3, )
@test size(out_ccipca1[3]) == (300, 3)
@test size(out_ccipca1[4]) == (99, 3)
@test size(out_ccipca1[5]) == ()
@test size(out_ccipca1[6]) == ()
#####################################

#####################################
println("####### CCIPCA (Command line) #######")
run(`$(julia) $(joinpath(bindir, "ccipca")) --input $(joinpath(tmp, "Data.zst")) --outdir $(tmp) --dim 3 --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))
#####################################

#####################################
println("####### RSGD (Julia API) #######")
out_rsgd1 = rsgd(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="robbins-monro",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

out_rsgd2 = rsgd(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="momentum",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

out_rsgd3 = rsgd(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="nag",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

out_rsgd4 = rsgd(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="adagrad",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

@test size(out_rsgd1[1]) == (99, 3)
@test size(out_rsgd1[2]) == (3, )
@test size(out_rsgd1[3]) == (300, 3)
@test size(out_rsgd1[4]) == (99, 3)
@test size(out_rsgd1[5]) == ()
@test size(out_rsgd1[6]) == ()

@test size(out_rsgd2[1]) == (99, 3)
@test size(out_rsgd2[2]) == (3, )
@test size(out_rsgd2[3]) == (300, 3)
@test size(out_rsgd2[4]) == (99, 3)
@test size(out_rsgd2[5]) == ()
@test size(out_rsgd2[6]) == ()

@test size(out_rsgd3[1]) == (99, 3)
@test size(out_rsgd3[2]) == (3, )
@test size(out_rsgd3[3]) == (300, 3)
@test size(out_rsgd3[4]) == (99, 3)
@test size(out_rsgd3[5]) == ()
@test size(out_rsgd3[6]) == ()

@test size(out_rsgd4[1]) == (99, 3)
@test size(out_rsgd4[2]) == (3, )
@test size(out_rsgd4[3]) == (300, 3)
@test size(out_rsgd4[4]) == (99, 3)
@test size(out_rsgd4[5]) == ()
@test size(out_rsgd4[6]) == ()
#####################################

#####################################
println("####### RSGD (Command line) #######")
run(`$(julia) $(joinpath(bindir, "rsgd")) --input $(joinpath(tmp, "Data.zst")) --outdir $(tmp) --dim 3 --scheduling robbins-monro --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))

run(`$(julia) $(joinpath(bindir, "rsgd")) --input $(joinpath(tmp, "Data.zst"))  --outdir $(tmp) --dim 3 --scheduling momentum --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))

run(`$(julia) $(joinpath(bindir, "rsgd")) --input $(joinpath(tmp, "Data.zst"))  --outdir $(tmp) --dim 3 --scheduling nag --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))

run(`$(julia) $(joinpath(bindir, "rsgd")) --input $(joinpath(tmp, "Data.zst"))  --outdir $(tmp) --dim 3 --scheduling adagrad --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))
#####################################

#####################################
println("####### SVRG (Julia API) #######")
out_svrg1 = svrg(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="robbins-monro",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

out_svrg2 = svrg(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="momentum",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

out_svrg3 = svrg(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="nag",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

out_svrg4 = svrg(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="adagrad",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

@test size(out_svrg1[1]) == (99, 3)
@test size(out_svrg1[2]) == (3, )
@test size(out_svrg1[3]) == (300, 3)
@test size(out_svrg1[4]) == (99, 3)
@test size(out_svrg1[5]) == ()
@test size(out_svrg1[6]) == ()

@test size(out_svrg2[1]) == (99, 3)
@test size(out_svrg2[2]) == (3, )
@test size(out_svrg2[3]) == (300, 3)
@test size(out_svrg2[4]) == (99, 3)
@test size(out_svrg2[5]) == ()
@test size(out_svrg2[6]) == ()

@test size(out_svrg3[1]) == (99, 3)
@test size(out_svrg3[2]) == (3, )
@test size(out_svrg3[3]) == (300, 3)
@test size(out_svrg3[4]) == (99, 3)
@test size(out_svrg3[5]) == ()
@test size(out_svrg3[6]) == ()

@test size(out_svrg4[1]) == (99, 3)
@test size(out_svrg4[2]) == (3, )
@test size(out_svrg4[3]) == (300, 3)
@test size(out_svrg4[4]) == (99, 3)
@test size(out_svrg4[5]) == ()
@test size(out_svrg4[6]) == ()
#####################################

#####################################
println("####### SVRG (Command line) #######")
run(`$(julia) $(joinpath(bindir, "svrg")) --input $(joinpath(tmp, "Data.zst")) --outdir $(tmp) --dim 3 --scheduling robbins-monro --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))

run(`$(julia) $(joinpath(bindir, "svrg")) --input $(joinpath(tmp, "Data.zst"))  --outdir $(tmp) --dim 3 --scheduling momentum --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))

run(`$(julia) $(joinpath(bindir, "svrg")) --input $(joinpath(tmp, "Data.zst"))  --outdir $(tmp) --dim 3 --scheduling nag --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))

run(`$(julia) $(joinpath(bindir, "svrg")) --input $(joinpath(tmp, "Data.zst"))  --outdir $(tmp) --dim 3 --scheduling adagrad --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))
#####################################

#####################################
println("####### RSVRG (Julia API) #######")
out_rsvrg1 = rsvrg(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="robbins-monro",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

out_rsvrg2 = rsvrg(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="momentum",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

out_rsvrg3 = rsvrg(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="nag",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

out_rsvrg4 = rsvrg(input=joinpath(tmp, "Data.zst"),
	dim=3, scheduling="adagrad",
	stepsize=1.0e-15, numepoch=1,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

@test size(out_rsvrg1[1]) == (99, 3)
@test size(out_rsvrg1[2]) == (3, )
@test size(out_rsvrg1[3]) == (300, 3)
@test size(out_rsvrg1[4]) == (99, 3)
@test size(out_rsvrg1[5]) == ()
@test size(out_rsvrg1[6]) == ()

@test size(out_rsvrg2[1]) == (99, 3)
@test size(out_rsvrg2[2]) == (3, )
@test size(out_rsvrg2[3]) == (300, 3)
@test size(out_rsvrg2[4]) == (99, 3)
@test size(out_rsvrg2[5]) == ()
@test size(out_rsvrg2[6]) == ()

@test size(out_rsvrg3[1]) == (99, 3)
@test size(out_rsvrg3[2]) == (3, )
@test size(out_rsvrg3[3]) == (300, 3)
@test size(out_rsvrg3[4]) == (99, 3)
@test size(out_rsvrg3[5]) == ()
@test size(out_rsvrg3[6]) == ()

@test size(out_rsvrg4[1]) == (99, 3)
@test size(out_rsvrg4[2]) == (3, )
@test size(out_rsvrg4[3]) == (300, 3)
@test size(out_rsvrg4[4]) == (99, 3)
@test size(out_rsvrg4[5]) == ()
@test size(out_rsvrg4[6]) == ()
#####################################

#####################################
println("####### RSVRG (Command line) #######")
run(`$(julia) $(joinpath(bindir, "rsvrg")) --input $(joinpath(tmp, "Data.zst")) --outdir $(tmp) --dim 3 --scheduling robbins-monro --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))

run(`$(julia) $(joinpath(bindir, "rsvrg")) --input $(joinpath(tmp, "Data.zst"))  --outdir $(tmp) --dim 3 --scheduling momentum --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))

run(`$(julia) $(joinpath(bindir, "rsvrg")) --input $(joinpath(tmp, "Data.zst"))  --outdir $(tmp) --dim 3 --scheduling nag --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))

run(`$(julia) $(joinpath(bindir, "rsvrg")) --input $(joinpath(tmp, "Data.zst"))  --outdir $(tmp) --dim 3 --scheduling adagrad --stepsize 1.0e-15 --numepoch 1 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))
#####################################

#####################################
println("####### Orthogonal Iteration (Julia API) #######")
out_orthiter1 = orthiter(input=joinpath(tmp, "Data.zst"),
	dim=3, numepoch=3,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

@test size(out_orthiter1[1]) == (99, 3)
@test size(out_orthiter1[2]) == (3, )
@test size(out_orthiter1[3]) == (300, 3)
@test size(out_orthiter1[4]) == (99, 3)
@test size(out_orthiter1[5]) == ()
@test size(out_orthiter1[6]) == ()
#####################################

#####################################
println("####### Orthogonal Iteration (Command line) #######")
run(`$(julia) $(joinpath(bindir, "orthiter")) --input $(joinpath(tmp, "Data.zst")) --outdir $(tmp) --dim 3 --numepoch 3 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))
#####################################

#####################################
println("####### Randomized Block Krylov Iteration (Julia API) #######")
out_rbkiter1 = rbkiter(input=joinpath(tmp, "Data.zst"),
	dim=3, numepoch=3,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

@test size(out_rbkiter1[1]) == (99, 3)
@test size(out_rbkiter1[2]) == (3, )
@test size(out_rbkiter1[3]) == (300, 3)
@test size(out_rbkiter1[4]) == (99, 3)
@test size(out_rbkiter1[5]) == ()
@test size(out_rbkiter1[6]) == ()
#####################################

#####################################
println("####### Randomized Block Krylov Iteration (Command line) #######")
run(`$(julia) $(joinpath(bindir, "rbkiter")) --input $(joinpath(tmp, "Data.zst")) --outdir $(tmp) --dim 3 --numepoch 3 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))
#####################################

#####################################
println("####### ARNOLDI (Julia API) #######")
out_arnoldi1 = arnoldi(input=joinpath(tmp, "Data.zst"),
	dim=3,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"))

@test size(out_arnoldi1[1]) == (99, 3)
@test size(out_arnoldi1[2]) == (3, )
@test size(out_arnoldi1[3]) == (300, 3)
@test size(out_arnoldi1[4]) == (99, 3)
@test size(out_arnoldi1[5]) == ()
@test size(out_arnoldi1[6]) == ()
#####################################

#####################################
println("####### ARNOLDI (Command line) #######")
run(`$(julia) $(joinpath(bindir, "arnoldi")) --input $(joinpath(tmp, "Data.zst")) --outdir $(tmp) --dim 3 --numepoch 3 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv"))`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))
#####################################

#####################################
println("####### LANCZOS (Julia API) #######")
out_lanczos1 = lanczos(input=joinpath(tmp, "Data.zst"),
	dim=3,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"))

@test size(out_lanczos1[1]) == (99, 3)
@test size(out_lanczos1[2]) == (3, )
@test size(out_lanczos1[3]) == (300, 3)
@test size(out_lanczos1[4]) == (99, 3)
@test size(out_lanczos1[5]) == ()
@test size(out_lanczos1[6]) == ()
#####################################

#####################################
println("####### LANCZOS (Command line) #######")
run(`$(julia) $(joinpath(bindir, "lanczos")) --input $(joinpath(tmp, "Data.zst")) --outdir $(tmp) --dim 3 --numepoch 3 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv"))`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))
#####################################

#####################################
println("####### HALKO (Julia API) #######")
out_halko = halko(input=joinpath(tmp, "Data.zst"),
	dim=3,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

@test size(out_halko[1]) == (99, 3)
@test size(out_halko[2]) == (3, )
@test size(out_halko[3]) == (300, 3)
@test size(out_halko[4]) == (99, 3)
@test size(out_halko[5]) == ()
####################################

#####################################
println("####### HALKO (Command line) #######")
run(`$(julia) $(joinpath(bindir, "halko")) --input $(joinpath(tmp, "Data.zst")) --outdir $(tmp) --dim 3 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))
#####################################

#####################################
println("####### ALGORITHM971 (Julia API) #######")
out_algorithm971 = algorithm971(input=joinpath(tmp, "Data.zst"),
	dim=3,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

@test size(out_algorithm971[1]) == (99, 3)
@test size(out_algorithm971[2]) == (3, )
@test size(out_algorithm971[3]) == (300, 3)
@test size(out_algorithm971[4]) == (99, 3)
@test size(out_algorithm971[5]) == ()
#####################################

#####################################
println("####### ALGORITHM971 (Command line) #######")
run(`$(julia) $(joinpath(bindir, "algorithm971")) --input $(joinpath(tmp, "Data.zst")) --outdir $(tmp) --dim 3 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))
#####################################

#####################################
println("####### SINGLEPASS (Julia API) #######")
out_singlepass = singlepass(input=joinpath(tmp, "Data.zst"),
	dim=3,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

@test size(out_singlepass[1]) == (99, 3)
@test size(out_singlepass[2]) == (3, )
@test size(out_singlepass[3]) == (300, 3)
@test size(out_singlepass[4]) == (99, 3)
@test size(out_singlepass[5]) == ()
####################################

#####################################
println("####### SINGLEPASS (Command line) #######")
run(`$(julia) $(joinpath(bindir, "singlepass")) --input $(joinpath(tmp, "Data.zst")) --outdir $(tmp) --dim 3 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))
#####################################

#####################################
println("####### SINGLEPASS2 (Julia API) #######")
out_singlepass2 = singlepass2(input=joinpath(tmp, "Data.zst"),
	dim=3,
	rowmeanlist=joinpath(tmp, "Feature_FTTMeans.csv"),
	logdir=tmp)

@test size(out_singlepass2[1]) == (99, 3)
@test size(out_singlepass2[2]) == (3, )
@test size(out_singlepass2[3]) == (300, 3)
@test size(out_singlepass2[4]) == (99, 3)
@test size(out_singlepass2[5]) == ()
####################################

#####################################
println("####### SINGLEPASS2 (Command line) #######")
run(`$(julia) $(joinpath(bindir, "singlepass2")) --input $(joinpath(tmp, "Data.zst")) --outdir $(tmp) --dim 3 --rowmeanlist $(joinpath(tmp, "Feature_FTTMeans.csv")) --logdir $(tmp)`)

testfilesize(true,
	joinpath(tmp, "Eigen_vectors.csv"),
	joinpath(tmp, "Eigen_values.csv"),
	joinpath(tmp, "Loadings.csv"),
	joinpath(tmp, "Scores.csv"))
#####################################


text:
@testset "setters" begin
    @testset "Mag" begin
        @test π <
              Float64(Arblib.set!(Mag(), π)) ==
              Float64(Arblib.set!(Arblib.radref(Arb()), π)) ==
              Float64(Arblib.set!(Arblib.mag_struct(), π)) <
              3.15

        @test Arblib.set!(Mag(), 1) == Arblib.set!(Mag(), UInt(1))
        @test Arblib.set!(Mag(), 3, 4) == Arblib.set!(Mag(), 3 * 2^4)
    end

    @testset "Arb" begin
        # MagLike and BigFloat
        @test Arblib.set!(Arb(), BigFloat(1.0)) == one(Arb)
        @test Arblib.set!(Arblib.realref(Acb()), BigFloat(1.0)) == one(Arb)
        @test Arblib.equal(Arblib.set!(Arblib.arb_struct(), BigFloat(1.0)), one(Arb))
        @test Arblib.set!(Arb(), one(Mag)) == one(Arb)
        @test Arblib.set!(Arblib.realref(Acb()), one(Mag)) == one(Arb)
        @test Arblib.equal(Arblib.set!(Arblib.arb_struct(), one(Mag).mag), one(Arb))
        # Check that aliasing works
        x = Arb()
        Arblib.set!(Arblib.radref(x), 1)
        @test Arblib.set!(x, Arblib.radref(x)) == one(Arb)

        # Rational
        @test Arblib.set!(Arb(), 1 // 2) == one(Arb) / 2
        @test Arblib.set!(Arblib.realref(Acb()), 1 // 2) == one(Arb) / 2
        @test Arblib.equal(Arblib.set!(Arblib.arb_struct(), 1 // 2), one(Arb) / 2)

        # Irrationals
        @test Arb(3) < Arblib.set!(Arb(), π) < Arb(4)
        @test Arblib.equal(Arblib.set!(Arb(), π), Arblib.set!(Arblib.realref(Acb()), π))
        @test Arblib.equal(Arblib.set!(Arb(), π), Arblib.set!(Arblib.arb_struct(), π))

        @test Arb(2) < Arblib.set!(Arb(), ℯ) < Arb(3)
        @test Arblib.equal(Arblib.set!(Arb(), ℯ), Arblib.set!(Arblib.realref(Acb()), ℯ))
        @test Arblib.equal(Arblib.set!(Arb(), ℯ), Arblib.set!(Arblib.arb_struct(), ℯ))

        @test Arb(0) < Arblib.set!(Arb(), MathConstants.γ) < Arb(1)
        @test Arblib.equal(
            Arblib.set!(Arb(), MathConstants.γ),
            Arblib.set!(Arblib.realref(Acb()), MathConstants.γ),
        )
        @test Arblib.equal(
            Arblib.set!(Arb(), MathConstants.γ),
            Arblib.set!(Arblib.arb_struct(), MathConstants.γ),
        )

        @test Arb(0) < Arblib.set!(Arb(), MathConstants.catalan) < Arb(1)
        @test Arblib.equal(
            Arblib.set!(Arb(), MathConstants.catalan),
            Arblib.set!(Arblib.realref(Acb()), MathConstants.catalan),
        )
        @test Arblib.equal(
            Arblib.set!(Arb(), MathConstants.catalan),
            Arblib.set!(Arblib.arb_struct(), MathConstants.catalan),
        )

        @test Arblib.overlaps(
            Arb(BigFloat(MathConstants.φ)),
            Arblib.set!(Arb(), MathConstants.φ),
        )
        @test Arblib.equal(
            Arblib.set!(Arb(), MathConstants.φ),
            Arblib.set!(Arblib.realref(Acb()), MathConstants.φ),
        )
        @test Arblib.equal(
            Arblib.set!(Arb(), MathConstants.φ),
            Arblib.set!(Arblib.arb_struct(), MathConstants.φ),
        )

        # Intervals: MagLike, ArfLike, BigFloat
        @test Arblib.overlaps(
            Arblib.set!(Arb(), (Mag(1), Mag(2))),
            Arblib.set!(Arb(), (Arf(1), Arf(2))),
        )

        @test Arblib.equal(
            Arblib.set!(Arb(), (Arf(1), Arf(2))),
            Arblib.set!(Arb(), (BigFloat(1), BigFloat(2))),
        )

        @test Arblib.radref(Arblib.set!(Arb(), (Arf(1), Arf(3)))) >= Mag(1)
        @test Arblib.midref(Arblib.set!(Arb(), (Arf(1), Arf(3)))) >= Arf(1)

        @test Arblib.radref(Arblib.set!(Arb(prec = 64), (BigFloat(π), BigFloat(π)))) >
              Mag(0)
        @test iszero(
            Arblib.radref(
                Arblib.set!(
                    Arb(prec = 64),
                    (BigFloat(π), BigFloat(π)),
                    prec = precision(BigFloat),
                ),
            ),
        )

        @test_throws ArgumentError Arblib.set!(Arb(), (Arf(2), Arf(1)))

        # Intervals: General
        @test Arblib.contains(Arblib.set!(Arb(), (1, π)), Arb(π))
        @test Arblib.contains(Arblib.set!(Arb(), (ℯ, Arb(4))), Arb(π))
        @test Arblib.equal(
            Arblib.set!(Arb(), (-2, 2)),
            Arblib.set!(Arb(), (Arf(-2), Arf(2))),
        )

        @test Mag(1e-20) < Arblib.radref(Arblib.set!(Arb(prec = 64), (π, π))) < Mag(1e-10)
        @test Mag(1e-80) <
              Arblib.radref(Arblib.set!(Arb(prec = 64), (π, π), prec = 256)) <
              Mag(1e-70)

        @test_throws ArgumentError Arblib.set!(Arb(), (2, 1))
    end

    @testset "Acb" begin
        @test Arblib.set!(Acb(), BigFloat(1.0)) == one(Acb)
        @test Arblib.equal(Arblib.set!(Arblib.acb_struct(), BigFloat(1.0)), one(Acb))
        x = Arblib.set!(Acb(), ℯ)
        @test isequal(Arblib.realref(x), Arb(ℯ))
        @test Arblib.imagref(x) == zero(Arb)
        @test Arblib.equal(Arblib.set!(Arblib.acb_struct(), BigFloat(1.0)), one(Acb))
        @test Arblib.set!(Acb(), one(Arf).arf) == one(Acb)
        @test Arblib.equal(Arblib.set!(Arblib.acb_struct(), one(Arf).arf), one(Acb))
        @test Arblib.set!(Acb(0, 5), BigFloat(1.0)) == one(Acb)

        @test Arblib.set!(Acb(), BigFloat(1.0), 2) == Acb(1, 2)
        x = Arblib.set!(Acb(), π, Arf(4))
        @test isequal(Arblib.realref(x), Arb(π))
        @test isequal(Arblib.imagref(x), Arb(4))
        @test Arblib.set!(Acb(), Arb(1), Arf(2)) == Acb(1, 2)
        @test Arblib.set!(Acb(), Arb(1).arb, Arf(2).arf) == Acb(1, 2)

        @test Arblib.set!(Acb(), Complex(1, 2)) == Acb(1, 2)
        @test Arblib.set!(Acb(), Complex(BigFloat(1), BigFloat(2))) == Acb(1, 2)
        @test Arblib.equal(Arblib.set!(Arblib.acb_struct(), Complex(1, 2)), Acb(1, 2))

        x = Arblib.set!(Acb(), π)
        @test isequal(Arblib.realref(x), Arb(π))
        @test Arblib.imagref(x) == zero(Arb)
        x = Arblib.set!(Arblib.acb_struct(), π)
        @test Arblib.equal(Arblib.realref(x), Arb(π))
        @test Arblib.equal(Arblib.imagref(x), zero(Arb))

        @test Arblib.equal(Arblib.set!(Acb(), (1, 2)), Arblib.set!(Acb(), Arb((1, 2))))
        @test Arblib.equal(
            Arblib.set!(Acb(), (1, 2), (3, 4)),
            Arblib.set!(Acb(), Arb((1, 2)), Arb((3, 4))),
        )
    end

    @testset "Set Rational" begin
        z = Arb(prec = 64)
        z[] = 1 // 2
        @test z == 0.5

        z = Acb(prec = 64)
        z[] = 1 // 2
        @test z == 0.5

        z = Acb(prec = 64)
        z[] = 1 // 2 + 3 // 4 * im
        @test z == 0.5 + 0.75im
    end

    @testset "Set Int128" begin
        x = 227725055589944414706309
        @test Arb("-227725055589944414706309") == Arb(-x)
        @test Acb("-227725055589944414706309") == Acb(-x)
        @test Arb("227725055589944414706309") == Arb(UInt128(x))
        @test Acb("227725055589944414706309") == Acb(UInt128(x))
    end

end


text:
using Test
using IPFitting, ProgressMeter, JuLIP, SHIPs
using JuLIP: read_dict
using JuLIP.Testing: test_fio
using JuLIP.MLIPs: IPSuperBasis
Fit = IPFitting
DB = IPFitting.DB
Data = Fit.Data

##

function rand_data(sym, N, configtype="rand")
   cubic = (N == 1)
   at = bulk(sym, cubic=cubic) * N
   rattle!(at, 0.1)
   F = (N == 1) ? nothing : rand(3, length(at))
   return Dat(at, configtype; E = rand(), F = F) # , V = rand(3,3))
end

##
println("Double-Check (de-)dictionisation of basis: ")
basis1 = rpi_basis(species = :Ti, N = 2, maxdeg = 8)
basis2 = rpi_basis(species = :Ti, N = 3, maxdeg = 6)
B = IPSuperBasis(basis1, basis2)
println(@test all(test_fio(basis1)))
println(@test all(test_fio(basis2)))
println(@test all(test_fio(B)))

println("Double-Check (de-)dictionisation of Dat: ")
data1 = [ rand_data(:Ti, 3, "md") for n = 1:10 ]
data2 = [ rand_data(:Ti, 1, "cell") for n = 1:10 ]
println(@test read_dict.(write_dict.(data1)) == data1)
println(@test read_dict.(write_dict.(data2)) == data2)

##
println("Create a temporary database.")
dbpath = tempname()
data = [data1; data2]
db = nothing

try
   global db = DB.LsqDB(dbpath, B, data)
   println(@test true)
catch
   @info("...something went wrong...")
   println(@test false)
end

println("checking consistency of db")
println(@test DB.dbpath(db) == dbpath)
println(@test DB.kronfile(dbpath) == dbpath * "_kron.h5")
println(@test DB.infofile(dbpath) == dbpath * "_info.json")
println(@test db.basis == B)
println(@test db.configs == data)

println("re-load the database")
db1 = DB.LsqDB(dbpath)
println(@test DB.dbpath(db1) == dbpath)
println(@test db1.basis == db.basis)
println(@test db1.configs == db.configs)
println(@test db1.Ψ == db.Ψ)

##
println("re-load the database with mmap=true")
db2 = DB.LsqDB(dbpath, mmap=true)
println(@test DB.dbpath(db2) == dbpath)
println(@test db2.basis == db.basis)
println(@test db2.configs == db.configs)
println(@test db2.Ψ == db.Ψ)


text:
using CUDA: has_cuda
using OrderedCollections: OrderedDict

using Oceananigans: AbstractModel, AbstractOutputWriter, AbstractDiagnostic

using Oceananigans.Architectures: AbstractArchitecture
using Oceananigans.Distributed: MultiArch
using Oceananigans.Advection: CenteredSecondOrder
using Oceananigans.BuoyancyModels: validate_buoyancy, regularize_buoyancy, SeawaterBuoyancy
using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions
using Oceananigans.Fields: BackgroundFields, Field, tracernames, VelocityFields, TracerFields, PressureFields
using Oceananigans.Forcings: model_forcing
using Oceananigans.Grids: inflate_halo_size, with_halo, architecture
using Oceananigans.Solvers: FFTBasedPoissonSolver
using Oceananigans.TimeSteppers: Clock, TimeStepper, update_state!
using Oceananigans.TurbulenceClosures: with_tracers, DiffusivityFields, time_discretization, implicit_diffusion_solver
using Oceananigans.LagrangianParticleTracking: LagrangianParticles
using Oceananigans.Utils: tupleit
using Oceananigans.Grids: topology

import Oceananigans.Architectures: architecture

const ParticlesOrNothing = Union{Nothing, LagrangianParticles}

mutable struct NonhydrostaticModel{TS, E, A<:AbstractArchitecture, G, T, B, R, SD, U, C, Φ, F,
                                   V, S, K, BG, P, I, AF} <: AbstractModel{TS}

         architecture :: A        # Computer `Architecture` on which `Model` is run
                 grid :: G        # Grid of physical points on which `Model` is solved
                clock :: Clock{T} # Tracks iteration number and simulation time of `Model`
            advection :: V        # Advection scheme for velocities _and_ tracers
             buoyancy :: B        # Set of parameters for buoyancy model
             coriolis :: R        # Set of parameters for the background rotation rate of `Model`
         stokes_drift :: SD       # Set of parameters for surfaces waves via the Craik-Leibovich approximation
              forcing :: F        # Container for forcing functions defined by the user
              closure :: E        # Diffusive 'turbulence closure' for all model fields
    background_fields :: BG       # Background velocity and tracer fields
            particles :: P        # Particle set for Lagrangian tracking
           velocities :: U        # Container for velocity fields `u`, `v`, and `w`
              tracers :: C        # Container for tracer fields
            pressures :: Φ        # Container for hydrostatic and nonhydrostatic pressure
   diffusivity_fields :: K        # Container for turbulent diffusivities
          timestepper :: TS       # Object containing timestepper fields and parameters
      pressure_solver :: S        # Pressure/Poisson solver
    immersed_boundary :: I        # Models the physics of immersed boundaries within the grid
     auxiliary_fields :: AF       # User-specified auxiliary fields for forcing functions and boundary conditions
end

"""
    NonhydrostaticModel(;
                   grid,
           architecture = CPU(),
                  clock = Clock{eltype(grid)}(0, 0, 1),
              advection = CenteredSecondOrder(),
               buoyancy = nothing,
               coriolis = nothing,
           stokes_drift = nothing,
                forcing = NamedTuple(),
                closure = nothing,
    boundary_conditions = NamedTuple(),
                tracers = (),
            timestepper = :QuasiAdamsBashforth2,
      background_fields = NamedTuple(),
              particles = nothing,
             velocities = nothing,
              pressures = nothing,
     diffusivity_fields = nothing,
        pressure_solver = nothing,
      immersed_boundary = nothing,
       auxiliary_fields = NamedTuple(),
    )

Construct a model for a non-hydrostatic, incompressible fluid, using the Boussinesq approximation
when `buoyancy != nothing`. By default, all Bounded directions are rigid and impenetrable.

Keyword arguments
=================

  - `grid`: (required) The resolution and discrete geometry on which `model` is solved.
  - `architecture`: `CPU()` or `GPU()`. The computer architecture used to time-step `model`.
  - `advection`: The scheme that advects velocities and tracers. See `Oceananigans.Advection`.
  - `buoyancy`: The buoyancy model. See `Oceananigans.BuoyancyModels`.
  - `closure`: The turbulence closure for `model`. See `Oceananigans.TurbulenceClosures`.
  - `coriolis`: Parameters for the background rotation rate of the model.
  - `forcing`: `NamedTuple` of user-defined forcing functions that contribute to solution tendencies.
  - `boundary_conditions`: `NamedTuple` containing field boundary conditions.
  - `tracers`: A tuple of symbols defining the names of the modeled tracers, or a `NamedTuple` of
               preallocated `CenterField`s.
  - `timestepper`: A symbol that specifies the time-stepping method. Either `:QuasiAdamsBashforth2` or
                   `:RungeKutta3`.
"""
function NonhydrostaticModel(;    grid,
                                 clock = Clock{eltype(grid)}(0, 0, 1),
                             advection = CenteredSecondOrder(),
                              buoyancy = nothing,
                              coriolis = nothing,
                          stokes_drift = nothing,
                   forcing::NamedTuple = NamedTuple(),
                               closure = nothing,
       boundary_conditions::NamedTuple = NamedTuple(),
                               tracers = (),
                           timestepper = :QuasiAdamsBashforth2,
         background_fields::NamedTuple = NamedTuple(),
         particles::ParticlesOrNothing = nothing,
                            velocities = nothing,
                             pressures = nothing,
                    diffusivity_fields = nothing,
                       pressure_solver = nothing,
                     immersed_boundary = nothing,
                      auxiliary_fields = NamedTuple(),
    )

    arch = architecture(grid)

    if arch == GPU() && !has_cuda()
         throw(ArgumentError("Cannot create a GPU model. No CUDA-enabled GPU was detected!"))
    end

    tracers = tupleit(tracers) # supports tracers=:c keyword argument (for example)

    validate_buoyancy(buoyancy, tracernames(tracers))
    buoyancy = regularize_buoyancy(buoyancy)

    # Adjust halos when the advection scheme or turbulence closure requires it.
    # Note that halos are isotropic by default; however we respect user-input here
    # by adjusting each (x, y, z) halo individually.
    user_halo = grid.Hx, grid.Hy, grid.Hz
    required_halo = Hx, Hy, Hz = inflate_halo_size(user_halo..., topology(grid), advection, closure)
    if any(user_halo .< required_halo) # Replace grid
        @warn "Inflating model grid halo size to ($Hx, $Hy, $Hz) and recreating grid. " *
              "The model grid will be different from the input grid. To avoid this warning, " *
              "pass halo=($Hx, $Hy, $Hz) when constructing the grid."

        grid = with_halo((Hx, Hy, Hz), grid)
    end

    # Collect boundary conditions for all model prognostic fields and, if specified, some model
    # auxiliary fields. Boundary conditions are "regularized" based on the _name_ of the field:
    # boundary conditions on u, v, w are regularized assuming they represent momentum at appropriate
    # staggered locations. All other fields are regularized assuming they are tracers.
    # Note that we do not regularize boundary conditions contained in *tupled* diffusivity fields right now.
    #
    # First, we extract boundary conditions that are embedded within any _user-specified_ field tuples:
    embedded_boundary_conditions = merge(extract_boundary_conditions(velocities),
                                         extract_boundary_conditions(tracers),
                                         extract_boundary_conditions(pressures),
                                         extract_boundary_conditions(diffusivity_fields))

    # Next, we form a list of default boundary conditions:
    prognostic_field_names = (:u, :v, :w, tracernames(tracers)...)
    default_boundary_conditions = NamedTuple{prognostic_field_names}(FieldBoundaryConditions() for name in prognostic_field_names)

    # Finally, we merge specified, embedded, and default boundary conditions. Specified boundary conditions
    # have precedence, followed by embedded, followed by default.
    boundary_conditions = merge(default_boundary_conditions, embedded_boundary_conditions, boundary_conditions)
    boundary_conditions = regularize_field_boundary_conditions(boundary_conditions, grid, prognostic_field_names)

    # Ensure `closure` describes all tracers
    closure = with_tracers(tracernames(tracers), closure)

    # Either check grid-correctness, or construct tuples of fields
    velocities         = VelocityFields(velocities, grid, boundary_conditions)
    tracers            = TracerFields(tracers,      grid, boundary_conditions)
    pressures          = PressureFields(pressures,  grid, boundary_conditions)
    diffusivity_fields = DiffusivityFields(diffusivity_fields, grid, tracernames(tracers), boundary_conditions, closure)

    if isnothing(pressure_solver)
        pressure_solver = PressureSolver(arch, grid)
    end

    # Materialize background fields
    background_fields = BackgroundFields(background_fields, tracernames(tracers), grid, clock)

    # Instantiate timestepper if not already instantiated
    implicit_solver = implicit_diffusion_solver(time_discretization(closure), arch, grid)
    timestepper = TimeStepper(timestepper, grid, tracernames(tracers), implicit_solver=implicit_solver)

    # Regularize forcing for model tracer and velocity fields.
    model_fields = merge(velocities, tracers)
    forcing = model_forcing(model_fields; forcing...)

    model = NonhydrostaticModel(arch, grid, clock, advection, buoyancy, coriolis, stokes_drift,
                                forcing, closure, background_fields, particles, velocities, tracers,
                                pressures, diffusivity_fields, timestepper, pressure_solver, immersed_boundary,
                                auxiliary_fields)

    update_state!(model)

    return model
end

architecture(model::NonhydrostaticModel) = model.architecture

#####
##### Recursive util for building NamedTuples of boundary conditions from NamedTuples of fields
#####
##### Note: ignores tuples, including tuples of Symbols (tracer names) and
##### tuples of DiffusivityFields (which occur for tupled closures)
#####

extract_boundary_conditions(::Nothing) = NamedTuple()
extract_boundary_conditions(::Tuple) = NamedTuple()

function extract_boundary_conditions(field_tuple::NamedTuple)
    names = propertynames(field_tuple)
    bcs = Tuple(extract_boundary_conditions(field) for field in field_tuple)
    return NamedTuple{names}(bcs)
end

extract_boundary_conditions(field::Field) = field.boundary_conditions


text:
@testset "AABB" begin
    @testset "Construction" begin
        for T in [Float16, Float32, Float64]
            e = one(T)
            z = zero(T)
            a = AABB(SVector(-2e,-2e,-2e),SVector(-e,-e,-e))
            b = AABB(SVector(-2e,-2e,-2e),SVector(-3e,-3e,-3e))
            c = AABB(a, b)
            @test c.pmin == SVector(-3e,-3e,-3e)
            @test c.pmax == SVector(-e,-e,-e)
        end
    end
    @testset "Intersections" begin
        for T in [Float16, Float32, Float64]
            e = one(T)
            for ϵ in [eps(T), e/10]
                z = zero(T)

                finite_aabb = AABB(SVector(-e,-e,-e), SVector(e,e,e))

                r₁ = Ray(SVector(z, z, 2e), SVector(z, z, -e))
                @test intersect(r₁, finite_aabb)

                r₂ = Ray(SVector(z, z, 2e), SVector(z, z, e))
                @test !intersect(r₂, finite_aabb)

                r₃ = Ray(SVector(z, z, 2e), SVector(e,e,-(e+ϵ)))
                @test intersect(r₃, finite_aabb)

                r₄ = Ray(SVector(z, z, 2e), SVector(e,e,-(e-ϵ)))
                @test !intersect(r₄, finite_aabb)

                r₅ = Ray(SVector(z, e-ϵ, 2e), SVector(z,z,-e))
                @test intersect(r₅, finite_aabb)

                r₆ = Ray(SVector(z, e+ϵ, 2e), SVector(z,z,-e))
                @test !intersect(r₆, finite_aabb)
            end
        end
    end
end


text:
# Nonlinear curve fitting.

"""
# Nonlinear multi-variable least squares

Uses a Newton like algorithm to compute the least squares fit
of a model

## Parameters:

 * `x` an array where each colum is a variable
 * `fun` The function that should be fitted (more later on)
 * `a0` An initial guess for each fitting parameter
 * `eps` Convergence criterion
 * `maxiter` Maximum number of iterations

## Specifying the model

The model is specified in argument `fun` that should be callable according
to the following signature:

`r = fun(x, a)`

where both `x` and `a` are vectors

The way the algorithm is implemented allows for the function to be implicit. This means
that we are trying to fit the relationship `fun(x, a) = 0`.

## Initial guess

The initial guess is important and should be as close as possible to the expeted values. It should, at least, give an order of magnitude of each parameter. Zero values are not recommended since in this
case no order of magnitude is available and it is assumed. The initial step is assumed to be `a0/10`.

## Convergence criterion

The initial guess `a0` will probably not fit the data very well. The largest value
of `maxerr = maxabs(fun(x, a0))` (maximum error) is used as a reference and
every time a new approximation is computed, if the change in paramaters is compared to this
reference error. If it is small, `maxabs(da) < eps*maxerr`, the algorithm converged.

## Return values

A tuple containing:

 * The estimated coefficients
 * A `Bool` stating whether the algorithm converged
 * Number of iterations it took to converge.

"""
function nonlinear_fit(x, fun, a0, eps=1e-8, maxiter=200)

    na = length(a0)
    np = size(x, 1)
    nv = size(x, 2)
    for i in 1:na
        if a0[i] == 0
            a0[i] = 0.01
        end
    end

    da = a0 / 10
    a1 = zeros(na)
    a = zeros(na)
    for k = 1:na
        a1[k] = a0[k] - da[k]
    end
    xp = zeros(nv)
    A = zeros(np, na)
    r0 = zeros(np)
    r1 = zeros(np)
    for p = 1:np
        for k = 1:nv
            xp[k] = x[p,k]
        end
        r0[p] = fun(xp, a0)
        r1[p] = fun(xp, a1)
    end
    maxerr = maximum(abs, [maximum(abs, r0), maximum(abs, r1)])
    iter = 1
    convergence = false
    for i = 1:maxiter
        iter = i
        for p = 1:np
            for k = 1:nv
                xp[k] = x[p,k]
            end
            for k = 1:na
                aa = a1[k]
                a1[k] = a1[k] + da[k]
                r = fun(xp, a1)
                a1[k] = aa

                A[p,k] = -(r1[p] - r) / da[k]
            end
        end
        da = A \ r1
        for k = 1:na
            a1[k] = a1[k] - da[k]
        end
        for p = 1:np
            r0[p] = r1[p]
            for k = 1:nv
                xp[k] = x[p,k]
            end
            r1[p] = fun(xp, a1)
        end
        if maximum(abs, r1) < eps * maxerr
            convergence = true
            break
        end
    end

    return a1, convergence, iter
end

"""
   a = gauss_newton_fit(x, y, fun, ∇fun!, a0[[, eps,] maxiter])

Gauss-Newton nonlinear least squares. Given vectors `x` and `y`, the tries to fit parameters `a` to
a function `f` using least squares approximation:

 ``y = f(x, a₁, ..., aₙ)``

For more general approximations, see [`gauss_newton_fit`](@ref).

### Arguments:

 * `x` Vector with x values
 * `y` Vector with y values
 * `fun` a function that is called as `fun(x, a)` where `a` is a vector of parameters.
 * `∇fun!` A function that calculares the derivatives with respect to parameters `a`
 * `a0` Vector with the initial guesses of the parameters
 * `eps` Maximum residual for convergence
 * `maxiter` Maximum number of iterations for convergence

## Return value

A vector with the convrged array. If no convergence is achieved, the function throws an error.

## Specification of the fitting function

The function that should be fitted shoud be specified by Julia funcion with the following signature:

```julia
fun(x::T, a::AbstractVector{T}) where {T<:Number}
```

The derivatives with respect to each fitting parameter `a[i]` should have the following signature:

```julia
∇fun!(x::T, a::AbstractVector{T}, df::AbstractVector{T}) where {T<:Number}
```

No return value is expected and the derivatives are returned in argument `df`.

## Initial approximation (guess)

If the initial approximation is not good enough, divergence is possible.

**Careful** with parameters close to 0. The initial guess should never be 0.0 because the initial
value of the parameter is used as reference value for computing resiuduals.

## Convergence criteria

The argumento `maxiter` specifies the maximum number of iterations that should be carried out.
At each iteration,

``aₖⁿ⁺¹ = aₖⁿ + δₖ``

Convergence is achieved when

``|δᵏ / aₖ⁰| < ε``

## Example
```julia
x = 1.0:10.0
a = [3.0, 2.0, 1.0]
y = a[1] + a[2]*x + a[3]*x^2
fun(x, a) = a[1] + a[2]*x + a[3]*x^2

function ∇fun!(x, a, df)
    df[1] = 1.0
    df[2] = x
    df[3] = x^2
end

a = gauss_newton_fit(x, y, fun, ∇fun!, [0.5, 0.5, 0.5], 1e-8, 30)
```
"""
function gauss_newton_fit(x::AbstractVector{T}, y::AbstractVector{T}, fun, ∇fun!, a0::AbstractVector{T},
                          eps=1e-8, maxiter=200) where {T<:Number}
    P = length(x) # Number of points
    N = length(a0) # Number of parameters

    xi = zero(T)
    df = zeros(T, N)
    a = zeros(T, N)
    for i in 1:N
        a[i] = a0[i]
        if a[i] == 0
            a[i] = 0.01
        end
    end

    A = zeros(T, N, N)
    b = zeros(T, N)

    δref = abs.(a)
    maxerr = zero(T)
    for iter = 1:maxiter
        A .= zero(T)
        b .= zero(T)
        for i in 1:P
            xi = x[i]
            yi = y[i]
            f = fun(xi, a) - yi
            ∇fun!(xi, a, df)

            # Assemble LHS
            for k in 1:N
                for j in 1:N
                    A[j,k] += df[j] * df[k]
                end
            end
            # Assemble RHS
            for j in 1:N
                b[j] -= f * df[j]
            end
        end

        δ = A\b
        a .+= δ

        # Verify convergence:
        maxerr = maximum(abs, δ./δref)
        if maxerr < eps
            return(a)
        end

    end

    error("gauss_newton_fit failed to converge in $maxiter iterations with relative residual of $maxerr !")

    return(a)
end

"""
   a = gauss_newton_generic_fit(x, y, fun, ∇fun!, a0[[, eps,] maxiter])

Gauss-Newton nonlinear least squares. Given matrix `x` the function tries to fit parameters `a` to
an implicit function `f` using least squares approximation:

``f(x₁,..., xₘ, a₁, ..., aₙ) = 0``

This function is very similar to [`gauss_newton_fit`](@ref) but more generic. It doesn't limit the number of
independent variables and doesn't require a `y` LHS.

### Arguments:

 * `x` Matrix where each column is a variable
 * `
 * `fun` a function that is called as `fun(x, a)` where `a` is a vector of parameters.
 * `∇fun!` A function that calculares the derivatives with respect to parameters `a`
 * `a0` Vector with the initial guesses of the parameters
 * `eps` Maximum residual for convergence
 * `maxiter` Maximum number of iterations for convergence

## Return value

A vector with the convrged array. If no convergence is achieved, the function throws an error.

## Specification of the fitting function

The function that should be fitted shoud be specified by Julia funcion with the following signature:

```julia
fun(x::T, a::AbstractVector{T}) where {T<:Number}
```

The derivatives with respect to each fitting parameter `a[i]` should have the following signature:

```julia
∇fun!(x::T, a::AbstractVector{T}, df::AbstractVector{T}) where {T<:Number}
```

No return value is expected and the derivatives are returned in argument `df`.

## Initial approximation (guess)

If the initial approximation is not good enough, divergence is possible.

**Careful** with parameters close to 0. The initial guess should never be 0.0 because the initial
value of the parameter is used as reference value for computing resiuduals.

## Convergence criteria

The argumento `maxiter` specifies the maximum number of iterations that should be carried out.
At each iteration,

``aₖⁿ⁺¹ = aₖⁿ + δₖ``

Convergence is achieved when

``|δᵏ / aₖ⁰| < ε``

## Example
```julia
x = 1.0:10.0
a = [3.0, 2.0, 1.0]
y = a[1] + a[2]*x + a[3]*x^2
fun(x, a) = a[1] + a[2]*x + a[3]*x^2

function ∇fun!(x, a, df)
    df[1] = 1.0
    df[2] = x
    df[3] = x^2
end

a = gauss_newton_fit(x, y, fun, ∇fun!, [0.5, 0.5, 0.5], 1e-8, 30)
```

"""
function gauss_newton_generic_fit(x::AbstractMatrix{T}, fun, ∇fun!, a0::AbstractVector{T},
                          eps=1e-8, maxiter=200) where {T<:Number}
    P = size(x, 1) # Number of points
    M = size(x, 2) # Number of columns (variables)
    N = length(a0) # Number of parameters

    xi = zeros(T, M)
    df = zeros(T, N)
    a = zeros(T, N)
    for i in 1:N
        a[i] = a0[i]
        if a[i] == 0
            a[i] = 0.01
        end
    end

    A = zeros(T, N, N)
    b = zeros(T, N)

    δref = abs.(a)
    maxerr = zero(T)
    for iter = 1:maxiter
        A .= zero(T)
        b .= zero(T)
        for i in 1:P
            for k in 1:M
                xi[k] = x[i,k]
            end
            f = fun(xi, a)

            ∇fun!(xi, a, df)
            # Assemble LHS
            for k in 1:N
                for j in 1:N
                    A[j,k] += df[j] * df[k]
                end
            end
            # Assemble RHS
            for j in 1:N
                b[j] -= f * df[j]
            end
        end

        δ = A\b
        a .+= δ

        # Verify convergence:
        maxerr = maximum(abs, δ./δref)
        if maxerr < eps
            return(a)
        end

    end

    error("gauss_newton_fit failed to converge in $maxiter iterations with relative residual of $maxerr !")

    return(a)
end

"""
   a = secant_nls_fit(x, y, fun, ∇fun!, a0[[, eps,] maxiter])

Secant/Gauss-Newton nonlinear least squares. DOESN'T NEED A DERIVATIVE FUNCTION. Given vectors `x` and `y`, the tries to fit parameters `a` to
a function `f` using least squares approximation:

 ``y = f(x, a₁, ..., aₙ)``

For more general approximations, see [`gauss_newton_fit`](@ref).

### Arguments:

 * `x` Vector with x values
 * `y` Vector with y values
 * `fun` a function that is called as `fun(x, a)` where `a` is a vector of parameters.
 * `∇fun!` A function that calculares the derivatives with respect to parameters `a`
 * `a0` Vector with the initial guesses of the parameters
 * `eps` Maximum residual for convergence
 * `maxiter` Maximum number of iterations for convergence

## Return value

A vector with the convrged array. If no convergence is achieved, the function throws an error.

## Specification of the fitting function

The function that should be fitted shoud be specified by Julia funcion with the following signature:

```julia
fun(x::T, a::AbstractVector{T}) where {T<:Number}
```

The derivatives with respect to each fitting parameter `a[i]` should have the following signature:

```julia
∇fun!(x::T, a::AbstractVector{T}, df::AbstractVector{T}) where {T<:Number}
```

No return value is expected and the derivatives are returned in argument `df`.

## Initial approximation (guess)

If the initial approximation is not good enough, divergence is possible.

**Careful** with parameters close to 0. The initial guess should never be 0.0 because the initial
value of the parameter is used as reference value for computing resiuduals.

## Convergence criteria

The argumento `maxiter` specifies the maximum number of iterations that should be carried out.
At each iteration,

``aₖⁿ⁺¹ = aₖⁿ + δₖ``

Convergence is achieved when

``|δᵏ / aₖ⁰| < ε``

## Example
```julia
x = 1.0:10.0
a = [3.0, 2.0, 1.0]
y = a[1] + a[2]*x + a[3]*x^2
fun(x, a) = a[1] + a[2]*x + a[3]*x^2

a = secant_nls_fit(x, y, fun, ∇fun!, [0.5, 0.5, 0.5], 1e-8, 30)
```
"""

function secant_nls_fit(x::AbstractVector{T}, y::AbstractVector{T}, fun, aguess::AbstractVector{T},
                          eps=1e-8, maxiter=200) where {T<:Number}
    P = length(x) # Number of points
    N = length(aguess) # Number of parameters

    xi = zero(T)
    df = zeros(T, N)
    a = zeros(T, N)
    for i in 1:N
        a[i] = aguess[i]
        if a[i] == 0
            a[i] = 0.01
        end
    end

    δ = a .* (one(T)/20)
    f1 = zeros(T, P)
    a .+= δ
    A = zeros(T, N, N)
    b = zeros(T, N)

    δref = abs.(a)
    maxerr = zero(T)
    for iter = 1:maxiter

        A .= zero(T)
        b .= zero(T)
        f1 .= fun.(x, Ref(a))
        for i in 1:P
            xi = x[i]
            yi = y[i]
            f = f1[i] - yi
            for k in 1:N
                a[k] -= δ[k]
                df[k] = (f1[i] - fun(xi, a)) / δ[k]
                a[k] += δ[k]
            end
            # Assemble LHS
            for k in 1:N
                for j in 1:N
                    A[j,k] += df[j] * df[k]
                end
            end
            # Assemble RHS
            for j in 1:N
                b[j] -= f * df[j]
            end
        end
        δ = A\b
        a .+= δ
        # Verify convergence:
        maxerr = maximum(abs, δ./δref)
        if maxerr < eps
            return(a)
        end

    end

    error("gauss_newton_fit failed to converge in $maxiter iterations with relative residual of $maxerr !")

    return(a)
end


text:
module NNLSSolverInterface

export NNLSSolver

using LinearAlgebra
using NNLS: QPWorkspace, load!, solve!
import MathProgBase.SolverInterface
const MSI = SolverInterface

struct NNLSSolver <: MSI.AbstractMathProgSolver
end

mutable struct NNLSModel <: MSI.AbstractLinearQuadraticModel
    workspace::QPWorkspace{Float64, Int}
    lb::Vector{Float64}
    ub::Vector{Float64}
    A::Matrix{Float64}
    constr_lb::Vector{Float64}
    constr_ub::Vector{Float64}
    Q::Matrix{Float64}
    q::Vector{Float64}
    sense::Symbol
    solution::Vector{Float64}
    duals::Vector{Float64}

    function NNLSModel()
        m = new(QPWorkspace{Float64, Int}(0, 0))
        resize!(m, 0, 0)
        m
    end
end

function Base.resize!(m::NNLSModel, q::Integer, n::Integer)
    resize!(m.workspace, q, n)
    m.lb = zeros(n)
    m.ub = zeros(n)
    m.A = zeros(q, n)
    m.constr_lb = zeros(q)
    m.constr_ub = zeros(q)
    m.Q = zeros(n, n)
    m.q = zeros(n)
    m.sense = :Min
    m.solution = fill(NaN, n)
    m.duals = fill(NaN, q)
end

MSI.LinearQuadraticModel(s::NNLSSolver) = NNLSModel()

MSI.numvar(work::QPWorkspace) = size(work.G, 2)
MSI.numconstr(work::QPWorkspace) = size(work.G, 1)

MSI.numvar(m::NNLSModel) = MSI.numvar(m.workspace)
MSI.numconstr(m::NNLSModel) = MSI.numconstr(m.workspace)

function MSI.loadproblem!(m::NNLSModel, A, lb, ub, obj, constr_lb, constr_ub, sense)
    q, n = size(A)
    if (q != size(m.A, 1)) || (n != size(m.A, 2))
        resize!(m, q, n)
    end
    m.sense = sense
    m.A .= A
    m.lb .= lb
    m.ub .= ub
    m.constr_lb .= constr_lb
    m.constr_ub .= constr_ub
    m.q .= obj
    m.Q .= 0
    m.solution .= NaN
    m.duals .= NaN
    m.workspace.status = :Unsolved
end

MSI.setquadobj!(m::NNLSModel, Q) = m.Q = Q
function MSI.setquadobj!(m::NNLSModel, rowidx, colidx, quadval)
    m.Q .= 0
    for i in 1:length(rowidx)
        x = m.Q[rowidx[i], colidx[i]] + quadval[i]
        m.Q[rowidx[i], colidx[i]] = x
        m.Q[colidx[i], rowidx[i]] = x
    end
end

function MSI.optimize!(m::NNLSModel)
    nvars = size(m.A, 2)
    nrows = size(m.A, 1)
    @assert nrows == length(m.constr_lb)
    @assert nrows == length(m.constr_ub)
    @assert nvars == length(m.lb)
    @assert nvars == length(m.ub)
    @assert nvars == length(m.q)

    nconstr = 2 * nrows + 2 * nvars

    if nvars != MSI.numvar(m.workspace) || nconstr != MSI.numconstr(m.workspace)
        resize!(m.workspace, nconstr, nvars)
        m.solution = fill(NaN, nvars)
    end

    m.workspace.L .= m.Q
    m.workspace.c .= m.q

    m.workspace.G .= 0

    for j in 1:nvars, i in 1:nrows
        m.workspace.G[i, j] = m.A[i, j]
    end
    for i in 1:nrows
        m.workspace.g[i] = m.constr_ub[i]
    end
    for j in 1:nvars, i in 1:nrows
        m.workspace.G[i + nrows, j] = -m.A[i, j]
    end
    for i in 1:nrows
        m.workspace.g[i + nrows] = -m.constr_lb[i]
    end
    for j in 1:nvars
        m.workspace.G[2 * nrows + j, j] = 1
    end
    for i in 1:nvars
        m.workspace.g[2 * nrows + i] = m.ub[i]
    end
    for i in 1:nvars
        m.workspace.G[2 * nrows + nvars + i, i] = -1
    end
    for i in 1:nvars
        m.workspace.g[2 * nrows + nvars + i] = -m.lb[i]
    end

    for i in 1:length(m.workspace.g)
        if isinf(m.workspace.g[i])
            m.workspace.G[i, :] .= 0
            m.workspace.g[i] = 0
        end
    end

    m.workspace.status = :Unsolved

    m.solution, λ = solve!(m.workspace)
    for i in 1:nrows
        if abs(λ[i]) > abs(λ[i + nrows])
            m.duals[i] = λ[i]
        else
            m.duals[i] = λ[i + nrows]
        end
    end
end

MSI.getsolution(m::NNLSModel) = m.solution
MSI.getobjval(m::NNLSModel) = 0.5 * m.solution' * m.Q * m.solution + m.q' * m.solution
MSI.getconstrduals(m::NNLSModel) = m.duals

MSI.status(m::NNLSModel) = m.workspace.status

end


text:
"""
    Selector

Abstract supertype for all selectors.

Selectors are wrappers that indicate that passed values are not the array indices,
but values to be selected from the dimension index, such as `DateTime` objects for
a `Ti` dimension.

Selectors provided in DimensionalData are:

- [`At`](@ref)
- [`Between`](@ref)
- [`Near`](@ref)
- [`Where`](@ref)
- [`Contains`](@ref)

"""
abstract type Selector{T} end

val(sel::Selector) = sel.val
rebuild(sel::Selector, val) = basetypeof(sel)(val)

const SelTuple = Tuple{<:Selector,Vararg{<:Selector}}

"""
    At <: Selector

    At(x, atol, rtol)
    At(x; atol=nothing, rtol=nothing)

Selector that exactly matches the value on the passed-in dimensions, or throws an error.
For ranges and arrays, every intermediate value must match an existing value -
not just the end points.

`x` can be any value or `Vector` of values.

`atol` and `rtol` are passed to `isapprox`.
For `Number` `rtol` will be set to `Base.rtoldefault`, otherwise `nothing`,
and wont be used.

## Example

```jldoctest
using DimensionalData

A = DimArray([1 2 3; 4 5 6], (X(10:10:20), Y(5:7)))
A[X(At(20)), Y(At(6))]

# output

5
```
"""
struct At{T,A,R} <: Selector{T}
    val::T
    atol::A
    rtol::R
end
At(val; atol=nothing, rtol=nothing) =
    At{typeof.((val, atol, rtol))...}(val, atol, rtol)

atol(sel::At) = sel.atol
rtol(sel::At) = sel.rtol

"""
    Near <: Selector

    Near(x)

Selector that selects the nearest index to `x`.

With [`Points`](@ref) this is simply the index values nearest to the `x`,
however with [`Intervals`](@ref) it is the interval _center_ nearest to `x`.
This will be offset from the index value for `Start` and
[`End`](@ref) loci.

## Example

```jldoctest
using DimensionalData

A = DimArray([1 2 3; 4 5 6], (X(10:10:20), Y(5:7)))
A[X(Near(23)), Y(Near(5.1))]

# output

4
```
"""
struct Near{T} <: Selector{T}
    val::T
end

"""
    Contains <: Selector

    Contains(x)

Selector that selects the interval the value is contained by. If the
interval is not present in the index, an error will be thrown.

Can only be used for [`Intervals`](@ref) or [`Categorical`](@ref).

## Example

```jldoctest
using DimensionalData

dims_ = X(10:10:20; mode=Sampled(sampling=Intervals())),
        Y(5:7; mode=Sampled(sampling=Intervals()))
A = DimArray([1 2 3; 4 5 6], dims_)
A[X(Contains(8)), Y(Contains(6.8))]

# output

3
```
"""
struct Contains{T} <: Selector{T}
    val::T
end

"""
    Between <: Selector

    Between(a, b)

Selector that retreive all indices located between 2 values,
evaluated with `>=` for the lower value, and `<` for the upper value.
This means the same value will not be counted twice in 2 adjacent
`Between` selections.

For [`Intervals`](@ref) the whole interval must be lie between the
values. For [`Points`](@ref) the points must fall between
the values. Different [`Sampling`](@ref) types may give different
results with the same input - this is the intended behaviour.

`Between` for [`Irregular`](@ref) intervals is a little complicated. The
interval is the distance between a value and the next (for `Start` locus)
or previous (for [`End`](@ref) locus) value.

For [`Center`](@ref), we take the mid point between two index values
as the start and end of each interval. This may or may not make sense for
the values in your indes, so use `Between` with `Irregular` `Intervals(Center())`
with caution.

## Example

```jldoctest
using DimensionalData

A = DimArray([1 2 3; 4 5 6], (X(10:10:20), Y(5:7)))
A[X(Between(15, 25)), Y(Between(4, 6.5))]

# output

1×2 DimArray{Int64,2} with dimensions:
  X: 20:10:20 Sampled: Ordered Regular Points,
  Y: 5:6 Sampled: Ordered Regular Points
 4  5
```
"""
struct Between{T<:Union{Tuple{Any,Any},Nothing}} <: Selector{T}
    val::T
end
Between(args...) = Between(args)

Base.first(sel::Between) = first(val(sel))
Base.last(sel::Between) = last(val(sel))

"""
    Where <: Selector

    Where(f::Function)

Selector that filters a dimension by any function that accepts
a single value from the index and returns a `Bool`.

## Example

```jldoctest
using DimensionalData

A = DimArray([1 2 3; 4 5 6], (X(10:10:20), Y(19:21)))
A[X(Where(x -> x > 15)), Y(Where(x -> x in (19, 21)))]

# output

1×2 DimArray{Int64,2} with dimensions:
  X: Int64[20] Sampled: Ordered Regular Points,
  Y: Int64[19, 21] Sampled: Ordered Regular Points
 4  6
```
"""
struct Where{T} <: Selector{T}
    f::T
end

val(sel::Where) = sel.f

# sel2indices ==========================================================================

# Converts Selectors to regular indices

@inline sel2indices(x, ls...) = sel2indices(dims(x), ls...)
@inline sel2indices(dims::Tuple, l1, ls...) = sel2indices(dims, (l1, ls...))
@inline sel2indices(dims::Tuple, lookup::Tuple) =
    map((d, l) -> sel2indices(d, l), dims, lookup)
@inline sel2indices(dims::Tuple, lookup::Tuple{}) = ()
@inline sel2indices(dim::Dimension, sel::Selector) = _sel2indices(dim, sel)
@inline sel2indices(dim::Dimension, sel::Val) = _sel2indices(dim, At(sel))
# Standard indices are just returned.
@inline sel2indices(::Dimension, sel::StandardIndices) = sel
@inline function sel2indices(::Dimension, sel)
    selstr = sprint(show, sel)
    throw(ArgumentError("Invalid index `$selstr`. Did you mean `At($selstr)`? Use stardard indices, `Selector`s, or `Val` for compile-time `At`."))
end

# First filter based on rough selector properties -----------------

# Vectors are mapped
@inline _sel2indices(dim::Dimension, sel::Selector{<:AbstractVector}) =
    [_sel2indices(mode(dim), dim, rebuild(sel, v)) for v in val(sel)]
@inline _sel2indices(dim::Dimension, sel::Selector) = _sel2indices(mode(dim), dim, sel)

# Where selector ==============================
# Yes this is everything. Where doesn't need mode specialisation
@inline _sel2indices(dim::Dimension, sel::Where) =
    [i for (i, v) in enumerate(index(dim)) if sel.f(v)]

# Then dispatch based on IndexMode -----------------
# Selectors can have varied behaviours depending on the index mode.

# Noindex Contains just converts the selector to standard indices. Implemented
# so the Selectors actually work, not because what they do is useful or interesting.
@inline _sel2indices(mode::NoIndex, dim::Dimension, sel::Union{At,Near,Contains}) = val(sel)
@inline _sel2indices(mode::NoIndex, dim::Dimension, sel::Union{Between}) =
    val(sel)[1]:val(sel)[2]
@inline _sel2indices(mode::Categorical, dim::Dimension, sel::Selector) =
    if sel isa Union{Contains,Near}
        _sel2indices(Points(), mode, dim, At(val(sel)))
    else
        _sel2indices(Points(), mode, dim, sel)
    end
@inline _sel2indices(mode::AbstractSampled, dim::Dimension, sel::Selector) =
    _sel2indices(sampling(mode), mode, dim, sel)

# For Sampled filter based on sampling type and selector -----------------

@inline _sel2indices(sampling::Sampling, mode::IndexMode, dim::Dimension, sel::At) =
    at(sampling, mode, dim, sel)
@inline _sel2indices(sampling::Sampling, mode::IndexMode, dim::Dimension, sel::Near) = begin
    span(mode) isa Irregular && locus(mode) isa Union{Start,End} && _nearirregularerror()
    near(sampling, mode, dim, sel)
end
@inline _sel2indices(sampling::Points, mode::IndexMode, dim::Dimension, sel::Contains) =
    _containspointserror()
@inline _sel2indices(sampling::Intervals, mode::IndexMode, dim::Dimension, sel::Contains) =
    contains(sampling, mode, dim, sel)
@inline _sel2indices(sampling::Sampling, mode::IndexMode, dim::Dimension, sel::Between{<:Tuple}) =
    between(sampling, mode, dim, sel)

@noinline _nearirregularerror() =
    throw(ArgumentError("Near is not implemented for Irregular with Start or End loci. Use Contains"))
@noinline _containspointserror() =
    throw(ArgumentError("`Contains` has no meaning with `Points`. Use `Near`"))

# Unaligned IndexMode ------------------------------------------

# unalligned2indices is callled directly from dims2indices

# We use the transformation from the first Transformed dim.
# In practice the others could be empty.
@inline unalligned2indices(dims::DimTuple, sel::Tuple) = sel
@inline unalligned2indices(dims::DimTuple, sel::Tuple{<:Dimension,Vararg{<:Dimension}}) =
    unalligned2indices(dims, map(val, sel))
@inline unalligned2indices(dims::DimTuple, sel::Tuple{<:Selector,Vararg{<:Selector}}) = begin
    coords = [map(val, sel)...]
    transformed = transformfunc(mode(dims[1]))(coords)
    map(_transform2int, sel, transformed)
end

_transform2int(::At, x) = convert(Int, x)
_transform2int(::Near, x) = round(Int, x)

# Selector methods

function hasselection(x, selectors::Union{DimTuple,SelTuple,Selector,Dimension})
    hasselection(dims(x), selectors)
end
hasselection(x::Nothing, selectors::Union{DimTuple,SelTuple,Selector,Dimension}) = false
function hasselection(dims::DimTuple, seldims::DimTuple)
    sorted = DD.dims(seldims, dims)
    hasselection(DD.dims(dims, sorted), map(val, sorted))
end
hasselection(dims::DimTuple, selectors::SelTuple) = all(map(hasselection, dims, selectors))
function hasselection(dims::DimTuple, selector::Dimension)
    hasselection(DD.dims(dims, selector), selector)
end
function hasselection(dims::DimTuple, selector::Selector)
    throw(ArgumentError("Cannot select from multiple Dimensions with a single Selector"))
end
hasselection(dim::Dimension, seldim::Dimension) = hasselection(dim, val(seldim))
hasselection(dim::Dimension, sel::At) = at(dim, sel; err=_False()) === nothing ? false : true
hasselection(dim::Dimension, sel::Contains) = contains(dim, sel; err=_False()) === nothing ? false : true
# Near an Between only fail on UnorderedIndex
# Otherwise Near returns the nearest index, and Between and empty range
hasselection(dim::Dimension, selnear::Near) = order(dim) isa UnorderedIndex ? false : true
hasselection(dim::Dimension, selnear::Between) = order(dim) isa UnorderedIndex ? false : true

# Selector methods

# at =============================================================================

struct _True end
struct _False end

at(dim::Dimension, sel::At; kw...) = at(sampling(mode(dim)), mode(dim), dim, sel; kw...)
function at(::Sampling, mode::IndexMode, dim::Dimension, sel::At; kw...)
    at(indexorder(dim), val(dim), dim, val(sel), atol(sel), rtol(sel); kw...)
end
#
function at(
    o::IndexOrder, a::AbstractArray{<:Union{Number,Dates.TimeType}}, dim::Dimension,
    selval, atol, rtol::Nothing;
    err=_True()
)
    x = unwrap(selval)
    i = searchsortedlast(a, x; order=_ordering(o))
    if i === 0 || !_is_at(x, a[i], atol)
        # Try again with the next index
        i1 = i + 1
        if checkbounds(Bool, a, i1) && _is_at(x, a[i1], atol)
            return relate(dim, i1)
        else
            return _selnotfound_or_nothing(err, dim, selval)
        end
    else
        return relate(dim, i)
    end
end

@inline _is_at(x, y, atol) = x == y
@inline _is_at(x::Real, y::Real, atol::Real) = abs(x - y) <= atol

# catch-all for an unordered or non-number index
function at(::IndexOrder, a, dim::Dimension, selval, atol, rtol::Nothing; err=_True())
    i = findfirst(x -> x == unwrap(selval), index(dim))
    if i === nothing
        return _selnotfound_or_nothing(err, dim, selval)
    else
        return relate(dim, i)
    end
end
# compile-time indexing
@generated function at(
    ::IndexOrder, ::Val{Index}, dim::Dimension, selval::Val{X}, atol, rtol::Nothing;
    err=_True()
) where {Index,X}
    i = findfirst(x -> x == X, Index)
    if i === nothing
        return :(_selnotfound_or_nothing(err, dim, selval))
    else
        return :(relate(dim, $i))
    end
end

_selnotfound_or_nothing(err::_True, dim, selval) = _selvalnotfound(dim, selval)
_selnotfound_or_nothing(err::_False, dim, selval) = nothing
@noinline _selvalnotfound(dim, selval) = throw(ArgumentError("$selval not found in $(name(dim))"))

# near ===========================================================================

# Finds the nearest point in the index, adjusting for locus if necessary.
# In Intevals we are finding the nearest point to the center of the interval.

near(dim::Dimension, sel::Near) = near(sampling(mode(dim)), mode(dim), dim, sel)
function near(::Sampling, mode::IndexMode, dim::Dimension, sel::Near)
    order = indexorder(dim)
    order isa UnorderedIndex && _nearunorderederror()
    locus = DD.locus(dim)

    v = _locus_adjust(locus, unwrap(val(sel)), dim)
    i = _inbounds(_searchorder(order)(order, dim, v), dim)
    i = if (order isa ForwardIndex ? (<=) : (>=))(i, _dimlower(order, dim))
        _dimlower(order, dim)
    else
        previ = _prevind(order, i)
        vl, vi = map(abs, (dim[previ] - v, dim[i] - v))
        # We have to use the right >/>= for Start/End locus
        _lt(locus)(vl, vi) ? previ : i
    end
    relate(dim, i)
end

@noinline _nearunorderederror() = throw(ArgumentError("`Near` has no meaning in an `Unordered` index"))

_locus_adjust(locus::Center, v, dim) = v
_locus_adjust(locus::Start, v, dim) = v - abs(step(dim)) / 2
_locus_adjust(locus::End, v, dim) = v + abs(step(dim)) / 2
_locus_adjust(locus::Start, v::DateTime, dim) = v - (v - (v - abs(step(dim)))) / 2
_locus_adjust(locus::End, v::DateTime, dim) = v + (v + abs(step(dim)) - v) / 2

# contains ================================================================================

# Finds which interval contains a point

function contains(dim::Dimension, sel::Contains; kw...)
    contains(sampling(mode(dim)), mode(dim), dim, sel; kw...)
end
# Points --------------------------------------
@noinline function contains(::Points, ::IndexMode, dim::Dimension, sel::Contains; err=_True())
    if err isa _True
        throw(ArgumentError("Points IndexMode cannot use 'Contains', use 'Near' instead."))
    else
        nothing
    end
end
# Intervals -----------------------------------
function contains(sampling::Intervals, mode::IndexMode, dim::Dimension, sel::Contains; kw...)
    contains(span(mode), sampling, indexorder(mode), locus(mode), dim, sel; kw...)
end
# Regular Intervals ---------------------------
function contains(span::Regular, ::Intervals, order, locus, dim::Dimension, sel::Contains;
    err=_True()
)
    _locus_checkbounds(locus, dim, sel) || return _boundserror_or_nothing(err)
    v = val(sel)
    absstep = abs(val(span))
    i = _whichsearch(locus, order)(order, dim, _maybeaddhalf(locus, absstep, v))
    # Check the value is in this cell.
    # It is always for AbstractRange but might not be for Val tuple or Vector.
    if (val(dim) isa AbstractRange) || _lt(locus)(v, dim[i] + absstep)
        return relate(dim, i)
    else
        return _notcontained_or_nothing(err, v)
    end
end
# Explicit Intervals ---------------------------
function contains(span::Explicit, ::Intervals, order, locus, dim::Dimension, sel::Contains;
    err=_True()
)
    x = val(sel)
    i = searchsortedlast(view(val(span), 1, :), x; order=_ordering(order))
    if i === 0 || val(span)[2, i] < x
        return _notcontained_or_nothing(err, x)
    else
        return relate(dim, i)
    end
end
# Irregular Intervals -------------------------
function contains(
    span::Irregular, ::Intervals, order::IndexOrder,
    locus::Locus, dim::Dimension, sel::Contains;
    err=_True()
)
    _locus_checkbounds(locus, dim, sel) || return _boundserror_or_nothing(err)
    i = _whichsearch(locus, order)(order, dim, val(sel))
    return relate(dim, i)
end
function contains(
    span::Irregular, ::Intervals, order::IndexOrder,
    locus::Center, dim::Dimension, sel::Contains;
    err=_True()
)
    _locus_checkbounds(locus, dim, sel) || return _boundserror_or_nothing(err)
    v = val(sel)
    i = _searchfirst(order, dim, v)
    i = if i <= firstindex(dim)
        firstindex(dim)
    elseif i > lastindex(dim)
        lastindex(dim)
    else
        interval = abs(dim[i] - dim[i - 1])
        distance = abs(dim[i] - v)
        _order_lt(order)(interval / 2, distance) ? i - 1 : i
    end
    return relate(dim, i)
end

_boundserror_or_nothing(err::_True) = throw(BoundsError())
_boundserror_or_nothing(err::_False) = nothing

_notcontained_or_nothing(err::_True, selval) = _notcontainederror(selval)
_notcontained_or_nothing(err::_False, selval) = nothing
@noinline _notcontainederror(v) = throw(ArgumentError("No interval contains $(v)"))

_whichsearch(::Locus, ::ForwardIndex) = _searchlast
_whichsearch(::Locus, ::ReverseIndex) = _searchfirst
_whichsearch(::End, ::ForwardIndex) = _searchfirst
_whichsearch(::End, ::ReverseIndex) = _searchlast

_maybeaddhalf(::Locus, s, v) = v
_maybeaddhalf(::Center, s, v) = v + s / 2

_order_lt(::ForwardIndex) = (<)
_order_lt(::ReverseIndex) = (<=)

# between ================================================================================

# Finds all values between two points, adjusted for locus where necessary

struct _Upper end
struct _Lower end

between(dim::Dimension, sel::Between) = between(sampling(mode(dim)), mode(dim), dim, sel)
function between(sampling::Sampling, mode::IndexMode, dim::Dimension, sel::Between)
    order = indexorder(dim)
    order isa UnorderedIndex && throw(ArgumentError("Cannot use `Between` with UnorderedIndex"))
    a, b = between(sampling, order, mode, dim, sel)
    relate(dim, a:b)
end
# Points ------------------------------------
function between(sampling::Points, o::IndexOrder, ::IndexMode, dim::Dimension, sel::Between)
    b1, b2 = _maybeflipbounds(o, _sorttuple(sel))
    s1, s2 = _maybeflipbounds(o, (_searchfirst, _searchlast))
    _inbounds((s1(o, dim, b1), s2(o, dim, b2)), dim)
end
# Intervals -------------------------
function between(sampling::Intervals, o::IndexOrder, mode::IndexMode, dim::Dimension, sel::Between)
    between(span(mode), sampling, o, mode, dim, sel)
end
# Regular Intervals -------------------------
function between(span::Regular, ::Intervals, o::IndexOrder, mode::IndexMode, dim::Dimension, sel::Between)
    b1, b2 = _maybeflipbounds(o, _sorttuple(sel) .+ _locus_adjust(mode))
    _inbounds((_searchfirst(o, dim, b1), _searchlast(o, dim, b2)), dim)
end
# Explicit Intervals -------------------------
function between(span::Explicit, ::Intervals, o::IndexOrder, mode::IndexMode, dim::Dimension, sel::Between)
    _inbounds(
        (searchsortedfirst(view(val(span), 1, :), first(val(sel)); order=_ordering(o), lt=<),
         searchsortedlast(view(val(span), 2, :), last(val(sel)); order=_ordering(o), lt=<)),
        dim
    )
end
# Irregular Intervals -----------------------
function between(span::Irregular, ::Intervals, o::IndexOrder, mode::IndexMode, d::Dimension, sel::Between)
    l, h = _sorttuple(sel)
    bl, bh = bounds(span)
    a = l <= bl ? _dimlower(o, d) : between(_Lower(), locus(mode), o, d, l)
    b = h >= bh ? _dimupper(o, d) : between(_Upper(), locus(mode), o, d, h)
    _maybeflipbounds(o, (a, b))
end
function between(x, locus::Union{Start,End}, o::IndexOrder, d::Dimension, v)
    _search(x, o, d, v) - _ordscalar(o) * (_locscalar(locus) + _endshift(x))
end
function between(x, locus::Center, o::IndexOrder, d::Dimension, v)
    r = _ordscalar(o); sh = _endshift(x)
    i = _search(x, o, d, v)
    interval = abs(d[i] - d[i-r])
    distance = abs(d[i] - v)
    # Use the right >/>= to match interval bounds
    _lt(x)(distance, (interval / 2)) ? i - sh * r : i - (1 + sh) * r
end

_locus_adjust(mode) = _locus_adjust(locus(mode), abs(step(span(mode))))
_locus_adjust(locus::Start, step) = zero(step), -step
_locus_adjust(locus::Center, step) = step/2, -step/2
_locus_adjust(locus::End, step) = step, zero(step)

_locscalar(::Start) = 1
_locscalar(::End) = 0
_endshift(::_Lower) = -1
_endshift(::_Upper) = 1
_ordscalar(::ForwardIndex) = 1
_ordscalar(::ReverseIndex) = -1

_search(x, order, dim, v) = _inbounds(_searchorder(order)(order, dim, v, _lt(x)), dim)

_lt(::_Lower) = (<)
_lt(::_Upper) = (<=)

_maybeflipbounds(m::IndexMode, bounds) = _maybeflipbounds(indexorder(m), bounds)
_maybeflipbounds(o::ForwardIndex, (a, b)) = (a, b)
_maybeflipbounds(o::ReverseIndex, (a, b)) = (b, a)

# Shared utils ============================================================================

_searchlast(o::IndexOrder, dim::Dimension, v, lt=<) =
    searchsortedlast(index(dim), unwrap(v); order=_ordering(o), lt=lt)
_searchlast(o::IndexOrder, dim::Dimension{<:Val{Index}}, v, lt=<) where Index =
    searchsortedlast(Index, unwrap(v); order=_ordering(o), lt=lt)

_searchfirst(o::IndexOrder, dim::Dimension, v, lt=<) =
    searchsortedfirst(index(dim), unwrap(v); order=_ordering(o), lt=lt)
_searchfirst(o::IndexOrder, dim::Dimension{<:Val{Index}}, v, lt=<) where Index =
    searchsortedfirst(Index, unwrap(v); order=_ordering(o), lt=lt)

_asfunc(::Type{typeof(<)}) = <
_asfunc(::Type{typeof(<=)}) = <=

# Return an inbounds index
_inbounds(is::Tuple, dim::Dimension) = map(i -> _inbounds(i, dim), is)
_inbounds(i::Int, dim::Dimension) =
    if i > lastindex(dim)
        lastindex(dim)
    elseif i <= firstindex(dim)
        firstindex(dim)
    else
        i
    end

_sorttuple(sel::Between) = _sorttuple(val(sel))
_sorttuple((a, b)) = a < b ? (a, b) : (b, a)

_lt(::Locus) = (<)
_lt(::End) = (<=)
_gt(::Locus) = (>=)
_gt(::End) = (>)

_locus_checkbounds(loc, dim::Dimension, sel::Selector) =  _locus_checkbounds(loc, bounds(dim), val(sel))
_locus_checkbounds(loc, (l, h)::Tuple, v) = !(_lt(loc)(v, l) || _gt(loc)(v, h))

_prevind(::ForwardIndex, i) = i - 1
_prevind(::ReverseIndex, i) = i + 1

_dimlower(o::ForwardIndex, d) = firstindex(d)
_dimlower(o::ReverseIndex, d) = lastindex(d)
_dimupper(o::ForwardIndex, d) = lastindex(d)
_dimupper(o::ReverseIndex, d) = firstindex(d)

_searchorder(::ForwardIndex) = _searchfirst
_searchorder(::ReverseIndex) = _searchlast

_ordering(::ForwardIndex) = Base.Order.ForwardOrdering()
_ordering(::ReverseIndex) = Base.Order.ReverseOrdering()