text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.1
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://discourse.julialang.org/t/ode-solvers-why-is-matlab-ode45-uncannily-stable/63052/15

# %%
using DifferentialEquations
using Plots

# %%
function wk5!(dP, P, params, t)

    #=

    The 5-Element WK with serial L
    set Ls = 0 for 4-Element WK parallel

    Formulation for DifferentialEquations.jl

    P:          solution vector (pressures p1 and p2)
    params: parameter tuple
                (Rc, Rp, C, Lp, Ls, I, q)

    I need to find a way to tranfer the function name as well
    for the time being we have to have the function in "I"

    =#

    # Split parameter tuple:
    Rc, Rp, C, Lp, Ls, I, q = params

    dP[1] = (
        -Rc / Lp * P[1]
        + (Rc / Lp - 1 / Rp / C) * P[2]
        + Rc * (1 + Ls / Lp) * didt(I, t, q)
        + I(t, q) / C
        )

    dP[2] = -1 / Rp / C * P[2] + I(t, q) / C

    return

end

# %%
# Generic Input Waveform
# max volume flow in ml/s
max_i = 425

# min volume flow in m^3/s
min_i = 0.0

T = 0.9

# Syst. Time in s
systTime = 2 / 5 * T

# Dicrotic notch time
dicrTime = 0.02

q_generic = (max_i, min_i, T, systTime, dicrTime)

function I_generic(t, q_generic)
    max_i, min_i, T, systTime, dicrTime = q_generic
    # implicit conditional using boolean multiplicator
    # sine waveform
    (
        (max_i - min_i) * sin(pi / systTime * (t % T))
        * (t % T < (systTime + dicrTime) )
        + min_i
    )
end

# %%
function didt(I, t, q)
    dt = 1e-3;
    didt = (I(t+dt, q) - I(t-dt, q)) / (2 * dt)
    return didt
end

# %%
plot(range(0, 2; length=2000), t -> didt(I_generic, t, q_generic); label="didt")

# %%
# Initial condition and time span
P0 = [0.0, 0.0]
tspan = (0.0, 30.0)

# Set parameters for Windkessel Model
Rc = 0.033
Rp = 0.6
C = 1.25
# L for serial model!
Ls = 0.01
# L for parallel
Lp = 0.02

I = I_generic
q = q_generic

p5 = (Rc, Rp, C, Lp, Ls, I, q)

problem = ODEProblem(wk5!, P0, tspan, p5)

dtmax = 1e-4

@time solutionTsit = solve(problem); GC.gc()
@time solutionTsitLowdt = solve(problem, Tsit5(), dtmax=dtmax); GC.gc()
@time solutionBS3 = solve(problem, BS3()); GC.gc()
@time solutionBS3Lowdt = solve(problem, BS3(), dtmax=dtmax); GC.gc()
@time solutionDP5 = solve(problem,DP5()); GC.gc()
@time solutionDP5Lowdt = solve(problem,DP5(), dtmax=dtmax); GC.gc()
@time solutionStiff = solve(problem, alg_hints=[:stiff]); GC.gc()

# %%
@time solutionTsit = solve(problem); GC.gc()
@time solutionTsitLowdt = solve(problem, Tsit5(), dtmax=dtmax); GC.gc()
@time solutionBS3 = solve(problem, BS3()); GC.gc()
@time solutionBS3Lowdt = solve(problem, BS3(), dtmax=dtmax); GC.gc()
@time solutionDP5 = solve(problem,DP5()); GC.gc()
@time solutionDP5Lowdt = solve(problem,DP5(), dtmax=dtmax); GC.gc()
@time solutionStiff = solve(problem, alg_hints=[:stiff]); GC.gc()

# %%
@show length(solutionTsit.t)
@show length(solutionTsitLowdt.t)
@show length(solutionBS3.t)
@show length(solutionBS3Lowdt.t)
@show length(solutionDP5.t)
@show length(solutionDP5Lowdt.t)
@show length(solutionStiff.t);

# %%
a, b = 0, 2

plot()
#plot!(t -> solutionTsit(t; idxs=1), a, b; label="Tsit")
plot!(t -> solutionTsitLowdt(t; idxs=1), a, b; label="TsitLowdt")
#plot!(t -> solutionBS3(t; idxs=1), a, b; label="BS3", ls=:dash)
plot!(t -> solutionBS3Lowdt(t; idxs=1), a, b; label="BS3Lowdt", ls=:dash)
#plot!(t -> solutionDP5(t; idxs=1), a, b; label="DP5", ls=:dashdot)
plot!(t -> solutionDP5Lowdt(t; idxs=1), a, b; label="DP5Lowdt", ls=:dashdot)
plot!(t -> solutionStiff(t; idxs=1), a, b; label="Stiff", ls=:dot, lw=1.5)

# %%
a, b = 0, 2

plot()
plot!(t -> solutionTsit(t; idxs=1), a, b; label="Tsit")
#plot!(t -> solutionTsitLowdt(t; idxs=1), a, b; label="TsitLowdt")
plot!(t -> solutionBS3(t; idxs=1), a, b; label="BS3", ls=:dash)
#plot!(t -> solutionBS3Lowdt(t; idxs=1), a, b; label="BS3Lowdt", ls=:dash)
plot!(t -> solutionDP5(t; idxs=1), a, b; label="DP5", ls=:dashdot)
#plot!(t -> solutionDP5Lowdt(t; idxs=1), a, b; label="DP5Lowdt", ls=:dashdot)
plot!(t -> solutionStiff(t; idxs=1), a, b; label="Stiff", ls=:dot, lw=1.5)

# %%
a, b = 28, 30

plot()
plot!(t -> solutionTsit(t; idxs=1), a, b; label="Tsit")
#plot!(t -> solutionTsitLowdt(t; idxs=1), a, b; label="TsitLowdt")
plot!(t -> solutionBS3(t; idxs=1), a, b; label="BS3", ls=:dash)
#plot!(t -> solutionBS3Lowdt(t; idxs=1), a, b; label="BS3Lowdt", ls=:dash)
plot!(t -> solutionDP5(t; idxs=1), a, b; label="DP5", ls=:dashdot)
#plot!(t -> solutionDP5Lowdt(t; idxs=1), a, b; label="DP5Lowdt", ls=:dashdot)
plot!(t -> solutionStiff(t; idxs=1), a, b; label="Stiff", ls=:dot, lw=1.5)

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using MAT

file = matopen("D:\\repo\\Test_Balthazar\\mp_data\\118nodes\\test_0.mat")
varnames = names(file)
#println(varnames)
data = read(file, "mpc")
println(data["gen"][:,1:2])
# println(data["branch"][:,3:4])
close(file)

ext: jl


text:
function transform_inf_expr(integrating_depvars, dict_depvar_input, dict_depvars, integrating_variables, transform)
    τs = Symbolics.variables(:τ, 1:length(integrating_variables))
    τs = Symbol.(τs)
    dict_transformation_vars = Dict()
    dict_depvar_input_ = Dict()
    integrating_var_transformation = []
    j = 1

    for depvar in integrating_depvars
        indvars = dict_depvar_input[depvar]
        ans = []

        for i in 1:length(indvars)
            if indvars[i] ∈ integrating_variables
                push!(ans, τs[j])
                push!(integrating_var_transformation, τs[j])
                dict_transformation_vars[indvars[i]] = transform(τs[j])
                j += 1
            else
                push!(ans, indvars[i])
            end
        end

        dict_depvar_input_[depvar] = ans
    end

    this_eq_pair = Dict(map(intvars -> dict_depvars[intvars] => dict_depvar_input_[intvars], integrating_depvars))
    this_eq_indvars = unique(vcat(values(this_eq_pair)...))

    return dict_transformation_vars, this_eq_indvars,integrating_var_transformation
end

function v_inf(t)
    return :($t ./ (1 .- $t.^2))
end

function v_semiinf(t , a , upto_inf)
    a = first(a)

    if upto_inf == true
        return :($a .+ $t ./ (1 .- $t))
    else
        return :($a .+ $t ./ (1 .+ $t))
    end
end

function get_inf_transformation_jacobian(integrating_variable, _inf, _semiup, _semilw)
    j = []
        for var in integrating_variable
            if _inf[1]
                append!(j, [:((1+$var^2)/(1-$var^2)^2)])
            elseif _semiup[1] || _semilw[1]
                append!(j, [:(1/(1-$var)^2)])
            end
        end

    return j
end

function transform_inf_integral(lb, ub, integrating_ex, integrating_depvars, dict_depvar_input, dict_depvars, integrating_variable, eltypeθ; dict_transformation_vars = nothing, transformation_vars = nothing)
    lb_ = Symbolics.tosymbol.(lb)
    ub_ = Symbolics.tosymbol.(ub)

    if -Inf in lb_ || Inf in ub_

        if !(integrating_variable isa Array)
            integrating_variable = [integrating_variable]
        end

        lbb = lb_ .=== -Inf
        ubb = ub_ .=== Inf
        _none = .!lbb .& .!ubb
        _inf = lbb .& ubb
        _semiup = .!lbb .& ubb
        _semilw = lbb  .& .!ubb

        function transform_indvars(t)
            if _none[1]
                return t
            elseif _inf[1]
                return v_inf(t)
            elseif _semiup[1]
                return v_semiinf(t , lb , 1)
            elseif _semilw[1]
                return v_semiinf(t , ub , 0)
            end
        end

        dict_transformation_vars, transformation_vars, integrating_var_transformation = transform_inf_expr(integrating_depvars, dict_depvar_input, dict_depvars, integrating_variable,transform_indvars)

        ϵ = 1/20 #cbrt(eps(eltypeθ))

        lb = 0.00.*_semiup + (-1.00+ϵ).*_inf + (-1.00+ϵ).*_semilw +  _none.*lb
        ub = (1.00-ϵ).*_semiup + (1.00-ϵ).*_inf  + 0.00.*_semilw  + _none.*ub

        j = get_inf_transformation_jacobian(integrating_var_transformation, _inf, _semiup, _semilw)

        integrating_ex = Expr(:call, :*, integrating_ex, j...)
    end

    return lb, ub, integrating_ex, dict_transformation_vars, transformation_vars
end


text:
module Qlab
	using Compat, PyPlot

	export load_data, load_latest_data, KT_estimation, pauliAlphabet, filter_records, comp_to_paulis, savitzky_golay, read_records, digitize, unwrap, unwrap!, plot2D, plot1D, hinton

	include("load_data.jl")
	include("KT_estimation.jl")
	include("comp_to_paulis.jl")
	include("SavitzkyGolay.jl")
	include("read_records.jl")
	include("digitize.jl")
	include("common_fits.jl")
	include("tomography.jl")
	include("plotting.jl")
	include("shortcuts.jl")
	include("resonator_fits.jl")
	include("utils.jl")
end


text:
# Miscellaneous functions

# Logarithm of multivariate gamma function
#
# See: https://en.wikipedia.org/wiki/Multivariate_gamma_function
#
function logmvgamma(p::Int, a::Real)
    # NOTE: one(a) factors are here to prevent unnecessary promotion of Float32
    res = p * (p - 1) * log(pi * one(a)) / 4
    for ii in 1:p
        res += lgamma(a + (1 - ii) * one(a)/ 2)
    end
    return res
end

# Remainder term after Stirling's approximation to the log-gamma function
# lstirling(x) = lgamma(x) + x - (x-0.5)*log(x) - 0.5*log2π
#              = 1/(12x) - 1/(360x^3) + 1/(1260x^5) + ...
#
# Asymptotic expansion from:
#
#   Temme, N. (1996) Special functions: An introduction to the classical
#   functions of mathematical physics, Wiley, New York, ISBN: 0-471-11313-1,
#   Chapter 3.6, pp 61-65.
#
# Relative error of approximation is bounded by
#   (174611/125400 x^-19) / (1/12 x^-1 - 1/360 x^-3)
# which is < 1/2 ulp for x >= 10.0
# total numeric error appears to be < 2 ulps
#
function lstirling_asym(x::Float64)
    t = 1.0/(x*x)
    @horner(t,
             8.33333333333333333e-2, #  1/12 x^-1
            -2.77777777777777778e-3, # -1/360 x^-3
             7.93650793650793651e-4, #  1/1260 x^-5
            -5.95238095238095238e-4, # -1/1680 x^-7
             8.41750841750841751e-4, #  1/1188 x^-9
            -1.91752691752691753e-3, # -691/360360 x^-11
             6.41025641025641026e-3, #  1/156 x^-13
            -2.95506535947712418e-2, # -3617/122400 x^-15
             1.79644372368830573e-1)/x #  43867/244188 x^-17
end

function lstirling_asym(x::Float32)
    t = 1f0/(x*x)
    @horner(t,
             8.333333333333f-2, #  1/12 x^-1
            -2.777777777777f-3, # -1/360 x^-3
             7.936507936508f-4, #  1/1260 x^-5
            -5.952380952381f-4, # -1/1680 x^-7
             8.417508417508f-4)/x #  1/1188 x^-9
end

lstirling_asym(x::Integer) = lstirling_asym(float(x))


text:
function render!(
    rendermodel::RenderModel,
    veh::Vehicle,
    color::Colorant=RGB(rand(), rand(), rand())
    )

    p = veh.state.posG
    add_instruction!(rendermodel, render_vehicle, (p.x, p.y, p.θ, veh.def.length, veh.def.width, color))
    return rendermodel
end
function render!(
    rendermodel::RenderModel,
    veh::Entity{VehicleState, BicycleModel, Int},
    color::Colorant=RGB(rand(), rand(), rand())
    )

    veh2 = Vehicle(veh.state, veh.def.def, veh.id)
    render!(rendermodel, veh2, color)
end


text:
using DataDeps

register(DataDep(
  "monks-problems-2",
  "http://archive.ics.uci.edu/ml/datasets/MONK%27s+Problems",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/monks-problems/monks-2.train",
  "7ebcc278005737caad0cd52e6aa31256a2c3174da8047faf4ac3a9d35edbdaa8",
  post_fetch_method=(path -> begin
    # First char of each row is whitespace, remove
    write(path, join(strip.(readlines(path)), "\n"))

    UCIData.process_dataset(path,
        target_index=1,
        feature_indices=2:7,
        id_indices=8,
        categoric_indices=2:7,
        delim=' ',
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
export
    initialize_XML,
    get_lane_length,
    handle_junctions,
    convert_roadway!,
    Connection,
    Junction,
    connect_two_lane_general!,
    connect_two_seg!,
    connect_two_seg_general!,
    add_connection!,
    add_junction!,
    gen_connected_lanes,
    gen_intersection,
    gen_loop_roadway,
    add_line!,
    in_lanes,
    get_max_curvature

include("roadway_generation.jl")
include("roadway_perception.jl")
include("roadway_intersection.jl")
include("convert2xodr.jl")


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module Constraint

import LinearAlgebra
import MathOptInterface
import OrderedCollections: OrderedDict
import SparseArrays

const MOI = MathOptInterface
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

include("bridge.jl")
include("function_conversion.jl")
include("map.jl")
include("set_map.jl")
include("single_bridge_optimizer.jl")

include("bridges/det.jl")
include("bridges/flip_sign.jl")
include("bridges/functionize.jl")
include("bridges/geomean_to_relentr.jl")
include("bridges/geomean.jl")
include("bridges/indicator_activate_on_zero.jl")
include("bridges/indicator_sos.jl")
include("bridges/interval.jl")
include("bridges/ltgt_to_interval.jl")
include("bridges/norm_spec_nuc_to_psd.jl")
include("bridges/norm_to_lp.jl")
include("bridges/quad_to_soc.jl")
include("bridges/relentr_to_exp.jl")
include("bridges/rsoc_soc.jl")
include("bridges/scalarize.jl")
include("bridges/semi_to_binary.jl")
include("bridges/slack.jl")
include("bridges/soc_rsoc.jl")
include("bridges/soc_to_nonconvex_quad.jl") # do not add these bridges by default
include("bridges/soc_to_psd.jl")
include("bridges/square.jl")
include("bridges/vectorize.jl")
include("bridges/zero_one.jl")

"""
    add_all_bridges(bridged_model, ::Type{T}) where {T}

Add all bridges defined in the `Bridges.Constraint` submodule to
`bridged_model`. The coefficient type used is `T`.
"""
function add_all_bridges(bridged_model, ::Type{T}) where {T}
    if T <: AbstractFloat
        MOI.Bridges.add_bridge(bridged_model, GreaterToIntervalBridge{T})
        MOI.Bridges.add_bridge(bridged_model, LessToIntervalBridge{T})
    end
    MOI.Bridges.add_bridge(bridged_model, GreaterToLessBridge{T})
    MOI.Bridges.add_bridge(bridged_model, LessToGreaterBridge{T})
    MOI.Bridges.add_bridge(bridged_model, NonnegToNonposBridge{T})
    MOI.Bridges.add_bridge(bridged_model, NonposToNonnegBridge{T})
    MOI.Bridges.add_bridge(bridged_model, ScalarizeBridge{T})
    MOI.Bridges.add_bridge(bridged_model, VectorizeBridge{T})
    MOI.Bridges.add_bridge(bridged_model, ScalarSlackBridge{T})
    MOI.Bridges.add_bridge(bridged_model, VectorSlackBridge{T})
    MOI.Bridges.add_bridge(bridged_model, ScalarFunctionizeBridge{T})
    MOI.Bridges.add_bridge(bridged_model, VectorFunctionizeBridge{T})
    MOI.Bridges.add_bridge(bridged_model, SplitIntervalBridge{T})
    MOI.Bridges.add_bridge(bridged_model, QuadtoSOCBridge{T})
    # We do not add `(R)SOCtoNonConvexQuad` because it starts with a convex
    # conic constraint and generate a non-convex constraint (in the QCP
    # interpretation).
    MOI.Bridges.add_bridge(bridged_model, NormInfinityBridge{T})
    MOI.Bridges.add_bridge(bridged_model, NormOneBridge{T})
    MOI.Bridges.add_bridge(bridged_model, GeoMeantoRelEntrBridge{T})
    MOI.Bridges.add_bridge(bridged_model, GeoMeanBridge{T})
    MOI.Bridges.add_bridge(bridged_model, RelativeEntropyBridge{T})
    MOI.Bridges.add_bridge(bridged_model, NormSpectralBridge{T})
    MOI.Bridges.add_bridge(bridged_model, NormNuclearBridge{T})
    MOI.Bridges.add_bridge(bridged_model, SquareBridge{T})
    MOI.Bridges.add_bridge(bridged_model, LogDetBridge{T})
    MOI.Bridges.add_bridge(bridged_model, RootDetBridge{T})
    MOI.Bridges.add_bridge(bridged_model, RSOCtoSOCBridge{T})
    MOI.Bridges.add_bridge(bridged_model, SOCtoRSOCBridge{T})
    # We do not add `SOCtoPSDBridge` as transforming the `SOC` to `RSOC` and
    # then to `PSD` produces a smaller SDP constraint.
    MOI.Bridges.add_bridge(bridged_model, RSOCtoPSDBridge{T})
    MOI.Bridges.add_bridge(bridged_model, IndicatorActiveOnFalseBridge{T})
    MOI.Bridges.add_bridge(bridged_model, IndicatorSOS1Bridge{T})
    MOI.Bridges.add_bridge(bridged_model, SemiToBinaryBridge{T})
    MOI.Bridges.add_bridge(bridged_model, ZeroOneBridge{T})
    return
end

end


text:
"""
function checkOptions{T}(n::Int64, x::Vector{T}, xScal::Vector{T},
    opt::OptionsNLEQ)

Checking of common input parameters and options.

T = Float64 or BigFloat

## Input parameters
-------------------
| Variable | Description             |
|:---------|:------------------------|
| n        | Size of the problem     |
| x        | Initial guess           |
| xScal*   | Initial scaling vector  |
| opt*     | Options set by the user |

(* marks inout parameters)

## Output parameters
--------------------
| Variable | Description                 |
|:---------|:----------------------------|
| retCode  | Exit code in case of errors |
"""
function checkOptions{T}(n::Int64, x::Vector{T}, xScal::Vector{T},
    opt::OptionsNLEQ)
    # Begin
    (epMach, small, great) = getMachineConstants(T)
    # Check the IO
    pIOwarn = getOption!(opt, OPT_PRINTIOWARN, STDOUT)
    if typeof(pIOwarn) != IOStream && pIOwarn != STDOUT
        retCode = 30
        write(STDOUT, "ERROR: Please provide a file stream for writing warnings or\n",
        "use the default option.\n")
        return retCode
    end

    pIOmon = getOption!(opt, OPT_PRINTIOMON, STDOUT)
    if typeof(pIOmon) != IOStream && pIOmon != STDOUT
        retCode = 30
        write(STDOUT, "ERROR: Please provide a file stream for writing iterations or\n",
        "use the default option.\n")
        return retCode
    end

    pIOsol = getOption!(opt, OPT_PRINTIOSOL, STDOUT)
    if typeof(pIOsol) != IOStream && pIOsol != STDOUT
        retCode = 30
        write(STDOUT, "ERROR: Please provide a file stream for writing solution or\n",
        "use the default option.\n")
        return retCode
    end
    # Check whether warnings need to be printed
    printWarn = getOption!(opt,OPT_PRINTWARNING,0)

    # Initialize return code to 0
    retCode = 0

    # Check dimensional parameter n
    if n <= 0
        retCode = 20
        write(pIOwarn,"ERROR: Bad input to dimensional parameter n supplied","\n",
            "Choose n positive, your input is: n = $n\n")
        return retCode
    end

    # Problem type specification by user
    nonLin = getOption!(opt,OPT_NONLIN,3)

    # Checking and conditional adaptation of user given RTOL
    # if RTOL is not set, set it to 1e-6
    if T == Float64
        rTol = getOption!(opt,OPT_RTOL,1e-6)
    elseif T == BigFloat
        rTol = getOption!(opt,OPT_RTOL,BigFloat(1e-6))
    end
    if rTol <= 0.0
        retCode = 21
        write(pIOwarn,"ERROR: Nonpositive $OPT_RTOL supplied\n")
        return retCode
    else
        tolMin = epMach*10.0*n
        if rTol < tolMin
            rTol = tolMin
            setOption!(opt,OPT_RTOL,rTol)
            if printWarn == 1
                write(pIOwarn,"WARNING: User prescribed $OPT_RTOL increased to a reasonable smallest value RTOL = $rTol\n")
            end
        end

        tolMax = 1.0e-1
        if rTol > tolMax
            rTol = tolMax
            setOption!(opt,OPT_RTOL,rTol)
            if printWarn == 1
                write(pIOwarn,"WARNING: User prescribed $OPT_RTOL decreased to a reasonable largest value RTOL = $rTol\n")
            end
        end
    end

    # Test user prescribed accuracy and scaling on proper values
    if nonLin >= 3
        defScal = rTol
    else
        defScal = 1.0
    end

    for i = 1:n
        # Scaling Values cannot be negative
        # Positive scaling values give scale invariance
        if xScal[i] < 0.0
            retCode = 22
            write(pIOwarn,"ERROR: Negative value in xScal[$i] supplied\n")
            return retCode
        end

        if xScal[i] == 0.0
            xScal[i] = defScal
        end
        # Avoid overflow due to division by xScal[i]
        if xScal[i] > 0.0 && xScal[i] < small
            if printWarn == 1
                write(pIOwarn,"WARNING: xScal[$i] = $xScal[i] too small, increased to $small\n")
            end
            xScal[i] = small
        end
        # Avoid underflow due to division by xScal[i]
        if xScal[i] > great
            if printWarn == 1
                write(pIOwarn,"WARNING: xScal[$i] = $xScal[i] too big, increased to $great\n")
            end
            xScal[i] = great
        end
    end

    # Assign the Jacobian depending on user input
    # By default Forward mode automatic differentiation is used
    jacGen = getOption!(opt,OPT_JACGEN,4)
    if jacGen == 1
        jacFcn = getOption!(opt,OPT_JACFCN,0)
        if jacFcn == 0
            retCode = 30
            write(pIOwarn,"ERROR: The Jacobian function OPT_JACFCN is not supplied. ",
            "Please supply a Jacobian function or use OPT_JACGEN = 2 or 3 for numerical differentiation based jacobian evaluation.\n")
            return retCode
        end
    end

    return retCode
end

"""
function initializeOptions·(opt::OptionsNLEQ, wk::OptionsNLEQ, n::Int64,
    m1::Int64, nBroy::Int64, qRank1::Bool, solver::Int64, T::DataType)

Initialization of options based on the solver input argument.

## Input parameters
-------------------
| Variables  | Description                                                    |
|:-----------|:---------------------------------------------------------------|
| opt*       | Options set by the user                                        |
| wk*        | Internal workspace specific to the solver                      |
| n          | Size of the problem                                            |
| m1         | In full mode = n and in band mode = 2*ml+mu+1                  |
| nBroy      | Maximum number of possible consecutive iterative Broyden steps |
| qRank1     | Decision parameter for Rank-1 updates                          |
| solver = 1 | Specifies that the solver is NLEQ1                             |
|        = 2 | Specifies that the solver is NLEQ2                             |
| T          | Data type of the input. Acceptable values Float64 or BigFloat  |
(* marks inout parameters)
"""
function initializeOptions(opt::OptionsNLEQ, wk::OptionsNLEQ, n::Int64,
    m1::Int64, nBroy::Int64, qRank1::Bool, solver::Int64, T::DataType)
    # Begin
    # Initialize options: OPT
    initOption!(opt, OPT_FCMIN,     zero(T))
    initOption!(opt, OPT_SIGMA,     zero(T))
    initOption!(opt, OPT_SIGMA2,    zero(T))
    initOption!(opt, OPT_NOROWSCAL, 0)

    # Initialize workspace: WK
    if solver == 1
        initOption!(wk, WK_A, zeros(T,m1,n))
        if qRank1
            initOption!(wk, WK_DXSAVE, zeros(T,n,nBroy))
        else
            initOption!(wk, WK_DXSAVE, zero(T))
        end
    elseif solver == 2
        initOption!(wk, WK_QU, zeros(T,n))
        initOption!(wk, WK_A, zeros(T,n,n))
        if nBroy != 0
            initOption!(wk, WK_QA_DXSAVE, zeros(T,n,n))
        end
    end

    initOption!(wk, WK_DX  , zeros(T,n))
    initOption!(wk, WK_DXQ , zeros(T,n))
    initOption!(wk, WK_XA  , zeros(T,n))
    initOption!(wk, WK_XWA , zeros(T,n))
    initOption!(wk, WK_F   , zeros(T,n))
    initOption!(wk, WK_FA  , zeros(T,n))
    initOption!(wk, WK_ETA , zeros(T,n))
    initOption!(wk, WK_XW  , zeros(T,n))
    initOption!(wk, WK_FW  , zeros(T,n))
    initOption!(wk, WK_DXQA, zeros(T,n))

    initOption!(wk, WK_SUMXA0, zero(T))
    initOption!(wk, WK_SUMXA1, zero(T))
    initOption!(wk, WK_FCMON,  zero(T))
    initOption!(wk, WK_FCA,    zero(T))
    initOption!(wk, WK_FCKEEP, zero(T))
    initOption!(wk, WK_FCPRI,  zero(T))
    initOption!(wk, WK_DMYCOR, zero(T))
    initOption!(wk, WK_SUMXS,  zero(T))

    initOption!(wk, STATS_NITER,  0)
    initOption!(wk, STATS_NCORR,  0)
    initOption!(wk, STATS_NFCN,   0)
    initOption!(wk, STATS_NFCNJ,  0)
    initOption!(wk, STATS_NJAC,   0)
    initOption!(wk, STATS_NREJR1, 0)
    initOption!(wk, STATS_NEW,    0)
    initOption!(wk, STATS_ICONV,  0)
    initOption!(wk, STATS_CONV,   zero(T))
    initOption!(wk, STATS_SUMX,   zero(T))
    initOption!(wk, STATS_DLEVF,  zero(T))
    initOption!(wk, STATS_RTOL,   zero(T))
    initOption!(wk, "P_TOLALL", Vector{T}())

    # Check whether these iteration variables should be stored or not
    initOption!(opt, OPT_STORE, 0)
    if opt.options[OPT_STORE] == 1
        initOption!(wk, "P_XITER", Vector{Vector{T}}())
        initOption!(wk, "P_SUMXALL", Vector{T}())
        initOption!(wk, "P_DLEVFALL", Vector{T}())
        initOption!(wk, "P_SUMXQALL", Vector{T}())
        initOption!(wk, "P_FCALL", Vector{T}())
    end

    return nothing
end

"""
function initializeOptions(opt::OptionsNLEQ, solver::Int64, T::DataType)

Initialization of options based on the solver input argument.

## Input parameters
-------------------
| Variables  | Description                                                    |
|:-----------|:---------------------------------------------------------------|
| opt*       | Options set by the user                                        |
| solver = 1 | Specifies that the solver is NLEQ1                             |
|        = 2 | Specifies that the solver is NLEQ2                             |
| T          | Data type of the input. Acceptable values Float64 or BigFloat  |

(* marks inout parameters)
"""
function initializeOptions(n::Int64, opt::OptionsNLEQ, solver::Int64, T::DataType)
    # Begin
    (epMach, _, _) = getMachineConstants(T)
    # First call or successive call
    qSucc   = Bool(getOption!(opt,OPT_QSUCC,0))
    qIniMon = (opt.options[OPT_PRINTITERATION] >= 1 && !qSucc)

    # Check if the Jacobian is Dense/Sparse or Banded matrix
    mStor = getOption!(opt,OPT_MSTOR,0)
    ml = getOption!(opt,"OPT_ML",0)
    mu = getOption!(opt,"OPT_MU",0)
    if mStor == 0
        m1 = n
        m2 = n
    elseif mStor == 1
        m1 = 2*ml + mu + 1
        m2 = ml + mu + 1
    end

    qRank1 = Bool(getOption!(opt, OPT_QRANK1, 0))
    qOrdi  = Bool(getOption!(opt, OPT_QORDI,  0))
    qSimpl = Bool(getOption!(opt, OPT_QSIMPL, 0))

    if qRank1
        nBroy = getOption!(opt,OPT_NBROY,0)
        if nBroy == 0
            nBroy = max(m2,10)
            setOption!(opt,OPT_NBROY, nBroy)
        end
    else
        nBroy = 0
    end

    # If in stepwise mode and this is the first call then clear the workspace
    initOption!(opt, OPT_MODE, 0)

    if opt.options[OPT_MODE] == 1 && !qSucc
        if solver == 1
            empty!(wkNLEQ1.options)
        elseif solver == 2
            empty!(wkNLEQ1.options)
        end
    end

    # Check if this is a first call or successive call to nleq1
    # If first call then reset the workspace and persistent variables
    if !qSucc
        if solver == 1
            empty!(wkNLEQ1.options)
            initializeOptions(opt, wkNLEQ1, n, m1, nBroy, qRank1, 1, T)
        elseif solver == 2
            empty!(wkNLEQ2.options)
            initializeOptions(opt, wkNLEQ2, n, m1, nBroy, qRank1, 2, T)
        end
    end

    # Check for non linear option
    nonLin = opt.options[OPT_NONLIN]
    initOption!(opt, OPT_BOUNDEDDAMP, 0)

    if opt.options[OPT_BOUNDEDDAMP] == 0
        qBDamp = nonLin == 4
    elseif opt.options[OPT_BOUNDEDDAMP] == 1
        qBDamp = true
    elseif opt.options[OPT_BOUNDEDDAMP] == 2
        qBDamp = false
    end

    # Initialize bounded damping strategy restriction factor
    initOption!(opt, OPT_FCBAND, zero(T))
    if qBDamp
        if opt.options[OPT_FCBAND] < 1.0
            setOption!(opt, OPT_FCBAND, 10.0*one(T))
        end
    end

    # Maximum permitted number of iteration steps
    nItmax = getOption!(opt, OPT_NITMAX, 50)
    if nItmax <= 0
        nItmax = 50
        setOption!(opt, OPT_NITMAX, nItmax)
    end

    if qIniMon
        if solver == 1
            printInitialization(n, opt.options[OPT_PRINTIOMON], opt.options[OPT_RTOL],
            opt.options[OPT_JACGEN], mStor, ml, mu, opt.options[OPT_NOROWSCAL],
            qRank1, nonLin, qBDamp, opt.options[OPT_FCBAND], qOrdi, qSimpl, nItmax)
        elseif solver == 2
            printInitialization(n, opt.options[OPT_PRINTIOMON], opt.options[OPT_RTOL],
            opt.options[OPT_JACGEN], 0, 0, 0, opt.options[OPT_NOROWSCAL], qRank1,
            nonLin, qBDamp, opt.options[OPT_FCBAND], false, false, nItmax)
        end
    end

    # Initial damping factor for highly nonlinear problems
    initOption!(opt, OPT_FCSTART, zero(T))
    qFcStart = opt.options[OPT_FCSTART] > 0.0
    if !qFcStart
        setOption!(opt, OPT_FCSTART, T(1.0e-2))
        if nonLin == 4
            setOption!(opt, OPT_FCSTART, T(1.0e-4))
        end
    end

    # Minimal permitted damping factor
    initOption!(opt,OPT_FCMIN,zero(T))
    if opt.options[OPT_FCMIN] <= 0.0
        setOption!(opt, OPT_FCMIN, T(1.0e-4))
        if nonLin == 4
            setOption!(opt, OPT_FCMIN, T(1.0e-8))
        end
    end
    fcMin = opt.options[OPT_FCMIN]

    # Rank1 decision parameter SIGMA
    initOption!(opt,OPT_SIGMA,zero(T))
    if opt.options[OPT_SIGMA] < 1.0
        setOption!(opt, OPT_SIGMA, 3.0*one(T))
    end
    if !qRank1
        setOption!(opt, OPT_SIGMA, 10.0/fcMin)
    end

    # Decision parameter about increasing too small predictor
    # to greater corrector value
    initOption!(opt,OPT_SIGMA2,zero(T))
    if opt.options[OPT_SIGMA2] < 1.0
        setOption!(opt, OPT_SIGMA2, 10.0/fcMin)
    end

    # Starting value of damping factor (fcMin <= fc <= 1.0)
    if nonLin <= 2 && !qFcStart
        # for linear or mildly nonlinear problems
        fc = one(T)
    else
        # for highly or extremely nonlinear problems
        fc = getOption(opt, OPT_FCSTART, zero(T))
    end

    # Simplified Newton iteration implies ordinary Newton iteration mode
    if qSimpl
        setOption!(opt, OPT_QORDI, 1)
    end

    # If ordinary Newton iteration, damping factor is always 1
    if opt.options[OPT_QORDI] == 1
        fc = one(T)
    end

    # Set starting damping factor
    setOption!(opt, OPT_FCSTART, fc)

    if solver == 2
        iRank = getOption!(opt, OPT_IRANK, 0)
        if iRank <= 0 || iRank > n
            iRank = n
            setOption!(opt, OPT_IRANK, iRank)
        end

        cond = getOption!(opt, OPT_COND, 1.0/epMach)
        if cond < 1.0
            cond = 1.0/epMach
            setOption!(opt, OPT_COND, cond)
        end
    end

    if opt.options[OPT_PRINTITERATION] >= 2 && !qSucc
        write(opt.options[OPT_PRINTIOMON],"\nINFO: ","Internal parameters:",
        "\n\tStarting value for damping factor ",
        @sprintf("OPT_FCSTART\t= %1.2e",opt.options[OPT_FCSTART]),
        @sprintf("\n\tMinimum allowed damping factor OPT_FCMIN\t= %1.2e",fcMin),
        "\n\tRank-1 updates decision parameter ",
        @sprintf("OPT_SIGMA\t= %1.2e\n",opt.options[OPT_SIGMA]))
        if solver == 2
            write(opt.options[OPT_PRINTIOMON], @sprintf("\n\tInitial Jacobian pseudo-rank iRank\t\t= %6i", iRank),
            @sprintf("\n\tMaximum permitted subcondition cond\t\t= %1.2e\n", cond))
        end
    end

    return (m1, m2, nBroy, qBDamp)
end

"""
function printInitialization{T}(n::Int64, printIOmon, rTol::T, jacGen::Int64,
    mStor::Int64, ml::Int64, mu::Int64, qNoRowScal::Int64, qRank1::Bool, nonLin::Int64,
    qBDamp::Bool, fcBand::T, qOrdi::Bool, qSimpl::Bool, nItmax::Int64)

Print a summary of the initialization.

T = Float64 or BigFloat

## Input parameters
-------------------
| Variables  | Description                                        |
|:-----------|:---------------------------------------------------|
| n          | Size of the problem                                |
| printIOmon | IO handle for printing                             |
| rTol       | Relative tolerance                                 |
| jacGen     | Method of Jacobian generation                      |
| mStor      | Dense or band mode storage of Jacobian             |
| ml         | Lower bandwidth in case of band storage            |
| mu         | Upper bandwidth in case of band storage            |
| qNoRowScal | Decision parameter for automatic row scaling       |
| qRank1     | Decision parameter for Rank-1 updates              |
| nonLin     | Problem type specification                         |
| qBDamp     | Decision parameter for bounded damping strategy    |
| fcBand     | Bounded damping strategy restriction factor        |
| qOrdi      | Decision parameter for ordinary Newton iteration   |
| qSimpl     | Decision parameter for simplified Newton iteration |
| nItmax     | Maximum permitted Newton iterations                |
"""
function printInitialization{T}(n::Int64, printIOmon, rTol::T, jacGen::Int64,
    mStor::Int64, ml::Int64, mu::Int64, qNoRowScal::Int64, qRank1::Bool, nonLin::Int64,
    qBDamp::Bool, fcBand::T, qOrdi::Bool, qSimpl::Bool, nItmax::Int64)
    # Begin
    message = ""
    write(printIOmon,"\nINFO: ","N = $n\n")
    write(printIOmon,"\nINFO: ","Prescribed relative precision ",
    "$rTol\n")
    if jacGen == 1
        message = "a user function"
    elseif jacGen == 2
        message = "numerical differentation (without feedback strategy)"
    elseif jacGen == 3
        message = "numerical differentation (feedback strategy included)"
    end
    write(printIOmon,"\nINFO: ","The Jacobian is supplied by $message\n")
    if mStor == 0
        message = "full"
    elseif mStor == 1
        message = "banded"
    end
    write(printIOmon,"INFO: ","The Jacobian will be stored in $message mode\n")
    if mStor == 1
        write(printIOmon,"INFO: ","Lower bandwidth : $ml \t",
        "Upper bandwidth : $mu\n")
    end
    if qNoRowScal == 1
        message = "inhibited"
    else
        message = "allowed"
    end
    write(printIOmon,"INFO: ",
    "Automatic row scaling of the jacobian is $message\n")

    if qRank1
        message = "allowed"
    else
        message = "inhibited"
    end
    write(printIOmon,"\nINFO: ","Rank-1 updates are $message\n")
    if nonLin == 1
        message = "linear"
    elseif nonLin == 2
        message = "mildly nonlinear"
    elseif nonLin == 3
        message = "highly nonlinear"
    elseif nonLin == 4
        message = "extremely nonlinear"
    end
    write(printIOmon,"INFO: ","Problem is specified as being $message\n")
    if qBDamp
        write(printIOmon,"INFO: ","Bounded damping strategy is active\n",
        "bounding factor is $fcBand\n")
    else
        write(printIOmon,"INFO: ","Bounded damping strategy is off\n")
    end
    if qOrdi
        write(printIOmon,"INFO: ","Special mode: ",
        "Ordinary Newton iteration will be done\n")
    end
    if qSimpl
        write(printIOmon,"INFO: ","Special mode: ",
        "Simplified Newton iteration will be done\n")
    end

    write(printIOmon,"INFO: ","Maximum permitted number of ",
    "iteration steps : $nItmax\n")
end

"""
function printStats(stats::Dict{AbstractString,Any}, printIOmon)

Print a summary of the statistics.

## Input parameters
-------------------
| Variables  | Description                                      |
|:-----------|:-------------------------------------------------|
| stats      | Dictionary variable containing solver statistics |
| printIOmon | IO handle for printing                           |
"""
function printStats(stats::Dict{AbstractString,Any}, printIOmon)
    write(printIOmon,"\n",
    @sprintf("*************   Statistics   ************\n"),
    @sprintf("***  Newton-iterations     : %7i  ***\n", (stats[STATS_NITER])),
    @sprintf("***  Corrector steps       : %7i  ***\n", (stats[STATS_NCORR])),
    @sprintf("***  Rejected Rank-1 steps : %7i  ***\n", (stats[STATS_NREJR1])),
    @sprintf("***  Jacobian evaluations  : %7i  ***\n", (stats[STATS_NJAC])),
    @sprintf("***  Function evaluations  : %7i  ***\n", (stats[STATS_NFCN])),
    @sprintf("***  ... for Jacobain eval : %7i  ***\n", (stats[STATS_NFCNJ])),
    @sprintf("*****************************************\n"))
end

"""
function nScal{T}(n::Int64, x::Vector{T}, xa::Vector{T}, xScal::Vector{T},
    iScal::Int64, mPr::Int64, printIO, xw::Vector{T})

To be used in connection with NLEQ1 and NLEQ2.
Computation of the internal scaling vector XW used for the
Jacobian matrix, the iterate vector and it's related
vectors - especially for the solution of the linear system
and the computations of norms to avoid numerical overflow.

T = Float64 or BigFloat

## Input parameters
-------------------
| Variables | Description                                                                |
|:----------|:---------------------------------------------------------------------------|
| n         | Size of the problem                                                        |
| x         | Current iterate                                                            |
| xa        | Previous iterate                                                           |
| xScal     | Scaling vector                                                             |
| iScal     | Decision parameter for scaling                                             |
| mPr       | Decision parameter for printing                                            |
| printIO   | IO handle for printing                                                     |

## Output parameters
-------------------
| Variables | Description                                                                |
|:----------|:---------------------------------------------------------------------------|
| xw        | Scaling vector computed by this routine<br>All components must be positive.|

"""
function nScal{T}(n::Int64, x::Vector{T}, xa::Vector{T}, xScal::Vector{T},
    iScal::Int64, mPr::Int64, printIO, xw::Vector{T})
    # Begin
    (_, small, _) = getMachineConstants(T)
    if iScal == 1
        xw[:] = xScal
    else
        for l1 = 1:n
            xw[l1] = max(xScal[l1],max((abs(x[l1])+abs(xa[l1]))*0.5,small))
        end
    end

    if mPr >= 6
        write(printIO,"\n\n",
        "+++++++++++++++++++++++++++++++++++++++++++++++++\n",
        "      x-components         scaling-components\n")
        for l1 = 1:n
            write(printIO,"  %18.10e   %18.10e\n",x[l1],xw[l1])
        end
        write(printIO,"+++++++++++++++++++++++++++++++++++++++++++++++++\n\n")
    end
    return nothing
end

"""
function nScrf{T}(m::Int64, n::Int64, a::Array{T,2}, fw::Vector{T})

Row scaling of a (m,n)-matrix in full storage mode

T = Float64 or BigFloat

## Input parameters
-------------------
| Variables | Description                    |
|:----------|:-------------------------------|
| m         | Number of rows of the matrix   |
| n         | Numer of columns of the matrix |
| a[m,n]*   | Matrix to be scaled            |

(* marks inout parameters)

## Output parameters
-------------------
| Variables | Description                    |
|:----------|:-------------------------------|
| fw        | Row scaling factors.           |
"""
function nScrf{T}(m::Int64, n::Int64, a::Array{T,2}, fw::Vector{T})
    # Begin
    if issparse(a)
        nza = nnz(a)
        (row,col) = findn(a)
        aout = sparse(row,col,zeros(nza),m,n)
        for j = 1:nza
            rj = row[j]
            aaj = abs(a[rj,col[j]])
            fw[rj] = max(fw[rj],aaj)
        end
        for k = 1:m
            if fw[k] > 0.0
                fw[k] = 1.0/fw[k]
            else
                fw[k] = 1.0
            end
        end
        for j = 1:nza
            aout[row[j],col[j]] = a[row[j],col[j]]*fw[row[j]]
        end
        a[:,:] = aout
    else
        aout = zeros(a)
        for k = 1:m
            s1 = maximum(abs(a[k,1:n]))
            if s1 > 0.0
                s1 = 1.0/s1
                fw[k] = s1
                aout[k,1:n] = a[k,1:n]*s1
            else
                fw[k] = 1.0
                aout[k,1:n] = a[k,1:n]
            end
        end
        a[:,:] = aout
    end
end

"""
function nScrb{T}(n::Int64, lda::Int64, ml::Int64, mu::Int64, a::Array{T,2},
    fw::Vector{T})

Row scaling of a (n,n)-matrix in band storage mode

T = Float64 or BigFloat

## Input parameters
-------------------
| Variables | Description                              |
|:----------|:-----------------------------------------|
| n         | Number of rows and columns of the matrix |
| lda       | Leading dimension of the matrix array    |
| ml        | Lower bandwidth of the matrix            |
| mu        | Upper bandwidth of the matrix            |
| a[lda,n]* | Matrix to be scaled                      |

(* marks inout parameters)

## Output parameters
-------------------
| Variables | Description                              |
|:----------|:-----------------------------------------|
| fw        | Row scaling factors.                     |
"""
function nScrb{T}(n::Int64, lda::Int64, ml::Int64, mu::Int64, a::Array{T,2},
    fw::Vector{T})
    # Begin
    aout = zeros(a)
    m2 = ml + mu + 1
    for k = 1:n
        s1 = zero(T)
        l2 = max(1,k-ml)
        l3 = min(n,k+mu)
        k1 = m2 + k
        for l1 = l2:l3
            s1 = max(s1,abs(a[k1-l1,l1]))
        end
        if s1 > 0.0
            s1 = 1.0/s1
            fw[k] = s1
            for l1 = l2:l3
                aout[k1-l1,l1] = a[k1-l1,l1]*s1
            end
        else
            fw[k] = one(T)
        end
    end
    a[:,:] = aout
end

"""
function nLvls{T}(n::Int64, dxq::Vector{T}, dx1::Vector{T},
    xw::Vector{T}, f::Vector{T}, qdscal::Bool)

Provides descaled solution, error norm and level functions
To be used in connection with NLEQ1 and NLEQ2.

T = Float64 or BigFloat

## Input parameters
-------------------
| Variables | Description                                   |
|:----------|:----------------------------------------------|
| n         | Number of parameters to be estimated          |
| dx1       | Scaled Newton correction                      |
| xw        | Vector of scaling values                      |
| f         | Residual vector                               |
| qdscal    | true of descaling of dx1 required, else false |

## Output parameters
-------------------
| Variables | Description                              |
|:----------|:-----------------------------------------|
| dxq       | Leading dimension of the matrix array    |
| conv      | Scaled maximum norm of Newton correction |
| sumX      | Scaled natural level function value      |
| dLevF     | Standard level function value            |

"""
function nLvls{T}(n::Int64, dxq::Vector{T}, dx1::Vector{T},
    xw::Vector{T}, f::Vector{T}, qdscal::Bool)
    # Begin
    if qdscal
        # ----------------------------------------------------------------------
        # 1.2 Descaling of solution dx1 (stored to dxq)
        dxq[:] = dx1.*xw
    end
    # --------------------------------------------------------------------------
    # 2 Evaluation of scaled natural level function sumx and scaled maximum
    # error norm conv
    conv = maximum(abs(dx1))
    sumx = sum(dx1.^2)
    # --------------------------------------------------------------------------
    # 3 Evaluation of (scaled) standard level function dlevf
    dlevf = sqrt(sum(f.^2)/n)
    return (conv,sumx,dlevf)
end

# The following function nPrv1 is a multiple dispatch function
# The first one corresponds to nleq1 and the second corresponds to nleq2
"""
function nPrv1{T}(dlevf::T, dlevx::T, fc::T, niter::Int64,
    newt::Int64, mPr::Int64, printIO, qMixIO::Bool)

Printing of intermediate values (Type 1 routine)

T = Float64 or BigFloat

## Parameters
-------------
For all the parameters check n1int
"""
function nPrv1{T}(dlevf::T, dlevx::T, fc::T, niter::Int64,
    newt::Int64, mPr::Int64, printIO, qMixIO::Bool)
    # Begin
    if qMixIO
        write(printIO,"  ******************************************************************",
        "\n");
        if mPr >= 3
            write(printIO,"        It       Normf           Normx                     New\n")
        end
        if mPr == 2
            write(printIO,"        It       Normf           Normx         Damp.Fct.   New\n")
        end
    end
    if mPr >= 3 || niter == 0
        write(printIO,@sprintf("      %4i     %10.3e      %10.3e                 %2i\n",niter,dlevf,dlevx,newt))
    end
    if mPr == 2 && niter != 0
        write(printIO,@sprintf("      %4i     %10.3e      %10.3e      %7.5f    %2i\n",niter,dlevf,dlevx,fc,newt))
    end
    if qMixIO
        write(printIO,"  ******************************************************************",
        "\n")
    end
    flush(printIO)
    return nothing
end

"""
function nPrv1{T}(dlevf::T, dlevx::T, fc::T, niter::Int64, newt::Int64,
    iRank::Int64, mPr::Int64, printIO, qMixIO::Bool, cond1::T)

Printing of intermediate values (Type 1 routine)

T = Float64 or BigFloat

Parameters
-------------
For all the parameters check n2int
"""
function nPrv1{T}(dlevf::T, dlevx::T, fc::T, niter::Int64, newt::Int64,
    iRank::Int64, mPr::Int64, printIO, qMixIO::Bool, cond1::T)
    # Begin
    if qMixIO
        write(printIO,"  ******************************************************************",
        "\n");
        if mPr >= 3
            write(printIO,"        It       Normf           Normx                     New      Rank        Cond\n")
        end
        if mPr == 2
            write(printIO,"        It       Normf           Normx         Damp.Fct.   New      Rank        Cond\n")
        end
    end
    if mPr >= 3 || niter == 0
        write(printIO,@sprintf("      %4i     %10.3e      %10.3e",niter,dlevf,dlevx),
        @sprintf("                 %2i      %4i          %10.3e\n",newt,iRank,cond1))
    end
    if mPr == 2 && niter != 0
        write(printIO,@sprintf("      %4i     %10.3e      %10.3e",niter,dlevf,dlevx),
        @sprintf("      %7.5f    %2i      %4i          %10.3e\n",fc,newt,iRank,cond1))
    end
    if qMixIO
        write(printIO,"  ******************************************************************",
        "\n")
    end
    flush(printIO)
    return nothing
end

"""
function nPrv2{T}(dlevf::T, dlevx::T, fc::T, niter::Int64,
    printIO, qMixIO::Bool, cmark::AbstractString)

Printing of intermediate values (Type 2 routine)

T = Float64 or BigFloat

## Input parameters
-------------------
| Variables | Description                                 |
|:----------|:--------------------------------------------|
| dlevf     | Standard level function value               |
| dlevx     | Standard level value                        |
| fc        | Current damping factor                      |
| niter     | Current number of Newton iterations         |
| qMixIO    | Decision parameter for printing             |
| cmark     | Marker character to be printed before dlevx |
"""
function nPrv2{T}(dlevf::T, dlevx::T, fc::T, niter::Int64,
    printIO, qMixIO::Bool, cmark::AbstractString)
    if qMixIO
        write(printIO,"  ******************************************************************",
        "\n");
        write(printIO,"        It       Normf           Normx         Damp.Fct.\n")
    end
    write(printIO,@sprintf("      %4i     %10.3e    %1s %10.3e      %7.5f\n",niter,dlevf,cmark,dlevx,fc))
    if qMixIO
        write(printIO,"  ******************************************************************",
        "\n")
    end
    flush(printIO)
    return nothing
end

"""
function nSout{T}(n::Int64, x::Vector{T}, mode::Int64, mPr::Int64, printIO,
    nIter::Int64, dLevF::T, sumX::T)

Printing of iterate (user customizable routine)

T = Float64 or BigFloat

## Input parameters
-------------------
| Variables | Description                                                     |
|:----------|:----------------------------------------------------------------|
| n         | Size of the problem                                             |
| x         | Iterate vector                                                  |
| mode = 1  | This routine is called before the first Newton iteration step   |
|      = 2  | This routine is called with an intermediate iterate x           |
|      = 3  | This is the last call with the solution vector x                |
|      = 4  | This is the last call with the final, but not solution vector x |
| mPr       | Decision parameter for printing                                 |
| printIO   | IO handle for printing                                          |
| nIter     | Current number of Newton iterations                             |
| dLevF     | Standard level function value                                   |
| sumX      | Scaled natural level function value                             |
"""
function nSout{T}(n::Int64, x::Vector{T}, mode::Int64, mPr::Int64, printIO,
    nIter::Int64, dLevF::T, sumX::T)
    # Begin
    if mode == 1
        write(printIO,@sprintf("\n%s\n%s%5i\n\n%s\n","  Start data:","  N =",n,
            "  Format: iteration-number, (x(i),i=1,...N) , Normf , Normx "))
        write(printIO,@sprintf("%s\n","  Initial data:"))
    elseif mode == 3
        write(printIO,@sprintf("%s\n","  Solution data:"))
    elseif mode == 4
        write(printIO,@sprintf("%s\n","  Final data:"))
    end
    write(printIO,@sprintf(" %5i\n",nIter))
    l2 = 0
    for l1 = 1:n
        write(printIO,@sprintf("%18.10e ",x[l1]))
        l2 += 1
        if l2 == 3
            write(printIO," \n")
            l2 = 0
        end
    end
    write(printIO,@sprintf("%18.10e %18.10e \n",dLevF,
        sqrt(sumX/n)))
    if mode == 1 && mPr >= 2
        write(printIO,@sprintf("%s\n","  Intermediate data:"))
    elseif mode >= 3
        write(printIO,@sprintf("%s\n","  End data:"))
    end
    flush(printIO)
    return nothing
end

"""
function wnorm{T}(n::Int64, z::Vector{T}, xw::Vector{T})

Return the norm to be used in exit (termination) criteria

T = Float64 or BigFloat

## Input parameters
-------------------
| Variables | Description                                                     |
|:----------|:----------------------------------------------------------------|
| n         | Size of the problem                                             |
| z         | The vector of which norm is to be computed                      |
| xw        | Scaling value of z                                              |

## Output
---------
The mean square root norm of z subject to the scaling values in xw.
"""
function wnorm{T}(n::Int64, z::Vector{T}, xw::Vector{T})
    # Begin
    return sqrt(sum((z./xw).^2)/n)
end


text:
function save(f::Stream{format"PLY_BINARY"}, msh::AbstractMesh)
    io = stream(f)
    points = decompose(Point{3, Float32}, msh)
    faces = decompose(GLTriangleFace, msh)

    n_points = length(points)
    n_faces = length(faces)

    # write the header
    write(io, "ply\n")
    write(io, "format binary_little_endian 1.0\n")
    write(io, "element vertex $n_points\n")
    write(io, "property float x\nproperty float y\nproperty float z\n")
    write(io, "element face $n_faces\n")
    write(io, "property list uchar int vertex_index\n")
    write(io, "end_header\n")

    # write the vertices and faces
    write(io, points)

    for f in faces
        write(io, convert(UInt8, 3))
        write(io, raw.(f)...)
    end
    close(io)
end

function save(f::Stream{format"PLY_ASCII"}, msh::AbstractMesh)
    io = stream(f)
    points = coordinates(msh)
    meshfaces = faces(msh)

    n_faces = length(points)
    n_points = length(meshfaces)

    # write the header
    write(io, "ply\n")
    write(io, "format ascii 1.0\n")
    write(io, "element vertex $n_faces\n")
    write(io, "property float x\nproperty float y\nproperty float z\n")
    write(io, "element face $n_points\n")
    write(io, "property list uchar int vertex_index\n")
    write(io, "end_header\n")

    # write the vertices and faces
    for v in points
        println(io, join(Point{3, Float32}(v), " "))
    end
    for f in meshfaces
        println(io, length(f), " ", join(raw.(ZeroIndex.(f)), " "))
    end
    close(io)
end

function load(fs::Stream{format"PLY_ASCII"}; facetype=GLTriangleFace, pointtype=Point3f0)
    io = stream(fs)
    n_points = 0
    n_faces = 0

    properties = String[]

    # read the header
    line = readline(io)

    while !startswith(line, "end_header")
        if startswith(line, "element vertex")
            n_points = parse(Int, split(line)[3])
        elseif startswith(line, "element face")
            n_faces = parse(Int, split(line)[3])
        elseif startswith(line, "property")
            push!(properties, line)
        end
        line = readline(io)
    end

    faceeltype = eltype(facetype)
    points = Array{pointtype}(undef, n_points)
    #faces = Array{FaceType}(undef, n_faces)
    faces = facetype[]

    # read the data
    for i = 1:n_points
        points[i] = pointtype(parse.(eltype(pointtype), split(readline(io)))) # line looks like: "-0.018 0.038 0.086"
    end

    for i = 1:n_faces
        line = split(readline(io))
        len = parse(Int, popfirst!(line))
        if len == 3
            push!(faces, NgonFace{3, faceeltype}(reinterpret(ZeroIndex{UInt32}, parse.(UInt32, line)))) # line looks like: "3 0 1 3"
        elseif len == 4
            push!(faces, convert_simplex(facetype, QuadFace{faceeltype}(reinterpret(ZeroIndex{UInt32}, parse.(UInt32, line))))...) # line looks like: "4 0 1 2 3"
        end
    end
    return Mesh(points, faces)
end

function load(fs::Stream{format"PLY_BINARY"}; facetype=GLTriangleFace, pointtype=Point3f0)
    io = stream(fs)
    n_points = 0
    n_faces = 0

    properties = String[]

    # read the header
    line = readline(io)

    while !startswith(line, "end_header")
        if startswith(line, "element vertex")
            n_points = parse(Int, split(line)[3])
        elseif startswith(line, "element face")
            n_faces = parse(Int, split(line)[3])
        elseif startswith(line, "property")
            push!(properties, line)
        end
        line = readline(io)
    end

    faceeltype = eltype(facetype)
    points = Array{pointtype}(undef, n_points)
    #faces = Array{FaceType}(undef, n_faces)
    faces = facetype[]

    # read the data
    for i = 1:n_points
        points[i] = pointtype(read(io, Float32), read(io, Float32), read(io, Float32))
    end

    for i = 1:n_faces
        len = read(io, UInt8)
        indices = reinterpret(ZeroIndex{UInt32}, [ read(io, UInt32) for _ in 1:len ])
        if len == 3
            push!(faces, NgonFace{3, faceeltype}(indices)) # line looks like: "3 0 1 3"
        elseif len == 4
            push!(faces, convert_simplex(facetype, QuadFace{faceeltype}(indices))...) # line looks like: "4 0 1 2 3"
        end
    end

    return Mesh(points, faces)
end


text:
# DataLoaders - Basic data loaders for training machine learning
# models
#
# Lucas Ondel 2020

module BasicDataLoaders

using BSON

#######################################################################
# Basic input for loading / saving data, models, ...

export load
export save

include("io.jl")

#######################################################################
# Abstract data loader

export AbstractDataLoader

"""
    abstract type AbstractDataLoader end

Base type for all the data loaders.
"""
abstract type AbstractDataLoader{T} end

# Subtypes should implement the Iteration and Indexing interfaces

#######################################################################
# Concrete data loaders

export DataLoader

include("dataloader.jl")

end

max_stars_repo_name: lucasondel/DataLoaders

max_issues_repo_name: lucasondel/DataLoaders

max_forks_repo_name: lucasondel/DataLoaders


text:
@testset "MvNormalData summary statistics and MLE" begin
    for _ in 1:100
        K = rand(5:10)
        N = K + rand(10:100)
        X = randn(N, K)
        mvn = MvNormalData(X)
        @test mvn.W == N
        @test size(mvn.m) == (K, )
        @test mvn.m ≈ vec(mean(X, 1))
        Xc = X .- mvn.m'
        @test mvn.S ≈ (Xc'*Xc)/N
        @test MLE(MvNormalModel(), mvn) == MvNormalParams(mvn.m, mvn.S)
        @test size(mvn) == (N, K)
        @test firstline(repr(mvn)) ==
            "Summary statistics for multivariate normal, $(N) × $(K) samples"
    end
end

@testset "MvNormalData summary statistics and MLE (weighted)" begin
    for _ in 1:100
        K = rand(5:10)
        N = K + rand(10:100)
        X = randn(N, K)
        wv = Weights(abs.(randn(N)))
        W = sum(wv)
        mvn = MvNormalData(X, wv)
        @test mvn.W == W
        @test size(mvn.m) == (K, )
        @test mvn.m ≈ vec(mean(X, wv, 1))
        Xc = X .- mvn.m'
        @test mvn.S ≈ ((values(wv).*Xc)'*Xc)/W
        @test MLE(MvNormalModel(), mvn) == MvNormalParams(mvn.m, mvn.S)
        @test size(mvn) == (W, K)
        @test firstline(repr(mvn)) ==
            "Summary statistics for multivariate normal, $(sum(wv)) × $(K) samples"
    end
end

@testset "MvNormalData log likelihood" begin
    for _ in 1:100
        K = rand(5:10)
        N = K + rand(10:100)
        X = randn(2*N, K)
        mvnX = MvNormalData(X)
        paramsX = MLE(MvNormalModel(), mvnX)
        Y = randn(N, K)
        mvnY = MvNormalData(Y)
        ℓ = sum(logpdf(MvNormal(paramsX.μ, paramsX.Σ), Y'))
        @test ℓ ≈ loglikelihood(MvNormalModel(), mvnY, paramsX)
    end
end


text:
using ACEatoms
using Test

@testset "ACEatoms.jl" begin

   # ----------------------
   #   pair potentials
   @testset "PolyPairBasis" begin include("pair/test_pair_basis.jl") end
   @testset "PolyPairPot" begin include("pair/test_pair_pot.jl") end
   @testset "RepulsiveCore" begin include("pair/test_repulsion.jl") end

   # ----------------------
   #   species 1p basis
   @testset "ACEatoms 1p Basis" begin include("test_1pbasis.jl")  end
   #   basic site energy calculator
   @testset "ACESitePotential" begin  include("test_siteenergy.jl") end

   # ----------------------
   #  special physics
   @testset "Electrostatics" begin include("test_electro.jl") end
end


text:
# This file was generated, do not modify it. # hide
using RDatasets
using ScientificTypes

boston = dataset("MASS", "Boston")
sch = schema(boston)

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
platforms_and_sources = Dict(
    Platform("x86_64", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/cudnn/secure/8.3.1/local_installers/10.2/cudnn-linux-x86_64-8.3.1.22_cuda10.2-archive.tar.xz",
                      "5982bb96c2a720268fa44b908feb5258d060ad47f1f6e6030e760d13195ea964")],
    Platform("x86_64", "windows") => [
        ArchiveSource("https://developer.nvidia.com/compute/cudnn/secure/8.3.1/local_installers/10.2/cudnn-windows-x86_64-8.3.1.22_cuda10.2-archive.zip",
                      "83797cac5c1e6c59bcf6bcd13ae23a3d00f909b19295819fb50063ea3356fadc")],
)

max_stars_repo_name: AbhimanyuAryan/Yggdrasil

max_issues_repo_name: AbhimanyuAryan/Yggdrasil

max_forks_repo_name: AbhimanyuAryan/Yggdrasil


text:
# Root finder functions for 1d func
# Newton's Method, Linear Approximation
# Note that you can also find root by applying optimizer on an Error function
# ϵ(x,...)

# Newton's Method
# f1 is δf╱δx
function solve_Newton(f,f1,par;nlim = 1e4,ϵ=1e-3)
    for i in 1:nlim
        par1 = par - f(par)/f1(par)
        if ϵ>abs(par1 -par)
            return(par1)
        end
        par=par1
    end
    println("Warning! unconverged after" * string(nlim) *"Iterations.")
    return(par1)
end

# Linear Approximation
function solve_Linear_Approx(f,par1,par2;nlim = 1e4,ϵ=1e-3)
    #must follow :x1 * x2<0
    #             par1<par2
    for i in 1:nlim
        x1 = f(par1)
        x2 = f(par2)
        global par3 = par1+(par2 - par1)*abs(x1/(x1-x2))
        if abs(par3-par1)<ϵ
            if abs(par3 - par2)<ϵ
                return par3
            end
        else
            if (par3*par1)>0
                par1=par3
            else (par3*par2)>0
                par2=par3
            end
        end
    end
    println("Warning! unconverged after"*string(nlim)*"iterations.")
    return par3
end

#----------------------------------------------------------------
using Plots
x = 1:.1:4
y = sin.(x)
plot(x,y)
plot!(x,zeros(length(x)))

f(x) = sin(x)
f1(x) = cos(x)

solve_Newton(f,f1,3)

solve_Linear_Approx(f,3,4)
#----------------------------------------------------------------

max_stars_repo_path: Root Solving.jl

max_issues_repo_path: Root Solving.jl

max_forks_repo_path: Root Solving.jl


text:
#
# server.jl --
#
# Implement XPA client methods.
#
#------------------------------------------------------------------------------
#
# This file is part of XPA.jl released under the MIT "expat" license.
# Copyright (C) 2016-2020, Éric Thiébaut (https://github.com/JuliaAstro/XPA.jl).
#

# We must make sure that the `send` and `recv` callbacks exist during the life
# of the server.  To that end, we use the following dictionary to maintain
# references to callbacks while they are used by an XPA server.
const _SERVERS = Dict{Ptr{Cvoid},Tuple{Union{SendCallback, Nothing},
                                       Union{ReceiveCallback, Nothing}}}()

"""
```julia
XPA.Server(class, name, help, send, recv) -> srv
```

yields an XPA server identified by `class` and `name` (both specified as two
strings).

Argument `help` is a string which is meant to be returned by a help request
from `xpaget`:

```sh
xpaget class:name -help
```

Arguments `send` and `recv` are callbacks which will be called upon a client
[`XPA.get`](@ref) or [`XPA.set`](@ref) respectively.  At most one callback may
be `nothing`.

The send callback will be called in response to an external request from the
`xpaget` program, the `XPAGet()` or `XPAGetFd()` C routines, or the
[`XPA.get`](@ref) Julia method.  This callback is used to send data to the
requesting client and is a combination of a function (`sfunc` below) and
private data (`sdata` below) as summarized by the following typical example:

```julia
# Method to handle a send request:
function sfunc(sdata::S, srv::XPA.Server, params::String, buf::XPA.SendBuffer)
    result = ... # build up the result of the request
    try
        XPA.store!(buf, result)
        return XPA.SUCCESS
    catch err
        error(srv, err)
        return XPA.FAILURE
    end
end

# A send callback combines a method and some contextual data:
send = XPA.SendCallback(sfunc, sdata)
```

Here `sdata` is the client data (of type `S`) of the send callback, `srv` is
the XPA server serving the request, `params` is the parameter list of the
`[XPA.get](@ref)` call and `buf` specifies the addresses where to store the
result of the request and its size (in bytes).

The receive callback will be called in response to an external request from the
`xpaset` program, the `XPASet()` or `XPASetFd()` C routines, or the
[`XPA.set`](@ref) Julia method.  This callback is used to process sent data to
the requesting client and is a combination of a function (`rfunc` below) and
private data (`rdata` below) as summarized by the following typical example:

```julia
# Method to handle a send request:
function rfunc(rdata::R, srv::XPA.Server, params::String, buf::XPA.ReceiveBuffer)
    println("receive: \$params")
    # Temporarily wrap the received data into an array.
    bytes = XPA.peek(Vector{UInt8}, buf; temporary=true)
    ... # process the received bytes
    return XPA.SUCCESS
end

# A receive callback combines a method and some contextual data:
send = XPA.ReceiveCallback(rfunc, rdata)
```

Here `rdata` is the client data (of type `R`) of the receive callback, `srv` is
the XPA server serving the request, `params` is the parameter list of the
[`XPA.set`](@ref) call and `buf` specifies the address and size of the data to
process.

The callback methods `sfunc` and/or `rfunc` should return [`XPA.SUCCESS`](@ref)
if no error occurs, or [`XPA.FAILURE`](@ref) to signal an error.  The Julia XPA
package takes care of maintaining a reference on the client data and callback
methods.

See also [`XPA.poll`](@ref), [`XPA.mainloop`](@ref), [`XPA.store!`](@ref),
[`XPA.SendCallback`](@ref), [`XPA.ReceiveCallback`](@ref) and
[`XPA.peek`](@ref).

"""
function Server(class::AbstractString,
                name::AbstractString,
                help::AbstractString,
                send::Union{SendCallback, Nothing},
                recv::Union{ReceiveCallback, Nothing})
    # Create an XPA server and a reference to the callback objects to make
    # sure they are not garbage collected while the server is running.
    server = Server(class, name, help,
                    _callback(send), _context(send), _mode(send),
	            _callback(recv), _context(recv), _mode(recv))
    _SERVERS[server.ptr] = (send, recv)
    return server
end

function Server(class::AbstractString, name::AbstractString,
                help::AbstractString,
                sproc::Ptr{Cvoid}, sdata::Ptr{Cvoid}, smode::AbstractString,
                rproc::Ptr{Cvoid}, rdata::Ptr{Cvoid}, rmode::AbstractString)
    ptr = ccall((:XPANew, libxpa), Ptr{Cvoid},
                (Cstring, Cstring, Cstring,
	         Ptr{Cvoid}, Ptr{Cvoid}, Cstring,
	         Ptr{Cvoid}, Ptr{Cvoid}, Cstring),
                class, name, help,
                sproc, sdata, smode,
                rproc, rdata, rmode)
    ptr != C_NULL || error("failed to create an XPA server")
    obj = finalizer(close, Server(ptr))
    (get_send_mode(obj) & MODE_FREEBUF) != 0 ||
        error("send mode must have `freebuf` option set")
    return obj
end

# The following methods are helpers to build instances of an XPA server.
_callback(::Nothing) = C_NULL
_callback(::SendCallback) = _SEND_REF[]
_callback(::ReceiveCallback) = _RECV_REF[]
_context(::Nothing) = C_NULL
_context(cb::Callback) = pointer_from_objref(cb)
_mode(::Nothing) = ""
_mode(cb::SendCallback) = "acl=$(cb.acl),freebuf=true"
_mode(cb::ReceiveCallback) =
    "acl=$(cb.acl),buf=true,fillbuf=true,freebuf=true"

# The following method is called upon garbage collection of an XPA server.
function Base.close(srv::Server)
    if (ptr = srv.ptr) != C_NULL
        srv.ptr = C_NULL # avoid closing more than once!
        ccall((:XPAFree, libxpa), Cint, (Ptr{Cvoid},), ptr)
        haskey(_SERVERS, ptr) && pop!(_SERVERS, ptr)
    end
    return nothing
end

"""
```julia
SendCallback(func, data=nothing; acl=true)
```

yields an instance of `SendCallback` for sending the data requested by a call
to [`XPA.get`](@ref) (or similar) to an XPA server.  Argument `func` is the
method to be called to process the request and optional argument `data` is some
associated contextual data.

Keyword `acl` can be used to specify whether access control is enabled (true by
default).

!!! note
    The `freebuf` option is not available because we are always assuming that
    the answer to a [`XPA.get`](@ref) request is a dynamically allocated buffer
    which is automatically deleted by the XPA library.  This is like imposing
    that the `freebuf` option is laways true.  This choice has been made
    because it would otherwise be difficult to warrant that data passed by a
    Julia send callback be not garbage collected before being fully transfered
    to the client.

See also [`XPA.Server`](@ref), [`XPA.store!`](@ref) and
[`XPA.ReceiveCallback`](@ref).

"""
function SendCallback(func::F,
                      data::T = nothing;
                      acl::Bool = true) where {T,F<:Function}
    return SendCallback{T,F}(func, data, acl)
end

"""
```julia
ReceiveCallback(func, data=nothing; acl=true)
```

yields an instance of `ReceiveCallback` for processing the data sent by a call
to [`XPA.set`](@ref) (or similar) to an XPA server.  Argument `rfunc` is the
method to be called to process the request and optional argument `data` is some
associated contextual data.

Keyword `acl` can be used to specify whether access control is enabled (true by
default).

!!! note
    The `buf`, `fillbuf` and `freebuf` options are not available because we are
    always assuming that the data buffer accompanying an [`XPA.set`](@ref)
    request is always provided as a dynamically allocated buffer by the XPA
    library.  This is like imposing that the `buf`, `fillbuf` and `freebuf`
    options are always true.  This choice has been made because it would
    otherwise be difficult to warrant that data passed to a Julia receive
    callback can be safely stealed by Julia.

Also see [`XPA.Server`](@ref), [`XPA.SendCallback`](@ref) and
[`XPA.set`](@ref).

"""
function ReceiveCallback(func::F,
                         data::T = nothing;
                         acl::Bool = true) where {T,F<:Function}
    return ReceiveCallback{T,F}(func, data, acl)
end

const _MINIMAL_SEND_MODE = MODE_FREEBUF
const _MINIMAL_RECEIVE_MODE = (MODE_BUF | MODE_FILLBUF | MODE_FREEBUF)

# The send callback is executed in response to an external request from the
# `xpaget` program, the `XPAGet()` routine, or `XPAGetFd()` routine.
function _send(clientdata::Ptr{Cvoid}, handle::Ptr{Cvoid}, params::Ptr{Byte},
               bufptr::Ptr{Ptr{Byte}}, lenptr::Ptr{Csize_t})::Cint
    # Check assumptions.
    srv = Server(handle)
    (get_send_mode(srv) & _MINIMAL_SEND_MODE) == _MINIMAL_SEND_MODE ||
        return error(srv, "send mode must have option `freebuf=true`")

    # Call actual callback providing the client data is the address of a known
    # SendCallback object.
    return _send(unsafe_pointer_to_objref(clientdata), srv,
                 (params == C_NULL ? "" : unsafe_string(params)),
                 SendBuffer(bufptr, lenptr))
end

_send(cb::SendCallback, srv::Server, params::String, buf::SendBuffer) =
    cb.send(cb.data, srv, params, buf)

# The receive callback is executed in response to an external request from the
# `xpaset` program, the `XPASet()` routine, or `XPASetFd()` routine.
function _recv(clientdata::Ptr{Cvoid}, handle::Ptr{Cvoid}, params::Ptr{Byte},
               buf::Ptr{Byte}, len::Csize_t)::Cint
    # Check assumptions.
    srv = Server(handle)
    (get_recv_mode(srv) & _MINIMAL_RECEIVE_MODE) == _MINIMAL_RECEIVE_MODE ||
        return error(srv, "receive mode must have options `buf=true`, `fillbuf=true` and `freebuf=true`")

    # Call actual callback providing the client data is the address of a known
    # ReceiveCallback object.
    return _recv(unsafe_pointer_to_objref(clientdata), srv,
                 (params == C_NULL ? "" : unsafe_string(params)),
                 ReceiveBuffer(buf, len))
end

_recv(cb::ReceiveCallback, srv::Server, params::String, buf::ReceiveBuffer) =
    cb.recv(cb.data, srv, params, buf)

# Addresses of callbacks cannot be precompiled so we set them at run-time in
# the __init__() method of the module.
const _SEND_REF = Ref{Ptr{Cvoid}}(0)
const _RECV_REF = Ref{Ptr{Cvoid}}(0)
function __init__()
    global _SEND_REF, _RECV_REF
    _SEND_REF[] = @cfunction(_send, Cint,
                             (Ptr{Cvoid},     # client_data
                              Ptr{Cvoid},     # call_data
                              Ptr{Byte},      # paramlist
                              Ptr{Ptr{Byte}}, # buf
                              Ptr{Csize_t}))  # len
    _RECV_REF[] = @cfunction(_recv, Cint,
                             (Ptr{Cvoid},     # client_data
                              Ptr{Cvoid},     # call_data
                              Ptr{Byte},      # paramlist
                              Ptr{Byte},      # buf
                              Csize_t))       # len
end

"""
```julia
error(srv, msg) -> XPA.FAILURE
```

communicates error message `msg` to the client when serving a request by XPA
server `srv`.  This method shall only be used by the send/receive callbacks of
an XPA server.

Also see: [`XPA.Server`](@ref), [`XPA.message`](@ref),
          [`XPA.SendCallback`](@ref), [`XPA.ReceiveCallback`](@ref).

"""
function Base.error(srv::Server, msg::AbstractString)
    ccall((:XPAError, libxpa), Cint, (Server, Cstring),
          srv, msg) == SUCCESS ||
              error("XPAError failed for message \"$msg\"");
    return FAILURE
end

"""
```julia
XPA.message(srv, msg)
```

sets a specific acknowledgment message back to the client. Argument `srv` is
the XPA server serving the client and `msg` is the acknowledgment message.
This method shall only be used by the receive callback of an XPA server.

Also see: [`XPA.Server`](@ref), [`XPA.error`](@ref),
          [`XPA.ReceiveCallback`](@ref).

"""
message(srv::Server, msg::AbstractString) =
    ccall((:XPAMessage, libxpa), Cint, (Server, Cstring), srv, msg)

"""
```julia
XPA.store!(buf, data)
```

or

```julia
XPA.store!(buf, ptr, len)
```

store into the send buffer `buf` a dynamically allocated copy of the contents
of `data` or of the `len` bytes at address `ptr`.

!!! warning
    This method is meant to be used in a *send* callback to store the result of
    an [`XPA.get`](@ref) request processed by an XPA server.  Memory leaks are
    expected if used in another context.

See also [`XPA.Server`](@ref), [`XPA.SendCallback`](@ref) and
[`XPA.get`](@ref).

"""
function store!(buf::SendBuffer, ptr::Ptr{Byte}, len::Integer)
    # Before calling the send callback (see xpa.c), the buffer is empty
    # (*bufptr = NULL and *lenptr = 0).  On return of the send callback with a
    # successful status, if there are any data (*bufptr != NULL and *lenptr >
    # 0), this data is sent to the client.  Then, whatever the status and if
    # freebuf is true, the data buffer is destroy with free() or any specific
    # function set with XPASetFree().
    #
    # This function is similar to `XPASetBuf` except that it makes a dynamic
    # copy of the data to send (because option `freebuf` is always true) and
    # destroys any buffer which could have been set before.
    if (tmp = unsafe_load(buf.bufptr)) != NULL && unsafe_load(buf.lenptr) > 0
        unsafe_store!(buf.lenptr, 0)
        unsafe_store!(buf.bufptr, NULL)
        _free(tmp)
    end
    if ptr != NULL
        len > 0 || error("invalid number of bytes ($len) for non-NULL pointer")
        unsafe_store!(buf.bufptr, _memcpy(_malloc(len), ptr, len))
        unsafe_store!(buf.lenptr, len)
    else
        len == 0 || error("invalid number of bytes ($len) for NULL pointer")
    end
    return nothing
end

store!(buf::SendBuffer, ::Nothing) = store!(buf, NULL, 0)

store!(buf::SendBuffer, val::Union{Symbol,AbstractString}) =
    store!(buf, String(val))

store!(buf::SendBuffer, str::String) =
    store!(buf.bufptr, buf.lenptr,
           Base.unsafe_convert(Ptr{Byte}, str), sizeof(str))

function store!(buf::SendBuffer, arr::DenseArray{T,N}) where {T, N}
    @assert isbitstype(T)
    store!(buf.bufptr, buf.lenptr,
           convert(Ptr{Byte}, pointer(arr)), sizeof(arr))
end

function store!(buf::SendBuffer, val::T) where {T}
    @assert isbitstype(T)
    if (tmp = unsafe_load(buf.bufptr)) != NULL && unsafe_load(buf.lenptr) > 0
        unsafe_store!(buf.lenptr, 0)
        unsafe_store!(buf.bufptr, NULL)
        _free(tmp)
    end
    len = sizeof(T)
    ptr = _malloc(len)
    unsafe_store!(convert(Ptr{T}, ptr), val)
    unsafe_store!(buf.bufptr, ptr)
    unsafe_store!(buf.lenptr, len)
    return nothing
end

function Base.copyto!(dst::AbstractArray{T,N},
                      buf::ReceiveBuffer) where {T,N}
    @assert isbitstype(T)
    sizeof(T)*length(dst) ≤ sizeof(buf) ||
        error("data buffer is too small for array")
    return copyto!(dst, unsafe_wrap(Array, Ptr{T}(pointer(buf)), size(dst),
                                    own=false))
end

"""
```julia
XPA.peek(T, buf, i=1) -> val
```

yields the `i`-th binary value of type `T` stored into receive buffer `buf`.
Bounds checking is performed unless `@inbounds` is active.

Another usage of the `XPA.peek` method is to *convert* the contents of the
receive buffer into an array:

```julia
XPA.peek(Vector{T}, [len,] buf) -> vec
XPA.peek(Array{T[,N]}, (dim1, ..., dimN), buf) -> arr
```

yield a Julia vector `vec` or array `arr` whose elements are of type `T` and
dimensions are `len` or `(dim1, ..., dimN)`.  For a vector, if the length is
unspecified, the vector of maximal length that fits in the buffer is returned.

If keyword `temporary` is true, then `unsafe_wrap` is called (with option
`own=false`) to wrap the buffer contents into a Julia array whose life-time
cannot exceeds that of the callback.  Otherwise, a copy of the buffer contents
is returned.

See also [`XPA.ReceiveCallback`](@ref).

"""
@inline function peek(::Type{T},
                      buf::ReceiveBuffer,
                      i::Int=1)::T where {T}
    @assert isbitstype(T)
    @boundscheck checkbounds(T, buf, i)
    return unsafe_load(Ptr{T}(pointer(buf)), j)
end

@inline @propagate_inbounds function peek(::Type{T},
                                          buf::ReceiveBuffer,
                                          i::Integer) where {T}
    return peek(T, buf, Int(i))
end

Base.checkbounds(::Type{T}, buf::ReceiveBuffer, i::Int=1) where {T} =
    1 ≤ i && i*sizeof(T) ≤ sizeof(buf) || error("out of range index")

function peek(::Type{Vector{T}},
              buf::ReceiveBuffer;
              temporary::Bool=false)::Vector{T} where {T}
    @assert isbitstype(T)
    local vec::Vector{T}
    len = div(sizeof(buf), sizeof(T))
    if temporary
        vec = unsafe_wrap(Array, Ptr{T}(pointer(buf)), len)
    else
        vec = Vector{T}(undef, len)
        nbytes = sizeof(T)*len
        nbytes > 0 && _memcpy(pointer(vec), pointer(buf), nbytes)
    end
    return vec
end

function peek(::Type{Vector{T}},
              len::Integer,
              buf::ReceiveBuffer;
              kwds...)::Vector{T} where {T}
    return peek(Vector{T}, _dimensions(len), buf; kwds...)
end

function peek(::Type{Array{T}},
              dims::NTuple{N,Integer},
              buf::ReceiveBuffer;
              kwds...)::Array{T,N} where {T,N}
    return peek(Array{T,N}, _dimensions(dims), buf; kwds...)
end

function peek(::Type{Array{T,N}},
              dims::NTuple{N,Integer},
              buf::ReceiveBuffer;
              kwds...)::Array{T,N} where {T,N}
    return peek(Array{T,N}, _dimensions(dims), buf; kwds...)
end

function peek(::Type{Array{T,N}},
              dims::NTuple{N,Int},
              buf::ReceiveBuffer;
              temporary::Bool=false)::Array{T,N} where {T,N}
    @assert isbitstype(T)
    local arr::Array{T,N}
    len = 1
    @inbounds for d in 1:N
        (dim = dims[d]) < 0 && error("invalid Array dimension")
        len *= dim
    end
    nbytes = sizeof(T)*len
    nbytes ≤ sizeof(buf) || error("data buffer is too small for array")
    if temporary
        arr = unsafe_wrap(Array, Ptr{T}(pointer(buf)), dims)
    else
        arr = Array{T,N}(undef, dims)
        nbytes > 0 && _memcpy(pointer(arr), pointer(buf), nbytes)
    end
    return arr
end

"""
```julia
XPA.poll(sec, maxreq)
```

polls for XPA events.  This method is meant to implement a polling event loop
which checks for and processes XPA requests without blocking.

Argument `sec` specifies a timeout in seconds (rounded to millisecond
precision).  If `sec` is positive, the method blocks no longer than this amount
of time.  If `sec` is strictly negative, the routine blocks until the occurence
of an event to be processed.

Argument `maxreq` specifies how many requests will be processed.  If `maxreq <
0`, then no events are processed, but instead, the returned value indicates the
number of events that are pending.  If `maxreq == 0`, then all currently
pending requests will be processed.  Otherwise, up to `maxreq` requests will be
processed.  The most usual values for `maxreq` are `0` to process all requests
and `1` to process one request.

The following example implements a polling loop which has no noticeable impact
on the consumption of CPU when no requests are emitted to the server:

```julia
const __running = Ref{Bool}(false)

function run()
    global __running
    __running[] = true
    while __running[]
        XPA.poll(-1, 1)
    end
end
```

Here the global variable `__running` is a reference to a boolean whose value
indicates whether to continue to run the XPA server(s) created by the process.
The idea is to pass the reference to the callbacks of the server (as their
client data for instance) and let the callbacks stop the loop by setting the
contents of the reference to `false`.

Another possibility is to use [`XPA.mainloop`](@ref) (which to see).

To let Julia performs other tasks, the polling method may be repeatedly called
by a Julia timer.  The following example does this.  Calling `resume` starts
polling for XPA events immediately and then every 100ms.  Calling `suspend`
suspends the processing of XPA events.

```julia
const __timer = Ref{Timer}()

ispolling() = (isdefined(__timer, 1) && isopen(__timer[]))

resume() =
    if ! ispolling()
        __timer[] = Timer((tm) -> XPA.poll(0, 0), 0.0, interval=0.1)
    end

suspend() =
    ispolling() && close(__timer[])
```

Also see: [`XPA.Server`](@ref), [`XPA.mainloop`](@ref).

"""
poll(sec::Real, maxreq::Integer) =
    ccall((:XPAPoll, libxpa), Cint, (Cint, Cint),
          (sec < 0 ? -1 : round(Cint, 1E3*sec)), maxreq)

"""
```julia
XPA.mainloop()
```

runs XPA event loop which handles the requests sent to the server(s) created by
this process.  The loop runs until all servers created by this process have
been closed.

In the following example, the receive callback function close the server when
it receives a `"quit"` command:

```julia
function rproc(::Nothing, srv::XPA.Server, params::String,
               buf::Ptr{UInt8}, len::Integer)
    status = XPA.SUCCESS
    if params == "quit"
        close(srv)
    elseif params == ...
        ...
    end
    return status
end
```

Also see: [`XPA.Server`](@ref), [`XPA.mainloop`](@ref).

"""
mainloop() =
    ccall((:XPAMainLoop, libxpa), Cint, ())


text:
@testset "parse single file" begin
    file_path = joinpath(
        @__DIR__,
        "data",
        "bids_root",
        "sub-subtest",
        "ses-1",
        "test",
        "sub-test_ses-1_run-001_modlbl.nii.gz"
    )

    my_file = File(file_path)

    metadata = OrderedDict(
            "key1" => "value1",
            "key2" => "value2"
        )
    entities = OrderedDict(
            "sub" => "test",
            "ses" => "1",
            "run" => "001",
            "modality" => "modlbl"
        )
    @test my_file.path == file_path
    @test my_file.metadata == metadata
    @test my_file.entities == entities

    metadata_path = joinpath(
        @__DIR__,
        "data",
        "bids_root",
        "sub-subtest",
        "ses-1",
        "test",
        "sub-test_ses-1_run-001_modlbl.json"
    )

    @test get_metadata_path(my_file) == metadata_path
    @test get_sub(my_file) == "test"
    @test get_ses(my_file) == "1"
    @test get_sub(
        dirname(my_file.path)*Base.Filesystem.path_separator, from_fname=false
    ) == "subtest"
    @test get_ses(
        dirname(my_file.path)*Base.Filesystem.path_separator, from_fname=false
    ) == "1"

    new_path = joinpath(dirname(file_path), "run-001_modlbl.nii.gz")
    mv(file_path, new_path)
    my_file = File(new_path, load_metadata=false)
    @test my_file.path == new_path
    @test my_file.metadata == OrderedDict{String,Any}()
    entities["sub"] = "subtest"
    @test my_file.entities == entities
    mv(new_path, file_path)
end

@testset "parse filename" begin
    fname = "sub-test_ses-1_run-001_modlbl.json"
    entities = OrderedDict(
        "sub" => "test",
        "ses" => "1",
        "run" => "001",
        "modality" => "modlbl"
    )
    @test parse_fname(fname) == entities

    # Fail tests
    fname_err = [
        (
            "sub-subtest_ses-1_run-001_key-nomoderror.nii.gz",
            AssertionError(
                """Got unexpected modality SubString{String}["key", "nomoderror"]"""
            )
        ),
        (
            "sub-subtest_ses-1-more-than-one-dash_run-001_mod.nii.gz",
            ErrorException(
                """Invalid BIDS file name sub-subtest_ses-1-more-than-one-dash_run-001_mod.nii.gz
                     (part ses-1-more-than-one-dash should have exactly one '-')"""
            )
        ),
        (
            "sub-subtest_key_without_val_mod.nii.gz",
            ErrorException(
                """Invalid BIDS file name sub-subtest_key_without_val_mod.nii.gz
                     (part key should have exactly one '-')"""
            )
        ),
        (
            "sub-subtest_key1-val1_key1-val2_mod.nii.gz",
            ErrorException("Invalid BIDS file name (key key1 occurs twice)")
        ),
        (
            "sub-subtest_-val1_mod.nii.gz",
            ErrorException("Empty key in pair -val1")
        )
    ]
    for (fname_test, e) in fname_err
        @test_throws e parse_fname(fname_test)
    end

    @test construct_fname(entities, ext="json") == fname
end


text:
VERSION < v"0.7.0-beta2.199" && __precompile__()
module Roots

if VERSION >= v"0.7-"
    using Printf
end

using Compat: @nospecialize, lastindex, range, Nothing

export fzero,
       fzeros,
       secant_method

export find_zero, find_zeros,
       Order0, Order1, Order2, Order5, Order8, Order16

export Bisection, FalsePosition

## load in files
include("utils.jl")
include("find_zero.jl")
include("bracketing.jl")
include("derivative_free.jl")
include("simple.jl")
include("find_zeros.jl")
include("newton.jl")
include("fzero.jl")

end


text:
## Exercise 8-8
## There is a built-in function called count that is similar to the function in Looping and Counting. Read the documentation of this function and use it to count the number of a’s in "banana".
println("Ans: ")

println("Please enter something: ")
input = strip(readline())

println(count("a", input))

println("End.")


text:
# DB related

getdb(dbname) = SQLite.DB(dbname)
get_qids(v) = @_ map(_.question_id, v)
clean(questions) = @_ filter("julia" in _.tags, questions.items) |> filter(_.question_id > 0, __)

function show_query(db, query; params = (), n = 1)
    stmt = SQLite.Stmt(db, query)
    res = DBInterface.execute(stmt, params)
    for (i, row) in enumerate(res)
        println(row)
        i >= n && break
    end
end

function print_href(io, text, href)
    print(io, "[", text, "](", href, ")")
end

function print_unixtime(io, tm)
    print(io, "**[")
    print(io, Dates.format(Dates.unix2datetime(tm), "dd u yyyy HH:MM:SS"))
    print(io, "]**")
end

function process_post(x, body = nothing)
    html_body = @_ parsehtml(x.body) |> __.root |> matchFirst(sel"body", __)
    io = IOBuffer()
    ownername = x.owner.display_name
    ownerlink = x.owner.link
    print_href(io, ownername, ownerlink)
    print(io, " asks: ")
    print_href(io, x.title, x.link)
    print(io, "\n")
    print_unixtime(io, x.creation_date)
    println(io, "\t**Score: $(x.score)\tAnswers: $(x.answer_count)**\n")
    if isnothing(body)
        return @_ process_post_body(html_body, io) |> String(take!(__))
    else
        print(io, body)
        return String(take!(io))
    end
end

function process_post_body(x::HTMLText, io, list)
    print(io, x.text)

    return io
end

function process_post_body(x, io = IOBuffer(), list = "")
    if tag(x) == :pre
        # It should be large code block
        print(io, "```\n", nodeText(x), "```\n")
    elseif tag(x) == :code
        # This is inline code
        print(io, "`", nodeText(x), "`")
    elseif tag(x) == :a
        # Link
        print(io, "[", nodeText(x), "](", x.attributes["href"], ")")
    elseif tag(x) == :p
        # Paragraph
        @_ foreach(process_post_body(_, io, list), x.children)
        print(io, "\n")
    elseif tag(x) == :ul
        @_ foreach(process_post_body(_, io, "ul"), x.children)
        print(io, "\n")
    elseif tag(x) == :li
        # In the future, here should be if/else for different kinds of list (numbered, unnumbered, nested)
        print(io, "+ ")
        @_ foreach(process_post_body(_, io, ""), x.children)
        print(io, "\n")
    elseif tag(x) == :blockquote
        # Quote
        print(io, "```quote\n")
        # Bad thing is we are losing all formatting, but let's hope that it is not going to be an issue in most cases
        print(io, nodeText(x))
        print(io, "\n```\n")
    elseif tag(x) == :strong
        print(io, "**")
        @_ foreach(process_post_body(_, io, list), x.children)
        print(io, "**")
    else
        # Not much left to do
        @_ foreach(process_post_body(_, io, list), x.children)
    end

    return io
end

function process_question(item, db; zulip = ZulipGlobal[], to = "stackoverflow", type = "stream")
    status, msg_id, title = qstatus(db, item)
    status == "known" && return

    zulip_msg = process_post(item)
    if length(zulip_msg) > 9900
        zulip_msg = process_post(item, "Message is too long, please read it on [stackoverflow.com]($(item.link))")
    end
    if status == "new"
        try
            res = sendMessage(zulip; to = to, type = type, content = zulip_msg, topic = item.title)
            if get(res, :result, "fail") == "success"
                addquestion!(db, item, res)
            else
                @error "Get bad response from zulip server: $res"
            end
        catch err
            @info "Received error from the server"
            @info zulip_msg
            @info item.title
            @info item
            rethrow()
        end
    else # status == "update"
        res = updateMessage(zulip, msg_id; to = to, type = type, content = zulip_msg, topic = title)
        if get(res, :result, "fail") == "success"
            updquestion!(db, item)
        else
            @error "Get bad response from zulip server: $res"
        end
    end
    @info "Processed question $(item.question_id)"
end

##########################################
# Answers

function msg(x::Answer, body = nothing)
    html_body = @_ parsehtml(x.body) |> __.root |> matchFirst(sel"body", __)
    io = IOBuffer()
    if x.is_accepted
        print(io, "✅ ")
    end
    print_href(io, x.owner.display_name, x.owner.link)
    print(io, " ")
    print_href(io, "answered", "https://stackoverflow.com/a/$(x.answer_id)")
    print(io, ":\n")
    print_unixtime(io, x.creation_date)
    println(io, " **Score: $(x.score)**\n")
    if isnothing(body)
        return @_ process_post_body(html_body, io) |> String(take!(__))
    else
        print(io, body)
        return String(take!(io))
    end
end

function process(answer::Answer, db; zulip = ZulipGlobal[], to = "stackoverflow", type = "stream")
    status, msg_id, title = astatus(db, answer)
    if (isempty(title)) & (status != "new")
        # This is some sort of error
        # It should never happen
        @error "Empty title of the question: " * string(answer.question_id)
        return nothing
    end
    status == "known" && return nothing

    zulip_msg = msg(answer)
    if length(zulip_msg) > 9900
        zulip_msg = msg(answer, "Message is too long, please read it on [stackoverflow.com](https://stackoverflow.com/a/$(answer.answer_id))")
    end
    if status == "new"
        res = sendMessage(zulip; to = to, type = type, content = zulip_msg, topic = title)
        if get(res, :result, "fail") == "success"
            add!(db, answer, res)
        else
            @error "Get bad response from zulip server: $res"
        end
    else # status == "update"
        res = updateMessage(zulip, msg_id; to = to, type = type, content = zulip_msg, topic = title)
        if get(res, :result, "fail") == "success"
            update!(db, answer)
        else
            @error "Get bad response from zulip server: $res"
        end
    end
    @info "Processed answer_id $(answer.answer_id)"
end

function process_answers(answers, db; zulip = ZulipGlobal[], to = "stackoverflow", type = "stream")
    for x in answers
        for answer in x.items
            process(answer, db, zulip = zulip, to = to, type = type)
        end
    end
end

########################################
# Cleanup section
########################################

function remove_question(db, qid; zulip = ZulipGlobal[])
    msgid = get_question(db, qid)
    answids = get_answers(db, qid)
    if isempty(msgid)
        @info "Question $qid is not found"
        return
    end
    msgid = msgid[1]
    try
        deleteMessage(zulip, msgid)
    catch err
        @error err
    end
    for answid in answids
        try
            deleteMessage(zulip, answid)
        catch err
            @error err
        end
    end
    drop_question!(db, qid)
end

function prune_questions(db, qs; zulip = ZulipGlobal[], tagword = "julia")
    qids = @_ map(_.qid, qs)
    soqs = getquestions(qids; site = "stackoverflow", pagesize = length(qids))
    allowed = @_ filter(tagword in _.tags, soqs.items) |> map(_.question_id, __)
    for q in qs
        if !(q.qid in allowed)
            @info "Removing question '$(q.title)'"
            remove_question(db, q.qid, zulip = zulip)
        end
    end
end


text:
##

function benchmark(optimizer, f::BBOBFunction, run_lengths, Ntrials, dimensions, Δf)

    @info("$(string(optimizer))\t $f")

    t = (T) -> NamedDimsArray{(:trial, :run_length, :dimension)}(
        zeros(T, Ntrials, length(run_lengths), length(dimensions))
    )

    reached_minium = t(Bool)
    distance_to_xopt = t(Float64)
    fmin = t(Float64)

    t = 0.0
    for j in 1:length(run_lengths), k in 1:length(dimensions)
        for i in 1:Ntrials
            try
                D = dimensions[k]
                t += @elapsed mfit = optimize(optimizer, f.f, D, run_lengths[j])

                reached_minium[i, j, k] = minimum(mfit) < Δf + f.f_opt
                fmin[i, j, k] = minimum(mfit) - f.f_opt
                distance_to_xopt[i, j, k] =  √sum(abs2.(minimizer(mfit) - f.x_opt[1:D]))

            catch err
                @warn(err)

                reached_minium[i, j, k] = false
                fmin[i, j, k] = NaN
                distance_to_xopt[i, j, k] = NaN
                @warn(string(optimizer, " failed :", err))
            end
        end
    end
    t /= Ntrials*length(run_lengths)*length(dimensions)

    (
     success_counts = sum(reached_minium, dims=1),
     distance_to_minimizer = mean(distance_to_xopt, dims=1),
     minimum = mean(fmin, dims=1),
     runtime = t
    )
end

#
function benchmark(optimizer, funcs::Vector, run_lengths, Ntrials, dimensions, Δf)

    res = [benchmark(optimizer, f, run_lengths, Ntrials, dimensions, Δf) for f in funcs]
    reduce_res(res, field, f=mean) = f(getfield(r, field) for r in res)

    (
        success_counts = reduce_res(res, :success_counts, sum),
        distance_to_minimizer = reduce_res(res, :distance_to_minimizer),
        minimum = reduce_res(res, :minimum),
        runtime = reduce_res(res, :runtime, sum),
    )
end

##


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Tar_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Tar")
JLLWrappers.@generate_main_file("Tar", UUID("9b64493d-8859-5bf3-93d7-7c32dd38186f"))
end  # module Tar_jll


text:
import UncertainData:
    UncertainDataset,
    resample,
    UncertainValue,
    StrictlyDecreasing,
    NoConstraint,
    TruncateLowerQuantile,
    TruncateUpperQuantile,
    TruncateQuantiles,
    TruncateMinimum,
    TruncateMaximum,
    TruncateRange,
    TruncateStd

import Distributions:
    Uniform,
    Normal

import Test

# Create some uncertain data with decreasing magnitude and zero overlap between values,
# so we're guaranteed that a strictly decreasing sequence through the dataset exists.
N = 10
u_timeindices = [ i <= N/2 ? CertainValue(float(i)) : UncertainValue(Normal, i, 1) for i = N:-1:1]
u = UncertainDataset(u_timeindices)
UI = UncertainIndexDataset(u_timeindices)
UV = UncertainValueDataset(u_timeindices)

# No further constraints other than the order constraint
x = resample(UI, StrictlyDecreasing())
@test x isa Vector{Float64}

n_realizations = 100
X = [resample(UI, StrictlyDecreasing()) for i = 1:n_realizations]

# We're getting vectors
@test all([x isa Vector{Float64} for x in X])

# All sequences
@test all([all(diff(x) .< 0) for x in X])

test_constraints = [
    NoConstraint(),
    TruncateLowerQuantile(0.2),
    TruncateUpperQuantile(0.2),
    TruncateQuantiles(0.2, 0.8),
    TruncateMaximum(50),
    TruncateMinimum(-50),
    TruncateRange(-50, 50),
    TruncateStd(1)
]

# First constrain using a single regular constraint, then apply the order constraint.
for i = 1:length(test_constraints)
    constraint = test_constraints[i]
    @test resample(UI, constraint, StrictlyDecreasing()) isa Vector{Float64}
    @test all([resample(UI, constraint, StrictlyDecreasing()) isa Vector{Float64} for k = 1:n_realizations])
end

#First element-wise apply a vector of regular constraints to each element in the dataset,
#then apply the order constraint.
for i = 1:length(test_constraints)
    constraints = [test_constraints[i] for k = 1:length(UI)]
    @test resample(UI, constraints, StrictlyDecreasing()) isa Vector{Float64}
    @test all([resample(UI, constraints, StrictlyDecreasing()) isa Vector{Float64} for k = 1:n_realizations])
end

iv = UncertainIndexValueDataset(UI, UV)
@test resample(iv, StrictlyDecreasing()) isa Tuple{Vector{Float64}, Vector{Float64}}

# First constrain using a single regular constraint, then apply the order constraint.
for i = 1:length(test_constraints)
    constraint = test_constraints[i]
    @test resample(iv, constraint, StrictlyDecreasing()) isa Tuple{Vector{Float64}, Vector{Float64}}
    @test resample(iv, constraint, constraint, StrictlyDecreasing()) isa Tuple{Vector{Float64}, Vector{Float64}}
    @test all([resample(iv, constraint, StrictlyDecreasing()) isa Tuple{Vector{Float64}, Vector{Float64}} for k = 1:n_realizations])
    @test all([resample(iv, constraint, constraint, StrictlyDecreasing()) isa Tuple{Vector{Float64}, Vector{Float64}} for k = 1:n_realizations])

    cs = [test_constraints[i] for k = 1:length(UI)]
    @test resample(iv, cs, cs, StrictlyDecreasing()) isa Tuple{Vector{Float64}, Vector{Float64}}
    @test resample(iv, cs, constraint, StrictlyDecreasing()) isa Tuple{Vector{Float64}, Vector{Float64}}
    @test resample(iv, constraint, cs, StrictlyDecreasing()) isa Tuple{Vector{Float64}, Vector{Float64}}
end