text:
# Borrowed from the ModernGL.jl example created by Simon Danisch
# Apparently, this might originally come from @yurivish (https://github.com/yurivish)
# The copy here has only been very slightly modified to add convience functions
# for generating frame and render buffers.

function glGetInfoLog(obj::GLuint)
	isShader = glIsShader(obj)
	getiv = isShader == GL_TRUE ? glGetShaderiv : glGetProgramiv
	getInfo = isShader == GL_TRUE ? glGetShaderInfoLog : glGetProgramInfoLog

	int = GLint[0]
	getiv(obj, GL_INFO_LOG_LENGTH, int)
	maxlength = int[1]

	if maxlength > 0
		buffer = zeros(GLchar, maxlength)
		sizei = GLsizei[0]
		getInfo(obj, maxlength, sizei, buffer)
		length = sizei[1]
		bytestring(pointer(buffer), length)
	else
		""
	end
end

function validateShader(shader)
	success = GLint[0]
	glGetShaderiv(shader, GL_COMPILE_STATUS, success)
	success[1] == GL_TRUE
end

function glErrorMessage()
	err = glGetError()
	err == GL_NO_ERROR ? "" :
	err == GL_INVALID_ENUM ? "GL_INVALID_ENUM: An unacceptable value is specified for an enumerated argument. The offending command is ignored and has no other side effect than to set the error flag." :
	err == GL_INVALID_VALUE ? "GL_INVALID_VALUE: A numeric argument is out of range. The offending command is ignored and has no other side effect than to set the error flag." :
	err == GL_INVALID_OPERATION ? "GL_INVALID_OPERATION: The specified operation is not allowed in the current state. The offending command is ignored and has no other side effect than to set the error flag." :
	err == GL_INVALID_FRAMEBUFFER_OPERATION ? "GL_INVALID_FRAMEBUFFER_OPERATION: The framebuffer object is not complete. The offending command is ignored and has no other side effect than to set the error flag." :
	err == GL_OUT_OF_MEMORY ? "GL_OUT_OF_MEMORY: There is not enough memory left to execute the command. The state of the GL is undefined, except for the state of the error flags, after this error is recorded." : "Unknown OpenGL error with error code $err."
end

function glCheckError(actionName="")
	message = glErrorMessage()
	if length(message) > 0
		if length(actionName) > 0
			error("Error ", actionName, ": ", message)
		else
			error("Error: ", message)
		end
	end
end

function glGenOne(glGenFn)
	id = GLuint[0]
	glGenFn(1, id)
	glCheckError("generating a buffer, array, or texture")
	id[1]
end

glGenBuffer() = glGenOne(glGenBuffers)
glGenFramebuffer() = glGenOne(glGenFramebuffers)
glGenRenderbuffer() = glGenOne(glGenRenderbuffers)
glGenVertexArray() = glGenOne(glGenVertexArrays)
glGenTexture() = glGenOne(glGenTextures)

global GLSL_VERSION = ""

function createcontextinfo()
	global GLSL_VERSION

	glsl = split(bytestring(glGetString(GL_SHADING_LANGUAGE_VERSION)), ['.', ' '])
	if length(glsl) >= 2
		glsl = VersionNumber(parse(Int, glsl[1]), parse(Int, glsl[2]))
		GLSL_VERSION = string(glsl.major) * rpad(string(glsl.minor),2,"0")
	else
		error("Unexpected version number string. Please report this bug! GLSL version string: $(glsl)")
	end

	glv = split(bytestring(glGetString(GL_VERSION)), ['.', ' '])
	if length(glv) >= 2
		glv = VersionNumber(parse(Int, glv[1]), parse(Int, glv[2]))
	else
		error("Unexpected version number string. Please report this bug! OpenGL version string: $(glv)")
	end

	dict = Dict{Symbol, Any}()
	dict[:glsl_version]	= glsl
	dict[:gl_version] = glv
	dict[:gl_vendor] = bytestring(glGetString(GL_VENDOR))
	dict[:gl_renderer]	= bytestring(glGetString(GL_RENDERER))
	dict
end

function get_glsl_version_string()
	if isempty(GLSL_VERSION)
		error("couldn't get GLSL version, GLUTils not initialized, or context not created?")
	end

	return "#version $(GLSL_VERSION)\n"
end

function createShader(source, typ)
	shader = glCreateShader(typ)::GLuint
	if shader == 0
		error("Error creating shader: ", glErrorMessage())
	end
	glShaderSource(shader, 1, convert(Ptr{Uint8}, pointer([convert(Ptr{GLchar}, pointer(source))])), C_NULL)
	glCompileShader(shader)
	!validateShader(shader) && error("Shader creation error: ", glGetInfoLog(shader))
	return shader
end

function createShaderProgram(f, vertexShader, fragmentShader)
	prog = glCreateProgram()
	if prog == 0
		error("Error creating shader program: ", glErrorMessage())
	end

	glAttachShader(prog, vertexShader)
	glCheckError("attaching vertex shader")

	glAttachShader(prog, fragmentShader)
	glCheckError("attaching fragment shader")
	f(prog)

	glLinkProgram(prog)
	status = GLint[0]

	glGetProgramiv(prog, GL_LINK_STATUS, status)

	if status[1] == GL_FALSE
		glDeleteProgram(prog)
		error("Error linking shader: ", glGetInfoLog(prog))
	end
	prog
end

createShaderProgram(vertexShader, fragmentShader) = createShaderProgram(prog->0, vertexShader, fragmentShader)

function removeAllShaders(shaderBank)
	for (n, p) in shaderBank
		glUseProgram(p)
		as = glGetAttachedShaders(p)
		for s in as
			glDetachShader(p, s)
			glDeleteShader(s)
		end
		glUseProgram(0)
		glDeleteProgram(p)
	end
end


text:
"""
    eachobs(data; kws...)

Return an iterator over `data`.

Supports the same arguments as [`DataLoader`](@ref).
The `batchsize` default is `-1` here while
it is `1` for `DataLoader`.

# Examples

```julia
X = rand(4,100)

for x in eachobs(X)
    # loop entered 100 times
    @assert typeof(x) <: Vector{Float64}
    @assert size(x) == (4,)
end

# mini-batch iterations
for x in eachobs(X, batchsize=10)
    # loop entered 10 times
    @assert typeof(x) <: Matrix{Float64}
    @assert size(x) == (4,10)
end

# support for tuples, named tuples, dicts
for (x, y) in eachobs((X, Y))
    # ...
end
```
"""
function eachobs(data; batchsize=-1, kws...)
    DataLoader(data; batchsize, kws...)
end

"""
    DataLoader(data; [batchsize, buffer, partial, shuffle, parallel, rng])

An object that iterates over mini-batches of `data`,
each mini-batch containing `batchsize` observations
(except possibly the last one).

Takes as input a single data array, a tuple (or a named tuple) of arrays,
or in general any `data` object that implements the [`numobs`](@ref) and [`getobs`](@ref)
methods.

The last dimension in each array is the observation dimension, i.e. the one
divided into mini-batches.

The original data is preserved in the `data` field of the DataLoader.

# Arguments

- `data`: The data to be iterated over. The data type has to be supported by
  [`numobs`](@ref) and [`getobs`](@ref).
- `buffer`: If `buffer=true` and supported by the type of `data`,
a buffer will be allocated and reused for memory efficiency.
You can also pass a preallocated object to `buffer`. Default `false`.
- `batchsize`: If less than 0, iterates over individual observations.
Otherwise, each iteration (except possibly the last) yields a mini-batch
containing `batchsize` observations. Default `1`.
- `partial`: This argument is used only when `batchsize > 0`.
  If `partial=false` and the number of observations is not divisible by the batchsize,
  then the last mini-batch is dropped. Default `true`.
- `parallel`: Whether to use load data in parallel using worker threads. Greatly
    speeds up data loading by factor of available threads. Requires starting
    Julia with multiple threads. Check `Threads.nthreads()` to see the number of
    available threads. **Passing `parallel = true` breaks ordering guarantees**.
    Default `false`.
- `shuffle`: Whether to shuffle the observations before iterating. Unlike
    wrapping the data container with `shuffleobs(data)`, `shuffle=true` ensures
    that the observations are shuffled anew every time you start iterating over
    `eachobs`. Default `false`.
- `collate`: Batching behavior. If `nothing` (default), a batch is `getobs(data, indices)`. If `false`, each batch is
    `[getobs(data, i) for i in indices]`. When `true`, applies [`batch`](@ref) to the vector of observations in a batch,
   recursively collating arrays in the last dimensions. See [`batch`](@ref) for more information and examples.
- `rng`: A random number generator. Default `Random.GLOBAL_RNG`

# Examples

```jldoctest
julia> Xtrain = rand(10, 100);

julia> array_loader = DataLoader(Xtrain, batchsize=2);

julia> for x in array_loader
         @assert size(x) == (10, 2)
         # do something with x, 50 times
       end

julia> array_loader.data === Xtrain
true

julia> tuple_loader = DataLoader((Xtrain,), batchsize=2);  # similar, but yielding 1-element tuples

julia> for x in tuple_loader
         @assert x isa Tuple{Matrix}
         @assert size(x[1]) == (10, 2)
       end

julia> Ytrain = rand('a':'z', 100);  # now make a DataLoader yielding 2-element named tuples

julia> train_loader = DataLoader((data=Xtrain, label=Ytrain), batchsize=5, shuffle=true);

julia> for epoch in 1:100
         for (x, y) in train_loader  # access via tuple destructuring
           @assert size(x) == (10, 5)
           @assert size(y) == (5,)
           # loss += f(x, y) # etc, runs 100 * 20 times
         end
       end

julia> first(train_loader).label isa Vector{Char}  # access via property name
true

julia> first(train_loader).label == Ytrain[1:5]  # because of shuffle=true
false

julia> foreach(println∘summary, DataLoader(rand(Int8, 10, 64), batchsize=30))  # partial=false would omit last
10×30 Matrix{Int8}
10×30 Matrix{Int8}
10×4 Matrix{Int8}
```
"""
struct DataLoader{T, R<:AbstractRNG, C<:Val}
    data::T
    batchsize::Int
    buffer::Bool
    partial::Bool
    shuffle::Bool
    parallel::Bool
    collate::C
    rng::R
end

function DataLoader(
        data;
        buffer = false,
        parallel = false,
        shuffle = false,
        batchsize::Int = 1,
        partial::Bool = true,
        collate = Val(nothing),
        rng::AbstractRNG = Random.GLOBAL_RNG)
    buffer = buffer isa Bool ? buffer : true
    collate = collate isa Val ? collate : Val(collate)
    if !(collate ∈ (Val(nothing), Val(true), Val(false)))
        throw(ArgumentError("`collate` must be one of `nothing`, `true` or `false`."))
    end
    return DataLoader(data, batchsize, buffer, partial, shuffle, parallel, collate, rng)
end

function Base.iterate(e::DataLoader)
    # Wrapping with ObsView in order to work around
    # issue https://github.com/FluxML/Flux.jl/issues/1935
    data = ObsView(e.data)

    data = e.shuffle ? shuffleobs(e.rng, data) : data
    data = e.batchsize > 0 ? BatchView(data; e.batchsize, e.partial, e.collate) : data

    iter = if e.parallel
        eachobsparallel(data; e.buffer)
    else
        if e.buffer
            buf = getobs(data, 1)
            (getobs!(buf, data, i) for i in 1:numobs(data))
        else
            (getobs(data, i) for i in 1:numobs(data))
        end
    end
    obs, state = iterate(iter)
    return obs, (iter, state)
end

function Base.iterate(::DataLoader, (iter, state))
    ret = iterate(iter, state)
    isnothing(ret) && return
    obs, state = ret
    return obs, (iter, state)
end

function Base.length(e::DataLoader)
    numobs(if e.batchsize > 0
        # Wrapping with ObsView in order to work around
        # issue https://github.com/FluxML/Flux.jl/issues/1935
        data = ObsView(e.data)

        BatchView(data; e.batchsize, e.partial)
    else
        e.data
    end)
end

Base.IteratorEltype(::DataLoader) = Base.EltypeUnknown()

## This causes error in some cases of `collect(loader)`
# function Base.eltype(e::DataLoader)
#     eltype(if e.batchsize > 0
#         BatchView(e.data; e.batchsize, e.partial)
#     else
#         e.data
#     end)
# end

@inline function _dataloader_foldl1(rf, val, e::DataLoader, data)
    if e.shuffle
        _dataloader_foldl2(rf, val, e, shuffleobs(e.rng, data))
    else
        _dataloader_foldl2(rf, val, e, data)
    end
end

@inline function _dataloader_foldl2(rf, val, e::DataLoader, data)
    if e.batchsize > 0
        _dataloader_foldl3(rf, val, e, BatchView(data; e.batchsize, e.partial))
    else
        _dataloader_foldl3(rf, val, e, data)
    end
end

@inline function _dataloader_foldl3(rf, val, e::DataLoader, data)
    if e.buffer > 0
        _dataloader_foldl4_buffered(rf, val, data)
    else
        _dataloader_foldl4(rf, val, data)
    end
end

@inline function _dataloader_foldl4(rf, val, data)
    for i in 1:numobs(data)
        @inbounds x = getobs(data, i)
        # TODO: in 1.8 we could @inline this at the callsite,
        #       optimizer seems to be very sensitive to inlining and
        #       quite brittle in its capacity to keep this type stable
        val = Transducers.@next(rf, val, x)
    end
    Transducers.complete(rf, val)
end

@inline function _dataloader_foldl4_buffered(rf, val, data)
    buf = getobs(data, 1)
    for i in 1:numobs(data)
        @inbounds x = getobs!(buf, data, i)
        val = Transducers.@next(rf, val, x)
    end
    Transducers.complete(rf, val)
end

@inline function Transducers.__foldl__(rf, val, e::DataLoader)
    e.parallel && throw(ArgumentError("Transducer fold protocol not supported on parallel data loads"))
    _dataloader_foldl1(rf, val, e, ObsView(e.data))
end

ext: jl

max_stars_repo_name: JuliaML/MLBase

max_issues_repo_name: JuliaML/MLBase

max_forks_repo_name: JuliaML/MLBase


text:
using Test, GL2017Replication, MAT

@testset "Testing the testing" begin
    @test hello("Julia") == "Hello, Julia"
    @test domath(2.0) ≈ 7.0
end

@testset "Structures" begin
    gl = ModelGL()
    @test typeof(gl.β) == Float64
end

@testset "initilize" begin
    gl = ModelGL()
    initilize!(gl)
    @test sum(gl.pr) ≈ 1.0
end

@testset "EGM" begin
    gl = ModelGL()
    initilize!(gl)
    EGM!(gl)
    @test minimum(gl.b_pol)>=gl.b_grid[1]
    @test maximum(gl.b_pol)<=gl.b_grid[end]
end

@testset "compute distribution" begin
    gl = ModelGL()
    initilize!(gl)
    EGM!(gl)
    compute_distribution!(gl)
    @test sum(gl.JD) ≈ 1
end

@testset "aggregate" begin
    gl = ModelGL()
    initilize!(gl)
    EGM!(gl)
    compute_distribution!(gl)
    aggregate!(gl)
    @test gl.D_4Y_actual < gl.B_4Y_actual
end

cd(@__DIR__)
steady_mat = matopen("steady.mat")
@testset "calibrate" begin
    gl = ModelGL()
    load_parameters_iss!(gl)
    calibrate!(gl)
    @testset "initial" begin
        Y1_mat =  read(steady_mat, "Y1")
        cpol1_mat = read(steady_mat, "c_pol1")
        @test abs(gl.Y_actual - Y1_mat) < 0.001
        @test abs(maximum(gl.c_pol) - maximum(cpol1_mat)) < 0.001
    end

    gl_tss = calibrate_terminal(gl,load_params = true)

    @testset "calibrate_terminal" begin
        Y2_mat =  read(steady_mat, "Y2")
        cpol2_mat = read(steady_mat, "c_pol2")
        @test abs(gl_tss.Y_actual - Y2_mat) < 0.001
        @test abs(maximum(gl_tss.c_pol) - maximum(cpol2_mat)) < 0.001
    end

    #Takes very long to run
    #Tgl = TransGL()
    #gl_trans = transition!(gl,gl_tss,Tgl)
    #
    #@testset "transition" begin
    #    @test abs(Tgl.D_4Y_t[1] - gl.D_4Y_actual) < 0.001
    #    @test abs(Tgl.D_4Y_t[end] - gl_tss.D_4Y_actual) < 0.001
    #end

end


text:
module DESincEig
#=
A module by Philippe Gaudreau,
Department of Mathematical & Statistical Sciences,
University of Alberta, 2014.

The primary function of this module (function SincEigen) computes the eigenvalues of
Singular Sturm-Liouville problems using the DESINC Method.
_________________________________________________________________________________
A general S-L problem has the following form:

P1 : (-D^2 + q(x) ) u(x) = λ ρ(x) u(x),  a < x < b,   with   u(a) = u(b) = 0,
     where -∞ ≦ a < b ≦ ∞ .

In the problem P1,
1. D is the differential operator, an
2. q(x) is a continuous bounded function defined on (a,b)
3. ρ(x) is a continuous positive function defined on (a,b)
4. u(x) are the eigenfunction corresponding to the eigenvalues λ.
_________________________________________________________________________________

In this algorithm, we use the transformation developed by Eggert et al. in reference [1].
This transformation results in a symmetric generalized eigenvalue problem defined as:

P2 : (-D^2 + qtilde(x))v(x) = λ ρtilde(x) v(x),  -∞ < x < ∞,  with   v(±∞) = 0

In the problem P2,
1. D is the differential operator, an
2. qtilde(x) is the resulting transformed function defined on (-∞,∞)
3. ρtilde(x) is the resulting transformed function defined on (-∞,∞)
4. v(x) are the eigenfunction corresponding to the eigenvalues λ.
   v(x) now have double expoenential decay at both infinities.
See reference [2] for more details of the form of qtilde(x) and ρtilde(x)
_________________________________________________________________________________

References:
1.   N. Eggert, M. Jarratt, and J. Lund. Sinc function computation of the eigenvalues of
          Sturm-Liouville problems. SIAM Journal of Computational Physics,
          69:209-229, 1987
2.   P. Gaudreau, R.M. Slevinsky, and H. Safouhi. The Double Exponential Sinc Collocation
          Method for Singular Sturm-Liouville Problems. arXiv:1409.7471v2, 2014
3.   R.M. Slevinsky and S. Olver. On the use of conformal maps for the acceleration of
          convergence of the trapezoidal rule and Sinc numerical methods.
          arXiv:1406:3320, 2014
4.   R.M. Corless, G.H. Gonnet, D.E.G. Hare, D.J. Jeffrey, D.E. Knuth, On the Lambert W function.
          Advances in Computational Mathematics, 5(1):329-359, 1996
5.   P. Gaudreau, and H. Safouhi. Centrosymmetric Matrices in the Sinc Collocation Method for
          Sturm-Liouville Problems. arXiv:1507.06709v1, 2015
=#

using SincFun, Optim

export SincEigen, SincEigenStop

include("lambertW.jl")

####################### Conformal mappings based on different domains #####################################
#=
In this first section of the DESincEig module, we will construct various conformal mappings to induce
a double expoential decay at both infinities for the solution of the transformed S-L problem.

These conformal mappings are based on the domain of the orginal problem ( see problem P1 above )
as well as the type of decay of the solution u(x) at the end-points of its domain.

The conformal mappings, ϕ(t), will be created as a composition of an outer map, ψ(t), with an inner map, H(t).
In other words, the conformal mappings have the form:

ϕ(t) = (ψ∘H)(t) = ψ(H(t)),
_____________________________________ DEFINITION OF OUTER MAPS________________________________________
The SincFun type Domain is used to select from the outer maps depending on the domain of the problem.
Concrete Domains have methods:
 1. ψ  : The outer map
 2. ψp : The derivative of ψ, (ψ prime)
 3. ψtilde : The first half of the resulting transformation proposed by Eggert et al.
            ψtilde =  - sqrt{ψp} D( 1/ψp D(\sqrt{ψp}) ),
            where D is the differential operator w.r.t the varibale t.

1. For S-L problems on a finite domain I=(a,b) with algebraic decay at the endpoints:
Domain = Finite(a,b).
2. For S-L problems on a infinite domain I=(-∞,∞) with algebraic decay at the endpoints:
Domain = Infinite1{Float64}().
3. For S-L problems on a infinite domain I=(-∞,∞) with single-exponential decay at the endpoints:
Domain = Infinite2{Float64}().
4. For S-L problems on a semi-infinite domain I=(0,∞) with single-exponential decay at infinity and algebraic decay at 0:
Domain = SemiInfinite1{Float64}().
=#

type SincResults{T}
     RESULTS::Matrix{T}
     All_Abs_Error_Approx
     hoptimal::Vector{T}
     N::Vector{T}
     MatrixSizes
end

export ψtilde
ψtilde(d::Finite,t) = 0*t+1
ψtilde(d::Infinite1,t) = 0.25.-0.75.*sech(t).^2
ψtilde(d::Infinite2,t) = 0*t
ψtilde(d::SemiInfinite1,t) = (2*exp(t).+1)./(2.+2.*exp(t)).^2
ψtilde(d::SemiInfinite2,t) = 0*t + 0.25

#=_____________________________________ DEFINITION OF INNER MAP________________________________________
Inner map is given by a SincFun type ConformalMap.
It stores the data for h(t) = u₀sinh(t) + u₁ + u₂t + u₃ t^2 + ⋯ + uₙt^(n-1).
For more information on the form of the inner map, please see reference [3].
=#

################################### The Main function: SincEigen ####################################################
#=
This function computes the eigenvalues of the the Sturm-Liouville (P1) using the double exponential Sinc-Collocation method.
For optimal results, it is worth performing an asymptotic analysis of the solution of P2.
The resulting analysis will lead to the following bounds for constants βL,βR,γL,γR>0:
 |v(x)| < AL exp(-βL exp ( γL |x| ) ) , on (-∞,0)
 |v(x)| < AR exp(-βR exp ( γR |x| ) ) , on (0, ∞)
The parameters βL, βR, γL, γR are used in the calculation of the mesh size, h, for the DESCM.

Another important parameter involved in the calculations of the mesh size for the DESCM is the width, d, of the strip D_{d}.
Let S denote the set of complex singularities of the functions qtilde(x) and ρtilde(x):
S = { z ∈ C : qtilde(z) or ρtilde(z) does not exist a.}
let s denote the positive imaginary part of nearest singularity to the real axis:
s = min | Im{S} |,
then d = min{ π/2max{γL,γR} , s } .

Input:
Necessary parameters
1. q(x):: Function,      The function in P1
2. ρ(x):: Function,      The function in P1
3. domain:: Domain,      Finite, Infinite1, Infinite2 or SemiInfinite1
4. βopt:: Vector{T},     [βL,βR]
5. γopt:: Vector{T},     [γL,γR]
6. d:: Number,           min{ π/2max{γL,γR} , s }
____________________________________________________________________________________________________________________
Extra parameters ( not necessary but can offer more options )
7. enum::Tuple{Int,T}         (n,tv): Default = (0,NaN)
                            n ∈ {0,1,2,3,...} :  Eigenvalue number
                            tv = λ_{n} : True Value for eigenvalue λ_{n}
                            If specified, the absolute error is computed instead of the absolute error approximation
                            of the DESCM w.r.t. to λ_{n} will be returned.

8. tol::T,                  tolerance level: Default = 5e-12
                            tol is used to find the optimal Matrix Size in order for the approximate eigenvalue
                            μ to have an an aprroximation to the absolute error less than tol.

9. Range::Vector{Integer},  [start:skip:end]: Default = [1:1:100]
                            Will run algorithm with N or M = "start" with index jump "skip" until "end".

10. u0::T                   Parameter used in the inner map construction H(t) presented above.
                            Default = one(T)

11. u::Vector{T}            Parameters used in the inner map construction H(t) presented above.
                            Default = [zero(T)]

12. Trace_Mesh::Bool         Default = false
If true, the mesh size will be computed by minimizing the trace of the matrix D^(-1)HD^(-1).
resulting from the DESCM. For even functions q(x) and ρ(x) and an infinite domain: Domain = Infinite1 or Infinite2,
one can minimize this functional to obtain an alternate mesh size: htilde. This alternate mesh-size has proven
to be better suited for highly-oscillatory functions q(x). This functional is minimized using the Julia package: Optim.

13. Centro::Bool              Default = false
As demonstrated in reference [5], if the function q(x) and ρ(x) are even on an infinite domain, the sinc matices are
also centrosymmetric. Hence, If Centro is set to true, the simpifications detailed in reference [5] are implemented.
_______________________________________________________________________________________________________________________
=#
function SincEigen{T<:Number,S<:Integer}(q::Function,ρ::Function,domain::Domain{T},β::Vector{T},γ::Vector{T},dopt::T; enum::Tuple{S,Any}=(0,NaN), tol::T=5e-12, Range::Vector{S}=collect(1:1:100), u0::T=one(T), u::Vector{T}=[zero(T)], Trace_Mesh::Bool=false,Centro::Bool=false)
    # Functions used in Matrix construction.
    H = [ConformalMap(u0,u)]
    for i=1:3 push!(H,H[end]') end
    ϕ(t) = ψ(domain,H[1][t])
    ϕp2(t) = (ψp(domain,H[1][t]).*H[2][t]).^2
    qtilde(t) = ψtilde(domain,H[1][t]).*H[2][t].^2 .+ (3/4).*(H[3][t]./H[2][t]).^2 .-(H[4][t]./2H[2][t]) .+ q(ϕ(t)).*ϕp2(t)
    rhotilde(t) = ρ(ϕ(t)).*ϕp2(t)
    # Determining step sizes and left and right collocation points based on asymptotic information.
    if γ[1] == γ[2] && β[1] == β[2]
          if Trace_Mesh == true
               diag_sinc_matrix(t) = (ψtilde(domain,H[1][t] .+ (3/4).*(H[3][t]./H[2][t].^2).^2 .-(H[4][t]./2H[2][t].^3))./ψp(domain,H[1][t]).^2 .+ q(ϕ(t)))./ρ(ϕ(t))
               hoptimal = [optimize(h->sum( diag_sinc_matrix(collect(-Range[i]:Range[i])*h).+ (pi^2/(3h^2))./ rhotilde(collect(-Range[i]:Range[i])*h) ),0.001,(3.0*u0+log(pi*dopt*γ[2]*i/β[2]))./(γ[2]*i)).minimum for i in collect(1:length(Range))]
          elseif Trace_Mesh == false
               hoptimal = lambertW(pi*dopt*γ[2]*Range/β[2])./(γ[2]*Range)
          end #if loop
    else
          hoptimal = max(lambertW(pi*dopt*γ[2]*Range/β[2])./(γ[2]*Range),lambertW(pi*dopt*γ[1]*Range/β[1])./(γ[1]*Range))
    end # if loop
     N = convert(Vector{T},round(log(-log(sqrt(hoptimal).*exp(-pi*dopt./hoptimal))/β[2])/γ[2]./hoptimal))
     M = convert(Vector{T},round(log(-log(sqrt(hoptimal).*exp(-pi*dopt./hoptimal))/β[1])/γ[1]./hoptimal))
    #  INITIAL CONDITIONS
    Length = length(Range)                                            # Length of the vectors in the following algorithm given the previous conditions
    if Centro==false
          MatrixSizes = N.+M.+1                                       # Array of all square matrix sizes for every iteration.
          All_Eigenvalues = zeros(round(Int,MatrixSizes[end]),Length) # Matrix used for storing all obtained eigenvalues as columns at every iteration.
          for i = 1:Length
               h = hoptimal[i]                                        # Mesh size used at every iteration.
               k = collect(-M[i]:N[i])                                # Vector of Mesh points.
               A = Symmetric(diagm(qtilde(k*h))-sinc(2,k'.-k)/h^2)    # Construct symmetric A matrix.
               D2 = Symmetric(diagm(rhotilde(k*h)))                   # Construct Diagonal pos. def. matrix D^2
               E = eigvals(A,D2)                                      # Solving Generalized eigenvalue problem.
               All_Eigenvalues[1:length(E),i] = E                     # Storing all eigenvalues in columns in matrix All_Eigenvalues.
          end
          # Calculating an Approximation to the Absolute Error for all Energy values
          if isnan(enum[2]) == true
               All_Abs_Error_Approx  = abs(All_Eigenvalues[:,2:end].-All_Eigenvalues[:,1:end-1])
          else
               All_Abs_Error_Approx = zeros(All_Eigenvalues[:,2:end])
               Index = [collect(1:enum[1]);collect((enum[1]+2):Length)]
               All_Abs_Error_Approx[Index,:] =  abs(All_Eigenvalues[Index,2:end].-All_Eigenvalues[Index,1:(end-1)])
               All_Abs_Error_Approx[enum[1]+1,:] =  abs(All_Eigenvalues[enum[1]+1,2:end] .- enum[2])
          end # if loop
          ## Ouputing the convergence analysis of the algorithm given the number of iterations and the tolerance level tol.
          RESULTS = Convergence_Analysis(All_Eigenvalues,tol,MatrixSizes,All_Abs_Error_Approx)
          return SincResults(RESULTS, All_Abs_Error_Approx, hoptimal, N, MatrixSizes)
    elseif Centro==true
          All_Eigenvalues_even = zeros(round(Int,Range[end]+1),Length)
          All_Eigenvalues_odd = zeros(round(Int,Range[end]),Length)
          for i = 1:Length
               h = hoptimal[i]
               index = convert(Int,N[i])
               k = collect(-index:-1).*1.0
               j = collect(-index:index).*1.0
               Temp = -sinc(2,k'.-j)./h^2
               A = diagm(qtilde(k*h)) .+ Temp[1:index,:]
               JC = flipdim(Temp[index+2:2index+1,:],1)
               qc = qtilde(0.0).+ (pi^2/(3h^2))
               x = sqrt(2).*Temp[index+1,:]
               # Solving Generalized eigenvalue problems.
               E_odd = eigvals(Symmetric(A-JC),Symmetric(diagm(rhotilde(k*h))))
               E_even = eigvals(Symmetric([qc  x' ; x A+JC]),Symmetric(diagm(rhotilde([0.0;k]*h))))
               All_Eigenvalues_even[1:length(E_even),i] = E_even
               All_Eigenvalues_odd[1:length(E_odd),i] = E_odd
          end
          # Calculating an Approximation to the Absolute Error for all Energy values
          All_Abs_Error_Approx_even  = abs(All_Eigenvalues_even[:,2:end].-All_Eigenvalues_even[:,1:end-1])
          All_Abs_Error_Approx_odd  = abs(All_Eigenvalues_odd[:,2:end].-All_Eigenvalues_odd[:,1:end-1])
          if isnan(enum[2]) == false
            if  isodd(enum[1])
                index=round(Int,(enum[1]+1)/2)
                All_Abs_Error_Approx_odd[index+1,:] =  abs(All_Eigenvalues_odd[index+1,2:end] .- enum[2])
            elseif iseven(enum[1])
                index=round(Int,(enum[1]+2)/2)
                All_Abs_Error_Approx_even[index+1,:] =  abs(All_Eigenvalues_even[index+1,2:end] .- enum[2])
            end # if loop
          end
          ## Ouputing the convergence analysis of the algorithm given the number of iterations and the tolerance level tol.
          RESULTS_odd = Convergence_Analysis(All_Eigenvalues_odd,tol,N+0.0,All_Abs_Error_Approx_odd)
          RESULTS_even = Convergence_Analysis(All_Eigenvalues_even,tol,N.+1.0,All_Abs_Error_Approx_even)
          RESULTS_odd[:,1] = 2.*RESULTS_odd[:,1].+1
          RESULTS_even[:,1] = 2.*RESULTS_even[:,1]
          RESULTS = sortrows([RESULTS_odd;RESULTS_even])
          return SincResults(RESULTS, (All_Abs_Error_Approx_even,All_Abs_Error_Approx_odd), hoptimal, N, (N+1,N))
     end
end

######################################################################################################################################

################################### Convergence Analysis ####################################################
#=
The following function performs a convergence analysis of the DESCM. More specifically, After running the function
SincEigen for the values "Range", the function Convergence_Analysis will find all eigenvalues which achieved
an approximation to the absolute error less than a prediscribed accuray "tol".
Input: All_Eigenvalues::Matrix{T},         Matrix containing all eigenvalues computed by running from Range[1] to Range[end]
       tol::T,                             Tolerance level
       MatrixSizes::Vector{T},             Vector of matrix size used in the algorithm SincEigen
       enum::Vector{T},                    Vector of given eigenvalue number and corresponding eigenvalue if known.
       All_Abs_Error_Approx::Matrix{T}     Matrix containing the sequence of approximations to the absolute error for every eigenvalue obtained in the SincEigen algorithm
                                           Row i, Column j = approximations to the absolute error for eigenvalue number "n=i-1" when the matrix size = MatrixSizes[j].

Output:
RESULTS=[Eig_number MatrixSizeOpt Eigenvalues Error]::Matrix{T},
    Matrix of eigenvalues wich achieved an approximation to the absolute error less than a prediscribed accuray "tol"
    First column: The eigenvalue number,
    Second column: The optimal matrix size needed to achieve "tol".
    Third column: Estimated eigenvalue
    Fourth column: Approximation to the absolute error
=#

function  Convergence_Analysis{T<:Number}(All_Eigenvalues::Matrix{T},tol::T,MatrixSizes::Vector{T},All_Abs_Error_Approx::Matrix{T})
    #Finding which Eigenvalues satisfy the condition Absolute Error Approxmiation < tol
    m = size(All_Abs_Error_Approx,1)
    First_Non_Zero = [findfirst(All_Abs_Error_Approx[i,:]) for i in collect(1:m)]
    Less_than_tol_matrix =  All_Abs_Error_Approx .< tol
    Eig_less_than_tol_posi = [findnext(Less_than_tol_matrix[i,:], true, First_Non_Zero[i]) for i in collect(1:m)]
    # list of eigenvalues  that satisfy the condition Relative Error Approxmiation<tol
    Posi_index = findn(Eig_less_than_tol_posi)#[1]
    # number of eigenvalues that satisfy the condition Error<tol
    num_eig_less_than_tol = length(Posi_index)
    # Find the optimal value of N for the eignvalues from the list Posi_index
    idx = [Eig_less_than_tol_posi[Posi_index[j]]+1 for j in collect(1:num_eig_less_than_tol)]
    # Find the value of the eigenvalues and Relative Error Approximation for the eignvalues from the list Posi_index
    Error = [All_Abs_Error_Approx[Posi_index[i],idx[i]-1] for i in collect(1:num_eig_less_than_tol)]
    Eigenvalues = [All_Eigenvalues[Posi_index[i],idx[i]] for i in collect(1:num_eig_less_than_tol)]
    # Display wanted results
    Eig_number = Posi_index.-1
    MatrixSizeOpt = MatrixSizes[idx]
    [Eig_number MatrixSizeOpt Eigenvalues Error]
end
######################################################################################################################################

function SincEigenStop{T<:Number}(q::Function , ρ::Function , domain::Domain{T} , β::Vector{T} ,γ::Vector{T} ,dopt::T; enum::Vector{T}=[zero;NaN],tol::T = convert(T,5.0e-12),u0::T = one(T),u::Vector{T}=[zero(T)])
# Functions used in Matrix construction.
    H = [ConformalMap(u0,u)]
    for i=1:3 push!(H,H[end]') end
    ϕ(t) = ψ(domain,H[1][t])
    ϕp2(t) = (ψp(domain,H[1][t]).*H[2][t]).^2
    qtilde(t) = ψtilde(domain,H[1][t]).*H[2][t].^2 .+ (3/4).*(H[3][t]./H[2][t]).^2 .-(H[4][t]./2H[2][t]) .+ q(ϕ(t)).*ϕp2(t)
    rhotilde(t) = ρ(ϕ(t)).*ϕp2(t)
# Determining step sizes and left and right collocation points based on asymptotic information.
    n = max( floor((enum[1]+1)/2) , 1 )
    h,k = h_and_k(n,β,γ,dopt)
    A = Symmetric(diagm(qtilde(k*h))-sinc(2,k'.-k)/h^2)
    D2 = Symmetric(diagm(rhotilde(k*h)))
    E = eigvals(A,D2)
    Eigenvalue_enum_old = E[enum[1]+1]

    n += 1

    h,k = h_and_k(n,β,γ,dopt)
    A = Symmetric(diagm(qtilde(k*h))-sinc(2,k'.-k)/h^2)
    D2 = Symmetric(diagm(rhotilde(k*h)))
    E = eigvals(A,D2)
    Eigenvalue_enum_new = E[enum[1]+1]
    Abs_Error_enum = isnan(enum[2]) ? abs(Eigenvalue_enum_new- Eigenvalue_enum_old) : abs(Eigenvalue_enum_new - enum[2])

    while Abs_Error_enum > tol && n < 200
    Eigenvalue_enum_old = Eigenvalue_enum_new
    n += 1
    h,k = h_and_k(n,β,γ,dopt)
    A = Symmetric(diagm(qtilde(k*h))-sinc(2,k'.-k)/h^2)
    D2 = Symmetric(diagm(rhotilde(k*h)))
    E = eigvals(A,D2)
    Eigenvalue_enum_new = E[enum[1]+1]
    Abs_Error_enum = isnan(enum[2]) ? abs(Eigenvalue_enum_new- Eigenvalue_enum_old) : abs(Eigenvalue_enum_new - enum[2])
    end #while loop

    v = eigvecs(A,D2)[:,enum[1]+1]
    return (enum[1],round(Int,length(k)),Eigenvalue_enum_new,Abs_Error_enum,v,h,k)

end #function

function h_and_k{T<:Number}(n::T, β::Vector{T} ,γ::Vector{T} ,dopt::T)
if γ[1]>γ[2]
    h = lambertW(pi*dopt*γ[1]*n/β[1])./(γ[1]*n)
    N = max( ceil( (γ[1]/γ[2]).*n .+ log(β[1]/β[2])./ (γ[2].*h) ) , 0 )
    return h,collect(-n:N)
elseif γ[2]>γ[1]
    h = lambertW(pi*dopt*γ[2]*n/β[2])./(γ[2]*n)
    M = max( floor( (γ[2]/γ[1]).*n .+ log(β[2]/β[1])./ (γ[1].*h) ) , 0 )
    return h, collect(-M:n)
elseif γ[1] == γ[2] && β[1] > β[2]
    h = lambertW(pi*dopt*γ[1]*n/β[1])./(γ[1]*n)
    N = ceil( (γ[1]/γ[2]).*n .+ log(β[1]/β[2])./ (γ[2].*h) )
    return h, collect(-n:N)
elseif γ[1] == γ[2] && β[1] < β[2]
    h = lambertW(pi*dopt*γ[2]*n/β[2])./(γ[2]*n)
    M = floor( (γ[2]/γ[1]).*n .+ log(β[2]/β[1])./ (γ[1].*h) )
    return h , collect(-M:n)
elseif γ[1] == γ[2] && β[1] == β[2]
    h = lambertW(pi*dopt*γ[2]*n/β[2])./(γ[2]*n)
    return h , collect(-n:n)
end  # if loop
end

end #Module


text:
"""
    lautat(surf, curfield, nsteps=500, dtstar=0.015, startflag=0,
writeflag=0, writeInterval=1000., delvort=delNone(), maxwrite=100,
nround=6)

Simulates potential flow for an airfoil undergoing unsteady motion
    using Large-Angle Unsteady Thin-Airfoil Theory.

    Inputs required:
    - `TwoDSurf` : based on airfoil geometry and kinematics
    - `TwoDFlowField` : initial flowfield

    `dtstar` should be in accordance with the kinematics; see `?find_tstep`.

    `nsteps` should be calculated according to total run time required and
    dtstar.

    `startflag=1` can be used to continue a simulation using the
    *resultsSummary* file in current directory.

    `writeflag=1` along with a suitable writeInterval is used to write flow
    data at intervals. These can be used to create plots and movies of the
    flowfield.

    `delvort` is used to control the vortex count. The simulation cost
    increases as O(n``^2``) with number of vortices in the flow field and so
    long run times can get very computationally expensive. Look at
    `?delVortDef` for algorthims to control vortex count.

    `maxwrite` is the maximum number of timestamps/directories that may be
    written.

    `nround` is the number of digits to which the output timestamps are
    rounded off to.

    The Large-Angle Unsteady Thin-Airfoil theory is described in: Ramesh,
    K. et al., "An unsteady airfoil theory applied to pitching motions
    validated against experiment and computation", Theor. Comput. Fluid
    Dyn. (2013) 27: 843. [Weblink](https://doi.org/10.1007/s00162-012-0292-8)

"""
function lautat(surf::TwoDSurf, curfield::TwoDFlowField, nsteps::Int64 = 500, dtstar::Float64 = 0.015, startflag = 0, writeflag = 0, writeInterval = 1000., delvort = delNone(); maxwrite = 50, nround=6, wakerollup=1)

    # If a restart directory is provided, read in the simulation data
    if startflag == 0
        mat = zeros(0, 8)
        t = 0.
    elseif startflag == 1
        dirvec = readdir()
        dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec)
        latestTime = maximum(dirresults)
        mat = DelimitedFiles.readdlm("resultsSummary",header = true)
        t = mat[end,1]
    else
        throw("invalid start flag, should be 0 or 1")
    end
    mat = mat'

    dt = dtstar*surf.c/surf.uref

    # if writeflag is on, determine the timesteps to write at
    if writeflag == 1
        writeArray = Int64[]
        tTot = nsteps*dt
        for i = 1:maxwrite
            tcur = writeInterval*real(i)
            if t > tTot
                break
            else
                push!(writeArray, Int(round(tcur/dt)))
            end
        end
    end

    #Intialise flowfield
    for istep = 1:nsteps
        #Udpate current time
        t = t + dt

        #Update kinematic parameters
        update_kinem(surf, t)

        #Update bound vortex positions
        update_boundpos(surf, dt)

        #Add a TEV with dummy strength
        place_tev(surf,curfield,dt)

        #Solve for TEV strength to satisfy Kelvin condition
        kelv = KelvinCondition(surf,curfield)
        soln = nlsolve(not_in_place(kelv), [-0.01])
        curfield.tev[length(curfield.tev)].s = soln.zero[1]

        #Update adot
        update_a2a3adot(surf,dt)

        #Set previous values of aterm to be used for derivatives in next time step
        surf.a0prev[1] = surf.a0[1]
        for ia = 1:3
            surf.aprev[ia] = surf.aterm[ia]
        end

        #Update rest of Fourier terms
        update_a2toan(surf)

        #Calculate bound vortex strengths
        update_bv(surf)

        # Delete or merge vortices if required
        controlVortCount(delvort, surf.bnd_x[Int(round(surf.ndiv/2))], surf.bnd_z[Int(round(surf.ndiv/2))], curfield)

        #Wake rollup if flag on
        if wakerollup == 1
            wakeroll(surf, curfield, dt)
        end

        # Calculate force and moment coefficients
        cl, cd, cm = calc_forces(surf, [curfield.u[1], curfield.w[1]])

        # write flow details if required
        if writeflag == 1
            if istep in writeArray
                dirname = "$(round(t, sigdigits=nround))"
                writeStamp(dirname, t, surf, curfield)
            end
        end

        # for writing in resultsSummary
        mat = hcat(mat,[t, surf.kinem.alpha, surf.kinem.h, surf.kinem.u, surf.a0[1], cl, cd, cm])

    end

    mat = mat'

    f = open("resultsSummary", "w")
    println(f,"#time \t","alpha(deg) \t","h/c \t", "u/uref \t", "A0 \t", "Cl \t", "Cd \t", "Cm ")
    DelimitedFiles.writedlm(f, mat)
    close(f)

    mat, surf, curfield

end

function ldvm(surf::TwoDSurf, curfield::TwoDFlowField, nsteps::Int64 = 500, dtstar::Float64 = 0.015, startflag = 0, writeflag = 0, writeInterval = 1000., delvort = delNone(); maxwrite = 50, nround=6)

    # If a restart directory is provided, read in the simulation data
    if startflag == 0
        mat = zeros(0, 8)
        t = 0.
    elseif startflag == 1
        dirvec = readdir()
        dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec)
        latestTime = maximum(dirresults)
        mat = DelimitedFiles.readdlm("resultsSummary",header = true)
        t = mat[end,1]
    else
        throw("invalid start flag, should be 0 or 1")
    end
    mat = mat'

    dt = dtstar*surf.c/surf.uref

    # if writeflag is on, determine the timesteps to write at
    if writeflag == 1
        writeArray = Int64[]
        tTot = nsteps*dt
        for i = 1:maxwrite
            tcur = writeInterval*real(i)
            if t > tTot
                break
            else
                push!(writeArray, Int(round(tcur/dt)))
            end
        end
    end

    # time stepping
    for istep = 1:nsteps
        #Udpate current time
        t = t + dt

        #Update external flowfield
        update_externalvel(curfield, t)

        #Update kinematic parameters
        update_kinem(surf, t)

        #Update bound vortex positions
        update_boundpos(surf, dt)

        #Add a TEV with dummy strength
        place_tev(surf,curfield,dt)

        #Update incduced velocities on airfoil
        update_indbound(surf, curfield)

        #Solve for TEV strength to satisfy Kelvin condition
        kelv = KelvinCondition(surf,curfield)
        soln = nlsolve(not_in_place(kelv), [-0.01])
        curfield.tev[length(curfield.tev)].s = soln.zero[1]

        #Update adot
        update_a2a3adot(surf,dt)

        lesp = surf.a0[1]

        #Check for LESP condition
        if abs(lesp)>surf.lespcrit[1]
            #2D iteration if LESP_crit is exceeded

            #Add a LEV with dummy strength
            place_lev(surf,curfield,dt)

            #Solve for TEV and LEV strengths to satisfy Kelvin condition and Kutta condition at leading edge
            kelvkutta = KelvinKutta(surf,curfield)
            soln = nlsolve(not_in_place(kelvkutta), [-0.01; 0.01])
            (curfield.tev[length(curfield.tev)].s, curfield.lev[length(curfield.lev)].s) = soln.zero[1], soln.zero[2]

            #set flag for levshedding=on
            surf.levflag[1] = 1
        else
            surf.levflag[1] = 0
        end

        #Update rest of Fourier terms
        update_a2toan(surf)

        #Set previous values of aterm to be used for derivatives in next time step
        surf.a0prev[1] = surf.a0[1]
        for ia = 1:3
            surf.aprev[ia] = surf.aterm[ia]
        end

        #Calculate bound vortex strengths
        update_bv(surf)

        # Delete or merge vortices if required
        controlVortCount(delvort, surf.bnd_x[Int(round(surf.ndiv/2))], surf.bnd_z[Int(round(surf.ndiv/2))], curfield)

        # free wake rollup
        wakeroll(surf, curfield, dt)

        # Calculate force and moment coefficients
        cl, cd, cm, cn = calc_forces(surf, [curfield.u[1], curfield.w[1]])

        # write flow details if required
        if writeflag == 1
            if istep in writeArray
                dirname = "$(round(t, sigdigits=nround))"
                writeStamp(dirname, t, surf, curfield)
            end
        end

        # for writing in resultsSummary
        mat = hcat(mat,[t, surf.kinem.alpha, surf.kinem.h, surf.kinem.u, surf.a0[1], cl, cd, cm])

    end

    mat = mat'

    f = open("resultsSummary", "w")
    println(f,"#time \t","alpha(deg) \t","h/c \t", "u/uref \t", "A0 \t", "Cl \t", "Cd \t", "Cm ")
    DelimitedFiles.writedlm(f, mat)
    close(f)

    mat, surf, curfield

end

function ldvm(surf::Vector{TwoDSurf}, curfield::TwoDFlowField, nsteps::Int64 = 500, dtstar::Float64 = 0.015, startflag = 0, writeflag = 0, writeInterval = 1000., delvort = delNone(); maxwrite = 50, nround=6)

    nsurf = length(surf)

    # If a restart directory is provided, read in the simulation data
    if startflag == 0
        mat = zeros(0, 7*nsurf+1)
        t = 0.
    elseif startflag == 1
        dirvec = readdir()
        dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec)
        latestTime = maximum(dirresults)
        mat = DelimitedFiles.readdlm("resultsSummary",header = true)
        t = mat[end,1]
    else
        throw("invalid start flag, should be 0 or 1")
    end
    mat = mat'

    dt = dtstar*minimum(map(q->q.c/q.uref, surf))

    # if writeflag is on, determine the timesteps to write at
    if writeflag == 1
        writeArray = Int64[]
        tTot = nsteps*dt
        for i = 1:maxwrite
            tcur = writeInterval*real(i)
            if t > tTot
                break
            else
                push!(writeArray, Int(round(tcur/dt)))
            end
        end
    end

    lesp = zeros(nsurf)
    cl = zeros(nsurf)
    cd = zeros(nsurf)
    cm = zeros(nsurf)

    # time stepping
    for istep = 1:nsteps
        #Udpate current time
        t = t + dt

        #Update external flowfield
        update_externalvel(curfield, t)

        for is = 1:nsurf
            #Update kinematic parameters
            update_kinem(surf[is], t)

            #Update bound vortex positions
            update_boundpos(surf[is], dt)
        end

        #Add TEVs with dummy strength
        place_tev(surf,curfield,dt)

        for is = 1:nsurf
            #Update induced velocities on airfoil
            update_indbound(surf[is], curfield)
            for js = 1:nsurf
                if is != js
                    #add influence of other surfaces to downwash
                    add_indbound_b(surf[is], surf[js])
                end
            end
        end

        #Solve for TEV strength to satisfy Kelvin condition
        kelv = KelvinConditionMult(surf,curfield)
        soln = nlsolve(not_in_place(kelv), ones(nsurf)*-0.01)

        for i = 1:nsurf
            curfield.tev[end-nsurf+i].s = soln.zero[i]

            #Update adot
            update_a2a3adot(surf[i],dt)

            lesp[i] = surf[i].a0[1]
        end

        shed_ind = Int[]

        #Check for LESP condition
        for i = 1:nsurf
            if (abs(lesp[i])>surf[i].lespcrit[1])

                push!(shed_ind, i)
            end
        end
        nshed = length(shed_ind)

        if nshed > 0
            #2D iteration if LESP_crit is exceeded
            #Add LEVs with dummy strength
            place_lev(surf,curfield,dt,shed_ind)

            #Solve for TEV and LEV strengths to satisfy Kelvin condition and Kutta condition at leading edge
            kelvkutta = KelvinKuttaMult(surf,curfield,shed_ind)
            soln = nlsolve(not_in_place(kelvkutta), -0.01*ones(nshed+nsurf))

            for i = 1:nsurf
                if i in shed_ind
                    if shed_ind != [1] && length(shed_ind)==1 #incase LEV only formed on 2nd surface
                    surf[i].levflag[1] = 1
                    i=i-1
                    curfield.tev[end-nsurf+i].s = soln.zero[i]
                    curfield.lev[end-nsurf+shed_ind[i]].s = soln.zero[nsurf+i]
                else
                    curfield.tev[end-nsurf+i].s = soln.zero[i]
                    curfield.lev[end-nsurf+shed_ind[i]].s = soln.zero[nsurf+i]
                    surf[i].levflag[1] = 1
                else
                    curfield.tev[end-nsurf+i].s = soln.zero[i]
                    surf[i].levflag[1] = 0
                end
            end
        else
            for i = 1:nsurf
                surf[i].levflag[1] = 0
            end
        end

        for i = 1:nsurf
            #Set previous values of aterm to be used for derivatives in next time step
            surf[i].a0prev[1] = surf[i].a0[1]
            for ia = 1:3
                surf[i].aprev[ia] = surf[i].aterm[ia]
            end
        end

        # Delete or merge vortices if required
        mean_bndx = sum(map(q->q.bnd_x[Int(round(q.ndiv/2))], surf))/nsurf
        mean_bndz = sum(map(q->q.bnd_z[Int(round(q.ndiv/2))], surf))/nsurf
        controlVortCount(delvort, mean_bndx, mean_bndz, curfield)

        # free wake rollup
        wakeroll(surf, curfield, dt)

        # Calculate force and moment coefficients
        for i = 1:nsurf
            cl[i], cd[i], cm[i] = calc_forces(surf[i], [curfield.u[1], curfield.w[1]])
        end

        # write flow details if required
        if writeflag == 1
            if istep in writeArray
                dirname = "$(round(t, sigdigits=nround))"
                writeStamp(dirname, t, surf, curfield)
            end
        end

        # for writing in resultsSummary
        matvect = [t;]
        for is = 1:nsurf
            matvect = [matvect;[surf[is].kinem.alpha, surf[is].kinem.h,
                                surf[is].kinem.u, surf[is].a0[1], cl[is], cd[is], cm[is]]]
        end
        mat = hcat(mat,matvect)

    end

mat = mat'

f = open("resultsSummary", "w")
println(f,"#time \t","alpha(deg) \t","h/c \t", "u/uref \t", "A0 \t", "Cl \t", "Cd \t", "Cm ")
DelimitedFiles.writedlm(f, mat)
close(f)

mat, surf, curfield

end

function ldvmLin(surf::TwoDSurf, curfield::TwoDFlowField, nsteps::Int64 = 500, dtstar::Float64 = 0.015, startflag = 0, writeflag = 0, writeInterval = 1000., delvort = delNone(); maxwrite = 50, nround=6)

    # If a restart directory is provided, read in the simulation data
    if startflag == 0
        mat = zeros(0, 8)
        t = 0.
    elseif startflag == 1
        dirvec = readdir()
        dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec)
        latestTime = maximum(dirresults)
        mat = DelimitedFiles.readdlm("resultsSummary",header=true)
        t = mat[end,1]
    else
        throw("invalid start flag, should be 0 or 1")
    end
    mat = mat'

    dt = dtstar*surf.c/surf.uref

    # if writeflag is on, determine the timesteps to write at
    if writeflag == 1
        writeArray = Int64[]
        tTot = nsteps*dt
        for i = 1:maxwrite
            tcur = writeInterval*real(i)
            if t > tTot
                break
            else
                push!(writeArray, Int(round(tcur/dt)))
            end
        end
    end

    vcore = 0.02*surf.c

    T1 = zeros(surf.ndiv)
    T2 = zeros(surf.ndiv)
    T3 = zeros(surf.ndiv)

    # time stepping
    for istep = 1:nsteps
        #Udpate current time
        t = t + dt

        #Update external flowfield
        update_externalvel(curfield, t)

        #Update kinematic parameters
        update_kinem(surf, t)

        #Update bound vortex positions
        update_boundpos(surf, dt)

        #Update induced velocities on airfoil
        update_indbound(surf, curfield)

        #Calculate downwash
        update_downwash(surf, [curfield.u[1],curfield.w[1]])

        #The integrals I1 and J1 are based on this downwash which doesnt
        #include vortives shed at current step (equal to term T1)
        #I1 has units if circulation, J1 is dimensionless
        T1[:] = surf.downwash[:]
        I1 = surf.c*simpleTrapz(T1.*(cos.(surf.theta) .- 1. ), surf.theta)
        J1 = -simpleTrapz(T1,surf.theta)/(surf.uref*pi)

        # T2 depends on recenetly shed TEV
        ntev = length(curfield.tev)

        if ntev == 0
            xloc_tev = surf.bnd_x[surf.ndiv] + 0.5*surf.kinem.u*dt
            zloc_tev = surf.bnd_z[surf.ndiv]
        else
            xloc_tev = surf.bnd_x[surf.ndiv]+(1. /3.)*(curfield.tev[ntev].x - surf.bnd_x[surf.ndiv])
            zloc_tev = surf.bnd_z[surf.ndiv]+(1. /3.)*(curfield.tev[ntev].z - surf.bnd_z[surf.ndiv])
        end

        for ib = 1:surf.ndiv
            xdist = surf.bnd_x[ib] - xloc_tev
            zdist = surf.bnd_z[ib] - zloc_tev
            distsq = xdist*xdist + zdist*zdist
            T2[ib] = (surf.cam_slope[ib]*zdist + xdist)/(2*pi*sqrt(distsq^2 + vcore^4))
        end

        #sig_prev = sum(map(q->q.s, curfield.tev)) + sum(map(q->q.s, curfield.lev))
        sig_prev = -surf.uref*surf.c*pi*(surf.a0prev[1] + surf.aprev[1]/2. )

        I2 = simpleTrapz(T2.*(cos.(surf.theta) .- 1. ), surf.theta)
        J2 = -simpleTrapz(T2, surf.theta)/(pi*surf.uref)

        tevstr = -(I1 + sig_prev)/(1 + I2)

        #Calc first 3 fourier coefficients and derivatives
        surf.a0[1] = J1 + J2*tevstr
        for ia = 1:surf.naterm
            surf.aterm[ia] = 2. *(simpleTrapz(T1.*cos.(ia*surf.theta), surf.theta) + tevstr*simpleTrapz(T2.*cos.(ia*surf.theta), surf.theta))/(pi*surf.uref)
        end

        #Calculate adot
        surf.a0dot[1] = (surf.a0[1] - surf.a0prev[1])/dt
        for ia = 1:surf.naterm
            surf.adot[ia] = (surf.aterm[ia]-surf.aprev[ia])/dt
        end

        #Check if LEV shedding is true
        if abs(surf.a0[1]) > surf.lespcrit[1]
            if (surf.a0[1] >= 0.)
                lesp_cond = surf.lespcrit[1]
            else
                lesp_cond = -surf.lespcrit[1]
            end

            # T3 depends on recenetly shed LEV
            nlev = length(curfield.lev)
            if surf.levflag[1] == 0
                le_vel_x = surf.kinem.u - surf.kinem.alphadot*sin(surf.kinem.alpha)*surf.pvt*surf.c + surf.uind[1]
                le_vel_z = -surf.kinem.alphadot*cos(surf.kinem.alpha)*surf.pvt*surf.c- surf.kinem.hdot + surf.wind[1]
                xloc_lev = surf.bnd_x[1] + 0.5*le_vel_x*dt
                zloc_lev = surf.bnd_z[1] + 0.5*le_vel_z*dt
            else
                xloc_lev = surf.bnd_x[1] + (1. /3.)*(curfield.lev[nlev].x - surf.bnd_x[1])
                zloc_lev = surf.bnd_z[1]+(1. /3.)*(curfield.lev[nlev].z - surf.bnd_z[1])
            end

            for ib = 1:surf.ndiv
                xdist = surf.bnd_x[ib] - xloc_lev
                zdist = surf.bnd_z[ib] - zloc_lev
                distsq = xdist*xdist + zdist*zdist
                T3[ib] = (surf.cam_slope[ib]*zdist + xdist)/(2*pi*sqrt(distsq^2 + vcore^4))
            end
            I3 = simpleTrapz(T3.*(cos.(surf.theta) .- 1. ), surf.theta)
            J3 = -simpleTrapz(T3, surf.theta)/(pi*surf.uref)

            det = J3*(I2 + 1. ) - J2*(I3 + 1. )

            tevstr = (-J3*(I1 + sig_prev) + (I3 + 1)*(J1 - lesp_cond))/det
            levstr = (J2*(I1 + sig_prev) - (I2 + 1)*(J1 - lesp_cond))/det

            #Recalculate required fourier terms
            surf.a0[1] = J1 + J2*tevstr + J3*levstr
            for ia = 1:3
                surf.aterm[ia] = 2. *(simpleTrapz(T1.*cos.(ia*surf.theta), surf.theta) +
                                      tevstr*simpleTrapz(T2.*cos.(ia*surf.theta), surf.theta) +
                                      levstr*simpleTrapz(T3.*cos.(ia*surf.theta), surf.theta))/(pi*surf.uref)
            end

            push!(curfield.tev, TwoDVort(xloc_tev, zloc_tev, tevstr, vcore, 0., 0.))
            push!(curfield.lev, TwoDVort(xloc_lev, zloc_lev, levstr, vcore, 0., 0.))

            for ia = 4:surf.naterm
                surf.aterm[ia] = 2. *(simpleTrapz(T1.*cos.(ia*surf.theta), surf.theta) +
                                      tevstr*simpleTrapz(T2.*cos.(ia*surf.theta), surf.theta) +
                                      levstr*simpleTrapz(T3.*cos.(ia*surf.theta), surf.theta))/(pi*surf.uref)
            end

            surf.levflag[1] = 1
        else
            push!(curfield.tev, TwoDVort(xloc_tev, zloc_tev, tevstr, vcore, 0., 0.))

            for ia = 4:surf.naterm
                surf.aterm[ia] = 2. *(simpleTrapz(T1.*cos.(ia*surf.theta), surf.theta) +
                                      tevstr*simpleTrapz(T2.*cos.(ia*surf.theta), surf.theta))/(pi*surf.uref)
            end

            surf.levflag[1] = 0
        end

#Set previous values of aterm to be used for derivatives in next time step
surf.a0prev[1] = surf.a0[1]
for ia = 1:3
    surf.aprev[ia] = surf.aterm[ia]
end

#Calculate bound vortex strengths
update_bv(surf)

# Delete or merge vortices if required
controlVortCount(delvort, surf.bnd_x[Int(round(surf.ndiv/2))], surf.bnd_z[Int(round(surf.ndiv/2))], curfield)

# free wake rollup
wakeroll(surf, curfield, dt)

# Calculate force and moment coefficients
cl, cd, cm = calc_forces(surf, [curfield.u[1], curfield.w[1]])

# write flow details if required
if writeflag == 1
    if istep in writeArray
        dirname = "$(round(t,sigdigits=nround))"
        writeStamp(dirname, t, surf, curfield)
    end
end

# for writing in resultsSummary
mat = hcat(mat,[t, surf.kinem.alpha, surf.kinem.h, surf.kinem.u, surf.a0[1], cl, cd, cm])

end

mat = mat'

f = open("resultsSummary", "w")
println(f,"#time \t","alpha(deg) \t","h/c \t", "u/uref \t", "A0 \t", "Cl \t", "Cd \t", "Cm ")
DelimitedFiles.writedlm(f, mat)
close(f)

mat, surf, curfield

end

function ldvm2DOF(surf::TwoDSurf, curfield::TwoDFlowField, strpar::TwoDOFPar, kinem::KinemPar2DOF, nsteps::Int64 = 500, dtstar::Float64 = 0.015, startflag = 0, writeflag = 0, writeInterval = 1000., delvort = delNone(); maxwrite = 50, nround=6)

    # If a restart directory is provided, read in the simulation data
    if startflag == 0
        mat = zeros(0, 8)
        t = 0.
    elseif startflag == 1
        dirvec = readdir()
        dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec)
        latestTime = maximum(dirresults)
        mat = DelimitedFiles.readdlm("resultsSummary",header=true)
        t = mat[end,1]
    else
        throw("invalid start flag, should be 0 or 1")
    end
    mat = mat'

    dt = dtstar*surf.c/surf.uref

    # if writeflag is on, determine the timesteps to write at
    if writeflag == 1
        writeArray = Int64[]
        tTot = nsteps*dt
        for i = 1:maxwrite
            tcur = writeInterval*real(i)
            if t > tTot
                break
            else
                push!(writeArray, Int(round(tcur/dt)))
            end
        end
    end

    cl = 0.
    cm = 0.

    #Intialise flowfield

    for istep = 1:nsteps
        #Udpate current time
        t = t + dt

        #Update external flowfield
        update_externalvel(curfield, t)

        #Update kinematic parameters (based on 2DOF response)
        update_kinem2DOF(surf, strpar, kinem, dt, cl, cm)

        #Update bound vortex positions
        update_boundpos(surf, dt)

        #Add a TEV with dummy strength
        place_tev(surf,curfield,dt)

        #Update incduced velocities on airfoil
        update_indbound(surf, curfield)

        #Solve for TEV strength to satisfy Kelvin condition
        kelv = KelvinCondition(surf,curfield)
        soln = nlsolve(not_in_place(kelv), [-0.01])
        curfield.tev[length(curfield.tev)].s = soln.zero[1]

        #Update adot
        #update_a2a3adot(surf,dt)
        update_atermdot(surf, dt)

        lesp = surf.a0[1]

        #Check for LESP condition
        if (abs(lesp)>surf.lespcrit[1])
            #2D iteration if LESP_crit is exceeded
            #Remove the previous tev

            #Add a LEV with dummy strength
            place_lev(surf,curfield,dt)

            #Solve for TEV and LEV strengths to satisfy Kelvin condition and Kutta condition at leading edge
            kelvkutta = KelvinKutta(surf,curfield)
            soln = nlsolve(not_in_place(kelvkutta), [-0.01; 0.01])
            (curfield.tev[length(curfield.tev)].s, curfield.lev[length(curfield.lev)].s) = soln.zero[1], soln.zero[2]

            #set flag for levshedding=on
            surf.levflag[1] = 1
        else
            surf.levflag[1] = 0
        end

        #Update rest of Fourier terms
        update_a2toan(surf)

        #Set previous values of aterm to be used for derivatives in next time step
        surf.a0prev[1] = surf.a0[1]
        for ia = 1:surf.naterm
            surf.aprev[ia] = surf.aterm[ia]
        end

        #Calculate bound vortex strengths
        update_bv(surf)

        # Delete or merge vortices if required
        controlVortCount(delvort, surf.bnd_x[Int(round(surf.ndiv/2))], surf.bnd_z[Int(round(surf.ndiv/2))], curfield)

        # free wake rollup
        wakeroll(surf, curfield, dt)

        # Calculate force and moment coefficients
        cl, cd, cm = calc_forces(surf, [curfield.u[1], curfield.w[1]])

        # write flow details if required
        if writeflag == 1
            if istep in writeArray
                dirname = "$(round(t,sigdigits=nround))"
                writeStamp(dirname, t, surf, curfield)
            end
        end

        # for writing in resultsSummary
        mat = hcat(mat,[t, surf.kinem.alpha, surf.kinem.h, surf.kinem.u, surf.a0[1], cl, cd, cm])

    end

    mat = mat'

    f = open("resultsSummary", "w")
    println(f,"#time \t","alpha(deg) \t","h/c \t", "u/uref \t", "A0 \t", "Cl \t", "Cd \t", "Cm ")
    DelimitedFiles.writedlm(f, mat)
    close(f)

    mat, surf, curfield
end


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ a38a9620-6bea-11eb-2293-c5469f2d5bdc
begin
	import Pkg; Pkg.add("Plotly")
	import Pkg; Pkg.add("PyPlot")
	import Pkg; Pkg.add("PlutoUI")
	import Pkg; Pkg.add("DataFrames")
	using PlutoUI
	using Plots
	using LinearAlgebra
	using SparseArrays
	using DataFrames
end

# ╔═╡ 05f72bc0-6c45-11eb-3602-c955811f9acf

md" Choose the **number of data points** Pascal is using in her statistics 👉 $(@bind n_data_points Slider(500:100:2000, default=1000, show_value = true)) "

# ╔═╡ 70049580-6bea-11eb-390b-d3013dd12235
md"Choose **probability** of the **ocean current** to lead to a deviation 👉 $(@bind ocean_current_prob Slider(0:0.05:0.3, show_value = true)) "

# ╔═╡ deb2ef40-6bea-11eb-3a13-c9642d262b93
md"Choose **probability** of **Pascal to fail on keeping the course** 👉 $(@bind pascal_fail_prob Slider(0:0.05:0.3, show_value = true)) "

# ╔═╡ 08bf49de-6bec-11eb-1b80-5f6e74cbb8d0
@bind direction_of_current Select([ "north" => "North ⬆", "east" => "East ➡", "south" => "South ⬇", "west" => "West ⬅"])

# ╔═╡ 18893b90-6c5c-11eb-051d-95eef7af0f0f

# ╔═╡ 335a84c0-6bf2-11eb-3249-fbfb0210657c
dd = Dict([(1, "north"), (2, "east"), (3, "south"), (4, "west")])

# ╔═╡ 26e7f8d0-6bed-11eb-09d6-a1597aeefc22
direction_of_current == "north"

# ╔═╡ c4582980-6bef-11eb-18e9-0b1732c8d932
direction = rand(1:4, n_data_points)

# ╔═╡ dc84d260-6bef-11eb-1e6c-23d28cf6b1d3
begin
	rr = rand(n_data_points)
	sail_deviation = zeros(1,n_data_points)
	#sail_correct = pascal_fail_prob/2 .< rr .< 1-pascal_fail_prob/2;
	sail_deviation = (-(rr .<= pascal_fail_prob/2)) + (rr.>= 1-pascal_fail_prob/2)

end

# ╔═╡ 71561790-6bf1-11eb-1acf-a1a46a8eb3d8
begin
	ss = rand(n_data_points)
	if direction_of_current == "north"
		left_drift = 2;
		right_drift = 4;
	elseif direction_of_current == "east"
		left_drift = 3 ;
		right_drift = 1;
	elseif direction_of_current == "south"
		left_drift = 4 ;
		right_drift = 2;
	else
		left_drift = 1;
		right_drift = 3;
	end

	added_deviation = (ss .<= ocean_current_prob) .*((-(direction .== left_drift)) + (direction .== right_drift))

end

# ╔═╡ 054eba90-6bf4-11eb-2214-ede88296a7b8
deviation = sail_deviation .+ added_deviation

# ╔═╡ 34da1880-6bf5-11eb-3e34-838a118066b4
table = [[sum(deviation[ direction .==i] .<= -1) for i=1:4] [sum(deviation[ direction .==i] .== 0) for i=1:4] [sum(deviation[ direction .==i] .>= 1) for i=1:4]]

# ╔═╡ 9d0230c0-6bf8-11eb-2445-777f55b92f4a
DataFrame([["North: "; "East:  "; "South: "; "West:  "] table], [:direction, :left_dev, :no_dev, :right_dev])

# ╔═╡ d83840b0-6c54-11eb-31d3-4d1d978c4e4a
table

# ╔═╡ 83cb7cd0-6bf6-11eb-2486-3f0ffc2c8416
(table[:,1] .+ table[:,3])./sum(table, dims=2)

# ╔═╡ 9b0e3050-6bf5-11eb-18df-c9ec17fc42be
sum(table,dims=1)

# ╔═╡ Cell order:
# ╟─a38a9620-6bea-11eb-2293-c5469f2d5bdc
# ╟─05f72bc0-6c45-11eb-3602-c955811f9acf
# ╟─70049580-6bea-11eb-390b-d3013dd12235
# ╟─deb2ef40-6bea-11eb-3a13-c9642d262b93
# ╟─08bf49de-6bec-11eb-1b80-5f6e74cbb8d0
# ╟─9d0230c0-6bf8-11eb-2445-777f55b92f4a
# ╠═18893b90-6c5c-11eb-051d-95eef7af0f0f
# ╠═d83840b0-6c54-11eb-31d3-4d1d978c4e4a
# ╠═83cb7cd0-6bf6-11eb-2486-3f0ffc2c8416
# ╟─335a84c0-6bf2-11eb-3249-fbfb0210657c
# ╠═26e7f8d0-6bed-11eb-09d6-a1597aeefc22
# ╠═c4582980-6bef-11eb-18e9-0b1732c8d932
# ╠═dc84d260-6bef-11eb-1e6c-23d28cf6b1d3
# ╠═71561790-6bf1-11eb-1acf-a1a46a8eb3d8
# ╠═054eba90-6bf4-11eb-2214-ede88296a7b8
# ╠═34da1880-6bf5-11eb-3e34-838a118066b4
# ╠═9b0e3050-6bf5-11eb-18df-c9ec17fc42be


text:
#-----------------------------------------------------------------------
module seval     # EVALUATION MODULE
#-----------------------------------------------------------------------

using LinearAlgebra , Statistics, Combinatorics, ForwardDiff
using DoubleFloats

include("type.jl")
include("metrics/WeakFieldIso.jl")
include("metrics/Minkowski.jl")
include("srl/FHC22.jl")
include("srl/RTC21.jl")
include("srl/mloc.jl")
include("outlier.jl")
include("geosol.jl")

Δψ0 = (π/180)*10    # Default angle margin

#-----------------------------------------------------------------------
#       STRUCTS
#-----------------------------------------------------------------------

"""
The TestCases datatype may be populated by the associated function of
the form:

    TestCases( par::NTuple , N::Int , ne::Int , X::Array{RealMtx,1} , Xtar::Array{RealVec,1} )

The variables are defined in the following way:

    par     ::  NTuple              # Parameter tuple
    N       ::  Int                 # Number of generated test cases
    ne      ::  Int                 # Number of emission points
    X       ::  Array{RealMtx,1}    # Emission points
    Xtar    ::  Array{RealVec,1}    # Target point

"""
struct TestCases                                    # TestCases datatype
    par     ::  NTuple              # Parameter tuple
    N       ::  Int                 # Number of cases
    ne      ::  Int                 # Number of emission points
    X       ::  Array{RealMtx,1}    # Emission points
    Xtar    ::  Array{RealVec,1}    # Target point
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
The TestData datatype may be populated by the associated function of the
form:

    TestData(   par::NTuple , N::Int    ,
                X::Array{RealMtx,1}     ,   Xtar::Array{RealVec,1}     ,
                Xc::Array{RealVec,1}    ,   Xsc::Array{RealVec,1}      ,
                erh::RealVec    ,   erv::RealVec    ,   err::RealVec   ,
                erhC::RealVec   ,   ervC::RealVec   ,   errC::RealVec  ,
                Xc2::Array{RealVec,1}   ,   Xsc2::Array{RealVec,1}     ,
                erh2::RealVec   ,   erv2::RealVec   ,   err2::RealVec  ,
                erhC2::RealVec  ,   ervC2::RealVec  ,   errC2::RealVec )

The quantities are defined in the following way:

    par     ::  NTuple              # Parameter tuple
    N       ::  Int                 # Number of cases
    X       ::  Array{RealMtx,1}    # Emission points
    Xtar    ::  Array{RealVec,1}    # Target point
    Xc      ::  Array{RealVec,1}    # Flat spacetime position
    Xsc     ::  Array{RealVec,1}    # Squirrel position
    erh     ::  RealVec             # Horizontal error relative to Xsc
    erv     ::  RealVec             # Vertical error relative to Xsc
    err     ::  RealVec             # Total error relative to Xsc
    erhC    ::  RealVec             # Horizontal error relative to Xc
    ervC    ::  RealVec             # Vertical error relative to Xc
    errC    ::  RealVec             # Total error relative to Xc
    Xc2     ::  Array{RealVec,1}    # Flat spacetime position     (Aux.)
    Xsc2    ::  Array{RealVec,1}    # Squirrel position      (Auxiliary)
    erh2    ::  RealVec             # Horizontal errors      (Auxiliary)
    erv2    ::  RealVec             # Vertical error         (Auxiliary)
    err2    ::  RealVec             # Total error            (Auxiliary)
    erhC2   ::  RealVec             # Horizontal errors      (Auxiliary)
    ervC2   ::  RealVec             # Vertical error         (Auxiliary)
    errC2   ::  RealVec             # Total error            (Auxiliary)

The quantities suffixed with `2` are auxiliary quantities which are
needed in the four-point case, in which the special relativistic
location algorithms generally suffer from the bifurcation problem.

"""
struct TestData                                      # TestData datatype
    par     ::  NTuple              # Parameter tuple
    N       ::  Int                 # Number of cases
    X       ::  Array{RealMtx,1}    # Emission points
    Xtar    ::  Array{RealVec,1}    # Target point
    Xc      ::  Array{RealVec,1}    # Flat spacetime position
    Xsc     ::  Array{RealVec,1}    # Squirrel position
    erh     ::  RealVec             # Horizontal errors
    erv     ::  RealVec             # Vertical error
    err     ::  RealVec             # Total error
    erhC    ::  RealVec             # Horizontal errors
    ervC    ::  RealVec             # Vertical error
    errC    ::  RealVec             # Total error
    Xc2     ::  Array{RealVec,1}    # Flat spacetime position     (Aux.)
    Xsc2    ::  Array{RealVec,1}    # Squirrel position      (Auxiliary)
    erh2    ::  RealVec             # Horizontal errors      (Auxiliary)
    erv2    ::  RealVec             # Vertical error         (Auxiliary)
    err2    ::  RealVec             # Total error            (Auxiliary)
    erhC2   ::  RealVec             # Horizontal errors      (Auxiliary)
    ervC2   ::  RealVec             # Vertical error         (Auxiliary)
    errC2   ::  RealVec             # Total error            (Auxiliary)
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    tcfl( tc::TestCases , tpfl::DataType=Float64 )

The `tcfl` function changes the floating point type of `tc` to `tpfl`.

"""
function tcfl( tc::TestCases , tpfl::DataType=Float64 )
    par2    = tpfl.(tc.par)
    N       = tc.N
    ne      = tc.ne
    X2      = [zeros(tpfl,4,ne) for _ in 1:N]
    Xtar2   = [zeros(tpfl,4) for _ in 1:N]

    for i=1:N
        X2[i]     = tpfl.(tc.X[i])
        Xtar2[i]  = tpfl.(tc.Xtar[i])
    end

    return TestCases( par2 , N , ne , X2 , Xtar2 )
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    tc2tup( tc::TestCases )

The `tc2tup` function changes `tc` from an object of type `TestCases` to
a tuple.

"""
function tc2tup( tc::TestCases )
    return ( tc.par , tc.N , tc.ne , tc.X , tc.Xtar )
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    td2tup( td::TestData )

The `td2tup` function changes `td` from an object of type `TestData` to
a tuple.

"""
function td2tup( td::TestData )
    return ( td.par , td.N , td.X , td.Xtar , td.Xc , td.Xsc , td.erh ,
             td.erv , td.err , td.erhC , td.ervC , td.errC , td.Xc2 ,
             td.Xsc2 , td.erh2 , td.erv2 , td.err2 , td.erhC2 ,
             td.ervC2 , td.errC2 )
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    tup2tc( tctup::Tuple )

The `tc2tup` function changes `tctup` from a Tuple to an object of type
`TestCases`.

"""
function tup2tc( tctup::Tuple )
    return TestCases(tctup[1],tctup[2],tctup[3],tctup[4],tctup[5])
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    tup2td( tdtup::Tuple )

The `td2tup` function changes `tdtuple` from a Tuple to an object of
type `TestData`.

"""
function tup2td( tdtup::Tuple )
    return TestData(tdtup[1] ,tdtup[2] ,tdtup[3] ,tdtup[4] ,tdtup[5] ,
                    tdtup[6] ,tdtup[7] ,tdtup[8] ,tdtup[9] ,tdtup[10],
                    tdtup[11],tdtup[12],tdtup[13],tdtup[14],tdtup[15],
                    tdtup[16],tdtup[17],tdtup[18],tdtup[19],tdtup[20])
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#
#       GENERAL FUNCTIONS
#
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    roteuler( tpfl::DataType=Float64 )

The `roteuler` function constructs a random rotation matrix `R` and its
inverse as a tuple `(R,inv(R))`.

"""
function roteuler( tpfl::DataType=Float64 )
    α = 2*pi*rand(tpfl)     # Euler angles
    β = 2*pi*rand(tpfl)
    γ = 2*pi*rand(tpfl)

    ( c1 , s1 ) = ( cos(α) , sin(α) )
    ( c2 , s2 ) = ( cos(β) , sin(β) )
    ( c3 , s3 ) = ( cos(γ) , sin(γ) )

    R = [     c2        -c3*s2           s2*s3          ;
            c1*s2   c1*c2*c3-s1*s3  -c3*s1-c1*c2*s3     ;
            s1*s2   c1*s3+c2*c3*s1   c1*c3-c2*s1*s3     ]

    return (R,inv(R))
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    angvec( u::RealVec , v::RealVec )

The `angvec` function computes the angle between the vectors `u` and
`v`.

"""
function angvec( u::RealVec , v::RealVec )
    return acos( dot(u,v)/(norm(u)*norm(v)) )
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    angmar( θ::Real , Δψ::Real )

The `angmar` function checks whether the angle ``θ``` falls in the range
``0<θ<π/2 - |Δψ|``.

"""
function angmar( θ::Real , Δψ::Real )
    if θ <= π/2 && θ >+ 0
        ϕ = abs( π/2 - θ )
        if ϕ > Δψ
            return true
        else
            return false
        end
    else
        return false
    end
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    slchk( X::RealMtx )

The `slchk` function checks whether the emission points in the matrix
``X`` are spacelike separated with respect to the Minkowski metric
``η_{μν}``.

"""
function slchk( X::RealMtx )
    tpfl = typeof(X[1])
    ne = size(X)[2]

    if ne == 1
        return true
    elseif ne > 1
        Z = [zeros(tpfl,4) for _ in 1:ne]

        for i=1:ne
            Z[i] = X[:,i]
        end

        W = collect(combinations(Z,2))      # Combinations
        nc = length(W)
        Y = zeros(tpfl,nc)

        b = true
        for i=1:nc
            Y[i] = ηdot( W[i][1] - W[i][2] , W[i][1] - W[i][2] )
            if Y[i] > zero(tpfl)
                b = b * true
            else
                b = b * false
            end
        end

        return b
    end
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#
#       POINT GENERATION
#
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
#       Generates single random 3-component directional vector
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
"""
    vrgen( scale::Real , Δψ::Real , u::RealVec=Float64[1;0;0] )

The `vrgen` function generates a vector of length `scale` which makes
an angle ``θ<π/2-|ΔΨ|``.

"""
function vrgen( scale::Real , Δψ::Real , u::RealVec=Float64[1;0;0] )
    # scale determines length, Δψ determines angle margins
    tpfl=typeof(scale)

    b = false

    v = zeros(tpfl,3)

    while b == false
        vx = 2*(rand(tpfl)-1/2)
        vy = 2*(rand(tpfl)-1/2)
        vz = 2*(rand(tpfl)-1/2)
        vl = sqrt( vx^2 + vy^2 + vz^2 )
        if vl>1
            b = false
        else
            b = true
        end
        v = tpfl[ vx ; vy ; vz ]/vl
        b = angmar(angvec(u,v),Δψ)
    end
    return scale*v
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    λiRscalc( xtar::RealVec , V::RealVec , Re::Real , Rs::Real )

The `λiRscalc` function computes the value of the affine parameter `λ`
for the test cases. Given the straight line `Y=Re*xtar+λ*V` passing
through the point `Re*xtar`, it computes the value of `λ` at which the
spatial components of the line `Y[2:4]` have a norm `norm(Y[2:4])=Rs`.

"""
function λiRscalc( xtar::RealVec , V::RealVec , Re::Real , Rs::Real )
return  (
          -(Re*(V[1]*xtar[1] + V[2]*xtar[2] + V[3]*xtar[3])) +
          sqrt(abs(
            4*Rs^2*(V[1]^2 + V[2]^2 + V[3]^2) +
            Re^2*(
                  4*(V[1]*xtar[1] + V[2]*xtar[2] + V[3]*xtar[3])^2 -
                  4*(V[1]^2 + V[2]^2 + V[3]^2)*(xtar[1]^2 +
                  xtar[2]^2 + xtar[3]^2)
                 )
            ))/2
        )/(V[1]^2 + V[2]^2 + V[3]^2)
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    x3gen(tpfl::DataType=Float64)

The `x3gen` function generates a random 3-component position vector of
length `rell(x)` is the `WGS84` Earth reference ellipsoid radius.

"""
function x3gen(tpfl::DataType=Float64)
    x   = tpfl[ 1 ; 0 ; 0 ]
    rot = roteuler(tpfl)
    x = rot[1]*x
    return x*rell(x)
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    tidc( x::RealVec , v::RealVec , λ::Real , gfunc::Function )

The `tidc` function generates initial data for null geodesics given a
position vector `x`, velocity vector `v`, affine parameter `λ`, and a
metric function `gfunc`.

"""
function tidc( x::RealVec , v::RealVec , λ::Real , gfunc::Function )
    tpfl=typeof(x[1])

    X = zeros(tpfl,4)
    V = zeros(tpfl,4)

    X[2:4] = x
    V[2:4] = v

    V[1] = -λ

    V = nullenforcerp( V , X , gfunc )

    return vcat(X,V)
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#
#       PROBLEM GENERATOR
#
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    pgen( Rs::Real , gfunc::Function , tol::Real , ne::Int=4 ,
          Δψ::Real=Δψ0 )

The `pgen` function generates a target point `Xtar` and `ne` emission
points in a `4×ne` matrix `X` such that the points in `X` lie on the
past light cone of `Xtar` with respect to the metric `gfunc`, and have
spatial radii values of `~Rs` (defined with respect to spatial origin).
The parameter `tol` is the tolerance parameter for the integration. This
function returns the tuple `(X,Xtar)`.

"""
function pgen( Rs::Real , gfunc::Function , tol::Real , ne::Int=4 ,
               Δψ::Real=Δψ0 )
    tpfl=typeof(Rs)

    Z   = zeros(tpfl,8,ne)
    Zi  = zeros(tpfl,8,ne)

    X   = zeros(tpfl,4,ne)

    Xc  = zeros(tpfl,4,4)

    xi  = x3gen(tpfl)
    vi  = zeros(tpfl,3,ne)
    λi  = zeros(tpfl,ne)

    for i=1:ne
        B = true
        while B
            vi[:,i] = vrgen( one(tpfl) , Δψ , xi )
            λi[i]   = λiRscalc(xi/rell(xi),vi[:,i],rell(xi),Rs)
            vi[:,i] = λi[i]*vi[:,i]

            Zi[:,i] = tidc( xi , vi[:,i] , λi[i] , gfunc )
            Z[:,i]  = solveZ( Zi[:,i] , gfunc , tol , tol
                              , AutoVern9(Rodas5()) )
            X[:,i]  = Z[1:4,i]

            # Check spacelike separation
            if i==1
                B = false
            elseif i>1 && i<4
                if slchk( X[:,1:i] )
                    B = false
                end
            elseif i>=4
                if slchk(X[:,1:i])
                    B = false
                end
            end
        end
    end

    return ( X , Zi[1:4,1] )
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#
#       TEST CASE GENERATOR
#
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    gen( N::Int , g::Function , ne::Int=4 , Δψ::Real=Δψ0
              , tpfl::DataType=Float64
              , REs ::  Real=1.4365276211950395e9       # Earth radius
              , RR  ::  Real=1.4365277e9    # Just above Earth surface
              , Rs  ::  Real=6e9            # Satellite radius
              , tolh::  Real=1e-14          # Tolerance for ODE solvers
              , tol ::  Real=1e-10          # Tolerance for ODE solvers
            )

The `gen` function generates `N` sets of target points `Xtar` and
emission points `X` with respect to the metric `gfunc`. This function
returns a quantity of the TestCases datatype.

"""
function gen( N::Int , g::Function , ne::Int=4 , Δψ::Real=Δψ0
              , tpfl::DataType=Float64
              , REs ::  Real=1.4365276211950395e9       # Earth radius
              , RR  ::  Real=1.4365277e9    # Just above Earth surface
              , Rs  ::  Real=6e9            # Satellite radius
              , tolh::  Real=1e-14          # Tolerance for ODE solvers
              , tol ::  Real=1e-10          # Tolerance for ODE solvers
            )

par = (REs,RR,Rs,tolh,tol)

tc = TestCases(  par , N , ne , [zeros(tpfl,4,ne) for _ in 1:N]
                        , [zeros(tpfl,4) for _ in 1:N] )

for i=1:N
    ZZ          = pgen( Rs , g , tolh , ne , Δψ )
    tc.X[i]     = ZZ[1]
    tc.Xtar[i]  = ZZ[2]
    print("\r$i")
end

return tc
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
#
#       EVALUATION
#
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    erH( Xtar::RealVec , ΔX::RealVec )

The `erH` function projects the vector `ΔX` in the direction orthogonal
to the vector `Xtar`. This is used to compute the horizontal error in
the terrestrial positioning problem.

"""
function erH( Xtar::RealVec , ΔX::RealVec )
    tpfl = typeof(Xtar[1])
    x    = copy(Xtar[2:4])
    Δx   = copy(ΔX[2:4])
    rv   = x/norm(x)
    δ    = one(tpfl)*(I(3))
    P    = δ - rv*transpose(rv)
    ΔE   = P*Δx

    return  norm(ΔE)
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    erV( Xtar::RealVec , ΔX::RealVec )

The `erV` function projects the vector `ΔX` in the direction along the
vector `Xtar`. This is used to compute the vertical error in the
terrestrial positioning problem.

"""
function erV( Xtar::RealVec , ΔX::RealVec )
    tpfl = typeof(Xtar[1])
    x    = copy(Xtar[2:4])
    Δx   = copy(ΔX[2:4])
    rv   = x/norm(x)

    return  dot(rv,Δx)
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
    main(  tc::TestCases , sloc::Function , g::Function , Nx::Int=-1 ,
           tpfl::DataType=Double64 , tol::Real=1e-10 , ξ::Real=2e1 ,
           nb::Int=24 , ne::Real=6 )

This is the main evaluation function. It takes as input a collection of
test cases contained in a quantity of datatype `TestCases`, and
evaluates the squirrel locator function `sloc` on each test case
relative to a given metric function `g`. The function `main` returns
results in a quantity of type `TestData`.

The variable `Nx` is a limiter for the number of test cases to run. The
variable `tpfl` specifies the floating point precision for the
calculations. variables `tol`, `ξ`, `nb`,  and `ne` are inputs for the
`sloc` function, and correspond respectively to the integration
tolerance, outlier detection threshold, number of steps for the Broyden
algorithm, and number of emission points.

"""
function main(  tc::TestCases , sloc::Function , g::Function ,
                Nx::Int=-1 , tpfl::DataType=Float64 , tol::Real=1e-10 ,
                ξ::Real=2e1 , nb::Int=24 , ne::Real=6 )
tc2     = tcfl(tc,tpfl)
par     = tc2.par
REs     = par[1]
RR      = par[2]
Rs      = par[3]
tolh    = par[4]

if Nx==-1
N   = tc2.N
else
N   = Nx
end
np      = tc2.ne

nc      = binomial(np,4)

X4      = zeros(tpfl,4,4)
ΔXsq    = zeros(tpfl,4)
ΔXcr    = zeros(tpfl,4)

ΔXs2    = [[zeros(tpfl,4)];[zeros(tpfl,4)]]
ΔXc2    = [[zeros(tpfl,4)];[zeros(tpfl,4)]]

d = size(tc2.X[1])[2]

td = TestData( par , N , tc2.X , tc2.Xtar
         , [zeros(tpfl,4) for _ in 1:N]
         , [zeros(tpfl,4) for _ in 1:N]
         , zeros(tpfl,N) , zeros(tpfl,N) , zeros(tpfl,N)
         , zeros(tpfl,N) , zeros(tpfl,N) , zeros(tpfl,N)
         , [zeros(tpfl,4) for _ in 1:N]
         , [zeros(tpfl,4) for _ in 1:N]
         , zeros(tpfl,N) , zeros(tpfl,N) , zeros(tpfl,N)
         , zeros(tpfl,N) , zeros(tpfl,N) , zeros(tpfl,N)
         )

if d >= 5 && ne >= 5
    for i=1:N
        td.Xc[i]    = mlocator( td.X[i] )
        td.Xsc[i]   = sloc( td.X[i] , g , tol , ne , nb , ξ , Double64 )
        ΔXsq = td.Xtar[i] - td.Xsc[i]
        ΔXcr = td.Xtar[i] - td.Xc[i]

        td.erh[i]   = abs( erH( td.Xtar[i] , ΔXsq ) )
        td.erv[i]   = abs( erV( td.Xtar[i] , ΔXsq ) )
        td.err[i]   = norm( ΔXsq )

        td.erhC[i]  = abs( erH( td.Xtar[i] , ΔXcr ) )
        td.ervC[i]  = abs( erV( td.Xtar[i] , ΔXcr ) )
        td.errC[i]  = norm( ΔXcr )
        print("\r$i")
    end
elseif d == 4 || (d >= 4 && ne == 4)
    for i=1:N
        (td.Xc[i] ,td.Xc2[i])   = locator4FHC22( td.X[i] )
        (td.Xsc[i],td.Xsc2[i])  = sloc( td.X[i] , g , tol , ne , nb ,
                                        ξ , Double64 )
        ΔXsq  = td.Xtar[i] - td.Xsc[i]
        ΔXcr  = td.Xtar[i] - td.Xc[i]
        ΔXsq2 = td.Xtar[i] - td.Xsc2[i]
        ΔXcr2 = td.Xtar[i] - td.Xc2[i]

        td.erh[i]   = abs( erH( td.Xtar[i] , ΔXsq ) )
        td.erv[i]   = abs( erV( td.Xtar[i] , ΔXsq ) )
        td.err[i]   = norm( ΔXsq )

        td.erh2[i]   = abs( erH( td.Xtar[i] , ΔXsq2 ) )
        td.erv2[i]   = abs( erV( td.Xtar[i] , ΔXsq2 ) )
        td.er2r[i]   = norm( ΔXsq2 )

        td.erhC[i]  = abs( erH( td.Xtar[i] , ΔXcr ) )
        td.ervC[i]  = abs( erV( td.Xtar[i] , ΔXcr ) )
        td.errC[i]  = norm( ΔXcr )

        td.erhC2[i]  = abs( erH( td.Xtar[i] , ΔXcr2 ) )
        td.ervC2[i]  = abs( erV( td.Xtar[i] , ΔXcr2 ) )
        td.errC2[i]  = norm( ΔXcr2 )
        print("\r$i")
    end
end

return td
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
end     # END EVALUATION MODULE
#-----------------------------------------------------------------------


text:
# ---
# title: 1089. Duplicate Zeros
# id: problem1089
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/duplicate-zeros/description/>
# hidden: true
# ---
#
# Given a fixed length array `arr` of integers, duplicate each occurrence of
# zero, shifting the remaining elements to the right.
#
# Note that elements beyond the length of the original array are not written.
#
# Do the above modifications to the input array **in place** , do not return
# anything from your function.
#
#
#
# **Example 1:**
#
#
#
#     Input: [1,0,2,3,0,4,5,0]
#     Output: null
#     Explanation: After calling your function, the **input** array is modified to: [1,0,0,2,3,0,0,4]
#
#
# **Example 2:**
#
#
#
#     Input: [1,2,3]
#     Output: null
#     Explanation: After calling your function, the **input** array is modified to: [1,2,3]
#
#
#
#
# **Note:**
#
#   1. `1 <= arr.length <= 10000`
#   2. `0 <= arr[i] <= 9`
#
#
## @lc code=start
using LeetCode

function duplicate_zeros!(arr::Vector{Int})
    len = length(arr)
    for i in reverse(eachindex(arr))
        iszero(arr[i]) && insert!(arr, i, 0)
    end
    resize!(arr, len)
    nothing
end

## add your code here:
## @lc code=end


text:
function metapop()
  # make metapopulation
  mp = ZTMetaPopulation{LanguageDynamics.MomentumSelector}(7, 2.5)

  # add speakers
  for i = 1:200
    inject!(LanguageDynamics.MomentumSelector(0.01, 0.02, 1.0, 4, 0.1), mp)
  end

  # Travel + rendezvous
  while mp.nactions < 0.5*200*100#_000
    travel!(mp)
    rendezvous!(mp)
  end
end


text:
mutable struct BlackBoxOptimizer <: HeuristicOptimizer
    lossfn
    optparams::AbstractDict
    optres
    BlackBoxOptimizer(lossfn::Function; kws...) = new(lossfn, kws, nothing)
end

BlackBoxOptimizer(m::InterpretableModel; lossfn = l2loss, kws...) =
    BlackBoxOptimizer(lossfn; kws...)

const DefaultBBOParams = Dict(
    :PopulationSize => 500,
    :MaxTime => 15.0,
)

function getp(key, d1, opt, dbackup = DefaultBBOParams)
    haskey(d1, key) && return d1[key]
    haskey(opt.optparams, key) && return opt.optparams[key]
    haskey(dbackup, key) && return dbackup[key]
    error("No value for parameter $key found!")
end

function fit!(m::InterpretableModel, o::BlackBoxOptimizer,
                X::AbstractMatrix{Float64}, y::AbstractVector{Float64}; kws...)

    @assert size(X, 1) >= length(y) - maxlag(m)

    fitnessfn(ps::Vector{Float64}) = begin
        yhat, lossvalue = predictandcalcloss(m, X, y, ps, o.lossfn)
        lossvalue # Loss value is the fitness so should be minimixed
    end

    o.optres = bboptimize(fitnessfn;
                SearchSpace = parambounds(m),
                PopulationSize = getp(:PopulationSize, kws, o),
                MaxTime = getp(:MaxTime, kws, o))

    setparams!(m, best_candidate(o.optres))
    m
end


text:
"""
    ZeroDGrid()

Constructs a placeholder grid object for "0D" problems (in other words, systems of ODEs).
"""
struct ZeroDGrid <: AbstractGrid end

"""
    OneDGrid(nx, Lx; x0=-Lx/2, nthreads=Sys.CPU_THREADS, effort=FFTW.MEASURE)

Constrcut a OneDGrid object with size `Lx`, resolution `nx`, and leftmost
position `x0`. FFT plans are generated for `nthreads` CPUs using
FFTW flag `effort`.
"""
struct OneDGrid{T} <: AbstractOneDGrid
  nx::Int
  nk::Int
  nkr::Int
  Lx::T
  dx::T
  x::Array{T,1}
  k::Array{T,1}
  kr::Array{T,1}
  invksq::Array{T,1}
  invkrsq::Array{T,1}

  fftplan::FFTW.cFFTWPlan{Complex{T},-1,false,1}
  rfftplan::FFTW.rFFTWPlan{T,-1,false,1}

  # Range objects that access the aliased part of the wavenumber range
  kalias::UnitRange{Int}
  kralias::UnitRange{Int}
end

function OneDGrid(nx, Lx; x0=-0.5*Lx, nthreads=Sys.CPU_THREADS, effort=FFTW.MEASURE)
  T = typeof(Lx)
  dx = Lx/nx
  x = Array{T}(range(x0, step=dx, length=nx))

  nk = nx
  nkr = Int(nx/2+1)

  i1 = 0:Int(nx/2)
  i2 = Int(-nx/2+1):-1
  k = Array{T}(2π/Lx*cat(i1, i2; dims = 1))
  kr = Array{T}(2π/Lx*cat(i1; dims = 1))

  invksq = @. 1/k^2
  invksq[1] = 0
  invkrsq = @. 1/kr^2
  invkrsq[1] = 0

  FFTW.set_num_threads(nthreads)
  fftplan   = plan_fft(Array{Complex{T},1}(undef, nx); flags=effort)
  rfftplan  = plan_rfft(Array{T,1}(undef, nx); flags=effort)

  # Index endpoints for aliased i, j wavenumbers
  iaL, iaR = Int(floor(nk/3))+1, 2*Int(ceil(nk/3))-1
  kalias  = iaL:iaR
  kralias = iaL:nkr

  OneDGrid(nx, nk, nkr, Lx, dx, x, k, kr, invksq, invkrsq,
           fftplan, rfftplan, kalias, kralias)
end

"""
    TwoDGrid(nx, Lx)
    TwoDGrid(nx, Lx, ny, Ly; x0=-Lx/2, y0=-Ly/2, nthreads=Sys.CPU_THREADS, effort=FFTW.MEASURE)

Constrcut a TwoDGrid object. The two-dimensional domain has size (Lx, Ly),
resolution (nx, ny) and bottom left corner at (x0, y0). FFT plans are generated
which use nthreads threads with the specified planning effort.
"""
struct TwoDGrid{T} <: AbstractTwoDGrid
  nx::Int
  ny::Int
  nk::Int
  nl::Int
  nkr::Int
  Lx::T
  Ly::T
  dx::T
  dy::T
  x::Array{T,2}
  y::Array{T,2}
  X::Array{T,2}
  Y::Array{T,2}
  k::Array{T,2}
  l::Array{T,2}
  kr::Array{T,2}
  K::Array{T,2}
  L::Array{T,2}
  Kr::Array{T,2}
  Lr::Array{T,2}
  KKsq::Array{T,2}      # K^2 + L^2
  invKKsq::Array{T,2}   # 1/KKsq, invKKsq[1, 1]=0
  KKrsq::Array{T,2}     # Kr^2 + Lr^2
  invKKrsq::Array{T,2}  # 1/KKrsq, invKKrsq[1, 1]=0

  fftplan::FFTW.cFFTWPlan{Complex{T},-1,false,2}
  rfftplan::FFTW.rFFTWPlan{T,-1,false,2}

  # Range objects that access the aliased part of the wavenumber range
  kalias::UnitRange{Int}
  kralias::UnitRange{Int}
  lalias::UnitRange{Int}
end

function TwoDGrid(nx, Lx, ny=nx, Ly=Lx; x0=-0.5*Lx, y0=-0.5*Ly, nthreads=Sys.CPU_THREADS, effort=FFTW.MEASURE,
                  T=Float64)
  dx = Lx/nx
  dy = Ly/ny
  nk = nx
  nl = ny
  nkr = Int(nx/2+1)

  # Physical grid
  x = Array{T}(reshape(range(x0, step=dx, length=nx), (nx, 1)))
  y = Array{T}(reshape(range(y0, step=dy, length=ny), (1, ny)))
  X = [ x[i] for i = 1:nx, j = 1:ny]
  Y = [ y[j] for i = 1:nx, j = 1:ny]

  # Wavenubmer grid
  i1 = 0:Int(nx/2)
  i2 = Int(-nx/2+1):-1
  j1 = 0:Int(ny/2)
  j2 = Int(-ny/2+1):-1

  k  = reshape(2π/Lx*cat(i1, i2, dims=1), (nk, 1))
  kr = reshape(2π/Lx*cat(i1, dims=1), (nkr, 1))
  l  = reshape(2π/Ly*cat(j1, j2, dims=1), (1, nl))

  K = [ k[i] for i = 1:nk, j = 1:nl]
  L = [ l[j] for i = 1:nk, j = 1:nl]
  Kr = [ kr[i] for i = 1:nkr, j = 1:nl]
  Lr = [ l[j]  for i = 1:nkr, j = 1:nl]

  KKsq  = @. k^2 + l^2
  invKKsq = 1 ./ KKsq
  invKKsq[1, 1] = 0

  KKrsq = @. kr^2 + l^2
  invKKrsq = 1 ./ KKrsq
  invKKrsq[1, 1] = 0

  # FFT plans
  FFTW.set_num_threads(nthreads)
  fftplan   = plan_fft(Array{Complex{T},2}(undef, nx, ny); flags=effort)
  rfftplan  = plan_rfft(Array{T,2}(undef, nx, ny); flags=effort)

  # Index endpoints for aliased i, j wavenumbers
  iaL, iaR = Int(floor(nk/3))+1, 2*Int(ceil(nk/3))-1
  jaL, jaR = Int(floor(nl/3))+1, 2*Int(ceil(nl/3))-1

  kalias  = iaL:iaR
  kralias = iaL:nkr
  lalias  = jaL:jaR

  TwoDGrid(nx, ny, nk, nl, nkr, Lx, Ly, dx, dy, x, y, X, Y,
           k, l, kr, K, L, Kr, Lr, KKsq, invKKsq, KKrsq, invKKrsq,
           fftplan, rfftplan, kalias, kralias, lalias)
end

"""
    gridpoints(g)

Returns the collocation points of the grid `g` in 2D arrays `X, Y`.
"""
function gridpoints(g)
  X = [ g.x[i] for i=1:g.nx, j=1:g.ny]
  Y = [ g.y[j] for i=1:g.nx, j=1:g.ny]
  X, Y
end

"""
    dealias!(a, g, kalias)

Dealias `a` on the grid `g` with aliased x-wavenumbers `kalias`.
"""
function dealias!(a, g::OneDGrid)
  kalias = size(a, 1) == g.nkr ? g.kralias : g.kalias
  dealias!(a, g, kalias)
  nothing
end

function dealias!(a, g::OneDGrid, kalias)
  @views @. a[kalias, :] = 0
  nothing
end

function dealias!(a, g::TwoDGrid)
  kalias = size(a, 1) == g.nkr ? g.kralias : g.kalias
  dealias!(a, g, kalias)
  nothing
end

function dealias!(a, g::TwoDGrid, kalias)
  @views @. a[kalias, g.lalias, :] = 0
  nothing
end

"""
    makefilter(K; order=4, innerK=0.65, outerK=1)

Returns a filter acting on the non-dimensional wavenumber K that decays exponentially
for K>innerK, thus removing high-wavenumber content from a spectrum it is multiplied with.
The decay rate is determined by order and outerK determines the outer wavenumber at which
the filter is smaller than machine precision.
"""
function makefilter(K; order=4, innerK=0.65, outerK=1)
  TK = typeof(K)
  K = Array(K)
  decay = 15*log(10) / (outerK-innerK)^order # decay rate for filtering function
  filt = @. exp( -decay*(K-innerK)^order )
  filt[real.(K) .< innerK] .= 1
  TK(filt)
end

function makefilter(g::AbstractTwoDGrid; realvars=true, kwargs...)
  K = realvars ?
      @.(sqrt((g.kr*g.dx/π)^2 + (g.l*g.dy/π)^2)) : @.(sqrt((g.k*g.dx/π)^2 + (g.l*g.dy/π)^2))
  makefilter(K; kwargs...)
end

function makefilter(g::AbstractOneDGrid; realvars=true, kwargs...)
  K = realvars ? g.kr*g.dx/π : @.(abs(g.k*g.dx/π))
  makefilter(K; kwargs...)
end

makefilter(g, T, sz; kwargs...) = ones(T, sz).*makefilter(g; realvars=sz[1]==g.nkr, kwargs...)


text:
"""
    ExpandTable(scale, freq) => Array{Int64, 1}
    expandtable(scale, frea) => Array{Int64, 1}

Expand raw score array from freqency table.

# Arguments

- `scale` Vector of Real which stands for raw score scale.
- `freq` Vector of Int which contains counted frequency.

# Example

```julia
sc = [0, 1, 2, 3]
fq = [10, 3, 15, 4]
ExpandTable(sc, fq)
```
"""
function ExpandTable(scale, freq)
    fill.(scale, freq) |> Iterators.flatten |> collect
end

expandtable = ExpandTable


text:
using Flase, Test

n_sheeps = 10
sheeps = DenseSheeps( 10, n_sheeps = n_sheeps )
@testset "Diffusion" begin
    @test sheeps.grid[Flase.getRandomSheep( sheeps )] != 0
    for (n, coords) in sheeps
        @test Flase.getNSheep( sheeps, Tuple(coords)... ) == n
    end # for
    (n, coords) = first(sheeps)
    @test Flase.setNSheep!( sheeps, Tuple(coords)..., sheeps.capacity + 1 ) == nothing
    @test sheeps.grid[coords] == sheeps.capacity + 1
    @test sheeps.diffusion_candidates == [Tuple(coords)]
    @test Flase.kickSheep!( sheeps ) == nothing
    @test sheeps.diffusion_candidates |> isempty == true
    @test sheeps.current_sheep[] == n_sheeps + sheeps.capacity - n + 1
end # testset