text:
#####
##### Sampling: high-level interface and building blocks
#####

"Significant digits to display for reporting."
const REPORT_SIGDIGITS = 3

####
#### parts unaffected by warmup
####

"""
$(TYPEDEF)

A log density bundled with an RNG and options for sampling. Contains the parts of the
problem which are not changed during warmup.

# Fields

$(FIELDS)
"""
struct SamplingLogDensity{R,L,O,S}
    "Random number generator."
    rng::R
    "Log density."
    ℓ::L
    """
    Algorithm used for sampling, also contains the relevant parameters that are not affected
    by adaptation. See eg [`NUTS`](@ref).
    """
    algorithm::O
    "Reporting warmup information and chain progress."
    reporter::S
end

####
#### warmup building blocks
####

###
### warmup state
###

"""
$(TYPEDEF)

Representation of a warmup state. Not part of the API.

# Fields

$(FIELDS)
"""
struct WarmupState{TQ <: EvaluatedLogDensity,Tκ <: KineticEnergy, Tϵ <: Union{Real,Nothing}}
    Q::TQ
    κ::Tκ
    ϵ::Tϵ
end

function Base.show(io::IO, warmup_state::WarmupState)
    @unpack κ, ϵ = warmup_state
    ϵ_display = ϵ ≡ nothing ? "unspecified" : "≈ $(round(ϵ; sigdigits = REPORT_SIGDIGITS))"
    print(io, "adapted sampling parameters: stepsize (ϵ) $(ϵ_display), $(κ)")
end

###
### warmup interface and stages
###

"""
$(SIGNATURES)

Return the *results* and the *next warmup state* after warming up/adapting according to
`warmup_stage`, starting from `warmup_state`.

Use `nothing` for a no-op.
"""
function warmup(sampling_logdensity::SamplingLogDensity, warmup_stage::Nothing, warmup_state)
    nothing, warmup_state
end

"""
$(SIGNATURES)

Helper function to create random starting positions in the `[-2,2]ⁿ` box.
"""
random_position(rng, N) = rand(rng, N) .* 4 .- 2

"Docstring for initial warmup arguments."
const DOC_INITIAL_WARMUP_ARGS =
"""
- `q`: initial position. *Default*: random (uniform [-2,2] for each coordinate).

- `κ`: kinetic energy specification. *Default*: Gaussian with identity matrix.

- `ϵ`: a scalar for initial stepsize, or `nothing` for heuristic finders.
"""

"""
$(SIGNATURES)

Create an initial warmup state from a random position.

# Keyword arguments

$(DOC_INITIAL_WARMUP_ARGS)
"""
function initialize_warmup_state(rng, ℓ; q = random_position(rng, dimension(ℓ)),
                                 κ = GaussianKineticEnergy(dimension(ℓ)), ϵ = nothing)
    WarmupState(evaluate_ℓ(ℓ, q), κ, ϵ)
end

"""
$(TYPEDEF)

Find a local optimum (using quasi-Newton methods).

It is recommended that this stage is applied so that the initial stepsize selection happens
in a region which is at least plausible.
"""
Base.@kwdef struct FindLocalOptimum{T}
    """
    Add `-0.5 * magnitude_penalty * sum(abs2, q)` to the log posterior **when finding the local
    optimum**. This can help avoid getting into high-density edge areas of the posterior
    which are otherwise not typical (eg multilevel models).
    """
    magnitude_penalty::T = 1e-4
    """
    Maximum number of iterations in the optimization algorithm. Recall that we don't need to
    find the mode, or even a local mode, just be in a reasonable region.
    """
    iterations::Int = 20
    # FIXME allow custom algorithm, tolerance, etc
    # FIXME verify >0 of parameters in constructor
end

function warmup(sampling_logdensity, local_optimization::FindLocalOptimum, warmup_state)
    @unpack ℓ, reporter = sampling_logdensity
    @unpack magnitude_penalty, iterations = local_optimization
    @unpack Q, κ, ϵ = warmup_state
    @unpack q = Q
    report(reporter, "finding initial optimum")
    fg! = function(F, G, q)
        # NOTE: Optim.optimize *minimizes*, so we *add* a penalty
        ℓq, ∇ℓq = logdensity_and_gradient(ℓ, q)
        if G ≠ nothing
            @. G = -∇ℓq + q * magnitude_penalty
        end
        -ℓq + (0.5 * magnitude_penalty * sum(abs2, q))
    end
    objective = NLSolversBase.OnceDifferentiable(NLSolversBase.only_fg!(fg!), q)
    opt = Optim.optimize(objective, q, Optim.LBFGS(),
                         Optim.Options(; iterations = iterations))
    q = Optim.minimizer(opt)
    nothing, WarmupState(evaluate_ℓ(ℓ, q), κ, ϵ)
end

function warmup(sampling_logdensity, stepsize_search::InitialStepsizeSearch, warmup_state)
    @unpack rng, ℓ, reporter = sampling_logdensity
    @unpack Q, κ, ϵ = warmup_state
    @argcheck ϵ ≡ nothing "stepsize ϵ manually specified, won't perform initial search"
    z = PhasePoint(Q, rand_p(rng, κ))
    ϵ = find_initial_stepsize(stepsize_search, local_acceptance_ratio(Hamiltonian(κ, ℓ), z))
    report(reporter, "found initial stepsize",
           ϵ = round(ϵ; sigdigits = REPORT_SIGDIGITS))
    nothing, WarmupState(Q, κ, ϵ)
end

"""
$(TYPEDEF)

Tune the step size `ϵ` during sampling, and the metric of the kinetic energy at the end of
the block. The method for the latter is determined by the type parameter `M`, which can be

1. `Diagonal` for diagonal metric (the default),

2. `Symmetric` for a dense metric,

3. `Nothing` for an unchanged metric.

# Results

A `NamedTuple` with the following fields:

- `chain`, a vector of position vectors

- `tree_statistics`, a vector of tree statistics for each sample

- `ϵs`, a vector of step sizes for each sample

# Fields

$(FIELDS)
"""
struct TuningNUTS{M,D}
    "Number of samples."
    N::Int
    "Dual averaging parameters."
    stepsize_adaptation::D
    """
    Regularization factor for normalizing variance. An estimated covariance matrix `Σ` is
    rescaled by `λ` towards ``σ²I``, where ``σ²`` is the median of the diagonal. The
    constructor has a reasonable default.
    """
    λ::Float64
    function TuningNUTS{M}(N::Integer, stepsize_adaptation::D,
                           λ = 5.0/N) where {M <: Union{Nothing,Diagonal,Symmetric},D}
        @argcheck N ≥ 20        # variance estimator is kind of meaningless for few samples
        @argcheck λ ≥ 0
        new{M,D}(N, stepsize_adaptation, λ)
    end
end

function Base.show(io::IO, tuning::TuningNUTS{M}) where {M}
    @unpack N, stepsize_adaptation, λ = tuning
    print(io, "Stepsize and metric tuner, $(N) samples, $(M) metric, regularization $(λ)")
end

"""
$(SIGNATURES)

Form a matrix from positions (`q`), with each column containing a position.
"""
position_matrix(chain) = reduce(hcat, chain)

"""
$(SIGNATURES)

Estimate the inverse metric from the chain.

In most cases, this should be regularized, see [`regularize_M⁻¹`](@ref).
"""
sample_M⁻¹(::Type{Diagonal}, chain) = Diagonal(vec(var(position_matrix(chain); dims = 2)))

sample_M⁻¹(::Type{Symmetric}, chain) = Symmetric(cov(position_matrix(chain); dims = 2))

"""
$(SIGNATURES)

Adjust the inverse metric estimated from the sample, using an *ad-hoc* shrinkage method.
"""
function regularize_M⁻¹(Σ::Union{Diagonal,Symmetric}, λ::Real)
    # ad-hoc “shrinkage estimator”
    (1 - λ) * Σ + λ * UniformScaling(max(1e-3, median(diag(Σ))))
end

function warmup(sampling_logdensity, tuning::TuningNUTS{M}, warmup_state) where {M}
    @unpack rng, ℓ, algorithm, reporter = sampling_logdensity
    @unpack Q, κ, ϵ = warmup_state
    @unpack N, stepsize_adaptation, λ = tuning
    chain = Vector{typeof(Q.q)}(undef, N)
    tree_statistics = Vector{TreeStatisticsNUTS}(undef, N)
    H = Hamiltonian(κ, ℓ)
    ϵ_state = initial_adaptation_state(stepsize_adaptation, ϵ)
    ϵs = Vector{Float64}(undef, N)
    mcmc_reporter = make_mcmc_reporter(reporter, N; tuning = M ≡ Nothing ? "stepsize" :
                                       "stepsize and $(M) metric")
    for i in 1:N
        ϵ = current_ϵ(ϵ_state)
        ϵs[i] = ϵ
        Q, stats = sample_tree(rng, algorithm, H, Q, ϵ)
        chain[i] = Q.q
        tree_statistics[i] = stats
        ϵ_state = adapt_stepsize(stepsize_adaptation, ϵ_state, stats.acceptance_rate)
        report(mcmc_reporter, i; ϵ = round(ϵ; sigdigits = REPORT_SIGDIGITS))
    end
    if M ≢ Nothing
        κ = GaussianKineticEnergy(regularize_M⁻¹(sample_M⁻¹(M, chain), λ))
        report(mcmc_reporter, "adaptation finished", adapted_kinetic_energy = κ)
    end
    ((chain = chain, tree_statistics = tree_statistics, ϵs = ϵs),
     WarmupState(Q, κ, final_ϵ(ϵ_state)))
end

"""
$(TYPEDEF)

A composite type for performing MCMC stepwise after warmup.

The type is *not* part of the API, see [`mcmc_steps`](@ref) and [`mcmc_next_step`](@ref).
"""
struct MCMCSteps{TR,TA,TH,TE}
    rng::TR
    algorithm::TA
    H::TH
    ϵ::TE
end

"""
$(SIGNATURES)

Return a value which can be used to perform MCMC stepwise, eg until some criterion is
satisfied about the sample. See [`mcmc_next_step`](@ref).

Two constructors are available:

1. Explicitly providing
    - `rng`, the random number generator,
    - `algorithm`, see [`mcmc_with_warmup`](@ref),
    - `κ`, the (adapted) metric,
    - `ℓ`, the log density callable (see [`mcmc_with_warmup`](@ref),
    - `ϵ`, the stepsize.

2. Using the fields `sampling_logdensity` and `warmup_state`, eg from
    [`mcmc_keep_warmup`](@ref) (make sure you use eg `final_warmup_state`).

# Example

```julia
# initialization
results = DynamicHMC.mcmc_keep_warmup(RNG, ℓ, 0; reporter = NoProgressReport())
steps = mcmc_steps(results.sampling_logdensity, results.final_warmup_state)
Q = results.final_warmup_state.Q

# a single update step
Q, tree_stats = mcmc_next_step(steps, Q)

# extract the position
Q.q
```
"""
mcmc_steps(rng, algorithm, κ, ℓ, ϵ) = MCMCSteps(rng, algorithm, Hamiltonian(κ, ℓ), ϵ)

function mcmc_steps(sampling_logdensity::SamplingLogDensity, warmup_state)
    @unpack rng, ℓ, algorithm = sampling_logdensity
    @unpack κ, ϵ = warmup_state
    mcmc_steps(rng, algorithm, κ, ℓ, ϵ)
end

"""
$(SIGNATURES)

Given `Q` (an evaluated log density at a position), return the next `Q` and tree statistics.
"""
function mcmc_next_step(mcmc_steps::MCMCSteps, Q::EvaluatedLogDensity)
    @unpack rng, algorithm, H, ϵ = mcmc_steps
    sample_tree(rng, algorithm, H, Q, ϵ)
end

"""
$(SIGNATURES)

Markov Chain Monte Carlo for `sampling_logdensity`, with the adapted `warmup_state`.

Return a `NamedTuple` of

- `chain`, a vector of length `N` that contains the positions,

- `tree_statistics`, a vector of length `N` with the tree statistics.
"""
function mcmc(sampling_logdensity, N, warmup_state)
    @unpack reporter = sampling_logdensity
    @unpack Q = warmup_state
    chain = Vector{typeof(Q.q)}(undef, N)
    tree_statistics = Vector{TreeStatisticsNUTS}(undef, N)
    mcmc_reporter = make_mcmc_reporter(reporter, N)
    steps = mcmc_steps(sampling_logdensity, warmup_state)
    for i in 1:N
        Q, tree_statistics[i] = mcmc_next_step(steps, Q)
        chain[i] = Q.q
        report(mcmc_reporter, i)
    end
    (chain = chain, tree_statistics = tree_statistics)
end

"""
$(SIGNATURES)

Helper function for constructing the “middle” doubling warmup stages in
[`default_warmup_stages`](@ref).
"""
function _doubling_warmup_stages(M, stepsize_adaptation, middle_steps,
                                 doubling_stages::Val{D}) where {D}
    ntuple(i -> TuningNUTS{M}(middle_steps * 2^(i - 1), stepsize_adaptation), D)
end

"""
$(SIGNATURES)

A sequence of warmup stages:

1. find the local optimum using `local_optimization`,

2. select an initial stepsize using `stepsize_search` (default: based on a heuristic),

3. tuning stepsize with `init_steps` steps

4. tuning stepsize and covariance: first with `middle_steps` steps, then repeat with twice
   the steps `doubling_stages` times

5. tuning stepsize with `terminating_steps` steps.

`M` (`Diagonal`, the default or `Symmetric`) determines the type of the metric adapted from
the sample.

This is the suggested tuner of most applications.

Use `nothing` for `local_optimization` or `stepsize_adaptation` to skip the corresponding
step.
"""
function default_warmup_stages(;
                               local_optimization = FindLocalOptimum(),
                               stepsize_search = InitialStepsizeSearch(),
                               M::Type{<:Union{Diagonal,Symmetric}} = Diagonal,
                               stepsize_adaptation = DualAveraging(),
                               init_steps = 75, middle_steps = 25, doubling_stages = 5,
                               terminating_steps = 50)
    (local_optimization, stepsize_search,
     TuningNUTS{Nothing}(init_steps, stepsize_adaptation),
     _doubling_warmup_stages(M, stepsize_adaptation, middle_steps, Val(doubling_stages))...,
     TuningNUTS{Nothing}(terminating_steps, stepsize_adaptation))
end

"""
$(SIGNATURES)

A sequence of warmup stages for fixed stepsize:

1. find the local optimum using `local_optimization`,

2. tuning covariance: first with `middle_steps` steps, then repeat with twice
   the steps `doubling_stages` times

Very similar to [`default_warmup_stages`](@ref), but omits the warmup stages with just
stepsize tuning.
"""
function fixed_stepsize_warmup_stages(;
                                      local_optimization = FindLocalOptimum(),
                                      M::Type{<:Union{Diagonal,Symmetric}} = Diagonal,
                                      middle_steps = 25, doubling_stages = 5)
    (local_optimization,
     _doubling_warmup_stages(M, FixedStepsize(), middle_steps, Val(doubling_stages))...)
end

"""
$(SIGNATURES)

Helper function for implementing warmup.

!!! note
    Changes may imply documentation updates in [`mcmc_keep_warmup`](@ref).
"""
function _warmup(sampling_logdensity, stages, initial_warmup_state)
    foldl(stages; init = ((), initial_warmup_state)) do acc, stage
        stages_and_results, warmup_state = acc
        results, warmup_state′ = warmup(sampling_logdensity, stage, warmup_state)
        stage_information = (stage = stage, results = results, warmup_state = warmup_state′)
        (stages_and_results..., stage_information), warmup_state′
    end
end

"Shared docstring part for the MCMC API."
const DOC_MCMC_ARGS =
"""
# Arguments

- `rng`: the random number generator, eg `Random.GLOBAL_RNG`.

- `ℓ`: the log density, supporting the API of the `LogDensityProblems` package

- `N`: the number of samples for inference, after the warmup.

# Keyword arguments

- `initialization`: see below.

- `warmup_stages`: a sequence of warmup stages. See [`default_warmup_stages`](@ref) and
  [`fixed_stepsize_warmup_stages`](@ref); the latter requires an `ϵ` in initialization.

- `algorithm`: see [`NUTS`](@ref). It is very unlikely you need to modify
  this, except perhaps for the maximum depth.

- `reporter`: how progress is reported. By default, verbosely for interactive sessions using
  the log message mechanism (see [`LogProgressReport`](@ref), and no reporting for
  non-interactive sessions (see [`NoProgressReport`](@ref)).

# Initialization

The `initialization` keyword argument should be a `NamedTuple` which can contain the
following fields (all of them optional and provided with reasonable defaults):

$(DOC_INITIAL_WARMUP_ARGS)
"""

"""
$(SIGNATURES)

Perform MCMC with NUTS, keeping the warmup results. Returns a `NamedTuple` of

- `initial_warmup_state`, which contains the initial warmup state

- `warmup`, an iterable of `NamedTuple`s each containing fields

    - `stage`: the relevant warmup stage

    - `results`: results returned by that warmup stage (may be `nothing` if not applicable,
      or a chain, with tree statistics, etc; see the documentation of stages)

    - `warmup_state`: the warmup state *after* the corresponding stage.

- `final_warmup_state`, which contains the final adaptation after all the warmup

- `inference`, which has `chain` and `tree_statistics`, see [`mcmc_with_warmup`](@ref).

- `sampling_logdensity`, which contains information that is invariant to warmup

!!! warning
    This function is not (yet) exported because the the warmup interface may change with
    minor versions without being considered breaking. Recommended for interactive use.

$(DOC_MCMC_ARGS)
"""
function mcmc_keep_warmup(rng::AbstractRNG, ℓ, N::Integer;
                          initialization = (),
                          warmup_stages = default_warmup_stages(),
                          algorithm = NUTS(),
                          reporter = default_reporter())
    sampling_logdensity = SamplingLogDensity(rng, ℓ, algorithm, reporter)
    initial_warmup_state = initialize_warmup_state(rng, ℓ; initialization...)
    warmup, warmup_state = _warmup(sampling_logdensity, warmup_stages, initial_warmup_state)
    inference = mcmc(sampling_logdensity, N, warmup_state)
    (initial_warmup_state = initial_warmup_state, warmup = warmup,
     final_warmup_state = warmup_state, inference = inference,
     sampling_logdensity = sampling_logdensity)
end

"""
$(SIGNATURES)

Perform MCMC with NUTS, including warmup which is not returned. Return a `NamedTuple` of

- `chain`, a vector of positions from the posterior

- `tree_statistics`, a vector of tree statistics

- `κ` and `ϵ`, the adapted metric and stepsize.

$(DOC_MCMC_ARGS)

# Usage examples

Using a fixed stepsize:
```julia
mcmc_with_warmup(rng, ℓ, N;
                 initialization = (ϵ = 0.1, ),
                 warmup_stages = fixed_stepsize_warmup_stages())
```

Starting from a given position `q₀` and kinetic energy scaled down (will still be adapted):
```julia
mcmc_with_warmup(rng, ℓ, N;
                 initialization = (q = q₀, κ = GaussianKineticEnergy(5, 0.1)))
```

Using a dense metric:
```julia
mcmc_with_warmup(rng, ℓ, N;
                 warmup_stages = default_warmup_stages(; M = Symmetric))
```

Disabling the optimization step:
```julia
mcmc_with_warmup(rng, ℓ, N;
                 warmup_stages = default_warmup_stages(; local_optimization = nothing,
                                                         M = Symmetric))
```

Disabling the initial stepsize search (provided explicitly, still adapted):
```julia
mcmc_with_warmup(rng, ℓ, N;
                 initialization = (ϵ = 1.0, ),
                 warmup_stages = default_warmup_stages(; stepsize_search = nothing))
```
"""
function mcmc_with_warmup(rng, ℓ, N; initialization = (),
                          warmup_stages = default_warmup_stages(),
                          algorithm = NUTS(), reporter = default_reporter())
    @unpack final_warmup_state, inference =
        mcmc_keep_warmup(rng, ℓ, N; initialization = initialization,
                         warmup_stages = warmup_stages, algorithm = algorithm,
                         reporter = reporter)
    @unpack κ, ϵ = final_warmup_state
    (inference..., κ = κ, ϵ = ϵ)
end


text:
import Distributions

function updateOutput!{T <: FloatingPoint}(
    layer::SpatialConv,
    input::Array{T},
)
    layer.output = resize!(vec(layer.output), layer.nOutputPlane*(size(input,2)-layer.kH+1)*(size(input,3)-layer.kW+1))
    layer.output = reshape(layer.output,layer.nOutputPlane, (size(input,2)-layer.kH+1), (size(input,3)-layer.kW+1))
    print(size(layer.output))
    for i = 1:layer.nOutputPlane
        output_i = sub(layer.output, :, :, i)
        for j = 1:layer.nInputPlane
            input_i = sub(input, :, :, i)
            kernel_i = sub(layer.weight, i, j, :, :)
            output_i += conv2(input_i, kernel_i)
        end
    end
    return output
end

function forward!{T <: FloatingPoint}(
    layer::Layer,
    input::Array{T},
)
    return updateOutput!(layer, input)
end

function updateGradInput!{T <: FloatingPoint}(
    layer::Layer,
    input::Array{T},
    gradOutput::Array{T}
)
    error("updateGradInput!(layer, input, gradOutput) not implemented abstractly")
end

function accGradParameters!{T <: FloatingPoint}(
    layer::Layer,
    input::Array{T},
    gradOutput::Array{T},
)
    error("accGradParameters!(layer, input, gradOutput) not implemented abstractly")
end

function SpatialConv{T <: FloatingPoint}(
                                         ::Type{T},
                                         nInputPlane::Integer,
                                         nOutputPlane::Integer,
                                         kernelHeight::Integer,
                                         kernelWidth::Integer,
                                         strideHeight::Integer,
                                         strideWidth::Integer
                                         )
    obj = SpatialConv(
                      nInputPlane,
                      nOutputPlane,
                      kernelHeight,
                      kernelWidth,
                      strideHeight,
                      strideWidth,
                      Array(T, nOutputPlane, nInputPlane, kernelHeight, kernelWidth),
                      Array(T, nOutputPlane),
                      Array(T, nOutputPlane, nInputPlane, kernelHeight, kernelWidth),
                      Array(T, nOutputPlane),
    Array(T,1,1,1),
    Array(T,1,1,1)
                      )
    # fill weight and bias with a Uniform distribution normalized by fanin
    std = 1/sqrt(kernelWidth*kernelHeight*nInputPlane)
    rand!(Distributions.Uniform(-std, std), obj.weight)
    rand!(Distributions.Uniform(-std, std), obj.bias)
    return obj
end

ext: jl


text:
# TODO compress Union{Int, Missing} is no missing present

# TODO allow int to uint
"""
	type_compress!(df, compress_float = false, verbose = false)

Compress a DataFrame by using types of smaller bit. If safe to do so, it will
"downgrade" `Int*` and `UInt*` if safe to do some. It will compress any
`CategoricalVector` with `DataFrames.compress`.

For `Vector{String}`, if the number unique values is less than 2^16 then it will
be converted to `CategoricalVector` and otherwise will be stored as
`WeakRefStrings.StringVector`.

If `compress_float = true` then `Float64` will be downgraded to `Float32`; but
beware that this means the calculation will be done with reduce precision.
"""
type_compress!(df::DataFrame; compress_float = false, verbose = false) = begin
	for n in names(df)
		if verbose
			println("Compressing $n")
		end
		if compress_float || (nonmissingtype(eltype(df[!,n])) != Float64)
			df[!,n] = type_compress(df[!,n])
		end
	end
	df
end

# for CSV.Column read from CSV
type_compress(v::CC) where CC <: CSV.Column = type_compress(Vector(v))

#
type_compress(v::Vector{T}) where T <: Union{Int128, Int64, Int32, Int16} = begin
	min1, max1 = extrema(v)

	if typemin(Int8) <= min1 && max1 <= typemax(Int8)
		return Int8.(v)
	elseif typemin(Int16) <= min1 && max1 <= typemax(Int16)
		return Int16.(v)
	elseif typemin(Int32) <= min1 && max1 <= typemax(Int32)
		return Int32.(v)
	elseif typemin(Int64) <= min1 && max1 <= typemax(Int64)
		return Int64.(v)
	end
end

type_compress(v::Vector{T}) where T <: Union{UInt128, UInt64, UInt32, UInt16} = begin
	max1 = maximum(v)

	if max1 <= typemax(UInt8)
		return UInt8.(v)
	elseif max1 <= typemax(UInt16)
		return UInt16.(v)
	elseif max1 <= typemax(UInt32)
		return UInt32.(v)
	elseif max1 <= typemax(UInt64)
		return UInt64.(v)
	end
end

type_compress(v::Vector{Union{Missing, T}}) where T <: Union{UInt128, UInt64, UInt32, UInt16} = begin
	max1 = maximum(skipmissing(v))

	if max1 <= typemax(UInt8)
		return Vector{Union{Missing, UInt8}}(v)
	elseif max1 <= typemax(UInt16)
		return Vector{Union{Missing, UInt16}}(v)
	elseif max1 <= typemax(UInt32)
		return Vector{Union{Missing, UInt32}}(v)
	elseif max1 <= typemax(UInt64)
		return Vector{Union{Missing, UInt64}}(v)
	end
end

type_compress(v::Vector{Union{Missing, T}}) where T <: Union{Int128, Int64, Int32, Int16} = begin
	min1, max1 = extrema(skipmissing(v))

	if typemin(Int8) <= min1 && max1 <= typemax(Int8)
		return Vector{Union{Missing, Int8}}(v)
	elseif typemin(Int16) <= min1 && max1 <= typemax(Int16)
		return Vector{Union{Missing, Int16}}(v)
	elseif typemin(Int32) <= min1 && max1 <= typemax(Int32)
		return Vector{Union{Missing, Int32}}(v)
	elseif typemin(Int64) <= min1 && max1 <= typemax(Int64)
		return Vector{Union{Missing, Int64}}(v)
	end
end

type_compress(v::Vector{Float64}) = Vector{Float32}(v)
type_compress(v::Vector{Union{Missing, Float64}}) = Vector{Union{Missing, Float32}}(v)

type_compress(v::Vector{String}) = begin
	# TODO recommend a PooledString if necessary
	# use some heuristic to decide whether to compress

	## if the number of unique elements is predicted to be less than typemax(UInt16) + 1
	## then we can use categorical array to compress them
	# sv = countmap(sample(v, 888))
	#
	# # Estimate the number of unique items
	# # https://stats.stackexchange.com/questions/19014/how-can-i-estimate-unique-occurrence-counts-from-a-random-sampling-of-data
	# u1 = length([key for (key, value) in sv if value == 1])
	# u = length(sv)
	# # if the estimated number of estimtes is less than about 32000 then it's
	# # worth putting it into a categorical array
	# if u + u1/888 * (length(v) - 888) < typemax(UInt16) + 1
	# 	return categorical(v)
	# 	# rlev = rle(categorical(v))
	# 	# rlev[2] .= cumsum(rlev[2])
	# 	# return RLEVector(rlev...)
	# end

	if length(Set(v)) <= 2^16
		return categorical(v)
	end

	# check if the string is more compact at RLE level
	rlev = rle(v)

	# firstly check if the rlev can be whittle down further
	offset_dict = Dict{String, UInt}()

	# an IOBuffer to build up the string buffer
	long_string_builder = IOBuffer()
	offset_array_builder = IOBuffer()
	string_len_builder = IOBuffer()

	values = rlev[1]
	offset_upto = 0

	# add the first entry
	offset_dict[values[1]] = offset_upto

	write(offset_array_builder, offset_upto)
	offset_update = write(long_string_builder, values[1])
	write(string_len_builder, offset_update)

	offset_upto += offset_update

	for v in values[2:end]
		if !haskey(offset_dict, v)
			write(offset_array_builder, offset_upto)
			offset_update = write(long_string_builder, v)
			write(string_len_builder, offset_update)

			# update the offset_upto for next one
			offset_dict[v] = offset_upto
			offset_upto += offset_update
			# println(v, " : ", offset_dict[v])
		else
			# If the string doesn't existing then obtain the offset from the
			# offset_dict.
			# There is no need to write it to the string aagain because it's
			# already there

			offset_old = offset_dict[v]
			# still write the lengths
			write(string_len_builder, sizeof(v))
			# still write the offset
			write(offset_array_builder, offset_old)
		end
	end

	#long_str =String(take!(long_string_builder))
	offsets = inverse_rle(reinterpret(Int, take!(offset_array_builder)), rlev[2])
	lengths = inverse_rle(reinterpret(Int, take!(string_len_builder  )), rlev[2])

	buffer = take!(long_string_builder)

	close(long_string_builder)
	close(offset_array_builder)
	close(string_len_builder)

	#return String(buffer), offsets, lengths
	#return (StringArray{String, 1}(buffer, offsets, lengths), String(buffer), offsets, lengths)
	return StringArray{String, 1}(buffer, offsets, lengths)

end

type_compress(v::Vector{Missing}) = v

type_compress(v::CategoricalVector) where {T, IntType}  = compress(v)

type_compress(v) = begin
	# println("The compression for $(typeof(v)) is not yet supported. No compression is performed. Submit an idea if you think JDF should support it: https://github.com/xiaodaigh/JDF.jl/issues")
	v
end

if false

end


text:
#
# macros.jl -
#
# Definition of private macros for Julia Tcl interface.
#

# The following macro is used to build a list out of given arguments and
# keywords by the `exec` and `list` methods.
macro __build_list(listptr, args, kwds)
    quote
        for arg in $(esc(args))
            __lappend($(esc(listptr)), arg)
        end
        for (key, val) in $(esc(kwds))
            __lappendoption($(esc(listptr)), key, val)
        end
    end
end

# The following macro is used to build a list by concatenating given arguments
# keywords by the `eval` and `concat` methods.
macro __concat_args(listptr, args)
    quote
        for arg in $(esc(args))
            __concat($(esc(listptr)), arg)
        end
    end
end


text:
"""
    PointEvalHandler(grid::Grid, points::AbstractVector{Vec{dim,T}}) where {dim, T}

The `PointEvalHandler` can be used for function evaluation in *arbitrary points* in the
domain -- not just in quadrature points or nodes.

The constructor takes a grid and a vector of coordinates for the points. The
`PointEvalHandler` computes i) the corresponding cell, and ii) the (local) coordinate
within the cell, for each point. The fields of the `PointEvalHandler` are:
 - `cells::Vector{Union{Int,Nothing}}`: vector with cell IDs for the points, with `nothing`
   for points that could not be found.
 - `local_coords::Vector{Union{Vec,Nothing}}`: vector with the local coordinates
   (i.e. coordinates in the reference configuration) for the points, with `nothing` for
   points that could not be found.

There are two ways to use the `PointEvalHandler` to evaluate functions:

 - [`get_point_values`](@ref): can be used when the function is described by
   i) a `dh::DofHandler` + `uh::Vector` (for example the FE-solution), or
   ii) a `p::L2Projector` + `ph::Vector` (for projected data).
 - Iteration with [`PointIterator`](@ref) + [`PointValues`](@ref): can be used for more
   flexible evaluation in the points, for example to compute gradients.
"""
PointEvalHandler

struct PointEvalHandler{G,dim,T<:Real}
    grid::G
    cells::Vector{Union{Nothing, Int}}
    local_coords::Vector{Union{Nothing, Vec{dim,T}}}
end

function Base.show(io::IO, ::MIME"text/plain", ph::PointEvalHandler)
    println(io, typeof(ph))
    println(io, "  number of points: ", length(ph.local_coords))
    n_missing = sum(x -> x === nothing, ph.cells)
    if n_missing == 0
        print(io, "  Found corresponding cell for all points.")
    else
        print(io, "  Could not find corresponding cell for ", n_missing, " points.")
    end
end

function PointEvalHandler(grid::AbstractGrid, points::AbstractVector{Vec{dim,T}}) where {dim, T}
    node_cell_dicts = _get_node_cell_map(grid)
    cells, local_coords = _get_cellcoords(points, grid, node_cell_dicts)
    return PointEvalHandler(grid, cells, local_coords)
end

function _get_cellcoords(points::AbstractVector{Vec{dim,T}}, grid::AbstractGrid, node_cell_dicts::Dict{C,Dict{Int, Vector{Int}}}) where {dim, T<:Real, C}

    # set up tree structure for finding nearest nodes to points
    kdtree = KDTree(reinterpret(Vec{dim,T}, getnodes(grid)))
    nearest_nodes, _ = knn(kdtree, points, 3, true) #TODO 3 is a random value, it shouldn't matter because likely the nearest node is the one we want

    cells = Vector{Union{Nothing, Int}}(nothing, length(points))
    local_coords = Vector{Union{Nothing, Vec{dim, T}}}(nothing, length(points))

    for point_idx in 1:length(points)
        cell_found = false
        for (CT, node_cell_dict) in node_cell_dicts
            geom_interpol = default_interpolation(CT)
            # loop over points
            for node in nearest_nodes[point_idx]
                possible_cells = get(node_cell_dict, node, nothing)
                possible_cells === nothing && continue # if node is not part of the fieldhandler, try the next node
                for cell in possible_cells
                    cell_coords = getcoordinates(grid, cell)
                    is_in_cell, local_coord = point_in_cell(geom_interpol, cell_coords, points[point_idx])
                    if is_in_cell
                        cell_found = true
                        cells[point_idx] = cell
                        local_coords[point_idx] = local_coord
                        break
                    end
                end
                cell_found && break
            end
            cell_found && break
        end
        if !cell_found
            @warn("No cell found for point number $point_idx, coordinate: $(points[point_idx]).")
        end
    end
    return cells, local_coords
end

# check if point is inside a cell based on physical coordinate
function point_in_cell(geom_interpol::Interpolation{dim,shape,order}, cell_coordinates, global_coordinate) where {dim, shape, order}
    converged, x_local = find_local_coordinate(geom_interpol, cell_coordinates, global_coordinate)
    if converged
        return _check_isoparametric_boundaries(shape, x_local), x_local
    else
        return false, x_local
    end
end

# check if point is inside a cell based on isoparametric coordinate
function _check_isoparametric_boundaries(::Type{RefCube}, x_local::Vec{dim, T}) where {dim, T}
    tol = sqrt(eps(T))
    # All in the range [-1, 1]
    return all(x -> abs(x) - 1 < tol, x_local)
end

# check if point is inside a cell based on isoparametric coordinate
function _check_isoparametric_boundaries(::Type{RefTetrahedron}, x_local::Vec{dim, T}) where {dim, T}
    tol = sqrt(eps(T))
    # Positive and below the plane 1 - ξx - ξy - ξz
    return all(x -> x > -tol, x_local) && sum(x_local) - 1 < tol
end

# TODO: should we make iteration params optional keyword arguments?
function find_local_coordinate(interpolation, cell_coordinates, global_coordinate)
    """
    currently copied verbatim from https://discourse.julialang.org/t/finding-the-value-of-a-field-at-a-spatial-location-in-juafem/38975/2
    other than to make J dim x dim rather than 2x2
    """
    dim = length(global_coordinate)
    local_guess = zero(Vec{dim})
    n_basefuncs = getnbasefunctions(interpolation)
    max_iters = 10
    tol_norm = 1e-10
    converged = false
    for _ in 1:max_iters
        N = Ferrite.value(interpolation, local_guess)

        global_guess = zero(Vec{dim})
        for j in 1:n_basefuncs
            global_guess += N[j] * cell_coordinates[j]
        end
        residual = global_guess - global_coordinate
        if norm(residual) <= tol_norm
            converged = true
            break
        end
        dNdξ = Ferrite.derivative(interpolation, local_guess)
        J = zero(Tensor{2, dim})
        for j in 1:n_basefuncs
            J += cell_coordinates[j] ⊗ dNdξ[j]
        end
        local_guess -= inv(J) ⋅ residual
    end
    return converged, local_guess
end

# return a Dict with a key for each node that contains a vector with the adjacent cells as value
function _get_node_cell_map(grid::AbstractGrid)
    cells = getcells(grid)
    C = eltype(cells) # possibly abstract
    cell_dicts = Dict{Type{<:C}, Dict{Int, Vector{Int}}}()
    ctypes = Set{Type{<:C}}(typeof(c) for c in cells)
    for ctype in ctypes
        cell_dict = cell_dicts[ctype] = Dict{Int,Vector{Int}}()
        for (cellidx, cell) in enumerate(cells)
            cell isa ctype || continue
            for node in cell.nodes
                v = get!(Vector{Int}, cell_dict, node)
                push!(v, cellidx)
            end
        end
    end
    return cell_dicts
end

"""
    get_point_values(ph::PointEvalHandler, dh::AbstractDofHandler, dof_values::Vector{T}, [fieldname::Symbol]) where T
    get_point_values(ph::PointEvalHandler, proj::L2Projector, dof_values::Vector{T}) where T

Return a `Vector{T}` (for a 1-dimensional field) or a `Vector{Vec{fielddim, T}}` (for a
vector field) with the field values of field `fieldname` in the points of the
`PointEvalHandler`. The `fieldname` can be omitted if only one field is stored in `dh`.
The field values are computed based on the `dof_values` and interpolated to the local
coordinates by the function interpolation of the corresponding `field` stored in the
`AbstractDofHandler` or the `L2Projector`.

Points that could not be found in the domain when constructing the `PointEvalHandler` will
have `NaN`s for the corresponding entries in the output vector.
"""
get_point_values

function get_point_values(ph::PointEvalHandler, proj::L2Projector, dof_vals::AbstractVector)
    get_point_values(ph, proj.dh, dof_vals)
end

function get_point_values(ph::PointEvalHandler, dh::AbstractDofHandler, dof_vals::AbstractVector{T},
                           fname::Symbol=find_single_field(dh)) where {T}
    fdim = getfielddim(dh, fname)
    npoints = length(ph.cells)
    # for a scalar field return a Vector of Scalars, for a vector field return a Vector of Vecs
    if fdim == 1
        out_vals = fill!(Vector{T}(undef, npoints), NaN * zero(T))
    else
        nanv = convert(Vec{fdim,T}, NaN * zero(Vec{fdim,T}))
        out_vals = fill!(Vector{Vec{fdim, T}}(undef, npoints), nanv)
    end
    func_interpolations = get_func_interpolations(dh, fname)
    get_point_values!(out_vals, ph, dh, dof_vals, fname, func_interpolations)
    return out_vals
end
function find_single_field(dh)
    ns = getfieldnames(dh)
    if length(ns) != 1
        throw(ArgumentError("multiple fields in DoF handler, must specify which"))
    end
    return ns[1]
end

# values in dof-order. They must be obtained from the same DofHandler that was used for constructing the PointEvalHandler
function get_point_values!(out_vals::Vector{T2},
    ph::PointEvalHandler,
    dh::MixedDofHandler,
    dof_vals::Vector{T},
    fname::Symbol,
    func_interpolations
    ) where {T2, T}

    # TODO: I don't think this is correct??
    length(dof_vals) == ndofs(dh) || error("You must supply values for all $(ndofs(dh)) dofs.")

    fdim = getfielddim(dh, fname)

    for fh_idx in eachindex(dh.fieldhandlers)
        ip = func_interpolations[fh_idx]
        if ip !== nothing
            dofrange = dof_range(dh.fieldhandlers[fh_idx], fname)
            cellset = dh.fieldhandlers[fh_idx].cellset
            _get_point_values!(out_vals, dof_vals, ph, dh, ip, cellset, Val(fdim), dofrange)
        end
    end
    return out_vals
end

function get_point_values!(out_vals::Vector{T2},
    ph::PointEvalHandler,
    dh::DofHandler,
    dof_vals::Vector{T},
    fname::Symbol,
    func_interpolations
    ) where {T2, T}

    # TODO: I don't think this is correct??
    length(dof_vals) == ndofs(dh) || error("You must supply values for all $(ndofs(dh)) dofs.")

    fdim = getfielddim(dh, fname)
    dofrange = dof_range(dh, fname)
    _get_point_values!(out_vals, dof_vals, ph, dh, func_interpolations[1], nothing, Val(fdim), dofrange)
    return out_vals
end

# function barrier with concrete type of interpolation
function _get_point_values!(
    out_vals::Vector{T2},
    dof_vals::Vector{T},
    ph::PointEvalHandler,
    dh::AbstractDofHandler,
    ip::Interpolation,
    cellset::Union{Nothing, Set{Int}},
    fdim::Val{fielddim},
    dofrange::AbstractRange{Int},
    ) where {T2,T,fielddim}

    # extract variables
    local_coords = ph.local_coords
    # preallocate some stuff specific to this cellset
    pv = PointScalarValuesInternal(first(local_coords), ip)
    first_cell = cellset === nothing ? 1 : first(cellset)
    cell_dofs = Vector{Int}(undef, ndofs_per_cell(dh, first_cell))

    # compute point values
    for pointid in eachindex(ph.cells)
        cellid = ph.cells[pointid]
        cellid === nothing && continue # next point if no cell was found for this one
        cellset !== nothing && (cellid ∈ cellset || continue) # no need to check the cellset for a regular DofHandler
        celldofs!(cell_dofs, dh, ph.cells[pointid])
        reinit!(pv, local_coords[pointid], ip)
        @inbounds @views dof_vals_reshaped = _change_format(fdim, dof_vals[cell_dofs[dofrange]])
        out_vals[pointid] = function_value(pv, 1, dof_vals_reshaped)
    end
    return out_vals
end

###################################################################################################################
# utils

# reshape dof_values based on fielddim
_change_format(::Val{1}, dof_values::AbstractVector{T}) where T = dof_values
_change_format(::Val{fielddim}, dof_values::AbstractVector{T}) where {fielddim, T} = reinterpret(Vec{fielddim, T}, dof_values)

get_func_interpolations(dh::DH, fieldname) where DH<:DofHandler = [getfieldinterpolation(dh, find_field(dh, fieldname))]
function get_func_interpolations(dh::DH, fieldname) where DH<:MixedDofHandler
    func_interpolations = Union{Interpolation,Nothing}[]
    for fh in dh.fieldhandlers
        j = findfirst(i -> i === fieldname, getfieldnames(fh))
        if j === nothing
            push!(func_interpolations, missing)
        else
            push!(func_interpolations, fh.fields[j].interpolation)
        end
    end
    return func_interpolations
end

# Iteration of PointEvalHandler
"""
    PointIterator(ph::PointEvalHandler)

Create an iterator over the points in the [`PointEvalHandler`](@ref).
The elements of the iterator are either a [`PointLocation`](@ref), if the corresponding
point could be found in the grid, or `nothing`, if the point was not found.

A `PointLocation` can be used to query the cell ID with [`cellid`](@ref), and can be used
to reinitialize [`PointValues`](@ref) with [`reinit!`](@ref).

# Examples
```julia
ph = PointEvalHandler(grid, points)

for point in PointIterator(ph)
    point === nothing && continue # Skip any points that weren't found
    reinit!(pointvalues, point)   # Update pointvalues
    # ...
end
```
"""
PointIterator

struct PointIterator{PH<:PointEvalHandler}
    ph::PH
    coords::Vector{Vec{2,Float64}}
end

function PointIterator(ph::PointEvalHandler{G}) where {D,C,T,G<:Grid{D,C,T}}
    n = nnodes_per_cell(ph.grid)
    coords = zeros(Vec{D,T}, n) # resize!d later if needed
    return PointIterator(ph, coords)
end

"""
    PointLocation

Element of a [`PointIterator`](@ref), typically used to reinitialize
[`PointValues`](@ref). Fields:
 - `cid::Int`: ID of the cell containing the point
 - `local_coord::Vec`: the local (reference) coordinate of the point
 - `coords::Vector{Vec}`: the coordinates of the cell
"""
struct PointLocation{V}
    cid::Int
    local_coord::V
    coords::Vector{V}
end

function Base.iterate(p::PointIterator, state = 1)
    if state > length(p.ph.cells)
        return nothing
    elseif p.ph.cells[state] === nothing
        return (nothing, state + 1)
    else
        cid = (p.ph.cells[state])::Int
        local_coord = (p.ph.local_coords[state])::Vec
        n = nnodes_per_cell(p.ph.grid, cid)
        cellcoords!(resize!(p.coords, n), p.ph.grid, cid)
        point = PointLocation(cid, local_coord, p.coords)
        return (point, state + 1)
    end
end

cellid(p::PointLocation) = p.cid

function reinit!(pv::PointValues, point::PointLocation)
    reinit!(pv, point.coords, point.local_coord)
    return pv
end


text:
using QuantumFluidSpectra
using Test

@testset "Analysis" begin
    # Velocity and Helmholtz tests
    n = 100
    L = (1,1)
    N = (n,n)
    X,K,dX,dK = xk_arrays(L,N)
    kx,ky = K

    ##
    ktest = K[1][2]
    ψ = @. exp(im*ktest*X[1]*one.(X[2]'))
    psi = Psi(ψ,X,K)

    ## flow only in x direction, of correct value
    vx,vy = velocity(psi)
    @test vx ≈ ktest*one.(vx)
    @test vy ≈ zero.(vy)

    ## helmholtz decomposition
    Vi,Vc = helmholtz(vx,vy,kx,ky)

    ## Orthogonality
    vidotvc = Vi[1].*Vc[1] .+ Vi[2].*Vc[2]
    @test maximum(abs.(vidotvc)) < 1e-10

    ## Projective
    @test Vi[1] .+ Vc[1] ≈ vx
    @test Vi[2] .+ Vc[2] ≈ vy

    ## autocorrelation
    x,y = X; kx,ky = K
    dx = x[2]-x[1];dk = kx[2]-kx[1]

    Natoms = sum(abs2.(ψ))*dx^2
    AC = auto_correlate(ψ,X,K)
    Nr,Nim = AC[101,101] .|> (real,imag)

    @test Natoms ≈ Nr
    @test Nim ≈ 0.0

    ## cross correlated reduces to autocorrelate correctly
    CC = cross_correlate(ψ,ψ,X,K)
    @test AC == AC

end


text:
using Random
using SumBassSolver

Random.seed!(121)
n=50
p=0.001
q=0.1
m=100
c0=rand(Float64,n).*m./2
alpha = 1
F = 10
T = 100
Pi = 2500

x,v = sumBassDiffStartSolver(p,q,m,c0,alpha,F,T,Pi)

println(alpha .* Bass_derivative.(x,p,q,m,c0))

pp=fill(p,n)
qq=fill(q,n)
mm=fill(m,n)
alphas=fill(alpha,n)
FF=fill(F,n)
TT=fill(T,n)

x1,v1 = sumBassSolver(pp,qq,mm,c0,alphas,FF,TT,Pi)
x2,v2 = sumBassSolver(pp,qq,mm,c0,alphas,FF,TT,Pi,alg="SB")
x3,v3 = sumBassSolver(pp,qq,mm,c0,alphas,FF,TT,Pi,alg="SB",localopt=true)

max_stars_repo_name: andres-alban/res-alloc-bass-demands

max_issues_repo_name: andres-alban/res-alloc-bass-demands

max_forks_repo_name: andres-alban/res-alloc-bass-demands


text:
using Test
using NeXLParticle
using Random
using StaticArrays
using CoordinateTransformations, Rotations

@testset "ICP32" begin
    r, npts, ty = MersenneTwister(0xBADFEED), 100, Float32
    xy = SVector.(zip(rand(r, ty, npts), rand(r, ty, npts)))
    @test eltype(eltype(xy)) == ty
    # Muddle the points up
    xyp = AffineMap(Angle2d{ty}(deg2rad(1.2f0)), [0.1f0, -0.1f0]).(xy) .+ 0.02f0 * SVector.(collect(zip(0.5f0 .- rand(r, ty, npts), 0.5f0 .- rand(r, ty, npts))))
    xyps = shuffle(r, xyp)
    @test eltype(eltype(xyps)) == ty
    @test isapprox(NeXLParticle.icperror(xy, xyps), 5.79987f0, atol=0.001f0)
    @test isapprox(NeXLParticle.icperror(xy[1:2:100], xyps),7.61191f0,atol=0.001f0)
    @test isapprox(NeXLParticle.icperror(xy, xyps[1:2:100]), 2.85196f0, atol=0.001f0)
    res=NeXLParticle.icp(xy, xyps)
    @test isapprox(NeXLParticle.icperror(xy, res), 0.13382f0, atol=0.001f0)
end

@testset "ICP64_1" begin
    r, npts, ty = MersenneTwister(0xBADF00D), 100, Float64
    xy = SVector.(zip(rand(r, ty, npts), rand(r, ty, npts)))
    @test eltype(eltype(xy)) == ty
    # Muddle the points up
    xyp = AffineMap(Angle2d{ty}(deg2rad(1.2f0)), [0.1f0, -0.1f0]).(xy) .+ 0.02f0 * SVector.(collect(zip(0.5f0 .- rand(r, ty, npts), 0.5f0 .- rand(r, ty, npts))))
    xyps = shuffle(r, xyp)[1:npts÷2] # Pick only fifty points to match
    @test eltype(eltype(xyps)) == ty
    @test isapprox(NeXLParticle.icperror(xy, xyps), 2.5320, atol=0.0001)
    @test isapprox(NeXLParticle.icperror(xy[1:2:100], xyps), 3.032293, atol=0.0001)
    res=icp(xy, xyps)
    @test isapprox(NeXLParticle.icperror(xy, res), 0.075434, atol=0.0001)
end

@testset "ICP64_2" begin
    r, npts, ty = MersenneTwister(0xBADADDA), 100, Float64
    xy = SVector.(zip(rand(r, ty, npts), rand(r, ty, npts)))
    @test eltype(eltype(xy)) == ty
    # Muddle the points up
    xyp = AffineMap(Angle2d{ty}(deg2rad(1.2f0)), [0.1f0, -0.1f0]).(xy) .+ 0.02f0 * SVector.(collect(zip(0.5f0 .- rand(r, ty, npts), 0.5f0 .- rand(r, ty, npts))))
    xyps = shuffle(r, xyp)[1:npts÷2] # Pick only fifty points to match
    @test eltype(eltype(xyps)) == ty
    @test isapprox(NeXLParticle.icperror(xy, xyps), 3.458189, atol=0.0001)
    @test isapprox(NeXLParticle.icperror(xy[1:2:100], xyps), 4.235188, atol=0.0001)
    res=icp(xy, xyps)
    @test isapprox(NeXLParticle.icperror(xy, res), 0.172820, atol=0.0001)
end

max_stars_repo_name: NicholasWMRitchie/NeXLParticle

max_issues_repo_name: NicholasWMRitchie/NeXLParticle

max_forks_repo_name: NicholasWMRitchie/NeXLParticle


text:
using ArnoldiMethod

function test_fd_derivatives(a, b, N, ::Type{B}, f::Function, g::Function, h::Function) where {B<:AbstractQuasiMatrix}
    L = b-a
    Δx = L/(N+1)
    j = (1:N) .+ round(Int, a/Δx)

    R = if B == RadialDifferences
        RadialDifferences(N, Δx, 1.0, 0.0)
    else
        B(j, Δx)
    end

    D = Derivative(axes(R,1))

    ∇ = R'D*R
    ∇² = R'D'D*R

    r = axes(R,1)

    fv = R \ f.(r)
    gv = similar(fv)
    hv = similar(fv)

    mul!(gv, ∇, fv)
    mul!(hv, ∇², fv)

    δg = gv - R \ g.(r)
    δh = hv - R \ h.(r)

    r,fv,gv,hv,δg,δh,step(R)
end

function error_slope(loghs,ϵ)
    # To avoid the effect of round-off errors on the order
    # estimation.
    i = argmin(abs.(ϵ)) - 1

    ([loghs[1:i] ones(i)] \ log10.(abs.(ϵ[1:i])))[1]
end

function compute_derivative_errors(a, b, Ns, ::Type{B}, f::Function, g::Function, h::Function) where {B<:AbstractQuasiMatrix}
    errors = map(Ns) do N
        r,fv,gv,hv,δg,δh,Δx = test_fd_derivatives(a, b, N, B, f, g, h)

        [norm(δg)*Δx norm(δh)*Δx]
    end |> e -> vcat(e...)

    ϵg = errors[:,1]
    ϵh = errors[:,2]

    loghs = log10.(1.0 ./ Ns)
    pg = error_slope(loghs, ϵg)
    ph = error_slope(loghs, ϵh)

    ϵg,ϵh,pg,ph
end

function derivative_test_functions(d)
    a,b = √d*[-1,1]

    # The functions need to vanish at the boundaries, for the
    # derivative approximation to be valid (only Dirichlet0 boundary
    # conditions implemented).
    f = x -> exp(-1/(d-x^2))
    g = x -> -2*exp(-1/(d-x^2))*x/((d-x^2)^2)
    h = x -> -2*exp(-1/(d-x^2))*(d^2 + 2*(d-1)*x^2-3x^4)/((d-x^2)^4)

    f,g,h,a,b
end

function norm_rot!(v)
    normalize!(v)
    vc = v.args[2]
    vc[:] *= sign(vc[1])
    v
end

struct ShiftInvert{M}
    A⁻¹::M
end

Base.size(S::ShiftInvert, args...) = size(S.A⁻¹, args...)
Base.eltype(S::ShiftInvert) = eltype(S.A⁻¹)

LinearAlgebra.mul!(y, S::ShiftInvert, x) =
    ldiv!(y, S.A⁻¹, x)

function diagonalize_fd_hamiltonian(H, R::B, nev, σ; method=:arnoldi_shift_invert) where {B<:AbstractQuasiMatrix}
    A,target = if method == :arnoldi_shift_invert
        ShiftInvert(factorize(H - σ*I)),LR()
    else
        H, SR()
    end
    schur,history = partialschur(A, nev=nev, which=target)
    println(history)

    ϕ = [norm_rot!(R*schur.Q[:,j]) for j = 1:nev]

    θ = diag(schur.R)
    λ = if method == :arnoldi_shift_invert
        σ .+ inv.(θ)
    else
        θ
    end

    λ,ϕ
end

function get_kinetic_operator(R::B) where {B<:AbstractQuasiMatrix}
    D = Derivative(Base.axes(R,1))
    apply(*, R', D', D, R) / -2
end

function test_fd_particle_in_a_box(::Type{B}, N, L, nev; kwargs...) where {B<:AbstractQuasiMatrix}
    R = B(1:N, L/(N+1))
    display(R)
    println()

    Tm = get_kinetic_operator(R)
    λ,ϕ = diagonalize_fd_hamiltonian(Tm, R, nev, 0.0; kwargs...)

    r̃ = FiniteDifferencesQuasi.locs(R)

    n = 1:nev
    δλ = λ - n.^2*π^2/(2L^2)

    # Could/should also test eigenvectors

    λ,ϕ,r̃,R,δλ
end

function test_singular_fd_scheme(::Type{B}, N, rₘₐₓ, Z, ℓ, nev; kwargs...) where {B<:AbstractQuasiMatrix}
    R = if B == RadialDifferences
        ρ = rₘₐₓ/(N+0.5)
        RadialDifferences(N, ρ, Z)
    else
        ρ = rₘₐₓ/(N+1)
        NumerovFiniteDifferences(1:N, ρ, true, Z)
    end
    display(R)
    println()

    r = axes(R,1)

    n = 1:nev
    λₐ = -inv.(2(n .+ ℓ).^2)

    Tm = get_kinetic_operator(R) + R'*QuasiDiagonal(ℓ*(ℓ+1)./2r.^2)*R
    V = R'*QuasiDiagonal(-inv.(r))*R
    H = Tm + V
    λ,ϕ = diagonalize_fd_hamiltonian(H, R, nev, 1.1λₐ[1]; kwargs...)

    r̃ = FiniteDifferencesQuasi.locs(R)

    δλ = λ - λₐ

    # Could/should also test eigenvectors

    λ,ϕ,r̃,R,δλ
end

function compute_diagonalization_errors(f::Function, ::Type{B}, Ns, args...; kwargs...) where B
    errors = map(Ns) do N
        println("N = $N")
        t = time()
        λ,ϕ,r̃,R,δλ = f(B, N, args...; kwargs...)
        elapsed = time()-t
        vcat(δλ,elapsed)'
    end |> e -> vcat(e...)

    errors,elapsed = errors[:,1:end-1],errors[:,end]

    loghs = log10.(inv.(Ns))
    slopes = map(1:size(errors,2)) do j
        error_slope(loghs, errors[:,j])
    end

    errors,slopes,elapsed
end


text:
# Submitted by Gustorn

const mnames = ["A", "B", "C", "D"]
const wnames = ["E", "F", "G", "H"]

const Preferences = Dict{String,Vector{String}}
const Pairs = Dict{String,String}

# Returns a name => preference list dictionary, in decreasing order of preference
function genpreferences(mannames::Vector{String}, womannames::Vector{String})
    men   = Dict(map(m -> (m, shuffle(womannames)), mannames))
    women = Dict(map(w -> (w, shuffle(mannames)), womannames))
    return men, women
end

# Returns if `person` prefers the `first` candidate over the `second` one.
# This translates to `first` appearing *sooner* in the preference list
prefers(prefs, person, first, second) =
    findfirst(m -> m == first, prefs[person]) <
    findfirst(m -> m == second, prefs[person])

isfree(person, pairs) = !haskey(pairs, person)

function galeshapely(men::Preferences, women::Preferences)
    mentowomen = Dict{String,String}()
    womentomen = Dict{String,String}()
    while true
        bachelors = [m for m in keys(men) if isfree(m, mentowomen)]
        if length(bachelors) == 0
            return mentowomen, womentomen
        end

        for bachelor in bachelors
            for candidate in men[bachelor]
                if isfree(candidate, womentomen)
                    mentowomen[bachelor] = candidate
                    womentomen[candidate] = bachelor
                    break
                elseif prefers(women, candidate, bachelor, womentomen[candidate])
                    delete!(mentowomen, womentomen[candidate])
                    mentowomen[bachelor] = candidate
                    womentomen[candidate] = bachelor
                    break
                end
            end
        end
    end
end

function isstable(men::Preferences, women::Preferences, mentowomen::Pairs, womentoman::Pairs)
    for (husband, wife) in mentowomen
        for candidate in men[husband]
            if candidate != wife &&
               prefers(men, husband, candidate, wife) &&
               prefers(women, candidate, husband, womentoman[candidate])
                return false
            end
        end
    end
    return true
end

function main()
    men, women = genpreferences(mnames, wnames)
    mentowomen, womentomen = galeshapely(men, women)
    println(mentowomen)
    println(isstable(men, women, mentowomen, womentomen) ? "Stable" : "Unstable")
end

main()


text:
using AccretionFormulae
using Plots

function r_hor(a, M)
    (M +  √(M^2 - (a*M)^2))/M
end

M = 10*1.99e30
spin_vals_pos = LinRange(0, 1, 1000)
spin_vals_neg = LinRange(-1, 0 , 1000)
ISCO_pos_vals = AccretionFormulae.r_isco.(spin_vals_pos, M) ./ M
ISCO_neg_vals = reverse(AccretionFormulae.r_isco.(spin_vals_neg, M)) ./ M
EH_pos_vals = r_hor.(spin_vals_pos, M)

plt = plot(
            spin_vals_pos,
            ISCO_pos_vals,
            legend=:topleft,
            label="\$r_{ISCO} \\; \\textrm{Prograde}\$",
            xlabel="a (M)",
            ylabel="R (M)",
            grid=false,
            framestyle=:box
            )
plot!(spin_vals_pos, ISCO_neg_vals, label="\$r_{ISCO} \\; \\textrm{Retrograde}\$")
plot!(spin_vals_pos, EH_pos_vals, label="\$r_{s}\$")
png(plt, "ISCO_EH_spin_comp.png")


text:
#=
Basic 3D reconstruction
=#

function plot_arrays(neuron::Neuron)

    x=[zeros(Float64,0) for i=1:length(neuron.secs)]
    y=[zeros(Float64,0) for i=1:length(neuron.secs)]

    for i=1:length(neuron.secs), j=1:length(neuron.secs[i].pt3d)
        push!(x[i],neuron.secs[i].pt3d[j].x)
        push!(y[i],neuron.secs[i].pt3d[j].y)
    end

    [hcat(x[i],y[i]) for i=1:length(x)]
end

function getxyz(n::Neuron)

    x=Array(Float64,0)
    y=Array(Float64,0)
    z=Array(Float64,0)

    secind=[Array(Int64,0) for i=1:length(n.secs)]
    count=1

    @inbounds for i=1:length(n.secs)
        for j in n.secs[i].pnode
            if n.nodes[j].internal == true
                (x1,y1,z1)=getxyz(n.secs[i].pt3d[n.nodes[j].pt3d])
                push!(x,x1[1])
                push!(y,y1[1])
                push!(z,z1[1])
                push!(x,x1[end])
                push!(y,y1[end])
                push!(z,z1[end])
                push!(secind[i],count)
                push!(secind[i],count+1)
                count+=2
            end
        end
    end
    (hcat(x,y,z),secind)
end

function getxyz(p::Array{JNeuron.Pt3d,1})

    x=Array(Float64,0)
    y=Array(Float64,0)
    z=Array(Float64,0)

    @inbounds for i=1:length(p)

        push!(x,p[i].x)
        push!(y,p[i].y)
        push!(z,p[i].z)

    end
    (x,y,z)
end

function copyn(xyz::Array{Float64,2})

    xyz1=zeros(Float64,size(xyz))
    resetn!(xyz,xyz1)
    xyz1
end

function resetn!(xyz::Array{Float64,2},xyz1::Array{Float64,2})
    @inbounds @simd for i=1:length(xyz)
        xyz1[i]=xyz[i]
    end
    nothing
end

function translate3d!(neuron::Neuron,x::Float64,y::Float64,z::Float64)

    for i=1:length(neuron.secs), j=1:length(neuron.secs[i].pt3d)
            neuron.secs[i].pt3d[j]=add(neuron.secs[i].pt3d[j],x,y,z)
    end
    nothing
end

function translate3d!(xyz::Array{Float64,2},x::Float64,y::Float64,z::Float64)

    @inbounds @simd for i=1:size(xyz,1)
        xyz[i,1]+=x
        xyz[i,2]+=y
        xyz[i,3]+=z
    end
    nothing
end

function rotate3d!(neuron::Neuron,theta::Float64,myaxis::Int64)

    rotmat=zeros(Float64,3,3)

    if myaxis==1

        rotmat[1,1]=1.0
        rotmat[2,2]=cos(theta)
        rotmat[2,3]=-sin(theta)
        rotmat[3,2]=sin(theta)
        rotmat[3,3]=cos(theta)

    elseif myaxis==2

        rotmat[1,1]=cos(theta)
        rotmat[1,3]=sin(theta)
        rotmat[2,2]=1.0
        rotmat[3,1]=-sin(theta)
        rotmat[3,3]=cos(theta)
    else

        rotmat[1,1]=cos(theta)
        rotmat[1,2]=-sin(theta)
        rotmat[2,1]=sin(theta)
        rotmat[2,2]=cos(theta)
        rotmat[3,3]=1.0
    end

    for i=1:length(neuron.secs)
        for j=1:length(neuron.secs[i].pt3d)

            xyz=rotmat*[neuron.secs[i].pt3d[j].x,neuron.secs[i].pt3d[j].y,neuron.secs[i].pt3d[j].z]

            neuron.secs[i].pt3d[j]=Pt3d(neuron.secs[i].pt3d[j],xyz)
        end
    end
end

function randomize_shape!(neuron::Neuron)

    for i=1:length(neuron.secs)

        if (neuron.secs[i].parent!=length(neuron.secs))&&(neuron.secs[i].parent!=0)

            rot_seg(neuron,neuron.secs[i].parent,i)
        end
    end
    nothing
end

function randomize_shape!(n::Neuron,xyz::Array{Float64,2},inds::Array{Array{Int64,1},1})

    @inbounds for i=1:length(n.secs)
        if (n.secs[i].parent!=length(n.secs))&&(n.secs[i].parent!=0)
            rot_seg(n,n.secs[i].parent,i,xyz,inds)
        end
    end
    nothing
end

function rot_3d(a,b,c,u,v,w,x,y,z,dcos,dsin)

    xnew=(a*(v*v+w*w)-u*(b*v + c*w - u*x - v*y - w*z))*(1-dcos)+x*dcos+(-c*v + b*w - w*y + v*z)*dsin
    ynew=(b*(u*u+w*w)-v*(a*u + c*w - u*x - v*y - w*z))*(1-dcos)+y*dcos+(c*u - a*w + w*x - u*z)*dsin
    znew=(c*(u*u+v*v)-w*(a*u + b*v - u*x - v*y - w*z))*(1-dcos)+z*dcos+(-b*u + a*v - v*x + u*y)*dsin

    (xnew,ynew,znew)
end

function rot_seg(neuron::Neuron,pind::Int64,cind::Int64)

    #last point of previous segment
    a=neuron.secs[pind].pt3d[end].x
    b=neuron.secs[pind].pt3d[end].y
    c=neuron.secs[pind].pt3d[end].z

    (uvw,mag)=pt3d_vec(neuron.secs[pind].pt3d[end],neuron.secs[pind].pt3d[end-1])

    theta=randn()*10

    dcos=cosd(theta)
    dsin=sind(theta)

    rot_child(neuron,cind,a,b,c,uvw,dcos,dsin)
    nothing
end

function rot_seg(n::Neuron,pind::Int64,cind::Int64,xyz::Array{Float64,2},inds::Array{Array{Int64,1},1})

    #last point of previous segment
    @inbounds a=xyz[inds[pind][end],1]
    @inbounds b=xyz[inds[pind][end],2]
    @inbounds c=xyz[inds[pind][end],3]

    @inbounds (vx,vy,vz,mag)=pt3d_vec(xyz[inds[pind][end-1],1],xyz[inds[pind][end-1],2],xyz[inds[pind][end-1],3],a,b,c)

    theta=randn()*10

    dcos=cosd(theta)
    dsin=sind(theta)

    rot_child(n,cind,a,b,c,vx,vy,vz,dcos,dsin,xyz,inds)
    nothing
end

function rot_child(neuron::Neuron,cind::Int64,a::Float64,b::Float64,c::Float64,uvw::Array{Float64,1},dcos::Float64,dsin::Float64)

    for i in neuron.secs[cind].child

        rot_child(neuron,i,a,b,c,uvw,dcos,dsin)
    end

    @inbounds for j=1:length(neuron.secs[cind].pt3d)

        x=neuron.secs[cind].pt3d[j].x
        y=neuron.secs[cind].pt3d[j].y
        z=neuron.secs[cind].pt3d[j].z

        xyz=rot_3d(a,b,c,uvw[1],uvw[2],uvw[3],x,y,z,dcos,dsin)

        neuron.secs[cind].pt3d[j]=Pt3d(xyz[1],xyz[2],xyz[3],neuron.secs[cind].pt3d[j].d,neuron.secs[cind].pt3d[j].arc)

    end
    nothing
end

function rot_child(n::Neuron,cind::Int64,a::Float64,b::Float64,c::Float64,vx::Float64,vy::Float64,vz::Float64,dcos::Float64,dsin::Float64,xyz::Array{Float64,2},inds::Array{Array{Int64,1},1})

    for i in n.secs[cind].child
        rot_child(n,i,a,b,c,vx,vy,vz,dcos,dsin,xyz,inds)
    end

    @inbounds @fastmath for j=1:length(inds[cind])

        x=xyz[inds[cind][j],1]
        y=xyz[inds[cind][j],2]
        z=xyz[inds[cind][j],3]

        xyz1=rot_3d(a,b,c,vx,vy,vz,x,y,z,dcos,dsin)

        xyz[inds[cind][j],1]=xyz1[1]
        xyz[inds[cind][j],2]=xyz1[2]
        xyz[inds[cind][j],3]=xyz1[3]

    end
    nothing
end

function findcind(n)
    cinds=Array(Int64,0)
    @inbounds for i=1:length(n.secs)
        for j in n.secs[i].pnode
            if n.nodes[j].internal == true
                push!(cinds,j)
            end
        end
    end
    cinds
end

max_stars_repo_name: paulmthompson/JNeuron

max_issues_repo_name: paulmthompson/JNeuron

max_forks_repo_name: paulmthompson/JNeuron


text:
"""
    get_highlights(target::String; concatenate::Bool=true) -> Vector{String}

Get values from the `Highlights` column if a CSV file is passed,
or get the highlights from a PDF file if it is passed, or get the highlights from all
PDFs found recursively if a directory is passed.

# Arguments
- `target::String`: a CSV file, or a PDF file, or a directory with PDF files

# Keywords
- `concatenate::Bool=true`: $(CONCATENATE_BOTH_KEYWORD)

# Returns
- `Vector{String}`: the highlights

# Throws
- [`NotCSVorPDForDir`](@ref): $(NOT_CSV_OR_PDF_OR_DIR_EXCEPTION)
- Exceptions from: [`_get_highlights_from_CSV`](@ref), [`_get_highlights_from_PDF`](@ref)

# Example
```jldoctest; output = false
using PDFHighlights
HEADER = PDFHighlights.Internal.Constants.HEADER

path_to_pdf_dir = joinpath(pathof(PDFHighlights) |> dirname |> dirname, "test", "pdf")
path_to_pdf = joinpath(path_to_pdf_dir, "TestPDF.pdf")

(get_highlights(path_to_pdf_dir) == get_highlights(path_to_pdf) == String[
    "Highlight 1",
    "Highlight 2 Highlight 3",
    "Highlight 4",
    "Highhighlight 5",
    "6th Highhigh light-",
    "High light 7",
    "8th Highlight-",
]) |> println

_file, io = mktemp()
println(io, HEADER, '\\n', "The world didn't stop spinning,,,,1")
flush(io)
file = _file * ".csv"
mv(_file, file)

get_highlights(file) == ["The world didn't stop spinning"]

# output

true
true
```
"""
function get_highlights(target::String; concatenate::Bool=true)::Vector{String}

    if isdir(target) || endswith(target, ".pdf")
        _get_highlights_from_PDF(target; concatenate)
    elseif endswith(target, ".csv")
        _get_highlights_from_CSV(target)
    else
        throw(NotCSVorPDForDir(target))
    end

end


text:
function templateTestTmapreduce1(arr::AbstractArray)::Bool
	init = randInit();
	return mapreduce(x->2x, +, arr; init = init) == tmapreduce(x->2x, +, arr; init = init);
end

function testTmapreduce1Stable()::Bool
	return templateTestTmapreduce1(randArray());
end

function testTmapreduce1Unstable()::Bool
	arr = randStrArray();
	str = lipsum();
	init = randstring(rand(1:1000));
	return mapreduce(x->str[1:x], *, arr; init = init) == tmapreduce(x->str[1:x], *, arr; init = init);
end

function testTmapreduce1Offset()::Bool
	return templateTestTmapreduce1(offsetArray());
end

function testTmapreduce1Empty()::Bool
	return templateTestTmapreduce1([]);
end

#####

function templateTestTmapreduce1Multiple(arr1::AbstractArray, arr2::AbstractArray)::Bool
	init = randInit();
	return mapreduce((x, y)->2x + 3y, +, arr1, arr2; init = init) == tmapreduce((x, y)->2x + 3y, +, arr1, arr2; init = init);
end

function testTmapreduce1StableMultiple()::Bool
	return templateTestTmapreduce1Multiple(rand2Arrays()...);
end

function testTmapreduce1UnstableMultiple()::Bool
	arr1, arr2 = rand2StrArrays();
	str = lipsum();
	init = randstring(rand(1:1000));
	return mapreduce((x, y)->str[1:x + y], *, arr1, arr2; init = init) == tmapreduce((x, y)->str[1:x + y], *, arr1, arr2; init = init);
end

function testTmapreduce1MultipleOffset()::Bool
	return templateTestTmapreduce1Multiple(offsetArray(), offsetArray());
end

function testTmapreduce1MultipleEmpty()::Bool
	return templateTestTmapreduce1Multiple([], []);
end

#####

function templateTestTmapreduce2(arr::AbstractArray)::Bool
	return mapreduce(x->2x, +, arr) == tmapreduce(x->2x, +, arr);
end

function testTmapreduce2Stable()::Bool
	return templateTestTmapreduce2(randArray());
end

function testTmapreduce2Unstable()::Bool
	arr = randStrArray();
	str = lipsum();
	return mapreduce(x->str[1:x], *, arr) == tmapreduce(x->str[1:x], *, arr);
end

function testTmapreduce2Offset()::Bool
	return templateTestTmapreduce2(offsetArray());
end

function testTmapreduce2Empty()
	tmapreduce(x->2x, (x, y)->x, []);
end

#####

function templateTestTmapreduce2Multiple(arr1::AbstractArray, arr2::AbstractArray)::Bool
	return mapreduce((x, y)->2x + 3y, +, arr1, arr2) == tmapreduce((x, y)->2x + 3y, +, arr1, arr2);
end

function testTmapreduce2StableMultiple()::Bool
	return templateTestTmapreduce2Multiple(rand2Arrays()...);
end

function testTmapreduce2UnstableMultiple()::Bool
	arr1, arr2 = rand2StrArrays();
	str = lipsum();
	return mapreduce((x, y)->str[1:x + y], *, arr1, arr2) == tmapreduce((x, y)->str[1:x + y], *, arr1, arr2);
end

function testTmapreduce2MultipleOffset()::Bool
	return templateTestTmapreduce2Multiple(offsetArray(), offsetArray());
end

function testTmapreduce2MultipleEmpty()
	tmapreduce((x, y)->2x + 3y, (x, y)->x, [], []);
end

#####

function templateTestMaptreduce1(arr::AbstractArray)::Bool
	init = randInit();
	return mapreduce(x->2x, +, arr; init = init) == maptreduce(x->2x, +, arr; init = init);
end

function testMaptreduce1()::Bool
	return templateTestMaptreduce1(randArray());
end

function testMaptreduce1Offset()::Bool
	return templateTestMaptreduce1(offsetArray());
end

function testMaptreduce1Empty()::Bool
	return templateTestMaptreduce1([]);
end

#####

function templateTestMaptreduce1Multiple(arr1::AbstractArray, arr2::AbstractArray)::Bool
	init = randInit();
	return mapreduce((x, y)->2x + 3y, +, arr1, arr2; init = init) == maptreduce((x, y)->2x + 3y, +, arr1, arr2; init = init);
end

function testMaptreduce1Multiple()::Bool
	return templateTestMaptreduce1Multiple(rand2Arrays()...);
end

function testMaptreduce1MultipleOffset()::Bool
	return templateTestMaptreduce1Multiple(offsetArray(), offsetArray());
end

function testMaptreduce1MultipleEmpty()::Bool
	return templateTestMaptreduce1Multiple([], []);
end

#####

function templateTestMaptreduce2(arr::AbstractArray)::Bool
	return mapreduce(x->2x, +, arr) == maptreduce(x->2x, +, arr);
end

function testMaptreduce2()::Bool
	return templateTestMaptreduce2(randArray());
end

function testMaptreduce2Offset()::Bool
	return templateTestMaptreduce2(offsetArray());
end

function testMaptreduce2Empty()
	maptreduce(x->2x, +, []);
end

#####

function templateTestMaptreduce2Multiple(arr1::AbstractArray, arr2::AbstractArray)::Bool
	return mapreduce((x, y)->2x + 3y, +, arr1, arr2) == maptreduce((x, y)->2x + 3y, +, arr1, arr2);
end

function testMaptreduce2Multiple()::Bool
	return templateTestMaptreduce2Multiple(rand2Arrays()...);
end

function testMaptreduce2MultipleOffset()::Bool
	return templateTestMaptreduce2Multiple(offsetArray(), offsetArray());
end

function testMaptreduce2MultipleEmpty()
	maptreduce((x, y)->2x + 3y, +, [], []);
end

#####

function templateTestTmaptreduce1(arr::AbstractArray)::Bool
	init = randInit();
	return mapreduce(x->2x, +, arr; init = init) == tmaptreduce(x->2x, +, arr; init = init);
end

function testTmaptreduce1()::Bool
	return templateTestTmaptreduce1(randArray());
end

function testTmaptreduce1Offset()::Bool
	return templateTestTmaptreduce1(offsetArray());
end

function testTmaptreduce1Empty()::Bool
	return templateTestTmaptreduce1([]);
end

###

function templateTestTmaptreduce1Multiple(arr1::AbstractArray, arr2::AbstractArray)::Bool
	init = randInit();
	return mapreduce((x, y)->2x + 3y, +, arr1, arr2; init = init) == tmaptreduce((x, y)->2x + 3y, +, arr1, arr2; init = init);
end

function testTmaptreduce1Multiple()::Bool
	return templateTestTmaptreduce1Multiple(rand2Arrays()...);
end

function testTmaptreduce1MultipleOffset()::Bool
	return templateTestTmaptreduce1Multiple(offsetArray(), offsetArray());
end

function testTmaptreduce1MultipleEmpty()::Bool
	return templateTestTmaptreduce1Multiple([], []);
end

#####

function templateTestTmaptreduce2(arr::AbstractArray)::Bool
	return mapreduce(x->2x, +, arr) == tmaptreduce(x->2x, +, arr);
end

function testTmaptreduce2()::Bool
	return templateTestTmaptreduce2(randArray());
end

function testTmaptreduce2Offset()::Bool
	return templateTestTmaptreduce2(offsetArray());
end

function testTmaptreduce2Empty()
	tmaptreduce(x->2x, +, []);
end

###

function templateTestTmaptreduce2Multiple(arr1::AbstractArray, arr2::AbstractArray)::Bool
	return mapreduce((x, y)->2x + 3y, +, arr1, arr2) == tmaptreduce((x, y)->2x + 3y, +, arr1, arr2);
end

function testTmaptreduce2Multiple()::Bool
	return templateTestTmaptreduce2Multiple(rand2Arrays()...);
end

function testTmaptreduce2MultipleOffset()::Bool
	return templateTestTmaptreduce2Multiple(offsetArray(), offsetArray());
end

function testTmaptreduce2MultipleEmpty()
	tmaptreduce((x, y)->2x + 3y, +, [], []);
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

module ComputeManagementClient

using Random
using Dates
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_AvailabilitySetsApi.jl")
include("api_ComputeOperationsApi.jl")
include("api_DedicatedHostApi.jl")
include("api_DedicatedHostGroupsApi.jl")
include("api_DedicatedHostsApi.jl")
include("api_ImagesApi.jl")
include("api_LogAnalyticsApi.jl")
include("api_ProximityPlacementGroupsApi.jl")
include("api_SshPublicKeysApi.jl")
include("api_UsageApi.jl")
include("api_VirtualMachineExtensionImagesApi.jl")
include("api_VirtualMachineExtensionsApi.jl")
include("api_VirtualMachineImagesApi.jl")
include("api_VirtualMachineScaleSetExtensionsApi.jl")
include("api_VirtualMachineScaleSetRollingUpgradesApi.jl")
include("api_VirtualMachineScaleSetVMExtensionsApi.jl")
include("api_VirtualMachineScaleSetVMsApi.jl")
include("api_VirtualMachineScaleSetsApi.jl")
include("api_VirtualMachineSizesApi.jl")
include("api_VirtualMachinesApi.jl")

# export models
export convert, AdditionalCapabilities
export convert, AdditionalUnattendContent
export convert, ApiEntityReference
export convert, ApiError
export convert, ApiErrorBase
export convert, AutomaticOSUpgradePolicy
export convert, AutomaticOSUpgradeProperties
export convert, AutomaticRepairsPolicy
export convert, AvailabilitySet
export convert, AvailabilitySetListResult
export convert, AvailabilitySetProperties
export convert, AvailabilitySetSkuType
export convert, AvailabilitySetUpdate
export convert, BillingProfile
export convert, BootDiagnostics
export convert, BootDiagnosticsInstanceView
export convert, Caching
export convert, CloudError
export convert, ComputeOperationListResult
export convert, ComputeOperationValue
export convert, ComputeOperationValueDisplay
export convert, CreateOption
export convert, DataDisk
export convert, DataDiskImage
export convert, DedicatedHost
export convert, DedicatedHostAllocatableVM
export convert, DedicatedHostAvailableCapacity
export convert, DedicatedHostGroup
export convert, DedicatedHostGroupInstanceView
export convert, DedicatedHostGroupListResult
export convert, DedicatedHostGroupProperties
export convert, DedicatedHostGroupUpdate
export convert, DedicatedHostInstanceView
export convert, DedicatedHostInstanceViewWithName
export convert, DedicatedHostLicenseType
export convert, DedicatedHostListResult
export convert, DedicatedHostProperties
export convert, DedicatedHostUpdate
export convert, DiagnosticsProfile
export convert, DiffDiskOption
export convert, DiffDiskPlacement
export convert, DiffDiskSettings
export convert, DisallowedConfiguration
export convert, DiskEncryptionSetParameters
export convert, DiskEncryptionSettings
export convert, DiskInstanceView
export convert, EvictionPolicy
export convert, HardwareProfile
export convert, HyperVGenerationType
export convert, Image
export convert, ImageDataDisk
export convert, ImageDisk
export convert, ImageListResult
export convert, ImageOSDisk
export convert, ImageProperties
export convert, ImageReference
export convert, ImageStorageProfile
export convert, ImageUpdate
export convert, InnerError
export convert, InstanceViewStatus
export convert, KeyVaultKeyReference
export convert, KeyVaultSecretReference
export convert, LinuxConfiguration
export convert, ListUsagesResult
export convert, LogAnalyticsInputBase
export convert, LogAnalyticsOperationResult
export convert, LogAnalyticsOutput
export convert, MaintenanceRedeployStatus
export convert, ManagedDiskParameters
export convert, NetworkInterfaceReference
export convert, NetworkInterfaceReferenceProperties
export convert, NetworkProfile
export convert, OSDisk
export convert, OSDiskImage
export convert, OSProfile
export convert, OrchestrationServiceStateInput
export convert, OrchestrationServiceSummary
export convert, PatchSettings
export convert, Plan
export convert, Priority
export convert, ProximityPlacementGroup
export convert, ProximityPlacementGroupListResult
export convert, ProximityPlacementGroupProperties
export convert, ProximityPlacementGroupUpdate
export convert, PurchasePlan
export convert, RecoveryWalkResponse
export convert, RequestRateByIntervalInput
export convert, Resource
export convert, RetrieveBootDiagnosticsDataResult
export convert, RollbackStatusInfo
export convert, RollingUpgradePolicy
export convert, RollingUpgradeProgressInfo
export convert, RollingUpgradeRunningStatus
export convert, RollingUpgradeStatusInfo
export convert, RollingUpgradeStatusInfoProperties
export convert, ScaleInPolicy
export convert, ScheduledEventsProfile
export convert, SecurityProfile
export convert, Sku
export convert, SshConfiguration
export convert, SshPublicKey
export convert, SshPublicKeyGenerateKeyPairResult
export convert, SshPublicKeyResource
export convert, SshPublicKeyResourceProperties
export convert, SshPublicKeyUpdateResource
export convert, SshPublicKeysGroupListResult
export convert, StorageAccountType
export convert, StorageProfile
export convert, SubResource
export convert, SubResourceReadOnly
export convert, SubResourceWithColocationStatus
export convert, TerminateNotificationProfile
export convert, ThrottledRequestsInput
export convert, UpdateResource
export convert, UpgradeOperationHistoricalStatusInfo
export convert, UpgradeOperationHistoricalStatusInfoProperties
export convert, UpgradeOperationHistoryStatus
export convert, UpgradePolicy
export convert, Usage
export convert, UsageName
export convert, VMScaleSetConvertToSinglePlacementGroupInput
export convert, VaultCertificate
export convert, VaultSecretGroup
export convert, VirtualHardDisk
export convert, VirtualMachine
export convert, VirtualMachineAgentInstanceView
export convert, VirtualMachineAssessPatchesResult
export convert, VirtualMachineCaptureParameters
export convert, VirtualMachineCaptureResult
export convert, VirtualMachineExtension
export convert, VirtualMachineExtensionHandlerInstanceView
export convert, VirtualMachineExtensionImage
export convert, VirtualMachineExtensionImageProperties
export convert, VirtualMachineExtensionInstanceView
export convert, VirtualMachineExtensionProperties
export convert, VirtualMachineExtensionUpdate
export convert, VirtualMachineExtensionUpdateProperties
export convert, VirtualMachineExtensionsListResult
export convert, VirtualMachineHealthStatus
export convert, VirtualMachineIdentity
export convert, VirtualMachineIdentityUserAssignedIdentities
export convert, VirtualMachineImage
export convert, VirtualMachineImageProperties
export convert, VirtualMachineImageResource
export convert, VirtualMachineInstanceView
export convert, VirtualMachineListResult
export convert, VirtualMachineProperties
export convert, VirtualMachineReimageParameters
export convert, VirtualMachineScaleSet
export convert, VirtualMachineScaleSetDataDisk
export convert, VirtualMachineScaleSetExtension
export convert, VirtualMachineScaleSetExtensionListResult
export convert, VirtualMachineScaleSetExtensionProfile
export convert, VirtualMachineScaleSetExtensionProperties
export convert, VirtualMachineScaleSetExtensionUpdate
export convert, VirtualMachineScaleSetIPConfiguration
export convert, VirtualMachineScaleSetIPConfigurationProperties
export convert, VirtualMachineScaleSetIdentity
export convert, VirtualMachineScaleSetInstanceView
export convert, VirtualMachineScaleSetInstanceViewStatusesSummary
export convert, VirtualMachineScaleSetIpTag
export convert, VirtualMachineScaleSetListOSUpgradeHistory
export convert, VirtualMachineScaleSetListResult
export convert, VirtualMachineScaleSetListSkusResult
export convert, VirtualMachineScaleSetListWithLinkResult
export convert, VirtualMachineScaleSetManagedDiskParameters
export convert, VirtualMachineScaleSetNetworkConfiguration
export convert, VirtualMachineScaleSetNetworkConfigurationDnsSettings
export convert, VirtualMachineScaleSetNetworkConfigurationProperties
export convert, VirtualMachineScaleSetNetworkProfile
export convert, VirtualMachineScaleSetOSDisk
export convert, VirtualMachineScaleSetOSProfile
export convert, VirtualMachineScaleSetProperties
export convert, VirtualMachineScaleSetPublicIPAddressConfiguration
export convert, VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings
export convert, VirtualMachineScaleSetPublicIPAddressConfigurationProperties
export convert, VirtualMachineScaleSetReimageParameters
export convert, VirtualMachineScaleSetSku
export convert, VirtualMachineScaleSetSkuCapacity
export convert, VirtualMachineScaleSetStorageProfile
export convert, VirtualMachineScaleSetUpdate
export convert, VirtualMachineScaleSetUpdateIPConfiguration
export convert, VirtualMachineScaleSetUpdateIPConfigurationProperties
export convert, VirtualMachineScaleSetUpdateNetworkConfiguration
export convert, VirtualMachineScaleSetUpdateNetworkConfigurationProperties
export convert, VirtualMachineScaleSetUpdateNetworkProfile
export convert, VirtualMachineScaleSetUpdateOSDisk
export convert, VirtualMachineScaleSetUpdateOSProfile
export convert, VirtualMachineScaleSetUpdateProperties
export convert, VirtualMachineScaleSetUpdatePublicIPAddressConfiguration
export convert, VirtualMachineScaleSetUpdatePublicIPAddressConfigurationProperties
export convert, VirtualMachineScaleSetUpdateStorageProfile
export convert, VirtualMachineScaleSetUpdateVMProfile
export convert, VirtualMachineScaleSetVM
export convert, VirtualMachineScaleSetVMExtensionsSummary
export convert, VirtualMachineScaleSetVMInstanceIDs
export convert, VirtualMachineScaleSetVMInstanceRequiredIDs
export convert, VirtualMachineScaleSetVMInstanceView
export convert, VirtualMachineScaleSetVMListResult
export convert, VirtualMachineScaleSetVMNetworkProfileConfiguration
export convert, VirtualMachineScaleSetVMProfile
export convert, VirtualMachineScaleSetVMProperties
export convert, VirtualMachineScaleSetVMProtectionPolicy
export convert, VirtualMachineScaleSetVMReimageParameters
export convert, VirtualMachineSize
export convert, VirtualMachineSizeListResult
export convert, VirtualMachineSoftwarePatchProperties
export convert, VirtualMachineStatusCodeCount
export convert, VirtualMachineUpdate
export convert, WinRMConfiguration
export convert, WinRMListener
export convert, WindowsConfiguration

# export operations
export convert, AvailabilitySetsApi, ComputeOperationsApi, DedicatedHostApi, DedicatedHostGroupsApi, DedicatedHostsApi, ImagesApi, LogAnalyticsApi, ProximityPlacementGroupsApi, SshPublicKeysApi, UsageApi, VirtualMachineExtensionImagesApi, VirtualMachineExtensionsApi, VirtualMachineImagesApi, VirtualMachineScaleSetExtensionsApi, VirtualMachineScaleSetRollingUpgradesApi, VirtualMachineScaleSetVMExtensionsApi, VirtualMachineScaleSetVMsApi, VirtualMachineScaleSetsApi, VirtualMachineSizesApi, VirtualMachinesApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
export lmlovo, gnlslovo

"""
    lmlovo(model::Function [, θ::Vector{Float64} = zeros(n)], data::Array{Float64, 2},
           n::Int, p::Int [; kwargs...])

    lmlovo(model::Function, gmodel!::Function [, θ::Vector{Float64} = zeros(n)],
           data::Array{Float64,2}, n::Int, p::Int [; MAXITER::Int=200,
           ε::Float64=10.0^-4])

Fit the `n`-parameter model `model` to the data given by matrix
`data`. The strategy is based on the LOVO function, which means that
only `p` (0 < `p` <= rows of `data`) points are trusted. The
Levenberg-Marquardt algorithm is implemented in this version.

Matriz `data` is the data to be fit. This matrix should be in the form

    x11 x12 ... x1N y1
    x21 x22 ... x2N y2
    :

where `N` is the dimension of the argument of the model
(i.e. dimension of `x`).

If `θ` is provided, then it is used as the starting point.

The signature of function `model` should be given by

    model(x::Union{Vector{Float64}, SubArray}, θ::Vector{Float64})

where `x` are the variables and `θ` is a `n`-dimensional vector of
parameters. If the gradient of the model `gmodel!`

    gmodel! = (g::SubArray, x::Union{Vector{Float64}, SubArray},
               θ::Vector{Float64})

is not provided, then the function ForwardDiff.gradient! is called to
compute it.  **Note** that this choice has an impact in the
computational performance of the algorithm. In addition, if
`ForwardDiff.jl` is being used, then one **MUST** remove the signature
of vector `θ` from function `model`.

The optional arguments are

  - `MAXITER`: maximum number of iterations
  - `ε`: tolerance for the gradient of the function

Returns a [`RAFFOutput`](@ref) object.

"""
function lmlovo(model::Function, gmodel!::Function, θ::Vector{Float64},
                data::Array{Float64,2}, n::Int, p::Int;
                MAXITER::Int=400, ε::Float64=1.0e-4)

    @assert(n > 0, "Dimension should be positive.")
    @assert(p >= 0, "Trusted points should be nonnegative.")

    npun, = size(data)

    with_logger(lm_logger.x) do

        @debug("Size of data matrix ", size(data))

    end

    # Counters for calls to F and its Jacobian
    nf = 0

    nj = 0

    (p == 0) && return RAFFOutput(1, θ, 0, p, 0.0, nf, nj, [1:npun;])

    # Main function - the LOVO function
    LovoFun = let

        npun_::Int = npun

        ind::Vector{Int} = Vector{Int}(undef, npun_)

        F::Vector{Float64} = Vector{Float64}(undef, npun_)

        p_::Int = p

        # Return a ordered set index and lovo value
        (θ) -> begin

            nf += 1

            @views for i = 1:npun_
                F[i] = (model(data[i,1:(end - 1)], θ) - data[i, end])^2
            end

            indF, orderedF = sort_fun!(F, ind, p_)

            return indF, sum(orderedF)
        end

    end

    # Residue and Jacobian of residue
    val_res::Vector{Float64} = Vector{Float64}(undef, p)

    jac_res::Array{Float64, 2} = Array{Float64}(undef, p, n)

    # This function returns the residue and Jacobian of residue
    ResFun! = let

        data_    = data
        model_   = model
        gmodel!_ = gmodel!

        (θ::Vector{Float64}, ind, r::Vector{Float64},
         rJ::Array{Float64, 2}) -> begin

             nj += 1

             residual_fg!(model_, gmodel!_, data_, θ, ind, r, rJ)

         end

    end

    # Levenberg-Marquardt algorithm
    # -----------------------------

    Id = Matrix(1.0I, n, n)

    # Status = 1 means success
    status = 1

    # Parameters

    λ_up      = 2.0
    λ_down    = 2.0
    λ         = 1.0
    maxoutind = min(p, 5)

    # Allocation
    θnew = Vector{Float64}(undef, n)
    d = Vector{Float64}(undef, n)
    y = Vector{Float64}(undef, n)
    G = Array{Float64, 2}(undef, n, n)
    grad_lovo = Vector{Float64}(undef, n)

    # Initial steps

    ind_lovo, best_val_lovo = LovoFun(θ)

    ResFun!(θ, ind_lovo, val_res, jac_res)

    BLAS.gemv!('T', 1.0, jac_res, val_res, 0.0, grad_lovo)

    ngrad_lovo = norm(grad_lovo, 2)

    safecount = 1

    # Main loop

    while (ngrad_lovo >= ε) && (safecount < MAXITER)

        with_logger(lm_logger.x) do

            @info("Iteration $(safecount)")
            @info("  Current value:   $(best_val_lovo)")
            @info("  ||grad_lovo||_2: $(ngrad_lovo)")
            @info("  Current iterate: $(θ)")
            @info("  Best indices (first $(maxoutind)): $(ind_lovo[1:maxoutind])")
            @info("  lambda: $(λ)")

        end

        G .= Id

        BLAS.gemm!('T', 'N', 1.0, jac_res, jac_res, λ, G)

        try

            # TODO: use a smart version of the QR factorization
            #F = qr!(G)
            F = cholesky!(G, Val(true))

            ldiv!(d, F, grad_lovo)

            d .*= -1.0

        catch

            with_logger(lm_logger.x) do

                @warn "Failed to solve the linear system. Will try new point."
                d .= - 1.0 .* grad_lovo
                θ .= rand(n)

            end

        end

        θnew .= θ .+ d

        ind_lovo, val_lovo = LovoFun(θnew)

        if  val_lovo <= best_val_lovo

            θ .= θnew

            best_val_lovo = val_lovo

            λ = λ / λ_down

            ResFun!(θ, ind_lovo, val_res, jac_res)

            BLAS.gemv!('T', 1.0, jac_res, val_res, 0.0, grad_lovo)

            ngrad_lovo = norm(grad_lovo, 2)

            with_logger(lm_logger.x) do

                @info("  Better function value found, lambda changed to $(λ).")

            end

        else

            λ = λ * λ_up

            with_logger(lm_logger.x) do

                @info("  No improvement, lambda changed to $(λ).")

            end

        end

        safecount += 1

    end

    if safecount == MAXITER
        with_logger(lm_logger.x) do

            @info("No solution was found in $(safecount) iterations.")

        end
        status = 0
    end

    # TODO: Create a test for this case
    if isnan(ngrad_lovo)
        with_logger(lm_logger.x) do

            @info("Incorrect value for gradient norm $(ngrad_lovo).")

        end
        status = 0
    end

    outliers = [1:npun;]
    setdiff!(outliers, ind_lovo)

    with_logger(lm_logger.x) do

        @info("""

        Final iteration (STATUS=$(status))
          Solution found:       $(θ)
          ||grad_lovo||_2:      $(ngrad_lovo)
          Function value:       $(best_val_lovo)
          Number of iterations: $(safecount)
          Outliers:             $(outliers)

        """)

    end

    return RAFFOutput(status, θ, safecount, p, best_val_lovo, nf, nj, outliers)

end

function lmlovo(model::Function, θ::Vector{Float64}, data::Array{Float64,2},
                n::Int, p::Int; kwargs...)

    # Define closures for derivative and initializations

    # 'x' is considered as global parameter for this function
    model_cl(θ) = model(x, θ)

    grad_model!(g, x_, θ) = begin

        global x = x_

        ForwardDiff.gradient!(g, model_cl, θ)

    end

    return lmlovo(model, grad_model!, θ, data, n, p; kwargs...)

end

lmlovo(model::Function, gmodel!::Function, data::Array{Float64,2}, n::Int, p::Int; kwargs...) =
    lmlovo(model, gmodel!, zeros(Float64, n), data, n, p; kwargs...)

lmlovo(model::Function, data::Array{Float64,2}, n::Int, p::Int; kwargs...) =
    lmlovo(model, zeros(Float64, n), data, n, p; kwargs...)

"""

    gnlslovo(model, gmodel!, θ, data::Array{T, 2}, n, p;
             ε::Number=1.0e-4, MAXITER=400, αls=2.0, dinc=2.0,
             MAXLSITER=100) where {T<:Float64}

    gnlslovo(model, θ::Vector{Float64}, data::Array{Float64,2},
             n::Int, p::Int; kwargs...)

    gnlslovo(model, gmodel!, data::Array{Float64,2}, n::Int,
             p::Int; kwargs...)

    gnlslovo(model, data::Array{Float64,2}, n::Int, p::Int; kwargs...)

LOVO Gauss-Newton with line-search described in

> R. Andreani, G. Cesar, R. M. Cesar-Jr., J. M. Martínez, and
> P. J. S. Silva, “Efficient curve detection using a {Gauss-Newton}
> method with applications in agriculture,” in Proc. 1st International
> Workshop on Computer Vision Applications for Developing Regions in
> Conjunction with ICCV 2007-CVDR-ICCV07, 2007.

Fit the `n`-parameter model `model` to the data given by matrix
`data`. The strategy is based on the LOVO function, which means that
only `p` (0 < `p` <= rows of `data`) points are trusted.

Matriz `data` is the data to be fit. This matrix should be in the form

    x11 x12 ... x1N y1
    x21 x22 ... x2N y2
    :

where `N` is the dimension of the argument of the model
(i.e. dimension of `x`).

If `θ` is provided, then it is used as the starting point.

The signature of function `model` should be given by

    model(x::Union{Vector{Float64}, SubArray}, θ::Vector{Float64})

where `x` are the variables and `θ` is a `n`-dimensional vector of
parameters. If the gradient of the model `gmodel!`

    gmodel! = (g::SubArray, x::Union{Vector{Float64}, SubArray},
               θ::Vector{Float64})

is not provided, then the function ForwardDiff.gradient! is called to
compute it.  **Note** that this choice has an impact in the
computational performance of the algorithm. In addition, if
`ForwardDiff.jl` is being used, then one **MUST** remove the signature
of vector `θ` from function `model`.

The optional arguments are

  - `MAXITER`: maximum number of iterations
  - `ε`: tolerance for the gradient of the function
  - `αls`: number >1 to increase/decrease the parameter `t` in
    line-search
  - `dinc`: number >1 to increase the diagonal of the J^T J
    matrix in order to escape from singularity
  - `MAXLSITER`: maximum number of Linear System increases in
    diagonal before exiting. Also defines the maximum number of
    Line Search trials to satisfy Armijo (but does not exit in such
    case)

Returns a [`RAFFOutput`](@ref) object.

"""
function gnlslovo(model, gmodel!, θ, data::Array{T, 2}, n, p;
                  ε::Number=1.0e-4, MAXITER=400, αls=2.0, dinc=2.0,
                  MAXLSITER=100) where {T<:Float64}

    @assert(n > 0, "Dimension should be positive.")
    @assert(p >= 0, "Trusted points should be nonnegative.")

    npun, = size(data)

    with_logger(lm_logger.x) do

        @debug("Size of data matrix ", size(data))

    end

    # Counters for calls to F and its Jacobian
    nf = 0

    nj = 0

    (p == 0) && return RAFFOutput(1, θ, 0, p, 0.0, nf, nj, [1:npun;])

    # Main function - the LOVO function
    LovoFun = let

        npun_::Int = npun

        ind::Vector{Int} = Vector{Int}(undef, npun_)

        F::Vector{Float64} = Vector{Float64}(undef, npun_)

        p_::Int = p

        # Return a ordered set index and lovo value
        (θ) -> begin

            nf += 1

            @views for i = 1:npun_
                F[i] = (model(data[i,1:(end - 1)], θ) - data[i, end])^2
            end

            indF, orderedF = sort_fun!(F, ind, p_)

            return indF, sum(orderedF)
        end

    end

    # Residue and Jacobian of residue
    rp::Vector{Float64} = Vector{Float64}(undef, p)

    Jrp::Array{Float64, 2} = Array{Float64}(undef, p, n)

    # This function returns the residue and Jacobian of residue
    ResFun! = let

        data_    = data
        model_   = model
        gmodel!_ = gmodel!

        (θ::Vector{Float64}, ind, r::Vector{Float64},
         rJ::Array{Float64, 2}) -> begin

             nj += 1

             residual_fg!(model_, gmodel!_, data_, θ, ind, r, rJ)

         end

    end

    # Gauss-Newton algorithm

    maxoutind = min(p, 5)
    dtnf      = Inf
    status    = 1

    # Allocation
    θnew = Vector{Float64}(undef, n)
    d    = Vector{Float64}(undef, n)
    y    = Vector{Float64}(undef, n)
    G    = Array{Float64, 2}(undef, n, n)
    ∇f  = Vector{Float64}(undef, n)

    # Initial information

    ip, fk = LovoFun(θ)
    ResFun!(θ, ip, rp, Jrp)

    BLAS.gemv!('T', 1.0, Jrp, rp, 0.0, ∇f)

    norm∇f = norm(∇f, 2)

    safecount = 1

    safelscnt = 0

    # Main loop

    while (norm∇f >= ε) && (safecount < MAXITER)

        with_logger(lm_logger.x) do

            @info("Iteration $(safecount)")
            @info("  Current value:   $(fk)")
            @info("  ||∇f||_2: $(norm∇f)")
            @info("  Current iterate: $(θ)")
            @info("  Best indices (first $(maxoutind)): $(ip[1:maxoutind])")

        end

        # Solve the system

        BLAS.gemm!('T', 'N', 1.0, Jrp, Jrp, 0.0, G)

        safelscnt = 0

        while true && (safelscnt < MAXLSITER)

            F = qr(G)

            serror = false
            dtnf   = 0.0
            try
                ldiv!(d, F, ∇f)
                dtnf = - dot(d, ∇f)
            catch
                serror = true
            end

            with_logger(lm_logger.x) do
                @debug("  Found solution for linear system", d, dtnf)
            end

            (!serror) && (abs(dtnf) >= 1.0e-8) && break

            with_logger(lm_logger.x) do
                @debug("  Unable to solve the GN system (dtnf = $(dtnf))." *
                       "  Will increase the diagonal elements.")
            end

            # Adjust the diagonal elements of G in case of failure
            for i in diagind(G)
                G[i] = max(dinc * G[i], 1.0)
            end

            safelscnt += 1

        end

        if safelscnt >= MAXLSITER
            with_logger(lm_logger.x) do
                @debug("Unable to compute a descent direction.")
            end

            break
        end

        d .*= -1.0

        # Line search
        θnew .= θ .+ d

        t = 1.0

        safelscnt = 0

        # Warning. We cannot store ip to save computation, due to the
        # way that LovoFun is implemented (ip is a pointer to a hidden
        # vector, actually). This will change in the future
        ip, fkpd = LovoFun(θnew)

        while (fkpd <= fk + t * dtnf) && (safelscnt < MAXLSITER)
            t         *= αls
            θnew      .= θ .+ t .* d
            ip, fkpf   = LovoFun(θnew)
            safelscnt += 1
        end

        if t > 1.0
            t       /= αls
            θnew    .= θ .+ t .* d
            ip, fkpf = LovoFun(θnew)
        end

        # TODO: implement quadratic interpolation
        while (fkpd > fk + t * dtnf) && (safelscnt < MAXLSITER)
            t         /= αls
            θnew      .= θ .+ t .* d
            ip, fkpd   = LovoFun(θnew)
            safelscnt += 1
        end

        with_logger(lm_logger.x) do

            @info("""

                \tLine search
                \t  t          = $(t)
                \t  f(θ + t d) = $(fkpd)
                \t     d^T ∇f = $(dtnf)
                """)

            if safelscnt == MAXLSITER
                @warn("Armijo condition was not satisfied within" *
                      " $(MAXLSITER) iterations.")
            end

        end

        # New iteration

        θ .= θnew
        fk = fkpd

        ResFun!(θ, ip, rp, Jrp)

        BLAS.gemv!('T', 1.0, Jrp, rp, 0.0, ∇f)

        norm∇f = norm(∇f, 2)

        safecount += 1

    end

    if (norm∇f >= ε) && (safelscnt >= MAXLSITER)
        with_logger(lm_logger.x) do

            @info("Unable to compute a descent solution in " *
                  "($MAXLSITER) attempts.")

        end
        status = 0
    end

    if (norm∇f >= ε) && (safecount == MAXITER)
        with_logger(lm_logger.x) do

            @info("No solution was found in $(safecount) iterations.")

        end
        status = 0
    end

    # TODO: Create a test for this case
    if isnan(norm∇f)
        with_logger(lm_logger.x) do

            @info("Incorrect value for gradient norm $(norm∇f).")

        end
        status = 0
    end

    safecount -= 1

    outliers = [1:npun;]
    setdiff!(outliers, ip)

    with_logger(lm_logger.x) do

        @info("""

        Final iteration (STATUS=$(status))
          Solution found:       $(θ)
          ||grad_lovo||_2:      $(norm∇f)
          Function value:       $(fk)
          Number of iterations: $(safecount)
          Outliers:             $(outliers)

        """)

    end

    return RAFFOutput(status, θ, safecount, p, fk, nf, nj, outliers)

end

function gnlslovo(model, θ::Vector{Float64}, data::Array{Float64,2},
                  n::Int, p::Int; kwargs...)

    # Define closures for derivative and initializations

    # 'x' is considered as global parameter for this function
    model_cl(θ) = model(x, θ)

    grad_model!(g, x_, θ) = begin

        global x = x_

        ForwardDiff.gradient!(g, model_cl, θ)

    end

    return gnlslovo(model, grad_model!, θ, data, n, p; kwargs...)

end

gnlslovo(model, gmodel!, data::Array{Float64,2}, n::Int, p::Int; kwargs...) =
    gnlslovo(model, gmodel!, zeros(Float64, n), data, n, p; kwargs...)

gnlslovo(model, data::Array{Float64,2}, n::Int, p::Int; kwargs...) =
    gnlslovo(model, zeros(Float64, n), data, n, p; kwargs...)


text:
# This file was generated, do not modify it. # hide
nrows = size(first(dfs), 1)
for (df, name) in zip(dfs, col_mean)
    df[!, name] = zeros(nrows)
    for (i, row) in enumerate(eachrow(df))
      df[i, name] = mean(row)
    end
end;


text:
using Distributed
using SharedArrays
addprocs(8)
nprocs()
@everywhere using FRGRealTime, DelimitedFiles,Dierckx,Plots
@everywhere kdata=readdlm("/home/tyy/Documents/CTP-fRG-Test/Real-Time-data/realtime_zero_k.dat")[:,1]|>reverse
@everywhere lamdata=readdlm("/home/tyy/Documents/CTP-fRG-Test/Real-Time-data/realtime_zero_lamdak.dat")[:,1]|>reverse
@everywhere m2data=readdlm("/home/tyy/Documents/CTP-fRG-Test/Real-Time-data/realtime_zero_m2k.dat")[:,1]|>reverse
@everywhere lamfun=Spline1D(kdata,lamdata)
@everywhere m2fun=Spline1D(kdata,m2data)

m2fun
p0out=SharedArray(collect(0.0001:0.0001:0.1))

p0max=1e-3
p0min=1e-5
psmax=30.0
psmin=1.0
p0out=SharedArray(exp.(collect(log(p0min):(log(p0max)-log(p0min))/20:log(p0max))))
psout=SharedArray(exp.(collect(log(psmin):(log(psmax)-log(psmin))/80:log(psmax))))

dyv1=[(i,j) for j in psout for i in p0out]
outvdy = pmap(
    pv -> FRGRealTime.fastpropImintqs_All(
        pv[1],
        pv[2],
        20.43287,
        1.0,
        800.0,
        4.0,
        m2fun,
        lamfun,
        rtol=1e-10,atol=1e-10,maxevals=20000000,initdiv=180,
    ),
    dyv1,
)

plot(p0out,outvdy)

writedlm(path,data)

plot(psout,-(outvdy)./psout.^2,xaxis=:log,yaxis=:log,legend=:topleft)
plot(p0out,outvdy)

dir="/home/wjfu1/tyy/FRG-RealTime-data/DynamicalCE"
writedlm("$dir/dyp0_T=21.43287.dat",p0out)

writedlm("$dir/dyps_T=20.43287.dat",psout)
writedlm("$dir/dyp0_T=20.43287.dat",p0out)

writedlm("$dir/dy_T=20.43287.dat",outvdy)

dyfun=Spline1D(p0out,outvdy)
plot(p0out,outvdy)

plot(x->derivative(dyfun,x),p0min,p0max)

p0out=SharedArray(collect(0.0001:0.0001:0.001))

psout=SharedArray(exp.(collect(log(0.002):(log(0.1)-log(0.002))/80:log(0.1))))

@everywhere myfun(p0, ps) =
    FRGRealTime.fastpropImintqs_All(
        p0,
        ps,
        145.0,
        1.0,
        400.0,
        4.0,
        m2fun,
        lamfun,
    ) + FRGRealTime.fastpropImintqs_All_dq0(
        p0,
        ps,
        145.0,
        1.0,
        400.0,
        4.0,
        m2fun,
        lamfun,rtol=1e-4,atol=1e-4
    )

propp0out4=pmap(ps->myfun.(p0out, ps),psout)

writedlm("/home/tyy/Documents/CTP-fRG-Test/Real-Time-data/Λ=400_p0_0.001-0.01_ps.dat",p0out)
writedlm("/home/tyy/Documents/CTP-fRG-Test/Real-Time-data/Λ=400_p0_ps_0.002-1.0.dat",psout)
writedlm("/home/tyy/Documents/CTP-fRG-Test/Real-Time-data/Λ=400_Gamma2p0_0.001-0.01_ps_0.002-1.0_all.dat",propp0out4)

outvIm = pmap(
    p0 -> FRGRealTime.fastpropImintqs_All(
        p0,
        0.000001,
        145.0,
        1.0,
        400.0,
        4.0,
        m2fun,
        lamfun,
    ),
    v1,
)

plot(v1,outvIm)

plot(v1,(p0->FRGRealTime.deltasumImprop1(
    p0,
    0.05,
    145.0,
    1.0,
    400.0,
    4.0,
    m2fun,
    lamfun,
)).(v1).+(p0->FRGRealTime.deltasumImprop2(
    p0,
    0.05,
    145.0,
    1.0,
    400.0,
    4.0,
    m2fun,
    lamfun,
)).(v1))

outv1|>maximum

plot(p0->FRGRealTime.flowpm(p0,10.0,20.0,-10.0,145.0),0.0,10.0)

plot(p0->FRGRealTime.flowpp_intcostheqs(p0,0.00,100.0,100.0,m2fun(100.0),145.0),1.0,200.0)