text:
export WaterWaves

"""
    WaterWaves(param; kwargs)

Define an object of type `AbstractModel` in view of solving the initial-value problem for
the water waves system (via conformal mapping, see [Zakharov, Dyachenko and Vasilyev](https://doi.org/10.1016/S0997-7546(02)01189-5)).

# Argument
`param` is of type `NamedTuple` and must contain
- dimensionless parameters `ϵ` (nonlinearity) and `μ` (dispersion);
- optionally, `ν` the shallow/deep water scaling factor. By default, `ν=1` if `μ≦1` and `ν=1/√μ` otherwise. Set the infinite-layer case if `ν=0`, or `μ=Inf`.
- numerical parameters to construct the mesh of collocation points, if `mesh` is not provided as a keyword argument.

## Optional keyword arguments
- `IL`: Set the infinite-layer case if `IL=true`, in which case `ϵ` is the steepness parameter. Default is `false`.
- `mesh`: the mesh of collocation points. By default, `mesh = Mesh(param)`;
- `method ∈ {1,2,3}`: method used to initialize the conformal mapping, as a fix-point problem `F(u)=u`
    - if `method == 1`, use standard contraction fix-point iteration;
    - if `method == 2`, use Newton algorithm with GMRES iterative solver to invert the Jacobian;
    - if `method == 3`, use Newton algorithm with direct solver to invert the Jacobian;
- `tol`: (relative) tolerance of the fix-point algorithm (default is `1e-16`);
- `maxiter`: the maximal number of iteration in the fix-point algorithm (default is `100`);
- `ktol`: tolerance of the low-pass Krasny filter (default is `0`, i.e. no filtering);
- `dealias`: dealiasing with Orlicz rule `1-dealias/(dealias+2)` (default is `0`, i.e. no dealiasing);
- `label`: a label for future references (default is `"water waves"`);
- `verbose`: prints information if `true` (default is `true`).

# Return values
Generate necessary ingredients for solving an initial-value problem via `solve!`:
1. a function `WaterWaves.f!` to be called in the explicit time-integration solver (also `WaterWaves.f1!` and `WaterWaves.f2!` for the symplectic Euler solver);
2. a function `WaterWaves.mapto` which from `(η,v)` of type `InitialData` provides the raw data matrix on which computations are to be executed;
3. a function `WaterWaves.mapfro` which from such data matrix returns the Tuple of real vectors `(η,v,x)`, where
    - `x` is a vector of collocation points (non-regularly spaced);
	- `η` is the values of surface deformation at collocation points `x`;
    - `v` is the derivative of the trace of the velocity potential at `x`.

"""
mutable struct WaterWaves <: AbstractModel

	label   :: String
	f!		:: Function
	f1!		:: Function
	f2!		:: Function
	mapto	:: Function
	mapfro	:: Function
	info	:: String

    function WaterWaves(param::NamedTuple;
						mesh = Mesh(param),
						IL	    = false,
						method  = 1,
						tol	    = 1e-15,
						maxiter = 100,
						dealias	= 0,
						ktol	= 0,
						label	= "water waves",
						verbose	= true)

		# Set up
		μ 	= param.μ
		ϵ 	= param.ϵ
		if !in(:ν,keys(param))
			if μ > 1
				ν = 1/sqrt(μ)
				nu = "1/√μ (deep water case)"
			else
				ν = 1
				nu = "1 (shallow water case)"
			end
		else
			ν = param.ν
			nu = "$ν"
		end
		if μ == Inf || ν==0 || IL == true # infinite layer case
			IL = true;  # IL (=Infinite layer) is a flag to be used in functions tanh,cotanh,xdcotanh
			μ = 1; ν = 1; # Then we should set μ=ν=1 in subsequent formula.
		end

		# Print information
		info = "Water waves system.\n"
		if IL == true
			info *= "├─Steepness parameter ϵ=$ϵ (infinite depth case).\n"
		else
			info *= "├─Shallowness parameter μ=$μ, nonlinearity parameter ϵ=$ϵ, \
					scaling parameter ν=$nu.\n"
		end
		info *= "├─Initial data built with method $method: "
		if method == 1 info *= "standard contraction fix-point iteration." end
		if method == 2 info *= "Newton algorithm with GMRES iterative solver to invert the Jacobian." end
		if method == 3 info *= "Newton algorithm with direct solver to invert the Jacobian." end
		info *= " Relative tolerance $tol and maximum $maxiter iterations.\n"
		if dealias == 0
			info *= "└─No dealiasing. "
		else
			info *= "└─Dealiasing with Orszag's rule adapted to power $(dealias + 1) nonlinearity. "
		end
		if ktol == 0
			info *= "No Krasny filter. "
		else
			info *= "Krasny filter with tolerance $ktol."
		end
		info *= "\nDiscretized with $(mesh.N) collocation points on [$(mesh.xmin), $(mesh.xmax)]."

		# Pre-allocate useful data

		x 	= mesh.x
		x₀ = x[1]
		k 	= mesh.k

    	∂ₓ	=  1im * k             # Differentiation Fourier multiplier
		if dealias == 0
			Π⅔ 	= ones(size(k)) # no dealiasing (Π⅔=Id)
		else
			K = (mesh.kmax-mesh.kmin)/(2+dealias)
			Π⅔ 	= abs.(k) .<= K # Dealiasing low-pass filter
		end

		fz = zeros(Complex{Float64}, mesh.N)
        z = zeros(Float64, mesh.N)
        phi = zeros(Float64, mesh.N)
		ξ = zeros(Float64, mesh.N)
		xv = zeros(Float64, mesh.N)
		Dxv = zeros(Float64, mesh.N)
		Dz = zeros(Float64, mesh.N)
		Dphi = zeros(Float64, mesh.N)
		J = zeros(Float64, mesh.N)
		M1 = zeros(Float64, mesh.N)
		M2 = zeros(Float64, mesh.N)
		q0 = 0.

		# some useful functions
		if IL == true
			cotanh   = x -> sign.(x)
			xdcotanh = x -> zero(x)
			mytanh   = x -> sign.(x)
		else
			cotanh   = x -> (x.!=0) ./ tanh.(x+(x.==0))
			xdcotanh = x -> -x ./ (sinh.(x+(x.==0)).^2)
			mytanh 	 = x -> tanh.(x)
		end

		function meanf(f) #computes the mean of a function, being given its FFT
			real(f[1])/length(f)
		end
		function mean(f) #computes the mean of a function
			sum(f)/length(f)
		end

		# Build raw data from physical data.
		# Use a conformal change of coordinate
		function mapto(data::InitialData)
		  if data.η(x)==zero(x) && data.v(x)==zero(x)  #useful when defining the type of initial data in some solvers
			  return zeros(Complex{Float64}, (mesh.N,2))
		  else
			# preallocate some variables to save memory use
			fη = zeros(Complex{Float64}, mesh.N)
			u0 = zeros(Complex{Float64}, mesh.N)
			nu0 = zeros(Complex{Float64}, mesh.N)
			z0 = zeros(Float64, mesh.N)
			v0 = zeros(Float64, mesh.N)
			dx = zeros(Float64, mesh.N)

			# We will solve F(u)=0
			function F( u )
				dx.=real.(ifft(u))
				fη .= fft(data.η(x+ϵ/ν*dx))
			    return u+ν*sqrt(μ)*1im*Π⅔ .*cotanh(sqrt(μ)*(1+ϵ*meanf(fη))*k).*fη
			end

			# The Newton alogorithm demands the Jacobian of F,
			# which involves the derivative of η
			hatdη=∂ₓ.*fft(data.η(x))  # Fourier coefficients of η'
	        function dη( x :: Vector{Float64} )
	            return real.(exp.(1im*(x.-x₀)*k')*hatdη/length(k))
	        end

			# Define the Jacobian as a linear map
			# to be used with GMRES iterative solver
			function JacF( u )
				dx.= real.(ifft(u))
				fη .= fft(data.η(x+ϵ/ν*dx))
				δ  = meanf(fη)
				dF(φ) = φ+ϵ*sqrt(μ)*1im*Π⅔ .* ( cotanh(sqrt(μ)*(1+ϵ*δ)*k).*fft(dη(x+ϵ/ν*dx) .* ifft(φ))
								+ϵ/(1+ϵ*δ)* (xdcotanh(sqrt(μ)*(1+ϵ*δ)*k).*fη ) * mean(dη(x+ϵ/ν*dx) .* ifft(φ)) )
				return LinearMap(dF, length(u); issymmetric=false, ismutating=false)
			end

			# Define the Jacobian as a matrix
			# to be used with direct solver
			if method ==3
				FFT = exp.(-1im*k*(x.-x₀)');
				IFFT = exp.(1im*k*(x.-x₀)')/length(x);
				Id = Diagonal(ones(size(x)));
				Mean = ones(size(x))'/length(x)
				M(v) = Diagonal( v )

				function JacFMat( u )
					dx.=real.(ifft(u))
					fη .= fft(data.η(x+ϵ/ν*dx))
					δ  = meanf(fη)
		        	return Id + ϵ*sqrt(μ)*1im*Π⅔ .* (M(cotanh(sqrt(μ)*(1+ϵ*δ)*k)) * FFT * M(dη(x+ϵ/ν*dx)) * IFFT
								+ ϵ/(1+ϵ*δ)* (xdcotanh(sqrt(μ)*(1+ϵ*δ)*k).* fη )  * Mean * M(dη(x+ϵ/ν*dx)) * IFFT )
				end
			end

			# The iterative map to solve F(u)=u
			function iterate( u )
				if method==1  		# by contraction fix point algorithm
			    	return -F(u)
				elseif method==2	# by Newton algorithm with GMRES iterative solver to invert the Jacobian
					return gmres( JacF(u) , -F(u) ; reltol = tol/100, verbose=verbose )
				elseif method==3	# by Newton algorithm with direct solver to invert the Jacobian
					return - JacFMat( u ) \ F(u)
				else
					@error("In the function `WaterWaves`, the argument `method` should be 1, 2 or 3.")
				end

			end

			# initiate the iterative argument
			fη.=fft(data.η(x))
			δ = meanf(fη)
			u0 .= -ν*sqrt(μ)*1im*Π⅔ .*cotanh(sqrt(μ)*k*(1+ϵ*δ)).*fη

			# perform the iterative loop
			norm0=norm(fη)
			normdiff=norm0
			niter=0
			while  normdiff>tol*norm0 && niter<maxiter
				nu0 = iterate(u0)
				u0 += nu0
				normdiff=norm(nu0)
			    niter+=1
				if verbose
					@info "Relative error  $(normdiff/norm0)  at step  $niter"
				end
			end
			if verbose
				if niter == maxiter
				    @warn "The iterative solver did not converge. Maybe try a different method."
				    @warn "Estimated normalized error : $(normdiff/norm0)"
				else
					@info "The iterative solver converged in $niter iterations."
					@info "Estimated normalized error : $(normdiff/norm0)"
				end
			end

			# Constructs relevant variables from u0 the solution to F(u)=0
			z0 .= data.η(x+ϵ/ν*real.(ifft(u0)))
			v0 .= data.v(x+ϵ/ν*real.(ifft(u0))) .* (1 .+ ϵ/ν*real.(ifft(Π⅔.*∂ₓ.*u0)) )

			U=[ Π⅔ .* fft(z0)   Π⅔ .* fft(v0) ]  # Π⅔ for dealiasing
			U[ abs.(U).< ktol ].=0     # applies Krasny filter if ktol>0
 			return U
		  end
		end

		# Reconstruct physical variables from raw data
		# Return `(η,v,x)`, where
		# - `η` is the surface deformation;
		# - `v` is the derivative of the trace of the velocity potential;
		# - `x` is the vector of collocation points
		function mapfro(U)
		   	ξ  .= real.(sqrt(μ)*(1+ϵ*meanf(U[:,1]))*k)
	       	xv .= real.(-1im*sqrt(μ)*ifft( Π⅔ .*cotanh(ξ) .*  U[:,1] ))

		   return real.( ifft(U[:,1]) ) , real.( ifft(U[:,2])./(1 .+ ϵ*real.(ifft(∂ₓ.*fft(xv)) )) ), x + ϵ*xv
		end

		# Water Waves equations are ∂t U = f(U)
		function f!(U)
		   	Dphi .= real.(ifft(U[:,2]))
			ξ .= sqrt(μ)*(1 .+ ϵ*meanf(U[:,1]))*k
			Dxv .= real.(sqrt(μ)*ifft( k .* cotanh(ξ) .* U[:,1]  ))
			Dz .= real.(ifft( ∂ₓ.*  U[:,1] ))

			J .= (1 .+ ϵ*Dxv).^2 + μ*(ϵ*Dz).^2
			M1 .= real.(-1im/sqrt(μ)*ifft( mytanh(ξ).* U[:,2] ))
			M2 .= real.( 1im*sqrt(μ)*ifft( cotanh(ξ) .* fft(M1./J) ))
			q0 = mean((1 .+ ϵ*Dxv).*M2 + ϵ*μ*Dz.*M1./J)

			U[:,2] .= -∂ₓ .* U[:,1] - ϵ* Π⅔ .* ∂ₓ .* fft( Dphi.*M2 + 1/2 .*(Dphi.^2 -μ*M1.^2)./J - q0*Dphi )/ν
			U[:,1] .=  Π⅔ .* fft( (1 .+ ϵ*Dxv).*M1./J - ϵ*Dz.*M2 + ϵ*q0*Dz )/ν

		end

		# Water waves equations are ∂t (U1,U2) = (f1(U1,U2) , f2(U1,U2))
		function f1!(U1,U2)
		   	Dphi .= real.(ifft(U2))
			ξ .= sqrt(μ)*(1 .+ ϵ*meanf(U1))*k
			Dxv .= real.(sqrt(μ)*ifft( k .* cotanh(ξ) .* U1 ))
			Dz .= real.(ifft( ∂ₓ.* U1))

			J .= (1 .+ ϵ*Dxv).^2 + μ*(ϵ*Dz).^2
			M1 = real.(-1im/sqrt(μ)*ifft( mytanh(ξ).* U2 ))
			M2 = real.( 1im*sqrt(μ)*ifft( cotanh(ξ) .* fft(M1./J) ))
			q0 = mean((1 .+ ϵ*Dxv).*M2 + ϵ*μ*Dz.*M1./J)

			U1 .= Π⅔ .* fft( (1 .+ ϵ*Dxv).*M1./J - ϵ*Dz.*M2 + ϵ*q0*Dz )/ν
		end
		function f2!(U1,U2)
		   	Dphi .= real.(ifft(U2))
			ξ .= sqrt(μ)*(1 .+ ϵ*meanf(U1))*k
			Dxv .= real.(sqrt(μ)*ifft( k .* cotanh(ξ) .* U1))
			Dz .= real.(ifft( ∂ₓ.* U1))

			J .= (1 .+ ϵ*Dxv).^2 + μ*(ϵ*Dz).^2
			M1 = real.(-1im/sqrt(μ)*ifft( mytanh(ξ).* U2 ))
			M2 = real.( 1im*sqrt(μ)*ifft( cotanh(ξ) .* fft(M1./J) ))
			q0 = mean((1 .+ ϵ*Dxv).*M2 + ϵ*μ*Dz.*M1./J)

			U2 .= -∂ₓ .* U1 - ϵ* Π⅔ .* ∂ₓ .* fft( Dphi.*M2 + 1/2 .*(Dphi.^2 -μ*M1.^2)./J - q0*Dphi )/ν
		end

		new(label, f!, f1!, f2!, mapto, mapfro, info )
    end
end


text:
## Ref values for test POAP Tests
c_dirac_objs = [0.150401, 0.15072, 0.150806, 0.151217, 0.15157, 0.152157, 0.152176, 0.152169, 0.152273, 0.152625, 0.15326, 0.153586, 0.154105, 0.153971, 0.155061, 0.156201, 0.156065, 0.157107, 0.158064, 0.15822, 0.159576, 0.158978, 0.160406, 0.161868, 0.160384, 0.159505, 0.161051, 0.159079, 0.158525, 0.158969, 0.158921, 0.160585, 0.159371, 0.161072, 0.16279, 0.164225, 0.164118, 0.163612, 0.165383, 0.167163, 0.168953, 0.167443, 0.169247, 0.17007, 0.170889, 0.169784, 0.17123, 0.169543, 0.17058, 0.172401, 0.172817, 0.171209, 0.173027, 0.174843, 0.176656, 0.178464, 0.179831, 0.176824, 0.176724, 0.176429, 0.170939, 0.172716, 0.165667, 0.167429, 0.169183, 0.17093, 0.172668, 0.1699, 0.171621, 0.173333, 0.172384, 0.171649, 0.173333, 0.16552, 0.165682, 0.162314, 0.163958, 0.162344, 0.16006, 0.161673, 0.163275, 0.164867, 0.166448, 0.164142, 0.160528, 0.157595, 0.159133, 0.158655, 0.159929, 0.158477, 0.158735, 0.160253, 0.161726, 0.163188, 0.163358, 0.158799, 0.154889, 0.152338, 0.153746, 0.155143, 0.156529, 0.157906, 0.159272, 0.150353, 0.151764, 0.153164, 0.154555, 0.155935, 0.15458, 0.153801, 0.148818, 0.146722, 0.148052, 0.149373, 0.147037, 0.142586, 0.142068, 0.14335, 0.144623, 0.145887, 0.147142, 0.148387, 0.144226, 0.14174, 0.142958, 0.144168, 0.145369, 0.143739, 0.144859, 0.145971, 0.147074, 0.148168, 0.149255, 0.150333, 0.151403, 0.152465, 0.152207, 0.151604, 0.152643, 0.153673, 0.154696, 0.154209, 0.155216, 0.156186, 0.154538, 0.155523, 0.156501, 0.157472, 0.153079, 0.154036, 0.153504, 0.154446, 0.155382, 0.156311, 0.157233, 0.158148, 0.159056, 0.159888, 0.160784, 0.161673, 0.162555, 0.158035, 0.158904, 0.1598, 0.156165, 0.157076, 0.157981, 0.158879, 0.159771, 0.160657, 0.161538, 0.162412, 0.16328, 0.154548, 0.155405, 0.156256, 0.157101, 0.157941, 0.158775, 0.159604, 0.158697, 0.156979, 0.157791, 0.158598, 0.1594, 0.15829, 0.159082, 0.156014, 0.156795, 0.157571, 0.158342, 0.156745, 0.157433, 0.154931, 0.155627, 0.156317, 0.157003, 0.157684, 0.158361, 0.159033, 0.159678, 0.160341, 0.161, 0.157912, 0.158563, 0.159209, 0.157959, 0.154465, 0.148083, 0.148712, 0.149337, 0.149958, 0.150575, 0.151188, 0.151797, 0.152402, 0.153003, 0.153509, 0.154102, 0.154692, 0.155278, 0.15586, 0.151213, 0.151788, 0.152359, 0.152927, 0.153491, 0.154051, 0.148745, 0.14926, 0.14981, 0.150356, 0.150899, 0.151439, 0.151975, 0.152509, 0.153038, 0.153565, 0.154088, 0.154608, 0.155125, 0.155639, 0.15615, 0.156657, 0.157162, 0.157663, 0.158162, 0.158657, 0.159149, 0.159639, 0.160125, 0.160609, 0.16109, 0.161568, 0.162043, 0.162515, 0.162984, 0.163451, 0.163915, 0.164376, 0.164835, 0.16529, 0.165744, 0.165026, 0.163336, 0.163782, 0.164224, 0.164665, 0.165102, 0.165537, 0.16597, 0.164885, 0.165219, 0.163283, 0.163705, 0.164126, 0.164544, 0.16496, 0.165373, 0.165784, 0.16612, 0.166526, 0.16693, 0.167332, 0.164788, 0.165186, 0.165581, 0.165974, 0.166365, 0.166753, 0.16714, 0.167524, 0.167907, 0.168287, 0.168665, 0.169041, 0.169415, 0.169787, 0.170157, 0.170525, 0.170891, 0.171255, 0.171617, 0.171977, 0.172336, 0.172692, 0.173046, 0.173399, 0.17375, 0.174099, 0.174446, 0.174791, 0.175134, 0.175476, 0.175815, 0.176153, 0.174177, 0.174511, 0.174844, 0.175175, 0.175504, 0.175832, 0.176158, 0.176482, 0.176805, 0.177126, 0.177445, 0.177763, 0.178079, 0.178393, 0.178706, 0.179017, 0.179327, 0.179635, 0.179942, 0.180247, 0.18055, 0.180853, 0.181153, 0.181407, 0.181705, 0.182001, 0.182295, 0.182589, 0.18288, 0.183171, 0.183459, 0.183747, 0.184033, 0.184318, 0.184601, 0.184883, 0.185163, 0.185443, 0.18572, 0.185997, 0.184566, 0.18484, 0.185065, 0.185336, 0.185606, 0.185875, 0.186142, 0.186322, 0.186587, 0.186851, 0.187113, 0.187374, 0.187634, 0.187893, 0.188089, 0.188345, 0.1886, 0.188854, 0.189107, 0.189359, 0.189609, 0.189858, 0.190107, 0.190353, 0.190599, 0.190844, 0.191088, 0.19133, 0.191572, 0.191812, 0.192051, 0.192289, 0.192526, 0.192762, 0.192997, 0.193231, 0.193464, 0.193696, 0.193926, 0.194156, 0.194385, 0.192969, 0.188519, 0.188744, 0.188969, 0.189193, 0.189415, 0.189637, 0.189858, 0.190078, 0.190297, 0.190515, 0.190732, 0.190948, 0.191163, 0.191377, 0.19159, 0.191802, 0.192014, 0.192224, 0.192434, 0.192521, 0.192729, 0.192936, 0.193142, 0.193301, 0.193506, 0.193709, 0.193912, 0.194113, 0.194314, 0.194514, 0.190315, 0.190513, 0.190711, 0.190907, 0.191103, 0.191298, 0.191493, 0.191686, 0.191879, 0.192071, 0.192262, 0.192452, 0.192642, 0.19283, 0.193018, 0.193206, 0.193392, 0.193578, 0.193763, 0.193947, 0.194131, 0.190433, 0.190615, 0.190796, 0.190977, 0.191157, 0.191336, 0.191514, 0.191692, 0.191869, 0.192045, 0.186506, 0.18668, 0.186759, 0.186933, 0.187105, 0.187278, 0.182689, 0.18286, 0.18303, 0.183199, 0.183368, 0.183536, 0.183704, 0.183871, 0.184037, 0.184202, 0.184367, 0.184532, 0.184695, 0.184858, 0.185021, 0.185182, 0.185344, 0.185504, 0.185664, 0.185823, 0.185982, 0.18614, 0.186298, 0.186455, 0.186611, 0.186767, 0.186922, 0.187077, 0.187231, 0.187384, 0.183704, 0.183856, 0.184008, 0.184159, 0.18431, 0.18446]
c_boxStein_objs = [0.23616115868423004,0.23025212639816284,0.24660937216515688,0.21514684415891228,0.22196111668092788,0.21946103469163075,0.23997446159033486,0.2377587301912953,0.23371997352923735,0.19857990846373003,0.20182456885781075,0.21310235040503722,0.18582168394868576,0.1966111821432121,0.19926635135137688,0.17592562876470824,0.17367368288694812,0.16504454066170676,0.17874162126029552,0.20022185286680982,0.19140400357368625,0.1878961030360683,0.18707313655098284,0.1742853850671333,0.19753325645556363,0.1976770702283252,0.18558476330782891,0.17529072137060753,0.18703526051680985,0.17305800949827926,0.17602412126998124,0.1902397029928109,0.19997361047737489,0.1954569177512014,0.2009068943292797,0.18695848261987302,0.19012037664221298,0.19561543359284972,0.19231034622379187,0.2009178095200156,0.20089460107467977,0.19666218220966386,0.19832943445434523,0.20323135918635551,0.20804200659568928,0.20232616740607778,0.20132622674865433,0.19642430201664413,0.20132622674865433,0.21162026868587594,0.22142411814989657,0.2219143106230976,0.2258358504087058,0.22253032234909942,0.22743224708110982,0.23816780047918198,0.24789907649841803,0.2449553350846665,0.2400534103526562,0.2351514856206459,0.23515148562064586,0.23515148562064586,0.23515148562064586,0.23515148562064586,0.2323483256931152,0.2328385181663163,0.23185813321991425,0.2264660160147029,0.22593085417325975,0.22593085417325975,0.22151912191445045,0.21661719718244007,0.21976544565894934,0.21967912895555686,0.22010814931051642,0.22010814931051642,0.22501007404252676,0.22451988156932567,0.21912776436411427,0.21912776436411427,0.20790515515275032,0.20790515515275032,0.20790515515275032,0.20790515515275032,0.20778412457288553,0.20773810157693995,0.20333917675735086,0.20333917675735086,0.20333917675735086,0.2033391767573508,0.20387067029609682,0.2087725950281072,0.20305755571553002,0.19771194045196466,0.19281001571995435,0.19281001571995435,0.19281001571995435,0.19281001571995435,0.19281001571995435,0.18897704807236212,0.17868300613514043]
c_oracleStein_objs = [0.157011, 0.157239, 0.157476, 0.157682, 0.157776, 0.158003, 0.158026, 0.158199, 0.158355, 0.158482, 0.158991, 0.159125, 0.159302, 0.159418, 0.159522, 0.159662, 0.159722, 0.159727, 0.15977, 0.15987, 0.15992, 0.159992, 0.159993, 0.160318, 0.160401, 0.16044, 0.160678, 0.160796, 0.160884, 0.16091, 0.160949, 0.161083, 0.161107, 0.161124, 0.161128, 0.161403, 0.161559, 0.161602, 0.162013, 0.162078, 0.162144, 0.162203, 0.162353, 0.162478, 0.162589, 0.162714, 0.162752, 0.162872, 0.162901, 0.163, 0.163061, 0.163261, 0.163295, 0.163366, 0.163426, 0.16352, 0.163532, 0.163671, 0.163718, 0.163921, 0.16407, 0.164134, 0.16427, 0.164532, 0.164577, 0.164942, 0.165063, 0.165222, 0.165325, 0.165546, 0.165693, 0.165904, 0.166125, 0.16634, 0.166608, 0.166787, 0.166926, 0.166989, 0.167213, 0.167445, 0.167584, 0.167838, 0.167996, 0.168118, 0.16821, 0.168387, 0.168506, 0.168652, 0.168845, 0.169012, 0.169018, 0.169071, 0.169156, 0.169073, 0.169125, 0.169186, 0.169265, 0.169445, 0.16949, 0.16953, 0.169696, 0.169757, 0.169759, 0.169825, 0.169868, 0.169794, 0.16987, 0.169933, 0.169942, 0.170006, 0.170049, 0.170086, 0.170114, 0.170175, 0.170281, 0.170338, 0.170294, 0.170367, 0.170438, 0.170514, 0.170542, 0.170496, 0.170559, 0.170671, 0.170703, 0.17061, 0.170656, 0.170745, 0.170712, 0.170819, 0.170863, 0.171018, 0.171065, 0.17116, 0.171213, 0.171276, 0.17134, 0.171101, 0.171263, 0.171129, 0.171244, 0.171322, 0.171447, 0.171597, 0.171359, 0.170962, 0.171095, 0.171248, 0.171025, 0.171078, 0.171147, 0.171205, 0.170761, 0.170835, 0.170534, 0.170586, 0.170642, 0.170448, 0.170486, 0.170317, 0.169951, 0.169729, 0.16977, 0.169893, 0.169516, 0.169322, 0.169425, 0.169473, 0.169567, 0.16928, 0.169391, 0.169462, 0.169509, 0.169606, 0.169701, 0.169726, 0.169442, 0.169008, 0.168995, 0.16905, 0.1691, 0.169037, 0.169106, 0.169137, 0.169077, 0.169059, 0.169089, 0.169093, 0.169092, 0.168833, 0.168916, 0.168787, 0.168753, 0.168736, 0.168756, 0.168765, 0.168759, 0.168559, 0.168303, 0.168304, 0.168328, 0.168281, 0.168214, 0.168103, 0.168038, 0.168018, 0.168035, 0.16799, 0.167954, 0.167626, 0.167496, 0.167514, 0.167547, 0.167595, 0.167444, 0.167186, 0.166938, 0.166836, 0.166423, 0.166302, 0.166073, 0.165774, 0.16563, 0.165657, 0.165498, 0.16548, 0.165449, 0.165483, 0.165416, 0.165402, 0.165239, 0.165306, 0.165076, 0.165024, 0.164909, 0.164864, 0.16474]
c_oracleReg_objs = [0.158126, 0.159274, 0.159761, 0.160121, 0.16058, 0.161054, 0.161829, 0.162641, 0.163409, 0.164157, 0.164823, 0.165457, 0.166086, 0.166643, 0.166982, 0.167315, 0.167636, 0.167957, 0.168263, 0.16854, 0.168791, 0.169041, 0.169302, 0.16958, 0.169854, 0.170103, 0.170334, 0.170546, 0.170738, 0.170917, 0.171088, 0.17125, 0.171401, 0.171545, 0.171713, 0.171884, 0.172047, 0.172202, 0.172353, 0.172497, 0.172636, 0.172767, 0.172893, 0.173012, 0.173125, 0.173237, 0.173344, 0.17345, 0.173552, 0.173652]
c_SteinReg_objs = [-0.0827889, 0.146248, 0.262713, 0.280625, 0.264368]
c_LOOReg_objs = [0.202822, 0.200359, 0.198748, 0.197225, 0.1987, 0.200685, 0.201104, 0.202436, 0.20414, 0.204973, 0.204724, 0.204304, 0.204596, 0.204852, 0.204625, 0.203521, 0.201984, 0.200353, 0.198855, 0.197438, 0.196068, 0.19482, 0.193879, 0.193049, 0.192295, 0.191564, 0.190859, 0.190214, 0.189622, 0.18908, 0.188561, 0.188057, 0.187582, 0.187133, 0.186723, 0.186345, 0.185983, 0.185641, 0.185318, 0.185003, 0.184709, 0.184443, 0.184207, 0.184002, 0.183799, 0.183615, 0.183453, 0.183295, 0.183164, 0.183036]

ext: jl


text:
"""
    set_level_and_version(level, version)

A converter to pass into [`readSBML`](@ref) that enforces certain SBML level
and version.
"""
set_level_and_version(level, version) =
    doc -> begin
        ccall(
            sbml(:SBMLDocument_setLevelAndVersion),
            Cint,
            (VPtr, Cint, Cint),
            doc,
            level,
            version,
        )
    end

"""
    libsbml_convert(conversion_options::Vector{Pair{String, Dict{String, String}}})

A converter that runs the SBML conversion routine, with specified conversion
options. The argument is a vector of pairs to allow specifying the order of
conversions.
"""
libsbml_convert(conversion_options::Vector{Pair{String,Dict{String,String}}}) =
    doc -> begin
        for (converter, options) in conversion_options
            props = ccall(sbml(:ConversionProperties_create), VPtr, ())
            opt = ccall(sbml(:ConversionOption_create), VPtr, (Cstring,), converter)
            ccall(sbml(:ConversionProperties_addOption), Cvoid, (VPtr, VPtr), props, opt)
            for (k, v) in options
                opt = ccall(sbml(:ConversionOption_create), VPtr, (Cstring,), k)
                ccall(sbml(:ConversionOption_setValue), Cvoid, (VPtr, Cstring), opt, v)
                ccall(
                    sbml(:ConversionProperties_addOption),
                    Cvoid,
                    (VPtr, VPtr),
                    props,
                    opt,
                )
            end
            ccall(sbml(:SBMLDocument_convert), Cint, (VPtr, VPtr), doc, props)
        end
    end

"""
    libsbml_convert(converter::String; kwargs...)

Quickly construct a single run of a `libsbml` converter from keyword arguments.

# Example
```
readSBML("example.xml", libsbml_convert("stripPackage", package="layout"))
```
"""
libsbml_convert(converter::String; kwargs...) =
    libsbml_convert([converter => Dict(string(k) => string(v) for (k, v) in kwargs)])

"""
    convert_simplify_math

Shortcut for [`libsbml_convert`](@ref) that expands functions, local
parameters, and initial assignments in the SBML document.
"""
convert_simplify_math = libsbml_convert(
    ["promoteLocalParameters", "expandFunctionDefinitions"] .=>
        Ref(Dict{String,String}()),
)


text:
using Base.Test
using Musical
using Musical.Notations # A
using Musical.SPN # note octave freq A4
using MIDI # MIDI.Note

@test A isa Musical.Notations.Notation
@test A4 isa Musical.SPN.NoteOctave
@test note(A4) == A
@test octave(A4) == 4
@test freq(A4) == 440.0
@test MIDI.Note(A4) == MIDI.Note(69, 0, 0, 0)


text:
#===============================
<-----------Lathe.jl----------->
Programmed by Emmett Boudreau
    <emmett@emmettboudreau.com>
        <http://emmettboudreau.com>
MIT General Open Source License
    (V 3.0.0)
        Free for Modification and
        Redistribution
=~=~=~=~=~=~=~=~=~=~=~=~=~=~=~=|
         CONTRIBUTORS
=~=~=~=~=~=~=~=~=~=~=~=~=~=~=~=|
        ~ emmettgb
        ~ stefanches7
/><><><><><><><><><><><><><><><><\
Thank you for your forks!
<-----------Lathe.jl----------->
38d8eb38-e7b1-11e9-0012-376b6c802672
#[deps]
DataFrames.jl
Random.jl
================================#
@doc """
      |====== Lathe - Easily ML =====\n
      |= = = = = v. 0.1.0 = = = = = |\n
      |==============================\n
      |__________Lathe.stats\n
      |__________Lathe.validate\n
      |__________Lathe.preprocess\n
      |__________Lathe.models\n
      |__________Lathe.data\n
      |______________________________\n
      Use ?(Lathe.package) for information!\n
      [uuid]\n
      38d8eb38-e7b1-11e9-0012-376b6c802672\n
      [deps]\n
      DataFrames.jl\n
      Random.jl\n
       """
module Lathe
# <------- PARTS ----->
include("Stats.jl")
include("Models.jl")
include("Preprocess.jl")
include("Data.jl")
include("Distributions.jl")
# <------- PARTS ----->
# <------- DEPS ----->
using DataFrames
using Random
# <------- DEPS ----->

#==
This is the end of the main
module, nothing is to be written
beyond here
==#
end


text:
"""
    gRNA_frequency_distribution(m,
                                sd,
                                l,
                                u,
                                n_gRNA_total;
                                normalize = true,
                                visualize = false)

Generates vector with frequencies in the combinatorial gRNA/Cas9 construct library for all gRNAs

***INPUT***
m: the average abundance of the gRNAs (in terms of absolute or relative frequency)
sd: the standard deviation on the gRNA abundances (in terms of absolute or relative frequency)
l: minimal gRNA abundance (in terms of absolute or relative frequency)
u: maximal gRNA abundance (in terms of absolute or relative frequency)
n_gRNA_total: the total number of gRNAs in the experiment
normalize: if set to "true", the gRNA abundances (absolute frequencies) are converted into relative frequencies
visualize: if set to "true", a histogram of all gRNA abundances is plotted

***OUTPUT***
p_gRNA_freq: vector with frequencies for all gRNAs in the construct library
"""
function gRNA_frequency_distribution(m::Real,
                                    sd::Real,
                                    l::Real,
                                    u::Real,
                                    n_gRNA_total::Int64;
                                    normalize = true,
                                    visualize = false)

    ### constraints on parameters
    @assert m > 0 && sd > 0 && l > 0 && u > 0
    @assert l < m < u
    @assert n_gRNA_total > 0

    ### generate distribution
    d_gRNA_freq = truncated(Normal(m, sd), l, u)  # gRNA frequency distribution
    p_gRNA_freq = collect(rand(d_gRNA_freq, n_gRNA_total))  # sample gRNA frequencies from distribution

    if normalize # convert into relative frequencies
        p_gRNA_freq /= sum(p_gRNA_freq)
    end

    if visualize
        return histogram(p_gRNA_freq, label="",
            xlabel="Number of reads per gRNA",
            linecolor="white",
            normalize=:probability,
            xtickfontsize=10,ytickfontsize=10,
            color=:mediumturquoise, size=(600,350), bins = 25,
            ylabel="Relative frequency",
            title="gRNA frequency distribution")
    else
        return p_gRNA_freq
    end
end

"""
    gRNA_edit_distribution(f_act,
                            ϵ_edit_act,
                            ϵ_edit_inact,
                            sd_act,
                            n_gRNA_total;
                            visualize = false)

Generates vector with genome editing efficiencies for all the gRNAs in the experiment.

***INPUT***
f_act: fraction of all gRNAs that is active
ϵ_edit_act: Average genome editing efficiency for active gRNAs - mean of the genome editing efficiency distribution for active gRNAs
ϵ_edit_inact: Average genome editing efficiency for inactive gRNAs - mean of the genome editing efficiency distribution for inactive gRNAs
sd_act: standard deviation of the genome editing efficiency distributions for active and inactive gRNAs
n_gRNA_total: the total number of gRNAs in the experiment
visualize: if set to "true", a histogram of all genome editing efficiency is plotted

***OUTPUT***
p_gRNA_edit: vector with genome editing efficiencies for all gRNAs
"""
function gRNA_edit_distribution(f_act::Real,
                                ϵ_edit_act::Real,
                                ϵ_edit_inact::Float64,
                                sd_act::Float64,
                                n_gRNA_total::Int64;
                                visualize=false)

    ### constraints on parameters
    @assert 0 < f_act <= 1
    @assert 0.5 < ϵ_edit_act <= 1
    @assert 0 < ϵ_edit_inact <= 0.5
    @assert sd_act > 0
    @assert n_gRNA_total > 0

    ### generate distribution
    d_act = Binomial(1, f_act) # there is a probability f_act that a gRNA is active
    d_high_act = truncated(Normal(ϵ_edit_act, sd_act), 0.01, 1)  # average genome editing efficiency for active gRNAs is equal to ϵ_edit_act
    d_low_act = truncated(Normal(ϵ_edit_inact, sd_act), 0.01, 1) # average genome editing efficiency for inactive gRNAs is equal to ϵ_edit_inact
    p_gRNA_edit = zeros(n_gRNA_total) # initialize vector with genome editing efficiencies for gRNAs

    for i in 1:n_gRNA_total
        if rand(d_act, 1) == [1]  # gRNA is active
            p_gRNA_edit[i] = rand(d_high_act, 1)[1]
        else  # gRNA is inactive
            p_gRNA_edit[i] = rand(d_low_act, 1)[1]
        end
    end

    if visualize
        return histogram(p_gRNA_edit,
                normalize = :probability,
                linecolor = "white",
                label="",
                color=:turquoise4,
                xtickfontsize=10,ytickfontsize=10, xlim = (0, 1),
                xticks=(0:0.1:1),
                bins = 150,
                xlabel="gRNA editing efficiency",
                ylabel="Relative frequency",
                title="gRNA genome editing effiency distribution")

    else
        return p_gRNA_edit
    end
end

"""
    simulate_Nₓ₁(x,
                    g,
                    r,
                    n_gRNA_total,
                    p_gRNA_freq,
                    p_gRNA_edit,
                    ϵ_KO;
                    iter = 500)

Computes the expected value and the standard deviation of the minimal plant library size for full coverage of all single gene knockouts (E[Nx,1] and σ[Nx,1]) using simulation

***INPUT***
x: number of target genes in the experiment
g: number of gRNAs designed per target gene
r: number of gRNA sequences per combinatorial gRNA/Cas construct (>= 1)
n_gRNA_total: total number of gRNAs in the experiment (should be equal to x * g)
p_gRNA_freq: vector with relative frequencies for all gRNAs in the construct library (normalized!)
p_gRNA_edit: vector with genome editing efficiencies for all gRNAs (each genome editing efficiency in ]0,1])
ϵ_KO: global knockout efficiency; fraction of mutations leading to effective gene knockout (]0, 1])
iter: number of CRISPR/Cas experiments that are simulated to obtain E[Nₓ₁] and σ[Nₓ₁]

***OUTPUT***
E_Nₓ₁: expected value of the plant library size for full coverage of all single gene knockouts
sd_Nₓ₁: standard deviation on the plant library size for full coverage of all single gene knockouts
"""
function simulate_Nₓ₁(x::Int64,
                    g::Int64,
                    r::Int64,
                    n_gRNA_total::Int64,
                    p_gRNA_freq::Vector{Float64},
                    p_gRNA_edit::Vector{Float64},
                    ϵ_KO::Real;
                    iter = 500)

    ### constraints on parameters
    @assert x > 0 && g > 0
    @assert x * g == n_gRNA_total
    @assert r >= 1
    @assert 1 - 1e-4 < sum(p_gRNA_freq) < 1 + 1e4
    @assert length(p_gRNA_freq) == n_gRNA_total
    @assert all(0 .< p_gRNA_edit .<= 1)
    @assert length(p_gRNA_edit) == n_gRNA_total
    @assert 0 < ϵ_KO <= 1
    @assert iter > 1

    Nₓ₁_vec = [] #stores number of plants to reach full coverage for each simulated experiment
        for i in 1:iter
            genes_vec = [] # Initialize vector to store single gene knockouts that are observed in plants
            Nₓ₁ = 0
            while genes_vec != collect(1:x) # check if all possible single gene knockouts are present: if no full coverage, sample an additional plant
                Nₓ₁ += 1 # count how many plants must be sampled to observe all single gene knockouts

                # sample combinatorial gRNA/Cas9 construct
                gRNA_indices_construct = findall((rand(Multinomial(r, p_gRNA_freq))) .!= 0)

                # execute mutations
                gRNA_indices_mutations = [gRNA for gRNA in gRNA_indices_construct if rand(Binomial(1, p_gRNA_edit[gRNA])) == 1]

                # effective gene knockout (loss-of-function) ?
                gRNA_indices_KO = [gRNA for gRNA in gRNA_indices_mutations if rand(Binomial(1, ϵ_KO)) == 1]

                # which genes are knocked out?
                genes_indices_KO = Int.(ceil.(gRNA_indices_KO / g))
                append!(genes_vec, genes_indices_KO)

                # update vector with observed gene knockouts
                genes_vec = Int.(sort(unique(genes_vec)))
            end
            push!(Nₓ₁_vec, Nₓ₁) # add plant library size for full coverage of current experiment to vector
        end

    # Calculate expected value and standard deviation
    E_Nₓ₁ = mean(Nₓ₁_vec);
    sd_Nₓ₁ = std(Nₓ₁_vec)

    return E_Nₓ₁, sd_Nₓ₁
end

"""
    BioCCP_Nₓ₁(x,
                g,
                r,
                n_gRNA_total,
                p_gRNA_freq,
                p_gRNA_edit,
                ϵ_KO)

Computes the expected value and the standard deviation of the minimal plant library size for
full coverage of all single gene knockouts (E[Nx,1] and σ[Nx,1]) using BioCCP

***INPUT***
x: number of target genes in the experiment
g: number of gRNAs designed per target gene
r: number of gRNA sequences per combinatorial gRNA/Cas construct
n_gRNA_total: total number of gRNAs in the experiment
p_gRNA_freq: vector with relative frequencies for all gRNAs in the construct library (normalized!)
p_gRNA_edit: vector with genome editing efficiencies for all gRNAs
ϵ_KO: global knockout efficiency; fraction of mutations leading to effective gene knockout

***OUTPUT***
E_Nₓ₁ : expected value of the plant library size for full coverage of all single gene knockouts
sd_Nₓ₁ : standard deviation on the plant library size for full coverage of all single gene knockouts
"""
function BioCCP_Nₓ₁(x::Int64,
                    g::Int64,
                    r::Int64,
                    n_gRNA_total::Int64,
                    p_gRNA_freq::Vector{Float64},
                    p_gRNA_edit::Vector{Float64},
                    ϵ_KO::Real)

    ### constraints on parameters
    @assert x > 0 && g > 0
    @assert x * g == n_gRNA_total
    @assert r >= 1
    @assert 1 - 1e-4 < sum(p_gRNA_freq) < 1 + 1e4
    @assert length(p_gRNA_freq) == n_gRNA_total
    @assert all(0 .< p_gRNA_edit .<= 1)
    @assert length(p_gRNA_edit) == n_gRNA_total
    @assert 0 < ϵ_KO <= 1

    # prepare input for BioCCP functions
    p_gRNAs = p_gRNA_freq .* p_gRNA_edit * ϵ_KO  # calculate probability for each gRNA to induce effective gene knockout
    p_genes = [sum(p_gRNAs[i:i+g-1]) for i in 1:g:n_gRNA_total]  # obtain probability of single gene knockout by summing up probability of all corresponding gRNAs to induce effective gene knockout

    # Apply BioCCP functions
    E_Nₓ₁ = expectation_minsamplesize(x; p=p_genes, r=r, normalize=false)
    sd_Nₓ₁ = std_minsamplesize(x; p=p_genes, r=r, normalize=false)

    return E_Nₓ₁, sd_Nₓ₁
end

"""
    simulate_Nₓ₂(x,
                    g,
                    r,
                    n_gRNA_total,
                    p_gRNA_freq,
                    p_gRNA_edit,
                    ϵ_KO;
                    iter=500)

Computes the expected value and the standard deviation of the minimal plant library size for full coverage of all pairwise combinations of gene knockouts
in a multiplex CRISPR/Cas experiment (E[Nx,2] and σ[Nx,2]) using simulation

***INPUT***
x: number of target genes in the experiment
g: number of gRNAs designed per target gene
r: number of gRNA sequences per combinatorial gRNA/Cas construct
n_gRNA_total: total number of gRNAs in the experiment
p_gRNA_freq: vector with relative frequencies for all gRNAs in the construct library (normalized!)
p_gRNA_edit: vector with genome editing efficiencies for all gRNAs
ϵ_KO: global knockout efficiency; fraction of mutations leading to effective gene knockout
iter: number of CRISPR/Cas experiments that are simulated to obtain E[Nₓ₂] and σ[Nₓ₂]

***OUTPUT***
E_Nₓ₂ : expected value of the plant library size for full coverage of all pairwise combinations of gene knockouts
sd_Nₓ₂ : standard deviation on the plant library size for full coverage of all pairwise combinations of gene knockouts
"""
function simulate_Nₓ₂(x::Int64,
                    g::Int64,
                    r::Int64,
                    n_gRNA_total::Int64,
                    p_gRNA_freq::Vector{Float64},
                    p_gRNA_edit::Vector{Float64},
                    ϵ_KO::Real;
                    iter = 500)

    ### constraints on parameters
    @assert x > 0 && g > 0
    @assert x * g == n_gRNA_total
    @assert r >= 2
    @assert 1 - 1e-4 < sum(p_gRNA_freq) < 1 + 1e4
    @assert length(p_gRNA_freq) == n_gRNA_total
    @assert all(0 .< p_gRNA_edit .<= 1)
    @assert length(p_gRNA_edit) == n_gRNA_total
    @assert 0 < ϵ_KO <= 1
    @assert iter > 1

    Nₓ₂_vec = [] #stores number of plants to reach full coverage of all pairwise combinations of gene knockouts for each simulated experiment

    for i in 1:iter
        X_interactions_count = zeros(x, x) # Initialize matrix to count pairwise interactions
        Nₓ₂ = 0
        while X_interactions_count != ones(x, x) # check if all pairwise combinations are present
            Nₓ₂ += 1 # count how many plants must be sampled to fill pairwise interaction matrix

            # sample combinatorial gRNA/Cas9 construct
            gRNA_indices_construct = findall((rand(Multinomial(r, p_gRNA_freq))) .!= 0)

            # execute mutations
            gRNA_indices_mutations = [gRNA for gRNA in gRNA_indices_construct if rand(Binomial(1, p_gRNA_edit[gRNA])) == 1]

            # effective gene knockout (loss-of-function) ?
            gRNA_indices_KO = [gRNA for gRNA in gRNA_indices_mutations if rand(Binomial(1, ϵ_KO)) == 1]

            # which genes are knocked out?
            genes_indices_KO = Int.(ceil.(gRNA_indices_KO / g))

            # which pairwise combinations are present?
            interactions = collect(combinations(genes_indices_KO, 2))

            # Store represented combinations in matrix
            for interaction in interactions
                j = interaction[1]; k = interaction[2]
                X_interactions_count[j,k] = 1; X_interactions_count[k,j] = 1; X_interactions_count[j,j] = 1; X_interactions_count[k,k] = 1
            end
        end

        push!(Nₓ₂_vec, Nₓ₂)
        end

    # Calculate expected value and standard deviation
    E_Nₓ₂ = mean(Nₓ₂_vec)
    sd_Nₓ₂ = std(Nₓ₂_vec)

    return E_Nₓ₂, sd_Nₓ₂
end

"""
    BioCCP_Nₓ₂(x,
                g,
                r,
                n_gRNA_total,
                p_gRNA_freq,
                p_gRNA_edit, ϵ_KO)

Computes the expected value and the standard deviation of the minimal plant library size for full coverage of all pairwise combinations of gene knockouts in a multiplex CRISPR/Cas experiment (E[Nx,2] and σ[Nx,2]) using BioCCP

    ***INPUT***
x: number of target genes in the experiment
g: number of gRNAs designed per target gene
r: number of gRNA sequences per combinatorial gRNA/Cas construct
n_gRNA_total: total number of gRNAs in the experiment
p_gRNA_freq: vector with relative frequencies for all gRNAs in the construct library (normalized!)
p_gRNA_edit: vector with genome editing efficiencies for all gRNAs
ϵ_KO: global knockout efficiency; fraction of mutations leading to effective gene knockout

***OUTPUT***
E_Nₓ₂: expected value of the plant library size for full coverage of all pairwise combinations of gene knockouts
sd_Nₓ₂: standard deviation on the plant library size for full coverage of all pairwise combinations of gene knockouts
"""
function BioCCP_Nₓ₂(x::Int64,
                    g::Int64,
                    r::Int64,
                    n_gRNA_total::Int64,
                    p_gRNA_freq::Vector{Float64},
                    p_gRNA_edit::Vector{Float64},
                    ϵ_KO::Real)

    ### constraints on parameters
    @assert x > 0 && g > 0
    @assert x * g == n_gRNA_total
    @assert r >= 2
    @assert 1 - 1e-4 < sum(p_gRNA_freq) < 1 + 1e4
    @assert length(p_gRNA_freq) == n_gRNA_total
    @assert all(0 .< p_gRNA_edit .<= 1)
    @assert length(p_gRNA_edit) == n_gRNA_total
    @assert 0 < ϵ_KO <= 1

    # how many pairwise combinations of gRNAs
    ind_combinations_gRNA = collect(combinations(1:n_gRNA_total, 2))
    n_combinations_gRNA = length(ind_combinations_gRNA)

    # calculate probability and activity of gRNA combinations
    p_combinations_gRNA_library = zeros(n_combinations_gRNA)
    p_combinations_gRNA_act = zeros(n_combinations_gRNA)
    for i in 1:n_combinations_gRNA
        p_combinations_gRNA_library[i] = p_gRNA_freq[ind_combinations_gRNA[i][1]] * p_gRNA_freq[ind_combinations_gRNA[i][2]]
        p_combinations_gRNA_act[i] = p_gRNA_edit[ind_combinations_gRNA[i][1]] * p_gRNA_edit[ind_combinations_gRNA[i][2]]
    end

    # normalize probability gRNA combinations
    p_combinations_gRNA_library /= sum(p_combinations_gRNA_library)

    # select pairwise gRNA combinations of which each component codes for different gene (goal is to study combinations of knockouts in different genes)
    p_combinations_gRNA_library_interest = []
    p_combinations_gRNA_act_interest = []
    ind_combinations_gRNA_interest = []
    for i in 1:n_combinations_gRNA
        if ceil(ind_combinations_gRNA[i][1]/g) != ceil(ind_combinations_gRNA[i][2]/g)
            push!(p_combinations_gRNA_library_interest, p_combinations_gRNA_library[i])
            push!(p_combinations_gRNA_act_interest, p_combinations_gRNA_act[i])
            push!(ind_combinations_gRNA_interest, ind_combinations_gRNA[i])
        end
    end

    n_combinations_gRNA_interest = length(p_combinations_gRNA_library_interest)
    p_combinations_gRNA = p_combinations_gRNA_library_interest .* p_combinations_gRNA_act_interest * ϵ_KO^2

    # sum up probabilities or gRNA combinations for corresponding gene knockout combinations
    p_genes_matrix = zeros(x, x)
    for i in 1:n_combinations_gRNA_interest
        gene1 = Int(ceil(ind_combinations_gRNA_interest[i][1]/g))
        gene2 = Int(ceil(ind_combinations_gRNA_interest[i][2]/g))
        p_genes_matrix[gene1, gene2] += p_combinations_gRNA[i]
    end
    p_genes = collect([p_genes_matrix[i, j] for j in 2:size(p_genes_matrix, 1) for i in 1:j-1])
    n_combinations_genes = length(p_genes)
    combinations_pp = length(collect(combinations(1:r, 2)))

    # Apply BioCCP functions
    E_Nₓ₂ = expectation_minsamplesize(n_combinations_genes; p=p_genes, r=combinations_pp, normalize=false)
    sd_Nₓ₂ = std_minsamplesize(n_combinations_genes; p=p_genes, r=combinations_pp, normalize=false)

    return E_Nₓ₂, sd_Nₓ₂
end

"""
    simulate_Nₓ₂_countKOs(x,
                                         g,
                                         r,
                                         n_gRNA_total,
                                         p_gRNA_freq,
                                         p_gRNA_edit, ϵ_KO; iter=100000)

Counts the number of knockouts per plant in the experiment.

***INPUT***
x: number of target genes in the experiment
g: number of gRNAs designed per target gene
r: number of gRNA sequences per combinatorial gRNA/Cas construct
n_gRNA_total: total number of gRNAs in the experiment
p_gRNA_freq: vector with relative frequencies for all gRNAs in the construct library (normalized!)
p_gRNA_edit: vector with genome editing efficiencies for all gRNAs
ϵ_KO: global knockout efficiency; fraction of mutations leading to effective gene knockout
iter: number of plants that are sampled

***OUTPUT***
n_KOs_vec: vector with the number of knockouts for each plant
"""
function simulate_Nₓ₂_countKOs(x,
                                g,
                                r,
                                n_gRNA_total,
                                p_gRNA_freq,
                                p_gRNA_edit,
                                ϵ_KO;
                                iter = 100000)

    @assert x * g == n_gRNA_total

    n_KOs_vec = []

    for j in 1:iter
        # sample combinatorial gRNA/Cas9 construct
        gRNA_indices_construct = findall((rand(Multinomial(r, p_gRNA_freq))) .!= 0)

        # execute mutations
        gRNA_indices_mutations = [gRNA for gRNA in gRNA_indices_construct if rand(Binomial(1, p_gRNA_edit[gRNA])) == 1]

        # effective gene knockout (loss-of-function) ?
        gRNA_indices_KO = [gRNA for gRNA in gRNA_indices_mutations if rand(Binomial(1, ϵ_KO)) == 1]

        # which genes are knocked out?
        genes_indices_KO = Int.(ceil.(gRNA_indices_KO / g))

        push!(n_KOs_vec, length(unique((genes_indices_KO))))
    end

    return n_KOs_vec
end

"""
    simulate_Nₓ₃(x,
                    g,
                    r,
                    n_gRNA_total,
                    p_gRNA_freq,
                    p_gRNA_edit,
                    ϵ_KO;
                    iter=500)

Computes the expected value and the standard deviation of the minimal plant library size for full coverage of all triple combinations of gene knockouts in
a multiplex CRISPR/Cas experiment (E[Nx,3] and σ[Nx,3]) using simulation

***INPUT***
x: number of target genes in the experiment
g: number of gRNAs designed per target gene
r: number of gRNA sequences per combinatorial gRNA/Cas construct
n_gRNA_total: total number of gRNAs in the experiment
p_gRNA_freq: vector with relative frequencies for all gRNAs in the construct library (normalized!)
p_gRNA_edit: vector with genome editing efficiencies for all gRNAs
ϵ_KO: global knockout efficiency; fraction of mutations leading to effective gene knockout
iter: number of CRISPR/Cas experiments that are simulated to obtain E[Nₓ₃] and σ[Nₓ₃]

***OUTPUT***
E_Nₓ₃: expecteded value of the plant library size for full coverage of all triple combinations of gene knockouts
sd_Nₓ₃: standard deviation on the plant library size for full coverage of all triple combinations of gene knockouts
"""
function simulate_Nₓ₃(x::Int64,
                    g::Int64,
                    r::Int64,
                    n_gRNA_total::Int64,
                    p_gRNA_freq::Vector{Float64},
                    p_gRNA_edit::Vector{Float64},
                    ϵ_KO::Real;
                    iter = 500)

    ### constraints on parameters
    @assert x > 0 && g > 0
    @assert x * g == n_gRNA_total
    @assert r >= 3
    @assert 1 - 1e-4 < sum(p_gRNA_freq) < 1 + 1e4
    @assert length(p_gRNA_freq) == n_gRNA_total
    @assert all(0 .< p_gRNA_edit .<= 1)
    @assert length(p_gRNA_edit) == n_gRNA_total
    @assert 0 < ϵ_KO <= 1
    @assert iter > 1

    Nₓ₃_vec = [] # stores number of plants required for each experiment

    for i in 1:iter
        X_interactions_count = zeros(x, x, x) # Initialize matrix to count triple interactions
        Nₓ₃ = 0

        while X_interactions_count != ones(x, x, x) # check if all triple combinations are present
            Nₓ₃ += 1 # count how many plants must be sampled to fill triple interaction matrix

            # sample combinatorial gRNA/Cas9 construct
            gRNA_indices_construct = findall((rand(Multinomial(r, p_gRNA_freq))) .!= 0)

            # execute mutations
            gRNA_indices_mutations = [gRNA for gRNA in gRNA_indices_construct if rand(Binomial(1, p_gRNA_edit[gRNA])) == 1]

            # effective gene knockout (loss-of-function) ?
            gRNA_indices_KO = [gRNA for gRNA in gRNA_indices_mutations if rand(Binomial(1, ϵ_KO)) == 1]

            # which genes are knocked out?
            genes_indices_KO = Int.(ceil.(gRNA_indices_KO / g))

            # which triple combinations are present?
            interactions = collect(combinations(genes_indices_KO, 3))

            # Store represented triple combinations in 3D-matrix
            for interaction in interactions
                j = interaction[1]
                k = interaction[2]
                l = interaction[3]
                X_interactions_count[j,k,l] = 1
                X_interactions_count[k,j,l] = 1
                X_interactions_count[l,j,k] = 1
                X_interactions_count[l,k,j] = 1
                X_interactions_count[j,l,k] = 1
                X_interactions_count[k,l,j] = 1

                X_interactions_count[:,l,l] .= 1
                X_interactions_count[:,k,k] .= 1
                X_interactions_count[:,j,j] .= 1
                X_interactions_count[l,:,l] .= 1
                X_interactions_count[k,:,k] .= 1
                X_interactions_count[j,:,j] .= 1

                X_interactions_count[j,j,:] .= 1
                X_interactions_count[k,k,:] .= 1
                X_interactions_count[l,l,:] .= 1
            end
        end
        push!(Nₓ₃_vec, Nₓ₃)
    end

    # calculate expected value and standard deviation
    E_Nₓ₃ = mean(Nₓ₃_vec); sd_Nₓ₃ = std(Nₓ₃_vec)

    return E_Nₓ₃, sd_Nₓ₃
end

"""
    BioCCP_Nₓ₃(x,
            g,
            r,
            n_gRNA_total,
            p_gRNA_freq,
            p_gRNA_edit,
            ϵ_KO)

Computes the expected value and the standard deviation of the minimal plant library size
for full coverage of all triple combinations of gene knockouts in a multiplex CRISPR/Cas experiment (E[Nx,3] and σ[Nx,3]) using BioCCP

***INPUT***
x: number of target genes in the experiment
g: number of gRNAs designed per target gene
r: number of gRNA sequences per combinatorial gRNA/Cas construct
n_gRNA_total: total number of gRNAs in the experiment
p_gRNA_freq: vector with relative frequencies for all gRNAs in the construct library (normalized!)
p_gRNA_edit: vector with genome editing efficiencies for all gRNAs
ϵ_KO: global knockout efficiency; fraction of mutations leading to effective gene knockout

***OUTPUT***
E_Nₓ₃: expecteded value of the plant library size for full coverage of all triple combinations of gene knockouts
sd_Nₓ₃: standard deviation on the plant library size for full coverage of all triple combinations of gene knockouts
"""
function BioCCP_Nₓ₃(x::Int64,
                    g::Int64,
                    r::Int64,
                    n_gRNA_total::Int64,
                    p_gRNA_freq::Vector{Float64},
                    p_gRNA_edit::Vector{Float64},
                    ϵ_KO::Real)

    ### constraints on parameters
    @assert x > 0 && g > 0
    @assert x * g == n_gRNA_total
    @assert r >= 3
    @assert 1 - 1e-4 < sum(p_gRNA_freq) < 1 + 1e4
    @assert length(p_gRNA_freq) == n_gRNA_total
    @assert all(0 .< p_gRNA_edit .<= 1)
    @assert length(p_gRNA_edit) == n_gRNA_total
    @assert 0 < ϵ_KO <= 1

    # how many triple combinations of gRNAs
    ind_combinations_gRNA = collect(combinations(1:n_gRNA_total, 3))
    n_combinations_gRNA = length(ind_combinations_gRNA)

    # calculate probability and activity of triple gRNA combinations
    p_combinations_gRNA_library = zeros(n_combinations_gRNA)
    p_combinations_gRNA_act = zeros(n_combinations_gRNA)
    for i in 1:n_combinations_gRNA
        p_combinations_gRNA_library[i] = p_gRNA_freq[ind_combinations_gRNA[i][1]] * p_gRNA_freq[ind_combinations_gRNA[i][2]] * p_gRNA_freq[ind_combinations_gRNA[i][3]]
        p_combinations_gRNA_act[i] = p_gRNA_edit[ind_combinations_gRNA[i][1]] * p_gRNA_edit[ind_combinations_gRNA[i][2]] * p_gRNA_edit[ind_combinations_gRNA[i][3]]
    end

    # normalize probability gRNA combinations
    p_combinations_gRNA_library /= sum(p_combinations_gRNA_library)

    # select triple gRNA combinations of which each component codes for different gene (goal is to study combinations of knockouts in different genes)
    p_combinations_gRNA_library_interest = []
    p_combinations_gRNA_act_interest = []
    ind_combinations_gRNA_interest = []
    for i in 1:n_combinations_gRNA
        if ceil(ind_combinations_gRNA[i][1]/g) != ceil(ind_combinations_gRNA[i][2]/g) && ceil(ind_combinations_gRNA[i][1]/g) != ceil(ind_combinations_gRNA[i][3]/g) && ceil(ind_combinations_gRNA[i][3]/g) != ceil(ind_combinations_gRNA[i][2]/g)
            push!(p_combinations_gRNA_library_interest, p_combinations_gRNA_library[i])
            push!(p_combinations_gRNA_act_interest, p_combinations_gRNA_act[i])
            push!(ind_combinations_gRNA_interest, ind_combinations_gRNA[i])
        end
    end

    n_combinations_gRNA_interest = length(p_combinations_gRNA_library_interest)
    p_combinations_gRNA = p_combinations_gRNA_library_interest .* p_combinations_gRNA_act_interest * ϵ_KO^3

    # sum up probabilities or gRNA combinations for corresponding gene knockout combinations
    p_genes_matrix = zeros(x, x, x)
    for i in 1:n_combinations_gRNA_interest
        gene1 = Int(ceil(ind_combinations_gRNA_interest[i][1]/g))
        gene2 = Int(ceil(ind_combinations_gRNA_interest[i][2]/g))
        gene3 = Int(ceil(ind_combinations_gRNA_interest[i][3]/g))
        p_genes_matrix[gene1, gene2, gene3] += p_combinations_gRNA[i]
    end

    combinations_genes = collect(combinations(1:x, 3))
    p_genes = []
        for combination in combinations_genes
            push!(p_genes, p_genes_matrix[combination[1], combination[2], combination[3]])
        end

    n_combinations_genes = length(p_genes)
    combinations_pp = length(collect(combinations(1:r, 3)))

    # apply BioCCP functions
    E_Nₓ₃ = expectation_minsamplesize(n_combinations_genes; p=p_genes, r=combinations_pp, normalize=false)
    sd_Nₓ₃ = std_minsamplesize(n_combinations_genes; p=p_genes, r=combinations_pp, normalize=false)

    return E_Nₓ₃, sd_Nₓ₃
end

"""
    BioCCP_Pₓ₁(x,
                N,
                g,
                r,
                n_gRNA_total,
                p_gRNA_freq,
                p_gRNA_edit,
                ϵ_KO)

Computes the probability of full coverage of all single gene knockouts (Px,1) for an experiment with given plant library size using BioCCP

***INPUT***
x: number of target genes in the experiment
N: plant library size
g: number of gRNAs designed per target gene
r: number of gRNA sequences per combinatorial gRNA/Cas construct
n_gRNA_total: total number of gRNAs in the experiment
p_gRNA_freq: vector with relative frequencies for all gRNAs in the construct library (normalized!)
p_gRNA_edit: vector with genome editing efficiencies for all gRNAs
ϵ_KO: global knockout efficiency; fraction of mutations leading to effective gene knockout

***OUTPUT***
Pₓ₁: probability of full coverage of all single gene knockouts

"""
function BioCCP_Pₓ₁(x::Int64,
                    N::Int64,
                    g::Int64,
                    r::Int64,
                    n_gRNA_total::Int64,
                    p_gRNA_freq::Vector{Float64},
                    p_gRNA_edit::Vector{Float64},
                    ϵ_KO::Real)

    ### constraints on parameters
    @assert N > 0 && x > 0 && g > 0
    @assert x * g == n_gRNA_total
    @assert r >= 1
    @assert 1 - 1e-4 < sum(p_gRNA_freq) < 1 + 1e4
    @assert length(p_gRNA_freq) == n_gRNA_total
    @assert all(0 .< p_gRNA_edit .<= 1)
    @assert length(p_gRNA_edit) == n_gRNA_total
    @assert 0 < ϵ_KO <= 1

    # prepare input
    p_gRNAs = p_gRNA_freq .* p_gRNA_edit * ϵ_KO
    p_genes = [sum(p_gRNAs[i:i+g-1]) for i in 1:g:n_gRNA_total]

    # apply BioCCP function
    Pₓ₁ = success_probability(x, N; p=p_genes, r=r, normalize=false)

    return Pₓ₁
end

"""
BioCCP_γₓ₁(x,
            N,
            g,
            r,
            n_gRNA_total,
            p_gRNA_freq,
            p_gRNA_edit,
            ϵ_KO)

Computes the expected coverage of all single gene knockouts (Px,1) for an experiment with given plant library size using BioCCP

***INPUT***
x: number of target genes in the experiment
N: plant library size
g: number of gRNAs designed per target gene
r: number of gRNA sequences per combinatorial gRNA/Cas construct
n_gRNA_total: total number of gRNAs in the experiment
p_gRNA_freq: vector with relative frequencies for all gRNAs in the construct library (normalized!)
p_gRNA_edit: vector with genome editing efficiencies for all gRNAs
ϵ_KO: global knockout efficiency; fraction of mutations leading to effective gene knockout

***OUTPUT***
γₓ₁: expected coverage of all single gene knockouts
"""
function BioCCP_γₓ₁(x::Int64,
                    N::Int64,
                    g::Int64,
                    r::Int64,
                    n_gRNA_total::Int64,
                    p_gRNA_freq::Vector{Float64},
                    p_gRNA_edit::Vector{Float64},
                    ϵ_KO::Real)

    ### constraints on parameters
    @assert N > 0 && x > 0 && g > 0
    @assert x * g == n_gRNA_total
    @assert r >= 1
    @assert 1 - 1e-4 < sum(p_gRNA_freq) < 1 + 1e4
    @assert length(p_gRNA_freq) == n_gRNA_total
    @assert all(0 .< p_gRNA_edit .<= 1)
    @assert length(p_gRNA_edit) == n_gRNA_total
    @assert 0 < ϵ_KO <= 1

    p_gRNAs = p_gRNA_freq .* p_gRNA_edit * ϵ_KO
    p_genes = [sum(p_gRNAs[i:i+g-1]) for i in 1:g:n_gRNA_total]

    # Apply BioCCP function
    γₓ₁ = expectation_fraction_collected(x, N; p=p_genes, r=r, normalize=false)

    return γₓ₁
end

"""
    BioCCP_Pₓ₂(x,
                N,
                g,
                r,
                n_gRNA_total,
                p_gRNA_freq,
                p_gRNA_edit,
                ϵ_KO)

Computes the probability of full coverage of all pairwise combinations of gene knockouts (Px,2) for an experiment with given plant library size using BioCCP

***INPUT***
x: number of target genes in the experiment
N: plant library size
g: number of gRNAs designed per target gene
r: number of gRNA sequences per combinatorial gRNA/Cas construct
n_gRNA_total: total number of gRNAs in the experiment
p_gRNA_freq: vector with relative frequencies for all gRNAs in the construct library (normalized!)
p_gRNA_edit: vector with genome editing efficiencies for all gRNAs
ϵ_KO: global knockout efficiency; fraction of mutations leading to effective gene knockout

***OUTPUT***
Pₓ₂: probability of full coverage of all pairwise combinations of gene knockouts
"""
function BioCCP_Pₓ₂(x::Int64,
                    N::Int64,
                    g::Int64,
                    r::Int64,
                    n_gRNA_total::Int64,
                    p_gRNA_freq::Vector{Float64},
                    p_gRNA_edit::Vector{Float64},
                    ϵ_KO::Real)

    ### constraints on parameters
    @assert N > 0 && x > 0 && g > 0
    @assert x * g == n_gRNA_total
    @assert r >= 2
    @assert 1 - 1e-4 < sum(p_gRNA_freq) < 1 + 1e4
    @assert length(p_gRNA_freq) == n_gRNA_total
    @assert all(0 .< p_gRNA_edit .<= 1)
    @assert length(p_gRNA_edit) == n_gRNA_total
    @assert 0 < ϵ_KO <= 1

    # how many pairwise combinations of gRNAs
    ind_combinations_gRNA = collect(combinations(1:n_gRNA_total, 2))
    n_combinations_gRNA = length(ind_combinations_gRNA)

    # calculate probability and activity of gRNA combinations
    p_combinations_gRNA_library = zeros(n_combinations_gRNA)
    p_combinations_gRNA_act = zeros(n_combinations_gRNA)
    for i in 1:n_combinations_gRNA
        p_combinations_gRNA_library[i] = p_gRNA_freq[ind_combinations_gRNA[i][1]] * p_gRNA_freq[ind_combinations_gRNA[i][2]]
        p_combinations_gRNA_act[i] = p_gRNA_edit[ind_combinations_gRNA[i][1]] * p_gRNA_edit[ind_combinations_gRNA[i][2]]
    end

    # normalize probability gRNA combinations
    p_combinations_gRNA_library /= sum(p_combinations_gRNA_library)

    # select pairwise gRNA combinations of which each component codes for different gene (goal is to study combinations of knockouts in different genes)
    p_combinations_gRNA_library_interest = []
    p_combinations_gRNA_act_interest = []
    ind_combinations_gRNA_interest = []
    for i in 1:n_combinations_gRNA
        if ceil(ind_combinations_gRNA[i][1]/g) != ceil(ind_combinations_gRNA[i][2]/g)
            push!(p_combinations_gRNA_library_interest, p_combinations_gRNA_library[i])
            push!(p_combinations_gRNA_act_interest, p_combinations_gRNA_act[i])
            push!(ind_combinations_gRNA_interest, ind_combinations_gRNA[i])
        end
    end

    n_combinations_gRNA_interest = length(p_combinations_gRNA_library_interest)
    p_combinations_gRNA = p_combinations_gRNA_library_interest .* p_combinations_gRNA_act_interest * ϵ_KO^2

    # sum up probabilities or gRNA combinations for corresponding gene knockout combinations
    p_genes_matrix = zeros(x, x)
    for i in 1:n_combinations_gRNA_interest
        gene1 = Int(ceil(ind_combinations_gRNA_interest[i][1]/g))
        gene2 = Int(ceil(ind_combinations_gRNA_interest[i][2]/g))
        p_genes_matrix[gene1, gene2] += p_combinations_gRNA[i]
    end

    p_genes = collect([p_genes_matrix[i, j] for j in 2:size(p_genes_matrix, 1) for i in 1:j-1])
    n_combinations_genes = length(p_genes)
    combinations_pp = length(collect(combinations(1:r, 2)))

    # Apply BioCCP function
    Pₓ₂ = success_probability(n_combinations_genes, N; p=p_genes, r=combinations_pp, normalize=false)

    return Pₓ₂
end

"""
    BioCCP_γₓ₂(x,
                N,
                g,
                r,
                n_gRNA_total,
                p_gRNA_freq,
                p_gRNA_edit,
                ϵ_KO)

Computes the expected coverage of all pairwise combinations of gene knockouts (γx,2) for an experiment with given plant library size using BioCCP

***INPUT***
x: number of target genes in the experiment
N: plant library size
g: number of gRNAs designed per target gene
r: number of gRNA sequences per combinatorial gRNA/Cas construct
n_gRNA_total: total number of gRNAs in the experiment
p_gRNA_freq: vector with relative frequencies for all gRNAs in the construct library (normalized!)
p_gRNA_edit: vector with genome editing efficiencies for all gRNAs
ϵ_KO: global knockout efficiency; fraction of mutations leading to effective gene knockout

***OUTPUT***
γₓ₂: expected coverage of all pairwise combinations of gene knockouts
"""
function BioCCP_γₓ₂(x::Int64,
                    N::Int64,
                    g::Int64,
                    r::Int64,
                    n_gRNA_total::Int64,
                    p_gRNA_freq::Vector{Float64},
                    p_gRNA_edit::Vector{Float64},
                    ϵ_KO::Real)

    ### constraints on parameters
    @assert N > 0 && x > 0 && g > 0
    @assert x * g == n_gRNA_total
    @assert r >= 2
    @assert 1 - 1e-4 < sum(p_gRNA_freq) < 1 + 1e4
    @assert length(p_gRNA_freq) == n_gRNA_total
    @assert all(0 .< p_gRNA_edit .<= 1)
    @assert length(p_gRNA_edit) == n_gRNA_total
    @assert 0 < ϵ_KO <= 1

    # how many pairwise combinations of gRNAs
    ind_combinations_gRNA = collect(combinations(1:n_gRNA_total, 2))
    n_combinations_gRNA = length(ind_combinations_gRNA)

    # calculate probability and activity of gRNA combinations
    p_combinations_gRNA_library = zeros(n_combinations_gRNA)
    p_combinations_gRNA_act = zeros(n_combinations_gRNA)
    for i in 1:n_combinations_gRNA
        p_combinations_gRNA_library[i] = p_gRNA_freq[ind_combinations_gRNA[i][1]] * p_gRNA_freq[ind_combinations_gRNA[i][2]]
        p_combinations_gRNA_act[i] = p_gRNA_edit[ind_combinations_gRNA[i][1]] * p_gRNA_edit[ind_combinations_gRNA[i][2]]
    end

    # normalize probability gRNA combinations
    p_combinations_gRNA_library /= sum(p_combinations_gRNA_library)

    # select pairwise gRNA combinations of which each component codes for different gene (goal is to study combinations of knockouts in different genes)
    p_combinations_gRNA_library_interest = []
    p_combinations_gRNA_act_interest = []
    ind_combinations_gRNA_interest = []
    for i in 1:n_combinations_gRNA
        if ceil(ind_combinations_gRNA[i][1]/g) != ceil(ind_combinations_gRNA[i][2]/g)
            push!(p_combinations_gRNA_library_interest, p_combinations_gRNA_library[i])
            push!(p_combinations_gRNA_act_interest, p_combinations_gRNA_act[i])
            push!(ind_combinations_gRNA_interest, ind_combinations_gRNA[i])
        end
    end

    n_combinations_gRNA_interest = length(p_combinations_gRNA_library_interest)
    p_combinations_gRNA = p_combinations_gRNA_library_interest .* p_combinations_gRNA_act_interest * ϵ_KO^2

    # sum up probabilities or gRNA combinations for corresponding gene knockout combinations
    p_genes_matrix = zeros(x, x)
    for i in 1:n_combinations_gRNA_interest
        gene1 = Int(ceil(ind_combinations_gRNA_interest[i][1]/g))
        gene2 = Int(ceil(ind_combinations_gRNA_interest[i][2]/g))
        p_genes_matrix[gene1, gene2] += p_combinations_gRNA[i]
    end

    p_genes = collect([p_genes_matrix[i, j] for j in 2:size(p_genes_matrix, 1) for i in 1:j-1])
    n_combinations_genes = length(p_genes)
    combinations_pp = length(collect(combinations(1:r, 2)))

    # Apply BioCCP function
    γₓ₂ = expectation_fraction_collected(n_combinations_genes, N; p=p_genes, r=combinations_pp, normalize=false)

    return γₓ₂
end

"""
    BioCCP_γₓ₃(x,
                N,
                g,
                r,
                n_gRNA_total,
                p_gRNA_freq,
                p_gRNA_edit,
                ϵ_KO)

Computes the expected coverage of all triple combinations of gene knockouts (γx,3) for an experiment with given plant library size using BioCCP

***INPUT***
x: number of target genes in the experiment
N: plant library size
g: number of gRNAs designed per target gene
r: number of gRNA sequences per combinatorial gRNA/Cas construct
n_gRNA_total: total number of gRNAs in the experiment
p_gRNA_freq: vector with relative frequencies for all gRNAs in the construct library (normalized!)
p_gRNA_edit: vector with genome editing efficiencies for all gRNAs
ϵ_KO: global knockout efficiency; fraction of mutations leading to effective gene knockout

***OUTPUT***
γₓ₃: expected coverage of all triple combinations of gene knockouts
"""
function BioCCP_γₓ₃(x::Int64,
                    N::Int64,
                    g::Int64,
                    r::Int64,
                    n_gRNA_total::Int64,
                    p_gRNA_freq::Vector{Float64},
                    p_gRNA_edit::Vector{Float64},
                    ϵ_KO::Real)

    ### constraints on parameters
    @assert N > 0 && x > 0 && g > 0
    @assert x * g == n_gRNA_total
    @assert r >= 3
    @assert 1 - 1e-4 < sum(p_gRNA_freq) < 1 + 1e4
    @assert length(p_gRNA_freq) == n_gRNA_total
    @assert all(0 .< p_gRNA_edit .<= 1)
    @assert length(p_gRNA_edit) == n_gRNA_total
    @assert 0 < ϵ_KO <= 1

    # how many triple combinations of gRNAs
    ind_combinations_gRNA = collect(combinations(1:n_gRNA_total, 3))
    n_combinations_gRNA = length(ind_combinations_gRNA)

    # calculate probability and activity of triple gRNA combinations
    p_combinations_gRNA_library = zeros(n_combinations_gRNA)
    p_combinations_gRNA_act = zeros(n_combinations_gRNA)
    for i in 1:n_combinations_gRNA
        p_combinations_gRNA_library[i] = p_gRNA_freq[ind_combinations_gRNA[i][1]] * p_gRNA_freq[ind_combinations_gRNA[i][2]] * p_gRNA_freq[ind_combinations_gRNA[i][3]]
        p_combinations_gRNA_act[i] = p_gRNA_edit[ind_combinations_gRNA[i][1]] * p_gRNA_edit[ind_combinations_gRNA[i][2]] * p_gRNA_edit[ind_combinations_gRNA[i][3]]
    end

    # normalize probability gRNA combinations
    p_combinations_gRNA_library /= sum(p_combinations_gRNA_library)

    # select triple gRNA combinations of which each component codes for different gene (goal is to study combinations of knockouts in different genes)
    p_combinations_gRNA_library_interest = []
    p_combinations_gRNA_act_interest = []
    ind_combinations_gRNA_interest = []
    for i in 1:n_combinations_gRNA
        if ceil(ind_combinations_gRNA[i][1]/g) != ceil(ind_combinations_gRNA[i][2]/g) && ceil(ind_combinations_gRNA[i][1]/g) != ceil(ind_combinations_gRNA[i][3]/g) && ceil(ind_combinations_gRNA[i][3]/g) != ceil(ind_combinations_gRNA[i][2]/g)
            push!(p_combinations_gRNA_library_interest, p_combinations_gRNA_library[i])
            push!(p_combinations_gRNA_act_interest, p_combinations_gRNA_act[i])
            push!(ind_combinations_gRNA_interest, ind_combinations_gRNA[i])
        end
    end

    n_combinations_gRNA_interest = length(p_combinations_gRNA_library_interest)
    p_combinations_gRNA = p_combinations_gRNA_library_interest .* p_combinations_gRNA_act_interest * ϵ_KO^3

    # sum up probabilities or gRNA combinations for corresponding gene knockout combinations
    p_genes_matrix = zeros(x, x, x)
    for i in 1:n_combinations_gRNA_interest
        gene1 = Int(ceil(ind_combinations_gRNA_interest[i][1]/g))
        gene2 = Int(ceil(ind_combinations_gRNA_interest[i][2]/g))
        gene3 = Int(ceil(ind_combinations_gRNA_interest[i][3]/g))
        p_genes_matrix[gene1, gene2, gene3] += p_combinations_gRNA[i]
    end

    combinations_genes = collect(combinations(1:x, 3))
    p_genes = []
        for combination in combinations_genes
            push!(p_genes, p_genes_matrix[combination[1], combination[2], combination[3]])
        end

    n_combinations_genes = length(p_genes)
    combinations_pp = length(collect(combinations(1:r, 3)))

    # Apply BioCCP function
    γₓ₃ = expectation_fraction_collected(n_combinations_genes, N; p=p_genes, r=combinations_pp, normalize=false)

    return γₓ₃
end

"""
    BioCCP_Pₓ₃(x,
                N,
                g,
                r,
                n_gRNA_total,
                p_gRNA_freq,
                p_gRNA_edit,
                ϵ_KO)

Computes the probability of full coverage of all triple combinations of gene knockouts (Px,3) for an experiment with given plant library size using BioCCP

***INPUT***
x: number of target genes in the experiment
N: plant library size
g: number of gRNAs designed per target gene
r: number of gRNA sequences per combinatorial gRNA/Cas construct
n_gRNA_total: total number of gRNAs in the experiment
p_gRNA_freq: vector with relative frequencies for all gRNAs in the construct library (normalized!)
p_gRNA_edit: vector with genome editing efficiencies for all gRNAs
ϵ_KO: global knockout efficiency; fraction of mutations leading to effective gene knockout

***OUTPUT***
Pₓ₃: probability of full coverage of all triple combinations of gene knockouts
"""
function BioCCP_Pₓ₃(x::Int64,
                    N::Int64,
                    g::Int64,
                    r::Int64,
                    n_gRNA_total::Int64,
                    p_gRNA_freq::Vector{Float64},
                    p_gRNA_edit::Vector{Float64},
                    ϵ_KO::Real)

    ### constraints on parameters
    @assert N > 0 && x > 0 && g > 0
    @assert x * g == n_gRNA_total
    @assert r >= 3
    @assert 1 - 1e-4 < sum(p_gRNA_freq) < 1 + 1e4
    @assert length(p_gRNA_freq) == n_gRNA_total
    @assert all(0 .< p_gRNA_edit .<= 1)
    @assert length(p_gRNA_edit) == n_gRNA_total
    @assert 0 < ϵ_KO <= 1

    # how many triple combinations of gRNAs
    ind_combinations_gRNA = collect(combinations(1:n_gRNA_total, 3))
    n_combinations_gRNA = length(ind_combinations_gRNA)

    # calculate probability and activity of triple gRNA combinations
    p_combinations_gRNA_library = zeros(n_combinations_gRNA)
    p_combinations_gRNA_act = zeros(n_combinations_gRNA)
    for i in 1:n_combinations_gRNA
        p_combinations_gRNA_library[i] = p_gRNA_freq[ind_combinations_gRNA[i][1]] * p_gRNA_freq[ind_combinations_gRNA[i][2]] * p_gRNA_freq[ind_combinations_gRNA[i][3]]
        p_combinations_gRNA_act[i] = p_gRNA_edit[ind_combinations_gRNA[i][1]] * p_gRNA_edit[ind_combinations_gRNA[i][2]] * p_gRNA_edit[ind_combinations_gRNA[i][3]]
    end

    # normalize probability gRNA combinations
    p_combinations_gRNA_library /= sum(p_combinations_gRNA_library)

    # select triple gRNA combinations of which each component codes for different gene (goal is to study combinations of knockouts in different genes)
    p_combinations_gRNA_library_interest = []
    p_combinations_gRNA_act_interest = []
    ind_combinations_gRNA_interest = []
    for i in 1:n_combinations_gRNA
        if ceil(ind_combinations_gRNA[i][1]/g) != ceil(ind_combinations_gRNA[i][2]/g) && ceil(ind_combinations_gRNA[i][1]/g) != ceil(ind_combinations_gRNA[i][3]/g) && ceil(ind_combinations_gRNA[i][3]/g) != ceil(ind_combinations_gRNA[i][2]/g)
            push!(p_combinations_gRNA_library_interest, p_combinations_gRNA_library[i])
            push!(p_combinations_gRNA_act_interest, p_combinations_gRNA_act[i])
            push!(ind_combinations_gRNA_interest, ind_combinations_gRNA[i])
        end
    end

    n_combinations_gRNA_interest = length(p_combinations_gRNA_library_interest)
    p_combinations_gRNA = p_combinations_gRNA_library_interest .* p_combinations_gRNA_act_interest * ϵ_KO^3

    # sum up probabilities or gRNA combinations for corresponding gene knockout combinations
    p_genes_matrix = zeros(x, x, x)
    for i in 1:n_combinations_gRNA_interest
        gene1 = Int(ceil(ind_combinations_gRNA_interest[i][1]/g))
        gene2 = Int(ceil(ind_combinations_gRNA_interest[i][2]/g))
        gene3 = Int(ceil(ind_combinations_gRNA_interest[i][3]/g))
        p_genes_matrix[gene1, gene2, gene3] += p_combinations_gRNA[i]
    end

    combinations_genes = collect(combinations(1:x, 3))
    p_genes = []
        for combination in combinations_genes
            push!(p_genes, p_genes_matrix[combination[1], combination[2], combination[3]])
        end

    n_combinations_genes = length(p_genes)
    combinations_pp = length(collect(combinations(1:r, 3)))

    # Apply BioCCP function
    Pₓ₃ = success_probability(n_combinations_genes, N; p=p_genes, r=combinations_pp, normalize=false)

    return Pₓ₃
end

max_stars_repo_name: kirstvh/MultiplexCrisprDOE

max_issues_repo_name: kirstvh/MultiplexCrisprDOE

max_forks_repo_name: kirstvh/MultiplexCrisprDOE


text:
using MixedBlockArrays
using Test

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


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

using DataFrames
using Query
using FGenerators

using Lazy: @> # or using DataConvenience: @>, but Lazy exports groupby
using DataFramesMeta: @where
using Statistics: mean
using Unicode:normalize

const _AbstractColumnName = Union{AbstractString, Signed, Symbol, Unsigned}

### R-like function names
# Returns true for true, false, 1, and 0.  Otherwise returns false
_isbool(x) = try Bool(x) catch e false end
# Returns true for a given column in a dataframe if that column is categorical data
_iscategorical(df::DataFrame, col::Symbol) = isa(df[col], CategoricalVector)
_iscategorical(V::AbstractVector) = isa(V, CategoricalVector)
# Similar to as.factor
_ascategorical!(df::DataFrame, col::Symbol) = categorical!(df, col)
DataFrames.categorical!(V::AbstractVector) = CategoricalVector(eltype(V)[v for v in V])
_ascategorical!(V::AbstractVector) = categorical!(V)
# Returns true if the input is boolean
_islogical(x::Union{Bool, Missing}) = true
function _islogical(s::AbstractString)
    s′ = normalize(s, casefold = true)
    if s′ ∈ ("true", "false", "na", "n/a", "missing") # sometimes better than `uppercase` or `lowercase` as it handles unicode characters
        return true
    else
        return false
    end
end
_islogical(x) = false
_islogical(X::AbstractArray) = all(_islogical.(X))
# Returns false for 0 and false.  Otherwise returns true
_aslogical(x::Bool) = Bool(x)
_aslogical(x::Number) = iszero(x) ? false : true
_aslogical(X::AbstractArray) = _aslogical.(X)
function _aslogical(s::AbstractString)
    s′ = normalize(s, casefold = true) # sometimes better than `uppercase` or `lowercase` as it handles unicode characters
    if s′ == "true"
        return true
    elseif s′ == "false"
        return false
    elseif s′ ∈ ("na", "n/a", "missing")
        return missing
    else
        # return nothing
        return missing
    end
end
_aslogical!(df::DataFrame, col::Symbol) = (df[col] = _aslogical(df[col]))
# Returns true for a numerical value
_isnumeric(x::Number) = true
_isnumeric(x) = false
_isnumeric(X::AbstractArray) = all(_isnumeric.(X))
_isnumeric(df::DataFrame, col::Symbol) = _isnumeric(df[col])
# Attempts to parse input as a numeric value
_asnumeric(x::Number) =
_asnumeric(x::Bool) = x ? 1 : 0
_asnumeric(X::AbstractArray) = _asnumeric.(X)
_asnumeric!(df::DataFrame, col::Symbol) = (df[col] = _asnumeric(df[col]))
# equivalent behaviour to `is.na`
Base.ismissing(x::Nothing) = nothing
Base.ismissing(X::AbstractArray) = ismissing.(X)
# equivalent behaviour to `is.null`
Base.isnothing(X::AbstractArray) = all(isnothing.(X))
### Functions defined in eyetrackingR
function _islogical2(x)
    if _islogical(x)
        return true
    elseif _isnumeric(x)
        return false
    else
        throw(error("""
            One of your columns could not be converted
            to the correct format (Bool).  Please do so manually.
        """))
    end
end

# a function from eyetrackingR, but do we need it?
# _asnumeric2(x::T) where T = _asnumeric(parse(T, string(x)))
# _asnumeric2(X::AbstractArray) = _asnumeric2.(X)

# `_check_then_convert` function which takes in a dataframe.
# The `check_function` and `convert_function` necessarily
# uses the dataframe as a parameter
function _check_then_convert(
    df::DataFrame,
    x,
    check_function::Function,
    convert_function::Function,
    colname::Union{String, Symbol};
    verbose::Bool = true
    )
    x′ = Symbol(x)
    if ! check_function(df, x′)
        verbose && println("Converting $(colname) to proper type...")
        x = convert_function(df, x′)
        verbose && println("...done")
    end

    if colname == "Trackloss" && any(ismissing(x))
        @warn "Found NAs in trackloss column.  These will be treated as TRACKLOSS = false."
        x = ismissing(x) ? false : x
    end

    return x
end

function _check_then_convert(
    x,
    check_function::Function,
    convert_function::Function,
    colname::Union{String, Symbol};
    verbose::Bool = true
    )

    if ! check_function(x)
        verbose && println("Converting $(colname) to proper type...")
        x = convert_function(x)
        verbose && println("...done")
    end

    if colname == "Trackloss" && any(ismissing(x))
        @warn "Found NAs in trackloss column.  These will be treated as TRACKLOSS = false."
        x = ismissing(x) ? false : x
    end

    return x
end

### Helper functions
@fgenerator function deepvalues(itr::Union{AbstractArray, AbstractDict, Tuple, NamedTuple})
    for v ∈ values(itr)
        @yieldfrom deepvalues(v)
    end
end
@fgenerator deepvalues(x) = @yield x
# Try to return named tuples or pairs!
@fgenerator function deepkeyvalues(d::AbstractDict)
    for (k, v) ∈ d
        if v isa AbstractVector
            for vi ∈ v
                @yield k => vi
            end
        else
            @yield k => v
        end
    end
end
# Returns each index of a list, but not for single values
_eachlistindex(X::AbstractArray) = eachindex(X)
_eachlistindex(x::Nothing) = [] # null values are not iterable in R
_eachlistindex(x) = 1
# Retuns an index of a list but not for julia-indexable values
_getlistindex(X::AbstractVector, n::Int) = getindex(X, n)
_getlistindex(x, n::Int) = isone(n) ? x : missing # or should we return missing?

### Main data cleaning function

## Note to consider: *_column as kwargs or normal args?
# on one hand, kwargs may imply that it would work without them
# which it would not.  On the other hand, kwargs make the
# function more readable, both in this code and in practice.
function make_clean_data(
    data::DataFrame;
    participant_column = nothing,
    trackloss_column = nothing,
    time_column = nothing,
    trial_column = nothing,
    aoi_columns = nothing,
    treat_non_aoi_looks_as_missing = nothing,
    item_columns = nothing
    )

    ## Data options:
    data_options = Dict{Symbol, Any}(
        :participant_column => participant_column,
        :trackloss_column => trackloss_column,
        :time_column => time_column,
        :trial_column => trial_column,
        :item_columns => item_columns,
        :aoi_columns => aoi_columns,
        :treat_non_looks_as_missing => treat_non_aoi_looks_as_missing
        )

    ## Check for reserved column name:
    if data_options[:time_column] == "Time"
        throw(error("""
            We apologise for the invonvenience, but your `time_column` cannot be called 'Time'.
            This name is a reserved name that EyeTrackingUtils uses.  Please rename.
        """))
    end

    if "Time" ∈ names(data)
        @warn """
        Your dataset has a column called 'time', but this column
        is reserved for EyeTrackingUtils.  Will rename to 'TimeOriginal'...
        """
        rename!(data, Dict(:Time => "TimeOriginal"))
    end

    ## Verify columns:
    out = copy(data)
    col_type_converter = Dict{Symbol, Function}(
        :participant_column => x -> _check_then_convert(x, _iscategorical, categorical!, "Participants"), # previously a data column
        :time_column => x -> _check_then_convert(x, _isnumeric, _asnumeric!, "Time"), # previously a data column
        :trial_column => x -> _check_then_convert(x, _iscategorical, categorical!, "Trial"), # previously a data column
        :trackloss_column => x -> _check_then_convert(x, _islogical2, _aslogical, "Trackloss"),
        :item_columns => x -> _check_then_convert(x, _iscategorical, categorical!, "Item"), # previously a data column
        :aoi_columns => x -> _check_then_convert(x, _islogical2, _aslogical, "AOI")
    )

    # for col in keys(col_type_converter)
    #     for i in _seq_along(data_options[col])
    #         if data_options[col] isa AbstractArray
    #             if isnothing(out[data_options[col][i]])
    #                 throw(error("Data are missing: $(col)"))
    #             end
    #             out[data_options[col][i]] = col_type_converter[col](out[data_options[col][i]])
    #         else
    #             if isnothing(out[data_options[col]])
    #                 throw(error("Data are missing: $(col)"))
    #             end
    #             out[data_options[col]] = col_type_converter[col](out[data_options[col]])
    #         end
    #     end
    # end
    for col in keys(col_type_converter)
        for i in _eachlistindex(data_options[col])
            if isnothing(out[_getlistindex(data_options[col], i)])
                throw(error("Data are missing at column $(_getlistindex(data_options[col], i))"))
            end
            println("Checking $col at $(data_options[col]) at $(_getlistindex(data_options[col], i))")
            out[_getlistindex(data_options[col], i)] = col_type_converter[col](out[_getlistindex(data_options[col], i)])
        end
    end

    ## Deal with non-AOI looks:
    if treat_non_aoi_looks_as_missing
        # any_aoi = sum(, dims=2) > 0 ? # second dimension is the sum of the rows
        any_aoi = sum(skipmissing.(out[data_options[aoi_columns]]), dims = 1) .> 0
        out[data_options[Symbol(trackloss_column)]][!any_aoi] .= true
    end

    ## Set All AOI rows with trackloss to NA:
    # this ensures that any calculations of proportion-looking will not include trackloss in the denominator
    for aoi in data_options[aoi_columns]
        out.aoi[out.data_options[trackloss_column]] = missing # or nothing?
    end

    # Check for duplicate values of Trial column within Participants
    duplicates = @> out begin
        groupby([:participant_column, :trial_column, :time_column])
        combine(nrow => :n)
        @where :n .> 1
    end

    if nrow(duplicates) > 0
        println(duplicates)
        throw(error("""
        It appears that `trial_column` is not unique within participants. See above for a summary
        of which participant*trials have duplicate timestamps. EyeTrackingUtils requires that each participant
        only have a single trial with the same `trial_column` value. If you repeated items in your experiment,
        use `item_column` to specify the name of the item, and set `trial_column` to a unique value
        (e.g., the trial index).
        """))
    end

    out = @> out begin
        groupby([:participant_column, :trial_column, :time_column])
    end

end

struct eR_data_eR_df
    out::DataFrame
    eR
end

## Assign attribute:
# class(out) <- c("eyetrackingR_data", "eyetrackingR_df", class(out))
# attr(out, "eyetrackingR") <- list(data_options = data_options)
# return(out)


text:
using MultivariateStats
using LinearAlgebra
using Test
import Statistics: mean, cov
import Random

@testset "Whitening" begin

    Random.seed!(34568)

    ## data

    d = 3
    n = 5

    X = rand(d, n)
    mval = vec(mean(X, dims=2))
    C = cov(X, dims=2)
    C0 = copy(C)

    emax = maximum(eigvals(C))
    rc = 0.01
    Cr = C + (emax * rc) * Matrix(I, d, d)

    ## cov_whitening

    cf = cholesky(Hermitian(C, :U))
    W = cov_whitening(cf)
    @test istriu(W)
    @test W'C * W ≈ Matrix(I, d, d)

    cf = cholesky(Hermitian(C, :L))
    W = cov_whitening(cf)
    @test istriu(W)
    @test W'C * W ≈ Matrix(I, d, d)

    W = cov_whitening(C)
    @test C == C0
    @test istriu(W)
    @test W'C * W ≈ Matrix(I, d, d)

    W = cov_whitening(C, 0.0)
    @test C == C0
    @test istriu(W)
    @test W'C * W ≈ Matrix(I, d, d)

    W = cov_whitening(C, rc)
    @test istriu(W)
    @test W'Cr * W ≈ Matrix(I, d, d)

    # Whitening from data

    f = fit(Whitening, X)
    W = f.W
    @test isa(f, Whitening{Float64})
    @test mean(f) === f.mean
    @test istriu(W)
    @test W'C * W ≈ Matrix(I, d, d)
    @test transform(f, X) ≈ W' * (X .- f.mean)

    f = fit(Whitening, X; regcoef=rc)
    W = f.W
    @test W'Cr * W ≈ Matrix(I, d, d)

    f = fit(Whitening, X; mean=mval)
    W = f.W
    @test W'C * W ≈ Matrix(I, d, d)

    f = fit(Whitening, X; mean=0)
    Cx = (X * X') / (n - 1)
    W = f.W
    @test W'Cx * W ≈ Matrix(I, d, d)

    # invsqrtm

    R = invsqrtm(C)
    @test C == C0
    @test R ≈ inv(sqrt(C))

    # mixing types
    X = rand(Float64, 5, 10)
    XX = convert.(Float32, X)

    M = fit(Whitening, X)
    MM = fit(Whitening, XX)

    # mixing types should not error
    transform(M, XX)
    transform(MM, X)

    # type consistency
    @test eltype(mean(M)) == Float64
    @test eltype(mean(MM)) == Float32
end


text:
######################################################################
# waypoint.jl
######################################################################

# TODO: figure out export strategy
export WaypointPolicy

# to simplify some notation in this file
#typealias VVF Vector{Vector{Float64}}
const VVF = Vector{Vector{Float64}}
export VVF

type WaypointPolicy <: Policy
	current_index::Int
	waypoints::VVF

	function WaypointPolicy(waypoints::VVF)
		new(2, waypoints)
	end
end

# TODO: actually require the third index
function action(m::SearchDomain, x::Vehicle, o::Float64, f::AbstractFilter, p::WaypointPolicy)
	dx = p.waypoints[p.current_index][1] - x.x
	dy = p.waypoints[p.current_index][2] - x.y
	p.current_index += 1

	# Just wrap around? Or maybe we should stay still?
	if p.current_index > length(p.waypoints)
		p.current_index = 1
	end

	return dx, dy, 0.0
end


text:
using MimiPAGE2009
using Test

m = MimiPAGE2009.get_model()
run(m)

# Test the values from this new components compared to the total cost calculation originally suggested by Fran

atol = 1e-8

@test all(isapprox.(m[:TotalCosts, :total_costs_percap_peryear],
    m[:EquityWeighting, :cons_percap_aftercosts] - m[:EquityWeighting, :rcons_percap_dis] + m[:EquityWeighting, :act_percap_adaptationcosts] + m[:EquityWeighting, :tct_percap_totalcosts_total],
    atol = atol))

@test all(isapprox.(m[:TotalCosts, :total_damages_percap_peryear],
    m[:EquityWeighting, :cons_percap_aftercosts] - m[:EquityWeighting, :rcons_percap_dis] + m[:EquityWeighting, :act_percap_adaptationcosts], # without abatement costs, just damages
    atol = atol))


text:
__precompile__()

module BinaryTwoStageDesigns

import JuMP
import QuadGK, Roots
import UnicodePlots
import Base.show, Base.print
import Base.convert, Base.writecsv
import Distributions, Distributions.minimum, Distributions.maximum,
    Distributions.pdf, Distributions.cdf, Distributions.rand,
    Distributions.quantile, Distributions.mean, Distributions.var
import DataFrames, DataFrames.DataFrame
import DataArrays
import MathProgBase
import CSV

# Parameters.jl
export
  Parameters,
    samplespace, maxsamplesize, isgroupsequential, allowsstoppingforefficacy, label, mcrv,
  NoParameters,
  PointAlternative,
    null, mtoer, alternative, mtter, mcrv,
  VagueAlternative,
    prior

# Design.jl
export
  Design,
    parameters, interimsamplesize, samplesize, criticalvalue, power,
    expectedpower, test, pdf, simulate, support, ispossible,
    stoppingforfutility, stoppingforefficacy, score, expectedpower, jeffreysprior,
    save, writepropertiescsv

# SampleSize.jl
export
  SampleSize

# SampleSpace.jl
export
  SampleSpace,
    interimsamplesizerange, maxsamplesize, possible

# IPModel.jl
export
  IPModel,
    extractsolution

# MESS.jl
export
  MESS,
    minconditionalpower, isgroupsequential

# LiuScore.jl
export
  LiuScore,
    ros, rup

export
  EB,
    expectedtransformedpower, expectedcost, expectedbenefit

export
  MBESS, MBESS2

# optimal.jl
export optimaldesign

# estimate.jl
export
  Estimator,
    design, estimate, pvalue, bias, rmse, mae, incompatibleoutcomes,
  MLEstimator, # include("Estimators/mle.jl")
  RBEstimator, # inlcude("Estimators/rbe.jl")
  OCEstimator

#
export
  ConfidenceInterval,
    limits, confidence, design, coverage, estimator, findinconsistencies,
    meanwidth, meaninterval,
  ECPInterval,
  CPInterval,
  MMWInterval

export adaptstageone, adaptstagetwo

export adapt_stage_one, adapt_stage_two

include("util.jl")

include("Parameters.jl")

include("Design.jl")

include("SampleSize.jl")

include("SampleSpace.jl")

include("IPModel.jl")

include("Parameters/MESS.jl")
include("Parameters/LiuScore.jl")
include("Parameters/MBESS.jl")
include("Parameters/MBESS2.jl")
include("Parameters/EB.jl")

include("optimaldesign.jl")

include("Estimator.jl")
include("Estimators/MLEstimator.jl")
include("Estimators/RBEstimator.jl")
include("Estimators/OCEstimator.jl")

include("ConfidenceInterval.jl")
include("ConfidenceIntervals/ECPInterval.jl")
include("ConfidenceIntervals/CPInterval.jl")
include("ConfidenceIntervals/MMWInterval.jl")

include("adapt.jl")

end # module

max_stars_repo_name: kkmann/binaryTwoStageDesigns

max_issues_repo_name: kkmann/binaryTwoStageDesigns

max_forks_repo_name: kkmann/binaryTwoStageDesigns


text:
module TestCopy

include("../preamble.jl")
import Tables
using Compat: Iterators
using DataFrames: DataFrame, eachrow
using StructArrays: StructVector, StructArray
using TypedTables: Table

@testset "$copy" for copy in [copy, tcopy, dcopy]
    @testset "$copy(_, ::$(prettytypeof(src)))" for src in Any[
        DataFrame(a=[1:4;], b=[5:8;]),
        StructVector(a=[1:4;], b=[5:8;]),
        Table(a=[1:4;], b=[5:8;]),
    ]
        if copy in (tcopy, dcopy) && src isa DataFrame
            @test_broken copy(Map(identity), src) ==ₜ src
            @test_broken copy(Map(identity), src; basesize=1) ==ₜ src
            continue
        end
        @test copy(Map(identity), src) ==ₜ src
        @test copy(eduction(Iterators.map(identity, src))) ==ₜ src
        if copy in (tcopy, dcopy)
            @test copy(Map(identity), src; basesize=1) ==ₜ src
            @test copy(Iterators.map(identity, src)) ==ₜ src
            @test copy(Iterators.map(identity, src); basesize = 1) ==ₜ src
        end
    end
    @testset "$copy(_, $T, ::$(prettytypeof(src)))" for src in Any[
            StructVector(a = [1:4;], b = [5:8;]),
            Table(a = [1:4;], b = [5:8;]),
        ],
        T in Any[DataFrame, StructArray, Table]

        @test copy(Map(identity), T, src) ==ₜ T(src)
        @test copy(T, eduction(Iterators.map(identity, src))) ==ₜ T(src)
        if copy in (tcopy, dcopy)
            @test copy(Map(identity), T, src; basesize = 1) ==ₜ T(src)
            @test copy(T, Iterators.map(identity, src)) ==ₜ T(src)
            @test copy(T, Iterators.map(identity, src); basesize = 1) ==ₜ T(src)
        end
    end
    @testset "$copy(_, eachrow(df))" begin
        df = DataFrame(a=[1:4;], b=[5:8;])
        if Tables.materializer(eachrow(df)) == DataFrame
            # requires https://github.com/JuliaData/DataFrames.jl/pull/2055
            if copy in (tcopy, dcopy)
                @test copy(Map(identity), eachrow(df); basesize=1) ==ₜ df
                @test copy(Iterators.map(identity, eachrow(df)); basesize = 1) ==ₜ df
                @test copy(eduction(Iterators.map(identity, eachrow(df))); basesize = 1) ==ₜ
                      df
            else
                @test copy(Map(identity), eachrow(df)) ==ₜ df
                @test copy(eduction(Iterators.map(identity, eachrow(df)))) ==ₜ df
            end
        else
            @test_broken copy(Map(identity), eachrow(df)) ==ₜ df
        end
    end
end

@testset "$collect" for collect in [collect, tcollect, dcollect]
    @testset "$collect(_, ::$(prettytypeof(src)))" for src in Any[
        1:64,
        StructVector(a=[1:4;], b=[5:8;]),
     ]
        @test collect(Map(identity), src) ==ₜ Base.collect(src)
        @test collect(IdentityTransducer(), src) ==ₜ Base.collect(src)
        @test collect(eduction(Iterators.map(identity, src))) ==ₜ Base.collect(src)
        if collect in (tcollect, dcollect)
            @test collect(Map(identity), src; basesize=1) ==ₜ Base.collect(src)
            @test collect(Iterators.map(identity, src)) ==ₜ Base.collect(src)
            @test collect(Iterators.map(identity, src); basesize = 1) ==ₜ Base.collect(src)
        end
    end
end

end  # module