text:
#Finds the representation of a permutation for a particular partition
function example1()
	example1(5,[3, 2],[2,3,1,4,5])
end

# Parameters:
#	N::Int
#	- the problem size
#	partition::Array{Int, 1}
#	- a partition of N
#	permutation::Array{Int, 1}
#	- a permutation of N
function example1(N::Int, partition::Array{Int, 1}, permutation::Array{Int, 1})

	P, WI = partitions(N) #Creates the partitions of i for i is 1 through N
	p = 1 #Find the index of the specified partitions
	while P[N][p] != partition
		p += 1
	end

	YS = ys_partition(N, partition) #Creates the Yamanouchi Symbols of the specified partition
	degree = length(YS)

	YOR, PT = yor(N) #Creates Young's Orthogonal Representations and the Pointer Tree for partition decomposistion
	YORnp = YOR[N][p] #Get the Representation of each Adjacent Transposition of the specified partition

	RM = yor_permutation(permutation, YORnp) #Creates the Representation Matrix of the specified permutation for the specified partition

	println("Finds the representation of a permutation for a particular partition")
	println("")
	println("Permutation: ")
	println(permutation_string(permutation))
	println("")
	println("Partition: ")
	println(partition_string(partition))
	println("")
	println("Degree: ")
	println(degree)
	println("")
	println("Standard Tableaux: ")
	for i = 1:degree #Converts the Yamanouchi Symbols into Standard Tableau and prints them
		print_ys(YS[i])
		println("")
	end
	println("Representation matrix of the permutation for this partition: ")
	print(full(round(RM,5)))
end

#If a and b are permutations and c = a * b, demonstrates that the representation of c is the representation of a multiplied by the representation of b
function example2()
	example2(5,[3,2],[1,2,4,3,5],[2,3,1,4,5])
end

# Parameters:
#	N::Int
#	- the problem size
#	partition::Array{Int, 1}
#	- a partition of N
#	p1::Array{Int, 1}
#	- the first permutation of N
#	p2::Array{Int, 1}
#	- the second permutation of N
function example2(N::Int, partition::Array{Int, 1}, p1::Array{Int, 1}, p2::Array{Int, 1})

	P, WI = partitions(N) #Creates the partitions of i for i is 1 through N
	p = 1 #Find the index of the specified partitions
	while P[N][p] != partition
		p += 1
	end

	pm = sn_multiply(p1, p2) #Find the product of the two permutations

	YOR, PT = yor(N) #Creates Young's Orthogonal Representations and the Pointer Tree for partition decomposistion
	YORnp = YOR[N][p] #Get the Representation of each Adjacent Transposition of the specified partition

	RM1 = yor_permutation(p1, YORnp) #Find the representation of the first permutation
	RM2 = yor_permutation(p2, YORnp) #Find the representation of the second permutation
	RMm = yor_permutation(pm, YORnp) #Find the representation of the product of the permutations
	error = RM1 * RM2 - RMm #Find the error in the process

	println("If a and b are permutations and c = a * b, demonstrates that the representation of c is the representation of a multiplied by the representation of b")
	println("")
	ST1 = permutation_string(p1)
	println("Representation of ", ST1)
	println(full(round(RM1, 5)))
	println("")
	ST2 = permutation_string(p2)
	println("Representation of ", ST2)
	println(full(round(RM2, 5)))
	println("")
	STm = permutation_string(pm)
	println(ST1," * ", ST2, " is ", STm)
	println("")
	println("Representation of ", STm)
	println(full(round(RMm, 5)))
	println("")
	println("Maximum error: ")
	println(maximum(abs(error)))
end

#Demonstrates how to put the values of a function on Sn in the order used to compute the FFT
function example3()

	PA = Array(Array{Int, 1}, 24) #Create an array to hold all of the permutations of 1,2,3,4
	PA[1] = [1,2,3,4]
	PA[2] = [1,2,4,3]
	PA[3] = [1,3,2,4]
	PA[4] = [1,3,4,2]
	PA[5] = [1,4,2,3]
	PA[6] = [1,4,3,2]
	PA[7] = [2,1,3,4]
	PA[8] = [2,1,4,3]
	PA[9] = [2,3,1,4]
	PA[10] = [2,3,4,1]
	PA[11] = [2,4,1,3]
	PA[12] = [2,4,3,1]
	PA[13] = [3,1,2,4]
	PA[14] = [3,1,4,2]
	PA[15] = [3,2,1,4]
	PA[16] = [3,2,4,1]
	PA[17] = [3,4,1,2]
	PA[18] = [3,4,2,1]
	PA[19] = [4,1,2,3]
	PA[20] = [4,1,3,2]
	PA[21] = [4,2,1,3]
	PA[22] = [4,2,3,1]
	PA[23] = [4,3,1,2]
	PA[24] = [4,3,2,1]

	VA = Array(Float64, 24) #Create the values of the function on Sn
	for i = 1:24
		VA[i] = Float64(i)
	end

	SNF = snf(4,PA,VA) #Put the values into the order used to compute the FFT

	println("Demonstrates how to put the values of a function on Sn in the order used to compute the FFT")
	println("")
	println("Start with any ordering: ")
	println("Permutation   Function value for that Permutation")
	for i = 1:24
		ST = permutation_string(PA[i])
		ST = string(ST, "  ", VA[i])
		println(ST)
	end
	println("")
	println("Is sorted into this order for SnFFT to use: ")
	println("Permutation   Function value for that Permutation")
	for i = 1:24
		ni = round(Int, SNF[i])
		ST = permutation_string(PA[ni])
		ST = string(ST, "  ", VA[ni])
		println(ST)
	end
end

#Demonstrates how to compute the FFT, the IFFT, and shows that the IFFT recovers the initial function
function example4()
	example4(4)
end

# Parameters:
#	N::Int
#	- the problem size
function example4(N::Int)

	YOR, PT = yor(N) #Creates Young's Orthogonal Representations and the Pointer Tree for Partition decomposistion
	SNF = 2 * rand(factorial(N)) #Create a random function on Sn
	FFT = sn_fft(N, SNF, YOR, PT) #Find the FFT
	iSNF = sn_ifft(N, FFT, YOR, PT) #Find the inverse FFT
	dif = SNF - iSNF #Find the error in the process

	P, WI = partitions(N) #Create the Partitions of 1:N
	Pn = P[N] #Get the partitions that act as the labels for the irreducibles

	println("Demonstrates how to compute the FFT, the IFFT, and shows that the IFFT recovers the initial function")
	println("")
	println("A random function on Sn: ")
	for i = 1:length(SNF)
		println(SNF[i])
	end
	println("")
	println("The Fourier Transform: ")
	for i = 1:length(FFT)
		println(partition_string(Pn[i]))
		println(round(FFT[i], 4))
		println("")
	end
	println("Maximum error in recovering the original function on Sn: ")
	println(maximum(abs(dif)))
	end

#Demonstrates how to compute a sparse FFT and shows that it produces the same results as the normal FFT
function example5()
	example5(4,0.8)
end

# Parameters:
#	N::Int
#	- the problem size
#	SC::Float64
#	- the portion of the function that is zero-valued
function example5(N::Int, SC::Float64)

	#create a sparse function on Sn
	SNF = rand(factorial(N))
	L = 0
	for(i = 1:length(SNF))
		if(SNF[i] < SC) #determine if the value will be zero
			SNF[i] = 0.0
		else #if non-zero, assign it a random value
			L += 1
			SNF[i] = 2 * rand()
		end
	end
	sSNF = Array(Float64, L)
	NZL = Array(Int, L)
	index = 1
	for i = 1:length(SNF)
		if SNF[i] != 0.0
			sSNF[index] = SNF[i]
			NZL[index] = i
			index += 1
		end
	end

	YOR, PT = yor(N) #Creates Young's Orthogonal Representations and the Pointer Tree for Partition decomposistion

	nFFT = sn_fft(N, SNF, YOR, PT) #Compute the FFT normally
	sFFT = sn_fft_sp(N, sSNF, NZL, YOR, PT) #Compute the sparse FFT

	P, WI = partitions(N) #Create the Partitions of 1:N
	Pn = P[N] #Get the partitions that act as the labels for the irreducibles

	println("Demonstrates how to compute a sparse FFT and shows that it produces the same results as the normal FFT")
	println("")
	println("A random sparse function on Sn")
	for i = 1:length(SNF)
		println(SNF[i])
	end
	println("")
	println("The sparse Fourier Transform: ")
	for i = 1:length(sFFT)
		println(partition_string(Pn[i]))
		println(round(sFFT[i], 4))
		println("")
	end
	println("Comparing the results of the two methods of computing the FFT using '==' evaluates to: ")
	println(sFFT == nFFT)
end

#Demonstrates the property of the FFT where the components of the FFT are the representations when the FFT is computed on a delta function
function example6()
	example6(5,[1,3,2,4,5])
end

# Parameters:
#	N::Int
#	- the problem size
#	Permutation::Array{Int, 1}
#	- a permutation of N
function example6(N::Int, Permutation::Array{Int, 1})

	YOR, PT = yor(N) #Creates Young's Orthogonal Representations and the Pointer Tree for Partition decomposistion
	YORn = YOR[N]

	#Create a delta function on Sn
	SNF = [1.0]
	NZL = [permutation_index(Permutation)]

	FFT = sn_fft_sp(N, SNF, NZL, YOR, PT) #Find the FFT, because delta functions are almost completely sparse the sparse FFT is used

	error = 0.0
	for i = 1:length(FFT)
		YORp = yor_permutation(Permutation, YORn[i]) #Find the representation for the selected permutation for this partition

		#Find the error
		error_p = maximum(abs(YORp - FFT[i]))
		if error_p > error
			error = error_p
		end
	end

	P, WI = partitions(N) #Create the Partitions of 1:N
	Pn = P[N] #Get the partitions that act as the labels for the irreducibles

	println("Demonstrates the property of the FFT where the components of the FFT are the representations when the FFT is computed on a delta function")
	println("")
	println("The permutation for which the funciton is non zero: ", permutation_string(Permutation))
	println("")
	println("The Fourier transform: ")
	for i = 1:length(FFT)
		println(partition_string(Pn[i]))
		println(round(FFT[i], 4))
		println("")
	end
	println("Maximum absolute difference between the Fourier components and the YOR representations for the permutation: ")
	println(error)
end

#Creates a random Band Limited funtion and returns both a compressed version (that FFT_BL() can use) and a full version (that FFT() can use)
#computing a bandlimited FFT and iFFT
function example7()
	example7(4,2)
end

# Parameters:
#	N::Int
#	- the problem size
#	K::Int
#	- the problem is homogenous at N-K
function example7(N::Int, K::Int)

	blSNF = 2 * rand(round(Int, factorial(N) / factorial(N - K)))
	nSNF = Array(Float64, factorial(N))
	index = 1
	for i = 1:length(blSNF)
		v = blSNF[i]
		for j = 1:factorial(N - K)
			nSNF[index] = v
			index += 1
		end
	end

	YOR, PT = yor(N)  #Creates Young's Orthogonal Representations and the Pointer Tree for Partition decomposistion
	nFFT = sn_fft(N, nSNF, YOR, PT) #compute the normal FFT

	YOR, PT, ZFI = yor_bl(N, K)  #Creates Young's Orthogonal Representations and the Pointer Tree for Partition decomposistion for bandlimited functions
	blFFT = sn_fft_bl(N, K, blSNF, YOR, PT, ZFI) #compute the band limited FFT

	iblSNF = sn_ifft_bl(N, K, blFFT, YOR, PT, ZFI) #compute the band limited iFFT

	P, ZFI, WI = partitions_bl(N,K) #Create the Partitions of 1:N
	Pn = P[N] #Get the partitions that act as the labels for the irreducibles

	error = 0.0
	shift = length(nFFT) - length(blFFT)
	for i = 1:length(blFFT)
		#Find the error
		error_p = maximum(abs(nFFT[i + shift] - blFFT[i]))
		if error_p > error
			error = error_p
		end
	end

	println("Creates a random Band Limited funtion and returns both a compressed version (that FFT_BL() can use) and a full version (that FFT() can use) computing a bandlimited FFT and iFFT")
	println("")
	println("The full representation of a bandlimited function on Sn:")
	for i = 1:length(nSNF)
		println(nSNF[i])
	end
	println("")
	println("The compressed representation of the same function used by SnFFT:")
	for i = 1:length(blSNF)
		println(blSNF[i])
	end
	println("")
	println("The band limited Fourier transform: ")
	for i = 1:length(blFFT)
		println(partition_string(Pn[i]))
		println(round(blFFT[i], 4))
		println("")
	end

	println("The maximum absolute difference between the full and bandlimited FFT: ")
	println(error)
	println("")
	println("The absolute error in recovering the band limited function on Sn from the inverse bandlimited FFT: ")
	println(maximum(abs(iblSNF - blSNF)))
end

#Calculating a partial iFFT
function example8()
	example8(4,3)
end

# Parameters:
#	N::Int
#	- the problem size
#	M::Int
#	- the number of the top frequencies of the Fourier transform to use
function example8(N::Int, M::Int)

	SNF = rand(factorial(N)) #create a random function on Sn
	YOR, PT = yor(N) #create YOR matrics and the Pointer Tree
	FFT = sn_fft(N, SNF, YOR, PT) #calculate the normal FFT (the sparse FFT also works)
	piSNF = sn_ifft_p(N, M, FFT, YOR, PT) #calculate the partial iFFT

	println("Demonstrates how to calculate a partial iFFT")
	println("")
	println("The original random function on Sn:")
	for i = 1:length(SNF)
		println(SNF[i])
	end
	println("")
	println("The recovered function from the partial iFFT considering the lowest ", M, " frequencies: ")
	for i = 1:length(piSNF)
		println(piSNF[i])
	end
end

#Uses sparcl (an R library) to perform sparse Hierarchical Clustering
#The output is saved into your home directory
function example_clustering()
	example_clustering(2,5,1.0,7)
end

# Parameters:
#	C::Int
#	- number of clusters
#	S::Int
#	- sizes of the clusters
#	F::Float64
#	- spreading factors for the distributions
#	N::Int
#	- the problem size
# Return Values:
# None - saves the data to CSV files and runs the clustering script
# Notes: For Julia V0.4
# - the TmStruct and related code can be replaced with "now()"
# - the appropriate version number needs to be used to define $loc
function example_clustering(C::Int,S::Int,F::Float64,N::Int)

	mkdir = "mkdir"
	Time = now()
	dir = string("ClusteringExample-",Time)
	R = "R"
	CMD = "CMD"
	BATCH = "BATCH"
	loc = string(homedir(),"/.julia/v0.4/SnFFT/src/sparcl_script.R")

	run(`$mkdir $dir`)
	cd(dir)
	clustering_setup(C,S,F,N)
	run(`$R $CMD $BATCH $loc`)
end

#########################################
# Functionality needed for the examples #
#########################################

###
# ClusteringExample
###

#Sets up the data for the clustering
function clustering_setup(C::Int64, S::Int, F::Float64, N::Int64)

	#Construct the distribution Centers
	Centers = Array(Array{Int64, 1}, C)
	for i = 1:C
		Centers[i] = sn_p(N)
		println(Centers[i])
	end

	#Construct the Distributions
	Permutations = Array(Array{Int64, 1}, C * S)
	Distributions = Array(Array{Float64, 1}, C * S)
	k = 1
	for i = 1:C
		center = Centers[i]
		for j = 1:S
			rt = sn_t(N)
			perm = sn_multiply(rt, center)
			Permutations[k] = perm
			Distributions[k] = mallowsdistribution(perm, F)
			println("Distribution ", k)
			k += 1
		end
	end

	#Calculate the Fourier Transforms
	FTA = Array(Array{Array{Float64, 2}, 1}, C * S)
	YOR, PT = yor(N)
	for i = 1:(C * S)
		FTA[i] = sn_fft(N, Distributions[i], YOR, PT)
		println("FFT ", i)
	end

	X = Array(Float64, C * S, factorial(N))
	for c = 1:(C * S)
		k = 1
		FT = FTA[c]
		for p = 1:length(FT)
			FC = FT[p]
			d = size(FC, 1)
			for i = 1:d
				for j = 1:d
					X[c, k] = FC[i,j]
					k += 1
				end
			end
		end
	end

	writecsv("centers.csv",Centers)
	writecsv("perms.csv",Permutations)
	writecsv("input-x.csv",X)
end

# Parameters:
#	P::Array{Int, 1}
#	- a permutation
#	Gamma::Float64
#	- the spreading factor
# Return Values:
#	MD::Array{Float64, 1}
#	- the mallows distribution with spreading factor Gamma centered at P
function mallowsdistribution(P::Array{Int64, 1}, Gamma::Float64)
	N = length(P)
	N_F = factorial(N)
	Q = preferencematrix(P)
	MD = Array(Float64, N_F)
	Gamma_neg = -1 * Gamma
	V = (1 - exp(Gamma_neg)) ^ N
	for j = 1:N
		V /= (1 - exp(j * Gamma_neg))
	end
	for i = 1:N_F
		Pi = index_permutation(N, i)
		Qi = preferencematrix(Pi) #Find the Preference Matrix of that permutation
		MD[i] = V * exp(Gamma_neg * kendalldistance(Qi, Q)) #Find the value of the Mallow Distribution for that permutation
	end
	return MD
end

# Parameters:
#	P::Array{Int, 1}
#	- a permutation
# Return Values:
#	Q::Array{Int, 2}
#	- the preference matrix for P
#	- Q[i,j] = 1 if and only if j precedes i in P
function preferencematrix(P::Array{Int64, 1})
	N = length(P)
	Q = zeros(Int64, N, N)
	for i = 1:N
		for j = 1:N
			p = P[j]
			if p == i
				break
			else
				Q[i, p] = 1
			end
		end
	end
	return Q
end

# Parameters:
#	Q1::Array{Int, 2}
#	- the preference matrix for the first permutation
#	Q2::Array{Int, 2}
#	- the preference matrix for the second permutation
# Return Values:
#	D::Int
#	- the Kendall Tau Distance between two permutations
#	- the the number of pairs (i, j) such that: P1[i] < P1[j] and P2[i] > P2[j]
function kendalldistance(Q1::Array{Int64, 2}, Q2::Array{Int64, 2})
	N = size(Q1, 1)
	D = 0
	for i = 1:N
		for j = 1:N
			D += abs(Q1[i,j] - Q2[i,j])
		end
	end
	D /= 2
	D = round(Int, D)
	return D
end

###
# Printing methods
###

# Parameters:
#	Permutation::Array{Int, 1}
#	- a permutation
# Return Values:
#	ST::String
#	- the string representation of permutation
function permutation_string(Permutation::Array{Int, 1})
	ST = "[ "
	for i = 1:length(Permutation)
		ST = string(ST, Permutation[i], " ")
	end
	ST = string(ST, "] ")
	return ST
end

# Parameters:
#	Partition::Array{Int, 1}
#	- a partition
# Return Values:
#	ST::String
#	- the string representation of partition
function partition_string(Partition::Array{Int, 1})
	ST = "( "
	for i = 1:length(Partition)
		ST = string(ST, Partition[i], " ")
	end
	ST = string(ST, ") ")
	return ST
end

# Parameters:
#	YS::Array{Int, 1}
#	- a yamanouchi symbol
# Return Values:
#	None - prints YS
function print_ys(YS::Array{Int, 1})
	N = length(YS)
	numRows = maximum(YS)
	ST = Array(AbstractString, numRows)
	IA = ones(Int, numRows)
	for i = 1:numRows
		ST[i] = ""
	end
	for n = 1:N
		row = YS[n]
		NI = IA[row]
		ST[row] = string(ST[row]," ", n)
		NI += 1
		IA[row] = NI
	end
	for i = 1:numRows
		println(ST[i])
	end
end


text:
# using Revise, Test
using BifurcationKit, LinearAlgebra, Setfield
const BK = BifurcationKit

function test_newton(x0)
	Ty = eltype(x0)

	F(x, p) = x.^3 .- 1
	Jac(x, p) = diagm(0 => 3 .* x.^2)

	opts = NewtonPar( tol = Ty(1e-8), verbose = true)
	sol, hist, flag, _ = newton(F, Jac, x0, nothing, opts, normN = x->norm(x,Inf), callback = BK.cbMaxNorm(100.0))
end
####################################################################################################
# we test the regular newton algorithm
# simple case
test_newton(ones(10) .+ rand(10) * 0.1)

# test types for newton
sol, = test_newton(Float16.(ones(10) .+ rand(10) * 0.1))
@test eltype(sol) == Float16
####################################################################################################
function test_newton_palc(x0, p0)
	Ty = eltype(x0)
	N = length(x0)

	θ = Ty(0.2)
	dotθ = BK.DotTheta(dot)

	F(x, p) = x.^3 .- p
	Jac(x, p) = diagm(0 => 3 .* x.^2)

	z0 = BorderedArray(x0, p0)
	τ0 = BorderedArray(rand(Ty, N), convert(typeof(p0), 0.2))
	zpred = BorderedArray(x0, convert(typeof(p0), 0.3))
	optn = NewtonPar{Ty, DefaultLS, DefaultEig}()
	optc = ContinuationPar{Ty, DefaultLS, DefaultEig}(newtonOptions = optn)
	sol, hist, flag, _ = newtonPALC(F, Jac, p0, (@lens _), z0, τ0, zpred, Ty(0.02), θ, optc, dotθ)
end

sol, = test_newton_palc(ones(10) .+ rand(10) * 0.1, 1.)

#test type
sol, = test_newton_palc(Float32.(ones(10) .+ rand(10) * 0.1), Float32(1.))
@test typeof(sol) == BorderedArray{Vector{Float32}, Float32}
####################################################################################################
# test of  deflated problems
_T = Float32
F4def = (x, p) -> @. (x-1) * (x-2)
J4def = (x, p) -> BK.finiteDifferences(z->F4def(z,p), x)
deflationOp = DeflationOperator(_T(2), dot, _T(1), [[_T(1)]])
@test firstindex(deflationOp) == 1
@test lastindex(deflationOp) == 1
@test eltype(deflationOp) == _T
@test deflationOp(rand(_T,1)) isa _T
defpb = DeflatedProblem(F4def, J4def, deflationOp)
@test defpb(rand(_T, 1), nothing) |> eltype == _T
@test defpb(rand(_T, 1), nothing, rand(_T, 1)) |> eltype == _T

push!(deflationOp, rand(_T,1))
@test deflationOp(zeros(_T, 1)) isa _T
@test deflationOp(rand(_T, 1), rand(_T, 1)) isa _T
copy(deflationOp)

#  test custom distance
deflationOp2 = DeflationOperator(_T(2), BifurcationKit.CustomDist((u,v)->norm(u-v)), _T(1), deflationOp.roots)
@test deflationOp2(zeros(_T, 1)) isa _T
length(deflationOp2)
deflationOp2(rand(_T, 1), rand(_T, 1))

ext: jl


text:
#############################################################################
#Copyright (c) 2016-2017 Charles Le Losq
#
#The MIT License (MIT)
#
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the #Software without restriction, including without limitation the rights to use, copy, #modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, #and to permit persons to whom the Software is furnished to do so, subject to the #following conditions:
#
#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, #INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
#
#############################################################################

"""
	rameau(paths::Tuple,switches::Tuple;input_properties=('\t',0),prediction_coef=[0.0059;0.0005],temperature=23.0,laser=532.0,lb_break=1600.,hb_start=2600.,roi_hf_external = [3000. 3100.; 3800. 3900.],basetype="gcvspline",mmap_switch=true)

	The RAMEAU function in Spectra is mostly experimental. For a robust use of internal calibrations, please see the rampy.rameau() class built in Python which is more up to date.

INPUTS:

    paths: Tuple{Strings}, it contains the following strings:

		in_liste: the relative path of your liste of spectra (see the description of this liste below), e.g. "./liste_2012.csv"

		in_path = the relative path of the repertory where your raw spectra are stored, e.g. "./raw/"

		out_path = the relative path of the repertory where you want to output the corrected spectra, e.g. "./treated/" (THIS FOLDER MUST EXIST PRIOR TO RUN)

		fig_path= the relative path of the repertory for saving the figures, e.g. "./figures/"

		rws_save_file = the relative path of the csv file where the results will be stored, e.g. "rws_calculated.csv". If calibration is set to "yes", it outputs an array containing the Rws, the provided water content and the calculated water content in columns 1, 2 and 3, respectively. If calibration is set to anything else, it outputs the rws and the water content predicted with the provided prediction_coef

		rws_save_fig = the relative path of the figure for the calibration, e.g. "rws_calculated.jpg". Note: It only works when calibration is on "yes", and Type_of_calibration is set on "internal". Just set it to "trash.jpg" if you are not going to use it.

	input_properties: Tuple, it contains the delimiter of the columns in your spectra files, and the number of starting lines to skip. For instance, use ('\t',1) if your files have one header line and the columns are separated  by tabulations, or (',',0) if you use CSV files without header.

	switches: Tuple, it contains the different switches as (Type_of_calibration, calibration?, experimental?, temperature_laser_correction?).

		Type_of_calibration: should be set to "internal" or "external". The external mode requires standards, as described in Behrens et al. (2006) and Thomas et al. (2008). If you use the external mode, you can leave blanks for the other switches. For instance, enter ("external","","","").

		calibration? :  use only with the "internal" mode, this should be equal to "yes" or "no". Enter yes if you are setting up a new calibration, or "no" if you want to use a previously determined calibration coefficient with providing the later in the prediction_coef variable (see options).

		experimental? : use only with the "internal" mode, this is an experimental code. For now, only the "double" feature is advised. It allows using different smoothing spline coefficients for the water and silicate bands, delimited by the lb_break and hb_start variables (see below).

		temperature_laser_correction? : use only with the "internal" mode, this should be equal to "yes" or "no". This asks if you want to use the temperature-laser wavelength correction as done in Le Losq et al. (2012). If you use the double baseline mode, this correction is applied after removing the background under the water peak.

OPTIONS:

	prediction_coef: Array{Float64}, this array contains the calibration coefficient with its error bar. Those values will be used in the predictions, if you use the predictive mode (i.g., calibration switch is not set to "yes"). Default = [0.0059;0.0005].

	temperature: Float64, the temperature for the temperature-laser wavelength correction in Celsius. Default = 23.0.

	laser: Float64, the laser wavelength in nm for the temperature-laser wavelength correction. Default = 532.0.

	lb_break: Float64, for double baseline correction, the breaking point before which the software will consider the BIRs in the low frequency region. Default = 2010.0.

	hb_start: Float64, for double baseline correction, the breaking point after which the software will consider the BIRs in the high frequency region. Default = 1000.0.

  roi_hf_external: Array{Float64}, the roi for fitting the linear baseline in the external calibration mode. Default = [3000. 3100.; 3800. 3900.].

	basetype: String, the type of baseline you want to fit. Corresponds to the "basetype" parameter of the baseline function.  Default = "gcvspline".

    mmap_switch: false or true, this allows to switch on or off the memory mapping in the `readcsv`/`readdlm` functions that `rameau` uses. Default = "true".

OUTPUTS:

	Rameau does not provide any outputs directly in Julia, but saves everything in the folders you indicate in the variable "paths".

"""
function rameau(paths::Tuple,switches::Tuple;input_properties=('\t',0),prediction_coef=[0.0062;0.0005],temperature=23.0,laser=532.0,lb_break=1600.,hb_start=2600.,roi_hf_external = [3000. 3100.; 3800. 3900.],basetype="gcvspline",mmap_switch=true,lambda=10.0^12,p=0.01)

	# some function definition
	calibration_model(x, p) = p[1].*x
	baseline_lfmodel(x, p) = p[1].*x + p[2].*x.^2 + p[3].*x.^3

	if switches[1] == "internal"

		scale = 1 # no scaling, keeping that for history

		# reading the list of data to get the names, smoothing factors and roi
		liste=readcsv(paths[1], skipstart=1,use_mmap = mmap_switch)
		names = liste[:,1]
		water = Array{Float64}(liste[:,3])
		smo_lf = Array{Float64}(liste[:,4])
		smo_hf = Array{Float64}(liste[:,5])
		roi = ones(Int((size(liste,2)-5)/2),2,size(liste,1))
		for j = 1:size(liste,1)
			roi[:,1,j] = (liste[j,6:2:end])
			roi[:,2,j] = (liste[j,7:2:end])
		end

		rws = ones(size(liste,1),3) # output array of ratios

		for i = 1:size(liste,1) # starting a loop other the data list\

			# reading the spectra
			spectra = Array{Float64}(readdlm(string(paths[2],names[i]),input_properties[1],skipstart=input_properties[2],use_mmap = mmap_switch))

			# Making sure that the frequencies are in increasing order
			if spectra[end,1] < spectra[1,1]
				spectra = flipdim(spectra,1)
			end

			# Assigning things, primary treatments
			x = spectra[:,1]
			spectra[:,2] = spectra[:,2] - minimum(spectra[:,2]) + 0.01 #we avoid any 0 values by putting the smallest value to 0.01
			y = spectra[:,2]./maximum(spectra[:,2]).*scale # intensity normalisation
			ratio_bkg = minimum(y)/maximum(y) # to keep a record of the ratio of maximum signal intensity over minimum background intensity

			#### PRELIMINARY STEP: FIRST WE GRAB THE GOOD SIGNAL IN THE ROI
			interest_index::Array{Int64} = find(roi[1,1,i] .<= x[:,1] .<= roi[1,2,i])
			if size(roi)[1] > 1
				for j = 2:size(roi)[1]
					interest_index = vcat(interest_index,  find(roi[j,1,i] .<= x[:,1] .<= roi[j,2,i]))
				end
			end
			interest_x = x[interest_index,1]
			interest_y = y[interest_index,1]

			if switches[3] == "double" # if user asks for using two different smoothing factors for the silicate and water regions

				if switches[4] == "yes"

					# If the Long correction is asked, then it will use the two smoothing spline factors and perform a "traditional" internal baseline fit with gcvspline
					x, y_long, ~ = tlcorrection([x[:] y[:]],temperature,laser) # Long correction
					y_long=y_long./maximum(y_long).*scale

					lb = roi[:,1,i]
					hb = roi[:,2,i]

					roi_hf =  [lb[hb.>hb_start] hb[hb.>hb_start]]
					roi_lf =  [lb[hb.<lb_break] hb[hb.<lb_break]]

					y_hf = y_long[x.>hb_start]
					y_lf = y_long[x.<lb_break]

					# A first baseline in signals between ~1300 and ~ 2000 cm-1
					y_calc_hf, baseline_hf = baseline(x,y_long,roi_hf,basetype,p=smo_hf[i])
					y_calc_lf, baseline_lf = baseline(x,y_long,roi_lf,basetype,p=smo_lf[i])

					# putting back the signals together
					y_calc2 = [y_calc_lf[x.<=lb_break];y_calc_hf[x.>lb_break]]
					bas2 = [baseline_lf[x.<=lb_break];baseline_hf[x.>lb_break]]

					# making a nice figure
					figure(figsize=(20,20))
					suptitle("LL2012 method, sp. $(names[i])")

					subplot(3,2,(1,2))
					plot(x,y,"black",label="Raw sp.")
					plot(x,y_long,"blue",label="Long corr. sp.")
					scatter(interest_x,y_long[interest_index,1],s=10.,color="red",label="ROI") # the ROI signals after first baseline and long correction
					plot(x,bas2,"red",label="Baseline")
					legend(loc="best",fancybox="true")

					xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
					ylabel("Intensity,a. u.",fontsize=18,fontname="Arial")

					subplot(3,2,3)
					plot(x[x.<1500],y_long[x.<1500],"blue",label="Long corr. sp.")
					scatter(interest_x[interest_x.<1500],y_long[interest_index[interest_x.<1500],1],s=10.,color="red",label="ROI") # the ROI signals after first baseline and long correction
					plot(x[x.<1500],bas2[x.<1500],"red",label="Baseline")
					xlim(0,1500)
					legend(loc="best",fancybox="true")
					xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
					ylabel("Intensity,a. u.",fontsize=18,fontname="Arial")

					subplot(3,2,4)
					plot(x[x.>2800],y_long[x.>2800],"blue",label="Long corr. sp.")
					scatter(interest_x[interest_x.>2800],y_long[interest_index[interest_x.>2800],1],s=10.,color="red",label="ROI") # the ROI signals after first baseline and long correction
					plot(x[x.>2800],bas2[x.>2800],"red",label="Baseline")
					xlim(2600,4000)
					legend(loc="best",fancybox="true")
					xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
					ylabel("Intensity,a. u.",fontsize=18,fontname="Arial")

					subplot(3,2,(5,6))
					plot(x,y_calc2,"cyan",label="Final sp.")
					xlabel(L"Raman shift, cm$^{-1}$")
					ylabel("Intensity, area normalised")
					legend(loc="best",fancybox="true")
					xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
					ylabel("Intensity,a. u.",fontsize=18,fontname="Arial")

					tight_layout()

				else
					# ALSO FULLY EXPERIMENTAL!!!! IT DOES NOT WORK!
					y_hf = y[x.>hb_start]
					y_lf = y[x.<lb_break]

					# A first linear baseline in signals between ~1300 and ~ 2000 cm-1
					y_calc_hf, baseline_hf = baseline(x,y,[2600. 2900.;3775. 3900.],"poly",p=3.0)

				    #A first linear baseline in signals between ~1300 and ~ 2000 cm-1
				   	y_calc1_part1, baseline1 = baseline(x,y,[1300 1800.],"poly",p=1.0)

					# we get a mean signal for background below 1300 cm-1
					X_low_lf = x[y .== minimum(y[100 .<x.<1400])]
					b_low_lf = mean(y[X_low_lf[1]-10 .< x .< X_low_lf[1]+10])

					baseline1[x.<1250.] = b_low_lf
					y_calc1_part1 = y - baseline1

					bas2 = [baseline1[x.<lb_break];baseline_hf[x.>lb_break]]
					y_calc2 = [y_calc1_part1[x.<lb_break];y_calc_hf[x.>lb_break]]

					figure()
					plot(x,y,color="black")
					plot(x,baseline_hf,color="red")
					plot(x,y_calc2,color="blue")
					tight_layout()
				end

			else # only a single baseline treatment is asked

				# we test if the Long correction is asked, and react accordingly for baseline substraction
				if switches[4] == "yes"
					x, y_long, ~ = tlcorrection([x[:] y[:]],temperature,laser) # Long correction
					y_long=y_long./maximum(y_long).*scale

					if basetype == "arPLS" || basetype == "whittaker"
						y_calc2, bas2 = baseline(x,y_long,roi[:,:,i],basetype,lambda = lambda, p = p)
					else
						y_calc2, bas2 = baseline(x,y_long,roi[:,:,i],basetype,p=smo_lf[i])
					end

					y_calc2 = y_calc2[:,1]./trapz(x[:],y_calc2[:,1]).*(scale./10) # area normalisation

					figure(figsize=(20,20))
					suptitle("LL2012 method, sp. $(names[i])")

					subplot(3,2,(1,2))
					plot(x,y,"black",label="Raw sp.")
					plot(x,y_long,"blue",label="Long corr. sp.")
					scatter(interest_x,y_long[interest_index,1],s=10.,color="red",label="ROI") # the ROI signals after first baseline and long correction
					plot(x,bas2,"red",label="Baseline")
					legend(loc="best",fancybox="true")

					xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
					ylabel("Intensity,a. u.",fontsize=18,fontname="Arial")

					subplot(3,2,3)
					plot(x[x.<1500],y_long[x.<1500],"blue",label="Long corr. sp.")
					scatter(interest_x[interest_x.<1500],y_long[interest_index[interest_x.<1500],1],s=10.,color="red",label="ROI") # the ROI signals after first baseline and long correction
					plot(x[x.<1500],bas2[x.<1500],"red",label="Baseline")
					xlim(0,1500)
					legend(loc="best",fancybox="true")
					xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
					ylabel("Intensity,a. u.",fontsize=18,fontname="Arial")

					subplot(3,2,4)
					plot(x[x.>2800],y_long[x.>2800],"blue",label="Long corr. sp.")
					scatter(interest_x[interest_x.>2800],y_long[interest_index[interest_x.>2800],1],s=10.,color="red",label="ROI") # the ROI signals after first baseline and long correction
					plot(x[x.>2800],bas2[x.>2800],"red",label="Baseline")
					xlim(2600,4000)
					legend(loc="best",fancybox="true")
					xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
					ylabel("Intensity,a. u.",fontsize=18,fontname="Arial")

					subplot(3,2,(5,6))
					plot(x,y_calc2,"cyan",label="Final sp.")
					xlabel(L"Raman shift, cm$^{-1}$")
					ylabel("Intensity, area normalised")
					legend(loc="best",fancybox="true")
					xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
					ylabel("Intensity,a. u.",fontsize=18,fontname="Arial")

					tight_layout()

				else
					if basetype == "arPLS" || basetype == "whittaker"
						y_calc2, bas2 = baseline(x,y_long,roi[:,:,i],basetype,lambda = lambda, p = p)
					else
						y_calc2, bas2 = baseline(x,y_long,roi[:,:,i],basetype,p=smo_lf[i])
					end

					y_calc2 = y_calc2[:,1]./trapz(x[:],y_calc2[:,1]).*scale # area normalisation

					figure()
					plot(x,y.*10,"black",label="Raw sp.")
					scatter(interest_x,interest_y.*10,s=10.,color="red",label="ROI")
					plot(x,bas2.*10,"red",label="Baseline")
					plot(x,y_calc2,"cyan",label="Final sp.")
					xlabel(L"Raman shift, cm$^{-1}$")
					ylabel("Intensity, area normalised")
					title("Baseline sub., NO Long corr.")
					legend(loc="best",fancybox="true")
					tight_layout()
				end

			end

			# saving the corrected spectrum
			writecsv(string(paths[3],names[i]),[x y_calc2 bas2])

			#Saving the generated figure (specified directory should exist)
			savefig(string(paths[4],"baseline_",names[i],".pdf"))
    		close()
			if switches[3] == "double"
				# calc=ulating the areas under silicate and water bands
				As = trapz(x[150 .< x .<1300],y_calc2[150 .< x .<1300])
				Aw = trapz(x[2800 .< x .<3800],y_calc2[2800 .< x .<3800])
			else
				# calc=ulating the areas under silicate and water bands
				As = trapz(x[150 .< x .<1300],y_calc2[150 .< x .<1300])
				Aw = trapz(x[2800 .< x .<3800],y_calc2[2800 .< x .<3800])
			end

			# recording them
			rws[i,1] = As;
			rws[i,2] = Aw;
			rws[i,3] = Aw./(As);
			println("Spectrum $(names[i]), the rws is $(rws[i,3])")
		end

		if switches[2] == "yes"
			fit = curve_fit(calibration_model, rws[:,3], water[:]./(100.0 -water[:]), [0.5])
			coef = fit.param
			sigma = estimate_errors(fit, 0.95)

			rws_calibration = collect(0:0.01:round(maximum(rws[:,3]),2))
			water_ratio_calibration = calibration_model(rws_calibration,coef)
			water_compare =  100.0 .*calibration_model(rws[:,3],coef)./(calibration_model(rws[:,3],coef)+1.0) # eq. 3 Le Losq et al. (2012)

			rmse_calibration = sqrt(1.0 ./(size(rws,1)-1.0).*sum((water_compare-water).^2))
			rmse_ratio_ws = rmse_calibration./100.0 .*(rmse_calibration./(100.0-rmse_calibration)+1.0)

			figure()
			scatter(rws[:,3],water[:]./(100.0 -water[:]))
			plot(rws_calibration,water_ratio_calibration,color="red",linewidth=2.0)
			plot(rws_calibration,water_ratio_calibration+rmse_ratio_ws,color="red",linestyle="--",linewidth=1.0)
			plot(rws_calibration,water_ratio_calibration-rmse_ratio_ws,color="red",linestyle="--",linewidth=1.0)
			plot(rws_calibration,water_ratio_calibration,color="red")
			xlim(0,maximum(rws_calibration)+1.0 ./4.0 .*maximum(rws_calibration))
			ylim(0,maximum(water_ratio_calibration)+1.0 ./4.0 .*maximum(water_ratio_calibration))

			xlabel(L"A$_{water}$/A$_{silicates}$, area ratio",fontsize=18,fontname="Arial")
			ylabel("Water/Glass, weight ratio",fontsize=18,fontname="Arial")
			annotate(L"R$_{ws}$=",xy=(0.3,0.9),xycoords="axes fraction",fontsize=18,fontname="Arial",horizontalalignment="center")
			annotate("$(round(coef[1],5)) +/- $(round(sigma[1],5))",xy=(0.3,0.8),xycoords="axes fraction",fontsize=18,fontname="Arial",horizontalalignment="center")
			annotate("Standard deviation\n= $(round(rmse_calibration,2)) wt%",xy=(0.7,0.3),xycoords="axes fraction",horizontalalignment="center",fontsize=18,fontname="Arial")
			tight_layout()
			savefig(paths[6])

			writecsv(string(paths[5]), ["spectrum" "product" "Rws" "Water input" "Water Raman";liste[:,1] liste[:,2] rws[:,3] water water_compare])
			#return [rws[:,3] water water_compare]
		else
			water_predicted = 100.0 .*(rws[:,3].*prediction_coef[1]./(rws[:,3].*prediction_coef[1]+1.0))
			water_predicted_high = 100.0 .*(rws[:,3].*sum(prediction_coef)./(rws[:,3].*sum(prediction_coef)+1.0))
			writecsv(string(paths[5]), ["spectrum" "product" "Rws" "Water Raman" "Error";liste[:,1] liste[:,2] rws[:,3] water_predicted water_predicted_high-water_predicted])
			#return [rws[:,3] water_predicted]
		end

	elseif switches[1] == "external"
		liste=readcsv(paths[1], skipstart=1,use_mmap = mmap_switch)

		reference = liste[:,1]
		# reference product names are in first column for record
		reference_density = liste[:,3]
		reference_water = liste[:,4]

		sample = liste[:,5]
		# sample product names are in column 6 for record
		sample_density = liste[:,7]

		water = ones(size(liste,1),1) # the array to record water content

		for i = 1:size(liste,1)

			#importing data
			reference_sp = readdlm(string(paths[2],reference[i]),input_properties[1],skipstart=input_properties[2],use_mmap = mmap_switch)
			sample_sp = readdlm(string(paths[2],sample[i]),input_properties[1],skipstart=input_properties[2],use_mmap = mmap_switch)

			if reference_sp[end,1] < reference_sp[1,1]
				reference_sp = flipdim(reference_sp,1)
			end

			if sample_sp[end,1] < sample_sp[1,1]
				sample_sp = flipdim(sample_sp,1)
			end

			# temperature - excitation line correction with the Hehlen formula
			x_ref, reference_sp_long, ~ = tlcorrection(reference_sp,temperature,laser,correction="hehlen",normalisation="no",density=reference_density[i])
			x_sp, sample_sp_long, ~ = tlcorrection(sample_sp,temperature,laser,correction="hehlen",normalisation="no",density=sample_density[i])

			# Linear baseline subtraction
			reference_sp_corr2, reference_baseline = baseline(x_ref,reference_sp_long,roi_hf_external,"poly",p=1.0)
			sample_sp_corr2, sample_baseline = baseline(x_sp,sample_sp_long,roi_hf_external,"poly",p=1.0)

			# Area calculation between 3100 and 3800 cm-1
			Area_reference = trapz(x_ref[roi_hf_external[1,2] .<x_ref[:,1].<roi_hf_external[2,1]],reference_sp_corr2[roi_hf_external[1,2] .<x_ref[:,1].<roi_hf_external[2,1]])
			Area_sample = trapz(x_sp[roi_hf_external[1,2] .<x_sp[:,1].<roi_hf_external[2,1]],sample_sp_corr2[roi_hf_external[1,2] .<x_sp[:,1].<roi_hf_external[2,1]])

			# With the intensity
			Intensity_reference = maximum(reference_sp_corr2[roi_hf_external[1,2] .<x_ref[:,1].<roi_hf_external[2,1]])
			Intensity_sample = maximum(sample_sp_corr2[roi_hf_external[1,2] .<x_ref[:,1].<roi_hf_external[2,1]])

			# Water calculation
			water[i] = (reference_water[i]./1.8*reference_density[i]) .*Area_sample./Area_reference # water is in mol/L
			#water[i] = (reference_water[i]./1.8*reference_density[i]) .*Intensity_sample./Intensity_reference # water is in mol/L

			water[i] = water[i] *1.8./sample_density[i] # water converted in wt%

			figure(figsize=(7.5,15))

			subplot(311)
			title("External calibration: $(round(water[i],2)) wt% water,\n sp. $(sample[i]) with ref. $(reference[i]) ")
			plot(reference_sp[:,1],reference_sp[:,2],color="black",label="Reference")
			plot(sample_sp[:,1],sample_sp[:,2],color="blue",label="Sample")
			xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
			legend(loc="best",frameon=false)

			subplot(312)
			plot(x_ref[:],reference_sp_long[:],color="black",label=L"T-$\nu$ corrected reference")
			plot(x_ref[:],reference_baseline,color="grey",label="Reference baseline")
			plot(x_sp[:],sample_sp_long[:],color="blue",label=L"T-$\nu$ corrected sample")
			plot(x_sp[:],sample_baseline,color="cyan",label="Sample baseline")
			xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
			ylabel("Normalized intensity, a. u.",fontsize=18,fontname="Arial")
			legend(loc="best",frameon=false)

			subplot(313)
			plot(x_ref,reference_sp_corr2,color="black",label="Ref. spectrum")
			plot(x_sp,sample_sp_corr2,color="blue",label="Sample spectrum")
			xlabel(L"Raman shift, cm$^{-1}$",fontsize=18,fontname="Arial")
			legend(loc="best",frameon=false)
			tight_layout()
			# saving the water content
			writecsv(string(paths[5]), ["Spectrum" "Sample Name" "Estimated water content";liste[:,5] liste[:,6] water])

			#Saving the generated figure (specified directory should exist)
			savefig(string(paths[4],"External_Std_",sample[i],".pdf"))

		end

	else # or error message
		error("The first switch indicates if you want to use an internal or external calibration mode. Please choose between both.")
	end
end


text:
function get_clusbools(mass_of_average_cell,thresh, max_cats::AbstractArray{T,4}, den) where T
    """
    This function takes in a physical threshold float value and the
    output of the maximum_signature function (4D-Array of signatures)
    to determine the cluster boolean filter.
    """

    # step 1. create a bool filter that is only true when
    # cluster signatures are larger than threshold
    test_bools = (max_cats[:,:,:,1] .> thresh)

    nx, ny, nz = size(max_cats)[1:3]

    # step 2. tag components
    components = Images.label_components(test_bools)
    max_component = maximum(components)

    print("Components: $(max_component)\n")

    # step 3. loop through each component, assess average overdensity
    masses = zeros(T, max_component)
    volume = zeros(T, max_component)

    compare = zeros(nx,nx,nx);

    for i in 1:nx
        for j in 1:nx
            for k in 1:nx
                if components[i,j,k] != 0
                    masses[ components[i,j,k] ] += den[i,j,k]
                    volume[ components[i,j,k] ] += 1.0
                end
            end
        end
    end

    compbool = (masses .> 5e3 / mass_of_average_cell)  .& (volume .> 1)
    println( "viri frac: ", sum( (volume .> 1) .& (masses ./ volume .> 370) ) / max_component)
    println( "mass frac: ", sum(masses[compbool])   / sum(den) )

    clusbool = zeros(Bool, nx, nx, nx)
    for i in 1:nx
        for j in 1:ny
            for k in 1:nz
                if components[i,j,k] != 0
                    if compbool[components[i,j,k]]
                        clusbool[i,j,k] = true
                    end
                end
            end
        end
    end

    clusbool
end

function get_sig_plot(sig_vec, f)
    """
    a function to calculate |dlog(M)/dlog(𝒮)| = ΔM^2 which is
    known as the mass change with respect to signature. This can be
    calculated for both filaments and clusters.
    """
    log10S = collect(-2:0.2:3) #this just makes an array [-2.0, -1.8, -1.6, ... , 2.8, 3.0]
    yy = zeros(Float64, size(log10S))
    for i in 1:size(yy,1)
        yy[i] = sum(
            f[sig_vec .> (10^log10S[i])]
            )
    end

    dydx = abs.( diff(yy.^2) ./ diff( log10S ) );
    midx = (log10S[1:size(log10S,1)-1] .+ log10S[2:size(log10S,1)]) ./ 2;
    return (10 .^ midx), dydx
end

function calc_structure_bools(mass_of_average_cell, cluster_threshold, combined_NEXUS::AbstractArray{T,4}, combined_NEXUSPLUS::AbstractArray{T,4}, den) where T
    """
    A function that returns the boolean filters for each structure based on physical criterion
    """

    N_cells = size(den,1) * size(den,2) * size(den,3)

    ####################################
    #     CLUSTER FILTER PROCESS
    ####################################

    clusbool = get_clusbools(mass_of_average_cell, cluster_threshold, combined_NEXUS, den);

    ####################################
    #    REMAINING FILTER PROCESS
    ####################################

    # we create 3 N^3 arrays with the signature corresponding to each structure
    cluster_signature = reshape(combined_NEXUS[:,:,:,1], N_cells)
    filament_signature = reshape(combined_NEXUSPLUS[:,:,:,2], N_cells)
    wall_signature = reshape(combined_NEXUSPLUS[:,:,:,3], N_cells);

    not_clus_flat = reshape( clusbool .== false, N_cells ); #bool array that is false where there are clusters
    filament_valid = filament_signature[ not_clus_flat ];
    flat_den_valid = reshape(den, N_cells)[ not_clus_flat ];

    S_fil, dM2_fil = get_sig_plot(filament_valid, flat_den_valid)
    max_dM, ind = findmax(dM2_fil) #find the peak of the mass change curve
    filament_thresh = S_fil[ind]

    wall_valid_filt = (not_clus_flat) .& (filament_signature .< filament_thresh)
    wall_valid = wall_signature[wall_valid_filt];
    wall_den_valid = reshape(den, N_cells)[wall_valid_filt];

    S_wall, dM2_wall = get_sig_plot(wall_valid, wall_den_valid)
    max_dM, ind = findmax(dM2_wall)
    wall_thresh = S_wall[ind]

    ####################################
    #        CREATE BOOL FILTERS
    ####################################

    cluster = den[clusbool]
    filbool = (combined_NEXUSPLUS[:,:,:,2] .> filament_thresh) .& (.!clusbool)
    filament = den[filbool]
    wallbool = (combined_NEXUSPLUS[:,:,:,3] .> wall_thresh) .& .!(filbool .| clusbool)
    wall = den[wallbool];

    return clusbool, filbool, wallbool, S_fil, dM2_fil, S_wall, dM2_wall

end


text:
"""
A crude method for finding integer-valued rows of an array
"""
function findintegerfeatures(X::AbstractArray, N=10, prec=3)
    mask = fill(false, size(X, 1))
    for (i, f) ∈ enumerate(eachrow(X))
        σ = std(f)
        fr = round.(f./σ, sigdigits=prec)
        mask[i] = length(unique(fr)) < N
    end
    return mask
end
export findintegerfeatures

_nointegerfeatures(X::AbstractArray, args...) = Y->Y[.!findintegerfeatures(X, args...), :]
nointegerfeatures(X, args...) = _nointegerfeatures(X, args...)(X)
export _nointegerfeatures, nointegerfeatures

"""
Can we remove any perfectly correlated features?
"""
function findredundantfeatures(X::AbstractArray)
    Σ = cor(X, dims=2)
    for r ∈ 1:size(Σ, 1)
        cfs = findall(r == 1.0) # Features perfectly correlated to this one
        cfs = cfs[cfs .!= r] # Don't count the one we're looking at
        for cf ∈ cfs
            Σ[cf, :] .= NaN
            Σ[:, cf] .= NaN
        end
    end
    idxs = isnan.(diag(Σ))
end
export findredundantfeatures

_noredundantfeatures(X::AbstractArray, args...) = Y->Y[.!findredundantfeatures(X, args...), :]
noredundantfeatures(X, args...) = _noredundantfeatures(X, args...)(X)
export _noredundantfeatures, noredundantfeatures


text:
using NumericalDataManipulation.Data.Storage
using NumericalDataManipulation.Data.Manipulation

using Crayons.Box

function merge_two_tables(task::MergeTwoTablesTask)
    @info RED_FG("=======================")
    @info MAGENTA_FG(UNDERLINE(BOLD("[Merge Two Tables Task]")))

    @info RED_FG("-----------------------")
    table_a = load_numerical_table(task.file_a; has_title = task.skip_title_a)
    @info "$(MAGENTA_FG(BOLD("Loaded table A:")))\n$table_a"
    table_b = load_numerical_table(task.file_b; has_title = task.skip_title_b)
    @info "$(MAGENTA_FG(BOLD("Loaded table B:")))\n$table_b"
    @info RED_FG("-----------------------")

    @info RED_FG(".......................")
    grid_a = DataGrid(table_a.argument)
    @info "$(MAGENTA_FG(BOLD("Argument grid of table A:")))\n$grid_a"
    grid_b = DataGrid(table_b.argument)
    @info "$(MAGENTA_FG(BOLD("Argument grid of table B:")))\n$grid_b"

    @info RED_FG("⋯ ⋯ ⋯ ⋯ ⋯")
    tail = task.grid_merge_parameters.tail_knots
    ϵʳᵉˡ = task.grid_merge_parameters.ϵʳᵉˡ
    meanΔ = task.grid_merge_parameters.meanΔ
    merged_grid_result = merge_grids(grid_a, grid_b; ntail = tail, ϵʳᵉˡ = ϵʳᵉˡ, meanΔ = meanΔ)
    @info "$(MAGENTA_FG(BOLD("Merged grid A + B:")))\n$merged_grid_result"
    @info RED_FG("⋯ ⋯ ⋯ ⋯ ⋯")
    @info "$(LIGHT_MAGENTA_FG(BOLD("The master grid is: ")))$(LIGHT_GRAY_FG(task.master == :a ? "A" : "B"))"
    master_grid = task.master == :a ? grid_a : grid_b
    slave_grid =  task.master == :a ? grid_b : grid_a
    merged_grid = merged_grid_result
    @info RED_FG(".......................")

    # ----
    master_table = task.master == :a ? table_a : table_b
    slave_table =  task.master == :a ? table_b : table_a
    target_table = merge_table_slave_to_master_sigmoid(
                        slave_table, master_table,
                        task.columns,
                        map(ivl -> (ivl[1], ivl[2]), task.slave_intervals),
                        (task.master_interval[1], task.master_interval[2]),
                        master_grid, slave_grid, merged_grid)
    # ----

    @info MAGENTA_FG(UNDERLINE(BOLD("[Merge Two Tables Task] - Done")))
    @info RED_FG("=======================")
end


text:
using DiffEqSensitivity, OrdinaryDiffEq, RecursiveArrayTools, DiffEqBase,
      ForwardDiff, Calculus, QuadGK, LinearAlgebra, Zygote
using Test

function fb(du,u,p,t)
  du[1] = dx = p[1]*u[1] - p[2]*u[1]*u[2]*t
  du[2] = dy = -p[3]*u[2] + p[4]*u[1]*u[2]
end
function foop(u,p,t)
  dx = p[1]*u[1] - p[2]*u[1]*u[2]*t
  dy = -p[3]*u[2] + p[4]*u[1]*u[2]
  [dx,dy]
end
function jac(J,u,p,t)
  (x, y, a, b, c) = (u[1], u[2], p[1], p[2], p[3])
  J[1,1] = a + y * b * -1 * t
  J[2,1] = y
  J[1,2] = b * x * -1 * t
  J[2,2] = c * -1 + x
end

f = ODEFunction(fb,jac=jac)
p = [1.5,1.0,3.0,1.0]; u0 = [1.0;1.0]
prob = ODEProblem(f,u0,(0.0,10.0),p)
sol = solve(prob,Tsit5(),abstol=1e-14,reltol=1e-14)
probb = ODEProblem(fb,u0,(0.0,10.0),p)
proboop = ODEProblem(foop,u0,(0.0,10.0),p)

solb = solve(probb,Tsit5(),abstol=1e-14,reltol=1e-14)
sol_end = solve(probb,Tsit5(),abstol=1e-14,reltol=1e-14,
          save_everystep=false,save_start=false)

sol_nodense = solve(probb,Tsit5(),abstol=1e-14,reltol=1e-14,dense=false)
soloop = solve(proboop,Tsit5(),abstol=1e-14,reltol=1e-14)
soloop_nodense = solve(proboop,Tsit5(),abstol=1e-14,reltol=1e-14,dense=false)

# Do a discrete adjoint problem
println("Calculate discrete adjoint sensitivities")
t = 0.0:0.5:10.0
# g(t,u,i) = (1-u)^2/2, L2 away from 1
function dg(out,u,p,t,i)
  (out.=2.0.-u)
end

_,easy_res = adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                                 reltol=1e-14)
_,easy_res2 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                 reltol=1e-14,
                                 sensealg=QuadratureAdjoint(abstol=1e-14,reltol=1e-14))
_,easy_res22 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,
                                  sensealg=QuadratureAdjoint(autojacvec=false,abstol=1e-14,reltol=1e-14))
_,easy_res23 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,
                                  sensealg=QuadratureAdjoint(abstol=1e-14,reltol=1e-14,compile=true))
_,easy_res3 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,
                                  sensealg=InterpolatingAdjoint())
_,easy_res32 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,
                                  sensealg=InterpolatingAdjoint(autojacvec=false))
_,easy_res4 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,
                                  sensealg=BacksolveAdjoint())
_,easy_res42 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                   reltol=1e-14,
                                   sensealg=BacksolveAdjoint(autojacvec=false))
_,easy_res43 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                   reltol=1e-14,
                                   sensealg=BacksolveAdjoint(autojacvec=false,checkpointing=false))
_,easy_res5 = adjoint_sensitivities(sol,Kvaerno5(nlsolve=NLAnderson(), smooth_est=false),
                                 dg,t,abstol=1e-12,
                                 reltol=1e-10,
                                 sensealg=BacksolveAdjoint())
_,easy_res6 = adjoint_sensitivities(sol_nodense,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,
                                  sensealg=InterpolatingAdjoint(checkpointing=true),
                                  checkpoints=sol.t[1:500:end])
_,easy_res62 = adjoint_sensitivities(sol_nodense,Tsit5(),dg,t,abstol=1e-14,
                                   reltol=1e-14,
                                   sensealg=InterpolatingAdjoint(checkpointing=true,autojacvec=false),
                                   checkpoints=sol.t[1:500:end])

# It should automatically be checkpointing since the solution isn't dense
_,easy_res7 = adjoint_sensitivities(sol_nodense,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,
                                  sensealg=InterpolatingAdjoint(),
                                  checkpoints=sol.t[1:500:end])

_,easy_res8 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                    reltol=1e-14,
                                    sensealg=InterpolatingAdjoint(autojacvec=DiffEqSensitivity.TrackerVJP()))
_,easy_res9 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                    reltol=1e-14,
                                    sensealg=InterpolatingAdjoint(autojacvec=DiffEqSensitivity.ZygoteVJP()))
_,easy_res10 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                    reltol=1e-14,
                                    sensealg=InterpolatingAdjoint(autojacvec=DiffEqSensitivity.ReverseDiffVJP())
                                    )
_,easy_res11 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                     reltol=1e-14,
                                     sensealg=InterpolatingAdjoint(autojacvec=DiffEqSensitivity.ReverseDiffVJP(true))
                                     )

adj_prob = ODEAdjointProblem(sol,QuadratureAdjoint(abstol=1e-14,reltol=1e-14),dg,t)
adj_sol = solve(adj_prob,Tsit5(),abstol=1e-14,reltol=1e-14)
integrand = AdjointSensitivityIntegrand(sol,adj_sol,QuadratureAdjoint(abstol=1e-14,reltol=1e-14))
res,err = quadgk(integrand,0.0,10.0,atol=1e-14,rtol=1e-12)

@test isapprox(res, easy_res, rtol = 1e-10)
@test isapprox(res, easy_res2, rtol = 1e-10)
@test isapprox(res, easy_res22, rtol = 1e-10)
@test isapprox(res, easy_res23, rtol = 1e-10)
@test isapprox(res, easy_res3, rtol = 1e-10)
@test isapprox(res, easy_res32, rtol = 1e-10)
@test isapprox(res, easy_res4, rtol = 1e-10)
@test isapprox(res, easy_res42, rtol = 1e-10)
@test isapprox(res, easy_res43, rtol = 1e-10)
@test isapprox(res, easy_res5, rtol = 1e-7)
@test isapprox(res, easy_res6, rtol = 1e-9)
@test isapprox(res, easy_res62, rtol = 1e-9)
@test all(easy_res6 .== easy_res7)  # should be the same!
@test isapprox(res, easy_res8, rtol = 1e-9)
@test_broken isapprox(res, easy_res9, rtol = 1e-9)
@test isapprox(res, easy_res10, rtol = 1e-9)
@test isapprox(res, easy_res11, rtol = 1e-9)

println("OOP adjoint sensitivities ")

_,easy_res = adjoint_sensitivities(soloop,Tsit5(),dg,t,abstol=1e-14,
                                 reltol=1e-14)
_,easy_res2 = adjoint_sensitivities(soloop,Tsit5(),dg,t,abstol=1e-14,
                                 reltol=1e-14,
                                 sensealg=QuadratureAdjoint(abstol=1e-14,reltol=1e-14))
@test_broken easy_res22 = adjoint_sensitivities(soloop,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,
                                  sensealg=QuadratureAdjoint(autojacvec=false,abstol=1e-14,reltol=1e-14))[1] isa AbstractArray
_,easy_res2 = adjoint_sensitivities(soloop,Tsit5(),dg,t,abstol=1e-14,
                                    reltol=1e-14,
                                    sensealg=QuadratureAdjoint(abstol=1e-14,reltol=1e-14,compile=true))
_,easy_res3 = adjoint_sensitivities(soloop,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,
                                  sensealg=InterpolatingAdjoint())
@test easy_res32 = adjoint_sensitivities(soloop,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,
                                  sensealg=InterpolatingAdjoint(autojacvec=false))[1] isa AbstractArray
_,easy_res4 = adjoint_sensitivities(soloop,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,
                                  sensealg=BacksolveAdjoint())
@test easy_res42 = adjoint_sensitivities(soloop,Tsit5(),dg,t,abstol=1e-14,
                                   reltol=1e-14,
                                   sensealg=BacksolveAdjoint(autojacvec=false))[1] isa AbstractArray
_,easy_res5 = adjoint_sensitivities(soloop,Kvaerno5(nlsolve=NLAnderson(), smooth_est=false),
                                 dg,t,abstol=1e-12,
                                 reltol=1e-10,
                                 sensealg=BacksolveAdjoint())
_,easy_res6 = adjoint_sensitivities(soloop_nodense,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,
                                  sensealg=InterpolatingAdjoint(checkpointing=true),
                                  checkpoints=soloop_nodense.t[1:5:end])
@test_broken easy_res62 = adjoint_sensitivities(soloop_nodense,Tsit5(),dg,t,abstol=1e-14,
                                   reltol=1e-14,
                                   sensealg=InterpolatingAdjoint(checkpointing=true,autojacvec=false),
                                   checkpoints=soloop_nodense.t[1:5:end])

_,easy_res8 = adjoint_sensitivities(soloop_nodense,Tsit5(),dg,t,abstol=1e-14,
                                   reltol=1e-14,
                                   sensealg=InterpolatingAdjoint(autojacvec=DiffEqSensitivity.TrackerVJP()))
_,easy_res9 = adjoint_sensitivities(soloop_nodense,Tsit5(),dg,t,abstol=1e-14,
                                   reltol=1e-14,
                                   sensealg=InterpolatingAdjoint(autojacvec=DiffEqSensitivity.ZygoteVJP()))
_,easy_res10 = adjoint_sensitivities(soloop_nodense,Tsit5(),dg,t,abstol=1e-14,
                                   reltol=1e-14,
                                   sensealg=InterpolatingAdjoint(autojacvec=DiffEqSensitivity.ReverseDiffVJP())
                                   )
_,easy_res11 = adjoint_sensitivities(soloop_nodense,Tsit5(),dg,t,abstol=1e-14,
                                     reltol=1e-14,
                                     sensealg=InterpolatingAdjoint(autojacvec=DiffEqSensitivity.ReverseDiffVJP(true))
                                     )

@test isapprox(res, easy_res, rtol = 1e-10)
@test isapprox(res, easy_res2, rtol = 1e-10)
@test isapprox(res, easy_res22, rtol = 1e-10)
@test isapprox(res, easy_res23, rtol = 1e-10)
@test isapprox(res, easy_res3, rtol = 1e-10)
@test isapprox(res, easy_res32, rtol = 1e-10)
@test isapprox(res, easy_res4, rtol = 1e-10)
@test isapprox(res, easy_res42, rtol = 1e-10)
@test isapprox(res, easy_res5, rtol = 1e-9)
@test isapprox(res, easy_res6, rtol = 1e-10)
@test isapprox(res, easy_res62, rtol = 1e-9)
@test isapprox(res, easy_res8, rtol = 1e-9)
@test isapprox(res, easy_res9, rtol = 1e-9)
@test isapprox(res, easy_res10, rtol = 1e-9)
@test isapprox(res, easy_res11, rtol = 1e-9)

println("Calculate adjoint sensitivities ")

_,easy_res8 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,
                                  save_everystep=false,save_start=false,
                                  sensealg=BacksolveAdjoint())
_,easy_res82 = adjoint_sensitivities(solb,Tsit5(),dg,t,abstol=1e-14,
                                   reltol=1e-14,
                                   save_everystep=false,save_start=false,
                                   sensealg=BacksolveAdjoint(checkpointing=false))

@test isapprox(res, easy_res8, rtol = 1e-9)
@test isapprox(res, easy_res82, rtol = 1e-9)

_,end_only_res = adjoint_sensitivities(sol_end,Tsit5(),dg,t,abstol=1e-14,
                                  reltol=1e-14,
                                  save_everystep=false,save_start=false,
                                  sensealg=BacksolveAdjoint())

@test isapprox(res, end_only_res, rtol = 1e-9)

println("Calculate adjoint sensitivities from autodiff & numerical diff")
function G(p)
  tmp_prob = remake(prob,u0=convert.(eltype(p),prob.u0),p=p)
  sol = solve(tmp_prob,Tsit5(),abstol=1e-14,reltol=1e-14,sensealg=DiffEqBase.SensitivityADPassThrough(),saveat=t)
  A = Array(sol)
  sum(((2 .- A).^2)./2)
end
G([1.5,1.0,3.0,1.0])
res2 = ForwardDiff.gradient(G,[1.5,1.0,3.0,1.0])
res3 = Calculus.gradient(G,[1.5,1.0,3.0,1.0])

import Tracker
res4 = Tracker.gradient(G,[1.5,1.0,3.0,1.0])[1]

import ReverseDiff
res5 = ReverseDiff.gradient(G,[1.5,1.0,3.0,1.0])

@test norm(res' .- res2) < 1e-7
@test norm(res' .- res3) < 1e-5
@test norm(res' .- res4) < 1e-6
@test norm(res' .- res5) < 1e-6

# check other t handling

t2 = [0.5, 1.0]
t3 = [0.0, 0.5, 1.0]
t4 = [0.5, 1.0, 10.0]

_,easy_res2 = adjoint_sensitivities(sol,Tsit5(),dg,t2,abstol=1e-14,
                                  reltol=1e-14)
_,easy_res3 = adjoint_sensitivities(sol,Tsit5(),dg,t3,abstol=1e-14,
                                  reltol=1e-14)
_,easy_res4 = adjoint_sensitivities(sol,Tsit5(),dg,t4,abstol=1e-14,
                                  reltol=1e-14)

function G(p,ts)
  tmp_prob = remake(prob,u0=convert.(eltype(p),prob.u0),p=p)
  sol = solve(tmp_prob,Tsit5(),abstol=1e-10,reltol=1e-10,sensealg=DiffEqBase.SensitivityADPassThrough(),saveat=ts)
  A = convert(Array,sol)
  sum(((2 .- A).^2)./2)
end
res2 = ForwardDiff.gradient(p->G(p,t2),[1.5,1.0,3.0,1.0])
res3 = ForwardDiff.gradient(p->G(p,t3),[1.5,1.0,3.0,1.0])
res4 = ForwardDiff.gradient(p->G(p,t4),[1.5,1.0,3.0,1.0])

@test easy_res2' ≈ res2
@test easy_res3' ≈ res3
@test easy_res4' ≈ res4

println("Adjoints of u0")

function dg(out,u,p,t,i)
  out .= 1 .- u
end

ū0,adj = adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                         reltol=1e-14)

_,adjnou0 = adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                        reltol=1e-14)

ū02,adj2 = adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                                    sensealg=BacksolveAdjoint(),
                                    reltol=1e-14)

ū022,adj22 = adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                                    sensealg=BacksolveAdjoint(autojacvec=false),
                                    reltol=1e-14)

ū023,adj23 = adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                                    sensealg=BacksolveAdjoint(autojacvec=false,checkpointing=false),
                                    reltol=1e-14)

ū03,adj3 = adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                                    sensealg=InterpolatingAdjoint(),
                                    reltol=1e-14)

ū032,adj32 = adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                                    sensealg=InterpolatingAdjoint(autojacvec=false),
                                    reltol=1e-14)

ū04,adj4 = adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                                    sensealg=InterpolatingAdjoint(checkpointing=true),
                                    checkpoints=sol.t[1:500:end],
                                    reltol=1e-14)

@test_nowarn adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                                      sensealg=InterpolatingAdjoint(checkpointing=true),
                                      checkpoints=sol.t[1:5:end],
                                      reltol=1e-14)

ū042,adj42 = adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                                    sensealg=InterpolatingAdjoint(checkpointing=true,autojacvec=false),
                                    checkpoints=sol.t[1:500:end],
                                    reltol=1e-14)

ū05,adj5 = adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                                    sensealg=QuadratureAdjoint(abstol=1e-14,reltol=1e-14),
                                    reltol=1e-14)

ū052,adj52 = adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                                    sensealg=QuadratureAdjoint(autojacvec=false,abstol=1e-14,reltol=1e-14),
                                    reltol=1e-14)

ū05,adj53 = adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                                    sensealg=QuadratureAdjoint(abstol=1e-14,reltol=1e-14,compile=true),
                                    reltol=1e-14)

ū0args,adjargs = adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                        save_everystep=false, save_start=false,
                        sensealg=BacksolveAdjoint(),
                        reltol=1e-14)

ū0args2,adjargs2 = adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                        save_everystep=false, save_start=false,
                        sensealg=InterpolatingAdjoint(),
                        reltol=1e-14)

res = ForwardDiff.gradient(prob.u0) do u0
  tmp_prob = remake(prob,u0=u0)
  sol = solve(tmp_prob,Tsit5(),abstol=1e-14,reltol=1e-14,saveat=t)
  A = convert(Array,sol)
  sum(((1 .- A).^2)./2)
end

@test ū0 ≈ res rtol = 1e-10
@test ū02 ≈ res rtol = 1e-10
@test ū022 ≈ res rtol = 1e-10
@test ū023 ≈ res rtol = 1e-10
@test ū03 ≈ res rtol = 1e-10
@test ū032 ≈ res rtol = 1e-10
@test ū04 ≈ res rtol = 1e-10
@test ū042 ≈ res rtol = 1e-10
@test ū05 ≈ res rtol = 1e-10
@test ū052 ≈ res rtol = 1e-10
@test adj ≈ adjnou0 rtol = 1e-10
@test adj ≈ adj2 rtol = 1e-10
@test adj ≈ adj22 rtol = 1e-10
@test adj ≈ adj23 rtol = 1e-10
@test adj ≈ adj3 rtol = 1e-10
@test adj ≈ adj32 rtol = 1e-10
@test adj ≈ adj4 rtol = 1e-10
@test adj ≈ adj42 rtol = 1e-10
@test adj ≈ adj5 rtol = 1e-10
@test adj ≈ adj52 rtol = 1e-10
@test adj ≈ adj53 rtol = 1e-10

@test ū0args ≈ res rtol = 1e-10
@test adjargs ≈ adj rtol = 1e-10
@test ū0args2 ≈ res rtol = 1e-10
@test adjargs2 ≈ adj rtol = 1e-10

println("Do a continuous adjoint problem")

# Energy calculation
g(u,p,t) = (sum(u).^2) ./ 2
# Gradient of (u1 + u2)^2 / 2
function dg(out,u,p,t)
  out[1]= u[1] + u[2]
  out[2]= u[1] + u[2]
end

adj_prob = ODEAdjointProblem(sol,QuadratureAdjoint(abstol=1e-14,reltol=1e-14),g,nothing,dg)
adj_sol = solve(adj_prob,Tsit5(),abstol=1e-14,reltol=1e-10)
integrand = AdjointSensitivityIntegrand(sol,adj_sol,QuadratureAdjoint(abstol=1e-14,reltol=1e-14))
res,err = quadgk(integrand,0.0,10.0,atol=1e-14,rtol=1e-10)

println("Test the `adjoint_sensitivities` utility function")
_,easy_res = adjoint_sensitivities(sol,Tsit5(),g,nothing,dg,abstol=1e-14,
                                 reltol=1e-14)
println("2")
_,easy_res2 = adjoint_sensitivities(sol,Tsit5(),g,nothing,dg,abstol=1e-14,
                                  reltol=1e-14,
                                  sensealg=InterpolatingAdjoint())
_,easy_res22 = adjoint_sensitivities(sol,Tsit5(),g,nothing,dg,abstol=1e-14,
                                   reltol=1e-14,
                                   sensealg=InterpolatingAdjoint(autojacvec=false))
println("23")
_,easy_res23 = adjoint_sensitivities(sol,Tsit5(),g,nothing,dg,abstol=1e-14,
                                   reltol=1e-14,
                                   sensealg=QuadratureAdjoint(abstol=1e-14,reltol=1e-14))
_,easy_res232 = adjoint_sensitivities(sol,Tsit5(),g,nothing,dg,abstol=1e-14,
                                  reltol=1e-14,
                                  sensealg=QuadratureAdjoint(abstol=1e-14,reltol=1e-14,compile=false))
_,easy_res24 = adjoint_sensitivities(sol,Tsit5(),g,nothing,dg,abstol=1e-14,
                                   reltol=1e-14,
                                   sensealg=QuadratureAdjoint(autojacvec=false,abstol=1e-14,reltol=1e-14))
println("25")
_,easy_res25 = adjoint_sensitivities(sol,Tsit5(),g,nothing,dg,abstol=1e-14,
                                   reltol=1e-14,
                                   sensealg=BacksolveAdjoint())
_,easy_res26 = adjoint_sensitivities(sol,Tsit5(),g,nothing,dg,abstol=1e-14,
                                   reltol=1e-14,
                                   sensealg=BacksolveAdjoint(autojacvec=false))
_,easy_res262 = adjoint_sensitivities(sol,Tsit5(),g,nothing,dg,abstol=1e-14,
                                   reltol=1e-14,
                                   sensealg=BacksolveAdjoint(autojacvec=false,checkpointing=false))
println("27")
_,easy_res27 = adjoint_sensitivities(sol,Tsit5(),g,nothing,dg,abstol=1e-14,
                                  reltol=1e-14,
                                  checkpoints=sol.t[1:500:end],
                                  sensealg=InterpolatingAdjoint(checkpointing=true))
_,easy_res28 = adjoint_sensitivities(sol,Tsit5(),g,nothing,dg,abstol=1e-14,
                                  reltol=1e-14,
                                  checkpoints=sol.t[1:500:end],
                                  sensealg=InterpolatingAdjoint(checkpointing=true,autojacvec=false))
println("3")
_,easy_res3 = adjoint_sensitivities(sol,Tsit5(),g,nothing,abstol=1e-14,
                                  reltol=1e-14,
                                  sensealg=InterpolatingAdjoint())
_,easy_res32 = adjoint_sensitivities(sol,Tsit5(),g,nothing,abstol=1e-14,
                                   reltol=1e-14,
                                   sensealg=InterpolatingAdjoint(autojacvec=false))
println("33")
_,easy_res33 = adjoint_sensitivities(sol,Tsit5(),g,nothing,abstol=1e-14,
                                   reltol=1e-14,
                                   sensealg=QuadratureAdjoint(abstol=1e-14,reltol=1e-14))
_,easy_res34 = adjoint_sensitivities(sol,Tsit5(),g,nothing,abstol=1e-14,
                                   reltol=1e-14,
                                   sensealg=QuadratureAdjoint(autojacvec=false,abstol=1e-14,reltol=1e-14))
println("35")
_,easy_res35 = adjoint_sensitivities(sol,Tsit5(),g,nothing,abstol=1e-14,
                                   reltol=1e-14,
                                   sensealg=BacksolveAdjoint())
_,easy_res36 = adjoint_sensitivities(sol,Tsit5(),g,nothing,abstol=1e-14,
                                   reltol=1e-14,
                                   sensealg=BacksolveAdjoint(autojacvec=false))
println("37")
_,easy_res37 = adjoint_sensitivities(sol,Tsit5(),g,nothing,abstol=1e-14,
                                  reltol=1e-14,
                                  checkpoints=sol.t[1:500:end],
                                  sensealg=InterpolatingAdjoint(checkpointing=true))
_,easy_res38 = adjoint_sensitivities(sol,Tsit5(),g,nothing,abstol=1e-14,
                                  reltol=1e-14,
                                  checkpoints=sol.t[1:500:end],
                                  sensealg=InterpolatingAdjoint(checkpointing=true,autojacvec=false))

@test norm(easy_res .- res) < 1e-8
@test norm(easy_res2 .- res) < 1e-8
@test norm(easy_res22 .- res) < 1e-8
@test norm(easy_res23 .- res) < 1e-8
@test norm(easy_res232 .- res) < 1e-8
@test norm(easy_res24 .- res) < 1e-8
@test norm(easy_res25 .- res) < 1e-8
@test norm(easy_res26 .- res) < 1e-8
@test norm(easy_res262 .- res) < 1e-8
@test norm(easy_res27 .- res) < 1e-8
@test norm(easy_res28 .- res) < 1e-8
@test norm(easy_res3 .- res) < 1e-8
@test norm(easy_res32 .- res) < 1e-8
@test norm(easy_res33 .- res) < 1e-8
@test norm(easy_res34 .- res) < 1e-8
@test norm(easy_res35 .- res) < 1e-8
@test norm(easy_res36 .- res) < 1e-8
@test norm(easy_res37 .- res) < 1e-8
@test norm(easy_res38 .- res) < 1e-8

println("Calculate adjoint sensitivities from autodiff & numerical diff")
function G(p)
  tmp_prob = remake(prob,u0=eltype(p).(prob.u0),p=p,
                    tspan=eltype(p).(prob.tspan))
  sol = solve(tmp_prob,Tsit5(),abstol=1e-14,reltol=1e-14)
  res,err = quadgk((t)-> (sum(sol(t)).^2)./2,0.0,10.0,atol=1e-14,rtol=1e-10)
  res
end
res2 = ForwardDiff.gradient(G,[1.5,1.0,3.0,1.0])
res3 = Calculus.gradient(G,[1.5,1.0,3.0,1.0])

@test norm(res' .- res2) < 1e-8
@test norm(res' .- res3) < 1e-6

# Buffer length test
f = (du, u, p, t) -> du .= 0
p = zeros(3); u = zeros(50)
prob = ODEProblem(f,u,(0.0,10.0),p)
sol = solve(prob,Tsit5(),abstol=1e-14,reltol=1e-14)
@test_nowarn _,res = adjoint_sensitivities(sol,Tsit5(),dg,t,abstol=1e-14,
                                 reltol=1e-14)

@testset "Checkpointed backsolve" begin
  using DiffEqSensitivity, OrdinaryDiffEq
  tf = 10.0
  function lorenz(du,u,p,t)
    σ, ρ, β = p
    du[1] = σ*(u[2]-u[1])
    du[2] = u[1]*(ρ-u[3]) - u[2]
    du[3] = u[1]*u[2] - β*u[3]
    return nothing
  end
  prob_lorenz = ODEProblem(lorenz, [1.0, 0.0, 0.0], (0, tf), [10, 28, 8/3])
  sol_lorenz = solve(prob_lorenz,Tsit5(),reltol=1e-6,abstol=1e-9)
  function dg(out,u,p,t,i)
    (out.=2.0.-u)
  end
  t = 0:0.1:tf
  _,easy_res1 = adjoint_sensitivities(sol_lorenz,Tsit5(),dg,t,abstol=1e-6,
                                    reltol=1e-9,
                                    sensealg=BacksolveAdjoint())
  _,easy_res2 = adjoint_sensitivities(sol_lorenz,Tsit5(),dg,t,abstol=1e-6,
                                    reltol=1e-9,
                                    sensealg=InterpolatingAdjoint())
  _,easy_res3 = adjoint_sensitivities(sol_lorenz,Tsit5(),dg,t,abstol=1e-6,
                                    reltol=1e-9,
                                    sensealg=BacksolveAdjoint(),
                                    checkpoints=sol_lorenz.t[1:10:end])
  _,easy_res4 = adjoint_sensitivities(sol_lorenz,Tsit5(),dg,t,abstol=1e-6,
                                    reltol=1e-9,
                                    sensealg=BacksolveAdjoint(),
                                    checkpoints=sol_lorenz.t[1:20:end])
  # cannot finish in a reasonable amount of time
  @test_skip adjoint_sensitivities(sol_lorenz,Tsit5(),dg,t,abstol=1e-6,
                                   reltol=1e-9,
                                   sensealg=BacksolveAdjoint(checkpointing=false))
  @test easy_res2 ≈ easy_res1 rtol=1e-5
  @test easy_res2 ≈ easy_res3 rtol=1e-5
  @test easy_res2 ≈ easy_res4 rtol=1e-4

  ū1,adj1 = adjoint_sensitivities(sol_lorenz,Tsit5(),dg,t,abstol=1e-6,
                                     reltol=1e-9,
                                     sensealg=BacksolveAdjoint())
  ū2,adj2 = adjoint_sensitivities(sol_lorenz,Tsit5(),dg,t,abstol=1e-6,
                                     reltol=1e-9,
                                     sensealg=InterpolatingAdjoint())
  ū3,adj3 = adjoint_sensitivities(sol_lorenz,Tsit5(),dg,t,abstol=1e-6,
                                     reltol=1e-9,
                                     sensealg=BacksolveAdjoint(),
                                     checkpoints=sol_lorenz.t[1:10:end])
  ū4,adj4 = adjoint_sensitivities(sol_lorenz,Tsit5(),dg,t,abstol=1e-6,
                                     reltol=1e-9,
                                     sensealg=BacksolveAdjoint(),
                                     checkpoints=sol_lorenz.t[1:20:end])
  # cannot finish in a reasonable amount of time
  @test_skip adjoint_sensitivities(sol_lorenz,Tsit5(),dg,t,abstol=1e-6,
                                      reltol=1e-9,
                                      sensealg=BacksolveAdjoint(checkpointing=false))
  @test ū2 ≈ ū1 rtol=1e-5
  @test adj2 ≈ adj1 rtol=1e-5
  @test ū2 ≈ ū3 rtol=1e-5
  @test adj2 ≈ adj3 rtol=1e-5
  @test ū2 ≈ ū4 rtol=1e-4
  @test adj2 ≈ adj4 rtol=1e-4

  # LQR Tests from issue https://github.com/SciML/DiffEqSensitivity.jl/issues/300
  x_dim = 2
  T = 40.0

  cost = (x, u) -> x'*x
  params =  [-0.4142135623730951, 0.0, -0.0, -0.4142135623730951, 0.0, 0.0]

  function dynamics!(du,u,p,t)
      du[1] = -u[1] + tanh(p[1]*u[1]+p[2]*u[2])
      du[2] = -u[2] + tanh(p[3]*u[1]+p[4]*u[2])
  end

  function backsolve_grad(sol, lqr_params, checkpointing)
      bwd_sol = solve(
          ODEAdjointProblem(
              sol,
              BacksolveAdjoint(checkpointing = checkpointing),
              (x, lqr_params, t) -> cost(x,lqr_params),
          ),
          Tsit5(),
          dense = false,
          save_everystep = false,
      )

      bwd_sol.u[end][1:end-x_dim]
      #fwd_sol, bwd_sol
  end

  x0 = ones(x_dim)
  fwd_sol = solve(
      ODEProblem(dynamics!, x0, (0, T), params),
      Tsit5(),abstol=1e-9, reltol=1e-9,
      u0 = x0,
      p = params,
      dense = false,
      save_everystep = true
  )

  backsolve_results = backsolve_grad(fwd_sol, params, false)
  backsolve_checkpointing_results = backsolve_grad(fwd_sol, params, true)

  @test backsolve_results != backsolve_checkpointing_results

  int_u0, int_p = adjoint_sensitivities(fwd_sol,Tsit5(),(x, params, t)->cost(x,params), nothing, sensealg=InterpolatingAdjoint())

  @test isapprox(-backsolve_checkpointing_results[1:length(x0)], int_u0, rtol=1e-10)
  @test isapprox(backsolve_checkpointing_results[(1:length(params)) .+ length(x0)], int_p', rtol=1e-10)
end

using Test
using LinearAlgebra, DiffEqSensitivity, OrdinaryDiffEq, ForwardDiff, QuadGK
@testset "Adjoint of differential algebric equations with mass matrix" begin
  function G(p, prob, ts, cost)
    tmp_prob_mm = remake(prob,u0=convert.(eltype(p),prob.u0),p=p)
    sol = solve(tmp_prob_mm,Rodas5(autodiff=false),abstol=1e-14,reltol=1e-14,saveat=ts)
    cost(sol)
  end
  alg = Rodas5(autodiff=false)
  @testset "Fully ranked mass matrix" begin
    @info "discrete cost"
    A = [1 2 3; 4 5 6; 7 8 9]
    function foo(du, u, p, t)
      mul!(du, A, u)
      du .= du .+ p
      du[2] += sum(p)
      return nothing
    end
    mm = -[1 2 4; 2 3 7; 1 3 41]
    u0 = [1, 2.0, 3]
    p = [1.0, 2.0, 3]
    prob_mm = ODEProblem(ODEFunction(foo, mass_matrix=mm), u0, (0, 1.0), p)
    sol_mm = solve(prob_mm, Rodas5(), reltol=1e-14, abstol=1e-14)

    ts = 0:0.01:1
    dg(out,u,p,t,i) = out .= -1
    _, res = adjoint_sensitivities(sol_mm,alg,dg,ts,abstol=1e-5,reltol=1e-5,sensealg=QuadratureAdjoint())
    reference_sol = ForwardDiff.gradient(p->G(p, prob_mm, ts, sum),vec(p))
    @test res' ≈ reference_sol rtol=1e-3

    _, res_interp = adjoint_sensitivities(sol_mm,alg,dg,ts,abstol=1e-5,reltol=1e-5,sensealg=InterpolatingAdjoint())
    @test res_interp ≈ res rtol = 1e-2
    _, res_interp2 = adjoint_sensitivities(sol_mm,alg,dg,ts,abstol=1e-5,reltol=1e-5,sensealg=InterpolatingAdjoint(checkpointing=true),checkpoints=sol_mm.t[1:10:end])
    @test res_interp2 ≈ res rtol = 1e-2

    # backsolve doesn't work
    _, res_bs = adjoint_sensitivities(sol_mm,alg,dg,ts,abstol=1e-5,reltol=1e-5,sensealg=BacksolveAdjoint(checkpointing=false))
    @test res_bs ≈ res rtol = 1e-2
    _, res_bs2 = adjoint_sensitivities(sol_mm,alg,dg,ts,abstol=1e-5,reltol=1e-5,sensealg=BacksolveAdjoint(checkpointing=true),checkpoints=sol_mm.t)
    @test res_bs2 ≈ res rtol = 1e-2

    @info "continuous cost"
    g_cont(u,p,t) = (sum(u).^2) ./ 2
    dg_cont(out,u,p,t) = out .= sum(u)
    _,easy_res_cont = adjoint_sensitivities(sol_mm,alg,g_cont,nothing,
                                     dg_cont,abstol=1e-10,reltol=1e-10,
                                     sensealg=QuadratureAdjoint())
    function G_cont(p)
      tmp_prob_mm = remake(prob_mm,u0=eltype(p).(prob_mm.u0),p=p,
                        tspan=eltype(p).(prob_mm.tspan))
      sol = solve(tmp_prob_mm,Rodas5(autodiff=false),abstol=1e-14,reltol=1e-14)
      res,err = quadgk((t)-> (sum(sol(t)).^2)./2,prob_mm.tspan...,atol=1e-14,rtol=1e-10)
      res
    end
    reference_sol_cont = ForwardDiff.gradient(G_cont, p)
    @test easy_res_cont' ≈ reference_sol_cont rtol=1e-3
  end

  @testset "Singular mass matrix" begin
    function rober(du,u,p,t)
      y₁,y₂,y₃ = u
      k₁,k₂,k₃ = p
      du[1] = -k₁*y₁+k₃*y₂*y₃
      du[2] =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
      du[3] =  y₁ + y₂ + y₃ - 1
      nothing
    end
    function rober(u,p,t)
      y₁,y₂,y₃ = u
      k₁,k₂,k₃ = p
      return [-k₁*y₁+k₃*y₂*y₃,
       k₁*y₁-k₂*y₂^2-k₃*y₂*y₃,
       y₁ + y₂ + y₃ - 1]
    end
    M = [1. 0  0
         0  1. 0
         0  0  0]
    for iip in [true, false]
      f = ODEFunction{iip}(rober,mass_matrix=M)
      p = [0.04,3e7,1e4]

      prob_singular_mm = ODEProblem(f,[1.0,0.0,0.0],(0.0,100),p)
      sol_singular_mm = solve(prob_singular_mm,Rodas5(),reltol=1e-8,abstol=1e-8)
      ts = [50, sol_singular_mm.t[end]]
      dg_singular(out,u,p,t,i) = (fill!(out, 0); out[end] = -1)
      _, res = adjoint_sensitivities(sol_singular_mm,alg,dg_singular,ts,abstol=1e-5,reltol=1e-5,sensealg=QuadratureAdjoint())
      reference_sol = ForwardDiff.gradient(p->G(p, prob_singular_mm, ts, sol->sum(last, sol.u)), vec(p))
      @test res' ≈ reference_sol rtol = 1e-2

      _, res_interp = adjoint_sensitivities(sol_singular_mm,alg,dg_singular,ts,abstol=1e-5,reltol=1e-5,sensealg=InterpolatingAdjoint())
      @test res_interp ≈ res rtol = 1e-2
      _, res_interp2 = adjoint_sensitivities(sol_singular_mm,alg,dg_singular,ts,abstol=1e-5,reltol=1e-5,sensealg=InterpolatingAdjoint(checkpointing=true),checkpoints=sol_singular_mm.t[1:10:end])
      @test res_interp2 ≈ res rtol = 1e-2

      # backsolve doesn't work
      _, res_bs = adjoint_sensitivities(sol_singular_mm,alg,dg_singular,ts,abstol=1e-5,reltol=1e-5,sensealg=BacksolveAdjoint(checkpointing=false))
      @test_broken res_bs ≈ res rtol = 1e-2
      _, res_bs2 = adjoint_sensitivities(sol_singular_mm,alg,dg_singular,ts,abstol=1e-5,reltol=1e-5,sensealg=BacksolveAdjoint(checkpointing=true),checkpoints=sol_singular_mm.t)
      @test_broken res_bs2 ≈ res rtol = 1e-2
    end
  end
end


text:
module TestIndexing

using Base.Test
using RLEVectors

@testset begin

# ind2run
@test ind2run(RLEVector([3,3,4,4,5,5,6,6,7,7]),5) == 3
@test ind2run(RLEVector([3,3,4,4,5,5,6,6,7,7]),6) == 3
@test ind2run(RLEVector([3,3,4,4,5,5,6,6,7,7]),7) == 4
@test ind2run(RLEVector([3,3,4,4,5,5,6,6,7,7]),7:9) == 4:5
@test ind2run(RLEVector([3,3,4,4,5,5,6,6,7,7]),[3, 5, 6]) == [2, 3, 3]

# ind2runcontext
@test ind2runcontext(RLEVector([3,3,4,4,5,5,6,6,7,7]),5) == (3,1,1)
@test ind2runcontext(RLEVector([3,3,4,4,5,5,6,6,7,7]),6) == (3,2,0)
@test ind2runcontext(RLEVector([3,3,4,4,5,5,6,6,7,7]),7) == (4,1,1)
@test ind2runcontext(RLEVector([3,3,4,4,5,5,6,6,7,7]),7:9) == (4,5,1,1)

## getindex for single position
x = RLEVectors.RLEVector([1,2,3,4],[2,4,6,8])
@test x[1] == 1
@test x[3] == 2
@test x[8] == 4
@test_throws BoundsError x[9]

## getindex on multiple positions
x = RLEVectors.RLEVector([1,2,3,4],[2,4,6,8])
@test x[ [1,2,3] ] == RLEVector([1,1,2])
@test x[ [8,2,4] ] == RLEVector([4,1,2])
@test x[ 2:5 ] == RLEVector([1,2,2,3])
@test x[ 5:-1:2 ] == RLEVector([3,2,2,1])
@test x[ 1:end ] == x

## getindex with logical
#x = RLEVectors.RLEVector([1,2],[2,4])
@test x[ [true, false, false, true] ] == [1, 2]

## setindex! for single position
y = RLEVectors.RLEVector([3,4,4,5,5])
y[1] = 20
@test collect(y) == [20,4,4,5,5]

# No change
x = RLEVectors.RLEVector([1,2,3,4],[2,4,6,8])
x[4] = 2
@test collect(x) == [1,1,2,2,3,3,4,4]
@test x.runvalues == [1,2,3,4]
@test x.runends == [2,4,6,8]

# left of run
x = RLEVectors.RLEVector([1,2,3,4],[2,4,6,8])
x[1] = 9
@test collect(x) == [9,1,2,2,3,3,4,4]
@test x.runvalues == [9,1,2,3,4]
@test x.runends == [1,2,4,6,8]

# right of run
x = RLEVectors.RLEVector([1,2,3,4],[2,4,6,8])
x[2] = 7
@test collect(x) == [1,7,2,2,3,3,4,4]
@test x.runvalues == [1,7,2,3,4]
@test x.runends == [1,2,4,6,8]

# left of run making repeat
x = RLEVectors.RLEVector([1,2,3,4],[2,4,6,8])
x[3] = 1
@test collect(x) == [1,1,1,2,3,3,4,4]
@test x.runvalues == [1,2,3,4]
@test x.runends == [3,4,6,8]

# right of run making repeat
x = RLEVectors.RLEVector([1,2,3,4],[2,4,6,8])
x[4] = 3
@test collect(x) == [1,1,2,3,3,3,4,4]
@test x.runvalues == [1,2,3,4]
@test x.runends == [2,3,6,8]

# middle of run
x = RLEVectors.RLEVector([1,2,3],[3,6,9])
x[5] = 5
@test collect(x) == [1,1,1,2,5,2,3,3,3]
@test x.runvalues == [1,2,5,2,3]
@test x.runends == [3,4,5,6,9]

# ends
x = RLEVector([1,2,3],[3,6,9])
x[1] = 10
@test collect(x) == [10,1,1,2,2,2,3,3,3]
x = RLEVector([1,2,3],[3,6,9])
x[9] = 10
@test collect(x) == [1,1,1,2,2,2,3,3,10]

# singleton ranges
x = RLEVector([1,2,2,3,3])
x[1] = 5
@test collect(x) == [5,2,2,3,3]
x = RLEVector([1,2,2,3,3])
x[1] = 2
@test collect(x) == [2,2,2,3,3]
x = RLEVector([1,1,2,3,3])
x[3] = 3
@test collect(x) == [1,1,3,3,3]
x = RLEVector([1,1,2,3,3])
x[3] = 1
@test collect(x) == [1,1,1,3,3]
x = RLEVector([3,1,2,1,3])
x[3] = 1
@test collect(x) == [3,1,1,1,3]

## range with scalar
# in middle, no match
x = RLEVectors.RLEVector([1,2,3,4],[4,8,12,16])
x[10:11] = 5
@test x.runvalues == [1,2,3,5,3,4]
@test x.runends == [4,8,9,11,12,16]

# just left, no match
x = RLEVectors.RLEVector([1,2,3,4],[4,8,12,16])
x[9:10] = 5
@test x.runvalues == [1,2,5,3,4]
@test x.runends == [4,8,10,12,16]

# just right, no match
x = RLEVectors.RLEVector([1,2,3,4],[4,8,12,16])
x[11:12] = 5
@test x.runvalues == [1,2,3,5,4]
@test x.runends == [4,8,10,12,16]

# just left, match
x = RLEVectors.RLEVector([1,2,3,4],[4,8,12,16])
x[9:10] = 2
@test x.runvalues == [1,2,3,4]
@test x.runends == [4,10,12,16]

# just right, match
x = RLEVectors.RLEVector([1,2,3,4],[4,8,12,16])
x[11:12] = 4
@test x.runvalues == [1,2,3,4]
@test x.runends == [4,8,10,16]

# span, no match
x = RLEVectors.RLEVector([1,2,3,4],[4,8,12,16])
x[9:12] = 5
@test x.runvalues == [1,2,5,4]
@test x.runends == [4,8,12,16]

# span, left match
x = RLEVectors.RLEVector([1,2,3,4],[4,8,12,16])
x[9:12] = 2
@test x.runvalues == [1,2,4]
@test x.runends == [4,12,16]

# span, right match
x = RLEVectors.RLEVector([1,2,3,4],[4,8,12,16])
x[9:12] = 4
@test x.runvalues == [1,2,4]
@test x.runends == [4,8,16]

# span, both match
x = RLEVectors.RLEVector([1,2,3,2],[4,8,12,16])
x[9:12] = 2
@test x.runvalues == [1,2]
@test x.runends == [4,16]

# reverse range with scalar
x = RLEVectors.RLEVector([1,2,3,4],[2,4,6,8])
x[4:-1:2] = 5
@test collect(x) == [1,5,5,5,3,3,4,4]
@test x.runvalues == [1,5,3,4]
@test x.runends == [1,4,6,8]

## range with vector
x = RLEVectors.RLEVector([1,2,3,4],[2,4,6,8])
x[2:4] = [5,6,7]
@test x[2:4] == RLEVector([5,6,7])

## range with range
x = RLEVectors.RLEVector([1,2,3,4],[2,4,6,8])
x[2:4] = 5:7
@test x[ [2,3,4] ] == [5,6,7]

# reverse range with vector
x = RLEVectors.RLEVector([1,2,3,4],[2,4,6,8])
x[4:-1:2] = [5,6,7]
@test x[4:-1:2] == [5,6,7]

# Colon
x = RLEVectors.RLEVector([1,2,3,4],[2,4,6,8])
@test x[:] == x[1:end]
x[:] = 4
@test x == RLEVector([4 for i in 1:8])

# Logical
x = RLEVectors.RLEVector([1,2],[2,4])
@test x[ [ true,true,true,false ] ] == x[ [1,2,3] ]
x[ [true,true,true,false] ] = 4
@test x == RLEVector([4,4,4,2])
x = RLEVectors.RLEVector([1,2],[2,4])
x[ [true,true,true,false] ] = [4,5,6]
@test x == RLEVector([4,5,6,2])

# head and tail
x = RLEVectors.RLEVector([1,2,3,4],[2,4,6,8])
@test head(x) == [1,1,2,2,3,3]
@test head(x,3) == [1,1,2]
@test tail(x) == [2,2,3,3,4,4]
@test tail(x,3) == [3,4,4]

end # testset

end # module


text:
_NI() = error("Not implemented")

# Interface for ordering
abstract type Ordering
end

get_var(ordering::Ordering, layer::Int) = _NI

# Simple ordering implementations

# No ordering: variable corresponds to layer
struct NoOrdering <: Ordering
end

function get_var(ordering::NoOrdering, layer::Int)::Int
    return layer
end

# Fixed ordering
struct FixedOrdering <: Ordering
    fixed_ordering::Array{Int, 1}
end

function FixedOrdering(fixed_ordering::Array{Int, 1})::FixedOrdering
    # Check if the ordering is a valid permutation
    if !isperm(fixed_ordering)
        error("Invalid fixed ordering")
    end
    return new(fixed_ordering)
end

function get_var(ordering::FixedOrdering, layer::Int)::Array{Int, 1}
    return ordering.fixed_ordering[layer]
end

max_stars_repo_name: ctjandra/ddopt-ecp

max_issues_repo_name: ctjandra/ddopt-ecp

max_forks_repo_name: ctjandra/ddopt-ecp


text:
"""
    sampling_frequency(x) -> Hz

Sampling frequency (in Hz) of all the data in the recording, regardless of
their type (e.g., 2400)
"""

"""
    power_line_frequenc(x) -> Hz

Frequency (in Hz) of the power grid at the geographical location of the MEG
instrument (i.e. 50 or 60).
"""